Convergence properties of fixed-point search

with general but equal phase shifts for any number of iterations

111The paper was supported by NSFC(Grants No. 60433050 and 60673034), the
basic research fund of Tsinghua university NO: JC2003043.

Dafa Lia222email<EMAIL_ADDRESS>Xiangrong Lib, Hongtao
Huangc, Xinxin Lid

a Department of mathematical sciences, Tsinghua University, Beijing 100084
CHINA

b Department of Mathematics, University of California, Irvine, CA 92697-3875,
USA

c Electrical Engineering and Computer Science Department

University of Michigan, Ann Arbor, MI 48109, USA

d Department of computer science, Wayne State University, Detroit, MI 48202,
USA

The correspondence Author Dafa Li,

Phone Number is (8610)62773561

Fax No. is (8610) 62785847

Abstract

Grover presented the fixed-point search by replacing the selective inversions
by selective phase shifts of $\pi/3$. In this paper, we investigate the
convergence behavior of the fixed-point search algorithm with general but
equal phase shifts for any number of iterations.

PACS number: 03.67.Lx

Keywords: Amplitude amplification, the fixed-point search, quantum computing.

##  1 Introduction

Grover’s quantum search algorithm is used to find a target state in an
unsorted database of size $N$[1][2]. The Grover’s quantum search algorithm can
be considered as a rotation of the state vectors in two-dimensional Hilbert
space generated by the start ($s$) and target ($t$) vectors[2]. The amplitude
of the target state increases monotonically towards its maximum and decreases
monotonically after reaching the maximum [3]. This search algorithm is called
the amplitude amplification algorithm. For the size $N=2^{n}$ of the database,
quantum search algorithm requires $O(\sqrt{N})$ steps to find the target
state. As mentioned in [4] [5], unless we stop when it is right at the target
state, it will drift away. A fixed-point search algorithm was presented in [4]
to avoid drifting away from the target state. The fixed-point search algorithm
obtained by replacing the selective inversions by selective phase shifts of
$\pi/3$, converges to the target state irrespective of the number of
iterations. The main advantage of the fixed-point search with equal phase
shifts of $\pi/3$ is that it performs well for small but unknown initial error
probability and the fixed-point behavior leads to robust quantum search
algorithms [4]. However, the target state is the limit state when the number
of iterations tends to the infinite.

For readability, we introduce the fixed-point search algorithm as follows. In
[4] the transformation $UR_{s}^{\pi/3}U^{+}R_{t}^{\pi/3}U$, where $U$ is any
unitary operator, was applied to the start state $|s\rangle$,

$\displaystyle R_{s}^{\pi/3}$ $\displaystyle=$ $\displaystyle
I-[1-e^{i\frac{\pi}{3}}]|s\rangle\langle s|,$ $\displaystyle R_{t}^{\pi/3}$
$\displaystyle=$ $\displaystyle I-[1-e^{i\frac{\pi}{3}}]|t\rangle\langle t|,$
(1)

where $|t\rangle$ stands for the target state. The transformation
$UR_{s}^{\pi/3}U^{+}R_{t}^{\pi/3}U$ is denoted as Grover’s the Phase-$\pi/3$
search algorithm in [6].

Let us consider the fixed-point search algorithm with general but equal phase
shifts as follows.

$\displaystyle R_{s}^{\theta}$ $\displaystyle=$ $\displaystyle
I-[1-e^{i\theta}]|s\rangle\langle s|,$ $\displaystyle R_{t}^{\theta}$
$\displaystyle=$ $\displaystyle I-[1-e^{i\theta}]|t\rangle\langle t|.$ (2)

The transformation $UR_{s}^{\theta}U^{+}R_{t}^{\theta}U$ was called as the
Phase-$\theta$ search algorithm and studied in [7]. It is enough to let
$\theta$ be in $[0,\pi]$.

Note that if we apply $U$ to the start state $|s\rangle$, then the amplitude
of reaching the target state $|t\rangle$ is $U_{ts}$[2], where
$\left|\left|U_{ts}\right|\right|^{2}=1-\epsilon$. As indicated in [2], in the
case of database search, $|U_{ts}|$ is almost $1/\sqrt{N}$, where $N$ is the
size of the database. Thus, $\epsilon$ is almost $1-1/N$ and $\epsilon$ is
close to $1$ for the large size of database.

Apply the operations $U$, $R_{s}^{\theta}$, $U^{+}$, $R_{t}^{\theta}$, and $U$
to the start $|s\rangle$ and let $D(\theta)$ be the deviation of the state
$UR_{s}^{\theta}U^{+}R_{t}^{\theta}U|s\rangle$ from the $t$ state for any
phase shifts of $\theta$. The deviation $D(\theta)$ was reduced in [7] and is
rewritten as follows.

$D(\theta)=4(1-\cos\theta)^{2}\epsilon(\epsilon-d)^{2},$ (3)

where $d=\frac{1-2\cos\theta}{2(1-\cos\theta)}$. It was shown that $D(\theta)$
is between $0$ and $1$ in [7]. For the Phase-$\pi/3$ search algorithm,
$D(\pi/3)=\epsilon^{3}$[4].

In [8], we explored the performance of the fixed-point search with general but
different phase shifts for one iteration. In [7], we discussed the performance
of the fixed-point search with general but equal phase shifts for one
iteration.

In this paper, we investigate convergence behavior of the fixed-point search
with general but equal phase shifts for any number of iterations. It is useful
for designing fixed-point search algorithms for different choices of the phase
shift parameter $\theta$. The following results are established in Section 2.

(1). The fixed-point search with equal phase shifts of $\theta\leq\pi/2$
converges to the target state.

(2). The fixed-point search with equal phase shifts of $\theta$, where
$\pi/2<\theta\leq\arccos(-1/4)$, converges the target state with the
probability of at least $80\%.$

(3). The fixed-point search with equal phase shifts of $\theta$, where
$\arccos(-1/4)<\theta\leq 2\pi/3$, converges the target state with the
probability of among $66.6\%$ and $80\%.$

(4). The fixed-point search with equal phase shifts of $\theta$, where
$2\pi/3<\theta\leq\pi$, does not converge.

In section 3, we analyze the convergence rate for different values of
$\theta$. It is demonstrated that the Phase-$\pi/3$ is not always optimal and
the convergence rate can be improved by choosing $\theta>\pi/3$. In section 4,
we show that for the size $N=2^{n}$ of the database, $O(n)$ iterations of the
Phase-$\theta$ search can find the target state. However, as indicated in [4],
$O(n)$ iterations of the Phase-$\theta$ search involve the exponential
queries.

##  2 Convergence performance of the Phase-$\theta$ search for any number of
iterations

Let $\epsilon_{0}=\epsilon$ and $0<\epsilon<1$. Then, from Eq. (3) one can
obtain the following iteration equation

$\epsilon_{m+1}=4(1-\cos\theta)^{2}\epsilon_{m}(\epsilon_{m}-d)^{2}.$ (4)

In this section, we discuss the convergence behavior of the Phase-$\theta$
search for any number of iterations. For the Phase-$\pi/3$ search, after
recursive application of the basic iteration for $m$ times, the failure
probability $\epsilon_{m}=$ $\epsilon^{3^{m}}$ and the success probability
$\left|U_{m,ts}\right|=1-\epsilon^{3^{m}}$[4]. The $\epsilon_{m}$ in Eq. (4)
is the failure probability of the Phase-$\theta$ search algorithm after $m$
iterations.

From inference [11] in [7], Eq. (4) has the following fixed-points: $0$, $1$
($\theta\neq 0$), $a$, where $a=\cos\theta/(\cos\theta-1)$ ($\theta\neq 0$).
In other words, if the sequence $\\{\epsilon_{m}\\}$ in Eq. (4) has a limit
then the limit must be $0$, $1$ or $a$. Clearly $a<d$.

To study the convergence performance for any number of iterations, we need the
following results which are listed in the following paragraphs (A), (B), and
(C).

(A). From Eq. (4), we obtain the following,

$\epsilon_{m}-\epsilon_{m-1}=4\epsilon_{m-1}\left(\cos\theta-1\right)^{2}\left(1-\epsilon_{m-1}\right)(a-\epsilon_{m-1}).$
(5)

Eqs. (4) and (5) imply the following convergence property.

Property 1.

(1.1) If $\epsilon_{m}=d$, then $\epsilon_{m+l}=0$, for any $l>0$.

(1.2). if $\epsilon_{m-1}>a$ and $\epsilon_{m-1}\neq 0$,
$\epsilon_{m}<\epsilon_{m-1}$;

(1.3). If $\epsilon_{m-1}<a$ and $\epsilon_{m-1}\neq 0$,
$\epsilon_{m}>\epsilon_{m-1}$.

(B). When $\pi/2\leq\theta\leq\pi$, we have the following equation.

$\epsilon_{m}-a=4(\cos\theta-1)^{2}(\epsilon_{m-1}-a)(\epsilon_{m-1}-b)(\epsilon_{m-1}-c),$
(6)

where
$b=\frac{1}{2}-\frac{\sqrt{-\cos\theta(2-\cos\theta)}}{2(1-\cos\theta)}$, and
$c=\frac{1}{2}+\frac{\sqrt{-\cos\theta(2-\cos\theta)}}{2(1-\cos\theta)}$. When
$\epsilon_{i}=a$, $b$ or $c$, $\epsilon_{i+l}=a$ for any $l>0$. Note that
$b<d$, $a<d$, and $d<c$.

(C). Let

$f(x)=4(1-\cos\theta)^{2}x(x-d)^{2}.$ (7)

Then, the derivative of $f(x)$ is

