11institutetext:  Instituto de Física, Universidade Federal Fluminense, Av.
Gal. Milton Tavares de Souza s/n, Gragoatá, 24210-346 Niterói, Rio de Janeiro,
Brazil
Departamento de Física, Universidade Federal de São Carlos, P.O. Box 676,
13565-905, São Carlos, São Paulo, Brazil
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
Quantum mechanics Foundations of quantum mechanics Thermodynamics

# Charging power and stability of always-on transitionless driven quantum
batteries

Luiz F. C. Moraes∗ 11 Andreia Saguia† 11 Alan C. Santos‡ 22 Marcelo S.
Sarandy& 111122

###### Abstract

The storage and transfer of energy through quantum batteries are key elements
in quantum networks. Here, we propose a charger design based on transitionless
quantum driving (TQD), which allows for inherent control over the battery
charging time, with the speed of charging coming at the cost of the internal
energy available to implement the dynamics. Moreover, the TQD-based charger is
also shown to be locally stable, which means that the charger can be
disconnected from the quantum battery (QB) at any time after the energy
transfer to the QB, with no fully energy backflow to the charger. This
provides a highly charged QB in an always-on asymptotic regime. We illustrate
the robustness of the QB charge against time fluctuations and the full control
over the evolution time for a feasible TQD-based charger.

###### pacs:

03.65.-w

###### pacs:

03.65.Ta

###### pacs:

05.70.-a

##  1 Introduction

Quantum batteries (QBs) are quantum devices able to both temporarily store and
then transfer energy [1, 2, 3, 4]. They can potentially benefit from
entanglement and other quantum correlations to charge faster than conventional
classical batteries. QBs are potentially relevant as fuel for other quantum
devices and, more generally, for boosting the development of quantum networks.
A number of distinct experimental architectures have been proposed to
implement QBs, such as spin systems [5], quantum cavities [6, 7, 8, 9],
superconducting transmon qubits [10], quantum oscillators [11, 12] and spin
batteries in quantum dots [13].

A key challenge to implement QBs is to ensure the stability of the charging
and energy transfer process. A QB is defined as stable if it can be charged,
with respect to a reference state, with no energy backflow (spontaneous
discharge) at future times. A solution for the stability problem based on the
adiabatic theorem has been introduced in Ref. [14]. As a first step, it is
introduced the energy current (EC) operator, whose expectation value
quantifies the rate of internal energy transferred from the charger to the QB.
By taking the QB in pure states both at initial and final evolution times, the
expectation value of the EC operator at the final time will be the time
derivative of the extractable work from the QB [14]. Then, it has been shown
that the stability is achieved at the adiabatic time scale $\tau_{a}$ by
initially preparing the charger-QB system in an eigenstate $|p_{0}\rangle$ of
the total Hamiltonian acting on the charger-QB Hilbert space. Under such a
condition, the expectation value of the EC operator adiabatically vanishes,
preventing energy backflow to the charger [14]. Therefore, the QB can be
manufactured to charge or discharge its stored energy at any fixed total time
$\tau_{c}$ such that $\tau_{c}>\tau_{a}$. Even though this provides a general
adiabatic-based approach to yield stable QBs, which is independent of
architectures and corresponding Hamiltonians, it turns out to be still a
partial solution for the stability problem. First of all, it requires the
validity of the adiabatic approximation, which imposes constraints on the
energy gap dynamics and makes the system potentially more susceptible to
decoherence due to the competition between the adiabatic and the relaxation
time scales. Moreover, Ref. [14] solves a particular stability problem, which
we define here as the global stability problem. This means that stability is
ensured by manufacturing the QB to complete the charge at a specific global
time $\tau_{c}$ (such that $\tau_{c}>\tau_{a}$). The choice of $\tau_{c}$ is
arbitrary as long as it is greater than $\tau_{a}$. However, the Hamiltonian
dynamics is engineered to occur from the initial time $t=0$ up to
$t=\tau_{c}$, with no time fluctuation. Therefore, global stability does not
account for a possible energy backflow if the charger is still connected to
the QB after the time $\tau_{c}$. The ability of keeping the QB charged with
no energy revivals in the charger in an always-on asymptotic regime is what we
define as local stability, which is not dealt with in the solution proposed in
Ref. [14].

In this work, we propose a solution to the adiabatic time constraint as well
as the local stability problem. The adiabatic hypothesis can be relaxed here
by implementing transtionless quantum driving (TQD) [15, 16, 17]. TQD is a
useful technique to mimic adiabatic quantum tasks at arbitrary finite time,
providing a shortcut to adiabaticity. It has been applied for speeding up
quantum gate Hamiltonians [18, 19, 20, 21], heat engines in quantum
thermodynamics [22], quantum information processing [23, 24, 25], among other
applications (see e.g., Ref. [26] and references there in). Naturally, there
is a commitment between the speed of the evolution and the internal energy
required to drive the dynamics. Moreover, the implementation of TQD usually
requires many-body interaction terms in the Hamiltonian, whose number may
exponentially grow with the size of the system. As we will show, the
Hamiltonian complexity is tractable here, requiring up to three-body
interactions. Concerning the local stability problem, we solve it through a
suitable non-linear Hamiltonian interpolation, which yields local stability
either in the adiabatic or in the TQD scenario. These non-linear
interpolations keep a non-vanishing gap throughout the evolution while
ensuring an asymptotically fueled QB. The extractable energy will be
considered from the point of view of the ergotropy of the QB, which turns out
to show robustness against time variation.

##  2 Ergotropy and QB stability

Let us begin by discussing some preliminary concepts, namely, the work
extractable from a charger and the distinct kinds of stability behaviors in
the charging and discharging energy transfer processes.

###  2.1 Ergotropy and internal energy

The usefulness of the energy tranferred by a charger can be quantified in
terms of the ergotropy, which is defined as the maximum energy that can be
extractable as work from a quantum system through unitary operations $V$.
Ergotropy has been proposed by Allahverdyan et al. [27], reading

$\displaystyle\mathcal{E}=\text{Tr}\left\\{\rho
H_{0}\right\\}-\min_{V\in\mathcal{V}}\text{Tr}\left\\{V\rho
V^{\dagger}H_{0}\right\\},$ (1)

where $\rho$ denotes the charger state, $V$ is a unitary operator belonging to
the set $\mathcal{V}$ of all possible unitary operations over the system, and
$H_{0}$ is the reference Hamiltonian that defines the internal battery
structure (defining the charger energy states). Without loss of generality, we
can write $H_{0}$ by using its spectral decomposition
$H_{0}\\!=\\!\sum\nolimits_{n=1}^{\mathcal{D}}\epsilon_{n}|\epsilon_{n}\rangle\langle\epsilon_{n}|$,
with $\mathcal{D}$ being the dimension of the Hilbert space and
$\epsilon_{1}\\!\leq\\!\epsilon_{2}\\!\leq\\!\cdots\\!\leq\\!\epsilon_{\mathcal{D}}$.
Therefore, the empty battery energy is taken as $\epsilon_{1}$ and the full
energy value as $\epsilon_{\mathcal{D}}$, so that the charger has an energy
capacity given by
$\mathcal{E}_{\text{max}}\\!=\\!\epsilon_{\mathcal{D}}-\epsilon_{1}$. Let us
now consider the spectral decomposition of the charger density operator
$\rho\\!=\\!\sum\nolimits_{n=1}^{\mathcal{D}}\varrho_{n}|\varrho_{n}\rangle\langle\varrho_{n}|$,
where $|\varrho_{n}\rangle$ are eigenvectors of $\rho$ with eigenvalues
$\varrho_{n}$, so that
$\varrho_{1}\\!\geq\\!\varrho_{2}\\!\geq\\!\cdots\\!\geq\\!\varrho_{\mathcal{D}}$.
By taking this arrangement order, the minimization in Eq. (1) is attained for
$V=\sum\nolimits_{n}|\epsilon_{n}\rangle\langle\varrho_{n}|$, yielding

$\mathcal{U}_{0}=\min_{V\in\mathcal{V}}\text{Tr}\left\\{V\rho
V^{\dagger}H_{0}\right\\}=\sum_{n=1}^{\mathcal{D}}\varrho_{n}\epsilon_{n}.$
(2)

Notice that $\mathcal{U}_{0}$ arises from the operation $V$ that drives $\rho$
to a diagonal form in the energy eigenbasis. Then, Eq. (1) can be rewritten in
a more convenient way as [27, 28]

$\displaystyle\mathcal{E}=\sum_{i=1}^{\mathcal{D}}\sum_{n=1}^{\mathcal{D}}\varrho_{n}\epsilon_{i}\left(|\langle\varrho_{n}|\epsilon_{i}\rangle|^{2}-\delta_{ni}\right).$
(3)

From Eq. (3), it is worth emphasizing that not all the initial internal energy
can be effectively transferred as work (see, e.g., Refs. [29, 30]). In fact,
by writing the internal energy $\mathcal{U}(\rho,H_{0})$ for the reference
Hamiltonian $H_{0}$, we have

$\displaystyle\mathcal{U}(\rho,H_{0})=\text{Tr}\left\\{\rho
H_{0}\right\\}=\sum_{i=1}^{\mathcal{D}}\sum_{n=1}^{\mathcal{D}}\varrho_{i}\epsilon_{n}|\langle\varrho_{n}|\epsilon_{i}\rangle|^{2},$
(4)

which means that

$\displaystyle\mathcal{E}=\mathcal{U}({\rho,H_{0}})-\mathcal{U}_{0}.$ (5)

Therefore, since ergotropy is positive for any process (by construction), the
internal energy of the system needs to satisfy
$\mathcal{U}({\rho,H_{0}})\\!>\\!\mathcal{U}_{0}$ for an amount of extractable
work to be available in the QB setup.

###  2.2 Stability criteria

A quantum system exhibiting discrete constant eigenvalues obeys the quantum
version of the recurrence theorem of Poincaré [31], which establishes that, if
$|\psi(t_{0})\rangle$ is the system state vector at time $t_{0}$ and
$\epsilon$ is any positive number, then a future time $T$ will exist such that
$\|\,|\psi(T)\rangle-|\psi(t_{0})\rangle\,\|^{2}\\!<\\!\epsilon$, with
$\|\,|\phi\rangle\,\|$ denoting the vector norm
$\|\,|\phi\rangle\,\|\\!=\\!\sqrt{\langle\phi|\phi\rangle}$. Therefore, any QB
obeying the quantum recurrence theorem will be unstable, ı.e. subjected to
spontaneous charging and discharging processes [10, 14, 30]. As a strategy to
avoid spontaneous energy transfer between a charger and the QB, we may
consider suitable time-dependent evolutions, such as the adiabatic dynamics
[32, 33]. Indeed, this leads to the possibility of designing stable QBs [10,
14].

In this scenario, Ref. [14] introduced an adiabatic-based solution for a
specific stability problem, which we define here as the global stability
problem. Global stability requires the QB to complete the charge at a freely
programmed global time $\tau_{c}$, which is required to be large enough to
obey the adiabatic approximation [34, 35]. The dynamics is then assumed to
globally occur from the initial time $t=0$ up to $t=\tau_{c}$, with absolutely
no time fluctuation. Therefore, global stability does not account for a
possible energy backflow if the charger is still connected to the QB after the
time $\tau_{c}$.

