# Ion-acoustic cnoidal waves in a quantum plasma

S. Mahmood∗1,2 and F. Haas1
1Physics Institute, Federal University of Rio Grande do Sul,
915051-970, Porto Alegre, RS, Brazil
2Theoretical Physics Division (TPD), PINSTECH
P .O. Nilore Islamabad 44000, Pakistan

###### Abstract

Nonlinear ion-acoustic cnoidal wave structures are studied in an unmagnetized
quantum plasma. Using the reductive perturbation method, a Korteweg-de Vries
equation is derived for appropriate boundary conditions and nonlinear periodic
wave solutions are obtained. The corresponding analytical solution and
numerical plots of the ion-acoustic cnoidal waves and solitons in the phase
plane are presented using the Sagdeev pseudo-potential approach. The
variations in the nonlinear potential of the ion-acoustic cnoidal waves are
studied at different values of quantum parameter $H_{e}$ which is the ratio of
electron plasmon energy to electron Fermi energy defined for degenerate
electrons. It is found that both compressive and rarefactive ion-acoustic
cnoidal wave structures are formed depending on the value of the quantum
parameter. The dependence of the wavelength and frequency on nonlinear wave
amplitude is also presented.

###### pacs:

52.35.Fp, 52.35.Sb, 67.10.Db

##  I Introduction

The study of nonlinear wave propagation in quantum plasmas has gain importance
due its application in understanding the particle or energy transport
phenomenon on short scale lengths i.e., in micro and nano scale electronic
devices and in dense compact stars r2 ; r3 ; r4 ; r5 ; r6 ; r7 ; r8 .
Typically, the quantum effects in plasmas become important when the Fermi
temperature, which is related to the particles density, becomes equal or
greater than the system’s thermal temperature or the inter-particle distance
becomes smaller or of the same order of the particle’s de Broglie thermal
wavelength. In order to study the dynamics in quantum plasmas, the quantum
hydrodynamic (QHD) model is frequently used r2 ; r3 ; r7 . The QHD model
consists of a set of equations describing the transport of charge, momentum
and energy in a quantum charged particle system interacting through a self
consistent electromagnetic field. In QHD model, the quantum effects appear
through the quantum statistical (Fermi) pressure and the Bohm potential (due
to quantum diffraction or tunneling effects). The QHD is useful to study
collective effects on short scale lengths and has its limitation for systems
that are large compared to the Fermi Debye lengths of the species in the
system. Quantum ion-acoustic waves were investigated by Haas et al. r2 ; r3
using the QHD model. They derived a Korteweg-de Vries (KdV) equation in the
weakly nonlinear amplitude wave limit for studying the propagation of ion-
acoustic solitons in a quantum plasma. It is reported that the compressive or
rarefactive soliton solution depends on a quantum parameter ($H_{e}$) defined
for degenerate electrons, which is the ratio of the electron plasmon energy to
the Fermi energy. In the fully nonlinear regime, the existence of periodic
traveling wave patterns were reported for ion-acoustic waves in quantum
plasmas. The arbitrary amplitude ion-acoustic solitary waves using the QHD
model in electron-ion quantum plasmas with the Sagdeev potential approach has
also been investigated r9 . Nonlinear electrostatic wave structures such as
solitons, envelope and shocks have been studied in quantum electron-ion (EI)
plasmas ali  but no one has reported the propagation of cnoidal wave
structures in quantum plasmas. The purpose of the present work is to
investigate the formation of ion-acoustic cnoidal waves in a quantum EI plasma
using the well known QHD model.

In classical plasmas, a lot of research work has been done in studying
nonlinear ion-acoustic wave structures such as solitons, cnoidal waves and
envelopes. Solitons are single pulse structures which are formed due to the
balance between nonlinearity and dispersion effects in the system r10  and
they consist of isolated hump or dip like wave profile with no rapid
oscillations inside the packet. Envelope structure contains both fast and slow
oscillations, obtained when nonlinearity balances the wave group dispersion
effects. The envelope is a localized modulated wave packet whose dynamics is
governed by the nonlinear Schrödinger (NLS) equation r11 ; r12 . The periodic
(cnoidal) wave is the exact nonlinear periodic wave solution of the KdV
equation with appropriate boundary conditions. These solutions of the KdV
equation are also termed as cnoidal waves because they are written in terms of
Jacobian elliptic-function cn. In general, the nonlinear periodic waves are
expressed in terms of Jacobian elliptic-functions such as dn, sn or cn, and
the nonlinear dn waves are believed to be generated in the defocusing region
of the ionospheric plasmas r15 ; r16 ; r17 ; r18 . The ion-acoustic soliton
and double layer structures are observed in auroral and magnetospheric plasmas
and also nonlinear periodic wave signals appear frequently in these
observations r19 . The periodic signals are also observed at the edge of the
tokamak plasma, which can be described by cnoidal waves r20 . Kono et al. r21
studied the higher order contributions in the reductive perturbation theory
for the nonlinear ion-acoustic wave propagation under the periodic boundary
condition. The nonlinear periodic wave solution for small amplitude Langmuir
waves in electron-ion plasmas was studied by Schamel r22 . Jovanovic and
Shukla r23  presented a solution in the form of a cnoidal wave provided the
minimum value of the electrostatic potential remain finite in studying
coherent electric field structures. Prudskikh r231  studied the ion-acoustic
nonlinear periodic waves in dusty plasmas. The ion-acoustic cnoidal wave and
the associated nonlinear ion flux in dusty plasmas was studied by Jain et al.
r232 . They derived the coupled evolution equations for the first and second
order potentials for ion-acoustic waves in unmagnetized dusty plasmas using
reductive perturbation method with appropriate boundary conditions. Kaladze at
al. r24  investigated acoustic cnoidal waves and solitons in unmagnetized
pair-ion (PI) plasmas consisting of the same mass ion species with different
temperatures. They reported the formation of both compressive and rarefactive
cnoidal wave structures in PI plasmas which depends on the temperature ratio
of PI species. Recently, Kaladze and Mahmood r25  studied the effect of
positrons density and nonthermal parameter kappa on the propagation of the
ion-acoustic cnoidal waves in electron-positron-ion plasmas. Saha and
Chatterjee studied electron acoustic periodic and solitary wave solutions in
unmagnetized r251  and magnetized r252 ; r253  quantum plasmas. They derived a
KdV equation using the reductive perturbation method and investigated the
associated nonlinear structures using bifurcation theory.

The manuscript is organized in the following way. In the next section, the
model and set of dynamic equations for studying nonlinear ion-acoustic waves
in unmagnetized quantum plasmas is presented. Using the reductive perturbation
method, the KdV equation is also derived with appropriate boundary conditions.
In Section III, the stationary wave solution is obtained for ion-acoustic
cnoidal waves using the Sagdeev potential approach. In Section IV, the
numerical analysis and plots are presented for degenerate plasma cases at
different plasma densities chosen from literature and the conclusion is drawn
in the final Section V.

##  II Basic Model and Derivation of Korteweg-de Vries Equation

In order to study the electrostatic nonlinear periodic (cnoidal) waves in
unmagnetized electron-ion (EI) quantum plasmas, we will derive a KdV equation
using the reductive perturbation method. A KdV equation for quantum ion-
acoustic waves has already been derived by Haas et al. r3  with emphasis on
localized solutions obtained under decaying boundary conditions. For cnoidal
waves, periodic boundary conditions are more appropriate, hence we will derive
again the KdV equation for quantum ion-acoustic waves. The set of dynamic
equations for ion-acoustic wave using QHD model is described as follows.

The ion continuity and momentum equations for ion fluid are given by

$\frac{\partial n_{i}}{\partial t}+\frac{\partial}{\partial x}(n_{i}u_{i})=0,$
(1) $\frac{\partial u_{i}}{\partial t}+u_{i}\frac{\partial u_{i}}{\partial
x}=-\frac{e}{m_{i}}\frac{\partial\phi}{\partial x}.$ (2)

The dynamic equation for the inertialess electron quantum fluid is described
by

$0=e\frac{\partial\phi}{\partial x}-\frac{1}{n_{e}}\frac{\partial
p_{e}}{\partial x}+\frac{\hslash^{2}}{2m_{e}}\frac{\partial}{\partial
x}\left(\frac{1}{\sqrt{n_{e}}}\frac{\partial^{2}}{\partial
x^{2}}\sqrt{n_{e}}\right).$ (3)

The Poisson equation is written as

$\frac{\partial^{2}\phi}{\partial
x^{2}}=\frac{e}{\varepsilon_{0}}(n_{e}-n_{i}),$ (4)

where $\phi$ is the electrostatic potential. The ion fluid density and
velocity are represented by $n_{i}$ and $u_{i}$ respectively, while $n_{e}$ is
the electron fluid density. Also, $m_{e}$ and $m_{i}$ are the electron and ion
masses, $-e$ is the electronic charge, $\varepsilon_{0}$ and $\hslash$ are the
dielectric and scaled Planck’s constants. In equilibrium, we have
$n_{e0}=n_{i0}=n_{0}$ (say). Here $p_{e}$ is the electron pressure and
$p_{e}(n_{e})$ is obtained from the equation of state for the electron fluid.
The electrons are assumed to obey the equation of state pertaining to one-
dimensional zero-temperature Fermi gas r3 ; r4 , which is
$p_{e}=m_{e}v_{Fe}^{2}n_{e}^{3}/3n_{0}^{2}$ where $n_{0}$ is the equilibrium
plasma density. Here $v_{Fe}$ is the Fermi velocity of electron, connected to
the Fermi temperature by $m_{e}v_{Fe}^{2}/2=k_{B}T_{Fe}$ and $k_{B}$ is the
Boltzmann constant. The last term on the right hand side of the momentum
equations for electrons quantum fluid is the quantum force, which arises due
to the quantum Bohm potential and gives quantum diffraction or quantum
tunneling effects due to the wave-like nature of the charged particles. The
quantum effects due to ions are ignored in the model as they have large
inertia in comparison with the electrons.