$f^{\prime}(x)=12(1-\cos\theta)^{2}(x-d)(x-d/3).$ (8)

From Eq. (8), (1). $f^{\prime}(x)=0$ at $r$ and $d$, where $r=d/3$; (2).$\
f^{\prime}(x)<0$ when $r<x<d$; (3). $f^{\prime}(x)>0$ when $x<r$ or $x>d$;
(4). When $\theta>\pi/3$, $f(x)$ has a relative maximum
$g=\frac{2(1-2\cos\theta)^{3}}{27(1-\cos\theta)}$ at $r$ and a relative
minimum $0$ at $d$.

###  2.1 When $0<\theta\leq\pi/2$, for any $\epsilon_{0}\in(0,1)$, the
Phase-$\theta$ search converges to the target state.

Note that $0$ is an attractive fixed-point when $0<\theta<\pi/2$ and $0$ is
also a semi-attractive fixed-point when $\theta=\pi/2$. See inference [11] in
[7].

(1). $0<\theta\leq\pi/3$

For this case, $d\leq 0$ and $a<0$. In Eq. (4), $d=0$ means
$\epsilon_{m+1}=\epsilon_{m}^{3}$, which is Grover’s Phase-$\pi/3$ search.
From $d<0$ and Eq. (4), $\epsilon_{m}>0$. By property (1.2), always
$\epsilon_{m}<\epsilon_{m-1}$ when $0<\theta\leq\pi/3$. That is, the sequence
$\\{\epsilon_{m}\\}$ in Eq. (4) decreases monotonically.Therefore, for any
$\epsilon_{0}$ in $(0$, $1)$ $\lim_{m\rightarrow\infty}\epsilon_{m}=0$.

(2). $\pi/3<\theta\leq\pi/2$

For this case, $a\leq 0$, $0<d\leq 1/2$. Hence, from Eq. (4)
$0\leq\epsilon_{i}<1$. By property (1.2), always
$\epsilon_{m}\leq\epsilon_{m-1}$ when $\pi/3<\theta\leq\pi/2$. That is, the
sequence $\\{\epsilon_{m}\\}$ in Eq. (4) decreases. Factually, the sequence
$\\{\epsilon_{m}\\}$ in Eq. (4) decreases monotonically and $\epsilon_{m}>0$,
or is of the form $\epsilon_{0}>\epsilon_{1}>...>\epsilon_{k}=0$ and
$\epsilon_{l}=0$ for any $l>k$. Therefore, for any $\epsilon_{0}$ in $(0$,
$1)$ $\lim_{m\rightarrow\infty}\epsilon_{m}=0$.

Example 1. For the Phase-$\pi/2$ search,
$\epsilon_{m}=\epsilon_{m-1}(2\epsilon_{m-1}-1)^{2}$. Let
$\epsilon_{0}=0.99999$. See Fig. 1.

$\epsilon_{1}=\allowbreak 0.999\,95$, $\ \ \epsilon_{2}=\allowbreak
0.999\,75$, $\epsilon_{3}=\allowbreak 0.998\,75$, $\ \
\epsilon_{4}=\allowbreak 0.993\,76$,

$\epsilon_{5}=\allowbreak 0.969\,11$, $\ \ \epsilon_{6}=\allowbreak
0.853\,07$, $\epsilon_{7}=\allowbreak 0.425\,37$, $\ \
\epsilon_{8}=\allowbreak 9.\,\allowbreak 476\,6\times 10^{-3}$.

###  2.2 When $\pi/2<\theta<\arccos(-1/4)$, the Phase-$\theta$ search
converges the target state with the success probability of $(1-a)>80\%$.

For the Phase-$\theta$ search, $a<g<r<b<d<c$. Note that $a$ is an attractive
fixed-point. See inference [11] in [7]. From Eq. (6), we have the following
property.

Property 2.

(2.1). $a<\epsilon_{m}\leq g$ whenever $a<\epsilon_{m-1}<b$;

(2.2). $0\leq\epsilon_{m}<a$ whenever $b<\epsilon_{m-1}<c$ or
$\epsilon_{m-1}<a$.

The convergence region of the Phase-$\theta$ search

(A). When $\epsilon_{0}\in(0,c]$ and $\epsilon_{0}\neq d$, the deviation from
the target state converges to the fixed-point $a$.

There are four cases. The argument is the following.

Case 1. When $\epsilon_{0}=a$ or $b$ or $c$, it is trivial by Eq. (6).

Case 2. $\epsilon_{0}<a$. By property (2.2), $0<\epsilon_{m}<a$ for any $m$.
By property 1, the sequence $\\{\epsilon_{m}\\}$ increases monotonically.
Hence, the sequence $\\{\epsilon_{m}\\}$ converges to $a$ from below.

Case 3. $a<\epsilon_{0}<b$. By property (2.1), always $a<\epsilon_{m}\leq g$
for any $m>0$, and by property 1, the sequence $\\{\epsilon_{m}\\}$ decreases
monotonically. Hence, the sequence $\\{\epsilon_{m}\\}$ converges to $a$ from
above.

Case 4. $b<\epsilon_{0}<c$ and $\epsilon_{0}\neq d$. By property (2.2),
$0<\epsilon_{1}<a$. Then, it turns to case 2.

Conclusively, when $\epsilon_{0}\in(0,c]$ and $\epsilon_{0}\neq d$, from the
above four cases, $\epsilon_{m}\neq d$, hence
$\lim_{m\rightarrow\infty}\epsilon_{m}=a$.

(B). When $\epsilon_{0}\in(c,1)$, the deviation from the target state
converges to the fixed-points $a$ or $0$.

By property (1.2),
$\epsilon_{0}>...>\epsilon_{j^{\ast}-1}(>c)>\epsilon_{j^{\ast}\text{ }}(\leq
c)$. If $\epsilon_{j^{\ast}\text{ }}=d$, then $\epsilon_{m}=0$ for any
$m>j^{\ast}$. Otherwise, $\lim_{m\rightarrow\infty}\epsilon_{m}=a$ by the
above (A).

###  2.3 Phase-$\arccos(-1/4)$ search converges the target state with the
success probability of $80\%$.

For the Phase-$\arccos(-1/4)$ search, $a=1/5$ is an attractive fixed-point,
see inference [11] in [7]. $b=a=1/5$, $d=3/5$, and $c=4/5$. The iteration
equation is $\epsilon_{m}=\epsilon_{m-1}(\allowbreak
5\epsilon_{m-1}-3)^{2}/4$. Eq. (6) becomes the following.

$\epsilon_{m}-1/5=\allowbreak\frac{25}{4}\left(\epsilon_{m-1}-4/5\right)\left(\epsilon_{m-1}-1/5\right)^{2}.$
(9)

From Eq. (9) we have the following property.

Property 3.

(3.1). $\epsilon_{m}<1/5$ when $\epsilon_{m-1}<4/5$ and $\epsilon_{m-1}\neq
1/5$.

(3.2). $\epsilon_{m}>1/5$ when $\epsilon_{m-1}>4/5$.

The convergence region of the Phase-$\arccos(-1/4)$ search

(A). When $\epsilon_{0}\in(0,4/5]$ and $\epsilon_{0}\neq 3/5$, the deviation
from the target state converges to the fixed-point $1/5$.

When $\epsilon_{0}=1/5$ or $4/5$, it is trivial by Eq. (6). When
$\epsilon_{0}\in(0,4/5)$ and $\epsilon_{0}\neq 1/5$, always $\epsilon_{m}<1/5$
for $m>0$ by property (3.1) and the sequence $\\{\epsilon_{m}\\}$ increases
monotonically from $m>0$ by property (1.3). Therefore, the sequence
$\\{\epsilon_{m}\\}$ converges to $1/5$ from below.

(B) When $\epsilon_{0}\in(4/5,1)$, the deviation from the target state
converges to the fixed-points $1/5$ or $0$.

By property (1.2), $\epsilon_{0}>\epsilon_{1}>....>\epsilon_{m}(\leq 4/5)$.
Case 1. If $\epsilon_{m}=3/5$, then $\epsilon_{i}=0$ for any $i>m$. Case 2\.
Otherwise, by the above (A), $\lim_{m\rightarrow\infty}\epsilon_{m}=1/5$.

Example 2. Let $\epsilon_{0}=0.9999;$

$\epsilon_{1}=\allowbreak 0.999\,4$, $\ \ \epsilon_{2}=\allowbreak 0.996\,4$,
$\ \ \ \ \ \ \ \ \epsilon_{3}=\allowbreak 0.978\,55$, $\
\epsilon_{4}=\allowbreak 0.876\,41$,

$\epsilon_{5}=\allowbreak 0.418\,50$, $\ \epsilon_{6}=\allowbreak
8.\,\allowbreak 616\,5\times 10^{-2}$, $\ \epsilon_{7}=\allowbreak 0.142\,19$,
$\ \epsilon_{8}=\allowbreak 0.186\,26$,

$\epsilon_{9}=\allowbreak 0.199\,28$, $\ \epsilon_{10}=\allowbreak 0.2$.

###  2.4 When $\arccos(-1/4)<\theta\leq 2\pi/3$, the Phase-$\theta$ search
converges the target state with the success probability of $(1-a)$, where
$66\%\leq(1-a)<80\%$.

For the Phase-$\theta$ search, $b<r<a<g<d<c$. Note that $a$ is an attractive
fixed-point when $\arccos(-1/4)<\theta<2\pi/3$ and $1/3$ is a semi-attractive
fixed-point when $\theta=2\pi/3$. See inference [11] in [7]. From Eq. (6), we
have the following property.

