# The Volume of Healthy Red Blood Cells is Optimal for Advective Oxygen
Transport in Arterioles

Lucas Amoudruz  Computational Science and Engineering Laboratory, School of
Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138,
United States  Computational Science and Engineering Laboratory, ETH Zürich,
CH-8092, Switzerland  Athena Economides  Institute of Neuropathology,
University of Zurich, CH-8091 Zurich, Switzerland  Petros Koumoutsakos
Computational Science and Engineering Laboratory, School of Engineering and
Applied Sciences, Harvard University, Cambridge, MA 02138, United States
Corresponding author<EMAIL_ADDRESS>
###### Abstract

Red blood cells (RBCs) are vital for transporting oxygen from the lungs to the
body’s tissues through the intricate circulatory system. They achieve this by
binding and releasing oxygen molecules to the abundant hemoglobin within their
cytosol. The volume of RBCs affects the amount of oxygen they can carry, yet
whether this volume is optimal for transporting oxygen through the circulatory
system remains an open question. This study explores, through high-fidelity
numerical simulations, the impact of RBC volume on advectve oxygen transport
efficiency through arterioles which form the area of greatest flow resistance
in the circulatory system. The results show that, strikingly, RBCs with
volumes similar to those found in vivo are most efficient to transport oxygen
through arterioles. The flow resistance is related to the cell-free layer
thickness, which is influenced by the shape and the motion of the RBCs: at low
volumes RBCs deform and fold while at high volumes RBCs collide and follow
more diffuse trajectories. In contrast, RBCs with a healthy volume maximize
the cell-free layer thickness, resulting in a more efficient advectve
transport of oxygen.

Blood flow plays a vital role in sustaining life as it enables the delivery of
oxygen and nutrients to every tissue and organ in the body while also removing
harmful waste like carbon dioxide. RBC are a crucial component of blood,
making up almost half of the total blood volume in humans [37]. These
specialized cells are responsible for carrying oxygen from the lungs to the
rest of the body, where it is used to fuel cellular metabolism. They are
formed by a visco-elastic membrane that surrounds the cytosol of the cell. The
cytosol contains a high concentration of hemoglobin, which is responsible for
the transportation of oxygen [37, 10]. The red blood cell (RBC) membrane is
composed of a lipid bilayer anchored to a cytoskeleton and allows for
significant deformations [30]. The lipid bilayer of the membrane behaves like
an incompressible liquid crystal with localized bending resistance, while the
cytoskeleton provides resistance to local shear and dilation [22]. According
to current models, the unstressed shape of the cytoskeleton is an oblate
spheroid of reduced volume approximately $0.95$, compared to a sphere with the
same membrane area [30, 22, 3]. However, at equilibrium healthy RBCs adopt a
biconcave shape with a reduced volume ($v$) of approximately $0.65$ [14]. In
this paper, we study the effect of this reduced volume on the oxygen transport
efficiency in straight tubes with circular cross-section that have
representative characteristics of arterioles.

The biconcave shape and reduced volume of RBCs as well as metrics for their
optimality remain a subject of debate in the scientific community. Several
researchers proposed that the biconcave shape allows for a large surface-to-
volume ratio for RBCs so that they maximize oxygen exchange between the RBC
interior and their environment [29, 50]. The exchange of oxygen between the
cytoplasm and the environment is higher for cells with lower thickness because
cytoplasmic diffusion may be more limiting than membrane permeability, as
shown by Richardson et al. [41]. Other studies suggest that the biconcave
shape allows for a more efficient flow through capillaries and kidney tubules,
due to the RBCs’ deformability [30, 36, 51, 35]. Uzoigwe [50] proposed that
RBCs have a biconcave shape to maximize their moment of inertia, thereby
reducing shear stresses in blood flow and decreasing blood flow resistance in
arteries, but did not present any quantitative results. Guo et al. [23] showed
that the volume of cells affects their stiffness, while other studies found
that osmotic pressure and cell stiffness affect the viscosity of whole blood
[48, 47, 39, 43]. Alterations in RBC volume can significantly impact the
function of the circulatory system and individual health [50]. Farutin et al.
[15] studied numerically the cell transport efficiency in channels and
observed that the optimal hematocrit depends on the RBCs’ volume.

In this study, we investigate the effect of the reduced volume of RBCs on the
advectve oxygen transport efficiency along the blood flow direction. To this
end, we study the transport of a fixed quantity of hemoglobin by RBCs with
various resting shapes, parameterized by the reduced volume of the cells. We
investigate the effect of the RBC reduced volume on the oxygen flux through a
straight tube with circular cross-section. This metric has been used
experimentally to determine the optimal hematocrit for oxygen transport [47,
15]. We choose tubes with dimensions and pressure gradient similar to those
found in arterioles, where the flow resistance is the highest in the
circulatory system [7, 8, 33, 34].

The study relies on a recent and appropriately validated RBC model, comprising
visco-elastic membranes enclosing the cytosol and suspended in the blood
plasma. The model was extensively calibrated in Amoudruz et al. [3] and
validated against experimental data in various flow conditions. The evolution
of the solvent and the cytosol are described by dissipative particle dynamics
(DPD) [25, 13]. The numerical simulations are performed with Mirheo, a high-
performance software for blood flow and microfluidics [1].

## Methods

We model blood with RBCs composed of visco-elastic membranes surrounding their
viscous cytosol, and suspended in the blood plasma. The RBC membrane deforms
from forces that arise from bending resistance of the lipid-bilayer, as well
as due to the shear and dilation elasticity of the cytoskeleton with respect
to its stress-free state (SFS), and membrane viscosity. The resistance to
bending is described by the energy

$U_{bending}=2\kappa_{b}\oint{H^{2}dA},$ (1)

where the integral is taken over the membrane surface, $\kappa_{b}$ is the
bending modulus and $H$ is the local mean curvature. The in-plane elastic
energy is given by

$\begin{split}U_{in-
plane}=&\frac{K_{\alpha}}{2}\oint{\left(\alpha^{2}+a_{3}\alpha^{3}+a_{4}\alpha^{4}\right)dA_{0}}\\\
&+\mu\oint{\left(\beta+b_{1}\alpha\beta+b_{2}\beta^{2}\right)dA_{0}},\end{split}$

where the integral is taken over the SFS surface, $\alpha$ and $\beta$ are the
local dilation and shear strain invariants, respectively, $K_{\alpha}$ is the
dilation elastic modulus, $\mu$ is the shear elastic modulus and the
coefficients $a_{3}$, $a_{4}$, $b_{1}$ and $b_{2}$ are parameters that control
the non-linearity of the membrane elasticity for large deformations [31].

Each membrane is composed of 2562 particles located at the corners of a
triangulated mesh and evolve according to Newton’s law of motion. The bending
energy is discretized following references [26, 6], and the in-plane energy is
computed as described in Lim et al. [30]. The forces acting on the particles
are the negative gradients of the discretized energy terms, with respect to
the particle positions. Dissipation on the membrane is modeled by particles
sharing an edge in the triangle mesh that exert a pairwise force as described
in Fedosov [16], proportional to the membrane viscosity $\eta_{m}$. Finally,
the area of the membrane and the volume of the cytosol are constrained through
energy penalization terms,

$U_{area}=k_{A}\frac{\left(A-A_{0}\right)^{2}}{A_{0}},\quad
U_{volume}=k_{V}\frac{\left(V-V_{0}\right)^{2}}{V_{0}},$

where $A_{0}$ and $V_{0}$ are the area and volume of the cell at rest and $A$
and $V$ are the area and volume of the cell, respectively. These penalization
terms are required as the DPD method does not impose incompressibility of the
cytosol, and the membrane elastic energies alone do not conserve the membrane
area. The parameters of the model take physiological values calibrated from
experimental data [3] and are listed in the Appendix.

The RBC cytosol and surrounding plasma are represented with particles that
evolve through DPD interactions [25, 13, 3] (see Appendix for more details).
We emphasize that the model parameters correspond to the 1:5 ratio of
viscosities for the plasma ($1.2\text{\,}\mathrm{Pa}\text{\,}\mathrm{s}$ [53])
and the RBC cytosol ($6\text{\,}\mathrm{Pa}\text{\,}\mathrm{s}$ [52]). To
model the no-slip and no-flux boundary conditions on walls, particles are
bounced-back from the walls’ surface. Furthermore, the particles interact with
a layer of frozen particles that are inside the walls through DPD interactions
[40]. In addition, DPD particles are bounced-back from membrane surfaces, and
interact with the membrane particles through the dissipative and stochastic
parts of the DPD interactions only [16]. The cytosol and plasma particles
interact with each other only through the conservative part of the DPD forces.
The DPD parameters are chosen from the macroscopic properties of the fluids as
described in ref. [2].

## Results

Figure 1: Mean flow velocity in the tube against the reduced volume of the
RBCs for a fixed pressure gradient, with different tube radii, hematocrits and
pressure gradients.
: $R=$30\text{\,}\mathrm{\SIUnitSymbolMicro m}$$, $\mathrm{Hct}=40\%$, $\nabla
p=$1\text{\times}{10}^{-3}\text{\,}\mathrm{\text{cmH}_{2}\text{O}}\text{\,}{\mathrm{\SIUnitSymbolMicro
m}}^{-1}$$; : $R=$40\text{\,}\mathrm{\SIUnitSymbolMicro m}$$,
$\mathrm{Hct}=40\%$, $\nabla
p=$1\text{\times}{10}^{-3}\text{\,}\mathrm{\text{cmH}_{2}\text{O}}\text{\,}{\mathrm{\SIUnitSymbolMicro
m}}^{-1}$$; : $R=$40\text{\,}\mathrm{\SIUnitSymbolMicro m}$$,
$\mathrm{Hct}=45\%$, $\nabla
p=$1\text{\times}{10}^{-3}\text{\,}\mathrm{\text{cmH}_{2}\text{O}}\text{\,}{\mathrm{\SIUnitSymbolMicro
m}}^{-1}$$; : $R=$40\text{\,}\mathrm{\SIUnitSymbolMicro m}$$,
$\mathrm{Hct}=50\%$, $\nabla
p=$1\text{\times}{10}^{-3}\text{\,}\mathrm{\text{cmH}_{2}\text{O}}\text{\,}{\mathrm{\SIUnitSymbolMicro
m}}^{-1}$$; : $R=$40\text{\,}\mathrm{\SIUnitSymbolMicro m}$$,
$\mathrm{Hct}=45\%$, $\nabla
p=$2\text{\times}{10}^{-3}\text{\,}\mathrm{\text{cmH}_{2}\text{O}}\text{\,}{\mathrm{\SIUnitSymbolMicro
m}}^{-1}$$.  Figure 2: Snapshots of RBCs flowing near the walls of the tube
with $R=$40\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and $\mathrm{Hct}=0.45$.
The flow direction is from left to right. Cells with different reduced volumes
respond differently to the flow shear (from left to right: $v=0.35$, $0.65$,
$0.95$).