In order to find the nonlinear ion-acoustic periodic waves in a quantum
plasma, the set of nonlinear dynamic equations are written in a normalized
form as follows,

$\frac{\partial\tilde{n}_{i}}{\partial\tilde{t}}+\frac{\partial}{\partial\tilde{x}}(\tilde{n}_{i}\tilde{u}_{i})=0,$
(5)
$\frac{\partial\tilde{u}_{i}}{\partial\tilde{t}}+\tilde{u}_{i}\frac{\partial}{\partial\tilde{x}}\tilde{u}_{i}=-\frac{\partial\Phi}{\partial\tilde{x}},$
(6)
$0=\frac{\partial\Phi}{\partial\tilde{x}}-\frac{1}{2}\frac{\partial\tilde{n}_{e}^{2}}{\partial\tilde{x}}+\frac{H_{e}^{2}}{2}\frac{\partial}{\partial\tilde{x}}\left(\frac{1}{\sqrt{\tilde{n}_{e}}}\frac{\partial^{2}}{\partial\tilde{x}^{2}}\sqrt{\tilde{n}_{e}}\right),$
(7)
$\frac{\partial^{2}\Phi}{\partial\tilde{x}^{2}}=\tilde{n}_{e}-\tilde{n}_{i}.$
(8)

The normalization of space, time, ion velocity and electrostatic potential is
defined by $\tilde{x}\rightarrow\omega_{{}_{pi}}x/c_{s}$,
$\tilde{t}\rightarrow\omega_{{}_{pi}}t$, $\tilde{u}_{i}\rightarrow
u_{i}/c_{s}$ and $\Phi=e\phi/2k_{B}T_{Fe}$ respectively, where the ion plasma
frequency and ion-acoustic speed are
$\omega_{pi}=\left(n_{0}e^{2}/\varepsilon_{0}m_{i}\right)^{1/2}$ and$\
c_{s}=\left(2k_{B}T_{Fe}/m_{i}\right)^{1/2}$ respectively, and the non-
dimensional quantum parameter for electrons is defined as
$H_{e}=\hbar\omega_{pe}/2k_{B}T_{Fe}$ i.e., the ratio of electron plasmon
energy to the Fermi energy, here
$\omega_{pe}=\left(n_{0}e^{2}/\varepsilon_{0}m_{e}\right)^{1/2}$ is the
electron plasma frequency. The normalization of electron and ion fluid density
is defined as $\tilde{n}_{j}=n_{j}/n_{0}$ ($j=e,i$). In the following, for
simplicity we will not use the tilde sign.

In order to find an nonlinear evolution equation, a stretching of independent
variables $x,$ $t$ is defined as follows r10 ,

$\xi=\varepsilon^{1/2}(x-V_{0}t)\,,\quad\tau=\varepsilon^{3/2}t,$

where $\varepsilon$ is a small parameter and $V_{0}$ is the phase velocity of
the wave to be determined later on. The perturbed quantities can be expanded
in the powers of $\varepsilon$,

$n_{i1}=1+\varepsilon n_{i1}+\varepsilon^{2}n_{i2}+...,$ $n_{e1}=1+\varepsilon
n_{e1}+\varepsilon^{2}n_{e2}+...,$ $u_{i}=\varepsilon
u_{i1}+\varepsilon^{2}u_{i2}+...,$
$\Phi=\varepsilon\Phi_{1}+\varepsilon^{2}\Phi_{2}+...$ (9)

Moreover, $\partial/\partial x=\varepsilon^{1/2}\,\partial/\partial\xi$ and
$\partial/\partial t=\varepsilon^{3/2}\partial/\partial\tau-
V_{0}\varepsilon^{1/2}\partial/\partial\xi$.

From ion continuity and momentum equations the lowest order
$(\sim\varepsilon^{3/2})$ terms gives

$-V_{0}\frac{\partial n_{i1}}{\partial\xi}+\frac{\partial
u_{i1}}{\partial\xi}=0,$ (10) $-V_{0}\frac{\partial
u_{i1}}{\partial\xi}+\frac{\partial\Phi_{1}}{\partial\xi}=0,$ (11)
$\frac{\partial\Phi_{1}}{\partial\xi}-\frac{\partial n_{e1}}{\partial\xi}=0.$
(12)

The lowest order $(\sim\varepsilon)$ term of Poisson equation gives

$n_{i1}=n_{e1}.$ (13)

Multiplying Eq. (10) by $V_{0}$ and adding with Eq. (11), we have

$\frac{\partial
n_{i1}}{\partial\xi}=\frac{1}{V_{0}^{2}}\frac{\partial\Phi_{1}}{\partial\xi}.$
(14)

Using Eqs. (12), (13) and (14), we have

$V_{0}=\pm 1.$ (15)

which is the normalized phase velocity of the ion-acoustic wave. From now on,
we set $V_{0}=1$ without loss of generality.

Now integrating Eqs. (11), (12), (14) and using (13) we have

$n_{i1}=n_{e1}=\Phi_{1}+c_{1}(\tau)$ (16)

and

$u_{i1}=\Phi_{1}+c_{2}(\tau),$ (17)

where $c_{1}(\tau)$ and $c_{2}(\tau)$ are at this point arbitrary functions of
$\tau$ only.

Now collecting the next higher order terms from ion dynamic equations, we have

$-\frac{\partial n_{i2}}{\partial\xi}+\frac{\partial
u_{i2}}{\partial\xi}+\frac{\partial
n_{i1}}{\partial\tau}+\frac{\partial}{\partial\xi}(n_{i1}u_{i1})=0,$ (18)
$-\frac{\partial
u_{i2}}{\partial\xi}+\frac{\partial\Phi_{2}}{\partial\xi}+\frac{\partial
u_{i1}}{\partial\tau}+u_{i1}\frac{\partial u_{i1}}{\partial\xi}=0,$ (19)
$\frac{\partial\Phi_{2}}{\partial\xi}-\frac{\partial
n_{e2}}{\partial\xi}-\frac{1}{2}\frac{\partial}{\partial\xi}n_{e1}^{2}+\frac{H_{e}^{2}}{4}\frac{\partial^{3}n_{e1}}{\partial\xi^{3}}=0.$
(20)

The next higher order ($\sim\varepsilon^{2}$) term of Poisson equation gives

$\frac{\partial^{2}\Phi_{1}}{\partial\xi^{2}}=n_{e2}-n_{i2}.$ (21)

Adding Eqs. (18) and (19), we obtain

$\frac{\partial n_{i2}}{\partial\xi}=\frac{\partial
n_{i1}}{\partial\tau}+\frac{\partial}{\partial\xi}(n_{i1}u_{i1})+\frac{\partial
u_{i1}}{\partial\tau}+u_{i1}\frac{\partial}{\partial\xi}u_{i1}+\frac{\partial\Phi_{2}}{\partial\xi},$
(22)

and from Eq. (20) we have

$\frac{\partial
n_{e2}}{\partial\xi}=\frac{\partial\Phi_{2}}{\partial\xi}-\frac{1}{2}\frac{\partial}{\partial\xi}n_{e1}^{2}+\frac{H_{e}^{2}}{4}\frac{\partial^{3}n_{e1}}{\partial\xi^{3}}.$
(23)

Using $n_{i1}$, $u_{i1}$ and $n_{e1}$ from Eqs. (16) and (17) in Eqs. (22) and
(23), we obtain

$\frac{\partial
n_{i2}}{\partial\xi}=2\frac{\partial\Phi_{1}}{\partial\tau}+3\Phi_{1}\frac{\partial\Phi_{1}}{\partial\xi}+\left(2c_{2}+c_{1}\right)\frac{\partial\Phi_{1}}{\partial\xi}+\frac{\partial\Phi_{2}}{\partial\xi}$
(24)

and

$\frac{\partial
n_{e2}}{\partial\xi}=-\Phi_{1}\frac{\partial\Phi_{1}}{\partial\xi}-c_{2}\frac{\partial\Phi_{1}}{\partial\xi}+\frac{H_{e}^{2}}{4}\frac{\partial^{3}\Phi_{1}}{\partial\xi^{3}}+\frac{\partial\Phi_{2}}{\partial\xi}.$
(25)

Applying periodic boundary conditions we get $\partial
c_{1}/\partial\tau=\partial c_{2}/\partial\tau=0$, so that the functions
$c_{1}$ and $c_{2}$ become independent of both $\xi$ and $\tau$ and are from
now on constants.

Differentiating Eq. (21) and using Eqs. (24) and (25), and after some
simplifications, we have the KdV equation for the nonlinear dynamics of ion-
acoustic waves in a quantum plasma as follows,

$\frac{\partial\Phi}{\partial\tau}+2\Phi\frac{\partial\Phi}{\partial\xi}+c_{0}\frac{\partial\Phi}{\partial\xi}+D\frac{\partial^{3}\Phi}{\partial\xi^{3}}=0,$
(26)