In order to prevent energy revivals in the charger in the always-on time
regime, we also introduce here the concept of local stability. By taking
$\epsilon$ as a positive number, local stability is defined through the
requirement that the QB ergotropy $\mathcal{E}_{{QB}}(t)$ is such that
$\mathcal{E}_{{QB}}(t)>\epsilon$ after a finite charging time $\tau_{c}$.
Notice then that local stability refers to the prevention of full energy
backflow to the charger locally in time, as it would be expected from the
quantum recurrence theorem. Naturally, the large the value of $\epsilon$ the
higher the performance of the locally stable QB. Let us consider that
systematic errors or an inherent longer turned-on charger-QB connection could
lead us to an evolution time $\Delta t_{\text{real}}\\!=\\!\tau_{c}+\Delta t$.
Then, the battery ergotropy can deviate from its maximum value as
$\mathcal{E}_{\text{real}}(\Delta
t)\\!=\\!\mathcal{E}_{\text{max}}-|\Delta\mathcal{E}_{\text{dev}}(\Delta t)|$.
By rearranging this expression, we have

$\mathcal{E}_{\text{real}}(\Delta
t)/\mathcal{E}_{\text{max}}=1-\eta_{\text{ls}}(\Delta t),$ (6)

where $\eta_{\text{ls}}(\Delta t)\equiv|\Delta\mathcal{E}_{\text{dev}}(\Delta
t)|/\mathcal{E}_{\text{max}}\\!\in\\![0,1]$ can be understood as a local
stability coefficient for the QB. There are situations in which the charge
$\mathcal{E}_{\text{max}}$ is not achieved, but we have an asymptotic amount
of ergotropy $\mathcal{E}_{\text{asy}}\\!<\\!\mathcal{E}_{\text{max}}$ [30].
In this case, Eq. (6) can be generalized by replacing
$\mathcal{E}_{\text{max}}$ for $\mathcal{E}_{\text{asy}}$
($\mathcal{E}_{\text{max}}$ is merely a chosen normalization). As an example,
consider the case where one can use dark states to stabilize open system QBs
[30]. Due to the effect of decoherence, the maximum ergotropy is not achieved.
However, we have an asymptotic value $\mathcal{E}_{\text{asy}}\\!>\\!0$ for
which, after some instant $t\\!=\\!\tau_{\text{c}}$, one gets
$\eta_{\text{ls}}(\Delta t)\\!=\\!0$, for any $\Delta t\\!>\\!0$, since
$\Delta\mathcal{E}_{\text{dev}}(\Delta t>0)\\!=\\!0$. In this situation, one
has an optimal locally stable charging process. It is important mentioning
that there are situations for which stability is achieved only for some finite
interval $\Delta t$, such as recently reported in non-Markovian [29] and
Floquet engineered QBs [36].

##  3 Stability of the adiabatic ergotropy transfer

In this section, we will discuss the energy transfer process of QBs driven by
adiabatic Hamiltonians. Let us consider the spectral decomposition of a time-
dependent Hamiltonian
$H_{\text{ad}}(t)=\sum\nolimits_{n=1}^{\mathcal{D}}E_{n}(t)|n(t)\rangle\langle
n(t)|$. We prepare the initial state $|\psi(0)\rangle$ of the system in a
single non-degenerate eigenstate $|n(0)\rangle$ of $H_{\text{ad}}(0)$. By
assuming an adiabatic evolution, the evolved state reads
$|\psi(t)\rangle\\!=\\!e^{i\theta_{\text{ad}}(t)}|n(t)\rangle$, where the
adiabatic phase is $\theta_{\text{ad}}(t)\\!=\\!-E_{n}(t)/\hbar+i\langle
n(t)|\dot{n}(t)\rangle$. The density matrix is then
$\rho_{\text{ad}}(t)\\!=\\!|n(t)\rangle\langle n(t)|$, with no contribution of
the adiabatic phase.

The charger model adopted here is composed of a two-qubit cell (allowing for
entanglement), whose internal structure is defined by the Hamiltonian
$H_{0}\\!=\\!\hbar\omega(\sigma_{z}^{(1)}+\sigma_{z}^{(2)})$, with energy
basis given by $\sigma_{z}|n\rangle\\!=\\!(-1)^{1-n}|n\rangle$ for
$n\in\\{0,1\\}$. Under this configuration, each qubit has a storable energy
capacity $\mathcal{E}_{\text{max}}^{\text{qubit}}\\!=\\!2\hbar\omega$, with
the corresponding full charge state $|1\rangle$. Consequently, one can define
the charger empty and full charge states as
$|\text{emp}\rangle\\!=\\!|00\rangle$ and
$|\text{full}\rangle\\!=\\!|11\rangle$, respectively, so that for the charger
in the state $|\text{full}\rangle$ the ergotropy reads
$\mathcal{E}_{\text{max}}^{\text{cell}}\\!=\\!2\mathcal{E}_{\text{max}}^{\text{qubit}}\\!=\\!4\hbar\omega$.
The initial charger state is taken as in Ref. [14], which is given by the
singlet state

$\displaystyle|\psi(0)\rangle_{{\mathcal{C}}}=\frac{|0\rangle_{{\mathcal{C}}_{1}}|1\rangle_{{\mathcal{C}}_{2}}-|1\rangle_{{\mathcal{C}}_{1}}|0\rangle_{{\mathcal{C}}_{2}}}{\sqrt{2}},$
(7)

with ${\mathcal{C}}_{i}$ denoting the $i$-th qubit of the charger
($i\in\\{1,2\\}$). The maximally entangled state in Eq. (7) is responsible for
the quantumness of the charger, allowing for a more efficient energy transfer
to a QB than a factorizable state and for the design of switchable chargers
[14]. The process we will specifically consider here is the ergotropy transfer
from the charger to the QB. The QB will be taken here a single qubit system,
initially in the state $|0\rangle_{\text{{QB}}}$. Then, after some time
$\tau_{\text{c}}$, we expect to transfer the ergotropy initially stored in the
charger to the QB. By denoting
${\mathcal{C}}\\!=\\!{\mathcal{C}}_{1}\otimes{\mathcal{C}}_{2}$ as the charger
Hilbert space, the dynamics is then expected to lead the system from the
initial state
$|\phi(0)\rangle\\!=\\!|\psi(0)\rangle_{{\mathcal{C}}}|0\rangle_{\text{{QB}}}$
to the final state
$|\phi(\tau_{\text{c}})\rangle\\!=\\!|0\rangle_{{\mathcal{C}}_{1}}|0\rangle_{{\mathcal{C}}_{2}}|1\rangle_{\text{{QB}}}$
through a stable charging process of the QB. To this aim, we drive the system
by a piece-wise time-dependent Hamiltonian given by

$\displaystyle
H_{\text{ad}}(t)=\left[1-f(t)\right]H_{\text{ini}}+f(t)\left[1-f(t)\right]H_{\text{inter}}+g(t)H_{\text{fin}},$
(8)

where $f(t)$ and $g(t)$ are real functions satisfying
$f(0)\\!=\\!g(0)\\!=\\!0$ and $f(\tau)\\!=\\!g(\tau)\\!=\\!1$, for some total
evolution time $\tau$, and each time-independent Hamiltonian reads

$\displaystyle H_{\text{ini}}$
$\displaystyle=\hbar\Omega\left[\sigma_{x}^{({\mathcal{C}}_{1})}\sigma_{x}^{({\mathcal{C}}_{2})}+\sigma_{y}^{({\mathcal{C}}_{1})}\sigma_{y}^{({\mathcal{C}}_{2})}\right]\mathbbm{1}_{\text{{QB}}},$
(9) $\displaystyle H_{\text{inter}}$
$\displaystyle=\hbar\Omega\mathbbm{1}^{({\mathcal{C}}_{1})}\left[\sigma_{x}^{({\mathcal{C}}_{2})}\sigma_{x}^{\text{{QB}}}+\sigma_{y}^{({\mathcal{C}}_{2})}\sigma_{y}^{\text{{QB}}}\right],$
(10) $\displaystyle H_{\text{fin}}$
$\displaystyle=\hbar\Omega\left[\sigma_{z}^{({\mathcal{C}}_{1})}\mathbbm{1}^{({\mathcal{C}}_{2})}\sigma_{z}^{\text{{QB}}}+\mathbbm{1}^{({\mathcal{C}}_{1})}\sigma_{z}^{({\mathcal{C}}_{2})}\sigma_{z}^{\text{{QB}}}\right].$
(11)

Figure 1: Ergotropy $\mathcal{E}^{\text{ad}}_{\text{{QB}}}(\tau)$ (as a
multiple of $\mathcal{E}_{\text{max}}^{\text{qubit}}$) of the QB at the end of
the adiabatic evolution (at time $t\\!=\\!\tau$). Global stability is achieved
for $\tau$ at the adiabatic scale.

Each contribution to the Hamiltonian plays a specific role in the energy
transfer process. $H_{\text{ini}}$ allows us to encode the state
$|\phi(0)\rangle$ as one of its eigenstates. The final desired state
$|\phi(\tau_{\text{c}})\rangle$ can then be correctly addressed because it is
an eigenstate of $H_{\text{fin}}$ with identical energy and parity as the
initial state [14]. Concerning $H_{\text{inter}}$, it is introduced here in
order to ensure a non-crossing ground state energy gap.

It is worth mentioning that our choice of the function $f(t)$ may affect the
ergotropy transfer rate and the stability of the battery, since the minimum
instantaneous gap depends on the function $f(t)$. As an illustration, we will
consider three different functions given by
$f_{\text{lin}}(t)\\!=\\!g_{\text{lin}}(t)\\!=\\!t/\tau$,
$f_{\text{sin}}(t)\\!=\\!g_{\text{sin}}(t)\\!=\\!\sin\left(\pi
t/2\tau\right)$, and $f^{3}_{1/3}(t)\\!=\\!g_{1/3}(t)\\!=\\!t/\tau$. The
performance of the functions $f(t)$ and $g(t)$ can now be promptly verified.
By driving the system by the Hamiltonian $H_{\text{ad}}(t)$ and computing the
ergotropy $\mathcal{E}^{\text{ad}}_{\text{{QB}}}(\tau)$ of the QB at the end
of the evolution (at time $t\\!=\\!\tau$), we can analyze the global stability
of the energy transfer process. To this end, the Hamiltonian and reduced
density matrix used in Eq. (5) are given by
$H_{0}^{\text{{QB}}}\\!=\\!\hbar\omega\sigma_{z}$ and
$\rho^{\text{ad}}_{\text{{QB}}}(t)\\!=\\!\text{Tr}_{{\mathcal{C}}}[\rho_{\text{ad}}(t)]$,
respectively. The result is shown in Fig. 1. Notice that, for a given total
time $\tau$, the QB ergotropy is approximately equivalent for any
interpolation scheme, with some advantage for the linear interpolation in most
of regions except the region near $\Omega\tau=10$, where the interpolation via
$f_{1/3}(t)$ and $g_{1/3}(t)$ is a better option. Concerning stability, all
the three interpolations are globally stable at the adiabatic time scale. In
fact, the power of adiabatic chargers is dictated by the energy gap structure
of $H_{\text{ad}}(t)$, while its stability is achieved whenever the adiabatic
dynamics takes place and the initial state is encoded in a single eigenstate
of $H_{\text{ad}}(t)$ [10, 14].