We examine the equilibrated flow of RBCs inside a circular tube of radius
$R\in\\{$30\text{\,}\mathrm{\SIUnitSymbolMicro
m}$,$40\text{\,}\mathrm{\SIUnitSymbolMicro m}$\\}$. These sizes are typical
for arterioles [21]. The volume fraction of RBCs, or tube hematocrit, is set
to $\mathrm{Hct}\in\\{40\%,45\%,50\%\\}$ and the blood suspension is subjected
to a pressure gradient $\nabla
p\in\\{$1\text{\times}{10}^{-3}\text{\,}\mathrm{\text{cmH}_{2}\text{O}}\text{\,}{\mathrm{\SIUnitSymbolMicro
m}}^{-1}$,$2\text{\times}{10}^{-3}\text{\,}\mathrm{\text{cmH}_{2}\text{O}}\text{\,}{\mathrm{\SIUnitSymbolMicro
m}}^{-1}$\\}$, typical for arterioles of this size [21]. These conditions
correspond to mean flow velocities of the order of
$10\text{\,}\mathrm{mm}\text{\,}{\mathrm{s}}^{-1}$, consistent with values
reported in the literature [37]. Mean flow velocities of all simulations are
reported on fig. 1. We study different cases where all RBC membranes have the
same area and visco-elastic properties (calibrated from experimental data [3])
but have a different volume $V_{0}$ parameterized by the reduced volume
$v=V_{0}/V_{s}$, where $V_{s}$ is the volume of a sphere with the same area as
the cell’s membrane. All initial triangle mesh were produced by minimizing the
membrane energy with a target volume $V_{0}$. The number of RBCs is adapted to
keep the hematocrit constant, $N=\lfloor\pi R^{2}L\mathrm{Hct}/V_{0}\rfloor$.
We study a periodic domain with a length
$L=$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ along the flow direction. No-
slip and no-flux boundary conditions are applied on the side walls. All
quantities are reported from equilibrated blood flows.

We find that the reduced volume affects the flow patterns and RBCs
deformability (fig. 2). Specifically, for low reduced volumes, the cells
deform much more than at high values of $v$. In the limit $v\rightarrow 1$,
only spherical cells are allowed, and as the area is constant they cannot
deform. RBC with a physiological reduced volume do not exhibit significant
deformations compared to those with lower reduced volumes. Instead, RBCs with
$v=0.65$ seem to keep a relatively flat shape similar to tank-treading RBCs, a
type of motion where the membrane rotates around a steady elongated shape
similar to tank treads [18].

We quantify the amount of transported oxygen by computing the flux of cytosol
in the tube for each value of $v$,

$Q_{c}=\frac{V_{0}}{L}\sum\limits_{i=1}^{N}U_{i},$

where $U_{i}$ is the time-averaged velocity of the $i^{\text{th}}$ RBC’s
center of mass. We assume that the hemoglobin concentration inside the cytosol
is constant, and therefore the cytosol flux is proportional to the oxygen
flux. Figure 3 shows the flux of cytosol against the reduced volume of the
RBCs for different tube radii and hematocrits. The flux of cytosol is
normalized by the flux of plasma with no RBCs, $Q_{p}=|\nabla p|\pi
R^{4}/8\eta$, where $\eta$ is the dynamic viscosity of the plasma. In these
conditions, $Q_{c}$ reaches a maximum at a reduced volume $v\approx 0.65$.
Strikingly, the oxygen flux is maximized at the physiological reduced volume
of RBCs.

Figure 3: Cytosol flux $Q_{c}$, normalized by that of the plasma at zero
hematocrit, $Q_{p}$, against the reduced volume of RBCs $v$ flowing in the
tube of different radii and hematocrits. Same labels as in fig. 1.  Figure 4:
Snapshots of RBCs flowing near the walls (top horizontal line) of the tube.
The flow direction is from left to right. The slice of the RBCs in the
vertical plane is represented with black lines.  Figure 5: Cell-free layer
thickness $\delta$ against the reduced volume of RBCs $v$ flowing in the tube.
Same labels as in fig. 1.

To elucidate the mechanisms that reduce the cytosol flux when the volume of
the RBCs deviates from its physiological value, we examine the cell-free layer
(CFL) near the walls of the tube (fig. 4). We report its thickness, $\delta$,
with respect to the reduced volume $v$ on fig. 5. The viscosity is lower in
the CFL than in the blood suspension, hence the shear rate is higher in the
CFL region. Therefore, for a given pressure difference, a larger CFL thickness
causes a higher flux of hemoglobin (figs. 3 and 5) [9]. We suggest that the
change in the CFL thickness depends on the reduced volume of the RBCs. It has
been reported in the literature that cells that are close to the walls in a
Poiseuille flow tend to migrate towards the center of the pipe due to their
deformability [44]. In contrast, a large volume fraction of RBCs may cause the
cells to move closer to the walls due to collisions between the cells. In the
remainder of this study we investigate the deformations of RBCs and their
trajectories and correlate these observations to the CFL thickness. The
deformation of the cells are quantified by the rim angle (see below) and
bending energy of each cell. We distinguish tumbling-like motion, due to the
external shear stresses, from tank-treading motion, where the membrane rotates
around a steady shape. The cell trajectories are quantified by the diffusion
coefficient along the radial direction. We fix
$R=$40\text{\,}\mathrm{\SIUnitSymbolMicro m}$$, $\nabla
p=$1\text{\times}{10}^{-3}\text{\,}\mathrm{\text{cmH}_{2}\text{O}}\text{\,}{\mathrm{\SIUnitSymbolMicro
m}}^{-1}$$ and $\mathrm{Hct}=0.45$ for clarity, but similar results were
observed for the other conditions listed in fig. 1.

#### Rim angle

Figure 6: Top left: The rim (white) of a RBC at rest. Bottom left: The normal
of the rim plane (dark arrow) and the normal of the cell plane (red arrow)
form the rim angle $\theta$. Right: Violin plots illustrating the variation in
rim angle distribution $\theta$ (in degrees) of flowing cells across different
reduced volumes $v$. Crosses indicate the mean value of $\theta$.

In a sufficiently high shear, the RBC membrane rotates around the cytosol
[18]. In this situation, the region of the membrane that forms the rim at
equilibrium differs from the principal disc of the deformed cell, unlike RBCs
at rest. We define the rim angle as the angle between the normal to the
principal disc of the cell and the normal to the region corresponding to the
rim of the stress-free state of the cell, see fig. 6 and ref. [12]. For a
tumbling RBC, this angle is typically low while for tank-treading RBCs the
angle takes values in the whole range, $\theta\in[0,\pi/2]$. Figure 6 shows
the distribution of rim angles $\theta$ of cells flowing in the tube for
different values of $v$. On average, the rim angles are smaller at low values
of $v$. Furthermore, the density of $\theta$ is more uniform at high values of
$v$ than at low reduced volume values. This suggests that RBCs with a large
reduced volume tank-tread while their motion is closer to tumbling at lower
values of $v$.

#### Bending energy

Figure 7:  Left: Violin plots illustrating the variation in bending energy
distribution of flowing cells across different reduced volumes $v$. Crosses
indicate the mean value of the bending energy, shaded regions indicate the
maximum density. The dashed line corresponds to the bending energy of the
cells at rest. $\kappa_{b}$ is the bending modulus. Right: Time sequences of
typical RBCs in the tube (time increasing from left to right). From top to
bottom: $v=0.4$, $v=0.65$ and $v=0.95$. The white stripe represents the rim of
the cell at rest (see fig. 6).

The deformations of the cells are characterized by the bending energy
described by eq. 1. This energy is high when the cell deforms (e.g. folds) and
therefore characterizes complex dynamics and shape changes. The bending energy
distribution of the cells is shown on fig. 7 for different values of $v$. The
mean bending energy decreases when $v$ increases, consistent with the fact
that the minimum bending energy is achieved for a spherical shape.
Nevertheless, we observe differences in the distributions of the bending
energies. For $v\geq 0.6$, the bending energies are concentrated close to a
minimal value. Instead, when $v<0.6$, we observe a peak of bending energy well
above the minimal values observed at each reduced volume. These high values
suggest that RBCs undergo large deformations compared to their resting shape.
In particular, we observe cells that are “folding” periodically when $v$ is
relatively small (fig. 7), as opposed to cells that have a large reduced
volume. These large deformations cause a large bending energy density in the
bulk, or, equivalently, a large pressure compared to configurations with
undeformed cells. The larger deformations associated with low reduced volumes
may thus contribute to reducing the CFL thickness in this regime. Furthermore,
the deformation of RBCs is accompanied by an additional dissipation due to the
membrane viscosity and the recirculation of the cytosol inside the cell [18],
which may further decrease the blood flux.

#### Tank-treading frequency

Figure 8:  Left: TTF of the RBCs against their time averaged radial position,
for $v=0.4$ (circles), $v=0.7$ (triangles) and $v=0.9$ (squares). Right:
Radial diffusion $D$ of the RBCs against the reduced volume $v$.

In linear shear flows, RBCs tank-tread when the shear rate is large enough
[19, 11]. The tank-treading frequency (TTF) is reported against the time
averaged radial position of the cells on fig. 8 (see appendix for details).
For low values of $v$, most cells do not tank-tread hence the value estimated
for the TTF is low. For larger values of $v$, the TTF increases linearly with
the radial position. This increase is consistent with the nearly linear
relationship between the TTF and the shear rate for single cells [19], and
assuming that we have a nearly parabolic velocity profile along the radial
direction of the tube, hence a linear increase of the shear rate with the
radial position $r$. Furthermore, we observe that the TTF of cells is larger
when $v$ increases. This is also expected since spheres in shear flows rotate
with a larger frequency than the TTF of RBCs [19]. We suggest that, for
relatively thin cells, tank-treading is beneficial for cells to pass each
other without deforming and to minimize the effect of their collisions, thus
contributing to a lower value of $\delta$. It was shown that tank-treading
provides less hindrance to the motion of other particles compared to other
motions like tumbling [28]. This may be an indication that tank-treading cells
contribute to a lower value of $\delta$ compared to other motions, present
notably at low values of $v$.

#### Radial diffusion

RBC in blood flow undergo shear-induced diffusion, which is a result of the
cell-cell interactions. These interactions influence the radial distribution
of the cells and the CFL thickness [24, 46]. We thus measure the radial
diffusion coefficient of the RBCs in the tube. The diffusion coefficient is
estimated from the mean squared displacement of the cells along the radial
direction:

$D=\langle\left(r(T)-r(0)\right)^{2}\rangle/2T,$

where $r$ is the radial position of the cell’s center of mass and $T$ is a
time large enough to collect statistics but small enough to remain in the
linear regime of the mean squared displacement against time.

The theory of shear-induced diffusion suggests that, at a given volume
fraction, $D\propto\dot{\gamma}a^{2}$, where $\dot{\gamma}$ is the shear rate
and $a$ is the size of the particle [42]. Assuming that the shear rate is
constant across the simulations and $a\propto v^{1/3}$, we expect $Dv^{3/2}$
to be approximately constant against $v$. This quantity is reported on fig. 8.
We observe that $Dv^{3/2}$ remains constant within a 30% deviation, suggesting
that increasing $v$ increases the shear-induced diffusion. The deviations from
the theory may come from the cells deformability, the definition of the cell
radius from $v$ and the non-constant shear rate along the radial direction.
The shear-induced diffusion is higher for larger values of $v$ and may thus
contribute to reducing the CFL thickness for $v>0.65$ [46], hence decreasing
the transport efficiency of the cytosol.