where

$\displaystyle c_{0}$ $\displaystyle=$ $\displaystyle c_{1}+c_{2},$
$\displaystyle D$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(1-\frac{H_{e}^{2}}{4}\right).$ (27)

Here $\Phi_{1}$ has been replaced by $\Phi$.

In the above KdV equation (26), the term containing an arbitrary constant
$c_{0}$ can be removed with a Galilean transformation. Hence we can set
$c_{0}=0$ without loss of generality. It can be noticed easily from dispersive
coefficient $D$ that the cnoidal wave solution exist only when $H_{e}\neq 2$,
so that the dispersive coefficient does not disappear to balance the
nonlinearity.

##  III Nonlinear Periodic Wave Solutions

In order to find the steady state cnoidal and solitary waves solutions of the
KdV Eq. (26) for quantum ion-acoustic waves, we follow the same procedure as
already done in Refs.r24 ; r25 . Assume a solution $\Phi(\eta)$, where
$\eta=\xi-u\tau$ and $u$ is the velocity of the nonlinear structure moving
with the frame. Therefore, Eq. (26) can be written as

$D\frac{d^{3}\Phi}{d\eta^{3}}+\frac{d}{d\eta}\left(\Phi^{2}-u\Phi\right)=0.$
(28)

As said, the arbitrary constant $c_{0}$ has been ignored, since it gives just
a shift in the velocity of the nonlinear structure.

After integration of Eq. (28), we get the equation of a conservative nonlinear
oscillator i.e.,

$\frac{d^{2}\Phi}{d\eta^{2}}=-\frac{dW}{d\Phi},$ (29)

where its potential energy $W=W(\Phi)$ is defined as

$W(\Phi)=\frac{1}{3D}\Phi^{3}-\frac{u}{2D}\Phi^{2}+\rho_{0}\Phi.$ (30)

Here $\rho_{0}$ is an integration constant. The potential function $W(\Phi)$
has two points of extremum for $\Phi$ i.e. $\Phi=\tilde{\Phi}_{1,2}$ defined
by $\partial W/\partial\Phi=0$, which are given by

$\tilde{\Phi}_{1,2}=\frac{u}{2}\pm\sqrt{\frac{u^{2}}{4}-D\rho_{0}}.$ (31)

Thus, there are two equilibrium states. One of them defines a saddle point
while the other one represents a center point i.e., the bottom of the
potential well r24 . Moreover $u^{2}/4-D\rho_{0}>0$ must hold for real values.
The zero’s of the potential energy (30) i.e., $\Phi=z_{1},z_{2},z_{3}$ are
given as follows,

$z_{1}=0\,,\quad
z_{2,3}=\frac{3}{2}\left(\frac{u}{2}\pm\sqrt{\frac{u^{2}}{4}-\frac{4}{3}D\rho_{0}}\right).$
(32)

To get the shape of the real potential well $u^{2}/4-4D\rho_{0}/3>0$ must
holds. Note that the potential well having the center for $\Phi>0$ (positive)
defines the compressive cnoidal waves and solitons, while in the case of
$\Phi<0$ (negative) the potential well defines the rarefactive cnoidal waves
and solitons r24 . The shape of the potential well strongly depends on the
sign of dispersive coefficient $D$ (see Eq. (30)). As described in Eq. (27),
the value of the dispersive coefficient $D$ in the KdV equation is positive
for electron quantum parameter $H_{e}<2$ case, while it becomes negative for
$H_{e}>2$. Therefore, the formation of compressive or rarefactive ion-acoustic
nonlinear structure depends on the value of the quantum parameter $H_{e}$. The
amplitude of the nonlinear structures is defined by the width of the potential
well, which is the length between the last zero of the potential well and the
saddle point (see Eqs. (31) and (32)).

The energy first integral associated to (29) is

$\frac{1}{2}\left(\frac{d\Phi}{d\eta}\right)^{2}+W(\Phi)=\frac{E_{0}^{2}}{2},$
(33)

where $E_{0}^{2}$ is the integration constant (assumed positive definite in
order to access cnoidal wave solutions).

Using Eq. (30) in (33), we have

$\left(\frac{d\Phi}{d\eta}\right)^{2}=E_{0}^{2}-\frac{2}{3D}\Phi^{3}+\frac{u}{D}\Phi^{2}-2\rho_{0}\Phi.$
(34)

Let us consider the initial conditions $\Phi(0)=\Phi_{0}$ and
$d\Phi(0)/d\eta=0$. Then we can define

$E_{0}^{2}=\frac{2}{3D}\Phi_{0}^{3}-\frac{u}{D}\Phi_{0}^{2}+2\rho_{0}\Phi_{0}.$
(35)

Using Eq. (35) in Eq. (34) and after factorization, we have

$\left(\frac{d\Phi}{d\eta}\right)^{2}=\frac{2}{3D}(\Phi_{0}-\Phi)(\Phi-\Phi_{1})(\Phi-\Phi_{2}),$
(36)

where

$\Phi_{1,2}=\frac{3}{2}\left[\frac{u}{2}-\frac{\Phi_{0}}{3}\pm\sqrt{\frac{1}{3}(b_{1}-\Phi_{0})(\Phi_{0}-b_{2})}\right],$
(37)

and

$b_{1,2}=\frac{u}{2}\pm 2\sqrt{\frac{u^{2}}{4}-D\rho_{0}}.$ (38)

In addition, the following inequalities should be kept: $b_{2}\leq\Phi_{0}\leq
b_{1}$ or $b_{1}\leq\Phi_{0}\leq b_{2}$. From Eqs. (34)-(36), we have the
following relation,

$u=\frac{2}{3}(\Phi_{0}+\Phi_{1}+\Phi_{2}).$ (39)

The periodic (cnoidal) wave solution of Eq. (36) is given r26  by

$\Phi(\eta)=\Phi_{1}+(\Phi_{0}-\Phi_{1})\mathrm{cn}^{2}(R\eta,s),$ (40)

where cn is the Jacobian elliptic function, $s$ is the modulus defined as

$s^{2}=\frac{(\Phi_{0}-\Phi_{1})}{(\Phi_{0}-\Phi_{2})}<1,$ (41)

and the quantity $R=\sqrt{\frac{1}{6D}(\Phi_{0}-\Phi_{2})}$.

The amplitude $A$ of the cnoidal wave is defined from Eq. (40) as follows,

$A=\Phi_{0}-\Phi_{1}.$ (42)

As it is seen from the solution (40) at $\eta=0$, we have the initial
condition $\Phi(0)=\Phi_{0}$. In addition, the real numbers $\Phi_{i}$
$(i=0,1,2)$ are ordered as $\Phi_{0}>\Phi_{1}\geq\Phi_{2}$ and
$\Phi_{1}\leq\Phi\leq\Phi_{0}$ for $D>0$.

The modulus $0<s<1$ is a measure of the nonlinearity of the wave. The case
$s<<1$ corresponds to the weakly nonlinear oscillations near the bottom of the
potential well and the elliptic functions are close to trigonometric ones. At
$s\rightarrow 0$, the expression (40) passes to solution of linear equations
r26 .

The wavelength $\lambda$ of the cnoidal waves is defined as

$\lambda=4\sqrt{\frac{3D}{2(\Phi_{0}-\Phi_{2})}}K(s),$ (43)

where $K(s)$ is the complete elliptic integral of the first kind and the
corresponding frequency is $f=v_{1}/\lambda$ (where $v_{1}$ is the velocity of
the cnoidal waves). The velocity $v_{1}$ of the cnoidal waves in the
laboratory frame is equal to $v_{1}=V_{0}+u$, where the expression for the
frame velocity $u$ is given by,

$u=\frac{2}{3}\frac{(\Phi_{0}-\Phi_{1})}{s^{2}}(2-s^{2})+2\Phi_{2},$ (44)

which has been obtained using the expression of the modulus $s$ described by
Eq. (41).

The mean value of $\Phi$ can be expressed as

$\bar{\Phi}=\frac{1}{\lambda}\int_{0}^{\lambda}\Phi(\eta)d\eta=\Phi_{2}+(\Phi_{0}-\Phi_{2})\frac{E(s)}{K(s)},$
(45)

where $E(s)$ is the complete integral of the second kind.

The limiting case of the soliton i.e. $s=1$, can be obtained at
$\Phi_{0}\approx b_{1}$ or $\Phi_{0}\approx b_{2}$ (see Eq. (38)), so that

$\Phi_{1}\approx\Phi_{2}=\frac{u}{2}\mp\sqrt{\frac{u^{2}}{4}-D\rho_{0}}.$ (46)

Further, we take into account,

$K(s)\approx\frac{1}{2}\mathrm{ln}\left(\frac{16}{1-s^{2}}\right)\rightarrow\infty\,,\quad\mathrm{cn\,z}\rightarrow\frac{1}{\mathrm{ch\,z}},$
(47)

which imply that the wavelength of the cnoidal waves defined in (43) tends to
infinity and the solution (40) passes to a soliton-like shape r26  i.e.,

$\Phi(\eta)=\Phi_{1}+\frac{A}{\mathrm{ch}^{2}\left(\sqrt{\frac{1}{6D}(\Phi_{0}-\Phi_{1})}\eta\right)},$
(48)

where $\triangle=\sqrt{6\left|D/A\right|}$ is the width of the soliton and $A$
is its amplitude, defined by Eq. (42). As it follows from Eq. (44), the
propagation velocity of the solitons becomes