We also consider the instantaneous charge as a function of $t$, for a fixed
total time $\tau$. This is performed to study the local stability for
different choices of the interpolating functions. Again, the ergotropy is
computed for the QB as in Eq. (5). The results are shown in Fig. 2, where we
have set $\Omega\tau\\!=\\!10$. We highlight the intant in which $\Omega
t\\!=\\!10$ (vertical dotted line), so that we can see the efficiency in the
stability of each evolution considered. Notice that the interpolation via
$f_{1/3}(t)$ and $g_{1/3}(t)$ is able to asymptotically ensure local
stability, since it is designed to favor $H_{fin}$ for large $t$. For an
always-on connection, this interpolation prevents full backflow of energy from
the QB to the charger. From the behavior of the ergotropy, it is
straightforward to see that the efficiency $\eta_{\text{ls}}(\Delta t)$ can be
drastically affected by the choice of the interpolating function for
$H_{\text{ad}}(t)$.

##  4 Stability of TQD ergotropy transfer

As previously shown, adiabaticity provides an optimal strategy to achieve
global stability, even though local stability will depend on the interpolation
funtion adopted to drive the system. On the other hand, due to limitations
imposed by the validity conditions of adiabatic dynamics, we inevitably have a
loss of power in the charging/discharging stage. In order to provide an
alternative approach to recover the high-power performance, keeping stability,
we will use transitionless quantum driving (TQD) for speeding up the energy
transfer process.

Figure 2: Ergotropy $\mathcal{E}^{\text{ad}}_{\text{{QB}}}(\tau)$ of the QB
(as a multiple of $\mathcal{E}_{\text{max}}^{\text{qubit}}$) as a function of
$t$, for $\Omega\tau\\!=\\!10$. The horizontal dashed line denotes the maximal
charge, which has been normalized to 1. Local stability strongly depends on
the interpolation function adopted.

As originally proposed, the TQD approach for speeding up the adiabatic
evolution requires the inclusion of an additional term $H_{\text{cd}}(t)$,
called counter-diabatic Hamiltonian, to the original adiabatic Hamiltonian
$H_{\text{ad}}(t)$ so that the dynamics is governed by the TQD Hamiltonian
$H_{\text{tqd}}(t)=H_{\text{ad}}(t)+H_{\text{cd}}(t)$, where

$\displaystyle H_{\text{cd}}(t)=i\hbar\sum_{n}|\dot{n}(t)\rangle\langle
n(t)|+\langle\dot{n}(t)|n(t)\rangle|n(t)\rangle\langle n(t)|,$ (12)

with $|n(t)\rangle$ denoting the set of eigenstates of $H_{\text{ad}}(t)$.
Under this consideration, one can speed up the adiabatic dynamics by the
adiabatic path of $H_{\text{ad}}(t)$, since diabatic transitions are forbidden
due to presence of the term $H_{\text{cd}}(t)$. It is known that we can
arbitrarily speed up the dynamics, given that the $H_{\text{cd}}(t)$ can be
efficiently implemented (or simulated) [18].

In general, the stability of the battery is not a priori guaranteed for TQD
even in case where we have perfect control of the time evolution of
$H_{\text{tqd}}(t)$. In fact, the stability of the adiabatic approach is
obtained because the system state at instant $t\\!=\\!\tau$ is an eigenstate
of $H_{\text{ad}}(t)$, so that the system continuously to evolve to the
corresponding instantaneous eigenstate at latter times [14]. Therefore, we do
not need to turn off the charger-QB coupling. In the TQD case, due to term
$H_{\text{cd}}(t)$, the state of the system at instant $t\\!=\\!\tau$ is not
eigenstate of $H_{\text{tqd}}(t)$, then some dynamics is expected if we do not
turn off the charger-QB interaction. However, in cases where we use the
boundary conditions to get
$H_{\text{cd}}(t\\!=\\!0)\\!=\\!H_{\text{cd}}(t\\!=\\!\tau)\\!=\\!0$, the
stability is recovered. It is known that such boundary condition can be
obtained for a number of evolutions [37, 38, 26], but here we will adopt a
milder assumption, with no boundary conditions on time taken. As we shall see,
even under this assumption, we get considerable degree of stability while
speeding up the ergotropy transfer process.

To study the performance of TQD charging, we consider the shortcut to the
adiabatic evolution driven by the Hamiltonian in Eq. (8). First, let us
consider global stability. To this end, given the counter-diabatic term
$H_{\text{cd}}(t)$, we let the charger-QB evolve under $H_{\text{tqd}}(t)$ and
then we compute the QB ergotropy at $t\\!=\\!\tau$. In Fig. 3 we show the
ergotropy $\mathcal{E}^{\text{tqd}}_{\text{{QB}}}(\tau)$ of the QB after the
TQD evolution, where we choose two different interpolations, namely,
$f_{\text{lin}}(t)$ and $f_{1/3}(t)$ to get a comparison between the adiabatic
and TQD charging performance. While adiabatic charging requires a specific
minimum time $\tau$ to achieve maximum charge (or to get some non-zero amount
of ergotropy), the counter-diabatic energy transfer leads to a maximum charge
state at an arbitrary time interval. This is a consequence of the inverse
engineering approach, where a suitable counter-diabatic contribution allows us
to mimic the adiabatic dynamics at arbitrary $\tau$, as dictated by the
quantum speed limit for counter-diabatic dynamics [18]. On the other hand, it
is worth highlighting that such additional advantage comes at the well-known
cost of some intensity-field demand in the term $H_{\text{cd}}(t)$ in
$H_{\text{tqd}}(t)$ [18].

(a) $f_{\text{lin}}(t)$

(b) $f_{1/3}(t)$

Figure 3: Adiabatic and TQD ergotropy transfer from the charger to the QB (as
a multiple of $\mathcal{E}_{\text{max}}^{\text{qubit}}$) at instant
$t\\!=\\!\tau$, as a function of $\Omega\tau$, for distinct interpolation
functions, namely, $f_{\text{lin}}(t)$ [Fig. (3a)] and $f_{1/3}(t)$ [Fig.
(3b)].

Concerning local stability, we can also show that the conter-diabatic approach
provides a robust charge method even if fluctuations of the total time $\tau$
are taken into account. This is illustrated in Fig. 4, where both adiabatic
and TQD ergotropy transfer from the charger to the QB (as a multiple of
$\mathcal{E}_{\text{max}}$) are plotted as a function of $\Omega t$ for
different total times $\tau$ and for the interpolation function given by
$f_{1/3}(t)$ e $g_{1/3}(t)$. Notice that, by adopting $\Omega\tau=1$, the
adiabatic condition is violated. Again, we highlight the instant $\Omega
t\\!=\\!1$ as a vertical dotted line. In that case, the TQD approach is still
able to mimic adiabaticity and to provide a locally stable ergotropy transfer.
Even though we have the local stability coefficient $\eta_{\textrm{ls}}>0$,
the QB ergotropy is kept in a large value in the asymptotic regime. In the
inset of Fig. 4, we also consider a case of validity of adiabaticity by
setting $\Omega\tau=10$ (vertical dotted line denotes the instant $\Omega
t\\!=\\!10$). In this situation, both the adiabatic method and the TQD
approach provide a locally stable charging process from the charger to the QB.

Figure 4: Adiabatic and TQD ergotropy transfer from the charger to the QB (as
a multiple of $\mathcal{E}_{\text{max}}^{\text{qubit}}$) as a function of
$\Omega t$, for the interpolation function given by $f_{1/3}(t)$ e
$g_{1/3}(t)$. The ergotropy is shown for the case of $\Omega\tau=1$, where the
TQD approach is able to provide a stable charge, while the adiabatic method
fails. The horizontal dashed line denotes the maximal charge, which has been
normalized to 1. In the inset, a slower evolution is adopted by setting
$\Omega\tau\\!=\\!10$, so that both TQD and adiabatic dynamics achieve local
stability.

##  5 Energy cost of the charging process

In this section we briefly discuss the external energy cost of implementing
the TQD charging process in comparison with the adiabatic evolution. It is
known that the standard approach of TQD is more energetically costly than its
adiabatic counterpart [18]. However, as we shall see, it is possible to find a
range of values of the total evolution time for which the TQD charging
provides a beneficial setup. To this end, let us consider the measure
$\Sigma(\tau)$ of energy cost for implementing a time-dependent Hamiltonian
through the interval $t\in[0,\tau]$, which is provided by

$\displaystyle\Sigma(\tau)=\frac{1}{\tau}\int_{0}^{\tau}||H(t)||_{\text{HS}}dt,$
(13)

where $||A||_{\text{HS}}\\!=\\!(\text{Tr}\left\\{AA^{\dagger}\right\\})^{1/2}$
is the Hilbert-Schmidt norm of the operator $A$. We observe that
$\Sigma(\tau)$ has been used in experimental implementations of TQD in order
to quantify the energy cost of simulating adiabatic and TQD evolutions [39,
40]. Moreover, it has recently been shown in Ref. [41] that this quantity
allows for quantifying the thermodynamic cost of implementing the dynamics
driven by an arbitrary Hamiltonian $H(t)$. By explicitly evaluating the
adiabatic and TQD energy costs, we then obtain

$\displaystyle\Sigma_{\text{ad}}$
$\displaystyle=\frac{1}{\tau}\int_{0}^{\tau}||H_{\text{ad}}(t)||_{\text{HS}}dt=\frac{1}{\tau}\int_{0}^{\tau}\sqrt{\sum_{n}E_{n}^{2}(t)}dt,$
(14) $\displaystyle\Sigma_{\text{tqd}}$
$\displaystyle=\frac{1}{\tau}\int_{0}^{\tau}||H_{\text{tqd}}(t)||_{\text{HS}}dt$
$\displaystyle=\frac{1}{\tau}\int_{0}^{\tau}\sqrt{\sum_{n}E_{n}^{2}(t)+\hbar^{2}\mu_{n}(t)}dt,$

where $\mu_{n}(t)=\langle\dot{n}(t)|\dot{n}(t)\rangle-|\langle
n(t)|\dot{n}(t)\rangle|^{2}$ [19]. Therefore, given the parameters $\Omega$
and $\tau$, and the functions $f(t)$ and $g(t)$, the quantities above can be
directly computed. In our case, we consider the linear function
$f_{\text{lin}}(t)$ and perform the integration by numerical methods, with
results shown in Fig. 5. It is also convenient to define the ratio
$\Sigma_{\text{rel}}=\Sigma_{\text{tqd}}/\Sigma_{\text{ad}}$, which quantifies
how larger the energy cost is in the TQD dynamics in comparison with its
adiabatic counterpart. The result for $\Sigma_{\text{rel}}$ is shown in the
inset of Fig. 5.

Figure 5: Energy cost (as a multiple of $\hbar\Omega$) as a function of
$\Omega\tau$ for the implementation of the adiabatic and TQD charging of the
QB. Inset: The relative cost of the TQD process in comparison with the
adiabatic charging. Here we consider the interpolation function
$f_{\text{lin}}(t)$.

As expected, the faster the charging process is, the larger the energy cost to
implement it in a TQD process, as we can see in Fig. 5. However, it is
possible to find a range of values for the charging time $\tau$ for which the
TQD protocol presents better performance with respect to the adiabatic case.
In fact, in the inset of Fig. 5, we consider the time interval between
$\Omega\tau$ and $10\Omega\tau$. As previously shown in Fig. 3a, for this
interval, the TQD process allows us for reaching the full charged state of the
battery, while the adiabatic evolution cannot charge the battery efficiently (
with the maximum being around $80\%$ of
$\mathcal{E}_{\text{max}}^{\text{qubit}}$). Then, by taking into account the
cost of implementing the transitionless evolution, relatively to the adiabatic
dynamics, it is worth supporting the charging process via TQD. For example,
consider the range of values where
$\mathcal{E}^{\text{ad}}_{\text{QB}}<10\%\mathcal{E}_{\text{max}}^{\text{qubit}}$,
as highlighted in the yellow box in Fig. 5 (inset). For this interval of
values for $\Omega\tau$, the maximum value for
$\Sigma_{\text{tqd}}(\Omega\tau)$ corresponds to a mild TQD increasing cost,
with $\Sigma_{\text{tqd}}\\!=\Sigma_{\text{ad}}+\\!5\%\Sigma_{\text{ad}}$.