## Discussion

The physiological volume of RBCs, under the assumptions made in this study, is
the most efficient to transport oxygen through tubes that have the
characteristics of arterioles. When the reduced volume of the RBCs deviates
from its physiological value, we identify two factors that decrease the
transport efficiency of oxygen: the change of the CFL thickness, and the
dissipation due to the cytosol recirculation and the deformations of the
cells.

The viscosity of plasma is smaller than that of whole blood, thus a larger CFL
thickness contributes to a higher flux of blood. We find that the CFL
thickness is maximized at $v\approx 0.65$. Below this value, cells fold,
undergo large deformations and adopt a tumbling rather than tank-treading
motion. These deformations increase the effective thickness of the cells along
the radial direction, thus reducing the CFL thickness. In contrast, at larger
values of $v$, cells do not deform and have more diffusive trajectories, due
to collisions with their neighboring cells. These collisions cause the cells
to migrate perpendicularly to the flow direction, thus reducing the CFL
thickness. At intermediate values of $v$, the cells tank-tread, which
facilitate their motion relative to each other. This is supported by the
shear-thinning behavior of blood when the motion of RBCs transitions from
tumbling to tank-treading [17, 20]. The non-spherical shape of RBCs limits the
effect of collisions, and they form a maximal CFL thickness.

Figure 9: Cytosol flux against the CFL thickness. Circles correspond to cells
with $v<0.65$, triangles those with $v>0.65$ and the square is for $v=0.65$.

The CFL thickness alone does not explain the profile observed in fig. 3.
Indeed, for the same value of CFL thickness, the cytosol flux is lower for the
case with $v<0.65$ than when $v>0.65$ (fig. 9). We attribute this difference
to the additional dissipation due to the large deformations occurring for
cells with smaller reduced volumes, as membrane deformation was shown to be a
dissipative process [49]. In addition, we expect that tank-treading cells
dissipate more energy at lower reduced volumes, as shown for vesicles [27].

In this study we tackled the following question: what is the reduced volume
that corresponds to the highest flux of a given amount of hemoglobin? We
achieved this by adjusting the number of RBCs, assuming that the concentration
of hemoglobin in the cytosol is the same in every situation. However, we could
also consider a hemoglobin concentration that varies with $v$, e.g. a higher
concentration at low values of $v$ to model erythropoiesis [32], but this
approach would require to model the dependence of the cytosol viscosity with
the hemoglobin concentration. These different concentrations would change the
volume of oxygen transported by the system, and the viscosity of the cytosol.
In this study we did not explore these effects due to the large computational
cost of the simulations and made the approximation that the hemoglobin remains
constant across all cases, which corresponds to the simplest model. We leave
the more complex scenario of varying the concentration as well as the number
of RBCs to further studies.

The results of this study are valid for straight tubes, which is a simplistic
approximation of arterioles. The CFL thickness is known to vary around
bifurcations [38, 4, 54, 5] or curved geometries [45]. However, the
deformation of cells at low $v$ and the collisions of cells at larger values
of $v$ might still contribute to reduce the CFL thickness in those conditions
compared to the physiological values, and we thus expect the same trend as in
straight tubes. Nevertheless, the optimal value of $v$ might slightly differ
in those cases and further studies are required to provide a quantitative
description of these variations.

Finally, we remark that we considered the advective part of oxygen transport
in arterioles. The diffusive transport of oxygen, not studied in this work, is
responsible for transporting oxygen across the RBC membranes. This type of
transport may be affected by the shape, and thus the reduced volume, of RBCs
[29, 50, 41].

## Conclusion

We demonstrate, through high fidelity simulations, that RBCs with volumes
similar to those observed in vivo maximize the efficiency of advective oxygen
transport in straight tubes that have characteristics of arterioles, where the
vascular resistance is maximal. We qualitatively explain the variation of
advective oxygen transport efficiency with respect to the reduced volume of
the cells based on the CFL thickness and the dissipation due to the cells
deformability. At low reduced volumes, cells deform and fold, thus occupy a
larger effective volume on average. This causes the CFL thickness to decrease.
The deformations of cells induce additional dissipation due to the
recirculation of the cytosol and the highly viscous membranes. At large
reduced volumes, cells collide and migrate in the directions perpendicular to
the flow due to a high shear-induced diffusion. These trajectories also
contribute to reducing the CFL thickness and thus the oxygen transport is
lower.

The present findings provide valuable insights into the mechanisms of
advective oxygen transport in the body and could potentially have significant
implications for the advancement of therapies aimed at treating circulatory
disorders.

## Appendix

### Dissipative Particle Dynamics

The DPD method discretizes a fluid into $N$ particles with positions
$\mathbf{r}_{i}$, velocities $\mathbf{v}_{i}$ and mass $m$, $i=1,2,\dots,N$.
The particles evolve according to the Newton’s law of motion,

$\dot{\mathbf{r}}_{i}=\mathbf{v}_{i},\;\;\;\dot{\mathbf{v}}_{i}=\frac{1}{m}\sum\limits_{j=1}^{N}\mathbf{F}_{ij},$

where $\mathbf{F}_{ij}$ are pairwise forces that vanish after a cutoff
distance $r_{c}$. These interactions are formed by three terms [25, 13],

$\mathbf{F}_{ij}=aw(r_{ij})\mathbf{e}_{ij}-\gamma\left(\mathbf{e}_{ij}\cdot\mathbf{v}_{ij}\right)w_{D}(r_{ij})\mathbf{e}_{ij}+\sigma\xi_{ij}w_{R}(r_{ij})\mathbf{e}_{ij},$

where $\mathbf{v}_{ij}=\mathbf{v}_{i}-\mathbf{v}_{j}$,
$\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$, $r_{ij}$ =
$\|\mathbf{r}_{ij}\|$ and $\mathbf{e}_{ij}=\mathbf{r}_{ij}/r_{ij}$. The
coefficients $a$, $\gamma$ and $\sigma$ are the conservative, dissipative and
random force magnitudes, respectively. Furthermore, we use the standard DPD
conservative kernel

$w(r)=\begin{cases}1-r/r_{c},&r<r_{c},\\\ 0,&\text{otherwise}.\end{cases}$

We set $w_{D}=w^{1/4}$. $w_{R}$ satisfies the fluctuation-dissipation theorem,
$\sigma^{2}=2\gamma k_{B}T$ and $w_{D}=w_{R}^{2}$ [13], where $k_{B}T$ is the
temperature of the system in energy units.

### Parameters of the model

The values of the RBC parameters are listed in table 1. They correspond to the
mean value of the posterior distribution found in ref [3].

Parameter | Value
---|---
$\kappa_{b}$ | $2.10\text{\times}{10}^{-19}\text{\,}\mathrm{J}$
$\mu$ | $4.99\text{\,}\mathrm{\SIUnitSymbolMicro N}\text{\,}{\mathrm{m}}^{-1}$
$K_{\alpha}$ | $4.99\text{\,}\mathrm{\SIUnitSymbolMicro N}\text{\,}{\mathrm{m}}^{-1}$
$a_{3}$ | $-2$
$a_{4}$ | $8$
$b_{1}$ | $0.7$
$b_{2}$ | $1.84$
$\eta_{m}$ | $0.42\text{\times}{10}^{-6}\text{\,}\mathrm{Pa}\text{\,}\mathrm{s}\text{\,}\mathrm{m}$
$A_{0}$ | $135\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2}$
$V_{0}$ | varying
$k_{A}$ | $0.5\text{\,}\mathrm{J}\mathrm{/}\mathrm{m}^{2}$
$k_{V}$ | $7.23\text{\times}{10}^{5}\text{\,}\mathrm{J}\mathrm{/}\mathrm{m}^{3}$
Table 1: Parameters of the RBC model.

### Computation of the cell-free layer thickness

The cell-free layer thickness $\delta$ is computed from the vertices of the
RBC membranes. For a given time snapshot, we compute the maximum radial
position of vertices over bins placed at regular intervals along the flow
direction. The bins have a size of $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$.
The value of $\delta$ at a given time is then the difference between the
radius of the pipe and the average of this quantity over all bins. The values
reported in this work are the average of the cell-free layer thickness over
time, after equilibration of the flow.

### Computation of the tank-treading frequency

The TTF is estimated by computing the Fourier transform of the rim angle time
series of each cell and selecting the frequency of the largest mode, divided
by 2 since the rim angle of a tank-treading cell undergoes 2 revolutions per
tank-treading revolution. Note that we exclude cells having an average rim
angle $\bar{\theta}$ larger than $\pi/4$, and the tumbling cells, which we
characterize by $\bar{\theta}<\pi/8$.

## Author contributions

LA, AE and PK designed research. LA performed research. LA and AE analyzed
data. LA, AE and PK wrote the manuscript.

## Declaration of interests

The authors declare no competing interests.

## Acknowledgments

We would like to thank Xin Bian for the idea of using the rim angle to perform
the analysis. We acknowledge the computational resources granted by the Swiss
National Supercomputing Center (CSCS) under the project ID “s1160”.