$u=2\Phi_{1}+\frac{2}{3}A.$ (49)

From Eq. (48), we see that $\Phi_{1}$ defines the potential at
$\eta\rightarrow\pm\infty$.

Thus, at large values of $s$ i.e., $s\rightarrow 1$, $A=const$ the periodic
wave asymptotically approaches to the sequence of solitons having the
amplitude $A$ (relatively to the level $\Phi=\bar{\Phi}=\Phi_{1}$). By the
order of magnitude, the distance between them is equal to

$\lambda=\bigtriangleup\left|\mathrm{ln}(1-s^{2})\right|,$ (50)

where $\bigtriangleup$ is the width of the soliton already defined above.

##  IV Numerical Analysis

The numerical plots of the nonlinear wave potential and the phase plane plots
of the cnoidal wave structures and solitons are shown in Figs.1-3 at different
densities for a degenerate electron plasma cases, such as astrophysical
plasmas, laser plasmas and ultra-cold plasmas. For a completely degenerate
electron plasma, the electron Fermi energy and density are related as
$k_{B}T_{Fe}=\hslash^{2}(3\pi^{2}n_{0})^{2/3}/2m_{e}$ and the Fermi
temperature of degenerate electrons ($T_{Fe}$) should be much greater than the
thermal temperature $T$ of the system i.e., $T\ll T_{Fe}$. The quantum
parameter for electrons is related to density as $H_{e}\sim n_{0}^{-1/6}$,
which shows that in a completely degenerate electron plasma case, the value of
the quantum parameter decreases with the increase in the plasma density. So
the quantum diffraction effects tend to be less relevant in dense plasmas. In
case of astrophysical plasma conditions i.e., $n_{0}=10^{36}m^{-3}$, the
quantum parameter for degenerate electrons comes out to be $H_{e}=0.05$ and
the condition for thermal temperature becomes $T<<10^{9}K$, while for laser
plasmas we have $n_{0}=10^{32}m^{-3}$ then $H_{e}=0.24$ and $T<<10^{7}K$.
Further, for ultra-cold plasmas, we have $n_{0}=2.7\times 10^{26}m^{-3}$ for
which $H_{e}=2.002$ and $T<<1800K$ r28 ; r29 .

Figure 1: (a) The nonlinear potential $W(\Phi)$ is plotted for
$\rho_{0}=0.01$, $u=0.2$ and $H_{e}=0.05$. (b) The phase plane plots of the
compressive ion-acoustic cnoidal wave (solid bounded curve) and soliton
(dotted curve) are shown for the same numerical values as in Fig. 1a.
Dimensionless variables are used. Figure 2: (a) The nonlinear potential
$W(\Phi)$ is plotted for $\rho_{0}=0.01$, $u=0.2$ and $H_{e}=0.24$ (b) The
phase plane plots of the compressive ion-acoustic cnoidal wave (solid bounded
curve) and soliton (dotted curve) are shown for the same numerical values as
in Fig. 2a. Dimensionless variables are used. Figure 3: (a) The nonlinear
potential $W(\Phi)$ is plotted for $\rho_{0}=0.01$, $u=-0.05$ and
$H_{e}=2.002$ (b) The phase plane plots of the rarefactive ion-acoustic
cnoidal wave (solid bounded curve) and soliton (dotted curve) are shown for
the same numerical values as in Fig. 3a. Dimensionless variables are used.

The formation of ion-acoustic compressive (rarefactive) nonlinear structures
depends on the value of the quantum parameter for electrons i.e., $H_{e}<2$
($H_{e}>2$). It is also noticed that the velocity of the nonlinear structure
is positive i.e., $u>0$ for compressive cnoidal waves and solitons with
quantum parameter for degenerate electrons lies in the range $0\leq H_{e}<2$.
For the $H_{e}>2$ case, rarefactive ion-acoustic cnoidal waves and solitons
structures are formed and its solution exist only when the velocity of the
nonlinear structure is negative $u<0$ i.e., it moves in the backward
direction.

It can be seen from Figs. 1a and 2a that the Sagdeev potential $W(\Phi)$ are
formed for $\Phi>0$ for the electron quantum parameter values $H_{e}=0.05$ and
$0.24$. The corresponding compressive ion-acoustic cnoidal wave (solid curve)
and solitons (dotted curve) structures are shown in the phase plane plots of
Fig. 1b and 2b respectively. The cnoidal wave structure (solid bounded curve)
is formed inside the separatrix (dotted curve) which represent a soliton
structure as shown in the figures. The Sagdeev potential plot $W(\Phi)$ is
formed with $\Phi<0$ for the electron quantum parameter value $H_{e}=2.002$
($>2$) as shown in the Fig. 3a. The corresponding rarefactive ion-acoustic
cnoidal wave (solid bounded curve) and soliton (dotted curve) structures in
the phase plane plot are shown in the Fig. 3b. Rarefactive nonlinear ion-
acoustic cnoidal wave (solid bounded curve) structures are also formed inside
the separatrix (dotted curve), which represent the soliton.

The plots of the compressive ion-acoustic cnoidal wave from Eq. (40) and
solitons from Eq. (48) for $H_{e}<2$ case are shown in the Figs. 4a and 4b
respectively. It can be seen from the Fig. 4a that there is a little decrease
in the wavelength (frequency) of the compressive ion-acoustic cnoidal wave
case with the increase in the quantum parameter of the degenerate electrons.
The little increase in the width of the compressive ion-acoustic solitons with
the increase in the value of quantum parameter is shown in the Fig. 4b.
Similarly, the plots for the rarefactive ion-acoustic cnoidal wave and soliton
structures for $H_{e}>2$case are shown in the Figs. 5a and 5b respectively. It
can be seen from the Fig. 5a that wavelength (frequency) of the rarefactive
ion-acoustic cnoidal wave increases and its amplitude decreases significantly
with the little increase in the value of quantum parameter for $H_{e}>2$ case.
Also, the decrease in the wave amplitude as well as increase in the width of
the rarefactive ion-acoustic soliton with the increase in the value of quantum
parameter is shown in Fig. 5b.

Figure 4: (a) The plots of compressive ion-acoustic cnoidal waves (periodic
wave oscillations) from Eq. (40) and (b) solitons (single pulse) from Eq. (48)
are shown for the $H_{e}<2$ case i.e., with $H_{e}=0.05$ (solid) and
$H_{e}=0.24$ (dotted) respectively. Dimensionless variables are used. Figure
5: (a) The plots of compressive ion-acoustic cnoidal waves (periodic wave
oscillations) from Eq. (40) and (b) solitons (single pulse) from Eq. (48) are
shown for the $H_{e}>2$ case i.e., with $H_{e}=2.002$ (solid) and $H_{e}=2.01$
(dotted) respectively. Dimensionless variables are used.

The variations of the wavelength and frequency with wave amplitude of the
compressive ion-acoustic cnoidal waves at different quantum parameters are
shown in Fig. 6a and 6b, respectively. It can be seen from the figures that
the wavelength increases, while the frequency decreases with the increase in
the amplitude. Also, the decrease in the wavelength and increase in the wave
frequency is found with the increase in the electron quantum parameter for
compressive ($H_{e}<2$) ion-acoustic cnoidal waves case. However, this
decrease in the wavelength and increase in the wave frequency with the wave
amplitude for the ion-acoustic compressive cnoidal waves case seems to be very
small at the chosen degenerate plasma densities as indicated in Figs. 6a and
6b. The dependence of the wavelength and frequency on the wave amplitude in
case of rarefactive ion-acoustic cnoidal waves case i.e., at quantum
parameters with values $H_{e}>2$ are shown in Figs. 7a and 7b, respectively.
The wavelength increases while the frequency decreases with the increase in
the amplitude of the rarefactive ion-acoustic cnoidal waves in quantum
plasmas. On the other hand, the wavelength is found to be increased, while
frequency decreases for rarefactive ion-acoustic cnoidal waves case with the
increase in the value of the quantum parameter for degenerate electrons as
shown in the Figs. 7a and 7b.

Figure 6: (a) The dependence of wavelength $\lambda$ on compressive ion-
acoustic cnoidal wave amplitude $A$ is shown for degenerate electron quantum
parameter $H_{e}=0.05$ (solid curve) and $H_{e}=0.24$ (dotted curve) with
$v_{1}=0.1$. (b) The dependence of frequency $f$ on compressive ion-acoustic
cnoidal wave amplitude $A$ with same parameters as in Fig. 4a. Dimensionless
variables are used. Figure 7: (a) The dependence of wavelength $\lambda$ on
rarefactive ion-acoustic cnoidal wave amplitude $A$ is shown for degenerate
electron quantum parameter $H_{e}=2.002$ (solid curve) and $H_{e}=2.01$
(dotted curve). (b) The dependence of frequency $f$ on rarefactive ion-
acoustic cnoidal wave amplitude $A$ with the same parameters as in Fig. 6a.
Dimensionless variables are used.

##  V Conclusion