##  6 Conclusion

We have proposed a solution for the local stability of QBs as well as for the
inherent control of their charging power. This allows for the implementations
of always-on chargers that are able to charge a QB as fast as the energy
available for their internal dynamics. The charger-QB model proposed here is
composed by a set of independent three-qubit cells, which ensures the locality
of the interaction Hamiltonian. The energy resources can be further optimized
by implementing a more general approach of TQD technique [42]. In this
approach, the phase accompanying the state vector that mimics the adiabatic
evolution is set as a free parameter that can be used to adjust the
interaction Hamiltonian for a smoother TQD energy requirement. This is a
certainly interesting point for further analysis. Moreover, it is also
potentially fruitful the investigation of the relationship between the
correlations present in the charger and its performance with respect to the
energy storage and transfer. This may bring insights about possible
generalizations of the charger model to favor the relevant correlations. These
topics are left for future research.

It is worth mentioning that we recently became aware of a work in which the
authors also consider TQD for speeding up the adiabatic charging process of
three-level quantum batteries [43].

###### Acknowledgements.

L. F. C. M. acknowledges financial support from the Coordenação de
Aperfeiçoamento de Pessoal de Nível Superior (CAPES). A.C.S. acknowledges
financial support through the research grant from the São Paulo Research
Foundation (FAPESP) (Grant No 2019/22685-1). M.S.S. acknowledges financial
support from the Conselho Nacional de Desenvolvimento Científico e Tecnológico
(CNPq) (No. 307854/2020-5). This research is also supported in part by CAPES
(Finance Code 001) and by the Brazilian National Institute for Science and
Technology of Quantum Information [CNPq INCT-IQ (465469/2014-0)].

## References

  * [1] Alicki R. Fannes M. Phys. Rev. E872013042123.
  * [2] Hovhannisyan K. V., Perarnau-Llobet M., Huber M. Acín A. Phys. Rev. Lett.1112013240401.
  * [3] Campaioli F., Pollock F. A., Binder F. C., Céleri L., Goold J., Vinjanampathy S. Modi K. Phys. Rev. Lett.1182017150601.
  * [4] Andolina G. M., Keck M., Mari A., Campisi M., Giovannetti V. Polini M. Phys. Rev. Lett.1222019047702.
  * [5] Le T. P., Levinsen J., Modi K., Parish M. M. Pollock F. A. Phys. Rev. A972018022106.
  * [6] Binder F. C., Vinjanampathy S., Modi K. Goold J. New J. Phys.172015075015.
  * [7] Fusco L., Paternostro M. De Chiara G. Phys. Rev. E942016052122.
  * [8] Zhang Y.-Y., Yang T.-R., Fu L. Wang X. Phys. Rev. E992019052106.
  * [9] Ferraro D., Campisi M., Andolina G. M., Pellegrini V. Polini M. Phys. Rev. Lett.1202018117702.
  * [10] Santos A. C., Çakmak B., Campbell S. Zinner N. T. Phys. Rev. E1002019032107.
  * [11] Andolina G. M., Farina D., Mari A., Pellegrini V., Giovannetti V. Polini M. Phys. Rev. B982018205423.
  * [12] Andolina G. M., Keck M., Mari A., Giovannetti V. Polini M. Phys. Rev. B992019205437.
  * [13] Long W., Sun Q.-F., Guo H. Wang J. Appl. Phys. Lett.8320031397.
  * [14] Santos A. C., Saguia A. Sarandy M. S. Phys. Rev. E1012020062114.
  * [15] Demirplak M. Rice S. A. J. Phys. Chem. A10720039937.
  * [16] Demirplak M. Rice S. A. J. Phys. Chem. B10920056838.
  * [17] Berry M. J. Phys. A: Math. Theor.422009365303.
  * [18] Santos A. C. Sarandy M. S. Sci. Rep.5201515775.
  * [19] Santos A. C., Silva R. D. Sarandy M. S. Phys. Rev. A932016012311.
  * [20] Coulamy I. B., Santos A. C., Hen I. Sarandy M. S. Frontiers in ICT3201619.
  * [21] Hegade N. N., Paul K., Ding Y., Sanz M., Albarrán-Arriagada F., Solano E. Chen X. Phys. Rev. Applied152021024038.
  * [22] Beau M., Jaramillo J. del Campo A. Entropy182016168.
  * [23] Herrera M., Sarandy M. S., Duzzioni E. I. Serra R. M. Phys. Rev. A892014022323.
  * [24] Chen Z., Chen Y., Xia Y., Song J. Huang B. Sci. Rep.6201622202.
  * [25] Chen X., Lizuain I., Ruschhaupt A., Guéry-Odelin D. Muga J. G. Phys. Rev. Lett.1052010123003.
  * [26] Guéry-Odelin D., Ruschhaupt A., Kiely A., Torrontegui E., Martínez-Garaot S. Muga J. G. Rev. Mod. Phys.912019045001.
  * [27] Allahverdyan A. E., Balian R. Nieuwenhuizen T. M. Europhys. Lett.672004565.
  * [28] Çakmak B. Phys. Rev. E1022020042111.
  * [29] Kamian F. H., Tabesh F. T., Salimi S., Kheirandish F. Santos A. C. New J. Phys.222020083007.
  * [30] Quach J. Q. Munro W. J. Phys. Rev. Applied142020024092.
  * [31] Bocchieri P. Loinger A. Phys. Rev.1071957337.
  * [32] Born M. Fock V. Zeitschrift für Physik A Hadrons and Nuclei511928165.
  * [33] Kato T. Journal of the Physical Society of Japan51950435.
  * [34] Messiah A. Quantum Mechanics Quantum Mechanics (North-Holland Publishing Company) 1962.
  * [35] Sarandy M. S., Wu L.-A. Lidar D. A. Quantum Information Processing32004331.
  * [36] Bai S.-Y. An J.-H. Phys. Rev. A1022020060201.
  * [37] Ibáñez S., Chen X., Torrontegui E., Muga J. G. Ruschhaupt A. Phys. Rev. Lett.1092012100403.
  * [38] Torrontegui E., Ibáñez S., Martínez-Garaot S., Modugno M., del Campo A., Guéry-Odelin D., Ruschhaupt A., Chen X. Muga J. G. Chapter 2 \- shortcuts to adiabaticity in Advances in Atomic, Molecular, and Optical Physics, edited by Arimondo E., Berman P. R. Lin C. C. Vol. 62 of _Advances In Atomic, Molecular, and Optical Physics_ (Academic Press) 2013 pp. 117 – 169.
  * [39] Santos A. C., Nicotina A., Souza A. M., Sarthour R. S., Oliveira I. S. Sarandy M. S. EPL (Europhysics Letters)129202030008.
  * [40] Hu C.-K., Cui J.-M., Santos A. C., Huang Y.-F., Sarandy M. S., Li C.-F. Guo G.-C. Opt. Lett.4320183136.
  * [41] Deffner S. arXiv e-prints2021arXiv:2102.05118.
  * [42] Santos A. C. Sarandy M. S. J. Phys. A: Math. Theor.512018025301.
  * [43] Hu H., Qi S. Jing J. arXiv e-prints2021arXiv:2104.12143.


# Effects of Pauli blocking and in-medium nucleon-nucleon cross sections on
the stopping power at low-intermediate energy heavy ion collisions

Xiang Chen  China Institute of Atomic Energy, Beijing 102413, China  Yingxun
Zhang<EMAIL_ADDRESS>China Institute of Atomic Energy, Beijing 102413, China
Zhuxia Li  China Institute of Atomic Energy, Beijing 102413, China

###### Abstract

Three typical algorithms of Pauli blocking in the quantum molecular dynamics
type models are investigated in the nuclear matter, the nucleus and the heavy
ion collisions. The calculations in nuclear matter show that the blocking
ratios obtained with the three algorithms are underestimated 13-25% compared
to the analytical values of blocking ratios. For the finite nucleus, the
spurious collisions occur around the surface of the nucleus owing to the
defects of Pauli blocking algorithms. In the simulations of heavy ion
collisions, the uncertainty of stopping power from different Pauli blocking
algorithms is less than 5%. Furthermore, the in-medium effects of nucleon-
nucleon ($NN$) cross sections on the nuclear stopping power are discussed. Our
results show that the transport models calculations with free $NN$ cross
sections result in the stopping power decreasing with the beam energy at the
beam energy less than 300 MeV/u. To increase or decrease the values of
stopping power, an enhanced or suppressed model dependent in-medium $NN$ cross
section is required.

###### pacs:

Valid PACS appear here

##  I Introduction

Heavy ion collisions (HICs) provide crucial insights into the features of the
nuclear equation of state (EOS) and the in-medium $NN$ cross sections for a
wide range of densities, temperatures, and neutron-proton asymmetries.
However, the transient state of compressed/expanded nuclear matter during the
reaction, such as the pressure and density, can not be directly measured due
to spatial-temporal scale of the reaction system which is beyond the
capability of measurement. To extract the EOS or in-medium $NN$ cross
sections, the transport models which are used to simulate the heavy ion
collisions are indispensable.

Many transport codes have been developed to extract the isospin asymmetric
nuclear EOS and in-medium $NN$ cross sections Aichlin1991PR ; AOno1992POTP ;
SABass1998PPNP ; MPapa2005JOCP ; BALi2008PR ; Zhang2020FOP ; Aichlin2020PRC ;
Bertsch1988PR ; JSu2011PRC ; ZQFeng2011PRC ; Cozma2013PRC ; Nara1999PRC ;
Nara2016PRC ; Niita1995PRC ; Ogawa2015PRC ; QFLi2005PRC ; CMKo1988PRC ;
Colonna2013PRL ; Danie2000NPA . However, the model dependence of the
constraints of symmetry energy Tsang04 ; LWChen05 ; Tsang09 ; Rizzo08 ; Xiao09
; ZQFeng10 ; Xie13 ; JHong14 ; TSong15 ; Cozma17 ; YYLiu21  and in-medium $NN$
cross sections Zhang2007PRC ; QFLi2010MPLA ; QFLi2011PRC ; YJWang2014PRC ;
Henri2020PRC ; Basrak2016PRC  become apparent since the different conclusions
could be drawn from the same data. The situations led to the idea of a
systematic comparison and evaluation of transport codes under controlled
conditionsJXu2016PRC ; Zhang2018PRC ; AOno2019PRC ; Maria21 . Previous studies
along this direction were dedicated to the comparison of transport model
predictions for Au+Au collisionsJXu2016PRC , and for benchmarking the
treatment of nucleon-nucleon collisions and Pauli blockingZhang2018PRC  and
$\Delta$ productionAOno2019PRC  in box calculations with cascade mode. The
comparisons of the mean field have been done in box calculation with Vlasov
modeMaria21 . Currently, the observed differences in the reaction path and
corresponding observables are mainly resulted from differences in the
initialization of the systems and in the treatment of Pauli blocking effects.
The latter, i.e., Pauli blocking, describes the statistical ability to
populate the final states in the fermionic system in the gain(loss) term of
transport equation for a fermionic system and it is crucial for simulating the
low-intermediate energy heavy ion collisions in transport models.