## References

  * [1] Dmitry Alexeev, Lucas Amoudruz, Sergey Litvinov, and Petros Koumoutsakos.  Mirheo: High-performance mesoscale simulations for microfluidics.  Computer Physics Communications, 254:107298, 2020.
  * [2] Lucas Amoudruz.  Simulations and Control of Artificial Microswimmers in Blood.  PhD thesis, ETH Zurich, 2022.
  * [3] Lucas Amoudruz, Athena Economides, Georgios Arampatzis, and Petros Koumoutsakos.  The stress-free state of human erythrocytes: Data-driven inference of a transferable RBC model.  Biophysical Journal, 122(8):1517–1525, apr 2023.
  * [4] Christian Bächer, Alexander Kihm, Lukas Schrack, Lars Kaestner, Matthias W Laschke, Christian Wagner, and Stephan Gekle.  Antimargination of microparticles and platelets in the vicinity of branching vessels.  Biophysical journal, 115(2):411–425, 2018.
  * [5] Peter Balogh and Prosenjit Bagchi.  The cell-free layer in simulated microvascular networks.  Journal of fluid mechanics, 864:768–806, 2019.
  * [6] Xin Bian, Sergey Litvinov, and Petros Koumoutsakos.  Bending models of lipid bilayer membranes: Spontaneous curvature and area-difference elasticity.  Comput. Method. Appl. M., 359:112758, feb 2020.
  * [7] H Glenn Bohlen, Robert W Gore, and Phillip M Hutchins.  Comparison of microvascular pressures in normal and spontaneously hypertensive rats.  Microvascular research, 13(1):125–130, 1977.
  * [8] WILLIAM M Chilian, CHARLES L Eastham, and MELVIN L Marcus.  Microvascular distribution of coronary vascular resistance in beating left ventricle.  American Journal of Physiology-Heart and Circulatory Physiology, 251(4):H779–H788, 1986.
  * [9] Giles R Cokelet and Harry L Goldsmith.  Decreased hydrodynamic resistance in the two-phase flow of blood through small vertical tubes at low flow rates.  Circulation research, 68(1):1–17, 1991.
  * [10] Giuseppe Di Caprio, Chris Stokes, John M Higgins, and Ethan Schonbrun.  Single-cell measurement of red blood cell oxygen affinity.  Proceedings of the National Academy of Sciences, 112(32):9984–9989, 2015.
  * [11] Jules Dupire, Marius Socol, and Annie Viallat.  Full dynamics of a red blood cell in shear flow.  Proceedings of the National Academy of Sciences, 109(51):20808–20813, 2012.
  * [12] Eva Athena Economides.  Data informed, predictive simulations of blood microfluidics.  PhD thesis, ETH Zurich, 2020.
  * [13] Pep Espanol and Patrick Warren.  Statistical mechanics of dissipative particle dynamics.  EPL (Europhysics Letters), 30(4):191, 1995.
  * [14] Evan Evans and Yuan Cheng Fung.  Improved measurements of the erythrocyte geometry.  Microvascular Research, 4(4):335–347, 1972.
  * [15] Alexander Farutin, Zaiyi Shen, Gael Prado, Vassanti Audemar, Hamid Ez-Zahraouy, Abdelilah Benyoussef, Benoit Polack, Jens Harting, Petia M Vlahovska, Thomas Podgorski, et al.  Optimal cell transport in straight channels and networks.  Physical review fluids, 3(10):103603, 2018.
  * [16] Dmitry A Fedosov.  Multiscale modeling of blood flow and soft matter.  PhD thesis, Citeseer, 2010.
  * [17] Dmitry A Fedosov, Wenxiao Pan, Bruce Caswell, Gerhard Gompper, and George E Karniadakis.  Predicting human blood viscosity in silico.  Proceedings of the National Academy of Sciences, 108(29):11772–11777, 2011.
  * [18] Th M Fischer.  On the energy dissipation in a tank-treading human red blood cell.  Biophysical journal, 32(2):863–868, 1980.
  * [19] Thomas M Fischer.  Tank-Tread Frequency of the Red Cell Membrane : Dependence on the Viscosity of the Suspending Medium.  Biophysical Journal, 93(7):2553–2561, 2007.
  * [20] Alison M Forsyth, Jiandi Wan, Philip D Owrutsky, Manouk Abkarian, and Howard A Stone.  Multiscale approach to link red blood cell dynamics, shear viscosity, and atp release.  Proceedings of the National Academy of Sciences, 108(27):10986–10991, 2011.
  * [21] Yuan-Cheng Fung.  Biomechanics: circulation.  Shock, 9(2):155, 1998.
  * [22] Nadeeshani Maheshika Geekiyanage, Marie Anne Balanant, Emilie Sauret, Suvash Saha, Robert Flower, Chwee Teck Lim, and YuanTong Gu.  A coarse-grained red blood cell membrane model to study stomatocyte-discocyte-echinocyte morphologies.  PLoS One, 14(4):e0215447, 2019.
  * [23] Ming Guo, Adrian F Pegoraro, Angelo Mao, Enhua H Zhou, Praveen R Arany, Yulong Han, Dylan T Burnette, Mikkel H Jensen, Karen E Kasza, Jeffrey R Moore, et al.  Cell volume change through water efflux impacts cell stiffness and stem cell fate.  Proceedings of the National Academy of Sciences, 114(41):E8618–E8627, 2017.
  * [24] John M Higgins, David T Eddington, Sangeeta N Bhatia, and Lakshminarayanan Mahadevan.  Statistical dynamics of flowing red blood cells by morphological image processing.  PLoS computational biology, 5(2):e1000288, 2009.
  * [25] PJ Hoogerbrugge and JMVA Koelman.  Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics.  EPL (Europhysics Letters), 19(3):155, 1992.
  * [26] Frank Jülicher.  The morphology of vesicles of higher topological genus: conformal degeneracy and conformal modes.  Journal de Physique II, 6(12):1797–1824, 1996.
  * [27] Martin Kraus, Wolfgang Wintz, Udo Seifert, and Reinhard Lipowsky.  Fluid vesicles in shear flow.  Physical Review Letters, 77(17):3685–3688, 1996.
  * [28] Timm Krüger, Markus Gross, Dierk Raabe, and Fathollah Varnik.  Crossover from tumbling to tank-treading-like motion in dense simulated suspensions of red blood cells.  Soft Matter, 9(37):9008–9015, 2013.
  * [29] JG Lenard.  A note on the shape of the erythrocyte.  Bulletin of Mathematical Biology, 36(1):55–58, 1974.
  * [30] Gerald Lim H. W., Michael Wortis, and Ranjan Mukhopadhyay.  Red Blood Cell Shapes and Shape Transformations: Newtonian Mechanics of a Composite Membrane.  Soft Matter, 4, 2008.
  * [31] Gerald Lim HW, Michael Wortis, and Ranjan Mukhopadhyay.  Stomatocyte–discocyte–echinocyte sequence of the human red blood cell: Evidence for the bilayer–couple hypothesis from membrane mechanics.  Proceedings of the National Academy of Sciences, 99(26):16766–16769, 2002.
  * [32] Y Matoth, R Zaizov, and I Varsano.  Postnatal changes in some red cell parameters.  Acta Paediatrica, 60(3):317–323, 1971.
  * [33] Gerald A Meininger.  Responses of sequentially branching macro-and microvessels during reactive hyperemia in skeletal muscle.  Microvascular research, 34(1):29–45, 1987.
  * [34] MJ Mulvany and C Aalkjaer.  Structure and function of small arteries.  Physiological reviews, 70(4):921–961, 1990.
  * [35] Arman Namvar, Adam J Blanch, Matthew W Dixon, Olivia MS Carmo, Boyin Liu, Snigdha Tiash, Oliver Looker, Dean Andrew, Li-Jin Chan, Wai-Hong Tham, et al.  Surface area-to-volume ratio, not cellular viscoelasticity, is the major determinant of red blood cell traversal through small channels.  Cellular microbiology, 23(1):e13270, 2021.
  * [36] Igor V Pivkin, Zhangli Peng, George E Karniadakis, Pierre A Buffet, Ming Dao, and Subra Suresh.  Biomechanics of red blood cells in human spleen and consequences for physiology and disease.  Proceedings of the National Academy of Sciences, 113(28):7804–7809, 2016.
  * [37] Aleksander S Popel.  Theory of oxygen transport to tissue.  Critical reviews in biomedical engineering, 17(3):257, 1989.
  * [38] Yazdan Rashidi, Greta Simionato, Qi Zhou, Thomas John, Alexander Kihm, Mohammed Bendaoud, Timm Krüger, Miguel O Bernabeu, Lars Kaestner, Matthias W Laschke, et al.  Red blood cell lingering modulates hematocrit distribution in the microcirculation.  Biophysical Journal, 122(8):1526–1537, 2023.
  * [39] WH Reinhart, M Singh-Marchetti, and PW Straub.  The influence of erythrocyte shape on suspension viscosities.  European journal of clinical investigation, 22(1):38–44, 1992.
  * [40] M Revenga, I Zuniga, P Espanol, and I Pagonabarraga.  Boundary models in dpd.  International Journal of Modern Physics C, 9(08):1319–1328, 1998\.
  * [41] Sarah L Richardson, Alzbeta Hulikova, Melanie Proven, Ria Hipkiss, Magbor Akanni, Noémi BA Roy, and Pawel Swietach.  Single-cell o2 exchange imaging shows that cytoplasmic diffusion is a dominant barrier to efficient gas transport in red blood cells.  Proceedings of the National Academy of Sciences, 117(18):10067–10078, 2020.
  * [42] Pierre Rognon and Matthew Macaulay.  Shear-induced diffusion in dense granular fluids.  Soft matter, 17(21):5271–5277, 2021.
  * [43] H Schmid-Schönbein and RE Wells Jr.  Rheological properties of human erythrocytes and their influence upon the “anomalous” viscosity of blood.  In Ergebnisse der Physiologie Reviews of Physiology, Volume 63, pages 146–219. Springer, 2010.
  * [44] Lingling Shi, Tsorng-Whay Pan, and Roland Glowinski.  Numerical simulation of lateral migration of red blood cells in poiseuille flows.  International Journal for Numerical Methods in Fluids, 68(11):1393–1408, 2012.
  * [45] Christian J Spieker, Gábor Závodszky, Clarisse Mouriaux, Max Van der Kolk, Christian Gachet, Pierre H Mangin, and Alfons G Hoekstra.  The effects of micro-vessel curvature induced elongational flows on platelet adhesion.  Annals of Biomedical Engineering, 49:3609–3620, 2021.
  * [46] A Srivastav, X Grandchamp, G Coupier, and T Podgorski.  Shear-induced diffusion in a red blood cell suspension.  Computer methods in biomechanics and biomedical engineering, 15(sup1):34–35, 2012.
  * [47] HO Stone, HK Thompson Jr, and K Schmidt-Nielsen.  Influence of erythrocytes on blood viscosity.  American Journal of Physiology-Legacy Content, 214(4):913–918, 1968\.
  * [48] Max M Strumia, Martha Phillips, Albert B Sample, M Elizabeth Burns, and Pasquale Mariano.  Effect of red cell factors on the relative viscosity of whole blood.  American Journal of Clinical Pathology, 39(5):464–474, 1963.
  * [49] Roger Tran-Son-Tay, SP Sutera, and PR Rao.  Determination of red blood cell membrane viscosity from rheoscopic observations of tank-treading motion.  Biophysical journal, 46(1):65–72, 1984.
  * [50] Chika Uzoigwe.  The human erythrocyte has developed the biconcave disc shape to optimise the flow properties of the blood in the large vessels.  Medical hypotheses, 67(5):1159–1163, 2006.
  * [51] Koohyar Vahidkhah, Peter Balogh, and Prosenjit Bagchi.  Flow of red blood cells in stenosed microvessels.  Scientific reports, 6(1):1–15, 2016.
  * [52] R Wells and H Schmid-Schönbein.  Red cell deformation and fluidity of concentrated cell suspensions.  Journal of Applied Physiology, 27(2):213–217, 1969.
  * [53] Roe E Wells Jr and Edward W Merrill.  Shear rate dependence of the viscosity of whole blood and plasma.  Science, 133(3455):763–764, 1961.
  * [54] Qi Zhou, Joana Fidalgo, Miguel O Bernabeu, Mónica SN Oliveira, and Timm Krüger.  Emergent cell-free layer asymmetry and biased haematocrit partition in a biomimetic vascular network of successive bifurcations.  Soft Matter, 17(13):3619–3633, 2021.

  *[RBC]: red blood cell
  *[RBCs]: red blood cell
  *[DPD]: dissipative particle dynamics
  *[SFS]: stress-free state
  *[CFL]: cell-free layer
  *[TTF]: tank-treading frequency