To conclude, we have studied for the first time the ion-acoustic cnoidal waves
and solitons in an unmagnetized quantum plasma. The KdV equation for ion-
acoustic waves in a quantum plasma was obtained using the reductive
perturbation method with periodic wave boundary conditions, appropriate to
study cnoidal waves. It is found that both compressive and rarefactive
nonlinear ion-acoustic cnoidal wave structures are formed in such a degenerate
plasma, which depends on the quantum parameter i.e., $H_{e}\gtrless 2$. The
dependence of wave frequency and wavelengths on the nonlinear ion-acoustic
wave amplitude is also investigated at different values of quantum parameters
with the degenerate plasma densities exist in astrophysical and laboratory
plasmas. It is found that the dependencies of wavelength and frequency on wave
amplitude at different quantum parameters for electrons behave differently for
compressive and rarefactive ion-acoustic cnoidal wave cases. The results are
useful to understand how nonlinear wave propagates in quantum plasmas.

## Acknowledgements

SM acknowledges CNPq (National Council for Scientific and Technological
Development) and TWAS (The World Academy of Sciences) for a CNPq-TWAS
postdoctoral fellowship. FH acknowledges CNPq for financial support.

## References

  * (1) F. Haas, “Quantum Plasmas: an Hydrodynamic Approach” (Springer New York, USA, 2011).
  * (2) F. Haas, L. G. Garcia, J. Goedert and G. Manfredi, Phys. Plasmas 10, 3858 (2003).
  * (3) G. Manfredi and F. Haas, Phys. Rev. B 64, 075316 (2001).
  * (4) P.A. Markowich, C. A. Ringhofer and C. Schmeiser, “Semiconductor Equations” (Springer, Vienna, 1990).
  * (5) M. Bonitz, N. Horing and P. Ludwig, “Introduction to Complex Plasmas” (Springer-Verlag, Berlin-Heidelberg, 2010) chap-10.
  * (6) R. E. Wyatt, “Quantum Dynamics with Trajectories: Introduction to Quantum Hydrodynamics” (Springer, New York, 2005).
  * (7) P.K. Shukla and B. Eliasson, Rev. Mod. Phys. 83, 885 (2011).
  * (8) S. Mahmood and A. Mushtaq, Phys. Lett. A 372, 3467 (2008).
  * (9) S. Ali, W. M. Moslem, P.K. Shukla and I. Kourakis, Phys. Lett. A 366, 606 (2007).
  * (10) R. C. Davidson, “Methods in Nonlinear Plasma Theory” (Academic Press New York,1972) p-15.
  * (11) A. Hasegawa, “Plasma Instabilities and Nonlinear Effects” (Springer-Verlag, Berlin-Heidelberg, 1975) p-190.
  * (12) S. V. Vladimirov, V. N. Tsytovich, S. I. Popel and F. K. Khakimov, “Modulational Interactions in Plasmas” (Kluwer Academic Publishers, Netherlands, 1995) p-87.
  * (13) H. Schamel, Plasma Phys. 14, 905 (1972); Phys. Plasmas 19, 020501 (2012).
  * (14) L. C. Lee and J. R. Khan, Phys. Fluids 24, 430 (1981).
  * (15) I. Kourakis and P.K. Shukla, Eur. Phys. J. D 29, 247 (2004).
  * (16) A. P. Misra and A. Roy Choudhury, Phys. Plasmas 13, 062307 (2006).
  * (17) A. V. Gurevich and L. Stenflo, Phys. Scr. 38, 855 (1998).
  * (18) U. Kauschke and H. Schlüter, Plasma Phys. Cont. Fusion 32, 1149 (1990); 33, 1309 (1991).
  * (19) K. Konno, T. Mitsuhashi and Y. H. Ichikawa, J. Phys. Soc. Jpn. 46, 1907 (1979).
  * (20) H. Schamel, J. Plasma Phys. 13, 139 (1975).
  * (21) D. Jovanovic and P.K. Shukla, Phys. Lett. 84, 4373 (2000).
  * (22) V. V. Prudskikh, Plasma Phys. Reports 35, 651(2009).
  * (23) S. L. Jain, R. S. Tiwari and M. K. Mishra, Phys. Plasmas 19, 103702 (2012).
  * (24) T. Kaladze, S. Mahmood and H. Ur-Rehman, Physica Scripta 86, 035506 (2012).
  * (25) T. Kaladze and S. Mahmood, Phys. Plasmas 21, 032306 (2014).
  * (26) A. Saha and P. Chatterjee, Astrophys. Space Sci. 349, 239 (2014).
  * (27) U. K. Samanta, A. Saha and P. Chatterjee, Phys. Plasmas 20, 052111 (2013).
  * (28) U. K. Samanta, A. Saha and P. Chatterjee, Astrophys. Space Sci. 347, 293 (2013).
  * (29) V. I. Karpman, “Nonlinear Waves in Dispersive Media” (Pergamon, Oxford, 1975).
  * (30) A. F. Alexandrov, L.S. Bogdankevich and A. A. Rukhadze, “Principles of Plasma Electrodynamics” (Springer-Verlag, Berlin-Heidelberg, 1984) Chapter 1.
  * (31) T. C. Killian, Science 316, 705 (2007).


# Pileup Correction on Higher-order Cumulants with Unfolding Approach

Yu Zhang<EMAIL_ADDRESS>Key Laboratory of Quark & Lepton Physics
(MOE) and Institute of Particle Physics, Central China Normal University,
Wuhan 430079, China  Lawrence Berkeley National Laboratory, CA 94720, USA
Yige Huang<EMAIL_ADDRESS>Key Laboratory of Quark & Lepton Physics
(MOE) and Institute of Particle Physics, Central China Normal University,
Wuhan 430079, China  Toshihiro Nonaka<EMAIL_ADDRESS>Tomonaga Center for the History of the Universe, University of Tsukuba,
Tsukuba, Ibaraki 305, Japan  Xiaofeng Luo<EMAIL_ADDRESS>Key Laboratory of
Quark & Lepton Physics (MOE) and Institute of Particle Physics, Central China
Normal University, Wuhan 430079, China

###### Abstract

Higher-order cumulants of conserved charge distributions are sensitive
observables to probe the critical fluctuations near QCD critical point in
heavy-ion collisions. Due to high interaction rate, pileup event can be one of
the major sources of background in the measurements of higher-order cumulants.
In this paper, we studied the effects of pileup events on higher-order
cumulants of proton multiplicity distributions using UrQMD model. It is found
that the proposed pileup correction fails if the correction parameters are
determined by the Glauber fitting of charged particle multiplicities, which is
usually done in the real heavy-ion experiment. To address this, we propose a
model independent unfolding approach to determine the parameters in the pileup
correction. This approach can be applied in the pileup correction for the
future measurement of higher-order cumulants in heavy-ion collision
experiment.

##  I Introduction

One of the main goals of heavy-ion collision experiments is to understand the
phase structure of QCD matter with respect to baryon chemical potential
($\mu_{\rm B}$) and temperature ($T$) Gupta  _et al._ (2011); Luo and Xu
(2017); Bzdak  _et al._ (2020). The most important question is the position of
a QCD critical end point, which is a connecting point between crossover and
first-order phase transition, in ($\mu_{\rm B}$, T) plane Stephanov (2004);
Fukushima and Hatsuda (2011); Shi  _et al._ (2014); Ding  _et al._ (2015); Gao
and Liu (2016); Bazavov  _et al._ (2017); Fischer (2019); Fu  _et al._ (2020);
Isserstedt  _et al._ (2019); Gao and Pawlowski (2020); Bernhardt  _et al._
(2021); Borsanyi  _et al._ (2021). In order to explore the QCD phase diagram,
higher-order fluctuations of net-baryon, net-charge, and net-strangeness
multiplicity distributions have been proposed Hatta and Stephanov (2003);
Ejiri  _et al._ (2006); Stephanov (2009); Asakawa  _et al._ (2009); Friman
_et al._ (2011) and measured in ALICE Arslandok (2021), HADES Adamczewski-
Musch  _et al._ (2020), NA61/SHINE Maćkowiak-Pawłowska (2021), and STAR
collaborations Aggarwal  _et al._ (2010); Adamczyk  _et al._ (2014a, b, 2018);
Adam  _et al._ (2021); Nonaka (2021); Abdallah  _et al._ (2021a); Luo (2016).
Recently the STAR experiment reported a non-monotonic beam energy dependence
of the fourth-order fluctuations of net-proton multiplicity distributions from
the Beam Energy Scan program (BES-I) at RHIC Adam  _et al._ (2021); Abdallah
_et al._ (2021b). This is consistent with the expectations of critical
fluctuations in QCD-inspired models Schaefer and Wagner (2012); Chen  _et al._
(2016); Shao  _et al._ (2018); Fan  _et al._ (2019); Li  _et al._ (2019);
Mroczek  _et al._ (2021); Fu  _et al._ (2021). More definitive physics
messages will be achieved by enhanced statistics in Beam Energy Scan program
phase II (BES-II, 2019-2021) bes . According to the model calculations, the
value of the fourth-order fluctuations decreases again to the statistical
baseline at further low collision energies Stephanov (2011). To confirm this
the STAR experiment has been carrying out the fixed-target experiment for
$3.0\leq\sqrt{s_{\rm NN}}\;\leq 7.7$ GeV.