Property 4.

(4.1). $a<\epsilon_{m}\leq g$ whenever $b<\epsilon_{m-1}<a$;

(4.2). $0\leq\epsilon_{m}<a$ whenever $a<\epsilon_{m-1}<c$ or
$\epsilon_{m-1}<b$.

The convergence region of the Phase-$\theta$ search

(A). When $\epsilon_{0}\in(0,c]$ and $\epsilon_{0}\neq d$, the deviation from
the target state converges to the fixed-point $a$.

There are seven cases. We argue them as follows.

Case 1. If $\epsilon_{0}=a$ or $b$ or $c$, then it is trivial.

Case 2. $a<\epsilon_{0}\leq g$. The proof is put in Appendix A.

Case 3. $\epsilon_{0}<b$. By property (1.3), $\epsilon_{j}$ increases
monotonically from $\epsilon_{0}$ until $\epsilon_{j^{\ast}-1}<b$ and
$b\leq\epsilon_{j^{\ast}}<f(b)=a$ since $f^{\prime}(x)>0$ when $x<b$. If
$\epsilon_{j^{\ast}}=b$, it is trivial. Otherwise, by property (4.1),
$\epsilon_{j^{\ast}+1}$ is in $(a,g]$. Now it turns to case 2.

Case 4. $b<\epsilon_{0}\leq r$. When $\epsilon_{0}=r$,
$\epsilon_{m}=f^{(m-1)}(g)$. From the proof of case 2,
$\lim_{m\rightarrow\infty}f^{(m-1)}(g)=a$. Next consider that
$b<\epsilon_{0}<r$. Since $f^{\prime}(x)>0$ when $b<x<r$,
$f(b)<f(\epsilon_{0})<f(r)$. That is, $a<\epsilon_{1}<g$. It turns to case 2.

Case 5. $r<\epsilon_{0}<a$. Since $f^{\prime}(x)<0$ when $r<x<a$,
$a<\epsilon_{1}<g$. It turns to case 2.

Case 6. $g<\epsilon_{0}<d$. Since $f^{\prime}(x)<0$ when $g<x<d$ and $a<g$,
$0<\epsilon_{1}<f(g)<a$. Then, it turns to cases 1, 3, 4, 5.

Case 7. $d<\epsilon_{0}<c$. Since $f^{\prime}(x)>0$ when $d<x<c$,
$0<\epsilon_{1}<a$. Then, it turns to cases 1, 3, 4, 5.

(B). When $\epsilon_{0}\in(c,1)$, the deviation from the target state
converges to the fixed-points $a$ or $0$.

When $\epsilon_{0}>c$, by property (1.2) the sequence $\\{\epsilon_{i}\\}$
decreases monotonically from $\epsilon_{0}$ to $\epsilon_{i^{\ast}\text{
}}\leq c$. Case 1, if $\epsilon_{i^{\ast}\text{ }}=d$, then $\epsilon_{i}=0$
for any $i>i^{\ast}$. Case 2. Otherwise, by the above (A)
$\lim_{m\rightarrow\infty}\epsilon_{m}=a$.

Example 3. For the Phase-$2\pi/3$ search, $a=1/3$. The iteration equation
becomes $\epsilon_{m}=\epsilon_{m-1}(3\epsilon_{m-1}-2)^{2}$. Let
$\epsilon_{0}=0.99999$. We have the following iterations. See Fig. 1.

$\epsilon_{1}=0.999\,93$, $\ \epsilon_{2}=0.999\,51$,
$\epsilon_{3}=0.996\,57$, $\ \epsilon_{4}=0.976\,17$,

$\epsilon_{5}=0.841\,59$, $\ \epsilon_{6}=0.231\,76$, $\ \
\epsilon_{7}=\allowbreak 0.394\,52$, $\ \epsilon_{8}=\allowbreak 0.263\,50$,

$\epsilon_{9}=\allowbreak 0.385\,47$, $\ \epsilon_{10}=\allowbreak 0.274\,32$,
$\epsilon_{11}=\allowbreak 0.380\,05$, $\ \epsilon_{12}=\allowbreak
0.280\,99$,

$\epsilon_{13}=\allowbreak 0.376\,17$.

###  2.5 $2\pi/3<\theta\leq\pi$, the Phase-$\theta$ search does not converge.

For the Phase-$\theta$ search, $b<r<a<d<c$. From Eq. (6), we have the
following property.

Property 5

(5.1). $a<\epsilon_{m}\leq g$ whenever $b<\epsilon_{m-1}<a$;

(5.2). $0\leq\epsilon_{m}<a$ whenever $a<\epsilon_{m-1}<c$ or
$\epsilon_{m-1}<b$.

For large $\epsilon$, by property (1.2), the sequence $\\{\epsilon_{i}\\}$
decreases monotonically from $\epsilon_{0}$ to $\epsilon_{i^{\ast}}(\leq c)$.
If $\epsilon_{i^{\ast}}=d$, then $\epsilon_{i}=0$ for any $i>i^{\ast}$. If
$\epsilon_{i^{\ast}}=a$, $b$, or $c$, then $\epsilon_{i}=a$ when $i>i^{\ast}$.
Otherwise, when $i>i^{\ast}$, $\epsilon_{i}$ oscillate around the fixed point
$a$ by property 1. However, the sequence $\\{\epsilon_{i}\\}$ does not
converges because $a$, $0$ and $1$ are repulsive fixed-points.

Example 4. For the Phase-$\pi$ search, the iteration equation becomes
$\epsilon_{m}=\epsilon_{m-1}(4\epsilon_{m-1}-3)^{2}$, $a=1/2$. Let
$\epsilon=0.99999$. We have the following iterations. See Fig. 1.

$\epsilon_{1}==0.999\,91$, $\ \ \epsilon_{2}=0.999\,19,$ $\
\epsilon_{3}=0.992\,73$, $\ \epsilon_{4}=0.935\,83,$

$\epsilon_{5}=0.517\,07$, $\ \ \ \epsilon_{6}=0.448\,87$,
$\epsilon_{7}=0.651\,25$, $\ \epsilon_{8}=0.101\,61,$

$\epsilon_{9}=0.683\,49$, $\ \ \ \epsilon_{10}=4.\,837\,6\times 10^{-2}$.

Clearly, the sequence $\\{\epsilon_{i}\\}$ monotonically decreases from
$\epsilon_{0}$ to $\epsilon_{6}$. Note that after the sixth iteration,
$\epsilon_{m}$ oscillate around the fixed point $1/2$.

##  3 A comparison of rates of convergence after any number of iterations

For the Phase-$\pi/3$ search, let the iteration equation be
$\epsilon_{m}(\pi/3)=(\epsilon_{m-1}(\pi/3))^{3}$, where $\epsilon_{m}(\pi/3)$
is the failure probability of the Phase-$\pi/3$ search algorithm after $m$
iterations. For the Phase-$\theta$ search, we can rewrite Eq. (4) as
$\epsilon_{m}(\theta)=4(1-\cos\theta)^{2}\epsilon_{m-1}(\theta)(\epsilon_{m-1}(\theta)-d)^{2}$,
where the $\epsilon_{m}(\theta)$ is the failure probability of the
Phase-$\theta$ search algorithm after $m$ iterations. We want to compare the
failure probability of the Phase-$\theta$ ($\neq\pi/3$) search algorithm with
the one of the Phase-$\pi/3$ search after $m$ iterations. It is known that the
less the failure probability is, the faster the algorithm converges. By
factoring,

$\displaystyle\epsilon_{m}(\theta)-\epsilon_{m}(\pi/3)$ $\displaystyle=$ (10)
$\displaystyle\epsilon_{m-1}(\theta)(2\cos\theta-1)(1-\epsilon_{m-1}(\theta))(3-2\cos\theta)\ast$
$\displaystyle(\epsilon_{m-1}(\theta)-\frac{1-2\cos\theta}{3-2\cos\theta})+\epsilon_{m-1}^{3}(\theta)-\epsilon_{m-1}^{3}(\pi/3).$

We have the following results.

(1). $\pi/3<\theta\leq\pi$

Case 1. For large $\epsilon$, the Phase-$\theta$ search converges faster than
the Phase-$\pi/3$ search for $m$ iterations until $\epsilon_{m-1}(\theta)<$
$\frac{1-2\cos\theta}{3-2\cos\theta}$.

In [7], we show if
$\epsilon_{0}(\theta)=\epsilon_{0}(\pi/3)=\epsilon>\frac{1-2\cos\theta}{3-2\cos\theta}$
then $\epsilon_{1}(\theta)<\epsilon_{1}(\pi/3)=\epsilon^{3}$. If
$\epsilon_{m-1}(\theta)>\frac{1-2\cos\theta}{3-2\cos\theta}$ and
$\epsilon_{m-1}(\theta)<\epsilon_{m-1}(\pi/3)$, then by Eq. (10)$\
\epsilon_{m}(\theta)<\epsilon_{m}(\pi/3)$. Thus,
$\epsilon_{i}(\theta)<\epsilon_{i}(\pi/3)$, where $i=1$, $2$, …, $m-1$, until
$\epsilon_{m-1}(\theta)<$ $\frac{1-2\cos\theta}{3-2\cos\theta}$. It says that
after $m$ iterations, the failure probability of the Phase-$\theta$ search is
less than the one of the Phase-$\pi/3$ search until $\epsilon_{m-1}(\theta)<$
$\frac{1-2\cos\theta}{3-2\cos\theta}$. It suggests us first to use the fixed-
point search with large phase shifts for the large size of database.