# System size dependence of thermodynamic variables at kinetic freeze-out in
high-energy collisions using the Tsallis distribution

Rishabh Sharma<EMAIL_ADDRESS>Krishan Gopal
<EMAIL_ADDRESS>Sharang Rav Sharma
<EMAIL_ADDRESS>Chitrasen Jena
<EMAIL_ADDRESS>Department of Physics, Indian Institute of Science
Education and Research (IISER) Tirupati, Tirupati 517507, Andhra Pradesh,
India

###### Abstract

We use a thermodynamically consistent form of Tsallis distribution to study
the dependence of various thermodynamic quantities on the system size in high-
energy collisions. The charged hadron spectra obtained in $p$+$p$, $p$+Pb,
Xe+Xe, and Pb+Pb collisions at LHC are used to determine the energy density,
pressure, particle density, entropy density, mean free path, Knudsen number,
heat capacity, isothermal compressibility, expansion coefficient, and speed of
sound at the kinetic freeze-out surface. These quantities are studied as a
function of the system size. Notably, the rate of increase (or decrease) in
these thermodynamic variables is found to be more rapid in small systems such
as $p$+$p$ and $p$+Pb collisions than in large systems such as Xe+Xe and Pb+Pb
collisions. This may be due to the small volume of the hadronic system in
small collision systems at kinetic freeze-out. It is observed that high-
multiplicity $p$+$p$ collisions produce similar thermodynamic conditions as
peripheral heavy-ion collisions at kinetic freeze-out.

Suggested keywords

††preprint: APS/123-QED

##  I Introduction

High-energy collisions provide a unique opportunity to study the behavior of
matter at extraordinarily high temperatures and densities. In the last few
decades, considerable progress has been made in the study of the properties of
a new state of matter made up of de-confined quarks and gluons called Quark-
Gluon Plasma (QGP). This state is believed to have formed when normal baryonic
matter undergoes a first-order or crossover phase transition at very high
temperature and/or high baryon density [1, 2, 3]. Previously, it was believed
that QGP is only formed in heavy-ion collisions; however, recent results show
that high-multiplicity $p$+$p$ collisions behave similarly to heavy-ion
collisions, posing the question of whether QGP droplets are formed in such
collisions [4, 5, 6, 7]. Consequently, there is a pressing need to gain deeper
insights into high-multiplicity $p$+$p$ collisions in order to understand the
properties of the produced system.

To better understand the behavior of matter resulting from high-energy
collisions, it is crucial to comprehensively understand its thermodynamic
properties [8, 9, 10, 11, 12, 13, 14, 15]. In this study, we have explored
various thermodynamic properties of the hadronic medium produced in both small
collision systems ($p$+$p$ and $p$+Pb) and large collision systems (Xe+Xe and
Pb+Pb) at LHC energies.

We examine a scenario in which the fireball resulting from high-energy
collisions cools down as it undergoes spacetime expansion. The final-state
particles decouple from the system after Kinetic Freeze-Out (KFO), which
represents the surface of the last elastic scattering of the produced hadrons.
The production of final-state particles in abundance warrants the use of
statistical methods to describe such systems. One such approach involves
studying the transverse momentum ($p_{T}$) spectra of the produced particles
in the framework of Tsallis statistics. Tsallis distribution was initially
introduced as a generalization of the Boltzmann-Gibbs distribution and it
accounts for deviations from the standard equilibrium statistical description
[16]. Over time, Tsallis distribution has proven to be a valuable tool for
describing the hadron spectra in high-energy collisions. The Tsallis
distribution function is defined as [17]:

$f(E,q,T,\mu)=\left[1+(q-1)\frac{E-\mu}{T}\right]^{\frac{-1}{q-1}},$ (1)

where $E$ is the energy of the particle, $q$ is the entropy index, $T$ is the
temperature, and $\mu$ is the chemical potential. $q$ is indicative of the
degree of deviation from equilibrium, with $q=1$ representing complete
equilibrium. Although various forms of the Tsallis distribution have been used
to study the hadron spectra, ensuring thermodynamic consistency in their
formulations poses a significan challenge [18, 19, 20, 21]. In the present
study, we have employed Eq. (1), which has been demonstrated to maintain
thermodynamic consistency [17]. The thermodynamic quantities of interest were
derived using Eq. (1) in accordance with the following relations [22, 23, 24]:

$\displaystyle\text{Entropy: }S$
$\displaystyle=sV=-gV\int\frac{d^{3}p}{(2\pi)^{3}}\left(\frac{f-f^{q}}{1-q}-f\right),$
(2) $\displaystyle\text{Number: }N$
$\displaystyle=nV=gV\int\frac{d^{3}p}{(2\pi)^{3}}f^{q},$ (3)
$\displaystyle\text{Energy: }E$ $\displaystyle=\epsilon
V=gV\int\frac{d^{3}p}{(2\pi)^{3}}Ef^{q},$ (4) $\displaystyle\text{Pressure:
}P$ $\displaystyle=g\int\frac{d^{3}p}{(2\pi)^{3}}\frac{p^{2}}{3E}f^{q},$ (5)

where $g$ is the degeneracy factor, $V$ is the Tsallis volume; and $s$, $n$,
and $\epsilon$ correspond to entropy density, number density, and energy
density, respectively. Tsallis volume, $V$, should not be identified with the
volume determined from femtoscopy [25]. Additionally, we have also studied
thermodynamic response functions, such as heat capacity ($C_{V}$), isothermal
compressibility ($\kappa_{T}$), and expansion coefficient ($\alpha$), along
with mean free path ($\lambda$), Knudsen number ($Kn$), and squared speed of
sound ($c_{s}^{2}$) in the produced hadronic medium. The primary objective of
this work is to understand the interplay between collision energy and system
size on various thermodynamic quantities, contrasting their behavior in small
and large systems.

The remainder of this paper is organized as follows: In Section II, we discuss
the Tsallis distribution that is used to fit the transverse momentum spectra
of the charged hadrons. In Section III, we calculate the various thermodynamic
quantities and discuss their dependence on charged particle multiplicity.
Finally, in Section IV, we summarize our findings.

##  II Transverse momentum spectra

The $p_{T}$-spectra of different particle species can be accurately described
using the Tsallis distribution. In this framework, the invariant yield of a
particle is defined as:

$E\frac{d^{3}N}{d^{3}p}=\frac{gV}{(2\pi)^{3}}E\left[1+(q-1)\frac{E-\mu}{T}\right]^{\frac{-q}{q-1}}.$
(6)

This expression can be represented $p_{T}$, transverse mass
($m_{T}=\sqrt{p_{T}^{2}+m^{2}}$, where $m$ is the particle mass), and rapidity
($y$) as follows:

$\displaystyle\frac{d^{2}N}{p_{T}dp_{T}dy}=\frac{gV}{(2\pi)^{2}}m_{T}\cosh(y)$
$\displaystyle\
\times\left[1+(q-1)\frac{m_{T}\cosh(y)-\mu}{T}\right]^{\frac{-q}{q-1}}.$ (7)

Assuming $\mu=0$ at LHC energies, the expression for the invariant yield at
mid-rapidity ($y=0$) becomes

$\frac{d^{2}N}{dp_{T}dy}=\frac{gV}{(2\pi)^{2}}m_{T}p_{T}\left[1+(q-1)\frac{m_{T}}{T}\right]^{\frac{-q}{q-1}}.$
(8)

In the present study, since we are studying charged hadron spectra in various
collision systems, therefore, a transformation from rapidity ($y$) to
pseudorapidity ($\eta$) is necessary. This is achieved using the following
relation:

$\frac{dN}{dp_{T}d\eta}=\frac{dN}{dp_{T}dy}\left(\sqrt{1-\frac{m^{2}}{m_{T}^{2}\cosh^{2}(y)}}\right),$
(9)

which simplifies to:

$\frac{dN}{dp_{T}d\eta}=\frac{dN}{dp_{T}dy}\left(\frac{p_{T}}{m_{T}}\right),$
(10)

at mid-rapidity.

The charged hadron spectra from various collision systems were fitted with the
sum of Tsallis distributions for pions, kaons, and protons. The resulting
expression of the fit function is as follows [8]:

$\frac{d^{2}N}{dp_{T}d\eta}=2\frac{V}{(2\pi)^{2}}p_{T}^{2}\sum_{i=1}^{3}g_{i}\left[1+(q-1)\frac{m_{T,i}}{T}\right]^{\frac{-q}{q-1}},$
(11)

where $i=\pi^{+},K^{+},p$, and $g_{\pi}=g_{K}=1$ and $g_{p}=2$ are the
degeneracy factors. Factor of $2$ accounts for the antiparticles. This
notation is followed throughout the paper.

We have used Eq. (11) to fit the charged hadron spectra obtained in Pb+Pb
collisions at $\sqrt{s_{NN}}$ = 2.76 and 5.02 TeV, Xe+Xe collisions at
$\sqrt{s_{NN}}$ = 5.44 TeV, $p$+Pb collisions at $\sqrt{s_{NN}}$ = 5.02 TeV,
and $p$+$p$ collisions at $\sqrt{s}$ = 5.02 and 13 TeV, measured by the ALICE
Collaboration [26, 27, 28]. The extracted fit parameters are used to calculate
various thermodynamic variables. To achieve a fair comparison of these
thermodynamic quantities across different collision systems, we have fitted
the spectra measured within the same pseudorapidity window of $|\eta|<$ 0.8.
Furthermore, we have chosen to limit the fits to $p_{T}<$ 5 GeV/c to
investigate bulk properties, as high-$p_{T}$ particles are generated by hard
processes. In a previous study ([8]), various thermodynamic quantities were
calculated by fitting the spectra in the full $p_{T}$ range in Pb+Pb
collisions at $\sqrt{s_{NN}}$ = 2.76 and 5.02 TeV using Tsallis distribution.
It should be noted that the selection of the $p_{T}$ fitting range impacts the
fit parameters and consequently influences the resulting thermodynamic
variables [29].

Figure 1: Tsallis fits to the $p_{T}$-spectra of charged hadrons in different
centrality classes of Pb+Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV (left)
and 5.02 TeV (right). The lower panel of each plot shows the ratio of the data
and fit function.

Figure 1 shows the Tsallis fits to the charged hadron spectra in various
centrality classes of Pb+Pb collisions at $\sqrt{s_{NN}}$ = 2.76 and 5.02 TeV
using Eq. (11). It is evident that the Tsallis fits provide a more accurate
description of the charged hadron spectra in peripheral Pb+Pb collisions
compared to central collisions for both collision energies [8].

Figure 2: Tsallis fits to the $p_{T}$-spectra of charged hadrons in different
collision centralities of Xe+Xe collisions (left) and 0-100% centrality of
$p$+Pb collisions (right) at $\sqrt{s_{NN}}$ = 5.02 TeV. The lower panel of
each plot shows the ratio of the data and fit function.