The algorithms of Pauli blocking in Boltzmann-Uehling-Uhlenbeck (BUU) approach
and quantum molecular dynamics (QMD) approach are different. In the BUU
approach, the Pauli blocking may be improved by increasing the number of test
particles to infinity. In the QMD approach, a fixed width of Gaussian wave
packet is used to represent a nucleon. It leads to a strong fluctuation, which
is important in physics for describing the cluster formation and
multifragmentation, but underestimates Pauli blocking effect as in Ref.
Zhang2018PRC . Consequently, one can expect that the successful collision rate
of $NN$ collision could be overestimated in the transport codes, and thus, the
extracted model dependent in-medium $NN$ cross sections deviate from its real
values. Another, owing to the difficulties on the descriptions of Pauli
blocking in QMD approach, different Pauli blocking algorithms have been
developed in QMD codes Zhang2018PRC .

The goal of this work is to learn the systematic deviation of Pauli blocking
ratio compared to the analytical values and evaluate the uncertainties of
nuclear stopping power caused by different Pauli blocking algorithms in the
simulations of HICs. All the calculations are performed within the framework
of the improved quantum molecular dynamics model (ImQMD) Zhang2014PLB ;
Zhang2015PLB ; Zhang2020PRC ; Zhang2020FOP , but the Pauli blocking algorithms
are replaced by the three typical algorithms. The paper is organized as
follows: In Sec. II, we briefly describe the three typical algorithms of Pauli
blocking in the market and the in-medium $NN$ cross sections we try to
analyzed. In Sec. III, the influences of different Pauli blocking algorithms
in the nuclear matter, the finite nucleus and the heavy ion collisions are
presented and discussed. Furthermore, a simple discussion on the influence of
in-medium $NN$ cross sections on the stopping power is also presented. We do
not compare the model calculations to the data for extracting the real values
of in-medium $NN$ cross sections in this work, because the algorithms of Pauli
blocking still need to be improved. A summary and outlook are given in Sec.
IV.

##  II Pauli Blocking and in-medium $NN$ cross sections in ImQMD model

In this section, we only mention the three typical algorithms of Pauli
blocking that are directly related to Uehling-Uhlenbeck factor in the
transport equation, and the in-medium $NN$ cross sections that we used. There
are also some efforts to improve the Pauli blocking in the transport codes by
adding additional constraint Aichelin1985PRC ; Cozma2013PRC ; Cozma2018EPJA ;
QFLi2011PRC , but we will not touch it in this work. More details about the
mean field potential and the treatment of nucleon-nucleon collision in the
ImQMD model can be found in Refs. Zhang2020FOP ; Zhang2014PLB ; Zhang2015PLB ;
Zhang2020PRC .

###  II.1 Algorithms of Pauli blocking

In the quantum molecular dynamics type models, each nucleon is represented by
a Gaussian wave packet,

$\psi_{i}(\mathbf{r})=\frac{1}{(2\pi\sigma_{r}^{2})^{3/4}}{\rm
e}^{-\frac{(\mathbf{r}-\mathbf{r}_{i})^{2}}{4\sigma_{r}^{2}}+{\rm
i}\mathbf{r}\cdot\mathbf{p}_{i}/\hbar},i=1,...,A\vspace*{0.5mm}$ (1)

here, $\sigma_{r}$ and $\mathbf{r}_{i}$ are the width and centroid of wave
packet, respectively. Its Wigner density reads

$f_{W,i}(\mathbf{r},\mathbf{p})=\frac{1}{(\pi\hbar)^{3}}e^{-\frac{(\mathbf{r}-\mathbf{r}_{i})^{2}}{2\sigma_{r}^{2}}-\frac{(\mathbf{p}-{\mathbf{p}_{i}})^{2}}{2\sigma_{p}^{2}}},$
(2)

where $\sigma_{r}\sigma_{p}=\hbar/2$. The Wigner density expresses the
probability density of the simultaneous values of $\mathbf{r}$ and
$\mathbf{p}$ for $i$th nucleon. When Eq. (2) is integrated with respect to
$\mathbf{p}$, the correct probability in coordinate space
$|\psi_{i}(\bm{r})|^{2}$ is given; if we integrate Eq. (2) with respect to
$\mathbf{r}$, the correct probability in momentum space
$|C_{i}(\mathbf{p})|^{2}$ can also be verified Zhang2020FOP .

In the treatment of collision of $i+j\to i^{\prime}+j^{\prime}$ at certain
time step in the code, the positions of particle $i$ and $j$ are kept same
before and after collisions, i.e., $\mathbf{r}_{i}=\mathbf{r}_{i^{\prime}}$,
$\mathbf{r}_{j}=\mathbf{r}_{j^{\prime}}$, while the momenta of particle $i$
and $j$ are changed, i.e., $p_{i}+p_{j}\to p^{\prime}_{i}+p^{\prime}_{j}$.
Thus, the probability of the final state $p^{\prime}_{i}$ occupied by other
nucleons, i.e., $P(p^{\prime}_{i})$, can be calculated based on $\sum_{j\neq
i}f_{W,j}(\mathbf{r},\mathbf{p})$ in phase space cell around $p^{\prime}_{i}$.
In the following discussions, we briefly mention the methods on the
calculation of $P(p^{\prime}_{i})$, such as PB-Wigner, PB-Husimi, PB-HSP,
which are adopted in the different QMD codes.

####  II.1.1 PB-Wigner

For the algorithm of PB-Wigner, the probability of final state
$p^{\prime}_{i}$ occupied by other particles is expressed as
$P_{\tau}(p^{\prime}_{i})=P_{\tau}(\mathbf{r}_{i},\mathbf{p}^{\prime}_{i})$,

$\displaystyle
P_{\tau}(\mathbf{r}_{i},\mathbf{p}^{\prime}_{i})=\frac{1}{2/h^{3}}\sum_{j\in\tau,j\neq
i}f_{W,j}(\mathbf{r}_{i},\mathbf{p}^{\prime}_{i})$ $\displaystyle=$
$\displaystyle 4\sum_{j\in\tau(j\neq
i)}^{A}\exp\left[-\frac{(\mathbf{r}_{i}-\mathbf{r}_{j})^{2}}{2\sigma_{r}^{2}}\right]\times\exp\left[-\frac{(\mathbf{p}^{\prime}_{i}-\mathbf{p}_{j})^{2}}{2\sigma_{p}^{2}}\right],$

with $\tau=n$ or $p$. The factor $2/h^{3}$ results from consideration of the
spin in the phase-space cell.

The $P_{\tau}(p^{\prime}_{i})$ could be larger than 1, because of the
fluctuation and the semi-classical transport equation. If the occupation
probability $P_{\tau}(\mathbf{r}_{i},\mathbf{p}^{\prime}_{i})$ is larger than
$1$, occupation probability $P_{\tau}(\mathbf{r}_{i},\mathbf{p}^{\prime}_{i})$
is replaced by $\min{(P_{\tau}(\mathbf{r}_{i},\mathbf{p}^{\prime}_{i}),1)}$.
This method are used in ImQMD Zhang2020FOP ; Zhang2014PLB ; Zhang2015PLB ;
Zhang2020PRC , IQMD-BNU JSu2011PRC , JAMNara1999PRC ; Nara2016PRC ,
JQMDNiita1995PRC ; Ogawa2015PRC , UrQMD QFLi2005PRC ; QFLi2011PRC ;
SABass1998PPNP . In the ImQMD and UrQMD model, additional criteria is also
adopted to enhance the Pauli blocking ratio Zhang2014PLB ; QFLi2011PRC , but
for convenience, the effect will not be discussed in this paper.

####  II.1.2 PB-Husimi

In the algorithm of PB-Husimi, the probability of final state $p^{\prime}_{i}$
occupied by other particles is expressed according to the Husimi function
Ikeno2020PRC . The Husimi function is ensured to have a good property as a
probability and extra smearing, thus, the distribution is broader than the
Wigner function. In detail, the Husimi phase-space distribution function,
i.e., $f_{H}(\mathbf{r},\mathbf{p})$, in the QMD is related to the Wigner
phase-space distribution $f_{W}(\mathbf{r},\mathbf{p})$ as in Ref.
RFOCon1984PRA ; ZXLi1991CPL ,

$f_{H,i}(\mathbf{r},\mathbf{p})=\int
W(\mathbf{r},\mathbf{p}\,|\mathbf{r}^{\prime},\mathbf{p}^{\prime})f_{W,i}(\mathbf{r}^{\prime},\mathbf{p}^{\prime})\mathrm{d}\mathbf{r}^{\prime}\mathrm{d}\mathbf{p}^{\prime},$
(4)

with

$W(\mathbf{r},\mathbf{p}\,|\mathbf{r}^{\prime},\mathbf{p}^{\prime})=\frac{1}{(\pi\hbar)^{3}}\exp\left[-\frac{(\mathbf{r}^{\prime}-\mathbf{r})^{2}}{2\sigma_{r}^{2}}-\frac{(\mathbf{p}^{\prime}-\mathbf{p})^{2}}{2\sigma_{p}^{2}}\right].$
(5)

By using Eq. (4) and Eq. (5), we obtain the Husimi phase-space distribution
function as follows,

$f_{H,i}(\mathbf{r},\mathbf{p})=\frac{1}{h^{3}}\exp\left[-\frac{(\mathbf{r}-\mathbf{r}_{i})^{2}}{4\sigma_{r}^{2}}-\frac{(\mathbf{p}-\mathbf{p}_{i})^{2}}{4\sigma_{p}^{2}}\right].$
(6)

Thus, the occupation probability at $\mathbf{p}=\mathbf{p}^{\prime}_{i}$ and
$\mathbf{r}=\mathbf{r}_{i}$ can be obtained as
$P_{\tau}(p^{\prime}_{i})=P_{\tau}(\mathbf{r}_{i},\mathbf{p}^{\prime}_{i})$,

$\displaystyle
P_{\tau}(\mathbf{r}_{i},\mathbf{p}^{\prime}_{i})=\frac{1}{2/h^{3}}\sum_{j\in\tau,j\neq
i}f_{H,j}(\mathbf{r}_{i},\mathbf{p}^{\prime}_{i})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{j\in\tau(j\neq
i)}^{A}\exp\left[-\frac{(\mathbf{r}_{i}-\mathbf{r}_{j})^{2}}{4\sigma_{r}^{2}}\right]\times\exp\left[-\frac{(\mathbf{p}^{\prime}_{i}-\mathbf{p}_{j})^{2}}{4\sigma_{p}^{2}}\right].$

Similarly, if the occupation probability
$P_{\tau}(\mathbf{r}_{i},\mathbf{p}^{\prime}_{i})$ is larger than $1$,
$P_{\tau}(\mathbf{r}_{i},\mathbf{p}^{\prime}_{i})=\min{(P_{\tau}(\mathbf{r}_{i},\mathbf{p}^{\prime}_{i}),1)}$.

####  II.1.3 PB-HSP

In the algorithm of PB-HSP, the occupation probability
$P_{\tau}(\mathbf{r}_{i},\mathbf{p}^{\prime}_{i})$ is calculated as
$P_{\tau}(p^{\prime}_{i})=P_{\tau}(\mathbf{r}_{i},\mathbf{p}^{\prime}_{i})$,

$P_{\tau}(\mathbf{r}_{i},\mathbf{p}^{\prime}_{i})=\sum_{j\in\tau(j\neq
i)}^{A}(O_{ij}^{(x)}/\frac{4}{3}\pi R_{x}^{3})(O_{ij}^{(p)}/\frac{4}{3}\pi
R_{p}^{3}),$ (8)