Case 2. For small $\epsilon$, the Phase-$\pi/3$ search converges faster than
the Phase-$\theta$ search for $m$ iterations until
$\epsilon_{m-1}(\theta)>\frac{1-2\cos\theta}{3-2\cos\theta}$.

In [7], we show if
$\epsilon_{0}(\theta)=\epsilon_{0}(\pi/3)=\epsilon<\frac{1-2\cos\theta}{3-2\cos\theta}$
then $\epsilon_{1}(\theta)>\epsilon_{1}(\pi/3)=\epsilon^{3}$. If
$\epsilon_{m-1}(\theta)<\frac{1-2\cos\theta}{3-2\cos\theta}$ and
$\epsilon_{m-1}(\theta)>\epsilon_{m-1}(\pi/3)$, then by Eq. (10)$\
\epsilon_{m}(\theta)>\epsilon_{m}(\pi/3)$.

(2). When $0<\theta<\pi/3$, the Phase-$\pi/3$ search converges faster than the
Phase-$\theta$ search for any $\epsilon$ for any number of iterations.

When $\epsilon_{0}(\theta)=\epsilon_{0}(\pi/3)=\epsilon$, in [7] we show
$\epsilon_{1}(\theta)>\epsilon_{1}(\pi/3)$. Assume that
$\epsilon_{m-1}(\theta)>\epsilon_{m-1}(\pi/3)$. From Eq. (10), it is easy to
see that $\epsilon_{m}(\theta)>\epsilon_{m}(\pi/3)$. Therefore,
$\epsilon_{m}(\theta)>\epsilon_{m}(\pi/3)$ for any $m$. Hence, when
$0<\theta<\pi/3$, the Phase-$\pi/3$ search converges faster than the
Phase-$\theta$ search for any $\epsilon$ for any number of iterations.

##  4 For any known $\epsilon$, $O(n)$ iterations can find the target state.

Assume that a database has $N=2^{n}$ states (items). Then a state (an item) is
found with the probability of $1/N$[2]. In other words, the failure
probability $\epsilon=1-1/N$. It is known that the Phase-$\pi/3$ search
converges the target state. In this section, we investigate how to use the
fixed-point search to find the target state in a database when $\epsilon$ is
known. As discussed in [4], the fixed-point search is a recursive algorithm,
therefore the number of queries grows exponentially with the number of
recursion levels. For example, the Phase-$\pi/3$ search at $i$-level recursion
involves $q_{i}=(3^{i}-1)/2$ queries [6]. This implies that $O(n)$ iterations
of the Phase-$\theta$ search involve the exponential queries.

###  4.1 When $\epsilon\leq 3/4$, only one iteration is needed to find the
target state.

When $0\leq\epsilon\leq\frac{3}{4}$,
$\left|1-\frac{1}{2(1-\epsilon)}\right|\leq 1$. Let
$\cos\theta=1-\frac{1}{2(1-\epsilon)}$. Then $D(\theta)=0$. Therefore, if
$\epsilon$ is fixed and $0\leq\epsilon\leq\frac{3}{4}$, then we choose
$\theta=\arccos[1-\frac{1}{2(1-\epsilon)}],$ which is in $(\pi/3$ ,$\pi]$, as
phase shifts. The Phase-$\arccos[1-\frac{1}{2(1-\epsilon)}]$ search will
obviously make the deviation vanish. It means that one iteration will reach
$t$ state if the $\theta$ is chosen as phase shifts. Ref. [7].

###  4.2 When $\epsilon>3/4$, $O(n)$ iterations can find the target state.

####  4.2.1 First use the Phase-$\pi/3$ search

For the Phase-$\pi/3$ search, $\epsilon_{n}=\epsilon^{3^{n}}$. There exists
the least natural number $n^{\ast}$ such that $\epsilon^{3^{n^{\ast}}}\leq
3/4$. By calculating,
$n^{\ast}=\lceil(\ln\ln\frac{4}{3}-\ln\ln\frac{1}{\epsilon})/\ln 3\rceil$.

Lemma 1. For the Phase-$\pi/3$ search, $n^{\ast}=$ $O(n)$.

Proof. In the case of database search, Let $N=2^{n}$. Then
$\epsilon=1-2^{-n}$, and
$\lim_{n\rightarrow+\infty}\frac{n^{\ast}}{n}=\frac{\ln 2}{\ln 3}$. Almost
$\frac{n\ln 2}{\ln 3}\approx\allowbreak\lceil 0.63n\rceil$. Let $N=10^{n}$.
Then $\epsilon=1-10^{-n}$, and
$\lim_{n\rightarrow+\infty}\frac{n^{\ast}}{n}=\frac{\ln 10}{\ln
3}=\frac{1}{\lg 3}$. Almost $\frac{n}{\lg 3}\approx 2n$. Thus,
$n^{\ast}=O(n)$. Hence, when $\epsilon>3/4$, after $n^{\ast}$ iterations of
the Phase-$\pi/3$ search the failure probability $\epsilon_{n^{\ast}}\leq$
$3/4$. Then, after one iteration of the
Phase-$\arccos[1-\frac{1}{2(1-\epsilon_{n^{\ast}})}]$ search by using the
result in section 4.1, it will reach $t$ state.

Example 5. Let $N=10^{4}$. Then $\epsilon=1-10^{-4}$, $n^{\ast}=8$,$\
\epsilon_{7}=\allowbreak 0.803\,32$, $\epsilon_{8}=\allowbreak 0.518\,4$. See
Fig.1. However, for this purpose, it only needs 4 iterations for the
Phase-$\pi$ search. See example 7.

Example 6. Let $N=2^{10}$. Then $\epsilon=1-2^{-10}$, $n^{\ast}=6$,
$\epsilon_{5}=\allowbreak 0.788\,56$, $\epsilon_{6}=\allowbreak 0.490\,35$.

####  4.2.2 First use the Phase-$\theta$ ($\neq\pi/3$) search

Let $\epsilon>3/4$. Then, by property (1.2), for the Phase-$\theta$ search,
there exists the least natural number $m^{\ast}(\theta)$ such that $\
\epsilon_{0}>\epsilon_{1}>...>\epsilon_{m^{\ast}(\theta)-1}(>3/4)>\epsilon_{m^{\ast}(\theta)}(\leq
3/4)$…. Thus, after $m^{\ast}(\theta)$ iterations of the Phase-$\theta$
search, the failure probability $\epsilon_{m^{\ast}(\theta)}\leq 3/4$. Then,
after one iteration for the
Phase-$\arccos[1-\frac{1}{2(1-\epsilon_{m(\theta)^{\ast}})}]$ search by using
the result in section 4.1, it will reach $t$ state.

Next let us calculate $m^{\ast}(\theta)$. Let $\delta=1-\epsilon$, where
$\delta$ is the success probability. When $\epsilon$ is close to $1$, $\delta$
is close to $0$. Then, for large $\epsilon$, by induction
$\epsilon_{l}=1-[1+4(1-\cos\theta)]^{l}\delta+O(\delta^{2})$. Thus,
$\epsilon_{l}\approx 1-[1+4(1-\cos\theta)]^{l}\delta$. By this approximate
formula of $\epsilon_{l}$, $m^{\ast}(\theta)$ $\approx
M^{\ast}(\theta)=\lceil\frac{-2\lg
2-\lg\delta}{\lg(1+4(1-\cos\theta))}\rceil$.

In the case of database search, let $N=2^{n}$. Then $\epsilon=1-2^{-n}$,
$\delta=2^{-n}$, and $m^{\ast}(\theta)$ $\approx$
$M^{\ast}(\theta)=\lceil\frac{(n-2)\lg 2}{\lg(1+4(1-\cos\theta))}\rceil$. For
the Phase-$\pi/3$ search, $M^{\ast}(\pi/3)=\lceil\frac{n\ln 2}{\ln
3}-\frac{2\ln 2}{\ln 3}\rceil$. Note that $\frac{2\ln 2}{\ln 3}=\allowbreak
1.\,\allowbreak 261\,9$. Therefore, when $n$ is large enough
$M^{\ast}(\pi/3)\approx$ $m^{\ast}(\pi/3)=n^{\ast}$. For the Phase-$\pi$
search, $m^{\ast}(\pi)\approx M^{\ast}(\pi)=\lceil(n-2)\lg 2/(2\lg
3)\rceil\approx(\lg 2)n$. See Table (I).

Let $N=10^{n}$. Then $\epsilon=1-10^{-n}$, $\delta=10^{-n}$, and
$m^{\ast}(\theta)\approx M^{\ast}(\theta)=\lceil\frac{n-2\lg
2}{\lg(1+4(1-\cos\theta))}\rceil$. For the Phase-$\pi/3$ search,
$M^{\ast}(\pi/3)=\lceil\frac{n}{\lg 3}-\frac{2\lg 2}{\lg 3}\rceil$. Note that
$\frac{2\lg 2}{\lg 3}=\frac{2\ln 2}{\ln 3}$. Therefore, when $n$ is large
enough $M^{\ast}(\pi/3)\approx$ $m^{\ast}(\pi/3)=n^{\ast}$. For the
Phase-$\pi$ search, $m^{\ast}(\pi)\approx M^{\ast}(\pi)=\lceil(n-2\lg 2)/(2\lg
3)\rceil\approx n$. See

Example 7. Let $N=10^{4}$. Then $\epsilon=1-10^{-4}$, $M^{\ast}(\pi)=4$,
$\epsilon_{4}=0.475\,32$. See Table (II).

Lemma 2. For the Phase-$\theta$ $(\neq\pi/3)$ search, $m^{\ast}(\theta)=$
$O(n)$.