Similary, Fig. 2 presents the Tsallis distribution fits to the charged hadron
spectra in different centrality classes of Xe+Xe collisions at $\sqrt{s_{NN}}$
= 5.44 TeV, and 0-100% centrality in $p$+Pb collisions at $\sqrt{s_{NN}}$ =
5.02 TeV. As was the case in Pb+Pb collisions, Tsallis fits provide a more
accurate description of the charged hadron spectra in peripheral Xe+Xe
collisions compared to central collisions. We also observe that the Tsallis
distribution provide a good description of the charged particle spectra in
$p$+Pb collisions.

Figure 3: Tsallis fits to the $p_{T}$-spectra of charged hadrons in different
multiplicity classes of $p$+$p$ collisions at $\sqrt{s_{NN}}$ = 5.02 TeV
(left) and 13 TeV (right). The lower panel of each plot shows the ratio of the
data and fit function.

Figure 3 shows Tsallis fits to the charged hadron spectra measured in
different multiplicity classes of $p$+$p$ collisions at $\sqrt{s}$ = 5.02 and
13 TeV. Notably, Tsallis fits provide a better description of the charged
hadron spectra in lower multiplicity classes in $p$+$p$ collisions compared to
higher multiplicity classes.

Figure 4: Tsallis parameters (volume, temperature, and $q$) as a function of
$\langle dN_{ch}/d\eta\rangle$, obtained by fitting the charged hadron
$p_{T}$-spectra in Pb+Pb collisions at $\sqrt{s_{NN}}$ = 2.76 and 5.02 TeV,
Xe+Xe collisions at $\sqrt{s_{NN}}$ = 5.44 TeV, $p$+Pb collisions at
$\sqrt{s_{NN}}$ = 5.02 TeV, and $p$+$p$ collisions at $\sqrt{s}$ = 5.02 and 13
TeV using Eq. (11).

Figure 4 shows the dependence of fit parameters ($V$, $T$, and $q$) on charged
particle multiplicity ($\langle dN_{ch}/d\eta\rangle$), which is an indicator
of the size of the hadronic system produced in high-energy collisions [27, 30,
28, 31]. The findings from Fig. 4 lead to the following conclusions:

  * •

The Tsallis volume, $V$, increases linearly with increasing $\langle
dN_{ch}/d\eta\rangle$ and hence is a good indicator of the system size.
However, this increase appears to be slower for $p$+$p$ collisions compared to
Pb+Pb and Xe+Xe collisions. Similar trend was also observed in femtoscopy
studies where small systems have smaller radii compared to large systems
extracted with two-pion and three-pion cumulants [32]. This suggests that at
the same $\langle dN_{ch}/d\eta\rangle$, the volume of the hadronic system
produced in $p$+$p$ collisions is smaller than in Pb+Pb and Xe+Xe collisions.
The single measurement in $p$+Pb collisions falls within the trend observed in
$p$+$p$ collisions at the corresponding $\langle dN_{ch}/d\eta\rangle$,
suggesting a potential convergence of behavior in small collision systems.

  * •

The temperature, $T$, increases with increasing $\langle
dN_{ch}/d\eta\rangle$. It is observed that $T$ rises sharply in Pb+Pb and
Xe+Xe collisions compared to $p$+$p$ collisions with increasing $\langle
dN_{ch}/d\eta\rangle$.

  * •

The entropy index, $q$, for $p$+$p$ collisions increases with increasing
$\langle dN_{ch}/d\eta\rangle$, suggesting a deviation from equilibrium in
high-multiplicity $p$+$p$ collisions. The $q$-parameter for $p$+Pb collisions
aligns with the trend observed in $p$+$p$ collisions at the corresponding
$\langle dN_{ch}/d\eta\rangle$. However, for Pb+Pb and Xe+Xe collisions,
$q$-parameter decreases with increasing $\langle dN_{ch}/d\eta\rangle$,
suggesting that the system approaches equilibrium more closely in central
collisions than in peripheral collisions. It is speculated that the initial
energy density induced hot spots caused by the Color Glass Condensate
formalism may generate significant temperature fluctuations especially in
small systems [33]. In high multiplicity $p$+$p$ collisions this will cause an
increase in the value of $q$, contrary to central heavy-ion collisions. We
also observe that $q$ increases with increasing center-of-mass energy in both
small and large collision systems, indicating a tendency for the system to
deviate further from equilibrium in collisions with higher energy.

The parameters obtained from the Tsallis fits are used to calculate various
thermodynamic quantities, which are explored in the following section.

##  III Thermodynamic Variables

###  III.1 Energy

The Tsallis parameters, $V$, $T$, and $q$ are used to compute the energy
density using the following equation:

$\epsilon=2\sum_{i=1}^{3}g_{i}\int\frac{d^{3}p}{(2\pi)^{3}}E_{i}\left(1+(q-1)\frac{E_{i}}{T}\right)^{\frac{-q}{q-1}},$
(12)

The left panel of Fig. 5 shows the energy density, $\epsilon$, and the right
panel shows the total energy, $E$ (= $\epsilon V$) of the hadronic medium at
kinetic freeze-out plotted as a function of $\langle dN_{ch}/d\eta\rangle$ in
various collision systems.

It is observed that the energy density and total energy increase with
increasing $\langle dN_{ch}/d\eta\rangle$. Furthermore, we find that the
energy density exhibits a steeper rise in small collision systems compared to
large collision systems, leading to a discontinuity at $\langle
dN_{ch}/d\eta\rangle$ $\approx$ 30\. This difference might arise because the
transverse overlap area in small systems is much smaller than that in large
systems for the same event multiplicities [33, 34]. This is further
corroborated by the fact that the discontinuity disappears when $E$ is plotted
against $\langle dN_{ch}/d\eta\rangle$ and a linear scaling between the two is
observed. This shows that high multiplicity $p$+$p$ collisions are capable of
producing energy densities that are comparable to the most peripheral heavy-
ion collisions [35].

Figure 5: Energy density (left panel) and total energy (right panel) at
kinetic freeze-out as a function of $\langle dN_{ch}/d\eta\rangle$ in Pb+Pb
collisions at $\sqrt{s_{NN}}$ = 2.76 and 5.02 TeV, Xe+Xe collisions at
$\sqrt{s_{NN}}$ = 5.44 TeV, $p$+Pb collisions at $\sqrt{s_{NN}}$ = 5.02 TeV,
and $p$+$p$ collisions at $\sqrt{s}$ = 5.02 and 13 TeV calculated using Eq
(12). The expression $a\langle dN_{ch}/d\eta\rangle+b$ is fitted to the total
energy and we find $a=1.463\pm 0.026$ GeV and $b=-1.163\pm 0.298$ GeV with
$\chi^{2}/ndf$ = 0.207.

###  III.2 Pressure

Figure 6: Pressure at kinetic freeze-out as a function of $\langle
dN_{ch}/d\eta\rangle$ in Pb+Pb collisions at $\sqrt{s_{NN}}$ = 2.76 and 5.02
TeV, Xe+Xe collisions at $\sqrt{s_{NN}}$ = 5.44 TeV, $p$+Pb collisions at
$\sqrt{s_{NN}}$ = 5.02 TeV, and $p$+$p$ collisions at $\sqrt{s}$ = 5.02 and 13
TeV calculated using Eq (13).  Figure 7: $P(\epsilon)$ equation of state of
charged hadrons in Pb+Pb collisions at $\sqrt{s_{NN}}$ = 2.76 and 5.02 TeV,
Xe+Xe collisions at $\sqrt{s_{NN}}$ = 5.44 TeV, $p$+Pb collisions at
$\sqrt{s_{NN}}$ = 5.02 TeV, and $p$+$p$ collisions at $\sqrt{s}$ = 5.02 and 13
TeV at kinetic freeze-out. The expression $c_{s}^{2}\epsilon$ is fitted to the
data and we find $c_{s}^{2}=0.262\pm 0.002$ with $\chi^{2}/ndf$ = 0.038.

Examining the pressure of hadrnoic medium produced in high energy collisions
can provide valuable insights into its equation of state. It can be evaluated
using the following equation:

$P=2\sum_{i=1}^{3}g_{i}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{p^{2}}{3E_{i}}\left(1+(q-1)\frac{E_{i}}{T}\right)^{\frac{-q}{q-1}}.$
(13)

Figure 6 shows the pressure of the hadronic medium at kinetic freeze-out as a
function of $\langle dN_{ch}/d\eta\rangle$, exhibiting an increase with
increasing $\langle dN_{ch}/d\eta\rangle$. It is also observed that the
pressure rises more rapidly in $p$+$p$ and $p$+Pb collisions than in Pb+Pb and
Xe+Xe collisions, resulting to a discontinuity at $\langle
dN_{ch}/d\eta\rangle$ $\approx$ 30\. This can possibly be attributed to the
higher initial densities in $p$+$p$ and $p$+Pb collisions, which may result in
a larger pressure exerted against the surrounding environment [33, 36, 37].

Within Landau’s hydrodynamic model framework, the equation of state for an
ultra-relativistic hadron gas is, $P=c_{s}^{2}\epsilon$, where $c_{s}^{2}$ is
the square of the speed of sound [38]. Figure 7 shows the linear correlation
between pressure and energy density. A linear fit to all the data points
results in $c_{s}^{2}$ = 0.262 $\pm$ 0.002. A more detailed discussion on
$c_{s}^{2}$ is presented later.

###  III.3 Particle number

Figure 8: Number density (left panel) and total number of particles (right
panel) at kinetic freeze-out as a function of $\langle dN_{ch}/d\eta\rangle$
in Pb+Pb collisions at $\sqrt{s_{NN}}$ = 2.76 and 5.02 TeV, Xe+Xe collisions
at $\sqrt{s_{NN}}$ = 5.44 TeV, $p$+Pb collisions at $\sqrt{s_{NN}}$ = 5.02
TeV, and $p$+$p$ collisions at $\sqrt{s}$ = 5.02 and 13 TeV calculated using
Eq (14). The expression $a\langle dN_{ch}/d\eta\rangle$ is fitted to the total
number of particles and we find $a=1.904\pm 0.030$ and $b=0.035\pm 0.448$ with
$\chi^{2}/ndf$ = 0.043.

The particle number density can be evaluated from the expression:

$n=2\sum_{i=1}^{3}g_{i}\int\frac{d^{3}p}{(2\pi)^{3}}\left(1+(q-1)\frac{E_{i}}{T}\right)^{\frac{-q}{q-1}},$
(14)

Figure 8 shows increase in the number density, $n$, and total number of
particles, $N$ (= $nV$), with increasing $\langle dN_{ch}/d\eta\rangle$. It is
observed that the increase in particle density is more pronounced for $p$+$p$
and $p$+Pb collisions than for Pb+Pb and Xe+Xe collisions, which could be
attributed to the smaller volume of the former collision systems. Moreover, a
linear increase of $N$ with $\langle dN_{ch}/d\eta\rangle$ further strengthens
the argument that particle number density is closely related to the volume of
the system.

###  III.4 Entropy

The entropy plays a major role in understanding the evolution of the system
produced in heavy-ion collisions. Entropy density can be evaluated using the
following expression:

$\displaystyle
s=2\sum_{i=1}^{3}g_{i}\int\frac{d^{3}p}{(2\pi)^{3}}\left[\frac{E_{i}}{T}\left(1+(q-1)\frac{E_{i}}{T}\right)^{\frac{-q}{q-1}}\right.$
$\displaystyle\left.+\left(1+(q-1)\frac{E_{i}}{T}\right)^{\frac{-1}{q-1}}\right].$
(15)

Figure 9: Entropy density (left panel) and total entropy (right panel) at
kinetic freeze-out as a function of $\langle dN_{ch}/d\eta\rangle$ in Pb+Pb
collisions at $\sqrt{s_{NN}}$ = 2.76 and 5.02 TeV, Xe+Xe collisions at
$\sqrt{s_{NN}}$ = 5.44 TeV, $p$+Pb collisions at $\sqrt{s_{NN}}$ = 5.02 TeV,
and $p$+$p$ collisions at $\sqrt{s}$ = 5.02 and 13 TeV calculated using Eq
(III.4). The expression $a\langle dN_{ch}/d\eta\rangle+b$ is fitted to total
entropy and we find $a=20.489\pm 0.368$ and $b=15.012\pm 4.939$ with
$\chi^{2}/ndf$ = 1.702.

In Fig. 9, we observe that the entropy density, $s$, and total entropy, $S$ (=
$sV$), increase with increasing $\langle dN_{ch}/d\eta\rangle$. This can be
attributed to an increase in both the particle number and the volume of the
system. Additionally, we note that the entropy density increases more rapidly
in $p$+$p$ and $p$+Pb collisions than in Pb+Pb and Xe+Xe collisions, leading
to a discontinuity at $\langle dN_{ch}/d\eta\rangle$ $\approx$ 30\. It can be
inferred from an almost linear relationship between $S$ and $\langle
dN_{ch}/d\eta\rangle$ that the dissimilarity in the rate of increase may be
ascribed to the volume of the hadronic system. We have also explicitly
verified that the thermodynamic relation, $\epsilon+P=Ts$, holds.

###  III.5 Mean free path and Knudsen number

Figure 10: Mean free path (left panel) and Knudsen number (right panel) at
kinetic freeze-out as a function of $\langle dN_{ch}/d\eta\rangle$ in Pb+Pb
collisions at $\sqrt{s_{NN}}$ = 2.76 and 5.02 TeV, Xe+Xe collisions at
$\sqrt{s_{NN}}$ = 5.44 TeV, $p$+Pb collisions at $\sqrt{s_{NN}}$ = 5.02 TeV,
and $p$+$p$ collisions at $\sqrt{s}$ = 5.02 and 13 TeV. Mean free path is
evaluated using Eq (16) and the definition of Knudsen number is given in the
text.

The mean free path ($\lambda$) of a system is defined as the average distance
that a particle travels between two successive collisions in the system.
Mathematically, the mean free path is expressed as [10]:

$\lambda=\frac{1}{n\sigma},$ (16)

where $n$ is the particle density, and $\sigma$ represents the scattering
cross-section. In the specific case of hard core hadron radius approximation,
with $r_{h}=0.3$ fm, $\sigma=4\pi r_{h}^{2}$ [10, 39, 40]. Furthermore,
Knudsen number ($Kn$) is a dimensionless parameter, defined as the ratio of
the mean free path ($\lambda$) to the radius, $R$ (= $(3V/4\pi)^{1/3}$), of
the hadronic system [14, 15]. A small value of $Kn$ signifies that the gas
system is in the continuum region, allowing the application of hydrodynamics.
Generally, when $Kn<0.5$, the fluid-dynamical criterion is satisfied,
indicating fluid-like behavior, whereas a value greater than unity indicates
that the particles are mostly free-streaming [41, 42, 43].

Figure 10, shows both $\lambda$ and $Kn$ decrease with increasing $\langle
dN_{ch}/d\eta\rangle$. The decrease in $\lambda$ results from the increase in
the particle number density with larger system sizes. Consequently, the
hadrons experience more frequent collisions, leading to a decrease in the
average distance between successive collisions. The decrease appears to be
more rapid in $p$+$p$ collisions because of the smaller volume of the hadronic
system compared to large collision systems at the same charged particle
multiplicity.

It was observed that, except for collisions with large $\langle
dN_{ch}/d\eta\rangle$, the Knudsen number ($Kn$) for the hadronic gas at KFO
was greater than 1. This suggests that the mean free path ($\lambda$) is
comparable to the system size ($R$), indicating that the hadronic system at
KFO is not in the continuum phase. Instead, the particles are predominantly
free-streaming.

###  III.6 Heat capacity

The heat capacity ($C_{V}$) of a system quantifies the amount of energy
required to increase its temperature by one unit. Mathematically, this can be
expressed as [10, 11]:

$C_{V}=\left.\frac{\partial\epsilon}{\partial T}\right|_{V}$ (17)

Figure 11 shows the variation in $C_{V}$ as a function of $\langle
dN_{ch}/d\eta\rangle$. We observe an increase in $C_{V}$ with increasing
$\langle dN_{ch}/d\eta\rangle$. This could be a consequence of the increase in
the particle number with $\langle dN_{ch}/d\eta\rangle$. As the number of
particles in a system increases, a higher amount of energy is required to
increase its temperature by one unit.

We also observe a steeper increase in $C_{V}$ with $\langle
dN_{ch}/d\eta\rangle$ in the case of $p$+$p$ collisions, leading to a
discontinuity at $\langle dN_{ch}/d\eta\rangle$ $\approx$ 30\. This might be
due to the higher value of the $q$-parameter in $p$+$p$ collisions. Further
the system is away from equilibrium, as indicated by a higher $q$-value, more
is the change in its internal energy as temperature rises, resulting in a
higher $C_{V}$ [10].

Figure 11: Heat capacity at kinetic freeze-out as a function of $\langle
dN_{ch}/d\eta\rangle$ in Pb+Pb collisions at $\sqrt{s_{NN}}$ = 2.76 and 5.02
TeV, Xe+Xe collisions at $\sqrt{s_{NN}}$ = 5.44 TeV, $p$+Pb collisions at
$\sqrt{s_{NN}}$ = 5.02 TeV, and $p$+$p$ collisions at $\sqrt{s}$ = 5.02 and 13
TeV calculated using Eq (17).

###  III.7 Isothermal compressiblity

Figure 12: Isothermal compressibility at kinetic freeze-out as a function of
$\langle dN_{ch}/d\eta\rangle$ in Pb+Pb collisions at $\sqrt{s_{NN}}$ = 2.76
and 5.02 TeV, Xe+Xe collisions at $\sqrt{s_{NN}}$ = 5.44 TeV, $p$+Pb
collisions at $\sqrt{s_{NN}}$ = 5.02 TeV, and $p$+$p$ collisions at $\sqrt{s}$
= 5.02 and 13 TeV calculated using Eq. (18).

The isothermal compressibility ($\kappa_{T}$) of a system is a measure of the
extent to which its volume changes in response to an external pressure. For an
ideal fluid that is both incompressible and non-viscous, $\kappa_{T}$ is $0$.
It can be evaluated from the expression [10, 12, 13]:

$\kappa_{T}=\frac{\partial n/\partial\mu}{n^{2}}.$ (18)

Figure 12 shows isothermal compressibility as a function of $\langle
dN_{ch}/d\eta\rangle$. We observe that $\kappa_{T}$ decreases with increasing
$\langle dN_{ch}/d\eta\rangle$, suggesting that the system approaches near
ideal behavior with increasing multiplicity. This can be explained by the fact
that higher multiplicity classes consist of a greater number of particles,
thus requiring higher pressure to achieve a small change in volume [12].

###  III.8 Expansion coefficient

Figure 13: Expansion coefficient at kinetic freeze-out as a function of
$\langle dN_{ch}/d\eta\rangle$ in Pb+Pb collisions at $\sqrt{s_{NN}}$ = 2.76
and 5.02 TeV, Xe+Xe collisions at $\sqrt{s_{NN}}$ = 5.44 TeV, $p$+Pb
collisions at $\sqrt{s_{NN}}$ = 5.02 TeV, and $p$+$p$ collisions at $\sqrt{s}$
= 5.02 and 13 TeV calculated using Eq. (19).

The expansion coefficient or isobaric expansivity ($\alpha$) of a system is
the measure of relative expansion per degree change in temperature at a
constant pressure. For the hadronic medium at KFO, $\alpha$ is calculated
using the following expression [10]:

$\alpha=-\frac{\partial n/\partial T}{n}.$ (19)

Figure 13 shows $\alpha$ as a function $\langle dN_{ch}/d\eta\rangle$. The
value of $\alpha$ is observed to be negative in the studied collision systems.
$\alpha$ increases with increasing $\langle dN_{ch}/d\eta\rangle$. This is a
consequence of the increasing particle density with $\langle
dN_{ch}/d\eta\rangle$. For a denser system, the increase in temperature
produces a small change in its volume, whereas for a less dense system, the
change in volume is comparatively higher.

###  III.9 Speed of sound

Figure 14: Squared speed of sound at kinetic freeze-out as a function of
$\langle dN_{ch}/d\eta\rangle$ in Pb+Pb collisions at $\sqrt{s_{NN}}$ = 2.76
and 5.02 TeV, Xe+Xe collisions at $\sqrt{s_{NN}}$ = 5.44 TeV, $p$+Pb
collisions at $\sqrt{s_{NN}}$ = 5.02 TeV, and $p$+$p$ collisions at $\sqrt{s}$
= 5.02 and 13 TeV calculated using Eq. (20). The dotted line represents the
value of $c_{s}^{2}$ for a non-interacting massless ideal gas.

The speed of sound in hadronic medium is closely related to its equation of
state. The squared speed of sound is given as [11, 12]:

$c_{s}^{2}=\left(\frac{\partial
P}{\partial\epsilon}\right)_{s}=\left(\frac{\partial P/\partial
T}{\partial\epsilon/\partial T}\right).$ (20)

For a non-interacting massless ideal gas, the expected value of $c_{s}^{2}$ is
1/3 times the speed of light squared [44]. Figure 14 shows $c_{s}^{2}$ plotted
as a function of $\langle dN_{ch}/d\eta\rangle$. $c_{s}^{2}$ shows a weak
dependence with $\langle dN_{ch}/d\eta\rangle$. We observe that for studied
collision systems and energies, on an average $c_{s}^{2}$ = 0.266 $\pm$ 0.001.

##  IV Summary

The charged hadron spectra measured by the ALICE Collaboration in Pb+Pb
collisions at $\sqrt{s_{NN}}$ = 2.76 and 5.02 TeV, Xe+Xe collisions at
$\sqrt{s_{NN}}$ = 5.44 TeV, $p$+Pb collisions at $\sqrt{s_{NN}}$ = 5.02 TeV,
and $p$+$p$ collisions at $\sqrt{s}$ = 5.02 and 13 TeV have been studied in
the framework of Tsallis statistics. We find that the Tsallis fits provide a
better description of charged hadron spectra in peripheral Pb+Pb and Xe+Xe
collisions than in more central collisions. In addition, we also observe an
improvement in the fitting for lower multiplicity classes in $p$+$p$
collisions and a reasonably good description of the spectra in $p$+Pb
collisions. Furthermore, the dependence of the fit parameters ($V$, $T$, and
$q$) on the $\langle dN_{ch}/d\eta\rangle$ was also noted. It is observed
that:

  * •