where, $O_{ij}^{(x)}(O_{ij}^{(p)})$ is the volume of the overlap region of
spheres with the radius $R_{x}(R_{p})$ of nucleons $i$ and $j$ in coordinate
(momentum) space. This method has been adopted in QMD Aichelin1985PRC , LQMD
ZQFeng2011PRC  and TuQMD (the additional constraints named as surface
correction is also adopted in QMD and TuQMD) Cozma2013PRC ; Cozma2018EPJA .
The $O_{ij}^{(x)}$ is calculated by

$O_{ij}^{(x)}=\left\\{\begin{array}[]{ll}0\,,&L_{ij}^{(x)}\geq 2R_{x}\\\
\frac{4}{3}\pi R_{x}^{3}-\pi
L_{ij}^{(x)}\left[R_{x}^{2}-\frac{1}{3}(\frac{L_{ij}^{(x)}}{2})^{2}\right],&L_{ij}^{(x)}<2R_{x}\end{array}\right.$
(9)

where $L_{ij}^{(x)}=|\mathbf{r}_{i}-\mathbf{r}_{j}|$. By replacing
$L_{ij}^{(x)}$ to $L_{ij}^{(p)}=|\mathbf{p}^{\prime}_{i}-\mathbf{p}_{j}|$ and
$R_{x}$ to $R_{p}$, one can obtain $O_{ik}^{(p)}$. Here, we take $R_{x}$=3.367
fm , $R_{p}$=112.5 MeV/$c$ as in Ref. Zhang2018PRC . This method corresponds
to the uniform phase-space density distribution in coordinate and momentum
space, i.e., one nucleon occupies the phase space cell with the size
$\frac{4}{3}\pi R_{x}^{3}\cdot\frac{4}{3}\pi R_{p}^{3}=\frac{h^{3}}{2}$.
Similarly, if the occupation probability
$P_{\tau}(\mathbf{r}_{i},\mathbf{p}^{\prime}_{i})$ is larger than $1$,
occupation probability should be replaced by
$\min(P_{\tau}(\mathbf{r}_{i},\mathbf{p}^{\prime}_{i}),1)$.

###  II.2 In-medium $NN$ cross sections

In the ImQMD, the isospin dependent $NN$ cross sections and the differential
cross sections in free space are taken from Ref. Cugnon1996NIM . The in-medium
$NN$ cross sections in the ImQMD model is named as $\sigma^{med}_{QMD}$, and
$\sigma^{med}_{QMD}=R*\sigma^{free}$, where the medium correction factor
$R=(1+\eta(E_{beam})\rho/\rho_{0})$. Here, we want to stress the model
dependent $\sigma^{med}_{QMD}$ is not exactly same as in true $\sigma^{med}$,
and we will show it in the nuclear matter calculations in Sec. III.1.

##  III Results and discussion

In this section, we will first evaluate the algorithms of Pauli blocking in
nuclear matter in cascade mode. Then, the effects of Pauli blocking on the
finite nucleus and on heavy ion collision observable, i.e., stopping power,
are presented and discussed. The uncertainties of different Pauli blocking
algorithms on the stopping power are obtained.

###  III.1 Pauli blocking in nuclear matter

During the heavy ion collisions, various effects interplay and propagate
during the reaction process, thus, it is hard to evaluate the Pauli blocking
only in the heavy ion collisions. To disentangle the interplay between mean
field and $NN$ collisions, we first analyze the nucleon number distribution of
final state $\mathbf{p}^{\prime}_{i}$ (i.e., $dN/dp^{\prime}$) and occupation
probability of final state $\mathbf{p}^{\prime}_{i}$, i.e.,
$P_{\tau}(p^{\prime}_{i})$. The calculations are performed in the cascade
mode, i.e., without mean-field potential, in the nuclear matter.

To simulate the nuclear matter, a box with imposed periodic boundary
conditions is adopted. The periodic boundary conditions are considered as same
as in Ref. Zhang2018PRC , i.e., the dimensions of the cubic box are
$L_{\alpha}$ = 20 fm, $\alpha\equiv x,y,z$, and the position of the center of
box is ($L_{x}/2,L_{y}/2,L_{z}/2$). The box is initialized with a finite
temperature in uniform nuclear matter. The density is $\rho=0.16$ fm-3 and
isospin asymmetry equals zero, which corresponds to 1280 nucleons (640
neutrons + 640 protons) in the box. In coordinate space, the positions of
nucleon are initialized randomly from $0$ to $L_{\alpha}$. In momentum space,
the momenta of nucleon are initialized according to the Fermi-Dirac
distribution, $f=1/\\{1+$exp$[(\epsilon-\mu)/T]\\}$, with $\epsilon=p^{2}/2m$,
nucleon mass $m$ = 939 MeV/$c^{2}$, chemical potential $\mu$ = 36.84 MeV, and
temperature $T$ = 5 MeV. The $NN$ collision cross sections are set as 40 mb.

Figure 1: (Color online) Panels (a), (b), and (c) are the momentum
distribution of final state of nucleon-nucleon collisions at t=1 fm/$c$.
Panels (d), (e), and (f) are occupation probability for different Pauli
blocking algorithms, i.e., PB-Wigner, PB-Husimi, PB-HSP. Red lines are for the
analytical values of occupation probability (see text for more details).

In the upper panels of Fig. 1, we present the $dN/dp^{\prime}$ after the
evolution of one time step, i.e., at $t=1$ fm/$c$, to understand how many
nucleons participate in the $NN$ collisions at certain momentum. Panels (a),
(b), and (c) are for the cases of PB-Wigner, PB-Husimi, and PB-HSP,
respectively. 10,000 events are performed at $T$ = 5 MeV. The $dN/dp^{\prime}$
of PB-Wigner, PB-Husimi, and PB-HSP are almost same, because the momentums of
nucleons in initial state are sampled within the same Fermi-Dirac distribution
and only one time step is considered. $dN/dp^{\prime}$ increases with momentum
and reaches maximum around 240 MeV/$c$ and then decreases.

Fig. 1 (d), (e), and (f) show the occupation probability for the Pauli
blocking algorithms of PB-Wigner, PB-Husimi, PB-HSP, respectively. The mean
values of the occupation probability, i.e., $\langle P(p^{\prime})\rangle$,
and their standard deviation, i.e. $\langle(P(p^{\prime})-\langle
P(p^{\prime})\rangle)^{2}\rangle^{1/2}$, are shown as the blue curves and the
blue error bars. The actual averaged occupation probabilities used in the
ImQMD calculations, i.e., $\langle$min$(P(p^{\prime}),1)\rangle$, are shown as
the black curves. The actual occupation probability
$\langle$min$(P(p^{\prime}),1)\rangle$ is always slightly lower than $\langle
P(p^{\prime})\rangle$ due to the truncation of $P(p^{\prime})$ by using
min$(P(p^{\prime}),1)$. Red lines are the analytical values of occupation
probability, i.e., the Fermi-Dirac occupation probability. Among these three
Pauli blocking algorithms, the PB-Husimi has a smallest standard deviations
due to the large width in Eq. (II.1.2). Furthermore, all those algorithms used
in the QMD codes deviate from the theoretical values as discussed in Ref.
Zhang2018PRC . It underestimates the Pauli blocking probability at lower
momentum region and overestimates the blocking probability at high momentum
region due to the strong fluctuation.

To quantitatively evaluate the algorithms of Pauli blocking, we present the
attempted collision rate $\frac{dN^{att.}_{coll}}{dt}$, successful collision
rate $\frac{dN^{suc.}_{coll.}}{dt}$, and Pauli blocking ratio $R_{block}$
defined as,

$R_{block}=1-\frac{dN^{suc.}_{coll.}}{dt}/\frac{dN^{att.}_{coll.}}{dt},$ (10)

as a function of temperature in Fig. 2 (a) and (b), respectively. The
temperature $T$ is from 2 MeV to 10 MeV, which corresponds to the low-
intermediate energy heavy ion collisions. The black line is the result for
attempted collisions rate which is calculated as in Ref. Zhang2018PRC , i.e.,

$\langle\frac{dN^{att.}_{coll.}}{dt}\rangle=\frac{1}{2}A\rho\langle
v_{rel}\sigma^{med}\rangle.$ (11)

Here, $A$ is the nucleon number, $\rho$ is the density, $v_{rel}$ is the
relative velocity between two colliding nucleons, $\sigma^{med}$ is the in-
medium $NN$ cross section and $\sigma^{med}=40$ mb. Eq.(11) corresponds to the
case of no Pauli blocking. The bracket $\langle\cdot\rangle$ means the
averaged values over the time interval 60-140 fm/$c$ and from 200 events. As
illustrated in Fig. 2 (a), the averaged attempted collision rate (black line)
increases with the temperature increasing.

Figure 2: (Color online) Panel (a): Attempted collision rate
$\langle\frac{dN^{att.}_{coll.}}{dt}\rangle$ and successful collision rate
$\langle\frac{dN^{suc.}_{coll.}}{dt}\rangle$ as functions of temperature. The
red lines with different symbols are from different Pauli blocking algorithms.
Panel (b): The averaged blocking ratio $\langle R_{block}\rangle$ for
different Pauli blocking (see text for more details).

The red lines with different symbols in Fig. 2 (a) are the averaged successful
collision rates obtained with PB-Wigner (square), PB-Husimi (circle), PB-HSP
(triangle). The green line represents the result of the Pauli blocker which is
fixed to the initialized Fermi-Dirac distribution for given temperatures,
i.e., PB-FD. It is used for evaluating how well the Pauli blocking algorithms
is in the nuclear matter. Similar as the results in Ref. Zhang2018PRC , the
averaged successful collision rates from PB-Wigner, PB-Husimi and PB-HSP are
larger than the analytical successful collision rates from PB-FD. Among them,
the results from PB-HSP are more closer to the PB-FD than others.

In Fig. 2 (b), we plot the averaged blocking ratio $\langle R_{block}\rangle$
as a function of $T$. The values of $\langle R_{block}\rangle$ from different
Pauli blocking algorithms decrease with temperature increasing. However, the
$\langle R_{block}\rangle$ obtained with PB-Wigner, PB-Husimi, and PB-HSP are
smaller about 13-25% than the analytical values, i.e., $\langle
R_{block}\rangle$ obtained with PB-FD. It means that the different Pauli
blocking algorithms used in the transport codes overestimate the successful
collision rate. Thus, to get the same successful collision rate from the
algorithms of Pauli blocking adopted in QMD and from the analytical Pauli
blocking, one can image that the $\sigma^{med}_{QMD}$ is smaller than its true
values, i.e., $\sigma^{med}_{QMD}<\sigma^{med}$. For example, the values of
$\xi=\sigma^{med}_{QMD}/\sigma^{med}\approx 37-50\%$ for the selected three
Pauli blocking algorithms at $T=10$ MeV.

###  III.2 Pauli blocking in a finite nucleus

One should notice that the Pauli blocking effect in the nuclear matter is not
exactly the same as that in the finite size system due to the boundary effect.
Thus, we further check it in a finite nucleus where the mean-field potential
is also considered for binding the nucleons together. In the following
calculations, the interaction parameter set of SkM* is adopted.

Fig. 3 (a) shows the number of successful $NN$ collision as a function of
radial distance $r$, where the collisions occurs, in ${}^{124}\mathrm{Sn}$.
The lines with different colors represent the results obtained with the PB-
Wigner (black line), PB-Husimi (red line), and PB-HSP (green line),
respectively. Similar to the calculations in nuclear matter, the numbers of
successful collisions are obtained at 1 fm/$c$ and the values of successful
collision rates are about 0.9-1.6 $c$/fm. Ideally, the successful collisions
rate of initial nucleus should be zero, but it happens owing to the defect of
Pauli blocking algorithms in the QMD models. Among the three kinds of Pauli
blocking, the calculation with PB-HSP gives the smallest successful collision
rate.

In addition, most of the successful $NN$ collisions occur around $r$=5.2 fm,
which is close to the surface of nucleus as found in the density distribution
plots in panel (b). There were some ad hoc methods to overcome this defect in
Refs. Aichelin1985PRC ; Cozma2018EPJA , but a consistent improvement of Pauli
blocking near the surface of nucleus or reaction system is still a theoretical
challenge.

Figure 3: (Color online) Panel (a): the successful collision number as a
function of radial distance for 124Sn, (b) density distribution of 124Sn.
Panels (c), (d) and (e), momentum distribution of scattered nucleons; (f),
(g), and (h): occupation probability for different Pauli blocking algorithms,
i.e., PB-Wigner, PB-Husimi, PB-HSP (see the text for details).

The panels (c), (d) and (e) show $dN/dp^{\prime}$ of the final state of $NN$
collisions in different Pauli blocking algorithms, i.e., PB-Wigner, PB-Husimi,
and PB-HSP, respectively, and it clearly illustrates the importance of Pauli
blocking within the Fermi momentum. The mean values of occupation probability
$\langle P(p^{\prime})\rangle$ and mean actual occupation probability
$\langle$min$(P(p^{\prime}),1)\rangle$ at 1 fm/$c$ are shown by the blue and
black curves in panel (f)-(h), respectively. Similar to the finding in the
nuclear matter, there is also $\langle\min(P(p^{\prime}),1)\rangle\leq\langle
P(p^{\prime})\rangle$ in finite nucleus. The standard deviation of the
occupation probability $P(p^{\prime})$ is indicated by blue error bars, and it
shows that the PB-Wigner has the largest fluctuation among the three Pauli
blocking algorithms. The values obtained with PB-HSP are closer to 1 at low
momentum than other two algorithms.

###  III.3 Pauli blocking on stopping power in heavy ion collisions

The nuclear stopping governs the amount of dissipated energy under the
competition between the mean field potential and nucleon-nucleon collisions,
and can be measured by the ratio between transverse and longitudinal
components of kinematical observables Reisdorf2004PRL ; Lehaut2010PRL ;
Lopez2014PRC . For example, the ratio of the variances of the transverse
rapidity distribution to that of the longitudinal rapidity distributions of
emitted particles, which is named as $vartl$ Reisdorf2004PRL ,

$vartl=\frac{\langle y_{t}^{2}\rangle}{\langle y_{z}^{2}\rangle}.$ (12)

The energy-based isotropy ratio $R_{E}$ Lehaut2010PRL  and the momentum-based
isotropy ratio $R_{p}$ Lopez2014PRC  are also used in experiments to measure
the stopping power. Those quantities measure the transfer of momentum from
entrance direction to transverse direction, and thus, they are closely related
to the successful $NN$ collision rate. So, a study of nuclear stopping
observable can provide an insight on the in-medium $NN$ cross sections and
Pauli blocking effects.

Before discussing the stopping power in HICs, we first present the free $NN$
cross sections Cugnon1996NIM  used in the ImQMD calculations and the
corresponding mean attempted collision rate in Fig. 4 (a) and (b). As shown in
panel (a), both the cross sections of nn/pp (black line) and np (red line)
decrease with the beam energy. More details, in the ImQMD model, we also set
$\sigma^{free}_{nn/pp}=60$ mb and $\sigma^{free}_{np}=180$ mb at $p_{lab}<0.3$
GeV/$c$ (or $E_{lab}<50$ MeV) to avoid the spurious low energy $NN$ collisions
in nuclear medium. Fig. 4 (b) shows that the mean attempted collision rate by
nucleons, i.e., $\langle\rho v_{rel}\sigma^{free}_{NN}\rangle$, for nn/pp and
np, in the uniform nuclear matter with $\rho=0.16$ fm-3. The mean attempted
collision rate for np decreases with the beam energy increasing, while the
mean attempted collision rate for nn/pp first decreases and then increases at
$E_{lab}>\sim 100$ MeV. It is different with the behaviors in Fig. 2 where a
constant $NN$ cross section $\sigma=40$ mb is used. Thus, the stopping power
may weakly depend on the beam energy or decrease with the beam energy
increasing if $\sigma^{free}_{nn/pp,np}$ is used and Pauli blocking is
switched off in the model.

Figure 4: (Color online). Panel (a), the cross section of nn/pp and np used in
ImQMD model; panel (b), the mean attempted collision rate obtained with the
free $NN$ cross sections in nuclear matter.

However, the Pauli blocking effect is indispensable for low-intermediate
energy heavy ion collisions. In Fig. 5 (a), the average blocking ratios
$\langle R_{block}\rangle$ over the time interval 0-400 fm/$c$ as a function
of beam energy for ${}^{112}Sn+^{124}Sn$ at b=1 fm are presented. The black
lines with different symbols are the results obtained with PB-Wigner
(squares), PB-Husimi (circles), and PB-HSP (triangles), in the case of
$\sigma^{free}$ is adopted in the ImQMD calculations. The values of $\langle
R_{block}\rangle$ obtained with PB-Husimi and PB-HSP are almost same, but the
PB-Wigner results in the largest $\langle R_{block}\rangle$ among three kinds
of Pauli blocking algorithms. It is opposite to the finding in the nuclear
matter, and can be understood from the reaction dynamics. At early stage of
reaction, the weakest Pauli blocking algorithm observed in the nuclear matter
calculations, i.e., PB-Wigner, results in more $NN$ collisions than that for
PB-Husim or PB-HSP. More $NN$ collisions provide a larger repulsion for
nucleons during the compressed stage, and make the system expand to larger
momentum space than that with less $NN$ collisions. Thus, the successful
collision rates for PB-Wigner become smaller than that for PB-Husimi or PB-HSP
after the compression stage. Thus, the largest $\langle
R_{block}\rangle=\langle
1-\frac{dN^{suc.}_{coll.}}{dt}/\frac{dN^{att.}_{coll.}}{dt}\rangle$ values for
PB-Wigner are obtained by averaging over the time interval 0-400 fm/$c$. This
effect becomes obvious at high beam energy where the $NN$ collisions are more
frequent than low beam energies.

Figure 5: (Color online). The beam energy dependence of $\langle
R_{block}\rangle$ (a), of the mean attempted and successful collision rates
(b), and of $vartl$ (c). The different symbols are for different Pauli
blocking, and different colors are for different in-medium cross sections.

To explore the influence of different algorithms of Pauli blocking on the
stopping power, we present the averaged attempted and successful collision
rates over time interval 0-400 fm/$c$ as a function of beam energy in Fig. 5
(b) and $vartl$ in panel (c). The values of successful collision rates are in
the range of 2-5 $c$/fm, and the different Pauli blocking algorithms lead to
about 11-15% difference on the successful collision rate. The averaged
attempted collision rates and the averaged successful collision rates decrease
with the beam energy increasing.

The stopping power $vartl$ is calculated from the rapidity distribution of
charged particles $Z=1-6$ weighted by their charge number to weaken the defect
on the cluster formation mechanism at low-intermediate energy heavy ion
collisions. The values of $vartl$ weakly depend on the algorithms of Pauli
blocking we used owing to the small number of successful collision rate, as
shown in panel (b). Furthermore, the calculations with $\sigma^{free}$ predict
that the values of $vartl$ decrease with the beam energy increasing. This
behavior is opposite to the observation in experimentsReisdorf2004PRL  where
the $vartl$ increases with the beam energy increasing for Au+Au at
$E_{beam}>100$ MeV/u. This difference illustrates that the in-medium
correction on the $NN$ cross sections is needed.

###  III.4 In-medium $NN$ cross sections on stopping power

To test the effect of in-medium $NN$ cross sections on stopping power, we
simply take the $\sigma^{med}_{QMD}$ as
$\sigma^{med}_{QMD}=(1+\eta(E_{beam})\rho/\rho_{0})\sigma^{free}$ in the code.
For $E_{beam}\leq 100$ MeV, $\eta$ is set as -0.2. At $E_{beam}=150,200,300$
MeV, $\eta=0.0,0.2,0.8$, respectively. The red lines with symbols in panels
(a), (b) and (c) are the results obtained with $\sigma^{med}_{QMD}$. As
illustrated in Fig. 5 (a), the values of $\langle R_{block}\rangle$ weakly
depend on the correction of in-medium $NN$ cross sections. But, as shown in
panel (b), the in-medium $NN$ cross sections reduce the attempted collision
and successful collision rates by $\sim$ 7% at $E_{beam}\leq 100$ MeV and
enhance the attempted and successful collision rates by $>7$% at
$E_{beam}>100$ MeV. Consequently, as in panel (c), the values of $vartl$ are
influenced and they are suppressed at $E_{beam}\leq 100$ MeV/u, and enhanced
at $E_{beam}>100$ MeV/u.

However, one should bear in mind that a reliable extraction of in-medium $NN$
cross section by comparing HIC data to transport model calculations requires
an accuracy method to treat the Pauli blocking in the simulation of HICs. The
improvement of Pauli blocking is still on the way.

##  IV Summary

In summary, we first evaluate the different Pauli blocking algorithms in the
nuclear matter in cascade mode, i.e., only with $NN$ collisions. Our
calculations show that averaged occupation probability obtained with PB-Husimi
and PB-HSP are closer to the analytical values than that with PB-Wigner, but
all the three algorithms used in the QMD codes underestimate the Pauli
blocking ratio by 13-25% at $T\leq 10$ MeV. This underestimation may leads to
the extracted in-medium $NN$ cross sections from QMD type models are 37-50%
smaller than the true in-medium $NN$ cross sections, and $\sigma^{med}=2\sim
2.7\sigma^{med}_{QMD}$.

Furthermore, we analyzed the different Pauli blocking algorithms in the finite
nucleus and heavy ion collisions, in which both mean field potential and $NN$
collisions are included. For the finite nucleus, the Pauli blocking ratios are
in the range of 69-83% for different Pauli blocking algorithms owing to the
defects of Pauli blocking in QMD model. The spurious $NN$ collision mainly
occurs around the surface of finite nucleus. There have been some effort to
improve the Pauli blocking, especially near the surface of nucleus or the
reaction systems, but a consistent treatment is still a big challenge for
many-body transport theory.

By using the algorithms of Pauli blocking in the current market, the
influences of different Pauli blocking algorithms on the excitation function
of stopping power in HICs are discussed, and the uncertainties of stopping
power with different Pauli blocking are less than 5%. If one would like to
produce the behaviors of the $vartl$ increase with the beam energy, a strong
enhancement of in-medium $NN$ cross sections is needed. For obtaining the true
values of in-medium $NN$ cross sections by comparing the heavy ion collision
data to the transport model calculations, a refined Pauli blocking algorithm
is urged in future.

## Acknowledgements

The authors thank the helpful discussions with Dr. Yongjia Wang, Dr. D. Cozma
on the Pauli blocking algorithms. This work was partly inspired by the
transport code comparison project, and it was supported by the National
Natural Science Foundation of China Nos. 11875323, 11705163, 11790320,
11790323, and 11961141003, the National Key R&D Program of China under Grant
No. 2018YFA0404404, the Continuous Basic Scientific Research Project (No.
WDJC-2019-13, BJ20002501) and the funding of China Institute of Atomic Energy.

## References

  * (1) J. Aichelin, Phys. Rep. 202, 233 (1991).
  * (2) B. A. Li, L. W. Chen, and C. M. Ko, Phys. Rep. 464, 113(2008).
  * (3) Y. X. Zhang, N. Wang, Q. F. Li, L. Ou, J. L. Tian, M. Liu, K. Zhao, X. Z. Wu, and Z. X. Li, Front. Phys. 15(5), 54301 (2020).
  * (4) A. Ono, H. Horiuchi, T. Maruyama, and A. Ohnishi, Prog. Theo. Phys. 87, 1185 (1992).
  * (5) M. Papa, G. Giuliani, and A. Bonasera, J. Comp. Phys. 208, 403 (2005).
  * (6) S. A. Bass, M. Belkacem, M. Bleicher, M. Brandstetter, L. Bravina, C. Ernst, L. Gerland, M. Hofmann, S. Hofmann, S. Hofmann, J. Konopka, G. Mao, L. Neise, S. Soff, C. Spieles, H. Weber, L. A. Winckelmann, H. Stöcker, W. Greiner, C. Hartnack, J. Aichelin, and N. Amelin, Prog. Part. Nucl. Phys. 41, 255 (1998).
  * (7) J. Aichelin, E. Bratkovskaya, A. L. Fèvre, V. Kireyeu, V. Kolesnikov, Y. Leifels, V. Voronyuk, and G. Coci, Phys. Rev. C 101, 044905 (2020).
  * (8) J. Su, F. S. Zhang, and B. A. Bian, Phys. Rev. C 83, 014608 (2011).
  * (9) Z. Q. Feng, Phys. Rev. C 84, 024610 (2011).
  * (10) M. D. Cozma, Y. Leifels, W. Trautmann, Q. Li, and P. Russotto, Phys. Rev. C 88, 044912 (2013).
  * (11) Y. Nara, N. Otuka, A. Ohnishi, K. Niita, and S. Chiba, Phys. Rev. C 61, 024901 (1999).
  * (12) N. Ikeno, A. Ono, Y. Nara, and A. Ohnishi, Phys. Rev. C 93, 044612 (2016).
  * (13) K. Niita, S. Chiba, T. Maruyama, T. Maruyama, H. Takada, T.Fukahori, Y. Nakahara, and A. Iwamoto, Phys. Rev. C 52, 2620 (1995).
  * (14) T. Ogawa, T. Sato, S. Hashimoto, D. Satoh, S. Tsuda, and K. Niita, Phys. Rev. C 92, 024614 (2015).
  * (15) Q. Li, Z. Li, S. Soff, M. Bleicher, and H. Stöcker, Phys. Rev. C 72, 034613 (2005).
  * (16) G. F. Bertsch, and S. D. Gupta, Phys. Rep. 160, 189 (1988).
  * (17) C. M. Ko, and Q. Li, Phys. Rev. C 37, 2270 (1988).
  * (18) M. Colonna, Phys. Rev. Lett. 110, 042701 (2013).
  * (19) P. Danielewicz, Nucl. Phys. A 673, 375 (2000).
  * (20) M. B. Tsang, T. X. Liu, L. Shi, P. Danielewicz, C. K. Gelbke, X. D. Liu, W. G. Lynch, W. P. Tan, G. Verde, A. Wagner, H. S. Xu, W. A. Friedman, L. Beaulieu, B. Davin, R. T. de Souza, Y. Larochelle, T. Lefort, R. Yanez, V. E. Viola Jr., R. J. Charity, and L. G. Sobotka, Phys. Rev. Lett. 92, 062701 (2004).
  * (21) Lie-Wen Chen, Che Ming Ko, and Bao-An Li, Phys. Rev. Lett. 94, 032701 (2005).
  * (22) M. B. Tsang, Y. Zhang, P. Danielewicz, M. Famiano, Z. Li, W. G. Lynch, and A.W. Steiner, Phys. Rev. Lett. 102, 122701(2009).
  * (23) J. Rizzo, M. Colonna, V. Baran, M. Di Toro, H.H.Wolter, M. Zielinska-Pfabe, Nucl. Phys. A 806, 79 (2008).
  * (24) Zhigang Xiao, Bao-An Li, Lie-Wen Chen, Gao-Chan Yong, and Ming Zhang, Phys. Rev. Lett. 102, 062502 (2009).
  * (25) Zhao-Qing Feng, Gen-Ming Jin, Phys. Lett. B 683, 140 (2010).
  * (26) Wen-Jie Xie, Jun Su, Long Zhu, Feng-Shou Zhang, Phys. Lett. B 718, 1510(2013).
  * (27) Jun Hong, P. Danielewicz, Phys. Rev. C 90, 024605 (2014).
  * (28) Taesoo Song and Che Ming Ko, Phys. Rev. C 91, 014901 (2015).
  * (29) M. D. Cozma, Phys. Rev. C 95, 014601 (2017).
  * (30) Yangyang Liu, Yongjia Wang, Ying Cui, Cheng-Jun Xia, Zhuxia Li, Yongjing Chen, Qingfeng Li, and Yingxun Zhang, Phys.Rev.C 103, 014616 (2021).
  * (31) Y. X. Zhang, Z. X. Li, and P. Danielewicz, Phys. Rev. C 75, 034615 (2007).
  * (32) Q. F. Li, C. W. Shen, and M. D. Toro, Mod. Phys. Lett. A 25, 669 (2010).
  * (33) Q. F. Li, C. W. Shen, C. C. Guo, Y. J. Wang, Z. X. Li, J. Lukasik, and W. Trautmann, Phys. Rev. C 83, 044617 (2011).
  * (34) Y. J. Wang, C. C. Guo, Q. F. Li, H. F. Zhang, Z. X. Li, and W. Trautmann, Phys. Rev. C 89, 034606 (2014).
  * (35) M. Henri, O. Lopez, D. Durand, B. Borderie, R.Bougault, A. Chbihi, Q. Fable, J. D. Frankland, E. Galichet, D. Gruyer, J. Lemarié, N. L. Neindre, I. Lombardo, M. Pârlog, J. Quicray, G. Verde, E. Vient, and M. Vigilante, (INDRA Collaboration), Phys. Rev. C 101, 064622 (2020).
  * (36) Z. Basrak, P. Eudes, and V. de la Mota, Phys. Rev. C 93, 054609 (2016).
  * (37) J. Xu, L. W. Chen, M. Y. B. Tsang, H. Wolter, Y. X. Zhang, J. Aichelin, M. Colonna, D. Cozma, P. Danielewicz, Z. Q. Feng, A. L. Fèvre, T. Gaitanos, C. Hartnack, K. Kim, Y. Kim, C. M. Ko, B. A. Li, Q. F. Li, Z. X. Li, P. Napolitani, A. Ono, M. Papa, T. Song, J. Su, J. L. Tian, N. Wang, Y. J. Wang, J. Weil, W. J. Xie, F. S. Zhang, and G. Q. Zhang, Phys. Rev. C 93, 044609 (2016).
  * (38) Y. X. Zhang, Y. J. Wang, M. Colonna, P. Danielewicz, A. Ono, M. B. Tsang, H. Wolter, J. Xu, L. W. Chen, D. Cozma, Z. Q. Feng, S. D. Gupta, N. Ikeno, C. M. Ko, B. A. Li, Q. F. Li, Z. X. Li, S. Mallik, Y. Nara, T. Ogawa, A. Ohnishi, D. Oliinychenko, M. Papa, H. Petersen, J. Su, T. Song, J. Weil, N. Wang, F. S. Zhang, and Z. Zhang, Phys. Rev. C 97, 034625 (2018).
  * (39) A. Ono, J. Xu, M. Colonna, P. Danielewicz, C. M. Ko, M. B. Tsang, Y. J. Wang, H. Wolter, Y. X. Zhang, L. W. Chen, D. Cozma, H. Elfner, Z. Q. Feng, N. Ikeno, B. A. Li, S. Mallik, Y. Nara, T. Ogawa, A. Ohnishi, D. Oliinychenko, J. Su, T. Song, F. S. Zhang, and Z. Zhang, Phys. Rev. C 100, 044617 (2019).
  * (40) Maria Colonna, Yingxun Zhang, Yongjia Wang, Dan Cozma, Pawel Danielewicz, Che Ming Ko, Akira Ono, Manyee Betty Tsang, et al., preparing for submission.
  * (41) Y. X. Zhang, M. B. Tsang, Z. X. Li, and H. Liu, Phys. Lett. B 732, 186 (2014).
  * (42) Y. X. Zhang, M. B. Tsang, Z. X. Li, Phys. Lett. B 749, 262 (2015).
  * (43) Y. X. Zhang, M. Liu, C. J. Xia, Z. X. Li, and S. K. Biswal, Phys. Rev. C 101, 034303 (2020).
  * (44) J. Aichelin, and G. Bertsch, Phys. Rev. C 31, 1730 (1985).
  * (45) M. D. Cozma, Eur. Phys. J. A 54, 40 (2018).
  * (46) N. Ikeno, A. Ono, Y. Nara, and A. Ohnishi, Phys. Rev. C 101, 034607 (2020).
  * (47) R. F. O’Connell, Lipo Wang, and H. A. Williams, Phys. Rev. A 30, 2187 (1984).
  * (48) Z. X. Li, X. Z. Wu, and Y. Z. Zhuo, Chin. Phys. Lett. 8, 384 (1991).
  * (49) J. Cugnon, D. L’Hôte, and J. Vandermeulen, Nucl. Instrum. Meth. B 111, 215 (1996).
  * (50) G. Lehaut, D. Durand, O. Lopez, E. Vient, A. Chbihi, J. D. Frankland, E. Bonnet, B. Borderie, R. Bougault, E. Galichet, D. Guinet, P. Lautesse, N. Le Neindre, P. Napolitani, M. Parlog, M. F. Rivet, and E. Rosato (INDRA Collaboration), Phys. Rev. Lett. 104, 232701 (2010).
  * (51) W. Reisdorf, A. Andronic, A. Gobbi, O. N. Hartmann, N. Herrmann, K. D. Hildenbrand, Y. J. Kim, M. Kirejczyk, P. Koczoń, T. Kress, Y. Leifels, A. Schüttauf, Z. Tymiński, Z. G. Xiao, J. P. Alard, V. Barret, Z. Basrak, N. Bastid, M. L. Benabderrahmane, R. Čaplar, P. Crochet, P. Dupieux, M. Dželalija, Z. Fodor, Y. Grishkin, B. Hong, J. Kecskemeti, M. Korolija, R. Kotte, A. Lebedev, X. Lopez, M. Merschmeyer, J. Mösner, W. Neubert, D. Pelte, M. Petrovici. F. Rami, B. de Schauenburg, Z. Seres, B. Sikora, K. S. Sim, V. Simion, K. S. Wilczyńska, V. Smolyankin, M. Stockmeier, G. Stoicea, P. Wagner, K. Wiśniewski, D. Wohlfarth, I. Yushmanov, and A. Zhilin, (FOPI Collaboration), Phys. Rev. Lett. 92, 232301 (2004).
  * (52) O. Lopez, D. Durand, G. Lehaut, B. Borderie, J. D. Frankland, M. F. Rivet, R. Bougault, A. Chbihi, E. Galichet, D. Guinet, M. L. Commara, N. L. Neindre, I. Lombardo, L. Manduci, P. Marini, P. Napolitani, M. Pârlog, E. Rosato, G. Spadaccini, E. Vient, and M. Vigilante,(INDRA Collaboration), Phys. Rev. C 90, 064602 (2014).