Proof. When $\pi/3<\theta\leq\pi$, as discussed in case 1 of (1) in Sec. 3,
$m^{\ast}(\theta)<n^{\ast}$. By lemma 1, this lemma holds. When
$0<\theta<\pi/3$, from the approximate formula of $m^{\ast}(\theta)$,
$m^{\ast}(\theta)=$ $O(n)$.

Remark. $M^{\ast}(\theta)$ monotonically decreases as $\theta$ increases from
$0$ to $\pi$, especially $\frac{M^{\ast}(\pi)}{n^{\ast}}\approx 1/2$.
Therefore, we suggest first to use Phase-$\pi$ search for $m^{\ast}(\pi)$
times to get the failure probability $\epsilon_{m^{\ast}}\leq 3/4$.

##  5 Summary

In this paper, we investigate convergence performance of the Phase-$\theta$
search for any number of iterations. We discuss the convergence region and
rate of the Phase-$\theta$ search and study the convergence behavior of the
Phase-$\theta$ search for different initial $\epsilon_{0}$.

Acknowledgement

We want to thank the reviewer of [7] for suggesting us to study the
convergence behavior of the fixed-point search with general but equal phase
shifts for any number of iterations.

##  6 Appendix A

Proof. Since $f^{\prime}(x)<0$ when $r<x<d$, $f(g)\leq\epsilon_{1}<a$. Note
that $r<f(g)$. Thus, $r<f(g)\leq\epsilon_{1}<a$. Let
$f^{(k)}(x)=f(f^{(k-1)}(x))$. Since $f^{\prime}(x)<0$, $a<\epsilon_{2}\leq
f^{(2)}(g)<f(r)=g$ and $r<f(g)<f^{(3)}(g)\leq\epsilon_{3}<a$. By induction,
generally $a<\epsilon_{2k}\leq f^{(2k)}(g)<f^{(2k-2)}(g)<...f^{(2)}(g)<g$ and
$r<f(g)<...<f^{(2k-1)}(g)<f^{(2k+1)}(g)\leq\epsilon_{2k+1}<a$. That is,
$\epsilon_{i}$ oscillate around the fixed point $a$ by property 1 and between
$f^{(2k)}(g)$ and $f^{(2k+1)}(g)$. It is plain that the sequence
$\\{f^{(2k)}(g)\\}$ decreases monotonically as $k$ increases while the
sequence $\\{f^{(2k+1)}(g)\\}$ increases monotonically as $k$ does. Hence, the
sequences $\\{f^{(2k)}(g)\\}$ and $\\{f^{(2k+1)}(g)\\}$ have limits. Let
$\lim_{k\rightarrow\infty}f^{(2k)}(g)=\alpha$ and
$\lim_{k\rightarrow\infty}f^{(2k+1)}(g)=\beta$. Clearly, $\alpha$, $\beta<d$.
From Eq. (4),
$f^{(2k)}(g)=4(1-\cos\theta)^{2})f^{(2k-1)}(g)(f^{(2k-1)}(g)-d)^{2}$ and
$f^{(2k+1)}(g)=4(1-\cos\theta)^{2})f^{(2k)}(g)(f^{(2k)}(g)-d)^{2}$. By taking
the limits, we obtain $\alpha=4(1-\cos\theta)^{2})\beta(\beta-d)^{2}$ and
$\beta=4(1-\cos\theta)^{2})\alpha(\alpha-d)^{2}$. By substituting,
$\beta=[4(1-\cos\theta)^{2})]^{2}\beta(\beta-d)^{2}(\alpha-d)^{2}$. By
cancelling, $[4(1-\cos\theta)^{2})]^{2}(\beta-d)^{2}(\alpha-d)^{2}=1$. Then,
there are two cases. Case 1. $4(1-\cos\theta)^{2})(d-\beta)(d-\alpha)=1$. By
solving this equation, $\alpha=\beta=1$ or $\alpha=\beta=a$. Since $\alpha$,
$\beta<d<1$, then $\alpha=\beta=a$. Case 2.
$4(1-\cos\theta)^{2})(d-\beta)(d-\alpha)=-1$. There is no solution because
$\alpha$, $\beta<d$. Therefore,
$\lim_{k\rightarrow\infty}f^{(2k)}(g)=\lim_{k\rightarrow\infty}f^{(2k+1)}(g)=a$.
Then, $\lim_{k\rightarrow\infty}f^{(k)}(g)=a$, and also
$\lim_{m\rightarrow\infty}\epsilon_{m}=a$. We finish the proof.

## References

  * [1] L.K.Grover, Phys. Rev. Lett. 79 (1997) 325.
  * [2] L.K.Grover, Phys. Rev. Lett. 80 (1998) 4329.
  * [3] D. Li et al., Theor. Math. Phys. 144(3) (2005) 1279-1287.
  * [4] L.K.Grover, Phys. Rev. Lett. 95 (2005) 150501.
  * [5] G. Brassard, Science 275 (1997) 627.
  * [6] T.Tulsi, L. Grover, and A. Patel, quant-ph/0505007. Also, Quant. Inform. and Comput. 6(6) (2006) 483–494.
  * [7] D. Li et al., Eur. Phys. J. D 45 (2007) 335-340.
  * [8] D. Li et al., Phys. Lett. A 362 (2007) 260-264. Also see quant-ph/0604062.


# A model-free procedure to correct for volume fluctuations in E-by-E analyses
of particle multiplicities

Anar Rustamov  Joachim Stroth  Romain Holzmann  GSI Helmholtzzentrum für
Schwerionenforschung GmbH, 64291 Darmstadt, Germany  Institut für Kernphysik,
Goethe-Universität, 60438 Frankfurt am Main, Germany  Helmholtz
Forschungsakademie Hessen für FAIR, Campus Frankfurt,
60438 Frankfurt am Main, Germany

(August 27, 2024)

###### Abstract

We develop an innovative and unbiased procedure, based on event mixing, to
account for unavoidable contributions from volume (or system size)
fluctuations to experimentally measured moments of particle multiplicity
distributions produced in relativistic nuclear collisions. Within the wounded-
nucleon model they are characterized by fluctuations of the number of wounded
nucleons, the latter usually referred to as participants. For the first time
we extract participant fluctuations directly from the data used for the
fluctuation analysis, i.e., without involving model calculations. To achieve
this we constructed a dedicated event-mixing algorithm that eliminates all
possible correlations between produced particles while preserving the volume
fluctuations. The procedure provides direct access to the cumulants of
wounded-nucleon distributions, which can be used to account for non-critical
contributions to the experimentally measured cumulants of multiplicity
distributions.

###### keywords:

heavy-ion collisions, event-by-event observables, critical phenomena, QCD
phase diagram

##  1 Introduction

Fluctuations of conserved charges such as baryon number ($B$), electric charge
($Q$), strangeness ($S$) is the focus of current experimental and theoretical
investigations of relativistic nuclear collisions worldwide [1, 2, 3, 4]. They
are identified as promising probes of the Equation of State (EoS) of a system
under the study through its response to infinitesimal changes in external
parameters, thereby providing direct access to the phase structure of the
system. At center-of-mass energies available at LHC down to about 12 GeV per
nucleon pair recent theoretical calculations imply that chiral symmetry is
restored in a smooth crossover transition [5, 6, 7]. In the region of high
net-baryon density, where systematic lattice QCD calculations are not directly
applicable, effective model calculations suggest that the strongly interacting
matter undergoes a first order chiral phase transition [8, 9, 1, 10, 11]. The
conjectured chiral critical end point (CEP), at which the matter exhibits a
second order phase transition, should terminate the anticipated first order
chiral phase transition line. Extensive experimental efforts to measure net-
proton number fluctuations are underway [12, 13, 14, 15, 16, 17, 18], in order
to verify these theoretical predictions of a smooth crossover and to locate
the conjectured CEP in the phase diagram of strongly interacting matter.

Within the Grand Canonical Ensemble (GCE) of statistical mechanics net-charge
fluctuations can be quantified with generalized susceptibilities. In a thermal
system of fixed volume $V$ at a temperature $T$ the $n^{th}$ order
susceptibilities $\chi_{q}^{n}$ are defined as derivatives of the logarithm of
partition function $Z$ with respect to the chemical potentials $\mu_{q}$
responsible for the conservation of the corresponding charge $q\in\\{B,Q,S\\}$
on average, and evaluated at vanishing values of chemical potentials [3, 4]:

$\chi_{n}^{q}\equiv\frac{1}{VT^{3}}\left.\frac{\partial^{n}lnZ(T,V,\mu_{B},\mu_{Q},\mu_{S})}{\partial\hat{\mu}_{q}^{n}}\right\rvert_{\vec{\mu}=0}=\frac{1}{VT^{3}}\kappa_{n}(N_{q})\;.$
(1)

Here $\hat{\mu}_{q}=\mu_{q}/T$ and $\kappa_{n}(N_{q})$ is the $n^{th}$ order
cumulant of net-charge number distribution $N_{q}$, which can be measured in
experiments. Other parameters of the system, such as pressure, entropy density
etc., can be calculated by taking appropriate derivatives of $ln(Z)$.

In principle Eq. 1 already establishes a direct link between experimentally
measured cumulants $\kappa_{n}(N_{q})$ and theoretically evaluated
susceptibilities $\chi_{n}^{q}$, thus allowing to probe the EoS of the system,
which is encoded in the partition function $Z$. Clearly the experimental
measurements are to be done in a subspace of the full phase space, yet keeping
in mind that in the full phase space there are no fluctuations of conserved
charges, i.e., they are conserved in each event. And this is exactly what is
missing in Eq. 1. The solution to this problem is provided by exploiting the
Canonical Ensemble of statistical mechanics in the full phase space to account
for event-by-event net-charge conservation effects inside the selected
subspace [19, 20, 21, 22, 23].