Tsallis volume of the hadronic medium increases with increasing $\langle
dN_{ch}/d\eta\rangle$, but the growth rate of volume in $p$+$p$ collisions is
slower compared to Pb+Pb and Xe+Xe collisions, while $p$+Pb collisions show a
similar trend to $p$+$p$ collisions. This suggests that small collision
systems produce smaller hadronic systems at KFO for the same $\langle
dN_{ch}/d\eta\rangle$ compared to large collision systems.

  * •

Temperature increases with increasing $\langle dN_{ch}/d\eta\rangle$ with a
slow increase in small collision systems and a sharp rise in large collision
systems.

  * •

The entropy index, $q$, increases with increasing $\langle
dN_{ch}/d\eta\rangle$ in $p+p$ collisions, suggesting a deviation from
equilibrium in high multiplicity collisions. $q$-parameter in $p$+Pb
collisions is observed to follow the same trend as $p+p$ collisions at the
corresponding $\langle dN_{ch}/d\eta\rangle$. However, for Pb+Pb and Xe+Xe
collisions, $q$ decreases with increasing $\langle dN_{ch}/d\eta\rangle$,
indicating that the system is closer to equilibrium in central collisions.

The thermodynamic properties of the hadronic medium in high-energy collisions,
focusing on the dependence of energy density, pressure, particle density,
entropy density, mean free path, Knudsen number, heat capacity, isothermal
compressibility, expansion coefficient, and speed of sound on $\langle
dN_{ch}/d\eta\rangle$ are studied. We find that:

  * •

Energy density increases with increasing $\langle dN_{ch}/d\eta\rangle$, with
small collision systems having a higher energy density in the common $\langle
dN_{ch}/d\eta\rangle$ region, possibly because of the smaller volumes of the
hadronic medium compared to large collision systems.

  * •

Pressure increases with increasing $\langle dN_{ch}/d\eta\rangle$, with small
collision systems exhibiting a faster increase. This might be explained by the
fact that less space for particles to move around in small collision systems
results into a more rapid momentum exchange.

  * •

Number density increases with increasing $\langle dN_{ch}/d\eta\rangle$, with
smaller systems exhibiting a faster increase owing to the smaller volume of
the produced hadronic system.

  * •

Entropy density increases with increasing $\langle dN_{ch}/d\eta\rangle$, with
smaller systems having a higher rate of increase.

  * •

The mean free path decreases with increasing $\langle dN_{ch}/d\eta\rangle$
due to the increased particle number density. Knudsen number was found to be
much greater than unity, suggesting that produced hadronic medium is in the
free streaming region.

  * •

Heat capacity increases with increasing $\langle dN_{ch}/d\eta\rangle$. This
is a consequence of increase in particle number with increasing $\langle
dN_{ch}/d\eta\rangle$. Increasing the temperature of a system by one unit will
require more energy when the system contains a greater number of particles.

  * •

Isothermal compressibility decreases with increasing $\langle
dN_{ch}/d\eta\rangle$ indicating near-ideal behavior with increasing
multiplicity.

  * •

Expansion coefficient is negative for the studied collision systems and it
increases with increasing $\langle dN_{ch}/d\eta\rangle$. This is due to the
higher particle density associated with greater $\langle
dN_{ch}/d\eta\rangle$, resulting in a smaller change in volume for the same
amount of pressure in a relatively denser system.

  * •

Squared speed of sound shows a weak dependence with $\langle
dN_{ch}/d\eta\rangle$. We observe an average value of $c_{s}^{2}$ = 0.266
$\pm$ 0.001 for the studied collisions systems and energies.

This study characterizes the thermodynamic properties of the produced hadronic
system in different collisions across a range of LHC energies. We observe that
high multiplicity $p$+$p$ collisions produce a hadronic medium with
thermodynamic properties similar to peripheral heavy-ion collisions. These
observations may contribute to the efforts to understand the properties of the
medium formed in high-multiplicity $p$+$p$ collisions.

###### Acknowledgements.

CJ acknowledges the financial support from DAE-DST Project No.
3015/I/2021/Gen/RD-I/13283.

## References

  * Ejiri [2008] S. Ejiri, Phys. Rev. D 78, 074507 (2008).
  * Aoki  _et al._ [2006] Y. Aoki, G. Endrődi, Z. Fodor, S. D. Katz, and K. K. Szabó, Nature 443, 675 (2006).
  * Adams  _et al._ [2005] J. Adams _et al._ (STAR Collaboration), Nucl. Phys. A 757, 102 (2005).
  * Adam  _et al._ [2017] J. Adam _et al._ (ALICE Collaboration), Nature Phys. 13, 535 (2017).
  * Bjorken  _et al._ [2013] J. D. Bjorken, S. J. Brodsky, and A. S. Goldhaber, Phys. Lett. B 726, 344 (2013).
  * Khachatryan  _et al._ [2017] V. Khachatryan _et al._ (CMS Collaboration), Phys. Lett. B 765, 193 (2017).
  * Khuntia  _et al._ [2019a] A. Khuntia, H. Sharma, S. K. Tiwari, R. Sahoo, and J. Cleymans, Eur. Phys. J. A 55 (2019a).
  * Azmi  _et al._ [2020] M. D. Azmi, T. Bhattacharyya, J. Cleymans, and M. Paradza, J. Phys. G: Nucl. Part. Phys 47, 045001 (2020).
  * Deb  _et al._ [2021] S. Deb, G. Sarwar, R. Sahoo, and J. e Alam, Eur. Phys. J. A 57 (2021).
  * Sahu and Sahoo [2021a] D. Sahu and R. Sahoo, Physics 3, 207 (2021a).
  * Khuntia  _et al._ [2016] A. Khuntia, P. Sahoo, P. Garg, R. Sahoo, and J. Cleymans, Eur. Phys. J. A 52 (2016).
  * Jain  _et al._ [2023] S. Jain, R. Gupta, and S. Jena, Universe 9, 170 (2023).
  * Khuntia  _et al._ [2019b] A. Khuntia, S. K. Tiwari, P. Sharma, R. Sahoo, and T. K. Nayak, Phys. Rev. C 100, 014910 (2019b).
  * Sarkar and Ghosh [2017] N. Sarkar and P. Ghosh, Phys. Rev. C 96, 044901 (2017).
  * Sarkar  _et al._ [2019] N. Sarkar, P. Deb, and P. Ghosh, arxiv:1905.06532v3 (2019).
  * Tsallis [1988] C. Tsallis, J. Stat. Phys. 52, 479 (1988).
  * Cleymans and Worku [2012a] J. Cleymans and D. Worku, J. Phys. G: Nucl. Part. Phys. 39, 025006 (2012a).
  * Hui  _et al._ [2018] J.-Q. Hui, Z.-J. Jiang, and D.-F. Xu, Adv. High Energy Phys. 2018, 7682325 (2018).
  * Bíró  _et al._ [2017] G. Bíró, G. Barnaföldi, T. Biró, K. Ürmössy, and Á. Takács, Entropy 19, 88 (2017).
  * Si  _et al._ [2018] R.-F. Si, H.-L. Li, and F.-H. Liu, Adv. High Energy Phys. 2018, 7895967 (2018).
  * Lao  _et al._ [2017] H.-L. Lao, F.-H. Liu, and R. A. Lacey, Eur. Phys. J. A 53 (2017).
  * Cleymans and Worku [2012b] J. Cleymans and D. Worku, Eur. Phys. J. A 48 (2012b).
  * Bhattacharyya  _et al._ [2016] T. Bhattacharyya, J. Cleymans, and S. Mogliacci, Phys. Rev. D 94, 094026 (2016).
  * Bhattacharyya  _et al._ [2017] T. Bhattacharyya, J. Cleymans, P. Garg, P. Kumar, S. Mogliacci, R. Sahoo, and S. Tripathy, J. Phys.: Conf. Ser. 878, 012016 (2017).
  * Cleymans and Paradza [2020] J. Cleymans and M. W. Paradza, Physics 2, 654 (2020).
  * Acharya  _et al._ [2018] S. Acharya _et al._ (ALICE Collaboration), J. High Energ. Phys. 2018, 13 (2018).
  * Acharya  _et al._ [2019a] S. Acharya _et al._ (ALICE Collaboration), Phys. Lett. B 788, 166 (2019a).
  * Acharya  _et al._ [2019b] S. Acharya _et al._ (ALICE Collaboration), Eur. Phys. J. C 79 (2019b).
  * Patra  _et al._ [2021] R. N. Patra, B. Mohanty, and T. K. Nayak, Eur. Phys. J. Plus 136 (2021).
  * Abelev  _et al._ [2013] B. Abelev _et al._ (ALICE Collaboration), Phys. Rev. C 88, 044910 (2013).
  * Acharya  _et al._ [2020] S. Acharya _et al._ (ALICE Collaboration), Phys. Lett. B 800, 135043 (2020).
  * Grosse-Oetringhaus [2014] J. F. Grosse-Oetringhaus, Nuclear Physics A 931, 22 (2014), qUARK MATTER 2014.
  * Liu  _et al._ [2023] L. Liu, Z.-B. Yin, and L. Zheng, Chinese Phys. C 47, 024103 (2023).
  * Loizides  _et al._ [2018] C. Loizides, J. Kamin, and D. d’Enterria, Phys. Rev. C 97, 054910 (2018).
  * Sahu and Sahoo [2021b] D. Sahu and R. Sahoo, J. Phys. G: Nucl. Part. Phys. 48, 125104 (2021b).
  * Shuryak and Zahed [2013] E. Shuryak and I. Zahed, Phys. Rev. C 88, 044915 (2013).
  * Castorina  _et al._ [2020] P. Castorina, A. Iorio, D. Lanteri, H. Satz, and M. Spousta, Phys. Rev. C 101, 054902 (2020).
  * Mohanty and Alam [2003] B. Mohanty and J.-e. Alam, Phys. Rev. C 68, 064903 (2003).
  * Kadam and Mishra [2015] G. P. Kadam and H. Mishra, Phys. Rev. C 92, 035203 (2015).
  * Bugaev  _et al._ [2013] K. A. Bugaev, D. R. Oliinychenko, A. S. Sorin, and G. M. Zinovjev, Eur. Phys. J. A 49 (2013).
  * Huovinen and Molnar [2009] P. Huovinen and D. Molnar, Phys. Rev. C 79, 014906 (2009).
  * Bouras  _et al._ [2010] I. Bouras, E. Molnár, H. Niemi, Z. Xu, A. El, O. Fochler, C. Greiner, and D. H. Rischke, Phys. Rev. C 82, 024910 (2010).
  * Gombeaud  _et al._ [2009] C. Gombeaud, T. Lappi, and J.-Y. Ollitrault, Phys. Rev. C 79, 054914 (2009).
  * Ferrer and Hackebill [2023] E. Ferrer and A. Hackebill, Nuclear Physics A 1031, 122608 (2023).