# Revisiting the Auction Algorithm for Weighted Bipartite Perfect Matchings

Megha Khosla and Avishek Anand

(L3S Research Center, Hannover, Germany
<EMAIL_ADDRESS>
###### Abstract

We study the classical _weighted perfect matchings problem_ for bipartite
graphs or sometimes referred to as the _assignment problem_ , i.e., given a
weighted bipartite graph $G=(U\cup V,E)$ with weights
$w:E\rightarrow\mathcal{R}$ we are interested to find the maximum matching in
$G$ with the minimum/maximum weight. In this work we present a new and
arguably _simpler_ analysis of one of the earliest techniques developed for
solving the assignment problem, namely the _auction algorithm_. Using our
analysis technique we present tighter and improved bounds on the runtime
complexity for finding an approximate minumum weight perfect matching in
$k$-left regular sparse bipartite graphs.

##  1 Introduction

Let $G=(U\cup V,E)$ be a bipartite graph on $n$ vertices and $m$ edges and let
$w:E\rightarrow\mathcal{R}$ be a weight function on edges. A perfect matching
of $G$ is a subset $M\subseteq E$ of the edges such that for every node $v\in
U\cup V$ there is exactly one incident edge $e\in M$. The weight of a matching
$M$ is given by the sum of the weights of its edges, i.e. $w(M):=\sum_{e\in
M}w(e)$. The minimum weight bipartite perfect matching problem (MWPM) is to
find for a given bipartite graph $G$ and a given weight function $w$ a perfect
matching of minimum weight which is sometimes also referred to as the
_assignment problem_. WLOG we can assume that $|U|\leq|V|=n$.

In this work we present a novel and simple analysis of _auction_ algorithm
originally proposed by Bertsekas [7] for solving the assignment problem. The
auction algorithm resembles a competitive bidding process whereby unassigned
persons (nodes on the left set $U$) bid simultaneously for objects (nodes on
the right set $V$), thereby raising their prices. On obtaining all bids, the
object $v\in V$ is assigned to the highest bidder.

There is a cost $w(u,v)$ for matching person $u$ with object $v$ and we want
to assign persons to objects so as to minimize the overall cost (for a
minimization objective). Let $E$ be the set of all pairs $(u,v)$ that can be
matched. A typical iteration of the auction algorithm consists of a bidding
and an assignment phase. To start off, each object $v$ is initialized with
some initial price $L(v)$. The bidding and assignment phases are as follows:

Bidding Phase : Let $I$ be the set of unassigned persons. Each person $u\in I$
finds an object $\bar{v}_{u}$ which optimizes for minimum cost, $w(u,v)+L(v)$
, that is,

$\bar{v}_{u}\,=\,\arg\min_{v:(u,v)\in E}\,w(u,v)+L(v)$

and computes a _bidding increment_ $\gamma_{u}$ for some parameter
$\varepsilon>0$

$\gamma_{u}=w^{2}_{u}-w^{1}_{u}+\varepsilon,$

where

$w^{1}_{u}=\min_{v:(u,v)\,\in\,E}w(u,v)+L(v)$

is the best object value and

$w^{2}_{u}=\min_{v:(u,v)\,\in\,E,\,v\neq\bar{v}_{u}}w(u,v)+L(v)$

is the second best object value.

Assignment Phase : Note that an object $v$ can be the best bid for multiple
persons $u\in I$. In such a case $v$ is assigned to the person with the
_highest bid_ and its price is raised by the highest bid, i.e.,

$L(v)=L(v)+\max_{u\in I:v=\bar{v}_{u}}\gamma_{u}.$

The person that was assigned to $v$ at the beginning of the iteration (if any)
becomes unassigned. The algorithm continues with a sequence of iterations
until all persons have an assigned object.

The auction algorithm finds an approximate solution for maximum/minimum weight
perfect matching which is bounded by $OPT+n\varepsilon$, which we also refer
as the _$\varepsilon$ \- optimal_ solution. Let $w^{max}$ and $w^{min}$ denote
the maximum and minimum edge weights respectively. The worst case running time
for complete bipartite graphs is $O({n^{2}w^{max}\over\varepsilon})$. In order
to improve the running time a procedure called _$\varepsilon$ -scaling_ is
employed. It essentially consists of executing the algorithm several times
starting with a large value of $\varepsilon$ and successively reducing
$\varepsilon$ after each run. The procedure terminates typically when the
ultimate value of $\varepsilon$ is less than some critical value (for example,
$1/n$, for integral weights). For integral weights, it has be shown that the
worst-case running time of the auction algorithm for finding optimal solutions
using $\varepsilon$-scaling is $O(nm\log(nw^{max}))$ [5, 6]. However for the
asymmetric problem where number of persons is less than the number of objects,
the prices would need to be initialized by $0$ and $\varepsilon$-scaling
cannot be used out of the box.

In order to benefit from $\varepsilon$-scaling an _reverse_ auction method is
used (in addition to the regular auction) in which the objects also compete
for persons by offering discounts [4].

Subsequently the auction method was extended to solve the classical linear
network flow problem and many of its special classes. In particular, [10] and
[2] propose an extension to the minimum cost network flow problem using
$\varepsilon$-relaxation. Auction algorithms for transportation problems [3]
and shortest paths [8] have also been proposed. We refer the interested reader
to [9] for a more comprehensive discussion on auction algorithms.

In this paper we focus on the sequential version of auction algorithm where a
single unassigned person bids at every iteration. This version is also known
as the _Gauss-Seidel_ version because of its similarity with Gauss-Seidel
methods for solving systems of nonlinear equations. Moreover, we restrict
ourselves to its application to finding minimum weight perfect matchings in
bipartite graphs. Our analysis technique is inspired by the LSA method [23]
which is used to construct large hash tables and finding maximum matchings in
unweighted bipartite graphs. In fact, the original motivation was to extend
the label based technique in LSA for solving the weighted version of the
matching problem in bipartite graphs. Though our proposed algorithm later
turned out to be a version of the auction algorithm, our analysis technique
allows for a simpler interpretation of prices (labels) as a function of
_shortest unweighted paths_ in the underlying bipartite graphs. This in turn
helps us to bound the runtime of the weighted version in terms of the shortest
unweighted paths in the underlying graphs (which also provides a bound on the
runtime of LSA) with weight range ($w^{max}-w^{min}$) in the multiplicative
factor.

In particular, we use the main result in [23] to show that the worst case
runtime of auction algorithm for sparse $k$\- regular bipartite graphs is no
more than $O(n\cdot{w^{max}-w^{min}\over\varepsilon})$ with high probability.
For complete bipartite graphs, we prove a slightly better runtime bound of
$O(n^{2}\cdot{w^{max}-w^{min}\over\varepsilon})$, which is an improvement over
the previous bound when both $w^{max}$ and $w^{min}$ are large.

###  1.1 More on Related Work

The first polynomial time algorithm for the assignment problem, the so called
Hungarian method, was given by [25, 26]; implementations of this algorithm had
a running time of $O(n^{3})$, which is still optimal for dense graphs. Edmonds
and Karp [14] and Tomizawa [33] independently observed that the assignment
problem is reducible to computing single-source shortest paths on non-negative
weights. Further, Fredman and Tarjan [16] showed that using Fibonacci heaps,
$n$ executions of Dijkstra’s [11] shortest path algorithm take $O(mn+n^{2}\log
n)$ time. On integer weighted graphs this algorithm can be implemented
slightly faster, in $O(mn+n^{2}\log\log n)$ time [20, 32] or $O(mn)$ time
(randomized) [1, 31], independent of the maximum edge weight. Gabow and Tarjan
[18] gave a scaling algorithm for the assignment problem running in
$O(m\sqrt{n}\log(nw^{max}))$ time for integral weights that is based on the
Hungarian method. Orlin and Ahuja [28] using the auction approach, and
Goldberg and Kennedy [19], obtain the same time bound. Recently the new
scaling algorithm by Duan et al. [12] matches the same bound for weighted
matchings in general graphs. Sankowski takes an algebraic approach by
proposing a randomized algorithm [30], that solves the assignment problem
using fast matrix multiplications in $O(w^{max}n^{\omega})$ time with high
probability, where $\omega$ is the exponent of square matrix multiplication.
Improved runtime bounds have also been achieved by reducing maximum weighted
matching to maximum cardinality matching problems ( see [22, 21, 29] and
references therein)

The other area of related work concerns with the algorithms for maximum
weighted matchings (MWM). It is important to note that though the two problems
MWM and MWPM are reducible to each other, the reductions do not work for their
approximate versions as the approximation may compromise perfection. Duan and
Pettie [13] presented a $(1-\varepsilon)$-MWM algorithm that runs in
$O(m\varepsilon^{-1}\log(\varepsilon^{-1}))$ on general graphs and
$O(m\varepsilon^{-1}\cdot\min\\{\log\varepsilon^{-1},\log(w^{max})\\})$ time
on integer-weighted general graphs. These results cannot be compared with the
bound that we provide for random sparse bipartite graphs as we require that
the matching returned should also be perfect. Please refer to [13] for more
details on other works for finding approximate maximum weight matching.

###  1.2 Our Contribution

We present a novel analysis of auction algorithm for finding approximate
minimum weight perfect matchings in bipartite graphs for arbitrary weights.
Our approach is inspired by the _local search allocation_ (LSA) method [23]
used to construct large hash tables and finding maximum matchings in
unweighted bipartite graphs. Using our analysis we could easily provide a
better runtime bound for a random sparse $k-$ left regular bipartite graphs.
In a random $k$-left regular bipartite graph, each vertex in the left-set
chooses $k$ neighbors from the right set independently and uniform at random.

From [27, 15, 17] we know that for such graphs, for all $k\geq 3$, there
exists a threshold density $c^{*}_{k}$ ( and is computable , for example for
$k=3$ it is close $0.91$) such that when $|U|/|V|<c^{*}_{k}$, there exists a
left perfect matching in $G$ with high probability, otherwise this is not the
case. We will show that for random $k$\- left regular bipartite graphs obeying
this threshold condition, one can find an $\varepsilon$\- optimal minimum
weight perfect matching in near linear time. In case a perfect matching does
not exist, we will show that with appropriate stopping criteria, the algorithm
will stop and return an approximate solution for minimum weight maximum
matching (though we do not analyse the approximation guarantee in detail). We
note that there exists no restriction on the weights and they can be
arbitrary. For complete bipartite graphs the runtime bound is
$O(n^{2}\cdot{w^{max}-w^{min}\over\varepsilon})$. We now state the main
result.

######  Theorem 1.

Let $G=(U\cup V,E)$ be a weighted bipartite graph with $|U|\leq|V|=n$ and
arbitrary weights. Let $OPT_{G}$ be the weight of the optimal minimum weight
maximum matching in $G$. Then for an arbitrary parameter $0<\varepsilon$ the
auction algorithm returns a $OPT_{G}+n\varepsilon$ minimum weight maximum
matching. The worst case runtime for sparse and complete bipartite graphs is
as follows.

  1. 1.

Sparse Graphs: For $k\geq 3$, let $G=(U\cup V,E)$ be such that each vertex in
$U$ chooses $k$ neighbors in $V$ independently and uniform at random. In
addition let $|U|<c^{*}_{k}n$. The worst case runtime in case a perfect
matching exists is $O(n\cdot{w^{max}-w^{min}\over\varepsilon})$ with
probability $1-o(1)$.

  2. 2.

Complete Graphs: For complete bipartite graphs the runtime bound is
$O(n^{2}\cdot{w^{max}-w^{min}\over\varepsilon})$.

The main idea behind the proof for the runtime bounds is to show that the
prices (which we will refer to as labels in the subsequent section) when
initialized with zero are increasing by at least $\varepsilon$ and are bounded
in terms of shortest distances in the underlying graph and the weight range.
We will then use the result from [23] which provides a linear bound in
expectation and with high probability for the sum of the shortest distances in
the special class of sparse bipartite graphs as considered in Theorem 1. For
the approximation guarantee, we will show that no alternating path allows for
a decrease in the existing matching weight by more than $n\varepsilon$. Note
that we do not provide any improved bound on approximation quality of the
solution in this paper. We start by giving a detailed description of the
algorithm in the next section.

##  2 The Auction Algorithm for Perfect Weighted Matchings

###  2.1 Notations.

Throughout the paper we use the following notations unless stated otherwise.
We denote the set of integers $\\{1,2,\ldots,n\\}$ by $[n]$. Let $G=(U\cup
V;E)$ denote a weighted bipartite graph where $|U|\leq|V|=n$. For any $e\in
E$, $w(e)$ denotes the weight on edge $e$. In addition let $w^{max}$ and
$w^{min}$ denote the maximum and minimum weight on any edge respectively. For
any vertex $u\in U$ we refer to set of neighbors of $u$ in $V$ as $N(u)$. For
any $v\in V$, $L(v)$ denotes the label of $v$. Let $M$ be the set of matched
edges in the optimal matching. For $u\in U$ and $v\in V$, we say that an edge
$e=(u,v)$ is assigned to $v$ if $e\in M$. We call a vertex $v\in V$ _free_ if
and only if no edge is assigned to it. The optimal minimum weight for a
perfect matching in $G$ is then equal to $\sum_{e\in M}w(e)$ and is denoted by
$OPT_{G}$.

###  2.2 The Algorithm and its Analysis

We first explain the auction algorithm in detail with respect to our
interpretation. Initially we are provided with all vertices from the right set
$V$. We assign labels to all vertices in $V$ and initialize them by $0$. The
vertices of the left set $U$ appear one by one together with the incident
edges and their weights. An edge $e=(u,v)$ is assigned to vertex $v\in V$ such
that the sum $L(v)+w(e)$ is minimum for all edges incident on $u$. We refer to
the rule for choosing a candidate edge as the _choice rule_. The original
algorithm aimed only to find prefect matchings. In case a perfect matching
does not exist, the algorithm might enter in an endless loop. We will argue in
the next section that the label values are bounded and therefore resolve this
issue by checking the minimum label with the maximum possible value of
$n/2(w^{max}-w^{min}+\varepsilon)$. In case for some vertex $u$ all its
neighbors have labels greater than the maximum value, $u$ is discarded and
never matched.

Let for vertex $v^{\prime}\in V\backslash\\{v\\}$ and
$e^{\prime}=(u,v^{\prime})$, the sum $L(v^{\prime})+w(e^{\prime})$ is the
minimum. For some $0<\varepsilon$, the label of $v$ after its new assignment
will be updated as follows.

$L(v)=L(v^{\prime})+w(e^{\prime})-w(e)+\varepsilon.$

We call the above rule as the _update rule_.

In case $v$ is not empty and was already assigned another edge $(u,v)$, the
edge $(u,v)$ is moved out of the matching and the process is repeated for
vertex $u$.

Let $\mathbf{L}=\\{L(v_{1}),\ldots,L(v_{n})\\}$ and
$\mathbf{T}=\\{T(v_{1}),\ldots,T(v_{n})\\}$ where $L(v_{i})$ denotes the label
of vertex $v_{i}$ and $T(v_{i})$ denotes the vertex matched to $v_{i}$. We
initialize $\mathbf{L}$ with all $0$s , i.e., all vertices are free. Auction
algorithm is described in Algorithm 1 which calls the Procedure 2 to match an
arbitrary vertex from the left set when it appears.

Algorithm 1 Auction algorithm ($U,V,E$)

1: for all $v\in V$ do

2: Set $L(v)=0$

3: Set $T(v)=\emptyset$

4: for all $u\in U$ do

5: CALL MatchVertex ($u,\mathbf{L},\mathbf{T}$)

Procedure 2 MatchVertex ($u,\mathbf{L},\mathbf{T}$)

1: Choose $v\in N(u)$ such that $L(v)+w(u,v)$ is the minimum $\rhd$Choice Rule

2: if $L(v)>n/2(w^{max}-w^{min}+\varepsilon)$ then

3: RETURN

4: $L(v)\leftarrow\min{(L(v^{\prime})+w(u,v^{\prime})|v^{\prime}\in
N(u)\setminus\\{v\\})}-w(u,v)+\varepsilon$ $\rhd$Update Rule

5: if $(T(v)\neq\emptyset)$ then

6: $y\leftarrow T(v)$ $\rhd$Move that moves an edge out of matching

7: $T(v)\leftarrow u$ $\rhd$Move that assigns a new edge or matches a new
vertex

8: $\mathbf{CALL}$ MatchVertex($y,\mathbf{L},\mathbf{T}$)

9: else

10: $T(v)\leftarrow u$ $\rhd$Move that assigns an edge or matches a vertex

We observe that if Auction algorithm does not enter an endless loop in the
Procedure 2, it would return a matched vertex for each of the vertices in $U$.
This implies that Auction algorithm will find a left perfect matching if it
terminates.

In the next section we will prove that the algorithm does not run in endless
loops, i.e., it terminates, by showing that the (1) labels are non decreasing
and (2) labels are bounded by a maximum value.

###  2.3 Bounding the Labels

We need some additional notation. In what follows a _move_ denotes either
assigning an edge to a free vertex or replacing a previously assigned edge.
Let $P$ be the total number of moves performed by the algorithm. For $p\in[P]$
we use $L_{p}(v)$ to denote the label of vertex $v$ at the end of the $p$th
move. Let $\mathcal{M}_{p}$ denote the set of matched edges at the end of the
$p$th move.

Let an edge $(u,v)$ is assigned to a vertex $v\in V$ in some move $p$. The
choice and the update rules can then be rewritten as follows.

$\displaystyle\textbf{Choice Rule : }L_{p-1}(v)+w(u,v)\leq\min_{v^{\prime}\in
N(u)\setminus\\{v\\}}(L_{p-1}(v^{\prime})+w(u,v^{\prime})).$ (1)
$\displaystyle\textbf{Update Rule : }L_{p}(v)=\min_{v^{\prime}\in
N(u)\setminus\\{v\\}}(L_{p-1}(v^{\prime})+w(u,v^{\prime}))-w(u,v)+\varepsilon.$
(2)

We will first show that the label of exactly one vertex increases by at least
$\varepsilon$ at the end of a move. The labels of other vertices remain
unchanged. Even though this fact is obvious from the algorithm description, we
prove it here for completeness.

######  Proposition 1.

For any $p\in[P]$ there exists exactly one vertex $v\in V$ such that
$L_{p}(v)\geq L_{p-1}(v)+\varepsilon$ and for all other $v^{\prime}\in
V\backslash\\{v\\}$, $L_{p}(v^{\prime})=L_{p-1}(v^{\prime})$.

###### Proof.

By definition of the move, the label of exactly one vertex is altered at the
end of a move. Let $v$ be assigned as edge $(u,v)$ is some move $p\in[P]$.
Therefore labels of all vertices except $v$ remain unchanged at the end of the
$p$th move, i.e.,

$\forall v^{\prime}\in
V\backslash\\{v\\}:L_{p}(v^{\prime})=L_{p-1}(v^{\prime})$

For vertex $v$ the new label at the end of the $p$th move is defined by (2)
(the update rule) as

$L_{p}(v)=\min_{v^{\prime}\in
N(u)\backslash\\{v\\}}(L_{p-1}(v^{\prime})+w(u,v^{\prime}))-w(u,v)+\varepsilon,$

which combined with (1) gives $L_{p}(v)\geq L_{p-1}(v)+\varepsilon$, thereby
concluding the proof.

∎

In the following proposition we bound the label of any vertex with respect to
the labels of the other neighbors of its matched vertex and the corresponding
edge weights. We will need this to bound the maximum label of any vertex at
the end of the $(P-1)$th move.

######  Proposition 2.

For all $p\in[P]$ and all $(u,v)\in\mathcal{M}_{p}$, the following holds.

$L_{p}(v)\leq\min_{v^{\prime}\in
N(u)\setminus\\{v\\}}(L_{p}(v^{\prime})+w(u,v^{\prime}))-w(u,v)+\varepsilon.$

###### Proof.

Let in some move $p\in P$ an edge $(u,v)$ is placed in the matched set
$\mathcal{M}_{p}$. Note that the labels of all vertices $v^{\prime}\in
V\setminus\\{v\\}$ remain unchanged at the end of the $p$th move. By update
rule we obtain

$\displaystyle L_{p+1}(v)=$ $\displaystyle\min_{v^{\prime}\in
N(u)\backslash\\{v\\}}(L_{p}(v^{\prime})+w(u,v^{\prime}))-w(u,v)+\varepsilon$
$\displaystyle\leq$ $\displaystyle\min_{v^{\prime}\in
N(u)\backslash\\{v\\}}(L_{p+1}(v^{\prime})+w(u,v^{\prime}))-w(u,v)+\varepsilon.$

The last inequality holds $L_{p+1}(v^{\prime})\geq L_{p}(v^{\prime})$ for all
$v^{\prime}$ and $p$ by Proposition 1. For any other edge
$(u^{\prime},v^{\prime})\in\mathcal{M}_{p}$ that was last assigned in some
move $p^{\prime}<p$ the following holds by update rule

$\displaystyle L_{p^{\prime}}(v^{\prime})=\min_{v^{\prime\prime}\in
N(u^{\prime})\backslash\\{v^{\prime}\\}}(L_{p^{\prime}-1}(v^{\prime\prime})+w(u^{\prime},v^{\prime\prime}))-w(u^{\prime},v^{\prime})+\varepsilon$
$\displaystyle\leq\min_{v^{\prime\prime}\in
N(u^{\prime})\backslash\\{v^{\prime}\\}}(L_{p}(v^{\prime\prime})+w(u^{\prime},v^{\prime\prime}))-w(u^{\prime},v^{\prime})+\varepsilon$

The last inequality holds as $L_{p}(v^{\prime\prime})\geq
L_{p^{\prime}-1}(v^{\prime\prime})$ for all $v^{\prime\prime}$ and all
$p>p^{\prime}-1$. Also as $v^{\prime}$ was not updated after $p^{\prime}$th
move, $L_{p^{\prime}}=L_{p}$. We therefore obtain

$L_{p}(v^{\prime})\leq\min_{v^{\prime\prime}\in
N(u^{\prime})\backslash\\{v^{\prime}\\}}(L_{p}(v^{\prime\prime})+w(u^{\prime},v^{\prime\prime}))-w(u^{\prime},v^{\prime})+\varepsilon,$

thereby completing the proof.

∎

####  2.3.1  Maximum Label and the Runtime

We want to prove a bound of the labels of the vertices in terms of their
shortest distances to the set of free vertices at the end of the $(P-1)$th
move. We start by considering an ordered alternating path of unmatched and
matched edges. Using Proposition 2 we will exploit the relationship between
(a)label of the last vertex (b) label of the first vertex and (c) edge weights
in this ordered path.

For an even $4\leq t\leq 2n-4$, let
$B_{t}=(v_{1},u_{2},v_{2},\cdots,u_{{t\over 2}+1},v_{{t\over 2}+1})$ denote an
alternating path of unmatched and matched edges at the end of some move $p$.
We assume that the first edge in the path is an unmatched edge. Let
$\mathcal{M}_{p}(B_{t})$ and $\mathcal{M}^{\prime}_{p}(B_{t})$ denote the set
of matched and unmatched edges in $B_{t}$ at the end of some $p$th move.

######  Lemma 1.

For all $t\leq 2n-4$ and all paths $B_{t}$ as defined above, the following
holds.

$L_{p}(v_{t\over 2})\leq
L_{p}(v_{1})+\sum_{e\in\mathcal{M}^{\prime}_{p}(B_{t})}w(e)-\sum_{e\in\mathcal{M}_{p}(B_{t})}w(e)+{t\varepsilon\over
2}$

###### Proof.

We will prove the lemma by induction on the path length $t$. For the case
$t=2$ we have $B_{2}=(v_{1},u_{2},v_{2})$, where
$(v_{1},u_{2})\in\mathcal{M}^{\prime}_{p}(B_{2})$ and
$(u_{2},v_{2})\in\mathcal{M}_{p}(B_{2})$. By Proposition 2 we obtain

$L_{p}(v_{2})\leq L_{p}(v_{1})+w(v_{1},u_{2})-w(u_{2},v_{2})+\varepsilon.$

Clearly the lemma follows for path length $2$. Now assume that the lemma is
true for path length $t=2t^{\prime}$, i.e.,

$L_{p}(v_{t^{\prime}})\leq
L_{p}(v_{1})+\sum_{e^{\prime}\in\mathcal{M}^{\prime}_{P}(B_{2t^{\prime}})}w(e^{\prime})-\sum_{e\in\mathcal{M}_{p}(B_{2t^{\prime}})}w(e)+{t^{\prime}\varepsilon}$

We will now prove the lemma for paths of length $2t^{\prime}+2$. Consider the
last vertex $v_{{t^{\prime}}+1}$ in $B_{2t^{\prime}+2}$. As it is matched to
$u_{t^{\prime}+1}$, by Proposition 2 we obtain

$L_{p}(v_{t^{\prime}+1})\leq
L_{p}(v_{t^{\prime}})+w(u_{t^{\prime}+1},v_{t^{\prime}})-w(u_{t^{\prime}+1},v_{t^{\prime}+1})+\varepsilon.$

Combining the above inequality with the induction hypothesis we obtain

$L_{p}(v_{t^{\prime}+1})\leq
L_{p}(v_{1})+\sum_{e^{\prime}\in\mathcal{M}^{\prime}_{P}(B_{2t^{\prime}+2})}w(e^{\prime})-\sum_{e\in\mathcal{M}_{p}(B_{2t^{\prime}+2})}w(e)+(t^{\prime}+1)\varepsilon,$

hence completing the proof. ∎

We obtain the following corollary about the labels at the end of $(P-1)th$
move.

######  Corollary 1.

Let $d_{P-1}(v)$ denote the length of the shortest unweighted path to any free
vertex in $G$ after move $P-1$ is completed. Then for all $v\in V$,
$L_{P-1}(v)\leq{d_{P-1}(v)\over 2}(w^{max}-w^{min}+\varepsilon)$.

###### Proof.

We assume that the given bipartite graph is connected, otherwise we run the
algorithm on connected components. We note that at the end of $P-1$th move,
there is at least one free vertex in $V$. Note that the label of any free
vertex is $0$. For any $v\in V$we use Lemma 1 considering its shortest
alternating path to some free vertex and bounding the weights of unmatched and
matched edges on this path by $w^{max}$ and $w^{min}$ respectively, we obtain
the desired result, i.e.,

$L_{P-1}(v)\leq 0+{d_{P-1}(v)\over 2}(w^{max}-w^{min}+\varepsilon)$

∎

Before we prove the runtime bounds we will describe the main result from [23]
which we use to bound the distance values $d_{P-1}(v)$. Khosla [23] considers
the problem of assigning $m$ items to $n$ locations, such that each of the $m$
items chooses $k$ locations independently and uniformly at random. Each item
needs to be assigned to one of its $k$ choices. It is easy to see that such an
assignment instance represents a $k$-left regular bipartite graph and a valid
assignment corresponds to a left-perfect matching (where all vertices from the
left set have been matched).

A label based approach, LSA is used to find such a perfect matching for the
case $m<c^{*}_{k}n$, where $c^{*}_{k}$ is the _threshold density_ (known from
[27, 15, 17] before). The fact that $m<c^{*}_{k}n$ ensures that a left perfect
matching exists with probability $1-o(1)$. The runtime of the algorithm is
bounded by the sum of labels at the end of the algorithm which are in turn
bounded by the shortest distances to the set of free vertices. Note that
$c^{*}_{k}<1$ and there always exist at least one free vertex at the end of
the algorithm. The distances are in turn bounded using some structural
properties of the corresponding graphs. She shows that the sum of shortest
distances to the set of free vertices is bounded by $O(n)$ with probability
$1-o(1)$ . In a subsequent extension Anand and Khosla [24] show that the
result also holds in expectation. We state here their main result adjusted to
the terminology and its requirement in this paper.

######  Theorem 2.

For $k\geq 3$, let $G=(U\cup V,E)$ be such that each vertex in $U$ chooses $k$
neighbors in $V$ independently and uniformly at random. In addition let
$|U|<c^{*}_{k}|V|$. Then for some $\delta>0$, $\sum_{v\in V}d_{P-1}(v)=O(n)$
with probability $1-n^{-\delta}$ and in expectation.

We refer the reader to [23, 24] for a complete proof of Theorem 2. We are now
ready to prove the runtime bounds of the auction algorithm.

###### Proof of runtime bounds (Theorem 1).

We recall that $P\geq\sum_{v\in V}L_{P-1}(v)+1.$ From Corollary 1 and Theorem
2 we conclude that for sparse random $k$-left regular bipartite graphs $P\leq
O(n\cdot(w^{max}-w^{min}+\varepsilon))$ with high probability. We know that
labels increase by at least $\varepsilon$ in each step. Further $k$
comparisons are required, in each move, to find the best and the second best
vertices. From these two observations we can conclude that the worst case
bound is $O(nk\cdot{w^{max}-w^{min}+\varepsilon\over\varepsilon})$ with high
probability.

For complete bipartite graphs, note that $d_{P-1}(v)\leq 2$ for all $v\in V$
(for the last free vertex, it is $0$ ) Note that in each move $n$ comparisons
are made which gives us the worst case bound of
$O\left(n^{2}\cdot{w^{max}-w^{min}\over\varepsilon}\right)$. ∎

It is easy to see that for arbitrary bipartite graphs and also for the case
where the perfect matching does not exist, the worst case runtime is bounded
by $O(n^{2}d\cdot{w^{max}-w^{min}+\varepsilon\over\varepsilon})$, where $d$ is
the maximum degree of the vertices in $U$. We note that in practice the
maximum label value will be much lower than what we estimate here as we assume
all unmatched edges to be with the largest weight and all matched edges to be
of smallest weight. Also, what appears like the worst case for the runtime
analysis, i.e., where each unmatched edge has weight $w^{max}$ and each of the
matched edge has weight $w^{min}$ and when the weight ranges are large, is in
fact a considerably easy case for the algorithm. The parameter $\varepsilon$
will not have much role to play in this case as for each label update since
the label will be increased by a high value $w^{max}-w^{min}$.

For completeness, we prove in the following section that the algorithm outputs
a $\varepsilon$-optimal solution for the case where the perfect matching
exists. We do not go into the detailed analysis for the case where a perfect
matching does not exist because of the, previously stated, bad worst case
runtime bound. We believe that a closer analysis will lead to the same
approximation guarantee for this case too. In the future, we hope to improve
the runtime bound by using scaling type approach, in which the maximum
possible value of the label is initially set to a very small value. This would
lead to discarding most of the vertices and would result in a small matching
size. In each scale we would increase the maximum value by some factor and
improve on the matching size obtained from the previous scale.

###  2.4 An $\varepsilon$-optimal solution

Using Lemma 1 it is now to easy to show that Auction algorithm outputs a near
optimal minimum weight perfect matching for all classes of bipartite graphs
and arbitrary weights provided a perfect matching exists.

######  Lemma 2.

For any $\varepsilon>0$ and given bipartite graph $G$, Auction algorithm
outputs a perfect matching with weight at most $OPT_{G}+n\varepsilon$

###### Proof.

Let us assume that Auction algorithm does not output an optimal answer. First
consider the case where $|U|<|V|=n$, such that there is a free vertex $v\in V$
at the end of the last move. As in Lemma 1 let
$B_{t}=(v,u_{1},v_{1},u_{2},v_{2},\cdots,u_{t},v_{t})$ be an augmenting path
of length $2t$ where $\mathcal{M}_{P}(B_{t})$ and
$\mathcal{M}^{\prime}_{P}(B_{t})$ denote the set of matched and unmatched
edges in $B_{t}$. By Lemma 1 we obtain

$L_{P}(v_{t})\leq
L_{P}(v)+\sum_{e\in\mathcal{M}^{\prime}_{P}(B_{t})}w(e)-\sum_{e\in\mathcal{M}_{P}(B_{t})}w(e)+t\varepsilon.$

As $L_{P}(v)=0$ and $L_{P}(v_{t})>0$ we obtain

$\sum_{e\in\mathcal{M}_{P}(B_{t})}w(e)<\sum_{e\in\mathcal{M}^{\prime}_{P}(B_{t})}w(e)+t\varepsilon$
(3)

We next consider an augmenting cycle $C_{2t}$ of length $2t$ such that

$C_{2t}=(u_{1},v_{1},u_{2},v_{2},\cdots,u_{t},v_{t},u_{1})$

where for all $i\in[t]$, $(u_{i},v_{i})\in\mathcal{M}_{P}$. WLOG we assume
that vertex $v_{t}$ was matched later than all other right set vertices in
$C_{2t}$ in some move, say $p$. Let
$B_{2t-4}=(v_{1},u_{2},v_{2},\cdots,u_{t-1},v_{t-1})$. Now by Lemma 1 we
obtain

$L_{P}(v_{t-1})\leq
L_{P}(v_{1})+\sum_{e\in\mathcal{M}_{P}(B_{2t-4})}w(e)-\sum_{e\in\mathcal{M}^{\prime}_{P}(B_{2t-4})}w(e)+(t-2)\varepsilon$

As $v_{1}$ and $v_{t-1}$ by assumption were not updated after the $(p-1)$th
move we can write the above inequality as

$L_{p-1}(v_{t-1})\leq
L_{p-1}(v_{1})+\sum_{e\in\mathcal{M}_{P}(B_{2t-4})}w(e)-\sum_{e\in\mathcal{M}^{\prime}_{P}(B_{2t-4})}w(e)+(t-2)\varepsilon,$

because $L_{P}(v_{1})=L_{p-1}(v_{1})$ and $L_{P}(v_{t-1})=L_{p-1}(v_{t-1}).$

By the choice rule we have

$\displaystyle L_{p-1}(v_{t})+w(u_{t},v_{t})\leq
L_{p-1}(v_{t-1})+w(u_{t},v_{t-1})$
$\displaystyle\stackrel{{\scriptstyle}}{{\leq}}L_{p-1}(v_{1})+\sum_{i=1}^{t-1}w(u_{i+1},v_{i})-\sum_{i=2}^{t-1}w(u_{i},v_{i})+(t-2)\varepsilon$
$\displaystyle\leq
L_{p-1}(v_{t})+w(u_{1},v_{t})-w(u_{1},v_{1})+\sum_{i=1}^{t-1}w(u_{i+1},v_{i})-\sum_{i=2}^{t-1}w(u_{i},v_{i})+(t-1)\varepsilon$
$\displaystyle\implies$
$\displaystyle\sum_{i=1}^{t}w(u_{i},v_{i})\leq\sum_{i=2}^{t}w(u_{i},v_{i-1})+w(u_{1},v_{t})+(t-1)\varepsilon$

The proof is completed by adding (3) and (2.4) for all vertex disjoint
augmenting paths and cycles. ∎

## References

  * [1] Arne Andersson, Torben Hagerup, Stefan Nilsson, and Rajeev Raman.  Sorting in linear time.  In Proceedings of the twenty-seventh annual ACM symposium on Theory of computing, pages 427–436. ACM, 1995.
  * [2] D. P. Bertsekas.  Distributed relaxation methods for linear network flow problems.  In Decision and Control, 1986 25th IEEE Conference on, volume 25, pages 2101–2106. IEEE, 1986.
  * [3] D. P. Bertsekas and D. A. Castañon.  The auction algorithm for the transportation problem.  Annals of Operations Research, 20(1):67–96, 1989.
  * [4] D. P. Bertsekas, D.A. Castañon, and Haralampos Tsaknakis.  Reverse auction and the solution of inequality constrained assignment problems.  SIAM Journal on Optimization, 3(2):268–297, 1993.
  * [5] D. P. Bertsekas and J. Eckstein.  Dual coordinate step methods for linear network flow problems.  Mathematical Programming, 42(1):203–243, 1988.
  * [6] D. P. Bertsekas and J. N. Tsitsiklis.  Parallel and distributed computation: numerical methods, volume 23.  1989\.
  * [7] D.P. Bertsekas.  A new algorithm for the assignment problem.  Mathematical Programming, 21(1):152–171, 1981.
  * [8] D.P. Bertsekas.  An auction algorithm for shortest paths.  SIAM Journal on Optimization, 1(4):425–447, 1991.
  * [9] D.P. Bertsekas.  Auction algorithms for network flow problems: A tutorial introduction.  Computational Optimization and Applications, 1(1):7–66, 1992.
  * [10] D.P. Bertsekas and J. Eckstein.  Distributed asynchronous relaxation methods for linear network flow problems.  IFAC Proceedings Volumes, 20(5):103 – 114, 1987.  10th Triennial IFAC Congress on Automatic Control - 1987 Volume VII, Munich, Germany, 27-31 July.
  * [11] E. W. Dijkstra.  A note on two problems in connexion with graphs.  Numerische mathematik, 1(1):269–271, 1959.
  * [12] R. Duan, S. Pettie, and Hsin-Hao Su.  Scaling algorithms for weighted matching in general graphs.  In Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA ’17, pages 781–800, 2017.
  * [13] Ran Duan and Seth Pettie.  Linear-time approximation for maximum weight matching.  Journal of the ACM (JACM), 61(1):1, 2014.
  * [14] Jack Edmonds and Richard M Karp.  Theoretical improvements in algorithmic efficiency for network flow problems.  Journal of the ACM (JACM), 19(2):248–264, 1972.
  * [15] N. Fountoulakis and K. Panagiotou.  Sharp load thresholds for cuckoo hashing.  Random Structures & Algorithms, 41(3):306–333, 2012.
  * [16] Michael L Fredman and Robert Endre Tarjan.  Fibonacci heaps and their uses in improved network optimization algorithms.  Journal of the ACM (JACM), 34(3):596–615, 1987.
  * [17] A. Frieze and P. Melsted.  Maximum matchings in random bipartite graphs and the space utilization of cuckoo hash tables.  Random Structures & Algorithms, 41(3):334–364, 2012.
  * [18] H. N. Gabow and R. E. Tarjan.  Faster scaling algorithms for network problems.  SIAM Journal on Computing, 18(5):1013–1036, 1989.
  * [19] Andrew V Goldberg and Robert Kennedy.  Global price updates help.  SIAM Journal on Discrete Mathematics, 10(4):551–572, 1997.
  * [20] Yijie Han.  Deterministic sorting in o (n log log n) time and linear space.  In Proceedings of the thiry-fourth annual ACM symposium on Theory of computing, pages 602–608. ACM, 2002.
  * [21] Chien-Chung Huang and Telikepalli K.  Efficient algorithms for maximum weight matchings in general graphs with small edge weights.  In Proceedings of the twenty-third annual ACM-SIAM symposium on Discrete Algorithms, pages 1400–1412, 2012.
  * [22] Ming-Yang Kao, Tak-Wah Lam, Wing-Kin Sung, and Hing-Fung Ting.  A decomposition theorem for maximum weight bipartite matchings.  SIAM Journal on Computing, 31(1):18–26, 2001.
  * [23] Megha Khosla.  Balls into bins made faster.  In European Symposium on Algorithms, pages 601–612. Springer, 2013\.
  * [24] Megha Khosla and Avishek Anand.  A faster algorithm for cuckoo insertion and bipartite matching in large graphs.  Algorithmica, 81(9):3707–3724, 2019.
  * [25] Harold W Kuhn.  The hungarian method for the assignment problem.  Naval research logistics quarterly, 2(1-2):83–97, 1955.
  * [26] Harold W Kuhn.  Variants of the hungarian method for assignment problems.  Naval Research Logistics Quarterly, 3(4):253–258, 1956.
  * [27] M. Lelarge.  A new approach to the orientation of random hypergraphs.  In Proceedings of the Twenty-Third Annual ACM-SIAM Symposium on Discrete Algorithms, SODA ’12, pages 251–264, 2012.
  * [28] James B Orlin and Ravindra K Ahuja.  New scaling algorithms for the assignment and minimum mean cycle problems.  Mathematical programming, 54(1-3):41–56, 1992.
  * [29] S. Pettie.  A simple reduction from maximum weight matching to maximum cardinality matching.  Information Processing Letters, 112(23):893–898, 2012.
  * [30] Piotr Sankowski.  Maximum weight bipartite matching in matrix multiplication time.  Theoretical Computer Science, 410(44):4480–4488, 2009.
  * [31] Mikkel Thorup.  Equivalence between priority queues and sorting.  In Foundations of Computer Science, 2002. Proceedings. The 43rd Annual IEEE Symposium on, pages 125–134. IEEE, 2002.
  * [32] Mikkel Thorup.  Integer priority queues with decrease key in constant time and the single source shortest paths problem.  In Proceedings of the thirty-fifth annual ACM symposium on Theory of computing, pages 149–158. ACM, 2003.
  * [33] N Tomizawa.  On some techniques useful for solution of transportation network problems.  Networks, 1(2):173–194, 1971.


eqsplit

$\begin{split}\BODY\end{split}$ (1)

# Charm balance function in relativistic heavy-ion collisions

Tribhuban Parida<EMAIL_ADDRESS>Department of Physical Sciences,
Indian Institute of Science Education and Research Berhampur, Transit Campus
(Govt ITI), Berhampur-760010, Odisha, India  Piotr Bożek
<EMAIL_ADDRESS>AGH University of Krakow, Faculty of Physics and
Applied Computer Science, aleja Mickiewicza 30, 30-059 Cracow, Poland  Sandeep
Chatterjee<EMAIL_ADDRESS>Department of Physical Sciences, Indian
Institute of Science Education and Research Berhampur, Transit Campus (Govt
ITI), Berhampur-760010, Odisha, India

###### Abstract

We calculate the balance function for charm in relativistic heavy-ion
collisions. The distribution of pairs of charm-anticharm quarks produced in
hard processes in the early stages of the nucleus-nucleus collision evolves in
the dense fireball formed in the collision. The evolution of the dense matter
is described using a relativistic viscous hydrodynamic model and the quark
diffusion with a Langevin equation. The evolution of the charm quark balance
function from the formation of the charm-anticharm pair up to the freeze-out
traces the partial thermalization of the heavy quarks in the dense matter. For
the balance function in azimuthal angle we reproduce the collimation effect
due to the transverse flow. The evolution in rapidity shows the thermalization
of the longitudinal velocity of the quark in the fluid. We provide predictions
for the one and two-dimensional balance functions for $D^{0}$-$\bar{D^{0}}$
mesons produced in ultarelativistic Pb+Pb collisions at $\sqrt{s_{NN}}=5.02$
TeV. The shape of the charm balance function in relative rapidity is sensitive
to the rescattering of heavy quarks in the early stages of the collision,
while the shape of the balance function in azimuthal angle is sensitive to the
rescattering in the latter stages.

ultra-relativistic nuclear collisions, correlations

##  I Introduction

Heavy quarks are produced in hard processes at the beginning of relativistic
nuclear collisions. The interaction of such particles produced at the early
stages with the surrounding medium makes them an excellent probe of the
dynamics of the dense matter created in the interaction region of the
collision Svetitsky (1988); Prino and Rapp (2016); Dong  _et al._ (2019); Dong
and Greco (2019); He  _et al._ (2023). Measurements of the nuclear
modification factor and the harmonic flow of open charm mesons provide an
insight into the partial thermalization of hard produced charm quarks.

The charge balance function has been proposed as a measure of the interaction
of charged quarks in the quark-gluon plasma Bass  _et al._ (2000). The two
opposite charges produced in the plasma would drift apart in relative rapidity
due to rescattering with the thermal medium. The charge balance function is
defined in a way to identify the specific correlated charge-anticharge pairs.
The balance function in the relative azimuthal angle is sensitive to the
amount of transverse flow at the time of the charge pair separation Bożek
(2005). The balance function involving the electric charge, strangeness and
baryon number can be measured Adams  _et al._ (2003); Aggarwal  _et al._
(2010); Alt  _et al._ (2005, 2007); Abelev  _et al._ (2013); Tumasyan  _et
al._ (2023). Charges carried by light quarks and hadrons can be produced
during the entire span of the collision. In particular, a large part of the
balance function for the electric charge receives contribution from local
charge conservation for charges produced at hadronization Jeon and Pratt
(2002). The full description of the charge balancing dynamics combines the
charge-anticharge pairs produced in the medium and at the hadronization Pratt
(2012a, b); Pratt and Plumberg (2019). Charm quarks are almost exclusively
produced in the initial hard scatterings. The two balancing charm-anticharm
quarks evolve in the medium. It has been noted that the charm-anticharm
correlations or the charm balance function could serve as a clean probe of
such interactions Zhu  _et al._ (2007); Basu  _et al._ (2021). For example,
the charm balance function could be sensitive to the formation of charmonium
bound states in the quark-gluon plasma.

In this paper, we study the charm-anticharm balance function in relative
azimuthal angle and relative rapidity. The heavy quarks produced in hard
processes evolve in the quark-gluon plasma, which changes their relative
momenta. In comparison to the dynamics of light quarks, the hard produced
heavy quarks have large initial momenta. In our calculation, we take charm and
anticharm quarks as produced in p+p collisions using the Pythia model Bierlich
_et al._ (2022). The quarks rescatter in the dense matter of the fireball. In
relative azimuthal angle the balance function becomes narrower. This is
equivalent to the effect of the increase of the correlation function of the
final $D^{0}$ and $\bar{D}^{0}$ mesons at small relative azimuthal angle Zhu
_et al._ (2007, 2008); Younus and Srivastava (2013); Tsiledakis  _et al._
(2011); Nahrgang  _et al._ (2014); Adare  _et al._ (2014). The balance
function in relative rapidity becomes broader after the diffusion of the heavy
quarks in the medium. The behavior of the charm balance in relative rapidity
provides an additional information, because it is sensitive to the early
stages of the reaction, unlike the balance function in azimuthal angle.

In the next section we describe the calculation of the balance function in the
model. Then, in Sec. III the balance function is calculated directly for the
charm quarks before and after rescattering in the medium. In Sec. IV we study
the time evolution of different contributions leading to the broadening of the
charm balance function in rapidity. Predictions for the balance function for
$D^{0}$-$\bar{D}^{0}$ mesons are provided in Sec. V. A simple estimate of
possible effects of a pre-hydrodynamic phase on the charm balance function is
given in Sec. VI. Conclusions about the charge balance functions in relative
azimuthal angle and relative rapidity are summarized in the last section.

##  II Calculation of the charm balance function in the hydrodynamic model

We have adopted the Optical Glauber model to set up the initial condition for
the hydrodynamic evolution. The initial three dimensional distribution of the
energy density ($\epsilon$) at a constant proper time ($\tau_{0}$) takes the
following form:

$\displaystyle\epsilon(x,y,\eta_{s})$ $\displaystyle=$
$\displaystyle\epsilon_{0}\left[\left(N_{+}(x,y)f_{+}(\eta_{s})+N_{-}(x,y)f_{-}(\eta_{s})\right)\right.$
(2)
$\displaystyle\left.\times\left(1-\alpha\right)+N_{coll}(x,y)\epsilon_{\eta_{s}}\left(\eta_{s}\right)\alpha\right]$

where, $N_{+}(x,y)$ and $N_{-}(x,y)$ are the participant densities from the
forward and backward going nuclei. $N_{coll}(x,y)$ is the density of the
binary collision sources evaluated at any transverse position (x, y). $\alpha$
is the hardness factor and $\epsilon_{0}$ is the scaling parameter of energy
density. $\epsilon_{\eta_{s}}\left(\eta_{s}\right)$ and $f_{+,-}(\eta_{s})$
are the profiles that determine the deposition of $\epsilon$ in space-time
rapidity ($\eta_{s}$). The forms of $\epsilon_{\eta_{s}}\left(\eta_{s}\right)$
and $f_{+,-}(\eta_{s})$ are given in Eq. 3 and Eq. 4 respectively,

$\epsilon_{\eta_{s}}(\eta_{s})=\exp\left(-\frac{\left(|\eta_{s}|-\eta_{0}\right)^{2}}{2\sigma_{\eta}^{2}}\theta(|\eta_{s}|-\eta_{0})\right)\
,$ (3)

$f_{+,-}(\eta_{s})=\epsilon_{\eta_{s}}(\eta_{s})\epsilon_{F,B}(\eta_{s})\ ,$
(4)

where

$\epsilon_{F}(\eta_{s})=\begin{cases}0,&\text{if }\eta_{s}<-\eta_{m}\\\
\frac{\eta_{s}+\eta_{m}}{2\eta_{m}},&\text{if
}-\eta_{m}\leq\eta_{s}\leq\eta_{m}\\\ 1,&\text{if
}\eta_{m}<\eta_{s}\end{cases}$ (5)

and

$\epsilon_{B}(\eta_{s})=\epsilon_{F}(-\eta_{s})\ .$ (6)

$f_{+,-}(\eta_{s})$ introduces an asymmetric deposition of energy density by a
participant source along $\eta_{s}$. This kind of initial condition generates
a tilted profile in the reaction plane (plane made by beam axis and impact
parameter) of the collision Bożek and Wyskiel (2010). The strength of the
initial tilt in the model can be controlled by the parameter $\eta_{m}$
present in Eq. 5. We have used $\alpha=0.14,\eta_{0}=1.8$, $\sigma_{\eta}=2.0$
and $\eta_{m}=3.2$ that describe the experimental data of the charged particle
yield and directed flow in Pb+Pb collisions at 5.02 TeV.

The viscous hydrodynamic evolution of the energy density profile has been
performed by the publicly available MUSIC code Schenke  _et al._ (2010, 2011);
Paquet  _et al._ (2016). We have taken shear viscosity to entropy density
ratio $\eta/s=0.04$ Alba  _et al._ (2018); Giacalone  _et al._ (2018) and have
used a lattice QCD based Equation of State Monnai  _et al._ (2019) during the
hydrodynamic evolution.

The hydrodynamic evolution provides us the relevant dynamic properties of the
background (local temperature $T$ and flow velocity $u^{\mu}$), in which the
charm and anticharm quarks evolve. The quark evolution stops once the local
background temperature falls below 150 MeV. The dynamics of the heavy quarks
in the medium is described using the Langevin equation

$\displaystyle\Delta{\bf r}_{i}$ $\displaystyle=$ $\displaystyle\frac{{\bf
p}_{i}}{E}\Delta t$ $\displaystyle\Delta{\bf p}_{i}$ $\displaystyle=$
$\displaystyle-\gamma{\bf p}_{i}\Delta t+\rho_{i}\sqrt{2D\Delta t}\ ,$ (7)

where $\Delta{\bf r_{i}}$ and $\Delta{\bf p_{i}}$ are the updates of the
coordinates in position and momentum of the heavy quark respectively within
the time interval $\Delta t$. $\gamma$ and $D$ are the drag and diffusion
coefficients, which effectively describe the interactions between the heavy
quark and the medium. We assume a diagonal form for the diffusion matrix van
Hees  _et al._ (2006); Cao and Bass (2011); Scardina  _et al._ (2017). The
noise term $\rho_{i}$ is randomly sampled from a normal distribution at every
time step, with $\langle\rho_{i}\rangle=0$ and
$\langle\rho_{i}\rho_{j}\rangle=\delta_{ij}$. In the post-point realization of
the stochastic evolution, we take $D=\gamma ET$, where
$E=\sqrt{p^{2}+m_{c}^{2}}$ is the energy of the heavy quark with mass $m_{c}$.
We consider $\gamma=0.4T$, that captures the suppression factor and the
elliptic flow of heavy quarks Chatterjee and Bożek (2018).

The balance function is defined as a conditional probability isolating
correlations between balancing charges Bass  _et al._ (2000); Jeon and Pratt
(2002); Schlichting and Pratt (2011)

$\displaystyle B(p_{1}|p_{2})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[\frac{N^{c\bar{c}}(p_{1},p_{2})}{N^{\bar{c}}(p_{2})}-\frac{N^{\bar{c}\bar{c}}(p_{1},p_{2})}{N^{\bar{c}}(p_{2})}\right.$
(8)
$\displaystyle\left.+\frac{N^{\bar{c}c}(p_{1},p_{2})}{N^{{c}}(p_{2})}-\frac{N^{{c}{c}}(p_{1},p_{2})}{N^{{c}}(p_{2})}\right]\
,$

where $N^{a}(p_{1})$ and $N^{ab}(p_{1},p_{2})$ are the one and the two-
particle distributions of particles of type $a,b$ with momenta $p_{1}$ and
$p_{2}$. The balance functions have been studied in models and in experiments
for the electric charge, strangeness and baryon number Adams  _et al._ (2003);
Aggarwal  _et al._ (2010); Alt  _et al._ (2005, 2007); Abelev  _et al._
(2013); Tumasyan  _et al._ (2023). Theses charges can be produced in the
initial nucleon-nucleon scattering, in the evolving hot and dense fireball,
and during the hadronization process. The understanding of the charge
balancing in heavy-ion collisions requires in that case a careful modeling of
the different ways of the charge-anticharge production Pratt (2012a, b); Pratt
_et al._ (2015); Pratt and Plumberg (2021); Pruneau  _et al._ (2023a, b);
Pratt  _et al._ (2022). However, owing to their large mass, the charm quarks
are expected to be produced only during the initial hard scatterings.

Assuming the independent evolution of $c$-$\bar{c}$ pairs created in different
hard process, one can write balance function using the correlation function
calculated for each $c$-$\bar{c}$ pair separately. In this work, we do not
take into account possible correlations due to the formation and dissolution
of charmonium bound states Basu  _et al._ (2021). The quark and antiquark are
evolved in the bulk matter according to Langevin equation. The two-particle
distribution $\rho^{c\bar{c}}(p_{1},p_{2})$ and the one-particle distributions
$\rho^{c}(p_{1})$, $\rho^{\bar{c}}(p_{1})$ of hadrons containing the charm
quark and antiquark can be obtained as an average over all the pairs. These
distributions are normalized to one:

$\int dp_{1}dp_{2}\rho(p_{1},p_{2})=\int dp_{1}\rho^{c}(p_{1})=\int
dp_{1}\rho^{\bar{c}}(p_{1})=1\ .$ (9)

The two-dimensional balance function can be easily constructed as

$\displaystyle B(\Delta y,\Delta\phi)=\int_{A}dp_{1}\int_{A}dp_{2}$
$\displaystyle\frac{\rho^{c\bar{c}}(p_{1},p_{2})\delta\left(\Delta
y-y_{1}+y_{2}\right)\delta\left(\Delta\phi-|\phi_{1}-\phi_{2}|\right)}{\int_{A}dp_{1}\rho^{c}(p_{1})}$
$\displaystyle\ ,$ (10)

where we assume a symmetric distribution between $c$ and $\bar{c}$ quarks and
use a symmetric choice of the acceptance region $A$ in azimuthal angle and
rapidity. In this paper we study the balance function in relative azimuthal
angle and/or in relative rapidity, where the acceptance region is given by
cuts on the transverse momentum of the particles. If the region $A$ involves
all particles containing $c$ and $\bar{c}$ quarks and all momenta, the balance
function is normalized

$\int_{-\infty}^{\infty}d(\Delta y)\int_{0}^{\pi}d(\Delta\phi)B(\Delta
y,\Delta\phi)=1\ .$ (11)

With limited kinematic acceptance and efficiency or if not all open charm and
charmonium hadrons are accounted for, the normalization is $<1$.

##  III Balance function for $c$-$\bar{c}$ quarks

In this section we study the balance functions for the $c-\bar{c}$ quarks to
understand the evolution of the balance function in the collision. The two-
dimensional balance function for quark pairs emitted from p+p collisions at
$\sqrt{s_{NN}}=5.02$ TeV from the Pythia model is shown in panel (a) of Fig.
1. The correlation has a broad distribution in $\Delta\phi$ with a maximum for
back-to-back pairs. In the panel (b), we show the balance function for
$c$-$\bar{c}$ pairs emitted after Langevin diffusion in the dense fireball
formed in Pb+Pb collisions with centrality $10-20\%$.

After diffusion in the bulk matter the quark-antiquark momenta pick up a
contribution of the transverse flow, as a result the azimuthal distribution
becomes more peaked for pairs flying in the same direction. The longitudinal
distribution in rapidity is broadened due to the rescattering in the medium.

Figure 1: The two dimensional balance function for $c$-$\bar{c}$ quarks in
relative rapidity $\Delta y$ and relative azimuthal angle $\Delta\phi$ for
$p_{T}>$ 2 GeV. The panel (a) represents simulation results from Pythia for
quarks emitted in p+p collisions and the panel (b) shows the balance function
after the evolution of the quark-antiquark pairs in the bulk matter formed in
Pb+Pb collisions with centrality $10-20\%$.  Figure 2: The balance function,
$B(\Delta y)$, for $c$-$\bar{c}$ quarks in relative rapidity $\Delta y$, for
p+p collisions (solid line) and Pb+Pb collisions with centrality $10-20\%$
(dashed line), $40-50\%$ (dashed-dotted line), and $60-70\%$ (dotted line).
Quarks with any momentum are included in the balance function.

In the following we use the one dimensional balance function in relative
rapidity,

$B(\Delta y)=2\int_{0}^{\pi}d(\Delta\phi)B(\Delta y,\Delta\phi)\ ,$ (12)

and relative azimuthal angle

$B(\Delta\phi)=\int_{-\infty}^{\infty}d(\Delta y)B(\Delta y,\Delta\phi)\ .$
(13)

Due to symmetry we define the one-dimensional balance functions $B(\Delta y)$
and $B(\Delta\phi)$ in the ranges $y\geq 0$ and $[0$,$\pi]$ respectively.

The balance function $B(\Delta y)$ is shown in Fig. 2. The plots indicate that
the balance function for quarks undergoing rescattering in the medium created
in Pb+Pb collisions become broader as compared to the balance function for
quarks emitted in p+p collisions. There is some dependence of the shape of the
balance function on collision centrality in Pb+Pb collisions. The balance
function in relative azimuthal angle $B(\Delta\phi)$ shows a qualitative
difference between p+p and Pb+Pb collisions (Fig. 3). The balance function is
peaked at $\Delta\phi=\pi$ in p+p collisions (back-to-back pairs), while in
Pb+Pb collisions it is peaked at $\Delta\phi=0$ (collimated pairs). The peak
of the balance function at $\Delta\phi=0$ is an effect of the transverse flow.
We recover the known collimation of the heavy quarks due to the rescattering
in the expanding medium Zhu  _et al._ (2007, 2008); Younus and Srivastava
(2013); Tsiledakis  _et al._ (2011); Nahrgang  _et al._ (2014); Adare  _et
al._ (2014).

The one-dimensional balance functions for quarks and antiquarks with $p_{T}>2$
GeV are shown in Figs. 4 and 5. The cut $p_{T}>2$ GeV makes the balance
function in relative rapidity narrower. The effect of the $p_{T}$ cut on the
balance functions in relative azimuthal angle is to make the function more
peaked around $\Delta\phi=\pi$ for p+p collisions and around $\Delta\phi=0$
for Pb+Pb collisions. Note that the normalization of the balance functions,
Eq. 11, measured in a restricted part of phase-space ($p_{T}>2$ GeV) is
modified:

$\displaystyle\int_{-\infty}^{\infty}d(\Delta
y)\int_{0}^{\pi}d(\Delta\phi)B(\Delta y,\Delta\phi)<1$ (14)

and similarly for the one-dimensional balance functions.

Figure 3: Same as in Fig. 2 but for the balance function, $B(\Delta\phi)$, in
relative azimuthal angle.  Figure 4: Same as in Fig. 2 but for quarks and
antiquarks with $p_{T}>2$ GeV.

For the study of the balance function for charm quarks, we use the root mean
square (RMS) width of the balance function in relative rapidity,
$\sigma_{\Delta y}=\langle\Delta y^{2}\rangle^{1/2}$ and not the simple width,
$\langle|\Delta y|\rangle$. This makes simpler the discussion of different
contributions to the total broadening of the balance function. The RMS width
in Pb+Pb collisions is larger than in p+p collisions (Fig. 6). The quark-
antiquark pair undergoes additional rescattering in the fireball widening the
distribution. The width shows a clear centrality dependence. The calculation
with reduced starting time for the hydrodynamic evolution of the fireball
($\tau_{0}=0.2$ fm/c instead of $0.4$ fm/c) shows slightly larger width of the
balance function. The mechanism responsible for the widening of the balance
function in rapidity is discussed in the following section.

The width of the balance function in relative azimuthal angle,
$\sigma_{\Delta\phi}=\langle\Delta\phi^{2}\rangle^{1/2}$ , shows a strong
centrality dependence (Fig. 7). The longer the $c$-$\bar{c}$ quark pair
interacts with the expanding fluid, the more collimated the momenta of two
quarks become.

Figure 5: Same as in Fig. 3 but for quarks and antiquarks with $p_{T}>2$ GeV.
Figure 6: RMS width of the balance function in relative rapidity for
$c$-$\bar{c}$ pairs produced in Pb+Pb collisions as a function of collision
centrality. The dotted line and dashed lines represent the results for the
initial time $\tau_{0}=0.2$ and $0.4$fm/c respectively. The horizontal solid
line denotes the width of the $c-\bar{c}$ distribution in p+p collisions. All
results are for quarks with $p_{T}>2$ GeV. Figure 7: Same as in Fig. 6, but
for the balance function in relative azimuthal angle.

##  IV Evolution of the $c$-$\bar{c}$ distribution

Initially, the charm quark has typically a large transverse momentum, as
compared to the background thermal motion. During the evolution, the quark
rescatters in the medium, losing part of its energy and transverse momentum
and the direction of the transverse momentum becomes more aligned with the
transverse collective flow of the fluid. As a result, the distribution in
relative azimuthal angle becomes peaked at $\Delta\phi=0$ Zhu  _et al._ (2007,
2008); Younus and Srivastava (2013); Tsiledakis  _et al._ (2011); Nahrgang
_et al._ (2014); Adare  _et al._ (2014).

Figure 8: The time evolution of the spread $\langle(y-Y)^{2}\rangle$ of the
relative rapidity of the heavy quark $y$ and the rapidity of the fluid for
Pb+Pb collisions with $10$-$20$% centrality (solid line). The dotted line
represents the thermal spreads of the relative rapidity for thermally
equilibrated longitudinal velocity of the quark in the fluid (Eq. 16).

When a heavy quark is initialized in the fluid, its rapidity with respect to
the fluid is zero (in the scenario when both the heavy quarks and the fluid
follow the Bjorken scaling flow in the longitudinal direction). The
rescattering with the thermal medium leads to a spread of the relative
rapidity distribution. The thermal distribution in the difference of the quark
rapidity $y$ and the fluid rapidity $Y$ is:

$f_{th}(y-Y)=\frac{(1+\sinh(y-Y))^{2}}{2\pi
E_{T}T}\exp\left(-\frac{\sinh(y-Y)^{2}E_{T}}{2T}\right),$ (15)

where $E_{T}=\sqrt{m_{c}^{2}+p_{T}^{2}}$ is the transverse mass of the quark.

The average spread of the relative rapidity during the evolution of heavy
quarks in the medium formed in a Pb+Pb collision with centrality $10$-$20$% is
shown in Fig. 8. The spread grows rapidly from the initial value towards the
average thermal spread:

$\sigma_{th}^{2}=\int_{-\infty}^{\infty}dyy^{2}f_{th}(y)\ ,$ (16)

shown with the dotted line in Fig. 8. Due to the Bjorken expansion of the
system, the longitudinal momentum distribution of heavy-quarks never reaches
the equilibrium distribution Florkowski  _et al._ (2013). During the
evolution, the average transverse momentum of the heavy quark decreases, also
the temperature of the surrounding medium decreases. The net effect is that
the thermal spread decreases in time. The final thermal spread at the freeze-
out is small and gives a very small contribution to the broadening of the
balance function in the medium.

Figure 9: The width of the charm balance function as a function of collision
centrality for the initialization with both quarks with the same rapidity and
same position. Figure 10: The time evolution of the width of the charm balance
function as in Pb+Pb collisions in centrality range $10$-$20$% (solid line).
The dotted line represents the broadening for two quarks initialized at the
same rapidity.

The random evolution of the quark-antiquark pair in the medium leads to a
diffusion in rapidity and to an increase of the average rapidity separation
between the two charges Bass  _et al._ (2000). To test this effect we perform
a simulation with both the $c$ and the $\bar{c}$ quarks initialized at the
same rapidity, $y=\frac{y_{1}+y_{2}}{2}$ and with the same space-time
rapidity, $\eta_{s}=y$. Here $y_{1}$ and $y_{2}$ are the rapidities of the
$c$-$\bar{c}$ pair as sampled from p+p collision in PYTHIA. Initially the
width of the balance function is zero and increases in time due to the
diffusion and thermal spread. The total broadening in Fig. 9 is much larger
than the thermal broadening plotted Fig. 8.

In Fig. 10 is plotted the time evolution of the width of the charm balance
function in the actual realistic simulation. Notice, that the width of the
balance function increases strongly in the first stage of the diffusion. It
shows, that the increase of the width of the balance function in rapidity is
sensitive to the early stage of the evolution. Further, Fig. 10 also shows the
broadening expected from the diffusion for two quarks initialized at the same
space-time rapidity. It is slightly smaller than the actual broadening in the
full calculation. In the full calculation the two quarks are initialized in
two different fluid cells, with different rapidities. During a three-
dimensional hydrodynamic evolution the fluid accelerates in the longitudinal
direction, which increases the rapidity of the two fluid elements, leading to
violation of the Bjorken scaling flow by about $3$-$5$%. This effect increases
the relative rapidity separation of the two fluid elements and of the two
quarks as well.

##  V Charm balance function for $D^{0}-\bar{D^{0}}$ mesons

The charm quarks hadronize into open charm mesons and baryons. In this work,
we perform the hadronization at the freeze-out hypersurface of the medium
which is obtained at a constant temperature of 150 MeV. We model the
hadronization using the fragmentation model of Peterson et al. Peterson  _et
al._ (1983). Constructing the charm balance function using all open charm
mesons could recover almost all of the balancing charm charges. In this paper
we present prediction for the balance function of $D^{0}$-$\bar{D}^{0}$ mesons
only. The probability of producing a $D^{0}$ meson from a charm quark is
around $40$% Acharya  _et al._ (2022). This means that only about $40$% of the
balancing anticharm would be recovered in $\bar{D}^{0}$ for each observed
$D^{0}$, assuming full acceptance and efficiency.

Figure 11:  The balance function in relative rapidity for
$D^{0}$-$\bar{D^{0}}$ mesons with $p_{T}>2$GeV in Pb+Pb collisions with
centrality range $10$-$20$% (dashed line), $40$-$50$% (dashed-dotted line),
and $60$-$70$% (dotted line). The solid line represents the balance function
for p+p collisions. Figure 12:  Same as in Fig. 11 but for the balance
function in relative azimuthal angle. Figure 13:  The width $\langle|\Delta
y|\rangle$ of the balance function for $D^{0}$-$\bar{D^{0}}$ mesons as a
function of collision centrality. The thin solid and dashed lines represent
the results for the initial time $\tau_{0}=0.2$ and $0.4$fm/c respectively.
The horizontal solid line denotes the width of the $D^{0}$-$\bar{D}^{0}$
distribution in p+p collisions. All results are for mesons with $p_{T}>2$ GeV.
Figure 14:  Same as in Fig. 13 but for the width $\langle|\Delta\phi|\rangle$
of the balance function in relative azimuthal angle.

The balance function measured for the open charm mesons $D^{0}$-$\bar{D}^{0}$
represents qualitatively the same features as the balance function for quarks.
The balance function in relative rapidity gets broader due to the diffusion of
charm quarks in the medium (Fig. 11). As expected, normalization of the
balance function is modified, as not all the charm quarks are identified by
measuring only the $D^{0}$ mesons. The balance in the relative azimuthal angle
indicates that a collimation of the momenta of the $D^{0}$ and $\bar{D}^{0}$
meson occurs due to the interaction of the heavy quarks with the expanding
medium (Fig. 12).

For the prediction of the balance function for $D^{0}$-$\bar{D}^{0}$ pairs, we
use the average width (as typically used by experimental groups), not the RMS
width as for the discussion of the charm quark balance function in previous
sections. The width of the balance function in relative rapidity increase as
compared to p+p collisions (Fig. 13). The width, $\langle|\Delta y|\rangle$,
shows some centrality dependence. A change in the starting time, $\tau_{0}$,
for the quark evolution in the quark-gluon plasma, has a sizable effect on the
width of $D^{0}$ meson balance function. These effects indicate a possible
sensitivity of the width of the charm meson balance function to the dynamics
of the dense matter in the collision. Especially interesting is the
sensitivity to the rescattering of heavy quarks in the early stages of the
dynamics. The balance function of $D^{0}$ mesons in relative rapidity could
serve as probe of the formation time of the dense matter.

The width of the balance function in relative azimuthal angle is a very
sensitive probe of the degree of thermalization of the direction of heavy
quarks in the medium (Fig. 14). The width depends strongly on the collision
centrality. The collimation of open charm mesons comes from the drag of the
heavy quarks in the expanding medium. This mechanism is active only in the
latter stages of the collision. In the early phase of the collision, when the
transverse collective flow is not yet present, the rescattering leads to an
isotropization in the relative azimuthal angle, but not a collimation.
Therefore, the width of the balance function is not sensitive to the staring
time of the dynamics, $\tau_{0}$. On the other hand, the width,
$\langle\Delta\phi\rangle$, could serve as a clock for the total lifetime of
the quark-gluon plasma.

##  VI Effect of early dynamics on the balance function

The formation time of the heavy quark pair can be estimated as $1/2m_{c}\simeq
0.07$ fm/c. In the previous section are presented results assuming the
rescattering of the heavy quarks starting at the initial time for
hydrodynamics $\tau_{0}=0.4$ or $0.2$ fm/c. However, the evolution of the
momenta of heavy quarks in the very early, preequilibrium phase could change
the charm balance as well Avramescu  _et al._ (2023); Boguslavski  _et al._
(2023). In the present manuscript, we estimate this effect assuming that the
dynamics of heavy quarks in the early phase is still given by the Langevin
equation and by extrapolating the preequilibrium evolution of the density of
the medium down to $\tau_{f}$. For the evolution of the energy density in the
earliest phase, $\tau<\tau_{0}$, we use the assumption :

$\epsilon(\tau,x,y,\eta_{s})=\epsilon(\tau_{0},x,y,\eta_{s})\left(\frac{\tau_{0}}{\tau}\right)^{\alpha}\
,$ (17)

with $\alpha=4/3$ (ideal hydrodynamics with conformal equation of state), $1$
(zero longitudinal pressure), or $0$ (longitudinal electric field) Vredevoogd
and Pratt (2009). For the flow velocity we take the Bjorken flow, with no
transverse flow for $\tau<\tau_{0}$.

Figure 15:  The width, $\langle|\Delta y|\rangle$, of the balance function for
$D^{0}$-$\bar{D^{0}}$ mesons as a function of collision centrality obtained
with different assumptions on the evolution of the energy density in the early
phase. Results assuming a hydrodynamic evolution starting at $\tau_{0}=0.2$
fm/c and a preequilibrium evolution of the energy density,
$\epsilon(\tau)=\epsilon(\tau_{0})\left(\frac{\tau_{0}}{\tau}\right)^{\alpha}$
from $\tau_{f}=0.07$fm/c to $\tau_{0}=0.2$ fm/c are plotted for $\alpha=4/3$
(dotted line), $1$ (dashed line), and $0$ (dashed-dotted line). Results
obtained with a hydrodynamic evolution only, starting at $\tau_{0}=0.2$ fm/c,
are represented with the thin solid line. Figure 16:  Same as in Fig. 15 but
for the width, $\langle|\Delta\phi|\rangle$, of the balance function in
relative azimuthal angle.

We notice that the early diffusion of heavy quarks has only an effect on the
width of the balance function in relative rapidity (Fig. 15). On the other
hand, the width of the balance function in the relative azimuthal angle shows
almost no sensitivity to the evolution the heavy quarks in the early phase
(Fig. 16). This observation is in agreement with the results in Sec. V. It
would be interesting to study such effect using specific models of the
evolution of the heavy quarks in the preequilibrium phase Sun  _et al._
(2019); Avramescu  _et al._ (2023); Boguslavski  _et al._ (2023) or including
the effect of electromagnetic fields Das  _et al._ (2017). Further, it will be
also interesting to understand how the breaking of the Bjorken scaling in the
initialization of the charm quarks as well as the flow velocity of the medium
affects the charm balance function in relative rapidity.

##  VII Summary

The dynamics of charm-anticharm balancing is studied for relativistic Pb+Pb
collisions. Pairs of charm and anticharm quarks are formed in the initial hard
collisions. The two charges evolve in the dense matter of the fireball and
finally hadronize into open charm mesons. The charm balance function allows to
measure directly the distribution of balancing charm hadrons corresponding to
the early produced charm quarks.

The charm balance function becomes broader in relative rapidity and narrower
in relative azimuthal angle of the quarks. The dominant contribution to the
broadening of the charm balance in rapidity comes from diffusion of the two
quarks in rapidity. The thermal spreading of the quark rapidity is small. We
observe a centrality dependence in the widening of the charm balance function
in relative rapidity, with the largest broadening occurring in central
collisions. The balance function in relative azimuthal angle becomes narrower
after rescattering of the heavy quarks in the expanding medium. The
collimation of the momenta of the two heavy quarks is due to transverse
collective flow of the fluid. This collimation effect has a strong centrality
dependence.

The measurements of the charm balance function in relative rapidity and in
relative azimuthal angle are complementary. The narrowing of the balance
function in relative azimuthal angle is sensitive to the rescattering of heavy
quarks in latter stages of the collisions, when the transverse flow is
substantial. The broadening of the balance function in relative rapidity is
sensitive to early phase of the dynamics, when the rescattering of heavy
quarks is the strongest. It would be interesting to compare the charm balance
functions in rapidity and azimuthal angle obtained in model calculations to
results of future, high efficiency experiments. Such observations could
constrain both the total lifetime and the formation time of the dense fireball
formed in relativistic heavy-ion collisions.

## Acknowledgment

P.B. acknowledges support from the National Science Centre grant
2018/29/B/ST2/00244. SC acknowledges IISER Berhampur for seed grant.

## References

  * Svetitsky (1988) B. Svetitsky, Phys. Rev. D 37, 2484 (1988).
  * Prino and Rapp (2016) F. Prino and R. Rapp, J. Phys. G 43, 093002 (2016), arXiv:1603.00529 [nucl-ex] .
  * Dong  _et al._ (2019) X. Dong, Y.-J. Lee, and R. Rapp, Ann. Rev. Nucl. Part. Sci. 69, 417 (2019), arXiv:1903.07709 [nucl-ex] .
  * Dong and Greco (2019) X. Dong and V. Greco, Prog. Part. Nucl. Phys. 104, 97 (2019).
  * He  _et al._ (2023) M. He, H. van Hees, and R. Rapp, Prog. Part. Nucl. Phys. 130, 104020 (2023), arXiv:2204.09299 [hep-ph] .
  * Bass  _et al._ (2000) S. A. Bass, P. Danielewicz, and S. Pratt, Phys. Rev. Lett. 85, 2689 (2000), arXiv:nucl-th/0005044 .
  * Bożek (2005) P. Bożek, Phys. Lett. B609, 247 (2005), arXiv:nucl-th/0412076 [nucl-th] .
  * Adams  _et al._ (2003) J. Adams _et al._ (STAR), Phys. Rev. Lett. 90, 172301 (2003), arXiv:nucl-ex/0301014 .
  * Aggarwal  _et al._ (2010) M. M. Aggarwal _et al._ (STAR), Phys. Rev. C 82, 024905 (2010), arXiv:1005.2307 [nucl-ex] .
  * Alt  _et al._ (2005) C. Alt _et al._ (NA49), Phys. Rev. C 71, 034903 (2005), arXiv:hep-ex/0409031 .
  * Alt  _et al._ (2007) C. Alt _et al._ (NA49), Phys. Rev. C 76, 024914 (2007), arXiv:0705.1122 [nucl-ex] .
  * Abelev  _et al._ (2013) B. Abelev _et al._ (ALICE), Phys. Lett. B 723, 267 (2013), arXiv:1301.3756 [nucl-ex] .
  * Tumasyan  _et al._ (2023) A. Tumasyan _et al._ (CMS), (2023), arXiv:2307.11185 [nucl-ex] .
  * Jeon and Pratt (2002) S. Jeon and S. Pratt, Phys. Rev. C65, 044902 (2002), arXiv:hep-ph/0110043 .
  * Pratt (2012a) S. Pratt, Phys. Rev. C 85, 014904 (2012a), arXiv:1109.3647 [nucl-th] .
  * Pratt (2012b) S. Pratt, Phys. Rev. Lett. 108, 212301 (2012b), arXiv:1203.4578 [nucl-th] .
  * Pratt and Plumberg (2019) S. Pratt and C. Plumberg, Phys. Rev. C 99, 044916 (2019), arXiv:1812.05649 [nucl-th] .
  * Zhu  _et al._ (2007) X. Zhu, M. Bleicher, S. L. Huang, K. Schweda, H. Stoecker, N. Xu, and P. Zhuang, Phys. Lett. B 647, 366 (2007), arXiv:hep-ph/0604178 .
  * Basu  _et al._ (2021) S. Basu, P. Christiansen, A. Ohlson, and D. Silvermyr, Eur. Phys. J. C 81, 1024 (2021), arXiv:2110.05134 [hep-ph] .
  * Bierlich  _et al._ (2022) C. Bierlich _et al._ ,  (2022), 10.21468/SciPostPhysCodeb.8, arXiv:2203.11601 [hep-ph] .
  * Zhu  _et al._ (2008) X. Zhu, N. Xu, and P. Zhuang, Phys. Rev. Lett. 100, 152301 (2008), arXiv:0709.0157 [nucl-th] .
  * Younus and Srivastava (2013) M. Younus and D. K. Srivastava, J. Phys. G 40, 065004 (2013), arXiv:1301.5762 [nucl-th] .
  * Tsiledakis  _et al._ (2011) G. Tsiledakis, H. Appelshauser, K. Schweda, and J. Stachel, Nucl. Phys. A 858, 86 (2011), arXiv:0908.0427 [nucl-ex] .
  * Nahrgang  _et al._ (2014) M. Nahrgang, J. Aichelin, P. B. Gossiaux, and K. Werner, Phys. Rev. C 90, 024907 (2014), arXiv:1305.3823 [hep-ph] .
  * Adare  _et al._ (2014) A. M. Adare, M. P. McCumber, J. L. Nagle, and P. Romatschke, Phys. Rev. C 90, 024911 (2014), arXiv:1307.2188 [nucl-th] .
  * Bożek and Wyskiel (2010) P. Bożek and I. Wyskiel, Phys. Rev. C81, 054902 (2010), arXiv:1002.4999 [nucl-th] .
  * Schenke  _et al._ (2010) B. Schenke, S. Jeon, and C. Gale, Phys. Rev. C 82, 014903 (2010), arXiv:1004.1408 [hep-ph] .
  * Schenke  _et al._ (2011) B. Schenke, S. Jeon, and C. Gale, Phys. Rev. Lett. 106, 042301 (2011), arXiv:1009.3244 [hep-ph] .
  * Paquet  _et al._ (2016) J.-F. Paquet, C. Shen, G. S. Denicol, M. Luzum, B. Schenke, S. Jeon, and C. Gale, Phys. Rev. C 93, 044906 (2016), arXiv:1509.06738 [hep-ph] .
  * Alba  _et al._ (2018) P. Alba, V. Mantovani Sarti, J. Noronha, J. Noronha-Hostler, P. Parotto, I. Portillo Vazquez, and C. Ratti, Phys. Rev. C 98, 034909 (2018), arXiv:1711.05207 [nucl-th] .
  * Giacalone  _et al._ (2018) G. Giacalone, J. Noronha-Hostler, M. Luzum, and J.-Y. Ollitrault, Phys. Rev. C 97, 034904 (2018), arXiv:1711.08499 [nucl-th] .
  * Monnai  _et al._ (2019) A. Monnai, B. Schenke, and C. Shen, Phys. Rev. C 100, 024907 (2019), arXiv:1902.05095 [nucl-th] .
  * van Hees  _et al._ (2006) H. van Hees, V. Greco, and R. Rapp, Phys. Rev. C73, 034913 (2006), arXiv:nucl-th/0508055 .
  * Cao and Bass (2011) S. Cao and S. A. Bass, Phys. Rev. C 84, 064902 (2011), arXiv:1108.5101 [nucl-th] .
  * Scardina  _et al._ (2017) F. Scardina, S. K. Das, V. Minissale, S. Plumari, and V. Greco, Phys. Rev. C 96, 044905 (2017), arXiv:1707.05452 [nucl-th] .
  * Chatterjee and Bożek (2018) S. Chatterjee and P. Bożek, Phys. Rev. Lett. 120, 192301 (2018), arXiv:1712.01189 [nucl-th] .
  * Schlichting and Pratt (2011) S. Schlichting and S. Pratt, Phys. Rev. C 83, 014913 (2011), arXiv:1009.4283 [nucl-th] .
  * Pratt  _et al._ (2015) S. Pratt, W. P. McCormack, and C. Ratti, Phys. Rev. C 92, 064905 (2015), arXiv:1508.07031 [nucl-th] .
  * Pratt and Plumberg (2021) S. Pratt and C. Plumberg, Phys. Rev. C 104, 014906 (2021), arXiv:2104.00628 [nucl-th] .
  * Pruneau  _et al._ (2023a) C. Pruneau, V. Gonzalez, B. Hanley, A. Marin, and S. Basu, Phys. Rev. C 107, 014902 (2023a).
  * Pruneau  _et al._ (2023b) C. Pruneau, V. Gonzalez, B. Hanley, A. Marin, and S. Basu, Phys. Rev. C 107, 054915 (2023b), arXiv:2211.10770 [hep-ph] .
  * Pratt  _et al._ (2022) S. Pratt, D. Oliinychenko, and C. Plumberg, Phys. Rev. C 106, 064911 (2022), arXiv:2210.03877 [nucl-th] .
  * Florkowski  _et al._ (2013) W. Florkowski, R. Ryblewski, and M. Strickland, Phys. Rev. C 88, 024903 (2013), arXiv:1305.7234 [nucl-th] .
  * Peterson  _et al._ (1983) C. Peterson, D. Schlatter, I. Schmitt, and P. M. Zerwas, Phys. Rev. D 27, 105 (1983).
  * Acharya  _et al._ (2022) S. Acharya _et al._ (ALICE), Phys. Rev. D 105, L011103 (2022), arXiv:2105.06335 [nucl-ex] .
  * Avramescu  _et al._ (2023) D. Avramescu, V. Băran, V. Greco, A. Ipp, D. I. Müller, and M. Ruggieri, Phys. Rev. D 107, 114021 (2023), arXiv:2303.05599 [hep-ph] .
  * Boguslavski  _et al._ (2023) K. Boguslavski, A. Kurkela, T. Lappi, F. Lindenbauer, and J. Peuron, (2023), arXiv:2303.12520 [hep-ph] .
  * Vredevoogd and Pratt (2009) J. Vredevoogd and S. Pratt, Phys. Rev. C79, 044915 (2009), arXiv:0810.4325 [nucl-th] .
  * Sun  _et al._ (2019) Y. Sun, G. Coci, S. K. Das, S. Plumari, M. Ruggieri, and V. Greco, Phys. Lett. B 798, 134933 (2019), arXiv:1902.06254 [nucl-th] .
  * Das  _et al._ (2017) S. K. Das, S. Plumari, S. Chatterjee, J. Alam, F. Scardina, and V. Greco, Phys. Lett. B768, 260 (2017), arXiv:1608.02231 [nucl-th] .


# Feasibility of Proof Of Authority as a Consensus Protocol Model

Shashank Joshi
Department of Computer Science and Engineering
SRM Institute Of Science And Technology
Kattankulathur, Tamil Nadu – 603203,
<EMAIL_ADDRESS>

###### Abstract

Blockchain is a type of decentralized distributed network which acts as an
immutable digital ledger. Despite the absence of any central governing
authority to validate the blocks in the ledger, it is considered secure and
immutable due to the consensus protocol among various nodes of the network. A
consensus algorithm is a mechanism that guarantees the reliability of the
blockchain and helps all connected nodes or peers to reach common ground
regarding the present state of the blockchain network thus an ideal consensus
algorithm must be secure, reliable, and fast. There are several different
algorithms to reach a consensus among the nodes thus this article seeks to
test the practicality of Proof of Authority in the blockchain network as a
consensus algorithm and its comparison with current mainstream consensus
algorithms.

_K_ eywords Consensus $\cdot$ Blockchain $\cdot$ Proof Of Authority $\cdot$
Proof of Work $\cdot$ Proof of stake

##  1 Introduction

The Concept of blockchain and the use of consensus algorithm was initially
proposed in 1982 [1], but a quantum leap was made with the introduction of
Bitcoin, a peer-to-peer electronic cash system [2] in 2009. The main idea
behind blockchain was to create a decentralized and growing list of blocks or
records which are integrated with the help of cryptography protocols. As
blockchain functions as a decentralized immutable ledger which records the
transaction in the real-time, it gives rise to a complexity to validate and
verify the transaction and make the system fault-tolerant, this can be
achieved with the help of a consensus algorithm which acts as an instrument
using which nodes in the network can coordinate and come to a common consensus
or single truth without a central authority, thus a consensus algorithm is an
essential element of a blockchain network that maintains the integrity and
security of these distributed computing systems. Many consensus algorithms
developed with the time, for instance initially classical consensus protocols
were the only known protocols that use all-to-all voting for obtaining single
truth, then Satoshi proposed longest-chain consensus protocols after that,
there were various variants of consensus protocols that came into existence
like the snow family of consensus protocols [3]. The primary goal of all these
developments in consensus protocols is to make the blockchain network more
secure and increase the transaction per second in the blockchain, as the
properties of an ideal consensus algorithm are:

  1. 1.

Scalability: A consensus algorithm should maintain a smooth operation of the
network by eliminating the possibility of slow processing time, system
bloating, lags, etc.

  2. 2.

Integrity: The result generated by all the non faulty nodes must be
indistinguishable.

  3. 3.

Termination: Every non faulty node should participate in deciding the output
[4].

  4. 4.

Cooperative: A consensus protocol should aim to achieve the maximum agreement
in the network [4].

  5. 5.

Egalitarian:A consensus protocol must give equal weightage to every vote in
the network [3].

As an immutable decentralized system blockchain’s nodes act as a cross-
platform runtime that provides services both as host and server and they are
required to share the data across all the nodes to reach a single truth to
carry out transactions. No permissions are required to be a node in a public
blockchain network thus there is no guarantee of integrity from the node thus
a node can alter the transaction data and hence it creates a fork in the
network. To avoid block tempering every node should agree on a single source
of truth thus a consensus algorithm is required[5]. Different sorts of
blockchains have different application scenarios. Therefore, blockchain
consensus algorithms should be appropriate and should suit the needs of
specific applications [6]. As the cryptocurrency space evolved many consensus
algorithms like PoW(Proof of Work ), PoS (Proof of Stake), etc. are developed
and become the mainstream consensus protocols, but there are some drawbacks of
these consensus algorithms thus to reduce these flaws PoA is a feasible
alternative consensus algorithm.

##  2 Proof of Authority

Proof of Authority is a regulated or permissioned consensus protocol, a family
of BFT algorithms [7] that aim to reach a consortium through authorized nodes
or validators. The validators are restricted to a fixed set of n authorized
nodes called sealers [8]. The working of this protocol can be summarized as:
the validating machines generate the new blocks of transaction which are
required to be added in the network. The mechanism of accepting and mining new
blocks in the network depends upon the configuration of the respective
blockchain network. In this protocol, the reputation of the sealer is at
stake, thus making malicious behaviour very damaging to the reputation of the
entity. Therefore, validating nodes of the blockchain network are responsible
for the security of the network. The whole architecture of PoA depends upon:

  1. 1.

Liveliness of the Network: As the framework of PoA depends on the reputation
of the sealer hence the liveliness of the network implies that the reputation
system of the network does not provide a possibility to reach a state at which
the bad behaviour of the sealer no longer damages their reputation [9].

  2. 2.

Tolerance: In a blockchain network with n sealers a maximum of p nodes where

$p<\frac{n}{2}$ (1)

can be dishonest or byzantine. Therefore to solve this byzantine consensus
problem [10] with n nodes where $n-p$ nodes are required to be honest and they
all should come to a single source of truth i.e. a unique block in the
network.

  3. 3.

Monotonicity: The reputation system associated with the PoA framework should
not be monotonic in nature in any one direction as it will decrease the trust
factor of the network [9].

  4. 4.

Visibility: The sealers in the network are required to prove and validate
their actual identity and the reputation stats of all the sealers must be
global and open to all the nodes in the blockchain network.

  5. 5.

Difficulty in becoming a sealer: The process of becoming a sealer should be
very selective and rigorous so that we can lower down any potential threat to
the network [3].The process of selection should be transparent and equal for
all.

In a PoA based blockchain network, the validating node or sealers should be
incentivised to retain their reputation or gain subsequently as this will
refrain sealers from associating with any malicious activity on the network.
The PoA consensus depends upon the mining rotation schema [11] i.e. the
responsibility of mining new blocks and keeping the network safe is fully
dependent and distributed among the validatory nodes, hence the action of a
sealer is strictly monitored by other sealers in the network.

###  2.1 Comparison With Proof of Work

Proof of Work (PoW) is another consensus model used in various blockchain
networks and cryptocurrencies like bitcoin etc [12]. According to this
consensus model, an amount of work or computing power is required to mine new
blocks in the blockchain network. This mechanism is designed to prevent any
malicious activity in the network and make it less susceptible to a potential
attack on the network. In PoW based protocols the mining nodes are required to
compete and solve an unknown complex mathematical challenge and obtain a
result less than the target value before adding a new block in the network.
This computational intensive approach is the core safety mechanism in this
protocol. As PoW is a computationally intensive approach it slows down the tps
(transaction per second) of the network and the time lag in accepting the new
block may lead another miner to find the optimum solution as a proof of work
thus creating a fork in the network which increase the time taken to reach on
a consensus [6]. PoA on the other hand is less resource-consuming and has a
far better throughput as compared to PoW due to the limited number of
validatory nodes in the network thus improving the scalability of the system.
PoA also provides an easy development and maintenance of decentralized
applications and a centralized fashion decreases the number of forks,
therefore, scaling down the number of possible attacks on the network.

###  2.2 Comparison With Proof of Stake

Proof of Stake (PoS) is a type of consensus model which was initially
developed as an energy-saving alternative to PoW. In PoS, the leaders are
selected based on their stakes in the blockchain network, i.e., the number of
digital tokens owned by the node in the network. In this model, the selection
of the validatory node to add new blocks in the network is pseudorandom in
nature. In this consensus model, instead of solving a complex mathematical
problem, users who want to add new blocks in the network must lock a share of
their stakes in the network. The probability $p_{i}$ of node i to be selected
as a validatory node or leader among N participatory nodes is given by:

$p_{i}=\frac{s_{i}}{\sum_{j=1}^{N}s_{j}}$ (2)

(Where $s_{i}$ is the stake of i th node) [13]. PoA in contrast is a
reputation-based consensus model that relies on a limited number of validating
nodes for adding new blocks in the network. PoA gives more importance to the
reputation of validating nodes rather than on their stakes in the network.
Thus, improving the scalability and throughput of the network. PoA also
presents a better solution for a permissioned blockchain network as compared
to PoS.

##  3 Advantages

  * •

PoA is less resource-consuming as compared to other consensus models as there
is no direct competition among the nodes to add new blocks because of the
limited number of validatory nodes in the network.

  * •

This consensus model can achieve more efficiency and throughput relative to
other consensus models due to its centralized reputation-based approach to
reach a consensus.

  * •

PoA based consensus model also solves the problem of market-based fluctuating
processing fees by maintaining it constant throughout the network.

  * •

Byzantine fault tolerance of PoA based consensus models is better than non-
distributed protocols of a typical centralized system [11].

  * •

Because of less number of forks generated in the network, the PoA based
consensus protocols have better transaction rate as compared to PoS or PoW.

  * •

The PoA based consensus model is well suited for the maintenance and
development of DApps and is ascendable as compared to other consensus
models[14].

##  4 Limitations

  * •

The PoA consensus model relinquishes some degree of decentralization because
of its mechanism to reach a consensus using validatory nodes thus we can
conclude PoA sacrifices decentralization in order to achieve high throughput
and scalability [7].

  * •

The identity of validating nodes is public across the network, thus knowing
the validators’ identity can cause third-party manipulation and interference
in the blockchain network.

  * •

PoA consensus models’ centralized approach introduces censorship and
blacklisting in the blockchain network.

  * •

Byzantine fault tolerance of PoA based consensus models is better than non-
distributed protocols of a typical centralized system [11].

  * •

PoA is a reputation-based consensus model, but if a validatory node in the
network is intended to do malicious activities in the network, the threat of
harming the reputation does not shield the network from malicious activities
and third-party involvement.

  * •

The PoA consensus model is prone to Sybil and cloning attacks [15] in which
one or two sealers can manipulate a large number of validatory nodes by
forging multiple identities.

##  5 Conclusion

In this paper, there is a detailed analysis of the feasibility of proof of
authority as a consensus algorithm in a blockchain network. The inherent
features of PoA make it a potential alternative to overcome the shortcoming of
PoW and PoS consensus protocols. PoA, as a consensus model, sacrifices
decentralization to achieve high throughput and scalability, which makes it a
suitable option for the DApps development and maintenance, private
blockchains, and notary-based applications which require users to trust
validators. On the other hand, it is still unproven for the public blockchain
networks like cryptocurrencies, which give more importance to
decentralization. Still, the cost-efficiency and performance of this construct
make it a potential alternative for mainstream consensus protocols and an
emerging blockchain solution for private or permissioned blockchain networks.

## References

  * [1] David L. Chaum.  Computer Systems Established, Maintained, and Trusted by Mutually Suspicious Groups.  PhD thesis, University Of California, 1982.
  * [2] Satoshi Nakamoto.  Bitcoin: A peer-to-peer electronic cash system.  Cryptography Mailing list at https://metzdowd.com, 03 2009.
  * [3] Sivleen Kaur, Sheetal Chaturvedi, Aabha Sharma, and Jayaprakash Kar.  A research survey on applications of consensus protocols in blockchain.  Security and Communication Networks, 2021:1–22, 2021.
  * [4] Yang Xiao, Ning Zhang, Jin Li, Wenjing Lou, and Y. Hou.  Distributed Consensus Protocols and Algorithms, pages 25–50.  03 2019.
  * [5] Xiaoqi Li, Peng Jiang, Ting Chen, Xiapu Luo, and Qiaoyan Wen.  A survey on the security of blockchain systems.  Future Generation Computer Systems, 107, 08 2017.
  * [6] Arvind Narayanan, Joseph Bonneau, Edward Felten, Andrew Miller, and Steven Goldfeder.  Bitcoin and cryptocurrency technologies.
  * [7] Xuefeng Liu, Gansen Zhao, Xinming Wang, Yixing Lin, Ziheng Zhou, Hua Tang, and Bingchuan Chen.  Mdp-based quantitative analysis framework for proof of authority.  pages 227–236, 10 2019.
  * [8] Parinya Ekparinya, Vincent Gramoli, and Guillaume Jourjon.  The attack of the clones against proof-of-authority, 02 2019.
  * [9] Omar Hasan, L. Brunie, and E. Bertino.  Privacy preserving reputation systems based on blockchain and other cryptographic building blocks: A survey.  2020\.
  * [10] M. Pease, R. Shostak, and L. Lamport.  Reaching agreement in the presence of faults.  Journal of the ACM, 27(2):228–234, 1980.
  * [11] Stefano De Angelis, Leonardo Aniello, Federico Lombardi, Andrea Margheri, and V. Sassone.  Pbft vs proof-of-authority: applying the cap theorem to permissioned blockchain.  01 2017.
  * [12] Andreas M Antonopoulos.  Mastering Bitcoin.  O’Reilly Media, Incorporated, 2017.
  * [13] Cong Nguyen, Hoang Dinh Thai, Diep Nguyen, Dusit Niyato, Huynh Nguyen, and Eryk Dutkiewicz.  Proof-of-stake consensus mechanisms for future blockchain networks: Fundamentals, applications and opportunities.  IEEE Access, PP:1–1, 06 2019.
  * [14] S. D. Angelis.  Assessing security and performances of consensus algorithms for permissioned blockchains.  ArXiv, abs/1805.03490, 2018.
  * [15] John R. Douceur.  The sybil attack.  In IPTPS, 2002.