However, for a proper comparison of experimental results with theory
predictions several additional contributions to the experimental measurements
need to be controlled with high precision [24]. Among those are fluctuations
due to the statistical nature of the detector response [25, 26, 27, 28, 29,
30], effects stemming from detecting only proxies (e.g. protons) of the true
baryon number fluctuations [31, 28], a mitigation of the critical signal due
to residual dynamics [32, 33, 34], and, last but not least, an insufficient
control of the volume associated with the particle-emitting source. In this
paper, we will address the rather complicated problem of volume fluctuations
for which there is so far no solution proven to work also at low collision
energies. Equation 1 is valid only if the system volume is fixed, a condition
which is hard, if not impossible, to achieve in experiments where the volume
is not directly observable. Commonly, in particular at high collision energies
where the colliding nuclei are sufficiently Lorentz contracted in the center-
of-mass, the volume is associated with the number of wounded nucleons [35],
usually referred to as participants [24]. In fact, experimental data are
analyzed in centrality percentiles corresponding to $n\%$ most central
collisions by introducing selection criteria e.g. in the energy deposited in a
forward calorimeter or the multiplicity of charged particles, either in the
full acceptance or in its sub-ranges etc. [36, 37]. For the latter, care must
be taken to ensure that the evaluated particles are not simultaneously used to
determine the critical fluctuations [16]. By fitting the so obtained
distributions with the Glauber Monte Carlo Model [38], where in addition it is
assumed that particles are produced from statistically independent sources,
the distributions of participants corresponding to a given centrality class
are obtained. Consequently, for each centrality class the number of
participants fluctuates from event to event, which prevents direct
exploitation of Eq. 1 [24].

To date, three conceptually distinct approaches have been developed to
circumvent this crucial experimental artifact in measurements of event-by-
event particle number fluctuations: (i) By evaluating contributions from
volume [39] or participant [24] fluctuations within the model of independent
sources for particle production. The latter approach, however, needs
probability distributions of participants as input. In [24] a detailed
procedure is worked out to account for participant fluctuations using the
Glauber Monte Carlo model. (ii) Introducing the so-called strongly intensive
quantities which, by a specific construction, eliminate participant
fluctuations [40]. This method also assumes statistically independent particle
production. (iii) Using an unfolding algorithm which also relies on model
calculations [41].

In this work we propose, as a radically new approach, to reconstruct the
participant fluctuations directly from the data used for the fluctuation
analysis. This is achieved by constructing a dedicated algorithm for event
mixing which removes all correlations between particles while keeping the
participant fluctuations unaffected. The paper is organized in the following
way. In section 2 we present definitions and notations. The main idea behind
the method is introduced in sections 3 and 4, in which the necessary analytic
expressions are given and the procedure for event mixing is introduced. The
results obtained are presented in section 5 and conclusions are given in
section 6.

##  2 Definitions and notations used in the paper

In the following we use the notations for cumulants introduced in [24]. The
$r^{th}$ central moment of a discrete random variable $X$, with probability
distribution $P(X)$, is defined as

$\mu_{r}\equiv\left<\left(X-\left<X\right>\right)^{r}\right>=\sum_{X}\left(X-\left<X\right>\right)^{r}P(X),\\\
$ (2)

where $\left<X\right>$ denotes the mean of the distribution

$\left<X\right>=\sum_{X}XP(X).\\\ $ (3)

In a similar way, we introduce moments about the origin, thereafter referred
to as raw moments

$\left<X^{r}\right>=\sum_{X}X^{r}P(X).\\\ $ (4)

The cumulants of $X$ are defined as the coefficients in the Maclaurin series
of the logarithm of the characteristic function of $X$. The first four
cumulants read

$\displaystyle\kappa_{1}=\left<X\right>,$
$\displaystyle\kappa_{2}=\mu_{2}=\left<X^{2}\right>-\left<X\right>^{2},$
$\displaystyle\kappa_{3}=\mu_{3}=\left<X^{3}\right>-3\left<X^{2}\right>\left<X\right>+2\left<X\right>^{3},$
(5)
$\displaystyle\kappa_{4}=\mu_{4}-3\mu_{2}^{2}=\left<X^{4}\right>-4\left<X^{3}\right>\left<X\right>-3\left<X^{2}\right>^{2}$
$\displaystyle+12\left<X^{2}\right>\left<X\right>^{2}-6\left<X\right>^{4}.$

For a random variable $X$, following the Poisson probability distribution
$P(X=k;\lambda)=e^{-\lambda}\lambda^{k}/k!$, any order cumulants of its
distribution are equal to its mean, that is $\kappa_{n}(X)=\lambda$.

##  3 Factorization approach

The cumulants $\kappa_{n}$ reconstructed from the multiplicity distributions
contain both a signal component and a background contribution due to the
fluctuations of the participants. The ultimate goal of the analysis is to
disentangle the signal part from non-critical contributions before they are
compared to the theoretical calculations using e.g. Eq. 1. Contributions from
fluctuations of participants are usually accounted for by exploiting
statistically independent sources for particle production [24], such as the
Wounded Nucleon Model (WNM) [35]. Within WNM, the analytical expressions for
the cumulants of (net-)particle distributions can be obtained. The first four
cumulants read [39, 24]:

$\displaystyle\kappa_{1}(N)=\left<N_{W}\right>\kappa_{1}(n)=\left<N_{W}\right>\left<n\right>,$
(6)
$\displaystyle\kappa_{2}(N)=\left<N_{W}\right>\kappa_{2}(n)+\left<n\right>^{2}\kappa_{2}(N_{W}),$
(7)
$\displaystyle\kappa_{3}(N)=\left<N_{W}\right>\kappa_{3}(n)+3\left<n\right>\kappa_{2}(n)\kappa_{2}(N_{W})+\left<n\right>^{3}\kappa_{3}(N_{W}),$
(8)
$\displaystyle\kappa_{4}(N)=\left<N_{W}\right>\kappa_{4}(n)+4\left<n\right>\kappa_{3}(n)\kappa_{2}(N_{W})$
(9)
$\displaystyle+3\kappa_{2}^{2}(n)\kappa_{2}(N_{W})+6\left<n\right>^{2}\kappa_{2}(n)\kappa_{3}(N_{W})+\left<n\right>^{4}\kappa_{4}(N_{W}),$

where $\kappa_{n}(n)$, $\kappa_{n}(N)$ denote the cumulants of particles from
a single wounded nucleon (or source) and those obtained after averaging over
the wounded-nucleon distribution, respectively. Within the same framework of
statistically independent particle sources we obtain the following expression
for the covariance between the multiplicities of two distinct particle species
$N_{1}$, $N_{2}$:

$cov(N_{1},N_{2})=\left<N_{W}\right>cov(n_{1},n_{2})+\left<n_{1}\right>\left<n_{2}\right>\kappa_{2}(N_{W}),$
(10)

with $cov(n_{1},n_{2})$ and $cov(N_{1},N_{2})$ denoting the covariances due to
emission from a single source and after averaging over the wounded-nucleon
distribution, respectively. Obviously, _a priory_ the values of $k_{n}(n)$ and
$cov(n_{1},n_{2})$ are unknown – they actually represent the physics we are
after. In order to isolate them, as it is evident from Eqs. 6\- 10, we need
additional information, namely distributions of wounded nucleons or rather the
cumulants of wounded-nucleon distributions $k_{n}(N_{W})$. That is where the
model dependence sets in. Indeed, nearly in all experiments the cumulants of
wounded nucleons are obtained by exploiting models, typically Glauber Monte
Carlo simulations or transport codes. The latter are known to introduce
biases, in particular at low collision energies [42]. On the other hand, at
beam energies below 3 $A$ GeV the build-up of transverse energy and momentum
in the early stage of the collision, when projectile- and target-like
spectator volumes are not yet separated from the participant zone, can lead to
participant formation outside the geometrical overlap volume. In that case, a
strict scaling of the number of participants, as the sources of particle
emission, with the size of the fireball volume cannot be assumed anymore.

##  4 Model independent extraction of participant fluctuations

Figure 1: The strategy for event mixing used to remove correlations between
particles while preserving participant fluctuations.

The method we develop acts differently. In a first step we mix the events in a
way to preserve the fluctuations of wounded nucleons but remove all
correlations between particles; the procedure is illustrated in Fig. 1. For a
given experimental event we count the number of all reconstructed charged
particles $n_{ch}$. Next, we select randomly one charged track from that
event, while the rest of the charged particles are each randomly taken from
the following $n_{ch}-1$ events, one per each event. In this way, the newly
created mixed event contains $n_{ch}$ particles, but the composition of
particle species in the mixed event is randomized. The procedure is repeated
event-by-event until a desired statistics of mixed events is reached. The
correlations between all particle species in mixed events are completely
eliminated, yet the original wounded-nucleon fluctuations are preserved. This
is because the total number of charged particles in real and mixed events is
kept to be the same. This also implies that in Eq. 10 the single-source
covariance term vanishes, yielding:

$cov(N_{1},N_{2})=\left<n_{1}\right>\left<n_{2}\right>\kappa_{2}(N_{W})=\left<N_{1}\right>\left<N_{2}\right>\frac{\kappa_{2}(N_{W})}{\left<N_{W}\right>^{2}},$
(11)

In a similar way, using Eq. 7 and realizing that for mixed events
$\kappa_{n}(n)=\left<n\right>$ (Poisson-distributed particles, because fully
uncorrelated), we obtain the following expression:

$\kappa_{2}(N)=\left<N_{W}\right>\left<n\right>+\left<n\right>^{2}\kappa_{2}(N_{W})=\left<N\right>+\left<N\right>^{2}\frac{\kappa_{2}(N_{W})}{\left<N_{W}\right>^{2}},$
(12)

Using Eg. 11, the second cumulants of participant fluctuations are obtained
from covariances between multiplicities of different particle species $N_{1}$
and $N_{2}$:

$\frac{\kappa_{2}(N_{W})}{\left<N_{W}\right>^{2}}=\frac{cov(N_{1},N_{2})}{\left<N_{1}\right>\left<N_{2}\right>}$
(13)

In a similar way, using Eq. 12, participant fluctuations can be extracted from
the reconstructed second order cumulants for particle type $N$:

$\frac{\kappa_{2}(N_{W})}{\left<N_{W}\right>^{2}}=\frac{\kappa_{2}(N)}{\left<N\right>^{2}}-\frac{1}{\left<N\right>}$
(14)

The covariances and cumulants entering the right hand sides of Eqs. 13 and 14
are to be evaluated using mixed events. Note that, by exploiting Eqs. 8 and 9,
higher-order cumulants of participant distributions can also be obtained.

##  5 Results and discussions

Figure 2: Fluctuation measure $\Sigma[A,B]$ reconstructed for 10 different
particle pairs, indicated with colored symbols, as a function of collision
centrality for original (left panel) and mixed (right panel) IQMD events.

The developed procedure is tested on events generated for Au-Au collisions at
$\sqrt{s_{\mathrm{NN}}}$ = 2.4 GeV with the IQMD transport model [43, 44].
Cumulants of multiplicity distributions are reconstructed for five different
particle species: proton (p), pion ($\pi$), deuteron (d), triton (t) and
nucleus of $He$ isotope (${}^{3}He$). Moreover, covariances are reconstructed
between multiplicites of all possible combinations of particle pairs: [p,
${}^{3}He$], [p, d], [p, t], [p, $\pi$], [${}^{3}He$, d], [${}^{3}He$, t],
[${}^{3}He$, $\pi$], [d, t], [d, $\pi$], [t, $\pi$]. The statistical
uncertainties are estimated using the subsampling approach [14, 15]. The
analysis is performed in four different centrality classes of 0-10%, 10-20%,
20-30%, and 30-40$\%$, selected by applying appropriate windows in the impact
parameter distribution provided by IQMD. Alternatively one could apply
different centrality selection criteria to meet experimental centrality
selection methods, however for the results presented in this section this is
immaterial. The method presented works for any kind of centrality selection
used in experiments.

A sample of mixed events is generated as prescribed in section 4, further
taking care that only events falling into a given centrality class are
used111We cross-checked, that mixing events from different centrality classes
does not introduce significant biases.. In addition, only charged particles
within kinematic region delimited as $|y|<0.4$ and $0.4<p_{T}<$ 1.6 GeV/c are
used, with $y$ and $p_{T}$ referring to the rapidity (using proton mass) and
transverse momentum of charged particles. In order to verify that in the mixed
events all correlations between particle multiplicities are lifted, we present
in Fig. 2 the centrality dependence of the strongly intensive fluctuation
measure $\Sigma[A,B]$, estimated for 10 different particle pairs [A, B]
mentioned above. Within WNM it depends neither on the mean number of
participants (centrality, system size or system volume) nor on its
fluctuations, and is defined as [40]:

$\Sigma[A,B]=\frac{\left<A\right>\omega(B)+\left<B\right>\omega(A)-2cov(A,B)}{\left<A+B\right>},$
(15)

where $\omega(A)$ and $\omega(B)$ are scaled variances of multiplicity
distributions of particle species $A$ and $B$, respectively, and $cov(A,B)$ is
the covariance between the corresponding multiplicity distributions of
particles $A$ and $B$, $\left<A\right>$ and $\left<B\right>$ in Eq. 15 stand
for event-averaged mean multiplicities. It is evident from Eq. 15, that in the
case of missing correlations between multiplicities of particles $A$ and $B$,
the fluctuation measure $\Sigma[A,B]$ becomes unity. Indeed in this case
$cov(A,B)$ vanishes and both scaled variances $\omega(A)$ and $\omega(B)$
become unity.

Figure 3: Normalized second order cumulants of participants as obtained from
covariances (left panel) and single particle cumulants (right panel). Figure
4: Normalized second order cumulants of participants obtained from single
particle cumulants averaged over five different particle species (red symbols)
and those extracted from the covariances averaged over all combinations
between all particle pairs (blue symbols).

The left panel of Fig. 2 shows the centrality dependence of $\Sigma[A,B]$ for
real (original) events while in the right panel similar results for the mixed
events are presented. Inspection of Fig. 2 clearly demonstrates that in mixed
events the values of $\Sigma[A,B]$ are all becoming unity, hence indicating
that any possible correlations between multiplicities of every particle pair
are eliminated. After this verification procedure we use mixed events to
extract the fluctuations of participants. Hence we turn back to quantities
which do depend on fluctuations of participants, i.e., cumulants of
multiplicity distributions and covariances between them. In the following,
only mixed events are used. We first reconstruct covariances between
multiplicity distributions of all particle pairs and using Eq. 13 extract the
normalized second order cumulant of participant distributions
$\kappa_{2}(N_{W})/\left<N_{W}\right>^{2}$. Obtained results as a function of
centrality are presented in the left panel of Fig. 3. Alternatively, we obtain
the same information (cf. right panel of Fig. 3) using Eq. 14 and
reconstructed single particle cumulants, $\kappa_{2}(A)$. Fig. 3 demonstrates,
that the $\kappa_{2}(N_{W})/\left<N_{W}\right>^{2}$ values extracted from
covariances and single particle cumulants are in reasonable agreement with a
mild dependence on particle species, thus demonstrating the robustness of the
method. This is clearly visible in Fig. 4, where the centrality dependence of
the $\kappa_{2}(N_{W})/\left<N_{W}\right>^{2}$ values obtained from single
particle cumulants averaged over five different particle species (red symbols)
and those extracted from the covariances averaged over all combinations
between all particle pairs (blue symbols) are presented.

##  6 Conclusion

In summary, we present a model-independent approach to disentangle non-
critical contributions to experimentally measured particle number
fluctuations. Specifically we demonstrate a possibility of eliminating
contributions from participant fluctuations based on measurements only, i.e.,
without involving the models. The method was tested on events generated at
$\sqrt{s_{\mathrm{NN}}}$ = 2.4 GeV with the IQMD transport model. The
participant fluctuations are extracted in two different ways: (i) Using second
order cumulants of five different particles species and (ii) From covariances
between multiplicity distributions of all possible pairs of particles. Both
approaches yield similar results, demonstrating that the method developed is
robust enough. The proposed method is essential for nuclear collisions at low
collisions energies performed at GSI/SIS and BNL/RHIC, but can be applied at
LHC energies as well. Furthermore, the method is general enough to allow
reconstructing any higher-order cumulant of the participant distribution.

##  7 Acknowledgements

AR acknowledges stimulating discussions with Peter Braun-Munzinger and Johanna
Stachel.