At low collision energies, the fixed-target experiment is more efficient than
collider experiment due to much less interaction rate for the latter Luo  _et
al._ (2020). However, the high interaction rate in the fixed-target experiment
increases the probability of pileup collisions. Pileup is a superposition of
two or more single-collision events occurring within a short time interval. It
is straightforward to correct a mean value of an observable for the pileup
effects using known value of the pileup probability. Unfortunately, this is
not the case for higher-order fluctuations. A correction method on higher-
order fluctuations was proposed in Ref. Nonaka  _et al._ (2020). The pileup
correction was tested in simple toy model with pre-defined correction
parameters. It was thus a closure test where the proposed method should always
work perfectly. In real experiments, however, correction parameters are
unknown and need to be determined by simulations. The accuracy of the pileup
correction strongly depends on how precisely the simulations can reflect the
real conditions of the experimental data. To address this, we propose a model
independent way of an unfolding approach to determine the correction
parameters. We study the effect of pileup events on higher-order cumulants in
UrQMD model. Pileup corrections are carried out in both closure test and
Glauber model fits. These results are shown as a function of the pileup
probability and compared to true values in UrQMD. Finally, the unfolding
approach is applied to UrQMD samples for pileup corrections.

This paper is organized as follows. In Sec. III, we introduce the methodology
of pileup corrections on cumulants. Pileup corrections are applied to UrQMD
model in Sec. IV In Sec. V we propose an unfolding approach to extract the
correction parameters in a model independent way. The validity of the method
is checked in UrQMD model. We then summarize the present study in Sec. VI

##  II Moments and cumulants

Let us consider a probability distribution of particle number P(N). The $r$th
order moment is defined by

$\langle N^{r}\rangle=\sum_{N}N^{r}P(N),$ (1)

where the bracket represents the event average. It is convenient to introduce
a moment generating function,

$G(\theta)=\sum_{N}e^{N\theta}P(N)=\langle e^{N\theta}\rangle.$ (2)

Then the $r$th order moment can be expressed as $r$th derivative of
$G(\theta)$

$\langle N^{r}\rangle=\frac{d^{r}}{d\theta^{r}}G(\theta)\Bigl{|}_{\theta=0}.$
(3)

Since moments drastically increase with increasing the order $r$, cumulants
are easier to handle rather than moments. A cumulant generating function is
defined as

$K(\theta)={\rm ln}G(\theta).$ (4)

The $r$th order cumulant is then given by

$C_{r}=\frac{d^{r}}{d\theta^{r}}K(\theta)\Bigl{|}_{\theta=0}.$ (5)

In the rest of paper, moments will mainly appear for the pileup correction,
while results will be discussed based on cumulants.

##  III Pileup correction

Let $P_{m}(N)$ be a probability distribution function to find one event with
$N$ particles at reference multiplicity $m$. Throughout this paper we suppose
that pileup events are formed by independent superposition of two single-
collision events with the probability $\alpha$. Then $P_{m}(N)$ can be
rewritten as

$\displaystyle
P_{m}(N)=(1-\alpha_{m})P_{m}^{\mathrm{t}}(N)+\alpha_{m}P^{\mathrm{pu}}_{m}(N).$
(6)

where $P_{m}^{\mathrm{t}}(N)$ and $P_{m}^{\mathrm{pu}}(N)$ are probability
distribution functions for single collision event and pileup event
respectively. Pileup events at multiplicity $m$ can be decomposed into sub-
pileup events whose multiplicity satisfies $m$ = $i$ \+ $j$. By looping all
possible combinations of $i$ and $j$ which satisfies $m$ = $i$ \+ $j$,
probability distribution function for pileup events is obtained and written as

$\displaystyle
P_{m}^{\mathrm{pu}}(N)=\sum_{i,j}\delta_{m,i+j}w_{i,j}P_{i,j}^{\mathrm{sub}}(N),$
(7)

where $w_{i,j}$ is the probability to observe a sub-pileup event among all
pileup events at multiplicity $m$. $i$ and $j$ commutes in $w_{i,j}$ which
gives $w_{i,j}$ = $w_{j,i}$. Exhausting all combinations of $i$ and $j$ there
should be

$\displaystyle\sum_{i,j}\delta_{m,i+j}w_{i,j}=1.$ (8)

We also consider a multiplicity distribution $T(m)$ used for centrality
determination. The pileup events at the $m$th multiplicity bin are then
decomposed into two sub-pileup events which satisfies $m=i+j$. Then we get
Nonaka  _et al._ (2020)

$\displaystyle w_{i,j}$ $\displaystyle=$ $\displaystyle\frac{\alpha
T(i)T(j)}{\sum_{i,j}\delta_{m,i+j}\alpha T(i)T(j)},$ (9)
$\displaystyle\alpha_{m}$ $\displaystyle=$
$\displaystyle\frac{\alpha\sum_{i,j}\delta_{m,i+j}T(i)T(j)}{(1-\alpha)T(m)+\alpha\sum_{i,j}\delta_{m,i+j}T(i)T(j)}.$
(10)

The first equation defines the weight of sub-pileup events having
multiplicities $i$ and $j$. The second equation represents the pileup fraction
at $m$th multiplicity bin. True moments are then expressed recursively in
terms of the measured moments at the lower multiplicity bins:

$\displaystyle\langle N^{r}\rangle_{m}^{\rm t}$ $\displaystyle=$
$\displaystyle\cfrac{\langle
N^{r}\rangle_{m}-\alpha_{m}\beta_{m}^{(r)}}{1-\alpha_{m}+2\alpha_{m}w_{m,0}},$
(11)

with

$\displaystyle\beta_{m}^{(r)}$ $\displaystyle=$
$\displaystyle\mu_{m}^{(r)}+\sum_{i,j>0}\delta_{m,i+j}w_{i,j}\langle
N^{r}\rangle_{i,j}^{\rm sub},$ (12)

and

$\displaystyle\mu_{m}^{(r)}=\begin{cases}\displaystyle
2w_{m,0}\sum^{r-1}_{k=0}\binom{r}{k}\langle N^{r-k}\rangle_{0}^{\rm t}\langle
N^{k}\rangle_{m}^{\rm t}&(m>0),\\\
\displaystyle\sum^{r-1}_{k=1}\binom{r}{k}\langle N^{r-k}\rangle_{0}^{\rm
t}\langle N^{k}\rangle_{0}^{\rm t}&(m=0),\end{cases}$ (13)

where $\langle N^{r}\rangle^{t}_{m}$ and $\langle N^{r}\rangle_{m}$ represent
$r$th order true and measured moments at $m$th multiplicity bin, respectively.
As can be seen from Eqs. (9)–(13), the necessary information to perform the
pileup corrections are true multiplicity distribution $T(m)$ and pileup
probability $\alpha$. In real experiments, we can only measure the
multiplicity distributions including pileup events, and therefore, some models
are needed to extract the true multiplicity distribution for single-collision
events. One naive way is to use Glauber and particle production model to fit
the multiplicity distribution as widely used for centrality determination.
Details will be discussed in Sec. IV.

##  IV Pileup correction in UrQMD model

In this section, we demonstrate pileup correction method using UrQMD model.
The package of version 3.4 is used and configured as cascade mode. Around 80
million events are generated with the impact parameter from 0 fm to 15 fm for
Au+Au collisions at $\sqrt{s_{\rm NN}}=$3 GeV. Protons in rapidity and
transverse momentum range of $-0.5<y<0$ and $0.4<p_{T}<2.0$ (GeV/c) are
measured for cumulant calculations of proton multiplicity distributions.
Collision centralities are defined by charged particle multiplicity at the
pseudorapidity range of $|\eta|<1.0$, where protons and antiprotons are
excluded to suppress self-correlations.

Cumulants up to the fourth order are calculated for two kinds of data sets.
One is a raw UrQMD data for single-collision events. This will be a reference
in the simulations without pileup events. The other one is a UrQMD data
including both single-collision and pileup events. They will be referred to as
”raw-UrQMD” and ”pileup-UrQMD” for clarity. In the simulation, the pileup
events are generated via a superposition of two single-collision events. Two
single-collision events in UrQMD are randomly added in terms of particle
multiplicities with probability $\alpha$. Pileup-UrQMD samples are prepared
for $0.1\%<\alpha<10\%$.

Pileup corrections are then performed for pileup-UrQMD data, with correction
parameters determined in two ways as follows.

Figure 1:  (Color online) Charged particle multiplicity distribution in Au+Au
collisions at $\sqrt{s_{\rm NN}}=$3 GeV from UrQMD model. The black circles
are for pileup-UrQMD data with $\alpha=0.5$%. The red squares show raw-UrQMD
data for single-collision events, and a blue line is a Glauber fit to pileup-
UrQMD data. The blue dashed line is pileup-Glauber distribution (pileup +
single-collision) which is simulated using the single-collision distribution
(blue line) extracted by Glauber fit.

Case 1

Closure test
The raw-UrQMD data is used to determine the correction parameters according to
Eqs. (9) and (10) with the know value of $\alpha$. This is a closure test as
the same data sets are used to calculate cumulants and to determine correction
parameters.

Case 2

Glauber fit
The multiplicity distribution in pileup-UrQMD data is fitted by Glauber and
particle production model. In each event, we first determine the number of
sources as:

$N_{\rm sc}=(1-x)N_{\rm part}+xN_{\rm coll},$ (14)

where $N_{\rm part}$ and $N_{\rm coll}$ represent the number of participant
nucleons and binary collisions, respectively, given by the Glauber model.
Final state multiplicity is then produced from $N_{\rm sc}$ according to
negative binomial distributions given by

$\displaystyle
P_{\mu,k}(N)=\cfrac{\Gamma(N+k)}{\Gamma(N+1)\Gamma(k)}\cdot\cfrac{(\mu/k)^{N}}{(\mu/k+1)^{N+k}},$
(15)

