# An application of node and edge nonlinear hypergraph centrality to a protein
complex hypernetwork

Sarah Lawson 1 Corresponding author Email addresses<EMAIL_ADDRESS>(Sarah Lawson<EMAIL_ADDRESS>(Diane Donovan<EMAIL_ADDRESS>(James Lefevre). School of Mathematics and Physics, ARC Centre of Excellence,
Plant Success in Nature and Agriculture, The University of Queensland,
Brisbane, 4072, Australia Diane Donovan 1 School of Mathematics and Physics,
ARC Centre of Excellence, Plant Success in Nature and Agriculture, The
University of Queensland, Brisbane, 4072, Australia James Lefevre  School of
Mathematics and Physics, ARC Centre of Excellence, Plant Success in Nature and
Agriculture, The University of Queensland, Brisbane, 4072, Australia

###### Abstract

The use of graph centrality measures applied to biological networks, such as
protein interaction networks, underpins much research into identifying key
players within biological processes. This approach however is restricted to
dyadic interactions and it is well-known that in many instances interactions
are polyadic. In this study we illustrate the merit of using hypergraph
centrality applied to a hypernetwork as an alternative. Specifically, we
review and propose an extension to a recently introduced node and edge
nonlinear hypergraph centrality model which provides mutually dependent node
and edge centralities. A Saccharomyces Cerevisiae protein complex hypernetwork
is used as an example application with nodes representing proteins and
hyperedges representing protein complexes. The resulting rankings of the nodes
and edges are considered to see if they provide insight into the essentiality
of the proteins and complexes. We find that certain variations of the model
predict essentiality more accurately and that the degree-based variation
illustrates that the centrality-lethality rule extends to a hypergraph
setting. In particular, through exploitation of the models flexibility, we
identify small sets of proteins densely populated with essential proteins.

###### keywords:

eigenvector centrality, protein complex, hypergraph

###  1 Introduction

## Chapter 0

The origins of graph theory date back to the mid 18th century [1] with the
field advancing over the next three centuries to become an area of significant
interest drawing on and influencing many other areas of mathematics. One of
the many notable advances was the introduction of hypergraphs, formally
introduced by Berge in 1970[2]. The (hyper)edges of a hypergraph can contain
an arbitrary number of vertices, as opposed to a graph whose edges consist of
at most two vertices, thus a hypergraph can model more complex systems. The
analysis of large and complex datasets represented by hypergraphs, also known
as hypernetworks in the applied setting, has led to numerous intriguing
results in recent years. Examples include the analysis of polyadic
interactions between proteins in humans [3] and multi-modal biomarkers for
Alzheimer’s disease [4].

Following the introduction of hypergraphs, researchers started to look at
extending graph theoretical concepts to provide insights into this new area.
One such area of interest is centrality measures which rank vertices and edges
by some defined ‘importance’ or ‘centrality’. This has been extensively
researched in the context of graphs with each measure having their own merit
based on what constitutes as ‘important’ in the given problem. For example,
degree centrality ranks vertices by their degree in descending order. Of
particular interest is eigenvector centrality which “accords each vertex a
centrality that depends on both the number and quality of its connections”
[5]. In 1972, Bonacich [6] introduced graph eigenvector centrality utilising
the eigenvectors of the adjacency matrix. When ordered by descending size, the
entries in the eigenvector give a ranking of the graph vertices. However,
although Bonacich is held to have introduced the measure formally, the concept
had appeared earlier in different disciplines. Edmund Landaus’s ranking of
chess players in 1895 [7] is thought to be the first documented usage of an
eigenvector centrality type measure.

Eigenvector centrality can also be applied to the ‘importance’ of edges and
further extended to the ‘importance’ of edges and vertices being mutually
dependent on each other. This latter approach was introduced by Bonacich
through “simultaneous eigenvector centrality scores”[8]. Faust considered
these centralities in the context of affiliation networks, a social network
which consists of actors and subsets of actors called events[9]. These
centrality scores were the foundation of linear eigenvector centrality for
uniform hypergraphs which was formally defined in 2004 [10]. Aiming to
overcome the limitations of a linear centrality, Benson [11] proposed a non-
linear eigenvector centrality for uniform hypergraphs utilising adjacency
tensors. In 2021, Tudisco and Higham [12] took this further and introduced a
model which not only combined the aforementioned centralities but extended
them to nonlinear centralities for nonuniform hypergraphs. It is this model
that we focus on in this research. For ease we refer to their model as the
Nonlinear Singular Vector Centrality (NSVC) Model aligning with the author’s
reference to their centrality vector solutions as ‘Nonlinear Perron Singular
Vectors’. As detailed in Section 2, the model incorporates a set of four
functions, $\func$, the definitions of which dictate different node and edge
centralities. (Both the terms vertex and node are in common use. The former
more often in the context of theoretical hypergraphs and the latter in applied
hypernetworks. They are used interchangeably in this paper.)

When defining an eigenvector type centrality measure, one of the key
considerations is the conditions under which a unique positive eigenvector
exists for the defined problem. Classical Perron Frobenius Theorem (PFT)
concerns itself with proving the existence of unique, positive eigenvectors
for a positive matrix as proposed by Perron in 1907 [13] and extended by
Frobenius in 1912 to positive eigenvectors of a nonnegative irreducible matrix
[14]. In recent years, PFT has been extended to many areas. Of particular
interest is its extension to tensors [15] and multihomogeneous mappings [16],
the former being used to evidence the uniqueness and existence of certain
tensor eigenvector centralities for uniform hypergraphs [11] and the latter
used by Tudisco and Highams to prove their existence and uniqueness theorem
under certain conditions. Reference to the existence and particularly the
uniqueness of a solution is surprisingly often omitted when eigenvector
centrality is applied to real-world scenarios.

One of the primary objectives of this work is to showcase the adaptability of
the NSVC model and it’s potential as an analytic tool to identify important
nodes and edges in hypergraphs representing biological networks. We apply the
model to a hypernetwork representing protein complexes in Saccharomyces
Cerevisiae, commonly known as brewers yeast, and demonstrate that the
resulting rankings can be used to predict essential proteins and complexes.
Essential proteins carry out essential cellular functions and the deletion of
such proteins leads to cell inviability. A protein complex is a group of
proteins interacting with each other to enact cellular functions. Definitions
of an essential protein complex vary but are often based on the number or
proportion of essential proteins within the complex. This is discussed further
in Section 3. Protein interactions are often represented as a dyadic Protein
Interaction Network (PIN) where the nodes represent proteins and the edges
represent the interaction between two proteins. In a 2001 article, Jeong et al
[17] stated “the most highly connected proteins in the cell are the most
important for its survival”. This is known as the centrality-lethality rule,
namely that nodes with higher degree centrality are more likely to be
essential. This short paper has been cited over 6000 times and has motivated
many researchers in the following decades to utilise centrality as a tool to
identify essential proteins in PINs. One of the earlier papers building on
this idea was by Estrada in 2006 [18] who applied the degree, closeness,
betweenness, eigenvector, information and subgraph centralities to a yeast PIN
and showed that each performed better than random selection in identifying
essential proteins. Extensive progress improving on these preliminary results
has been made over the years not only incorporating network topology but also
including biological information and machine learning techniques to classify
essential proteins. Numerous comprehensive reviews are available detailing
these approaches such as [19, 20, 21, 22]. Many of these approaches have
incorporated protein complex information into their centrality measures, such
as ECC [23], PCSD [24], LBCC [25], LIDC [26] and UC [27]. In this work we
propose the use of hypergraph centrality measures applied to a protein complex
hypernetwork (PCH) with nodes representing proteins and edges representing
complexes [28, 3, 29] as an alternative to PIN-based approaches. Using a PCH
retains information that would otherwise be lost in a PIN.

In this paper we first consider what it means for a node or edge to rank
higher in each of Tudisco and Higham’s proposed function sets for their model
and then we propose a further generalisation. We then apply this generalised
model to the S.Cervisiae PCH as an example of its use. We consider whether the
resulting rankings place nodes and edges representing essential proteins and
complexes higher than those representing nonessential proteins and complexes.
Ultimately, we show that the NSVC model can be customised to improve its
predictive capabilities and provide an alternative platform for future
researchers to build on.

###  1 Methods

## Chapter 1

####  1 Hypergraph definitions

A hypergraph is a tuple $H=(V,E)$ where $V$ is a set of vertices and $E$ is a
set of (hyper)edges. An edge in an undirected hypergraph is defined as a set
$e=\\{v_{1},\dots,v_{k}\\}\subseteq V$ for some $1\leq k\leq|V|$. A k-uniform
hypergraph is a hypergraph where every edge contains exactly $k$ vertices.

A hyperpath in a hypergraph $H$ is an ordered list of vertices and edges
$(v_{0},e_{1},v_{1},e_{2},...,e_{n},v_{n})$ such that $v_{i-1},v_{i}\in e_{i}$
for $1\leq i\leq n$ and $v_{j}\neq v_{k}\ \text{ for }0\leq j<k\leq n$. Edges
are not necessarily distinct. A hypergraph $H$ is connected if for all
$v_{i},v_{j}\in V$ there exists a hyperpath between $v_{i}$ and $v_{j}$ in
$H$.

Two vertices are adjacent if they belong to the same edge and two edges are
adjacent if they have a non-empty intersection. A vertex $i$ in an edge $e$ is
said to be incident to $e$. Similarly, an edge $e$ containing a node $i$ is
said to be incident to $i$. An incidence matrix for a hypergraph $H=(V,E)$
where $|V|=n$ and $|E|=m$ is an $n\times m$ matrix $B=[b_{ie}]$ such that:

$b_{ie}=\begin{cases}1,\text{ if }i\in e,\\\ 0,\text{ otherwise}.\\\
\end{cases}$

Given a hypergraph $H$, the weighted node degree $d^{V}_{i}$ of a vertex $i\in
V$ is $d^{V}_{i}=\sum\limits_{e\in E}b_{ie}w(e)$ [2]. If a hypergraph is
unweighted then $d^{V}_{i}=\sum\limits_{e\in E}b_{ie}$, the number of edges to
which vertex $i$ belongs. The weighted edge degree $d^{E}_{e}$ of an edge
$e\in E$ is $d^{E}_{e}=\sum\limits_{i\in V}b_{ie}v(i)$. If a hypergraph is
unweighted then $d^{E}_{e}=\sum\limits_{i\in V}b_{ie}$, the cardinality of the
edge.

####  2 Tudisco and Higham’s NSVC model [12]

Consider a hypergraph $H$ with incidence matrix $B$, functions
$\func:\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$ acting component wise and
diagonal matrices $N$ and $W$ containing hypergraph node and edge weights
respectively. The node and edge Nonlinear Singular Vector Centralities
$\bm{x}$, $\bm{y}>0$ respectively together with associated scalars
$\lambda,\mu>0$ are a solution of:

$\begin{cases}\lambda\bm{x}=g(BWf(\bm{y}))\\\
\mu\bm{y}=\psi(B^{T}N\varphi(\bm{x}))\\\ \end{cases}$ (1)

with $||\bm{x}||=||\bm{y}||=1$. The $i_{th}$ entry of $\bm{x}$ represents the
centrality $x_{i}$ of node $i$ and the $e_{th}$ entry of $\bm{y}$ represents
the centrality $y_{e}$ of edge $e$. A node and edge ranking results from this
model by ordering the nodes and edges by centrality in descending order. The
use of the 1-norm when normalising $\bm{x}$ and $\bm{y}$ gives the proportion
of centrality for each node or edge.

Tudisco and Higham calculate these centralities through a non-linear power
method. We adapt their code provided in [12] for our calculations.

######  Theorem 0.1.

Existence and uniqueness of solutions to the NSVC model [12]

Let $\func:\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$ be order preserving and
homogeneous functions of degrees $\alpha,\beta,\gamma,\delta$ respectively.
Define the coefficient $\eta$ = $|\alpha\beta\gamma\delta|$. If either:

P1. $\eta<1$,

or

P2. $\eta$=1, $\func$ are differentiable and positive functions and the
bipartite graph with adjacency matrix $\begin{bmatrix}0&BW\\\ B^{T}N&0\\\
\end{bmatrix}$ is connected,

then there exists a unique solution $\bm{x},\bm{y}>0$ (up to scaling) and
$\lambda,\mu>0$ for Equation 1.

Interested readers are encouraged to refer to [12] for the background
information and proof of Theorem 0.1.

Tudisco and Higham give three examples of function sets for their model. The
following summarises these sets including a brief review of their resulting
rankings. A generalised function set is subsequently proposed which
incorporates two of these three function sets and allows for further
variations of the model to better fit a proposed use.

##### A linear function set

$f(\bm{x})=\bm{x},\ g(\bm{x})=\bm{x},\ \varphi(\bm{x})=\bm{x},\
\psi(\bm{x})=\bm{x}$.

The resulting centrality measures after inputting these functions into
Equation 1 encompass the linear node and edge centralities based on the
incidence matrix for graphs proposed by Bonacich [8] and its extension to
hypergraphs detailed in Bonacich et al [10]:

$\lambda\bm{x}=BW\bm{y}\Rightarrow\lambda x_{i}=\sum_{\begin{subarray}{c}e\in
E\\\ i\in e\end{subarray}}w(e)y_{e},\quad\quad\lambda\neq 0,$ (2)

$\mu\bm{y}=B^{T}N\bm{x}\Rightarrow\mu y_{e}=\sum_{\begin{subarray}{c}i\in
e\end{subarray}}v(i)x_{i},\quad\quad\mu\neq 0.$ (3)

In this model variation, a node which belongs to more and/or higher ranked
and/or higher weighted edges ranks higher than one in fewer and smaller ranked
and weighted edges. It is easily seen that $\func$ are order preserving and
homogeneous. Since $\rho$ = $|\alpha\beta\gamma\delta|=1$ it does not meet P1,
however $\func$ are differentiable and positive and thus if the bipartite
graph with adjacency matrix $\begin{bmatrix}0&BW\\\ B^{T}N&0\\\ \end{bmatrix}$
is connected then P2 is met and it follows that there exists unique
$\bm{x},\bm{y}$ (up to scaling) and unique $\lambda,\mu>0$ solution of the
NSVC model.

##### A log-exp function set

$f(\bm{x})=\bm{x},\ \ g(\bm{x})=\bm{x}^{\frac{1}{p+1}},\ \
\varphi(\bm{x})=ln(\bm{x}),\ \psi(\bm{x})=e^{(\bm{x})}$ for
$p\in\mathbb{Z}^{+}$, $p>1$.

For a connected $k$-uniform hypergraph, this function set encompasses tensor
eigenvector centralities introduced in [11]. In particular a $Z$-tensor
eigenvector node centrality for $p=1$ and a $H$-tensor eigenvector node
centrality for $p=k-1$.

Substituting $\func$ into Equation 1 gives:

$\lambda\bm{x}=(BW\bm{y})^{\frac{1}{p+1}}\Rightarrow\lambda
x_{i}=(\sum_{\begin{subarray}{c}e\in E\\\ i\in
e\end{subarray}}w(e)y_{e})^{\frac{1}{p+1}},\quad\lambda\neq 0,$ (4)

$\mu\bm{y}=e^{(B^{T}Nln(\bm{x}))}\Rightarrow\mu y_{e}=e^{(\sum\limits_{j\in
e}v(j)ln(x_{j}))}=\prod\limits_{j\in e}x_{j}^{v(j)},\quad\mu\neq 0.$ (5)

As with the linear set, Equation 4 indicates that the more and/or higher
ranked and/or weighted edges a node belongs to, the higher the node will rank.
However, Equation 5 indicates that an edge will rank higher if it contains
fewer and/or higher ranked and/or lower weighted nodes due to the effect of
$||\bm{x}||=1$.

These functions are an interesting choice as the range of $\varphi$ for the
domain $\mathbb{R}^{+}$ does not fall within the specified co-domain
$\mathbb{R}^{+}$. Further, neither $\varphi$ or $\psi$ are homogeneous thus
the existence and uniqueness of a solution to the NSVC model with these
functions cannot be guaranteed using Theorem 0.1. However, it is known that
for connected $k$-uniform hypergraphs, positive H and Z tensor eigenvectors
exist and uniqueness is guaranteed for a H eigenvector [15, 30]. Further
research into alternative existence and uniqueness of solutions for functions
not meeting the required conditions would be of interest, particularly in the
context of non-uniform hypergraphs.

##### A max function set

$f(\bm{x})=(\bm{x}),\ \ g(\bm{x})=\bm{x},\ \varphi(\bm{x})=\bm{x}^{\omega},\
\psi(\bm{x})=\bm{x}^{\frac{1}{\omega}}$ for some
$\omega\in\mathbb{Z}^{+}\setminus\\{0\\}.$

The node centrality equations are given in Equation 2. Substituting $\varphi$
and $\psi$ into Equation 1 gives edge centralities:

$\displaystyle\mu\bm{y}$
$\displaystyle=(B^{T}N\bm{x}^{\omega})^{\frac{1}{\omega}}\Rightarrow\mu
y_{e}=(\sum\limits_{j\in
e}(v(j)^{\frac{1}{\omega}}x_{j})^{\omega})^{\frac{1}{\omega}},\quad\mu\neq 0,$
(6)

which for an unweighted hypergraph and for large enough $\omega$ gives $\mu
y_{e}\approx max\\{x_{j_{1}},\dots,x_{j_{k}}\\}$.

Again, the node centrality is linear, the more and/or higher ranked edges a
node belongs to, the higher the node will rank. But for large enough $\omega$,
an edges rank is only dependent on the maximum of the centralities of its
incident nodes.

The same conditions as the linear case apply for existence of a unique
$\bm{x},\bm{y}$ (up to scaling) and unique $\lambda,\mu>0$ solution of the
NSVC model.

##### Generalised power function set

We now propose the following generalisation:

$f(\bm{x})=\bm{x}^{a},\ g(\bm{x})=\bm{x}^{b},\ \varphi(\bm{x})=\bm{x}^{c},\
\psi(\bm{x})=\bm{x}^{d}$ for $a,c\in\mathbb{R}_{\geq 0}\ \text{ and
}b,d\in\mathbb{R}^{+}$.

We will refer to this set of functions as $P:(a,b,c,d)$. Substituting these
functions into the NSVC model we get:

$\lambda\bm{x}=(BW\bm{y}^{a})^{b}\Rightarrow\lambda
x_{i}=(\sum_{\begin{subarray}{c}e\in E\\\ i\in
e\end{subarray}}w(e)y_{e}^{a})^{b},\quad\quad\lambda\neq 0,$ (7)
$\mu\bm{y}=(B^{T}N\bm{x}^{c})^{d}\Rightarrow\mu
y_{e}=(\sum_{\begin{subarray}{c}i\in
e\end{subarray}}v(i)x_{i}^{c})^{d},\quad\quad\mu\neq 0.$ (8)

from which follows

$\mu^{a}\lambda x_{i}=(\sum_{\begin{subarray}{c}e\in E\\\ i\in
e\end{subarray}}w(e)(\sum_{\begin{subarray}{c}j\in
e\end{subarray}}v(j)x_{j}^{c})^{da})^{b},$ (9)

$\mu\lambda^{c}y_{e}=(\sum_{\begin{subarray}{c}i\in
e\end{subarray}}v(i)(\sum_{\begin{subarray}{c}f\in E\\\ i\in
f\end{subarray}}w(f)y_{f}^{a})^{bc})^{d}.$ (10)

The flexibility to choose variables $a,b,c,d$ allows the model to be
customised to suit the application. By Theorem 0.1 there exists a unique
solution to functions of this form if $|abcd|<1$, or if $|abcd|=1$ and the
conditions of P2 are met.

It is clear that this function set incorporates the linear and max functions
defined earlier with $P:(1,1,1,1)$ and $P:(1,1,\omega,1/\omega)$ respectively.
Significantly, it also encompasses the centralities based on the node and edge
degrees with $P:(0,1,0,1)$. This is of key interest in the application of the
model to the S.Cervisiae dataset due to degree centrality motivating the
investigation of centrality measures as a tool to predict protein essentiality
for PINS [18].

Although there are infinitely many possible sets of $a,b,c,d$ they can be
combined into families when considering using the centralities as a ranking
mechanism. For example, since node and edge degree centralities are not
dependent on each other it can easily be seen that any function set of the
form $P:(0,b,0,d)$ gives the same node and edge rankings as $P:(0,1,0,1).$
Further, if we are to just consider nodes then any variation of the form
$P:(0,b,c,d)$ results in the same node rankings. Equivalently, any function
sets of the form $P:(a,b,0,d)$ result in the same edge rankings. Another
example is those function sets of the form $P:(1,0,c,d)$ and $P:(a,b,1,0)$
which result in the same node and edge rankings respectively.

An interesting combination is P:($\omega,1/\omega,\omega,1/\omega$) which
gives:

$\lambda\bm{x}=(BW\bm{y}^{\omega})^{1/\omega}\Rightarrow\lambda
x_{i}=(\sum_{\begin{subarray}{c}e\in E\\\ i\in
e\end{subarray}}w(e)y_{e}^{\omega})^{1/\omega},\quad\quad\lambda\neq 0,$ (11)

$\mu\bm{y}=(B^{T}N\bm{x}^{\omega}))^{1/\omega}\Rightarrow\mu
y_{e}=(\sum_{\begin{subarray}{c}i\in
e\end{subarray}}v(i)x_{i}^{\omega})^{1/\omega},\quad\quad\mu\neq 0.$ (12)

For large enough $\omega$ and an unweighted hypergraph we get

$\mu y_{e}\approx max\\{x_{j_{1}},\dots,x_{j_{k}}\\}\ \text{ where
}e=\\{j_{1},\dots j_{k}\\},$ (13) $\lambda x_{i}\approx
max\\{y_{e_{1}},\dots,y_{e_{l}}\\}\text{ where }i\in\\{e_{1},\dots,e_{l}\\}.$
(14)

In our application we consider $\omega=95$ as this value leads to a close
approximation for the maximum values whilst still maintaining efficiency. As
$\omega$ decreases, the contributions of the smaller node and edge
centralities to the edge and node centralitites respectively in Equations 11
and 12 increase until the linear case when $\omega=1$ and all centralities
contribute their exact value. Further, as $\omega\rightarrow 0$ the node and
edge centralities tend to unit contributions to the edge and node
centralitites respectively equalling the ranking given by $P:(0,b,0,d)$. This
motivates the approach of defining a parameter space for $P:(a,b,c,d)$, for
example $a=c=\omega$ and $b=d=1/\omega$ where $\omega\in$(0, 95] for
P:($\omega,1/\omega,\omega,1/\omega$), and performing a search over this space
to find parameters $a,b,c,d$ which permit a better performing model for a
particular set of data.

####  3 NSVC application to S.Cerevisiae

##### S. Cerevisiae Dataset

The S.Cerevisiae protein complex dataset was obtained from [26]. The data
contained in the authors dataset, Complex_745, combines four existing protein
complex sets (CM270 [31], CM425 [31] [32] [33], CYC408 and CTC428 both from
CYC2008 [34], [35]). Proteins that do not belong to a protein complex are not
included in this analysis. Protein essentiality data was obtained from [25]
and originated from MIPS [31], SGD [33], DEG [36] and SGDP [37]. Three single
member complexes were removed from the original dataset with the resulting
dataset containing 742 unique complexes of cardinality 2 or greater and 2166
proteins. The proteins are represented by the hypergraph nodes and the protein
complexes by the hypergraph edges. To align with the existence and uniqueness
conditions detailed in Theorem 0.1, we restrict our analysis to the largest of
113 connected components. This component contains 1739 proteins, 80% of the
original proteins, and 558 protein complexes, 75% of the original complexes.

Unlike protein essentiality, there is no one defined measure of “complex
essentiality”. One approach has been to consider a (human) protein complex as
“potentially” essential if it contains at least one essential protein [3].
Alternatively, in [38], the analysis of the core proteins gives rise to
another definition, that a specific complex is essential if at least 60% of
the core proteins are essential. Influenced by the latter definition we
consider complexes as a whole which contain at least 60% essential nodes as
essential but recognise that this is a simplification. In subsequent works,
consideration of the complex core would likely yield interesting results.

The dataset contained 46 proteins with unknown essentiality, these were
excluded reducing the set to 1693 proteins and 558 complexes. The final set
contains 693 essential proteins and 230 essential complexes. For simplicity,
we refer to the hypergraph nodes representing essential proteins as essential
nodes and equivalently essential edges are those representing essential
complexes.

##### Evaluation methods

The S.Cerevisiae PCH incidence matrix $B$ along with specified $\func$ are
used to calculate centralitites through Equation 1. There are no defined
weights thus $N$ and $W$ are both identity matrices. The nodes and edges are
ranked in descending order of their resulting centralities.

To evaluate the performance of the different function sets on essential
protein identification, the $t$ top ranked nodes are first identified. The
value of $t$ is based on a percentage of the total number of essential
proteins, namely 5%, 10%, 25%, 50%, 75% and 100% of the total number of
essential proteins. The ranking is then analysed against the hypothesis that a
node $i$ with rank $r(i)$ is essential if $r(i)\leq t$. Given a dataset with
$n$ essential proteins, a perfect classification model would have $r(i)\leq n$
for all essential nodes and a nonessential node would have $r(i)>n$. If an
essential node is correctly classified as essential then this is considered a
true positive (TP). A false positive (FP) is when a node $i$ representing a
nonessential protein has $r(i)\leq t$ thus is classified as essential. An
essential node with $r(i)>t$ gives a false negative (FN) and a nonessential
node with $r(i)>t$ gives a true negative (TN). An analogous approach is taken
for analysing the performance of the chosen function sets when classifying
essential edges.

As we are predominantly interested in essential nodes and edges being ranked
higher we restrict our analysis to measures involving true positives. The
precision of the model is the ratio between the number of correctly classified
essential proteins or complexes and the total number of essential proteins or
complexes: $\frac{TP}{TP+FP}$. The recall is the ratio between the number of
correctly classified essential proteins or complexes and the total number of
proteins or complexes classified as essential: $\frac{TP}{TP+FN}$. We map
these measures against each other for all values of $1\leq t\leq 1693$
resulting in the precision recall curve.

The area under the precision-recall curve (PR-AUC) gives a single measure to
compare classifiers. The better a classifier performs, the higher the PR-AUC
will be. A perfect classifier will have a PR-AUC of 1. We also consider the
cumulative frequency of true positives for $1\leq t\leq 1693$. In both the
precision recall curve and cumulative frequency graphs the ratio of essential
to total number of proteins/complexes is used to define a baseline
representing a random classifier to compare the performance of the chosen
function sets.

These methods are similar to those used in papers such as [26] and its
successors when identifying essential proteins in PINs.

###  1 Results

## Chapter 2

First we analyse the performance of the model with selected function sets for
six thresholds. The overall performance is then considered in Section 1.
Finally the associated classifications of edges is analysed in Section 2.

####  1 Performance of the NSVC model for essential protein classification in
S.Cervisiae

##### Classification performance for six thresholds

The six thresholds used are $t\in\\{34,69,173,346,519,693\\}.$ These values
are based on 5%, 10%, 25%, 50%, 75% and 100% of 693, the number of essential
proteins in the dataset.

Set A: Base functions
To investigate the centrality-lethality rule we apply the node and edge degree
centralities using $P:(0,1,0,1)$ to the S.Cervisiae PCH. Tudisco and Higham’s
max function set, redefined as $P:(1,1,95,1/95)$, the linear function set,
redefined as $P:(1,1,1,1)$ and the log-exp function are also applied to
investigate if essential nodes match the criteria required to be ranked higher
using these functions. For ease we use $p=14$ in the log-exp function set,
giving $g(x)=\bm{x}^{1/15}$, as this ensures a less skewed distribution of the
centrality sizes towards the extremes of the interval [0,1] than for smaller
values of $p$. The function set $P:(1,1,46,29/46)$ is also included, the
motivation for the inclusion of this set is discussed within the following
results.

Figure 1 shows that using the function set $P:(0,1,0,1)$ predicts more
essential proteins than $P:(1,1,1,1)$ and the log-exp function sets for the
nominated thresholds. However the function set $P:(1,1,95,1/95)$ performs
similarly to $P:(0,1,0,1)$ for the top 69, 173 and 346 proteins thus we
consider building on this function set to surpass the performance of
$P:(0,1,0,1)$. As detailed in Section 2, $P:(1,1,95,1/95)$ provides a node
centrality that is proportional to a linear combination of the centralities of
the edges to which it belongs and an edge centrality based only on the maximum
centrality of its incident nodes. Utilising Equation 9 we see that
$P:(1,1,95,1/95)$ gives

$\mu\lambda x_{i}=\sum_{\begin{subarray}{c}e\in E\\\ i\in
e\end{subarray}}(\sum_{\begin{subarray}{c}j\in V\\\ j\in
e\end{subarray}}x_{j}^{95})^{1/95}.$

Since $\omega=95$ is large enough this gives:

$\mu\lambda x_{i}\approx\sum_{\begin{subarray}{c}e\in E\\\ i\in
e\end{subarray}}max\\{x_{j1},\dots,x_{jk}\\}\text{ where
}e=\\{x_{j1},\dots,x_{jk}\\}.$

However, if a node’s centrality is not based on the maximum ranked node in
each adjacent edge but on the centralities of a set of the top ranking nodes,
then this motivates the use of function sets of the form $P:(1,1,c,k/c)$. To
investigate this we considered $P:(1,1,c,k/c)$ where $c,k\in\\{0.1,0.2,\dots
0.9\\}\cup\\{1,2,\dots 95\\}$. The top performing function sets for each of
the six thresholds were identified and those which were also top performers
for predicting essential complexes were singled out. From this set, the
function set $P:(1,1,46,29/46)$ was highlighted as it had the highest node PR-
AUC. This set exceeds the number of essential proteins predicted by
$P:(0,1,0,1)$ and $P:(1,1,95,1/95)$ for $t\in\\{69,173,346\\}$ (Figure 1).

Figure 1: Number of essential proteins correctly classified for Set A function
sets for six thresholds.

Set B: Combined functions
A second approach taken was to combine and extend the degree, linear and max
functions through an exhaustive search of the combinations for $P:(a,b,c,d)$
with $a,c\in\\{0,1,1/95,95\\}$ and $b,d\in\\{1,1/95,95\\}$ to see if any of
these combinations performed better than the set A functions. Included in
these functions are the intuitive approaches summarised in Table 1. These
functions do not necessarily align with underlying characteristics of the PCH
dictating essentiality but they are being showcased here as important
considerations for other applications. For example, $P:(0,1,95,1/95)$ ranks
nodes by node degree centrality, utilising the number of edges a node belongs
to, whilst simultaneously ranking edges by max centrality, utilising the
maximum centrality of the nodes in the edge.

Table 1: A subset of functions included in Set B Function Set | Node centrality type | Edge centrality type
---|---|---
P:(95,1/95,95,1/95) | max | max
P:(95,1/95,0,1) | max | degree
P:(95,1/95,1,1) | max | linear
P:(0,1,95,1/95) | degree | max
P:(0,1,1,1) | degree | linear
P:(1,1,0,1) | linear | degree
P:(1,1,95,1/95) | linear | max

For functions of this form, $P:(1/95,95,1,1/95)$ was identified as being a top
performing function set using the same approach as that which identified
$P:(1,1,46,29/46)$ in set A. Figure 2 shows the results for the function sets
summarised in Table 1 together with $P:(1/95,95,1,1/95)$. As expected from the
discussion in Section 2, $P:(0,1,1,1)$ and $P:(0,1,95,1/95)$ produce the same
node rankings as $P:(0,1,0,1)$ as they are of the form $P:(0,b,c,d)$.

Figure 2: Number of essential proteins correctly classified for a subset of
Set B functions.

##### Overall classification performance

We now compare the overall performance of the top performers described in
Section 1 alongside $P:(0,1,0,1)$ and $P:(1,1,95,1/95)$.

Figure 3 shows clearly that $P:(0,1,0,1)$ and $P:(1,1,95,1/95)$ perform
approximately equal at first but then $P:(1,1,95,1/95)$ drops off
significantly and predicts less than or equal to random for larger values of
$t$. This later poorer performance for $P:(1,1,95,1/95)$ results in a PR-AUC
of 0.45, only marginally above the random PR-AUC of 0.40. The PR-AUC for
$P:(0,1,0,1)$ is higher at 0.51. Function set $P:(1/95,95,1,1/95)$ performs
approximately equal to $P:(0,1,0,1)$ reflected in an identical PR-AUC of 0.51.

The performance of $P:(1,1,46,29/46)$, measured by a PR-AUC of 0.52, is only
marginally better than $P:(0,1,0,1)$ however its performance for smaller
values of $t$ is significantly better (Figures 3 and 4.). It is this that sets
it aside from the others. Being able to identify small sets of predominantly
essential proteins could be useful in many small scale experiments.

The $P:(1,1,1,1)$ function set has the lowest PR-AUC of 0.38 with predictive
capabilities less than random in the upper thresholds. The overall poor
performance of this function set indicates that the network characteristics of
the essential proteins in the PCH do not align well with those influencing a
higher ranking under $P:(1,1,1,1)$. Specifically, the essentiality of a
protein cannot be accurately predicted solely by the number and importance of
the complexes it belongs to (Equation 2) or by the number and centrality of
the proteins in those complexes (Equation 9).

Figure 3: Cumulative frequency of the number of essential proteins correctly
classified via the top $t$ ranked nodes. Figure 4: Comparison of the precision
versus recall for essential protein classification.

####  2 Performance of the NSVC model for essential edge classification

The focus so far has been primarily on node classification however an
advantage of the NSVC model is that it simultaneously provides edge rankings.
In this section we briefly look at the performance of select functions
identified for node classification and interpret their performance when they
are applied to classify edges.

The six thresholds nominated are $t\in\\{11,23,57,115,172,230\\}$
corresponding to 5%, 10%, 25%, 50%, 75% and 100% of 230, the number of
essential complexes in the dataset. We see in Figure 5 that $P:(1,1,46,26/46)$
again performs best for smaller values of $t$.

Figure 5: Number of essential complexes correctly classified for nominated
function sets.

We can see in Figure 6 that the function set $P:(0,1,0,1)$ is the best
performer in the middle ranges and $P:(1/95,95,1,1/95)$ performs best for
larger values of $t$. Overall, the highest PR-AUC is 0.53 resulting from
$P:(1,1,46,29/46)$ followed by 0.47 for $P:(0,1,0,1)$ compared to a random PR-
AUC of 0.41. The precision recall curve, Figure 7, depicts the varying PR-AUC
results with $P:(1,1,46,29/46)$ clearly dominating.

Figure 6: Cumulative frequency of the number of essential complexes correctly
classified via the top $t$ ranked edges. Figure 7: Comparison of the precision
versus recall for essential protein complex classification.

Although it is likely that other variations of $P:(a,b,c,d)$ would yield
additional positive results it is motivating that $P:(1,1,46,29/46)$ predicts
both a higher number of essential proteins and complexes than the other
trialed function sets including the degree based $P:(0,1,0,1)$ for smaller
values of $t$.

###  1 Discussion

## Chapter 3

Research identifying essential proteins through their interactions has been
predominantly focused on using graph measures applied to PINs. Our research
has shown there is merit in using a PCH and applying newer hypergraph
centrality measures, specifically the NSVC model. In our novel application of
this model we identified function sets which performed better than random for
both the classification of essential proteins and essential protein complexes.
We highlighted the above average performance of the node-edge degree
centrality $P:(0,1,0,1)$ evidencing the centrality-lethality rule in a
hypergraph setting and extended the NSVC model to a generalised function set
$P:(a,b,c,d)$. This flexible approach resulted in the identification of
$P:(1,1,46,29/46)$ which gave the best classification performance for a range
of thresholds.

One of the significant advantages of using the NSVC model is that we gain
mutually dependent node and edge centralities. In this work we have mainly
focused on the node centralities and the associated classification of
essential proteins. The next step would be to broaden the results by
refocusing on edge centralities and protein complexes. This could include
alternative definitions of complex essentiality such as essentiality based on
the core of a protein complex as discussed in [38]. A growing body of research
indicates that the essentiality of a protein is linked to the essentiality of
the complex to which it belongs as opposed to being a characteristic solely of
the protein itself [38, 39, 40]. This further motivates shifting the focus
onto edges. Research also shows that essential proteins are more likely to
belong to complexes [39]. This is evident in the dataset used in this research
which is denser in essential proteins than the dataset underlying the PIN used
in [26] for example. It would therefore be useful to reapply graph centrality
measures, such as those in [26], to a PIN based only on proteins contained
within complexes as a comparative exercise.

Recently Tudisco and Higham proposed a core-periphery score vector that
assigns large values to nodes in the core of a hypergraph and small values to
those in the periphery [41]. Although out of scope of this work and requiring
input from proteomic experts, this measure would be interesting to analyse due
to proposed connections between hypergraph cores and essential proteins and
complexes such as those discussed in [29]. Other centralities which provide
insight into a nodes interaction with the edges to which it belongs could also
be considered for application to a PCH. In particular, subgraph centrality
[42], which assigns a node ranking based on closed walks originating from the
node, and vector centrality [43], which gives a vector communicating a nodes
interaction with edges of different cardinality.

Factors influencing the ‘importance’ of a node or edge in a biological
hypernetwork, such as the PCH studied in this work, aren’t necessarily
topological. With input from experts in the biological field, non-topological
information could be incorporated into the NSVC model. For instance, weight
functions could be utilised to represent quantifiable node and edge data or
directed hypergraphs could be used to model interactions more accurately.

In this paper we have showcased an alternative approach to identifying
essential proteins by proposing the use of the NSVC model applied to a PCH.
The potential of extending this application to other biological hypernetworks
is self-evident. It is hoped that the work initiated in this paper motivates
further research into using hypergraph centrality, particularly the NSVC
model, to analyse complex biological processes.

### Acknowledgments

Funding: This work was supported by The Australian Research Council,through
the Centre of Excellence for Plant Success in Nature and Agriculture
(CE200100015).

### Data Availability

Data and code used in this paper is publicly available from the original
sources or via the online repository
https://github.com/LawsonSarah/PCH_node_edge_centrality

## References

  * [1] Biggs NL, Lloyd EK, Wilson RJ.  Graph theory 1736-1936.  Oxford England: Clarendon Press; 1976.
  * [2] Berge C.  Graphs and hypergraphs.  [rev. ed.] / translated by Edward Minieka. ed. Amsterdam: North-Holland Pub. Co. ; New York : American Elsevier Pub. Co.; 1973.
  * [3] Klimm F, Deane CM, Reinert G.  Hypergraphs for predicting essential genes using multiprotein complex data.  Journal of complex networks. 2021;9(2).
  * [4] Shao W, Xiang S, Zhang Z, Huang K, Zhang J.  Hyper-graph based sparse canonical correlation analysis for the diagnosis of Alzheimer’s disease from multi-dimensional genomic data.  Methods (San Diego, Calif). 2021;189:86-94.
  * [5] Newman MEJ.  In: Mathematics of Networks. London: Palgrave Macmillan UK; 2018. p. 8525-33.
  * [6] Bonacich P.  Factoring and weighting approaches to status scores and clique identification.  The Journal of Mathematical Sociology. 1972;2(1):113-20.
  * [7] Landau E.  Zur relativen Wertbemessung der Turnierresultate.  Deutsches Wochenschach. 1895;11:366-9.
  * [8] Bonacich P.  Simultaneous group and individual centralities.  Social networks. 1991;13(2):155-68.
  * [9] Faust K.  Centrality in affiliation networks.  Social Networks. 1997 Apr;19(2):157-91.
  * [10] Bonacich P, Cody Holdren A, Johnston M.  Hyper-edges and multidimensional centrality.  Social networks. 2004;26(3):189-203.
  * [11] Benson AR.  Three Hypergraph Eigenvector Centralities.  SIAM journal on mathematics of data science. 2019;1(2):293-312.
  * [12] Tudisco F, Higham DJ.  Node and Edge Nonlinear Eigenvector Centrality for Hypergraphs.  Communications Physics. 2021 Jan;4:244.
  * [13] Perron O.  Zur Theorie der Matrices.  Math Ann. 1907;64:248–263.
  * [14] Frobenius G.  Ueber Matrizen aus nicht negativen Elementen.  Sitzungsber Königl Preuss Akad Wiss. 1912:456–477.
  * [15] Chang KC, Pearson K, Zhang T.  Perron-Frobenius Theorem for nonnegative tensors.  Communications in mathematical sciences. 2008;6(2):507-20.
  * [16] Gautier A, Tudisco F, Hein M.  The Perron-Frobenius theorem for multi-homogeneous mappings; 2021.  arXiv:1801.05034.  Available from: https://doi.org/10.48550/arXiv.1801.05034.
  * [17] Jeong H, Mason SP, Barabási AL, Oltvai ZN.  Lethality and centrality in protein networks.  Nature. 2001 May;411(6833):41-2.
  * [18] Estrada E.  Virtual identification of essential proteins within the protein interaction network of yeast.  Proteomics. 2006 Jan;6(1):35-40.
  * [19] Ashtiani M, Salehzadeh-Yazdi A, Razaghi-Moghadam Z, Hennig H, Wolkenhauer O, Mirzaie M, et al.  A systematic survey of centrality measures for protein-protein interaction networks.  BMC Systems Biology. 2018 Jul;12(1).
  * [20] Jalili M, Salehzadeh-Yazdi A, Gupta S, Wolkenhauer O, Yaghmaie M, Resendis-Antonio O, et al.  Evolution of Centrality Measurements for the Detection of Essential Proteins in Biological Networks.  Frontiers in Physiology. 2016 Aug;7.
  * [21] Li X, Li W, Zeng M, Zheng R, Li M.  Network-based methods for predicting essential genes or proteins: a survey.  Briefings in Bioinformatics. 2019 Feb;21(2):566-83.
  * [22] Rasti S, Vogiatzis C.  A survey of computational methods in protein–protein interaction networks.  Annals of Operations Research. 2018 Jun;276(1–2):35-87.
  * [23] Ren J, Wang J, Li M, Wang H, Liu B.  Prediction of Essential Proteins by Integration of PPI Network Topology and Protein Complexes Information.  In: Lecture Notes in Computer Science. Springer Berlin Heidelberg; 2011. p. 12-24.
  * [24] Lei X, Yang X.  A new method for predicting essential proteins based on participation degree in protein complex and subgraph density.  PloS one. 2018 Jun;13(6):e0198998.
  * [25] Qin C, Sun Y, Dong Y.  A New Method for Identifying Essential Proteins Based on Network Topology Properties and Protein Complexes.  PloS one. 2016 Aug;11(8):e0161042.
  * [26] Luo J, Qi Y.  Identification of Essential Proteins Based on a New Combination of Local Interaction Density and Protein Complexes.  PloS one. 2015 Jun;10(6):e0131418.
  * [27] Li M, Lu Y, Niu Z, Wu FX.  United Complex Centrality for Identification of Essential Proteins from PPI Networks.  IEEE/ACM Transactions on Computational Biology and Bioinformatics. 2017 Mar;14(2):370-80.
  * [28] Klamt S, Haus UU, Theis F.  Hypergraphs and Cellular Networks.  PLoS Computational Biology. 2009 May;5(5):e1000385.
  * [29] Ramadan E, Tarafdar A, Pothen A.  A hypergraph model for the yeast protein complex network.  In: 18th International Parallel and Distributed Processing Symposium, 2004. Proceedings. IEEE; 2004. p. p.2647-54.
  * [30] Chang KC, Pearson KJ, Zhang T.  Some variational principles for Z-eigenvalues of nonnegative tensors.  Linear algebra and its applications. 2013;438(11):4166-82.
  * [31] Mewes HW.  MIPS: analysis and annotation of proteins from whole genomes in 2005.  Nucleic Acids Research. 2006 Jan;34(90001):D169-72.
  * [32] Aloy P, Böttcher B, Ceulemans H, Leutwein C, Mellwig C, Fischer S, et al.  Structure-Based Assembly of Protein Complexes in Yeast.  Science. 2004 Mar;303(5666):2026-9.
  * [33] Cherry J.  SGD: Saccharomyces Genome Database.  Nucleic Acids Research. 1998 Jan;26(1):73-9.
  * [34] Pu S, Wong J, Turner B, Cho E, Wodak SJ.  Up-to-date catalogues of yeast protein complexes.  Nucleic Acids Research. 2008 Dec;37(3):825-31.
  * [35] Pu S, Vlasblom J, Emili A, Greenblatt J, Wodak SJ.  Identifying functional modules in the physical interactome of Saccharomyces cerevisiae.  Proteomics. 2007 Mar;7(6):944-60.
  * [36] Zhang R, Lin Y.  DEG 5.0, a database of essential genes in both prokaryotes and eukaryotes.  Nucleic Acids Research. 2009 Jan;37(Database):D455-8.
  * [37] Winzeler EA, Shoemaker DD, Astromoff A, Liang H, Anderson K, Andre Bea.  Functional Characterization of the S. cerevisiae Genome by Gene Deletion and Parallel Analysis.  Science. 1999 Aug;285(5429):901-6.
  * [38] Dezső Z, Oltvai ZN, Barabási AL.  Bioinformatics Analysis of Experimentally Determined Protein Complexes in the Yeast Saccharomyces cerevisiae.  Genome Research. 2003 Jan;13(11):2450-4.
  * [39] Hart GT, Lee I, Marcotte EM.  A high-accuracy consensus map of yeast protein complexes reveals modular nature of gene essentiality.  BMC Bioinformatics. 2007 Jul;8(1).
  * [40] Wang H, Kakaradov B, Collins SR, Karotki L, Fiedler D, Shales M, et al.  A Complex-based Reconstruction of the Saccharomyces cerevisiae Interactome.  Molecular and Cellular Proteomics. 2009 jun;8(6):1361-81.
  * [41] Tudisco F, Higham DJ.  Core-Periphery Detection in Hypergraphs.  SIAM Journal on Mathematics of Data Science. 2023 Jan;5(1):1-21.
  * [42] Estrada E, Rodríguez-Velázquez JA.  Subgraph centrality and clustering in complex hyper-networks.  Physica A: Statistical Mechanics and its Applications. 2006 may;364:581-94.
  * [43] Kovalenko K, Romance M, Vasilyeva E, Aleja D, Criado R, Musatov D, et al.  Vector centrality in hypergraphs.  Chaos, Solitons &; Fractals. 2022 Sep;162:112397.


# Heavy-quark potential in Gribov-Zwanziger approach
around deconfinement phase transition

Wan Wu, Guojun Huang, Jiaxing Zhao and Pengfei Zhuang  Physics Department,
Tsinghua University, Beijing 100084, China

###### Abstract

The interaction potential between a pair of heavy quarks is calculated with
resummed perturbation method in Gribov-Zwanziger approach at finite
temperature. The resummed loop correction makes the potential complex. While
the real part is, as expected, screened and becomes short-ranged in hot
medium, the strength of the imaginary part increases with temperature and is
comparable with the real part, which is very different from the previous
calculation in HTL approach. This means that, both the color screening and
Landau damping play important role in the dissociation of heavy flavor hadrons
in hot medium.

Considering the large mass $m$ and small velocity $v$, there exists an
hierarchy of energy scales $m\gg mv\gg mv^{2}$ for heavy quarks. When
integrating out the momentum larger than $m$ and $mv$ from the Quantum
Chromodynamics (QCD) respectively, one obtains the non-relativistic QCD
(NRQCD) and potential non-relativistic QCD (pNRQCD) theories Caswell:1985ui ;
Brambilla:1999xf ; Brambilla:2004jw  for the study of heavy quark systems. At
the leading order, the equation of motion in pNRQCD returns to the Schrödinger
equation in quantum mechanics, and the dynamics is fully described by an
interaction potential. The two-body Schrödinger equation with the Cornell
potential between a pair of heavy quarks successfully describes the quarkonium
properties in vacuum Zhao:2020jqu . The one-gluon exchange between two heavy
quarks gives rise to the Coulomb part of the potential, while the confinement
part should come from non-perturbative calculations through for instance
lattice simulations Koma:2006si ; Kawanai:2011jt .

In recent years heavy-flavor hadrons are widely considered as a probe of the
new state of matter - quark-gluon plasma (QGP) created in relativistic heavy-
ion collisions Matsui:1986dk ; Zhao:2020jqu ; Dong:2019unq ; Rothkopf:2019ipj
; Chapon:2020heu ; Zhao:2020nwy ; Zhao:2016ccp ; He:2014tga ; ExHIC:2017smd .
This extends the study of heavy-quark potential from vacuum to finite
temperature. In a hot medium of light quarks and gluons, the heavy-quark
potential is expected to become complex: the color screening by the
surrounding quarks and gluons reduces the Cornell potential Nadkarni:1986as ;
Wong:2004zr ; Satz:2005hx ; Kaczmarek:2004gv , analogy to the Debye screening
in electromagnetical systems, and the imaginary part is introduced by the
Landau damping or color singlet to octet transition Brambilla:2008cx . At
extremely high temperature, the Landau damping becomes dominant, and the
potential can be well described by Hard-Thermal Loop (HTL) resummed
perturbation Laine:2006ns ; Beraudo:2007ky . At finite temperature, the heavy-
quark potential can be extracted from the quarkonium spectral functions via
lattice QCD simulations Rothkopf:2011db ; Bazavov:2014kva ; Burnier:2014ssa ;
Burnier:2015tda ; Lafferty:2019jpr ; Bala:2019cqu . Recently, the machine
learning method Shi:2021qri  and spectral extraction strategy Larsen:2019zqv ;
Bala:2021fkm  are used to calculate the heavy-quark potential, both indicate a
much larger imaginary part compared with the previous perturbative
calculations Laine:2006ns ; Beraudo:2007ky .

Since Faddeev and Popov quantized the Yang-Mills theory by introducing a new
kind of particles called ghost through the introduction of a gauge
Faddeev:1967fc , the quantization has become a standard textbook item. Some
years later, however, Gribov discovered that the gauge fixing in the
quantization is not complete, there exist still gauge copies called Gribov
copies which could affect the infrared region of those gauge quantities such
as the gluon and ghost propagators Gribov:1977wm . The Gribov action arises
from the restriction of the domain of the integration in Euclidean space to
the Gribov region $\Omega$, which is defined as the set of all gauge field
configurations fulfilling a gauge (for instance the Landau gauge
$\partial^{\mu}A_{\mu}^{a}=0$) and for which the Faddeev-Popov operator
$M^{ab}=-\partial_{\mu}(\delta^{ab}\partial^{\mu}-gf^{abc}A^{\mu}_{c})$ is
strictly positive Gribov:1977wm ; Dokshitzer:2004ie . The original Gribov
Lagrangian includes a nonlinear term and hard to be calculated. Zwanziger used
the BRST method Zwanziger:1982na ; Zwanziger:1989mf ; Zwanziger:1992qr  to
modify the Lagrangian through the introduction of a set of auxiliary fields
and derived the Gribov-Zwanziger (GZ) Lagrangian in local form which is now
used widely Zwanziger:2006sc ; Cucchieri:2000hv ; Gracey:2009mj ;
Canfora:2015yia . From the GZ Lagrangian one can conveniently read off the
gluon propagator. For instance, it reads in Landau gauge Gribov:1977wm  (in
Coulomb gauge see Ref.Burgio:2008jr ),

$D_{\mu\nu}^{ab}(p)=\delta^{ab}\left(\delta_{\mu\nu}-{p_{\mu}p_{\nu}\over
p^{2}}\right){p^{2}\over p^{4}+m_{G}^{4}},$ (1)

where $m_{G}$ is the Gribov mass parameter. In comparison with the normal
gluon propagator in the limit of $m_{G}\to 0$, the GZ propagator is suppressed
in the infrared region by the complex poles at $p^{2}=\pm im_{G}^{2}$. This
structure does not allow us to attach the usual particle meaning to the gluon
propagator, invalidating the interpretation of gluons as excitation of the
physical spectrum. Gluons are confined by the Gribov condition. This means
that, the GZ approach successfully describes the confinement of gluons by the
non-physical singularities in the gluon propagator.

Considering the fact that the confinement is already reflected in the gluon
propagator (1) at lowest order, the GZ approach provides a possibility to
perturbatively calculate the heavy-quark potential. This can provide a way to
understand the physics of the complex potential extracted from the lattice
data. In fact, the calculation at one-loop level in vacuum shows a Coulomb
potential and a linear term modified by a logarithm Golterman:2012dx . To
understand the parton deconfinement at finite temperature, we study in this
paper the heavy quark potential in the GZ approach at finite temperature in
the frame of resummed perturbation theory. We will focus on the region around
the deconfinement temperature $T_{c}$ ($T\leq 2T_{c}$) which can be realized
in high energy nuclear collisions and where the normal perturbation theories
like HTL are not suitable. We calculate firstly the gluon loop, then the gluon
propagator in terms of the resummation of gluon and quark loops, and finally
the heavy quark potential through constructing the real-time Wilson loop. We
summarize in the end.

The Gribov mass parameter $m_{G}$ is not a free parameter of the theory. It is
a dynamical quantity, being determined in a self-consistent way through a gap
equation by minimizing the partition function of the system, $\partial\ln
Z/\partial m_{G}^{2}=0$. At lowest level it can be derived from the
contribution of a closed loop to the gluon self-energy. In Coulomb gauge it
reads Zwanziger:2006sc

$\int{d^{4}p\over(2\pi)^{4}}{1\over(p_{0}^{2}+{\bm{p}}^{2}){\bm{p}}^{2}+m_{G}^{4}}={3\over
2N_{c}g^{2}},$ (2)

where $g(T)$ as a function of temperature is the QCD running coupling
constant, and $N_{c}$ the number of color degrees of freedom. In the imaginary
time formalism of finite temperature field theory, the energy integration
becomes a Matsubara frequency summation $\int d^{4}p/(2\pi)^{4}=T\sum_{n}\int
d^{3}{\bm{p}}/(2\pi)^{3}$ with $p_{0}=i\omega_{n}$ and $\omega_{n}=2n\pi T$.
By evaluating the frequency summation and making the standard $\overline{MS}$
subtraction, the gap equation is simplified as

${1\over 4}\ln\left({e\over 2}{\mu^{2}\over
m_{G}^{2}}\right)+\int_{0}^{\infty}{dx\over u}{1\over
e^{m_{G}u/T}-1}={3\pi^{2}\over N_{c}g^{2}(T)},$ (3)

with $u=\sqrt{x^{2}+1/x^{2}}$, where $\mu$ is the renormalization scale that
controls the vacuum value of the Gribov mass.

The temperature dependence of the Gribov mass is characterized by not only the
statistical distribution but also the coupling constant
$g^{2}(T)=4\pi\alpha_{s}(T)$. At very high temperatures with $T\gtrsim
3T_{c}$, the IR and UV behaviours of the coupling are obtained through lattice
simulations Kaczmarek:2004gv , and the results are qualitatively in consistent
with the perturbative QCD calculation up to two loops Kaczmarek:2004gv , as
shown in Fig. 1. In vacuum the coupling can be extracted by fitting the
quarkonium spectra via potential model Zhao:2020jqu  or fitting the vacuum
potential with Cornell potential Burnier:2015tda , both give a coupling
$\alpha_{s}\approx 0.5$. For the gap around $T_{c}$ where charmonium and
bottonium states are expected to be dissociated and we are interested in in
this paper, the non-perturbative effect is strong, it is hard to find precise
calculation. We then take an interpolation to continuously connect the vacuum
and high-temperature values, see the solid line in Fig. 1.

Figure 1: The running coupling constant $\alpha_{s}(T)$ calculated by lattice
IR and UV Kaczmarek:2004gv , perturbative QCD up to two loops Kaczmarek:2004gv
, potential models Zhao:2020jqu ; Lafferty:2019jpr  and interpolation. The
deconfinement temperature $T_{c}$ is chosen to be $170$ MeV Lafferty:2019jpr .

With the running coupling extracted from the interpolation, we solve the gap
equation numerically and obtain the Gribov mass $m_{G}$ as a function of
temperature, shown in Fig. 2. The vacuum value is chosen to be $m_{G}(0)=0.55$
GeV, corresponding to the renormalization scale $\mu=10.92$ GeV in the gap
equation. We will see in the following that this value can reproduce the
Cornell potential well. The Gribov mass drops down rapidly in the beginning,
when the temperature is below the critical temperature, and then becomes
smooth in the deconfinement phase. In the limit of high temperature
$T\to\infty$, the solution of the gap equation

$m_{G}(T)={N_{c}\over 2^{3/2}3\pi}g^{2}(T)T$ (4)

approaches to a standard magnetic mass $m\sim g^{2}(T)T$ Zwanziger:2006sc ;
Fukushima:2013xsa . In the asymptotic free interval with $g^{2}\to 0$, the
Gribov mass disappears.

Figure 2: The Gribov mass $m_{G}(T)/T$.

To include loop correction to the gluon propagator (1), we calculate first the
gluon loop. Since the interaction potential between two heavy quarks is only
related to the component $D_{00}$ of the gluon propagator $D_{\mu\nu}$, we
consider in the following only the component $\Pi_{00}$ of the loop function
$\Pi_{\mu\nu}(p)$. Defining $\Pi_{G}(p)=-\Pi_{00}(p)/(g^{2}N_{c}{\bm{p}}^{2})$
and using the free gluon propagator (1) in Coulomb gauge, a direct calculation
leads to

$\displaystyle\Pi_{G}(p)$ $\displaystyle=$ $\displaystyle J_{1}(p)-J_{2}(p),$
$\displaystyle J_{1}(p)$ $\displaystyle=$
$\displaystyle{3p_{i}p_{j}\over{\bm{p}}^{2}}\int{d^{4}k\over(2\pi)^{4}}{T_{ij}({\bm{k}})\over({\bm{k}}^{2}+m_{G}^{4}/{\bm{k}}^{2}-k_{0}^{2})({\bm{k}}-{\bm{p}})^{2}},$
$\displaystyle J_{2}(p)$ $\displaystyle=$
$\displaystyle{1\over{\bm{p}}^{2}}\int{d^{4}k\over(2\pi)^{4}}{T_{ij}({\bm{k}})T_{ij}({\bm{k}}-{\bm{p}})\over{\bm{k}}^{2}+m_{G}^{4}/{\bm{k}}^{2}-k_{0}^{2}}$
(5) $\displaystyle\times$
$\displaystyle{{\bm{k}}^{2}+m_{G}^{4}/{\bm{k}}^{2}+k_{0}^{2}-k_{0}p_{0}\over({\bm{k}}-{\bm{p}})^{2}+m_{G}^{4}/({\bm{k}}-{\bm{p}})^{2}-(k_{0}-p_{0})^{2}},$

with $T_{ij}({\bm{k}})=\delta_{ij}-k_{i}k_{j}/{\bm{k}}^{2},\ i,j=1,2,3$. The
above two integrals are divergent in vacuum and need to be regularized.
Introducing a momentum-cutoff $\Lambda$ and using the renormalization scheme
$\overline{MS}$, the divergence can be attracted into the running coupling,

${1\over g^{2}(\Lambda)N_{c}}={1\over
48\pi^{2}}\left(11\log{\Lambda^{2}\over\Lambda_{\overline{MS}}^{2}}-{49\over
3}+22\log 2\right),$ (6)

where the momentum-cutoff in the scheme $\overline{MS}$ is taken as
$\Lambda_{\overline{MS}}\approx m_{G}/1.62$ with the quark contributions
Golterman:2012dx .

$J_{2}$ can further be separated into three parts,

$\displaystyle J_{2}(p)$ $\displaystyle=$ $\displaystyle
J_{2A}(p)+J_{2B}(p)+J_{2C}(p),$ (7) $\displaystyle J_{2A}(p)$ $\displaystyle=$
$\displaystyle{1\over{\bm{p}}^{2}}\int{d^{4}k\over(2\pi)^{4}}{T_{ij}({\bm{k}})T_{ij}({\bm{k}}-{\bm{p}})\over(k_{0}-p_{0})^{2}-\epsilon_{k-p}^{2}},$
$\displaystyle J_{2B}(p)$ $\displaystyle=$
$\displaystyle{1\over{\bm{p}}^{2}}\int{d^{4}k\over(2\pi)^{4}}{2\epsilon_{k}^{2}T_{ij}({\bm{k}})T_{ij}({\bm{k}}-{\bm{p}})\over(k_{0}^{2}-\epsilon_{k}^{2})((k_{0}-p_{0})^{2}-\epsilon_{k-p}^{2})},$
$\displaystyle J_{2C}(p)$ $\displaystyle=$
$\displaystyle{1\over{\bm{p}}^{2}}\int{d^{4}k\over(2\pi)^{4}}{-p_{0}k_{0}T_{ij}({\bm{k}})T_{ij}({\bm{k}}-{\bm{p}})\over(k_{0}^{2}-\epsilon_{k}^{2})((k_{0}-p_{0})^{2}-\epsilon_{k-p}^{2})},$

with effective gluon energy
$\epsilon_{k}=\sqrt{{\bm{k}}^{2}+m_{G}^{4}/{\bm{k}}^{2}}$. Considering the
fact that a statical potential is calculated at time $t\to\infty$ which
corresponds to $p_{0}\to 0$, $J_{2C}$ disappears in this limit. Taking into
account the Matsubara summations over gluon frequency $k_{0}$,

$\displaystyle T\sum_{k_{0}}{1\over k_{0}^{2}-\epsilon_{k}^{2}}=-{1\over
2\epsilon_{k}}\coth{\epsilon_{k}\over 2T},$ (8) $\displaystyle
T\sum_{k_{0}}{1\over(k_{0}^{2}-\epsilon_{k}^{2})((k_{0}-p_{0})^{2}-\epsilon_{k-p}^{2})}$
$\displaystyle=$
$\displaystyle-\sum_{\eta_{1},\eta_{2}=\pm}{\eta_{1}\eta_{2}\over
4\epsilon_{k}\epsilon_{k-p}}{f_{B}(\eta_{2}\epsilon_{k-p})-f_{B}(\eta_{1}\epsilon_{k})\over
p_{0}-\eta_{1}\epsilon_{k}+\eta_{2}\epsilon_{k-p}},$

where $f_{B}$ is the Bose-Einstein distribution function, $J_{1}$ and $J_{2}$
can be expressed as

$\displaystyle J_{1}(p)$ $\displaystyle=$
$\displaystyle{3p_{i}p_{j}\over{\bm{p}}^{2}}\int{d^{3}{\bm{k}}\over(2\pi)^{3}}{T_{ij}({\bm{k}})\over
2\epsilon_{k}({\bm{k}}-{\bm{p}})^{2}}\coth{\epsilon_{k}\over 2T},$
$\displaystyle J_{2A}(p)$ $\displaystyle=$
$\displaystyle-{1\over{\bm{p}}^{2}}\int{d^{3}{\bm{k}}\over(2\pi)^{3}}{T_{ij}({\bm{k}})T_{ij}({\bm{k}}+{\bm{p}})\over
2\epsilon_{k}}\coth{\epsilon_{k}\over 2T},$ $\displaystyle J_{2B}(p)$
$\displaystyle=$
$\displaystyle-{1\over{\bm{p}}^{2}}\int{d^{3}{\bm{k}}\over(2\pi)^{3}}\sum_{\eta_{1},\eta_{2}=\pm}{\eta_{1}\eta_{2}\epsilon_{k}T_{ij}({\bm{k}})T_{ij}({\bm{k}}-{\bm{p}})\over
2\epsilon_{k-p}}$ (9) $\displaystyle\times$
$\displaystyle{f_{B}(\eta_{2}\epsilon_{k-p})-f_{B}(\eta_{1}\epsilon_{k})\over
p_{0}-\eta_{1}\epsilon_{k}+\eta_{2}\epsilon_{k-p}}.$

We now take analytic extension of $p_{0}\to p_{0}+i\epsilon$ which leads to
the real and imaginary parts of the gluon loop function,

$\displaystyle\text{Re}\Pi_{G}(p)$ $\displaystyle=$ $\displaystyle
J_{1}(p)-J_{2A}(p)-J_{2B}(p),$ $\displaystyle\text{Im}\Pi_{G}(p)$
$\displaystyle=$
$\displaystyle-{\pi\over{\bm{p}}^{2}}\int{d^{3}{\bm{k}}\over(2\pi)^{3}}\sum_{\eta_{1}\eta_{2}=\pm}{\eta_{1}\eta_{2}\epsilon_{k}\over
2\epsilon_{k-p}}$ (10) $\displaystyle\times$
$\displaystyle\left[f_{B}(\eta_{2}\epsilon_{k-p})-f_{B}(\eta_{1}\epsilon_{k})\right](1+\cos^{2}\theta)$
$\displaystyle\times$
$\displaystyle\delta(p_{0}-\eta_{1}\epsilon_{k}+\eta_{2}\epsilon_{k-p}),$

where $\theta$ is the angle between the two momentum vectors ${\bm{k}}$ and
${\bm{k}}-{\bm{p}}$, and the $\delta$ function means the energy conservation
during the decay process from one gluon to two gluons.

The Gribov region $\Omega$ changes only the path integration of the gauge
field, the free quark propagator and in turn the quark loop function are not
affected by the Gribov condition and can be found in textbooks QFT ,

$\displaystyle\text{Re}\Pi_{Q}(p)$ $\displaystyle=$ $\displaystyle{N_{f}\over
2N_{c}{\bm{p}}^{2}}\Bigg{[}{T^{2}\over
6}+2\int{d^{3}{\bm{k}}\over(2\pi)^{3}}\sum_{\eta_{1}\eta_{2}=\pm}$
$\displaystyle\times$
$\displaystyle{\eta_{1}\eta_{2}{\bm{k}}\cdot({\bm{k}}-{\bm{p}})\over|{\bm{k}}||{\bm{k}}-{\bm{p}}|}{f_{F}(\eta_{2}|{\bm{k}}-{\bm{p}}|)-f_{F}(\eta_{1}|{\bm{k}}|)\over
p_{0}-\eta_{1}|{\bm{k}}|+\eta_{2}|{\bm{k}}-{\bm{p}}|}\Bigg{]},$
$\displaystyle\text{Im}\Pi_{Q}(p)$ $\displaystyle=$ $\displaystyle-{\pi
N_{f}\over
N_{c}{\bm{p}}^{2}}\int{d^{3}{\bm{k}}\over(2\pi)^{3}}\sum_{\eta_{1}\eta_{2}=\pm}{\eta_{1}\eta_{2}{\bm{k}}\cdot({\bm{k}}-{\bm{p}})\over|{\bm{k}}||{\bm{k}}-{\bm{p}}|}$
(11) $\displaystyle\times$
$\displaystyle\left[f_{F}(\eta_{2}|{\bm{k}}-{\bm{p}}|)-f_{F}(\eta_{1}|{\bm{k}}|)\right]$
$\displaystyle\times$
$\displaystyle\delta(p_{0}-\eta_{1}|{\bm{k}}|+\eta_{2}|{\bm{k}}-{\bm{p}}|),$

where $f_{F}$ is the Fermi-Dirac distribution function. While the quark loop
$\Pi_{Q}$ is not explicitly affected by the GZ approach, its renormalization
in vacuum is coupled to the Gribov mass $m_{G}$ through the momentum cutoff
$\Lambda_{\overline{MS}}$.

The total loop function contains both the gluon and quark loops,

$\Pi(p)=\Pi_{G}(p)+\Pi_{Q}(p).$ (12)

By summarizing over all gluon and quark loops on a chain QFT , one derives the
loop corrected gluon propagator,

$\displaystyle D_{00}(p)$ $\displaystyle=$
$\displaystyle{1\over{\bm{p}}^{2}}\left[1-\Pi_{00}(p)/{\bm{p}}^{2}+\left(-\Pi_{00}(p)/{\bm{p}}^{2}\right)^{2}+\cdots\right]$
(13) $\displaystyle=$ $\displaystyle{1\over{\bm{p}}^{2}}{1\over
1-g^{2}N_{c}\Pi(p)}.$

We now turn to the calculation of heavy quark potential via the gluon
propagator $D$. Aiming to an in-medium potential, we follow the strategy in
Ref. Laine:2006ns  to construct a real-time Wilson loop which characterizes
the propagation of two infinitely heavy quarks. The evolution of the Wilson
loop satisfies the Schrödinger equation where the potential to the first order
of $g^{2}$ reads Laine:2006ns

$\displaystyle V_{>}(t,r)$ $\displaystyle=$ $\displaystyle
g^{2}C_{F}\int{d^{3}{\bm{p}}\over(2\pi)^{3}}{2-e^{ip_{3}r}-e^{-ip_{3}r}\over
2}$ (14) $\displaystyle\times$
$\displaystyle\Bigg{\\{}{1\over{\bm{p}}^{2}\left(1-g^{2}N_{c}\Pi(0,{\bm{p}})\right)}+\int{dp_{0}\over\pi}f_{B}(p_{0})p_{0}$
$\displaystyle\times$
$\displaystyle\left(e^{p_{0}/T}e^{-ip_{0}t}-e^{ip_{0}t}\right)\Bigg{[}\left({1\over{\bm{p}}^{2}}-{1\over
p_{0}^{2}}\right)\rho_{E}(p)$ $\displaystyle-$
$\displaystyle\left({1\over{\bm{p}}^{2}}-{1\over
p_{3}^{2}}\right)\rho_{T}(p)\Bigg{]}\Bigg{\\}},$

with the constant $C_{F}=(N_{c}^{2}-1)/(2N_{c})$, where $\rho_{E}(p)$ and
$\rho_{T}(p)$ are the two spectral functions Laine:2006ns . Using the relation
$(e^{ip_{0}t}-e^{-ip_{0}t})/p_{0}=2\pi i\delta(p_{0})$ and the approximation
$f(p_{0})\approx T/p_{0}$ for thermalized gluons, we obtain the statical
potential in the limit of $t\to\infty$,

$\displaystyle V(r)$ $\displaystyle=$
$\displaystyle\lim_{t\to\infty}V_{>}(t,r)=V_{R}(r)+iV_{I}(r)$ (15)
$\displaystyle V_{R}(r)$ $\displaystyle=$ $\displaystyle-{C_{F}\over
N_{c}}\int{d^{3}{\bm{p}}\over(2\pi)^{3}}{1\over{\bm{p}}^{2}\text{Re}\Pi(0,{\bm{p}})}\left(1-e^{ip_{3}r}\right),$
$\displaystyle V_{I}(r)$ $\displaystyle=$ $\displaystyle{C_{F}\over
N_{c}}T\int{d^{3}{\bm{p}}\over(2\pi)^{3}}{{\bm{p}}^{2}F(0,{\bm{p}})\over|{\bm{p}}|({\bm{p}}^{2}\text{Re}\Pi(0,{\bm{p}}))^{2}}\left(1-e^{ip_{3}r}\right)$

with

$F(p)=2{|{\bm{p}}|\over p_{0}}\text{Im}\Pi(p).$ (16)

In vacuum, there is no Landau damping, the imaginary parts of the loop
function, propagator and potential disappear automatically, and the potential
is reduced to

$V(r)=C_{F}g^{2}\int{d^{3}{\bm{p}}\over(2\pi)^{3}}D_{00}(0,{\bm{p}})\left(1-e^{i{\bm{p}}\cdot{\bm{r}}}\right).$
(17)

Taking the Gribov mass $m_{G}=0.55$ GeV, which is approximately the value used
in Ref. Golterman:2012dx , the one-loop corrected potential can reproduce very
well the Cornell potential $V(r)=-\alpha/r+\sigma r$ with $\alpha=0.4105$ and
$\sigma=0.2~{}\text{GeV}^{2}$ Zhao:2020jqu , see the comparison in Fig. 3.

Figure 3: The vacuum potential $V(r)$ scaled by the condition $V=0$ at $r=1$
fm. The dashed and solid lines are calculated via Cornell potential and
Gribov-Zwanziger approach.  Figure 4: The real potential $V_{R}(r)$ at finite
temperature, again scaled by the condition $V_{R}=0$ at $r=$ 1 fm.  Figure 5:
The imaginary potential $V_{I}(r)$ at finite temperature.  Figure 6: The
comparison of the real potentials $V_{R}(r)$ calculated via Gribov-Zwanziger
approach (solid lines) and HTL approach (dashed lines).  Figure 7: The
comparison of the imaginary potentials $V_{I}(r)$ calculated via Gribov-
Zwanziger approach (solid lines) and HTL approach (dashed lines).

The heavy quark potential in medium contains real and imaginary parts, shown
in Figs. 4 and 5 at several temperatures. The real part $V_{R}(r)$ is, as
usually discussed in literatures Satz:2005hx , controlled by color screening.
In vacuum, the potential increases linearly with the distance and is never
saturated, which means the parton confinement. In hot medium, the potential
becomes saturated at a finite distance (Debye screening length $r_{D}$), due
to the color screening. In the saturated region, the interaction force ${\bf
F}=-{\bm{\nabla}}V=0$ indicates the vanishing color interaction at the
distance $r>r_{D}$. With increasing temperature, the screening length
decreases from infinity in vacuum to about 0.8, 0.6 and 0.4 fm at $T=0.15$,
$0.25$ and $0.35$ GeV. Considering that $r_{D}=0.4$ fm is already less than
the $J/\psi$ radius $r_{J/\psi}\sim 0.5$ fm, charmonia are expected to be
dissociated in a fireball with temperature $0.25<T<0.35$ GeV.

For a dynamically evaluating system created in high energy nuclear collisions,
quarkonia are dissociated by not only the color screening, but also the Landau
damping at finite temperature which leads to a decay width characterized by
the imaginary part $V_{I}(r)$ in potential models. Our calculation with GZ
approach shows a strong imaginary part of the potential. Especially at the
critical distance of charmonium radius $r\sim 0.5$ fm, the strengths of the
real and imaginary parts are at the same order, see Figs. 4 and 5.

It is clear that the one-loop corrected potential should approach to the HTL
result in high temperature limit Laine:2006ns , when the Gribov mass
$m_{G}(T)$ goes to zero,

$\displaystyle V(r)$ $\displaystyle=$ $\displaystyle-{C_{F}g^{2}\over
4\pi}\left[\left(m_{D}+{e^{-m_{D}r}\over r}\right)+iT\phi(m_{D}r)\right],$
$\displaystyle\phi(x)$ $\displaystyle=$ $\displaystyle
2\int_{0}^{\infty}dz{z\over(z^{2}+1)^{2}}\left(1-{\sin(zx)\over zx}\right),$
(18)

where $m_{D}=gT\sqrt{N_{c}/3+N_{f}/6}$ is the Debye mass in HTL approximation.
From the comparison between the two real potentials (15) with GZ approach and
(Heavy-quark potential in Gribov-Zwanziger approach around deconfinement phase
transition) with HTL approach, shown in Fig. 6, while the two calculations are
close to each other at high temperature, indicating the almost same color
screening strengths in the two approaches, the HTL approach fails to go back
to a confinement potential at low temperature, according to its definition. As
for the imaginary part, shown in Fig.7, the strength of the one with GZ
approach is much larger than the one with HTL approach. The physics is the
Gribov mass which behaviors like a magnetic mass $m_{G}\sim g^{2}T$ at high
temperature. The still large difference between the two imaginary parts at
temperature $T=0.35$ GeV $\sim 2T_{c}$ comes from two reasons: the temperature
is not high enough to satisfy the HTL condition $T\gg|{\bm{p}}|$, and the
Gribov mass is still large $m_{G}\sim 0.3$ GeV. The slowly dropping down of
the Gribov mass at temperature $T\gtrsim T_{c}$ is probably due to the gap
equation at the lowest level. The loop correction to the gap equation should
be considered in the future.

Confinement in vacuum and deconfinement in hot medium are non-perturbative
problems in QCD, but they cab be described by the gluon propagator at lowest
level in the Gribov-Zwanziger approach. This leads to an alternative way to
perturbatively calculate the heavy quark potential through the gluon
propagator. We calculated in this paper the loop-corrected heavy quark
potential via the resummation of gluon and quark loops in the strong coupling
region around the deconfinement phase transition temperature $T_{c}$. The loop
correction makes the potential complex, the saturation of the real part is
controlled by the color screening, and the imaginary part is characterized by
the Landau damping in hot medium. In comparison with the HTL approach of QCD,
the high temperature behavior of the real part is almost the same in the two
approaches, while the Gribov mass results in a much stronger imaginary part
due to effectively including the magnetic interaction in the GZ approach. For
the distance around the typical quarkonium size, the strength of the imaginary
part is comparable with the real part. This means that, different from the
previous quarkonium dissociation picture by color screening, Landau damping
plays also an important role.

Acknowledgement: The work is supported by the Guangdong Major Project of Basic
and Applied Basic Research No. 2020B0301030008 and the NSFC grants Nos.
11890712 and 12075129.

## References

  * (1) W. E. Caswell and G. P. Lepage, Phys. Lett. B 167, 437-442 (1986).
  * (2) N. Brambilla, A. Pineda, J. Soto and A. Vairo, Nucl. Phys. B 566, 275 (2000).
  * (3) N. Brambilla, A. Pineda, J. Soto and A. Vairo, Rev. Mod. Phys. 77, 1423 (2005).
  * (4) See for instance the recent review paper, J.Zhao, K.Zhou, S.Chen and P.Zhuang, Prog. Part. Nucl. Phys. 114, 103801(2020).
  * (5) Y. Koma, M. Koma and H. Wittig, Phys. Rev. Lett. 97, 122003 (2006).
  * (6) T. Kawanai and S. Sasaki, Phys. Rev. D 85, 091503 (2012).
  * (7) T. Matsui and H. Satz, Phys. Lett. B 178, 416-422 (1986).
  * (8) X. Dong and V. Greco, Prog. Part. Nucl. Phys. 104, 97-141 (2019).
  * (9) A. Rothkopf, Phys. Rept. 858, 1-117 (2020).
  * (10) E. Chapon, D. d’Enterria, B. Ducloue, M. G. Echevarria, P. B. Gossiaux, V. Kartvelishvili, T. Kasemets, J. P. Lansberg, R. McNulty and D. D. Price, et al. Prog. Part. Nucl. Phys. 122, 103906 (2022).
  * (11) J. Zhao, S. Shi and P. Zhuang, Phys. Rev. D 102, no.11, 114001 (2020).
  * (12) J. Zhao, H. He and P. Zhuang, Phys. Lett. B 771, 349-353 (2017).
  * (13) H. He, Y. Liu and P. Zhuang, Phys. Lett. B 746, 59-63 (2015).
  * (14) S. Cho et al. [ExHIC], Prog. Part. Nucl. Phys. 95, 279-322 (2017).
  * (15) S. Nadkarni, Phys. Rev. D 34, 3904 (1986).
  * (16) C. Y. Wong, Phys. Rev. C 72, 034906 (2005).
  * (17) H. Satz, J. Phys. G 32, R25 (2006).
  * (18) O. Kaczmarek, F. Karsch, F. Zantow and P. Petreczky, Phys. Rev. D 70, 074505 (2004) [erratum: Phys. Rev. D 72, 059903 (2005)].
  * (19) N. Brambilla, J. Ghiglieri, A. Vairo and P. Petreczky, Phys. Rev. D 78, 014017 (2008).
  * (20) M. Laine, O. Philipsen, P. Romatschke and M. Tassler, JHEP 03, 054 (2007).
  * (21) A. Beraudo, J. P. Blaizot and C. Ratti, Nucl. Phys. A 806, 312-338 (2008).
  * (22) A. Bazavov, Y. Burnier and P. Petreczky, Nucl. Phys. A 932, 117-121 (2014).
  * (23) A. Rothkopf, T. Hatsuda and S. Sasaki, Phys. Rev. Lett. 108, 162001 (2012).
  * (24) Y. Burnier, O. Kaczmarek and A. Rothkopf, Phys. Rev. Lett. 114, no.8, 082001 (2015).
  * (25) Y. Burnier, O. Kaczmarek and A. Rothkopf, JHEP 12, 101 (2015).
  * (26) D. Lafferty and A. Rothkopf, Phys. Rev. D 101, no.5, 056010 (2020).
  * (27) D. Bala and S. Datta, Phys. Rev. D 101, no.3, 034507 (2020).
  * (28) S. Shi, K. Zhou, J. Zhao, S. Mukherjee and P. Zhuang, Phys. Rev. D 105, no.1, 1 (2022).
  * (29) R. Larsen, S. Meinel, S. Mukherjee and P. Petreczky, Phys. Lett. B 800, 135119 (2020).
  * (30) D. Bala et al. [HotQCD], Phys. Rev. D 105, no.5, 054513 (2022).
  * (31) L. D. Faddeev and V. N. Popov, Phys. Lett. B 25, 29-30 (1967).
  * (32) V. N. Gribov, Nucl. Phys. B 139, 1 (1978).
  * (33) Y.Dokshitzer and D.Kharzeev, Ann. Rev. Nucl. Part. Sci. 54, 487(2004).
  * (34) D. Zwanziger, Nucl. Phys. B 209, 336-348 (1982).
  * (35) D. Zwanziger, Nucl. Phys. B 323, 513-544 (1989).
  * (36) D. Zwanziger, Nucl. Phys. B 399, 477-513 (1993).
  * (37) D. Zwanziger, Phys. Rev. D 76, 125014 (2007).
  * (38) A. Cucchieri and D. Zwanziger, Phys. Rev. D 65, 014002 (2001).
  * (39) J. A. Gracey, JHEP 02, 009 (2010).
  * (40) F. E. Canfora, D. Dudal, I. F. Justo, P. Pais, L. Rosa and D. Vercauteren, Eur. Phys. J. C 75, no.7, 326 (2015).
  * (41) G. Burgio, M. Quandt and H. Reinhardt, Phys. Rev. Lett. 102, 032002 (2009).
  * (42) M. Golterman, J. Greensite, S. Peris and A. P. Szczepaniak, Phys. Rev. D 85, 085016 (2012).
  * (43) K. Fukushima and N. Su, Phys. Rev. D 88, 076008 (2013).
  * (44) J. I. Kapusta and C. Gale, Finite-temperature field theory: Principles and applications (Cambridge University Press, Cambridge, 2009). M. Le Bellac, Thermal Field Theory (Cambridge University Press, Cambridge, 1996).