## References

  * [1] M. A. Stephanov, K. Rajagopal, E. V. Shuryak, Signatures of the tricritical point in QCD, Phys. Rev. Lett. 81 (1998) 4816–4819.  arXiv:hep-ph/9806219, doi:10.1103/PhysRevLett.81.4816.
  * [2] A. Bazavov, et al., Skewness, kurtosis, and the fifth and sixth order cumulants of net baryon-number distributions from lattice QCD confront high-statistics STAR data, Phys. Rev. D 101 (7) (2020) 074502.  arXiv:2001.08530, doi:10.1103/PhysRevD.101.074502.
  * [3] A. Rustamov, Deciphering the phases of QCD matter with fluctuations and correlations of conserved charges, (10 2022).  arXiv:2210.14810.
  * [4] P. Braun-Munzinger, A. Rustamov, J. Stachel, QCD under extreme conditions, (11 2022).  arXiv:2211.08819.
  * [5] A. Bazavov, et al., Chiral crossover in QCD at zero and non-zero chemical potentials, Phys. Lett. B 795 (2019) 15–21.  arXiv:1812.08235, doi:10.1016/j.physletb.2019.05.013.
  * [6] G. A. Almasi, B. Friman, K. Redlich, Baryon number fluctuations in chiral effective models and their phenomenological implications, Phys. Rev. D 96 (1) (2017) 014027.  arXiv:1703.05947, doi:10.1103/PhysRevD.96.014027.
  * [7] S. Borsanyi, Z. Fodor, J. N. Guenther, S. K. Katz, K. K. Szabo, A. Pasztor, I. Portillo, C. Ratti, Higher order fluctuations and correlations of conserved charges from lattice QCD, JHEP 10 (2018) 205.  arXiv:1805.04445, doi:10.1007/JHEP10(2018)205.
  * [8] M. Asakawa, K. Yazaki, Chiral Restoration at Finite Density and Temperature, Nucl. Phys. A 504 (1989) 668–684.  doi:10.1016/0375-9474(89)90002-X.
  * [9] K. Rajagopal, F. Wilczek, Static and dynamic critical phenomena at a second order QCD phase transition, Nucl. Phys. B 399 (1993) 395–425.  arXiv:hep-ph/9210253, doi:10.1016/0550-3213(93)90502-G.
  * [10] M. A. Stephanov, QCD critical point and complex chemical potential singularities, Phys. Rev. D 73 (2006) 094508.  arXiv:hep-lat/0603014, doi:10.1103/PhysRevD.73.094508.
  * [11] C. Sasaki, B. Friman, K. Redlich, Density fluctuations in the presence of spinodal instabilities, Phys. Rev. Lett. 99 (2007) 232301.  arXiv:hep-ph/0702254, doi:10.1103/PhysRevLett.99.232301.
  * [12] A. Rustamov, Net-baryon fluctuations measured with ALICE at the CERN LHC, Nucl. Phys. A 967 (2017) 453–456.  arXiv:1704.05329, doi:10.1016/j.nuclphysa.2017.05.111.
  * [13] S. Acharya, et al., Closing in on critical net-baryon fluctuations at LHC energies: cumulants up to third order in Pb$-$Pb collisions, (6 2022).  arXiv:2206.03343, doi:10.1016/j.physletb.2022.137545.
  * [14] S. Acharya, et al., Global baryon number conservation encoded in net-proton fluctuations measured in Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 2.76 TeV, Phys. Lett. B 807 (2020) 135564.  arXiv:1910.14396, doi:10.1016/j.physletb.2020.135564.
  * [15] J. Adamczewski-Musch, et al., Proton-number fluctuations in $\sqrt{s_{NN}}$ =2.4 GeV Au + Au collisions studied with the High-Acceptance DiElectron Spectrometer (HADES), Phys. Rev. C 102 (2) (2020) 024914.  arXiv:2002.08701, doi:10.1103/PhysRevC.102.024914.
  * [16] J. Adam, et al., Nonmonotonic Energy Dependence of Net-Proton Number Fluctuations, Phys. Rev. Lett. 126 (9) (2021) 092301.  arXiv:2001.02852, doi:10.1103/PhysRevLett.126.092301.
  * [17] H. Adhikary, Search for the critical point of strongly interacting matter, in: 15th Conference on Quark Confinement and the Hadron Spectrum, 2022, p. .  arXiv:2211.10504.
  * [18] A. Rustamov, Overview of fluctuation and correlation measurements, Nucl. Phys. A 1005 (2021) 121858.  arXiv:2005.13398, doi:10.1016/j.nuclphysa.2020.121858.
  * [19] M. Nahrgang, T. Schuster, M. Mitrovski, R. Stock, M. Bleicher, Net-baryon-, net-proton-, and net-charge kurtosis in heavy-ion collisions within a relativistic transport approach, Eur. Phys. J. C 72 (2012) 2143.  arXiv:0903.2911, doi:10.1140/epjc/s10052-012-2143-6.
  * [20] A. Bzdak, V. Koch, V. Skokov, Baryon number conservation and the cumulants of the net proton distribution, Phys. Rev. C 87 (1) (2013) 014901.  arXiv:1203.4529, doi:10.1103/PhysRevC.87.014901.
  * [21] P. Braun-Munzinger, B. Friman, K. Redlich, A. Rustamov, J. Stachel, Relativistic nuclear collisions: Establishing a non-critical baseline for fluctuation measurements, Nucl. Phys. A 1008 (2021) 122141.  arXiv:2007.02463, doi:10.1016/j.nuclphysa.2021.122141.
  * [22] P. Braun-Munzinger, A. Rustamov, J. Stachel, Experimental results on fluctuations of conserved charges confronted with predictions from canonical thermodynamics, Nucl. Phys. A 982 (2019) 307–310.  arXiv:1807.08927, doi:10.1016/j.nuclphysa.2018.09.074.
  * [23] P. Braun-Munzinger, A. Rustamov, J. Stachel, The role of the local conservation laws in fluctuations of conserved charges, (7 2019).  arXiv:1907.03032.
  * [24] P. Braun-Munzinger, A. Rustamov, J. Stachel, Bridging the gap between event-by-event fluctuation measurements and theory predictions in relativistic nuclear collisions, Nucl. Phys. A 960 (2017) 114–130.  arXiv:1612.00702, doi:10.1016/j.nuclphysa.2017.01.011.
  * [25] A. Bzdak, V. Koch, Local Efficiency Corrections to Higher Order Cumulants, Phys. Rev. C 91 (2) (2015) 027901.  arXiv:1312.4574, doi:10.1103/PhysRevC.91.027901.
  * [26] A. Bzdak, R. Holzmann, V. Koch, Multiplicity-dependent and nonbinomial efficiency corrections for particle number cumulants, Phys. Rev. C 94 (6) (2016) 064907.  arXiv:1603.09057, doi:10.1103/PhysRevC.94.064907.
  * [27] M. Arslandok, A. Rustamov, TIdentity module for the reconstruction of the moments of multiplicity distributions, Nucl. Instrum. Meth. A 946 (2019) 162622\.  arXiv:1807.06370, doi:10.1016/j.nima.2019.162622.
  * [28] X. Luo, N. Xu, Search for the QCD Critical Point with Fluctuations of Conserved Quantities in Relativistic Heavy-Ion Collisions at RHIC : An Overview, Nucl. Sci. Tech. 28 (8) (2017) 112.  arXiv:1701.02105, doi:10.1007/s41365-017-0257-0.
  * [29] A. Rustamov, M. I. Gorenstein, Identity Method for Moments of Multiplicity Distribution, Phys. Rev. C 86 (2012) 044906.  arXiv:1204.6632, doi:10.1103/PhysRevC.86.044906.
  * [30] M. Gazdzicki, M. I. Gorenstein, M. Mackowiak-Pawlowska, A. Rustamov, Particle-Set Identification method to study multiplicity fluctuations, Nucl. Phys. A 1001 (2020) 121915.  arXiv:1903.08103, doi:10.1016/j.nuclphysa.2020.121915.
  * [31] M. Kitazawa, M. Asakawa, Relation between baryon number fluctuations and experimentally observed proton number fluctuations in relativistic heavy ion collisions, Phys. Rev. C 86 (2012) 024904, [Erratum: Phys.Rev.C 86, 069902 (2012)].  arXiv:1205.3292, doi:10.1103/PhysRevC.86.024904.
  * [32] M. Asakawa, M. Kitazawa, Fluctuations of conserved charges in relativistic heavy ion collisions: An introduction, Prog. Part. Nucl. Phys. 90 (2016) 299–342.  arXiv:1512.05038, doi:10.1016/j.ppnp.2016.04.002.
  * [33] V. Vovchenko, V. Koch, C. Shen, Proton number cumulants and correlation functions in Au-Au collisions at sNN=7.7–200 GeV from hydrodynamics, Phys. Rev. C 105 (1) (2022) 014904.  arXiv:2107.00163, doi:10.1103/PhysRevC.105.014904.
  * [34] V. Vovchenko, V. Koch, Thermodynamic approach to proton number fluctuations in baryon-rich heavy-ion matter created at moderate collision energies, Phys. Lett. B 833 (2022) 137368.  arXiv:2204.00137, doi:10.1016/j.physletb.2022.137368.
  * [35] A. Bialas, M. Bleszynski, W. Czyz, Multiplicity Distributions in Nucleus-Nucleus Collisions at High-Energies, Nucl. Phys. B 111 (1976) 461–476.  doi:10.1016/0550-3213(76)90329-1.
  * [36] B. Abelev, et al., Centrality determination of Pb-Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV with ALICE, Phys. Rev. C 88 (4) (2013) 044909.  arXiv:1301.4361, doi:10.1103/PhysRevC.88.044909.
  * [37] J. Adamczewski-Musch, et al., Centrality determination of Au + Au collisions at 1.23A GeV with HADES, Eur. Phys. J. A 54 (5) (2018) 85.  arXiv:1712.07993, doi:10.1140/epja/i2018-12513-7.
  * [38] C. Loizides, J. Nagle, P. Steinberg, Improved version of the PHOBOS Glauber Monte Carlo, SoftwareX 1-2 (2015) 13–18.  arXiv:1408.2549, doi:10.1016/j.softx.2015.05.001.
  * [39] V. Skokov, B. Friman, K. Redlich, Volume Fluctuations and Higher Order Cumulants of the Net Baryon Number, Phys. Rev. C 88 (2013) 034911.  arXiv:1205.4756, doi:10.1103/PhysRevC.88.034911.
  * [40] M. I. Gorenstein, M. Gazdzicki, Strongly Intensive Quantities, Phys. Rev. C 84 (2011) 014904.  arXiv:1101.4865, doi:10.1103/PhysRevC.84.014904.
  * [41] S. Esumi, K. Nakagawa, T. Nonaka, Reconstructing particle number distributions with convoluting volume fluctuations, Nucl. Instrum. Meth. A 987 (2021) 164802\.  arXiv:2002.11253, doi:10.1016/j.nima.2020.164802.
  * [42] M. S. Abdallah, et al., Measurements of Proton High Order Cumulants in $\sqrt{s_{{}_{\mathrm{NN}}}}$ = 3 GeV Au+Au Collisions and Implications for the QCD Critical Point, Phys. Rev. Lett. 128 (20) (2022) 202303.  arXiv:2112.00240, doi:10.1103/PhysRevLett.128.202303.
  * [43] C. Hartnack, R. K. Puri, J. Aichelin, J. Konopka, S. A. Bass, H. Stoecker, W. Greiner, Modeling the many body dynamics of heavy ion collisions: Present status and future perspective, Eur. Phys. J. A 1 (1998) 151–169.  arXiv:nucl-th/9811015, doi:10.1007/s100500050045.
  * [44] P. B. Gossiaux, J. Aichelin, The Importance of initial final state correlations for the formation of fragments in heavy ion collisions, Phys. Rev. C 56 (1997) 2109–2131.  arXiv:nucl-th/9706033, doi:10.1103/PhysRevC.56.2109.