where $\mu$ is a mean value of particles generated from one source, and $k$
corresponds to an inverse of a width of the distribution. The values of
parameters employed are $\mu$ = 0.304, $k$ = 3.6 and $x$ = 0.11.

Figure 1 shows multiplicity distribution in UrQMD model. In top panel red
squares represent the raw-UrQMD data, and black circles are the pileup-UrQMD
data with $\alpha=0.5$%. A blue line is a Glauber fit to the pileup-UrQMD
data. Parameters for pileup corrections can be constructed from blue and red
lines in Fig. 1 according to Eqs. 9 and 10 with known value of $\alpha$. In
bottom panel of Figure 1, the black line represents a ratio of multiplicity
distributions between raw-UrQMD data and Glauber fit. With the single-
collision distribution extracted by Glauber fit, we then constructed pileup
events and scan pileup fraction $\alpha$ to find a best fit to Pileup-UrQMD
distribution. The best fit we found is shown as blue dashed line (pileup-
Glauber) in Fig. 1 with $\alpha=0.1$% which is underestimated because of
imperfect Glauber fit to pileup-UrQMD.

Figure 2 shows cumulants and their ratios up to the fourth order as a function
of pileup probability $\alpha$. The black points are results for raw-UrQMD and
black squares are for pileup-UrQMD data. It is found that the results of
pileup-UrQMD deviate from raw-UrQMD with increasing the pileup probability. It
is noteworthy that the value of $C_{4}/C_{2}$ changes sign from negative to
positive then goes above unity with increasing the pileup probability. Results
of pileup corrections for pileup-UrQMD data are shown in red stars for the
closure test, and blue circles for the Glauber fit. It is found that the
pileup correction works well for the wide range of pileup probability up to
10%. On the other hand, the pileup correction using the Glauber fit is not
sufficient with $\alpha>1.0$%. This would be because the fit quality is too
poor to correct cumulants. We have tried our best to tune the parameters of
Glauber and two-component models to fit the multiplicity distribution in UrQMD
data, but could not fit the distribution perfectly. Complicated processes of
particle productions implemented in UrQMD may not be described by the two-
component model.

Figure 2:  (Color online) Proton cumulants and ratios up to fourth-order as a
function of pileup level at top 5 % in Au+Au collisions at $\sqrt{s_{NN}}$=3
GeV from UrQMD model. The black points show results from pure UrQMD model
without any pileup simulation and the black squares are results from pileup
simulated UrQMD data but are not pileup corrected. The red stars and blue
circles are pileup corrected results which used true multiplicity distribution
or Glauber fitted distribution as inputs respectively.

##  V Unfolding

As discussed in Sec. IV, pileup corrections depend on how one can precisely
extract the true multiplicity distribution for single-collision events. The
issue is that the Glauber and particle production models, which are commonly
used for centrality determination, cannot fit even the UrQMD data perfectly. A
model independent way is necessary to make sure the quality of pileup
corrections.

In this section, we propose a model independent way of an unfolding approach
to determine the parameters for pileup corrections. The idea was originally
developed to reconstruct the particle multiplicity distributions in terms of
non-binomial detector efficiencies Esumi  _et al._ (2021). We found that the
similar methodology is applicable in pileup corrections.

Figure 3:  Flowcharts in unfolding to extract the true multiplicity
distribution. The dotted arrows show the procedures repeated for iterations.
Figure 4:  (Left) Correlations between two independence multiplicity
distribution in pileup events. (Right) Response matrices for 0th, 50th, and
100th iteration.  Figure 5:  (Top) Multiplicity distributions for UrQMD,
Glauber fit, and MC samples at 100th iteration. (Bottom) Difference between
UrQMD and MC samples as a function of multiplicity. Left-hand side panels are
for the true coordinates, while right-hand side panels are for the measured
coordinates. The range of the x-axis is limited from 10 to 70 for illustration
purpose.

Figure 3 depicts a flowchart for the unfolding procedure. In real experiments,
the multiplicity distribution is measured with the pileup events on top of the
single-collision events as shown in the left row in Fig. 3. ”True” and
”Measured” in the figure represent the multiplicity distributions for single-
collisions and inclusive distributions for both single-collisions and pileup
events. They are labeled as ”(a) UrQMD experiment true” and ”(b) UrQMD
experiment measured”, respectively. Both are related via a numerical process
to generate pileup events called ”Pileup Filter”, which is defined as an
independent superposition of two single-collision events with probability
$\alpha$ for simplicity. Similarly, we suppose Monte-Carlo samples labeled as
”(c) toy-MC true” and ”(d) toy-MC measured”. They are also related via the
same pileup filter between (a) and (b). In the rest of this paper, samples in
the top row will be referred to as ”true coordinates”, while the bottom row
will be ”measured coordinates”. The goal of the unfolding approach is to
reconstruct (a) starting from (c). Detailed procedures are shown below:

0

Generate (a) UrQMD-experiment and (b) UrQMD-measured samples. They correspond
to raw-UrQMD and pileup-UrQMD distributions in Sec. IV.

1

Generate a (c) toy-MC distribution based on the Glauber model discussed in
Sec. IV.

2

The pileup filter is applied to (c) to get (d) toy-MC measured distribution.

3

During the MC process from 1 to 2, we compute the reversed response matrices,
${\cal R}$, numerically as shown in Fig. 4. Note that any inversion procedure
is not necessary here.

4

The correction function is determined by subtracting (b) from (d). It
represents the difference between UrQMD-experiment and toy-MC distributions in
the measured coordinates. See lower panels in Fig. 5.

5

The response matrix ${\cal R}$ is multiplied to (f) to get (e) the correction
functions in the true coordinates.

6

By adding (e) to (c), the toy-MC distribution is modified to be closer to (a).

7

Repeat 1–6 until the correction functions become close enough to zero.

The response matrices in 3 are defined as

$T(i,j)=\sum_{i,j}\delta_{m;i,j}{\cal R}(i,j;m){\tilde{T}}(m),$ (16)

where ${\tilde{T}}(m)$ represents the probability distribution function of
multiplicity in pileup events at measured coordinates, and $T(i,j)$ is a
correlation between two multiplicities which forms pileup events. $P(i,j)$ and
R(i,j;m) are shown in Fig. 4. Distributions in Fig. 4-(b) are projections of
Fig. 4-(a) onto a diagonal plane for $m=10$, $50$, and $100$ with $m=i+j$. The
response matrices relate the multiplicity $m$ observed in pileup events at
measured coordinates and their original multiplicities from two single-
collision events, $i$ and $j$ at true coordinates. Note that the response
matrices are determined during the numerical process of the pileup filter for
each iteration. Figure 5 shows multiplicity distributions and correction
functions as a function of multiplicity for true and measured coordinates,
respectively. The initial distribution of MC samples are taken from the best
fit of the Glauber model to the UrQMD-experiment distribution in the true
coordinates. Nevertheless, there are large differences from the UrQMD-
experiment distribution as can be seen in the correction functions. After 100
iterations, the correction functions are found to be flat, which indicates
that the multiplicity distribution for MC samples are successfully unfolded to
UrQMD-experiment distributions. The resulting multiplicity distribution for
the true coordinates can be used to determine the parameters for pileup
corrections according to Eq. 9.

Figure 6: (Color online) Centrality dependence of proton cumulant ratios up to
fourth-order as a function of centrality in Au+Au collisions at $\sqrt{s_{\rm
NN}}=$3 GeV from UrQMD model. The label on x-axis from 0 to 8 represents 0-5%,
5-10%, 10-20%, …, 70-80% centralities, respectively. The black points are
results from raw-UrQMD model and black open squares are results from pileup-
UrQMD data. The red stars in upper panels are pileup corrected results which
utilized correction parameters extracted from the unfolding approach, and in
lower panels the red stars and black squares represent ratios of corrected
results and pileup-UrQMD results over raw-UrQMD results respectively.

In our simulations the preseted value of pileup probability $\alpha$ is used
for the unfolding approach. In real experiments, the pileup probability can be
basically calculated from the beam rates and thickness of the target material.
To determine this more precisely, the unfolding approach needs to be repeated
by varying the pileup probability to find the best parameter which yields the
smallest values of $\chi^{2}/{\rm NDF}$.

Let us then move to the pileup corrections on cumulants. The multiplicity
distribution for the true coordinate after 100 iterations is used to define
the parameters for pileup corrections. Results are shown in Fig. 6 for up to
the fourth-order cumulant as a function of centrality. Due to the effect of
pileup events, the results at the most central collisions deviate from the
true value of cumulants. The results of pileup correction using Glauber fits,
however, still deviate from the true cumulants, which is because the Glauber
fit is not perfect enough to describe the multiplicity distribution in UrQMD,
as discussed in Sec. IV. We then apply pileup corrections with correction
parameters determined by the unfolding approach. The results are consistent
with true values of cumulants in the most central collisions. Therefore, it is
concluded that our unfolding approach works well to determine the correction
parameters for pileup corrections. In this work, we simulated pileup events by
the superposition of two single-collision events, in fact the pileup events
merged from more than two single-collision events can be also studied. In the
unfolding approach the MC samples are taken from the best fit of the Glauber
model to the UrQMD-experiment distributions. In principle, MC samples can
start from any distributions like a flat distribution, but we propose to start
from the distribution close to the experimental data to avoid possible
systematics on the initial conditions of the MC samples.

##  VI Summary

We perform detailed studies on the pileup correction for higher-order
cumulants within UrQMD model. Two methods have been applied to determine the
correction parameters and the results of pileup corrections are shown as a
function of pileup probability. It is found that the pileup correction works
well if the correction parameters are determined by the true multiplicity
distribution of single-collision events and the known value of the pileup
probability whereas the correction failed once the Glauber model fit is used
to determine the parameters. This would be because the model of particle
productions are different from those implemented in UrQMD model, thus the fits
to the multiplicity distribution does not work perfectly. The study indicates
that the pileup correction in real experiments strongly depend on the model
and the robustness of the pileup corrections are sensitive to the multiplicity
distribution of single-collision events extracted from the measured
multiplicity distributions. To resolve the issue, we proposed an unfolding
approach to extract the multiplicity distribution of single-collision events
in a model independent way. The validity of this approach is verified in the
UrQMD simulation, in which the pileup correction can successfully recover the
true values of cumulants. Therefore, the model independent unfolding approach
can be used in pileup corrections for the future measurement of higher-order
cumulants in heavy-ion collisions.

##  VII Acknowledgement

We thank D. Cebra, X. Dong, S. Esumi, S. Heppelmann and N. Xu for stimulating
discussions. This work was supported by the National Key Research and
Development Program of China (Grant No. 2020YFE0202002 and 2018YFE0205201),
the National Natural Science Foundation of China (Grant No. 12122505, 11890711
and 11861131009), China Scholarship Council No. 201906770055 and JSPS KAKENHI
Grant No. 19H05598.

## References

  * Gupta  _et al._ (2011) S. Gupta _et al._ , Science 332, 1525 (2011), arXiv:1105.3934 [hep-ph] .
  * Luo and Xu (2017) X. Luo and N. Xu, Nucl. Sci. Tech. 28, 112 (2017), arXiv:1701.02105 [nucl-ex] .
  * Bzdak  _et al._ (2020) A. Bzdak _et al._ , Phys. Rept. 853, 1 (2020), arXiv:1906.00936 [nucl-th] .
  * Stephanov (2004) M. A. Stephanov, Prog. Theor. Phys. Suppl. 153, 139 (2004), arXiv:hep-ph/0402115 .
  * Fukushima and Hatsuda (2011) K. Fukushima and T. Hatsuda, Rept. Prog. Phys. 74, 014001 (2011), arXiv:1005.4814 [hep-ph] .
  * Shi  _et al._ (2014) C. Shi _et al._ , JHEP 07, 014 (2014), arXiv:1403.3797 [hep-ph] .
  * Ding  _et al._ (2015) H. Ding, F. Karsch, and S. Mukherjee, Int. J. Mod. Phys. E 24, 1530007 (2015), arXiv:1504.05274 [hep-lat] .
  * Gao and Liu (2016) F. Gao and Y. Liu, Phys. Rev. D 94, 076009 (2016), arXiv:1607.01675 [hep-ph] .
  * Bazavov  _et al._ (2017) A. Bazavov _et al._ , Phys. Rev. D 95, 054504 (2017), arXiv:1701.04325 [hep-lat] .
  * Fischer (2019) C. S. Fischer, Prog. Part. Nucl. Phys. 105, 1 (2019), arXiv:1810.12938 [hep-ph] .
  * Fu  _et al._ (2020) W. Fu, J. M. Pawlowski, and F. Rennecke, Phys. Rev. D 101, 054032 (2020), arXiv:1909.02991 [hep-ph] .
  * Isserstedt  _et al._ (2019) P. Isserstedt _et al._ , Phys. Rev. D 100, 074011 (2019), arXiv:1906.11644 [hep-ph] .
  * Gao and Pawlowski (2020) F. Gao and J. M. Pawlowski, Phys. Rev. D 102, 034027 (2020), arXiv:2002.07500 [hep-ph] .
  * Bernhardt  _et al._ (2021) J. Bernhardt _et al._ , Phys. Rev. D 104, 074035 (2021).
  * Borsanyi  _et al._ (2021) S. Borsanyi _et al._ , (2021), arXiv:2108.09213 [hep-lat] .
  * Hatta and Stephanov (2003) Y. Hatta and M. A. Stephanov, Phys. Rev. Lett. 91, 102003 (2003), [Erratum: Phys.Rev.Lett. 91, 129901 (2003)], arXiv:hep-ph/0302002 .
  * Ejiri  _et al._ (2006) S. Ejiri, F. Karsch, and K. Redlich, Phys. Lett. B 633, 275 (2006), arXiv:hep-ph/0509051 .
  * Stephanov (2009) M. A. Stephanov, Phys. Rev. Lett. 102, 032301 (2009), arXiv:0809.3450 [hep-ph] .
  * Asakawa  _et al._ (2009) M. Asakawa, S. Ejiri, and M. Kitazawa, Phys. Rev. Lett. 103, 262301 (2009), arXiv:0904.2089 [nucl-th] .
  * Friman  _et al._ (2011) B. Friman, F. Karsch, K. Redlich, and V. Skokov, Eur. Phys. J. C 71, 1694 (2011), arXiv:1103.3511 [hep-ph] .
  * Arslandok (2021) M. Arslandok (ALICE), Nucl. Phys. A 1005, 121979 (2021), arXiv:2002.03906 [nucl-ex] .
  * Adamczewski-Musch  _et al._ (2020) J. Adamczewski-Musch _et al._ (HADES), Phys. Rev. C 102, 024914 (2020), arXiv:2002.08701 [nucl-ex] .
  * Maćkowiak-Pawłowska (2021) M. Maćkowiak-Pawłowska (NA61/SHINE), Nucl. Phys. A 1005, 121753 (2021), arXiv:2002.04847 [nucl-ex] .
  * Aggarwal  _et al._ (2010) M. M. Aggarwal _et al._ (STAR), Phys. Rev. Lett. 105, 022302 (2010), arXiv:1004.4959 [nucl-ex] .
  * Adamczyk  _et al._ (2014a) L. Adamczyk _et al._ (STAR), Phys. Rev. Lett. 113, 092301 (2014a), arXiv:1402.1558 [nucl-ex] .
  * Adamczyk  _et al._ (2014b) L. Adamczyk _et al._ (STAR), Phys. Rev. Lett. 112, 032302 (2014b), arXiv:1309.5681 [nucl-ex] .
  * Adamczyk  _et al._ (2018) L. Adamczyk _et al._ (STAR), Phys. Lett. B 785, 551 (2018), arXiv:1709.00773 [nucl-ex] .
  * Adam  _et al._ (2021) J. Adam _et al._ (STAR), Phys. Rev. Lett. 126, 092301 (2021), arXiv:2001.02852 [nucl-ex] .
  * Nonaka (2021) T. Nonaka (STAR), Nucl. Phys. A 1005, 121882 (2021), arXiv:2002.12505 [nucl-ex] .
  * Abdallah  _et al._ (2021a) M. Abdallah _et al._ (STAR), Phys. Rev. Lett. 127 (2021a), 10.1103/PhysRevLett.127.262301.
  * Luo (2016) X. Luo, Nucl. Phys. A 956, 75 (2016), arXiv:1512.09215 [nucl-ex] .
  * Abdallah  _et al._ (2021b) M. Abdallah _et al._ (STAR), Phys. Rev. C 104, 024902 (2021b), arXiv:2101.12413 [nucl-ex] .
  * Schaefer and Wagner (2012) B. J. Schaefer and M. Wagner, Phys. Rev. D 85, 034027 (2012), arXiv:1111.6871 [hep-ph] .
  * Chen  _et al._ (2016) J.-W. Chen, J. Deng, H. Kohyama, and L. Labun, Phys. Rev. D 93, 034037 (2016), arXiv:1509.04968 [hep-ph] .
  * Shao  _et al._ (2018) G. Shao, Z. Tang, X. Gao, and W. He, Eur. Phys. J. C 78, 138 (2018), arXiv:1708.04888 [hep-ph] .
  * Fan  _et al._ (2019) W. Fan, X. Luo, and H. Zong, Chin. Phys. C 43, 033103 (2019), arXiv:1702.08674 [hep-ph] .
  * Li  _et al._ (2019) Z. Li, K. Xu, X. Wang, and M. Huang, Eur. Phys. J. C 79, 245 (2019), arXiv:1801.09215 [hep-ph] .
  * Mroczek  _et al._ (2021) D. Mroczek _et al._ , Phys. Rev. C 103, 034901 (2021), arXiv:2008.04022 [nucl-th] .
  * Fu  _et al._ (2021) W. Fu _et al._ , Phys. Rev. D 104, 094047 (2021).
  * (40) BES-II White Paper (STAR Note): https://drupal.star.bnl.gov/STAR/starnotes/public/sn0598.
  * Stephanov (2011) M. Stephanov, Phys. Rev. Lett. 107, 052301 (2011), arXiv:1104.1627 [hep-ph] .
  * Luo  _et al._ (2020) X. Luo, S. Shi, N. Xu, and Y. Zhang, Particles 3, 278 (2020), arXiv:2004.00789 [nucl-ex] .
  * Nonaka  _et al._ (2020) T. Nonaka, M. Kitazawa, and S. Esumi, Nucl. Instrum. Meth. A 984, 164632 (2020), arXiv:2006.15809 [physics.data-an] .
  * Esumi  _et al._ (2021) S. Esumi, K. Nakagawa, and T. Nonaka, Nucl. Instrum. Meth. A 987, 164802 (2021), arXiv:2002.11253 [physics.data-an] .