# Hadron momentum spectra from analytical solutions of relativistic
hydrodynamics

Mahammad Sabir Ali  School of Physical Sciences, National Institute of Science
Education and Research, An OCC of Homi Bhabha National Institute,
Jatni-752050, India  Deeptak Biswas  School of Physical Sciences, National
Institute of Science Education and Research, An OCC of Homi Bhabha National
Institute, Jatni-752050, India  Amaresh Jaiswal  School of Physical Sciences,
National Institute of Science Education and Research, An OCC of Homi Bhabha
National Institute, Jatni-752050, India  Institute of Theoretical Physics,
Jagiellonian University, ul. St. Łojasiewicza 11, 30-348 Krakow, Poland
Sushant K. Singh  Department of Physics & Astronomy, University of Florence,
Via G. Sansone 1, 50019 Sesto Fiorentino, Florence, Italy  Variable Energy
Cyclotron Centre, 1/AF, Bidhan Nagar , Kolkata 700064, India

###### Abstract

We present analytical solution of relativistic hydrodynamics for a system
having cylindrical symmetry with boost-invariant longitudinal expansion and
Hubble-like transverse expansion. We also consider analytical solution for
Hubble-like spherically expanding system. For these two cases, we calculate
analytical expression for transverse momentum spectra of hadrons, at constant
temperature freeze-out hypersurface using Cooper-Frye prescription. We compare
our results for transverse momentum spectra with experimental results from
Large Hadron Collider and CERN SPS where one expects cylindrical and spherical
geometry of the fireball, respectively. In the case of low-energy collisions
with spherical geometry, we calculate rapidity spectra and compare with the
results from CERN SPS.

##  I Introduction

The earliest hydrodynamic approach for describing nucleus-nucleus collisions
dates back to Landau’s work in 1953, where he focused on longitudinal
expansion along the collision axis [1]. Landau’s non-dissipative expansion
model provided an approximate analytical solution to ideal hydrodynamic
equations, known as Landau hydrodynamics, resulting in Gaussian-like rapidity
distributions of produced particles [1, 2]. Subsequent developments in
hydrodynamics included the Hwa-Bjorken framework, proposing a boost-invariant
scenario with a plateau-like rapidity distribution [3, 4]. However, the
applicability of boost-invariance symmetry was limited to the mid-rapidity
region in ultra-relativistic heavy-ion collisions. In contrast, the
transverse-momentum ($p_{T}$) integrated yield across the entire rapidity
region demonstrated an overall better agreement with the Gaussian structure
suggested by Landau [5, 6, 7, 8, 9, 2, 10]. Since Landau and Hwa-Bjorken
solutions, several different forms of analytical solutions were proposed
motivated by anticipated collective behaviour in relativistic heavy ion
collisions [11, 12, 13, 14, 15, 16]. Analytical solutions of relativistic
hydrodynamics in simple symmetric cases are important because they serve as
benchmarks for testing more realistic hydrodynamic simulation codes. In the
last two decades, relativistic hydrodynamic simulations have been successfully
applied to model the space-time evolution of the hot and dense matter produced
at the Relativistic Heavy Ion Collider (RHIC) at BNL and the Large Hadron
Collider (LHC) at CERN [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29,
30, 31].

Majority of studies in hydrodynamic modelling of the evolution of hot and
dense QCD matter, formed in relativistic heavy-ion collisions, rely on
numerical simulations. The hydrodynamic evolution equations, comprising of
energy-momentum and current conservation equations, are solved numerically
with appropriate initial conditions. The conversion from hydrodynamic fields
to the particle degree of freedom is achieved via the Cooper-Frye freezeout
prescription [32], which leads to the momentum distribution of observed
hadrons. Semi-analytical framework of hydrodynamics-inspired phenomenological
models provides another approach to study hadron momentum distribution [33,
34, 35, 36, 37, 38]. Due to their simplicity, these phenomenological models
have been employed extensively to analyze the momentum spectra of produced
hadrons, offering insights into the properties of the strongly interacting
matter at kinetic freeze-out [39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
51, 52, 53, 54, 55, 56, 57]. These phenomenological models, like the Blast-
wave model, primarily rely on the Cooper-Frye freezeout prescription with ad-
hoc parameterization of fluid variables at the freeze-out hypersurface. While
the parameterization of the fluid variables are inspired by hydrodynamic
solutions, a direct correspondence between them has not yet been established.
Therefore it would be highly beneficial to have a simple and fully analytical
hydrodynamic framework for computation of hadron momentum distribution.

In this work, we obtain an analytical solution for relativistic hydrodynamics
in a system characterized by cylindrical symmetry with boost-invariant
longitudinal expansion and Hubble-like transverse expansion. Additionally, we
explore an analytical solution for a Hubble-like spherically expanding system.
These straightforward profiles allow us to derive a simple analytical solution
for the temperature evolution relative to the proper time, facilitating a
connection with the freeze-out hypersurface. For both cylindrical and
spherical flow profiles, we compute the analytical expressions for the
transverse momentum spectra of hadrons on a constant temperature freeze-out
hypersurface, utilizing the Cooper-Frye prescription. We then compare these
results for transverse momentum spectra with experimental findings from the
Large Hadron Collider (LHC) and CERN SPS, where the fireball is expected to
exhibit cylindrical and spherical geometries, respectively. In the case of low
energy collisions with spherical expansion, we further calculate rapidity
spectra and compare them with results from CERN SPS. The paper is organized as
follows, in Sec. II we provide the derivation of the analytical hydrodynamic
solutions for spherical and cylindrical geometries. Subsequently, in Sec. III
we formulate the hadron momentum spectra using Cooper-Frye prescription for
spherical and cylindrical freezeout hypersurface. Further, in Sec. IV, we
apply the derived expressions for hadron momentum distribution to compare with
experimental results for hadron transverse momentum spectrum and rapidity
spectrum. Finally in Sec. V, we provide the summary and conclusion of the
present work. Throughout the text we employ natural units where
$\hbar=c=k_{B}=1$ and use $(+,-,-,-)$ for the metric signature. We use dot to
denote scalar products of both three- and four-vectors, i.e., $a\cdot
b=a^{0}b^{0}-\vec{a}\cdot\vec{b}$.

##  II Analytic solution of hydrodynamics

Hydrodynamic evolution of the hot and dense QCD matter is governed by the
conservation equations for particle charge current as well as energy and
momentum conservation. Assuming a system with no conserved charges, the
hydrodynamical equations are

$\partial_{\mu}T^{\mu\nu}(x)=0,$ (1)

where $T^{\mu\nu}(x)$ is the local energy-momentum tensor. In the following,
we consider non-dissipative evolution of the system where the energy-momentum
tensor takes the form

$T^{\mu\nu}=(\varepsilon+p)u^{\mu}u^{\nu}-pg^{\mu\nu},$ (2)

where $\varepsilon$ and $p$ denote the energy density and pressure of the
system, respectively. The local four-velocity is denoted by
$u^{\mu}(x)=\gamma(x)(1,\vec{v}(x))$, with
$\gamma(x)=1/\sqrt{1-\vec{v}^{2}(x)}$ being the local Lorentz factor. We use
the Minkowski metric for our calculations, given by $g^{\mu\nu}={\rm
diag}(1,-1,-1,-1)$.

Using the thermodynamic expressions $Ts=\varepsilon+p$ and $dp=s\,dT$, we can
express Eq. (1) in terms of temperature $T$ as

$u^{\mu}\partial_{\mu}(Tu^{\nu})-\partial^{\nu}T=0,$ (3)

where $s$ is the entropy density. The above equation can be re-written in the
following form

$\frac{\partial}{\partial
t}(T\gamma\vec{v})+\vec{\nabla}(T\gamma)=\vec{v}\times\vec{\nabla}\times(T\gamma\vec{v}).$
(4)

For non-dissipative fluid, the entropy conservation implies

$\partial_{\mu}(su^{\mu})=0,$ (5)

which can be re-written in the following form

$\frac{\partial}{\partial t}(s\gamma)+\vec{\nabla}\cdot(s\gamma\vec{v})=0.$
(6)

In the following, we will seek solutions to Eqs. (4) and (6) for a
hydrodynamically expanding system with spherical and cylindrical symmetry.

###  II.1 Spherical solution

In the case of radially expanding system, it is convenient to work in the
$(t,r,\theta,\phi)$ co-ordinates. Due to spherical symmetry, the hydrodynamic
variables depend only on time $t$ and radial distance $r$ from the origin. In
this case, the fluid velocity vector is given by $\vec{v}=(v_{r},0,0)$. For
such a system, Eqs. (4) and (6) become

$\displaystyle\frac{\partial}{\partial t}(T\gamma
v_{r})+\frac{\partial}{\partial r}(T\gamma)$ $\displaystyle=0$ (7)
$\displaystyle\frac{\partial}{\partial
t}(s\gamma)+\frac{1}{r^{2}}\frac{\partial}{\partial r}(r^{2}s\gamma v_{r})$
$\displaystyle=0.$ (8)

Above equations admit scaling type analytical solution with the fluid flow
velocity of the form $v_{r}=r/t$ [58, 59]. Assuming this flow profile, Eqs.
(7) and (8), become

$\displaystyle r\frac{\partial T}{\partial t}+t\frac{\partial T}{\partial r}$
$\displaystyle=0$ (9) $\displaystyle t\frac{\partial s}{\partial
t}+r\frac{\partial s}{\partial r}+3s$ $\displaystyle=0.$ (10)

We now require the equation of state relating energy density and pressure. We
assume an equation of state of the form $p=c_{s}^{2}\,\varepsilon$, with the
speed of sound $c_{s}$ kept constant.

By employing thermodynamic relations, we reformulate the equation of state to
relate entropy density and temperature as $s\propto T^{1/c_{s}^{2}}$.
Substituting this expression for $s$ in Eq. (10), we obtain

$t\frac{\partial T}{\partial t}+r\frac{\partial T}{\partial
r}+3\,c_{s}^{2}\,T=0.$ (11)

The above equation can be solved using the method of characteristics. The
equation for obtaining the characteristic curves is given by

$\frac{dt}{t}=\frac{dr}{r}=-\frac{dT}{3\,c_{s}^{2}\,T}.$ (12)

Solving the above equation, we obtain

$T(r,t)\sim t^{-3c_{s}^{2}}\,f\\!\left(\frac{r}{t}\right),$ (13)

The above equation implies that any functional form for $f(r/t)$ is a solution
of Eq. (11). This functional form can be fixed by employing Eq. (9).
Substituting the expression for $T(r,t)$ from Eq. (13) in Eq. (9), we get

$f^{\prime}(w)(1-w^{2})-3c_{s}^{2}wf(w)=0,$ (14)

where we have defined $w\equiv r/t$. The above equation can be easily
integrated to obtain

$f(w)\sim\frac{1}{(1-w^{2})^{3c_{s}^{2}/2}}.$ (15)

Using the above equation in the expression for $T(r,t)$ and $s(r,t)$, we
obtain the solution to Eqs. (7) and (8) as

$\displaystyle T(r,t)$
$\displaystyle=T_{0}\left[\frac{t_{0}^{2}-r_{0}^{2}}{t^{2}-r^{2}}\right]^{3c_{s}^{2}/2},$
(16) $\displaystyle s(r,t)$
$\displaystyle=s_{0}\left[\frac{t_{0}^{2}-r_{0}^{2}}{t^{2}-r^{2}}\right]^{3/2},$
(17) $\displaystyle v_{r}$ $\displaystyle=\frac{r}{t},$ (18)

where $T_{0}$ and $s_{0}$ are the initial temperature and entropy density,
respectively. We note that the scaling solution, given in Eqs. (16)-(18), is
in agreement with the family of solutions obtained in Ref. [12].

###  II.2 Cylindrical solution

In the case of a longitudinally and transversely expanding system, it is
convenient to work in the $(t,\rho,\phi,z)$ coordinates. Due to cylindrical
symmetry, the hydrodynamic variables depend only on time $t$, radial distance
$\rho$ from the origin in the transverse plane, and the longitudinal distance
$z$. In this case, the fluid velocity vector is given by
$\vec{v}=(v_{\rho},0,v_{z})$. For such a system, Eqs. (4) and (6) leads to

$\displaystyle\frac{\partial}{\partial t}(T\gamma
v_{\rho})+\frac{\partial}{\partial\rho}(T\gamma)-v_{z}\left[\frac{\partial}{\partial\rho}(T\gamma
v_{z})-\frac{\partial}{\partial z}(T\gamma v_{\rho})\right]$ $\displaystyle=0$
(19) $\displaystyle\frac{\partial}{\partial t}(T\gamma
v_{z})+\frac{\partial}{\partial
z}(T\gamma)-v_{\rho}\left[\frac{\partial}{\partial z}(T\gamma
v_{\rho})-\frac{\partial}{\partial\rho}(T\gamma v_{z})\right]$
$\displaystyle=0$ (20) $\displaystyle\frac{\partial}{\partial
t}(s\gamma)+\frac{1}{\rho}\frac{\partial}{\partial\rho}(\rho s\gamma
v_{\rho})+\frac{\partial}{\partial z}(s\gamma v_{z})$ $\displaystyle=0.$ (21)

In the following, we consider boost-invariance for the longitudinal expansion,
i.e., $v_{z}=\frac{z}{t}$. Bjorken’s notion of boost-invariance implies that
the expansion geometry of the fireball is identical in all longitudinally
boosted frames. Therefore, considering the $z=0$ slice, Eqs. (19)-(21) reduces
to

$\displaystyle\frac{\partial}{\partial t}(T\gamma
v_{\rho})+\frac{\partial}{\partial\rho}(T\gamma)$ $\displaystyle=0,$ (22)
$\displaystyle\frac{\partial}{\partial
t}(s\gamma)+\frac{\partial}{\partial\rho}(s\gamma
v_{\rho})+s\gamma\left(\frac{v_{\rho}}{\rho}+\frac{1}{t}\right)$
$\displaystyle=0.$ (23)

As done in the previous case of spherical expansion, we again assume an
equation of state of the form $p=c_{s}^{2}\,\varepsilon$ with constant
$c_{s}^{2}$.

In order to find the solution of Eqs. (22)-(23), we introduce the potential
$\Phi(T)$ and define $a_{\pm}$ as [59]

$\displaystyle d\Phi$ $\displaystyle=\frac{d\ln T}{c_{s}}=c_{s}\ d\ln s,$ (24)
$\displaystyle a_{\pm}$ $\displaystyle=e^{\Phi\pm\alpha},$ (25)

where $\alpha$ denotes the transverse rapidity, i.e.,
$\alpha=\tanh^{-1}(v_{\rho})$. The above equations can be inverted to write
the hydrodynamic variables in terms of these new variables $a_{\pm}$ as

$\displaystyle\Phi(T)$ $\displaystyle=\frac{1}{2}\ln a_{+}a_{-},$ (26)
$\displaystyle v_{\rho}$ $\displaystyle=\frac{a_{+}-a_{-}}{a_{+}+a_{-}}.$ (27)

In terms of the new variables $a_{\pm}$, Eqs. (22) and (23) can be rewritten
as [59]

$\displaystyle\frac{\partial a_{+}}{\partial
t}+\frac{v_{\rho}+c_{s}}{1+v_{\rho}\,c_{s}}\,\frac{\partial
a_{+}}{\partial\rho}+\frac{c_{s}a_{+}}{1+v_{\rho}\,c_{s}}\left(\frac{v_{\rho}}{\rho}+\frac{1}{t}\right)$
$\displaystyle=0,$ (28) $\displaystyle\frac{\partial a_{-}}{\partial
t}+\frac{v_{\rho}-c_{s}}{1-v_{\rho}\,c_{s}}\,\frac{\partial
a_{-}}{\partial\rho}+\frac{c_{s}a_{-}}{1-v_{\rho}\,c_{s}}\left(\frac{v_{\rho}}{\rho}+\frac{1}{t}\right)$
$\displaystyle=0.$ (29)

We intend to find a self-similar solution for the above set of equations.

We start with a simple scaling ansatz for the transverse component of the
fluid velocity, $v_{\rho}=\rho/t$. Using this ansatz in Eq. (27), the
variables $a_{\pm}$ become

$\displaystyle a_{+}(\rho,t)$
$\displaystyle=A(\rho,t)\left(1+\frac{\rho}{t}\right),$ (30) $\displaystyle
a_{-}(\rho,t)$ $\displaystyle=A(\rho,t)\left(1-\frac{\rho}{t}\right),$ (31)

where $A(\rho,t)$ is a constant of proportionality. Substituting the above
expressions in Eqs. (28) and (29), we get

$\displaystyle t(t+\rho c_{s})\frac{\partial A}{\partial
t}+t(\rho+tc_{s})\frac{\partial A}{\partial\rho}+Ac_{s}(3t-\rho)$
$\displaystyle=0$ (32) $\displaystyle t(t-\rho c_{s})\frac{\partial
A}{\partial t}+t(\rho-tc_{s})\frac{\partial A}{\partial\rho}+Ac_{s}(3t+\rho)$
$\displaystyle=0$ (33)

Assumption of the form of the scaling solution $v_{\rho}=\rho/t$, along with
Eqs. (28) and (29), puts constraints on $A(\rho,t)$. The form of $A(\rho,t)$
must be such that it satisfies the above differential equations. By adding and
subtracting Eqs. (32) and (33), we obtain a simpler set of equations given by

$\displaystyle\frac{\partial A}{\partial t}+\frac{\rho}{t}\frac{\partial
A}{\partial\rho}+\frac{3Ac_{s}}{t}=0,$ (34) $\displaystyle\frac{\partial
A}{\partial t}+\frac{t}{\rho}\frac{\partial A}{\partial\rho}-\frac{A}{t}=0.$
(35)

The above equations can be solved using the method of characteristics.

Let us first consider the partial differential equation in Eq. (34). The
characteristic curves are obtained by solving the equations

$\frac{dt}{t}=\frac{d\rho}{\rho}=-\frac{dA}{3\,c_{s}A}.$ (36)

It is interesting to note that the above equations for characteristic curve
are similar in structure to Eq. (12), suggesting that they belong to a broader
family of solutions [12]. The solution of the above equations is given by

$A(\rho,t)\sim t^{-3c_{s}}\,g\\!\left(\frac{\rho}{t}\right),$ (37)

where $g(\rho/t)$ is an arbitrary function. The form of $g(\rho/t)$ is
obtained by substituting Eq. (37) in Eq. (35). Defining $w\equiv\rho/t$, we
get

$\left(\frac{1}{w}-w\right)g^{\prime}(w)-(1+3c_{s})g(w)=0.$ (38)

The above equation can be easily integrated to obtain $g(w)$ which leads to

$A(\rho,t)\sim\left(1-\frac{\rho^{2}}{t^{2}}\right)^{-(1+3c_{s})/2}t^{-3c_{s}}.$
(39)

Using this result for $A(\rho,t)$ in Eqs. (26), (30) and (31), we get

$\Phi(T)\sim-3\,c_{s}\ln(t^{2}-\rho^{2}).$ (40)

Substituting the above expression in Eq. (24), we obtain the solutions for
temperature and entropy density as,

$\displaystyle T(\rho,t)$
$\displaystyle=T_{0}\left[\frac{t_{0}^{2}-\rho_{0}^{2}}{t^{2}-\rho^{2}}\right]^{3c_{s}^{2}/2},$
(41) $\displaystyle s(\rho,t)$
$\displaystyle=s_{0}\left[\frac{t_{0}^{2}-\rho_{0}^{2}}{t^{2}-\rho^{2}}\right]^{3/2}.$
(42)

Note that the solution given in Eqs. (41) and (42) has been derived for the
$z=0$ slice. It is then straightforward to show that

$\displaystyle T(\rho,z,t)$
$\displaystyle=T_{0}\left[\frac{t_{0}^{2}-\rho_{0}^{2}-z_{0}^{2}}{t^{2}-\rho^{2}-z^{2}}\right]^{3c_{s}^{2}/2},$
(43) $\displaystyle s(\rho,z,t)$
$\displaystyle=s_{0}\left[\frac{t_{0}^{2}-\rho_{0}^{2}-z_{0}^{2}}{t^{2}-\rho^{2}-z^{2}}\right]^{3/2},$
(44) $\displaystyle v_{\rho}(\rho,z,t)$ $\displaystyle=\frac{\rho}{t},$ (45)
$\displaystyle v_{z}(\rho,z,t)$ $\displaystyle=\frac{z}{t}.$ (46)

form a solution of Eqs. (19)-(21) for all $z$. We note that the scaling
solution, given in Eqs. (43)-(46), is in agreement with the family of
solutions obtained in Ref. [12].

##  III Freezeout and momentum spectra

The hadron momentum spectra can be derived by applying the Cooper-Frye
prescription for particle production at the freezeout hypersurface

$E\frac{dN}{d^{3}p}=\frac{g}{(2\pi)^{3}}\int p_{\mu}d\Sigma^{\mu}f(x,p),$ (47)

where $g$ is the particle degeneracy factor 111Not to be confused with $g(w)$
introduced in the previous section., $E$ is the relativistic particle energy,
$d\Sigma_{\mu}$ is the freeze-out hyper-surface, $p^{\mu}$ is the particle
momenta and $f(x,p)$ is the phase-space distribution function of the
particles. In this work, we consider the Maxwell-Jüttner form of the
equilibrium distribution function, $f(x,p)=\exp[-\beta\,u\cdot p]$, with
$\beta\equiv 1/T$ being the inverse temperature. To proceed further, we need
to specify the freeze-out hypersurface for the case of spherical and
cylindrical expansion. Defining $(3+1)$-dimensional proper time
$\tau_{3}=\sqrt{t^{2}-r^{2}}=\sqrt{t^{2}-\rho^{2}-z^{2}}$, we see from Eqs.
(16) and (43) that analytical solution for temperature, in case of both
spherical and cylindrical expansion, leads to

$T(\tau_{3})=T_{0}\left(\frac{\tau_{30}}{\tau_{3}}\right)^{3c_{s}^{2}/2}.$
(48)

From here, we can easily conclude that $T_{0}$ is the initial temperature at
initial proper time $\tau_{30}$. The above equation also suggests that the
scaling nature of the solution determines the initial temperature profile of
the fireball.

Kinetic freezeout in heavy-ion collision is considered to occur when the
constituent hadrons stop interacting and free stream to the detectors. In
hydrodynamic simulations, this is imposed when the temperature reaches the
kinetic freezeout temperature due to the cooling caused by the expansion of
the fireball. In the present work, we also assume that the freezeout
hypersurface is specified at a constant temperature $T_{f}$. From Eq. (48), we
see that the $T=T_{f}$ hypersurface uniquely corresponds to a constant proper
time hypersurface $\tau_{3}=\tau_{3f}$. In the following, we calculate this
freezeout hypersurface and the corresponding hadron momentum spectra for the
two cases of spherical and cylindrical expansion.

###  III.1 Spherical geometry

Points on a spherically symmetric hypersurface can be parametrized as [60]

$x^{\mu}=\left[t(\zeta),r(\zeta)\sin\theta\cos\phi,r(\zeta)\sin\theta\sin\phi,r(\zeta)\cos\theta\right],$
(49)

where $r,\,\theta,\,\phi$ are usual spherical polar coordinate variables. Here
$\zeta$ is the parameter which defines the spherical hypersurface with
$0<\zeta<1$ such that $r(0)=0$ and $r(1)=R$, with $R$ being the radius of the
fireball at freezeout. The surface element of the hypersurface can be obtained
as

$d\Sigma^{\mu}=\left(\frac{dr}{d\zeta},\frac{dt}{d\zeta}\sin\theta\cos\phi,\frac{dt}{d\zeta}\sin\theta\sin\phi,\frac{dt}{d\zeta}\cos\theta\right)r^{2}\sin\theta\
d\zeta\ d\theta\ d\phi.$ (50)

For radially expanding fireball, the spherically symmetric flow velocity
profile can be written as

$u^{\mu}=\gamma(\zeta)\left[1,v_{r}(\zeta)\sin\theta\cos\phi,v_{r}(\zeta)\sin\theta\sin\phi,v_{r}(\zeta)\cos\theta\right],$
(51)

Where we recall that we have Hubble-like solution for the fluid flow velocity,
$v_{r}(\zeta)=\frac{r(\zeta)}{t}$.

Using spherical polar co-ordinate variables in momentum space, the four-
momentum of a hadron, $p^{\mu}=(E,\,p^{x},\,p^{y},\,p^{z})$, can be written as

$p^{\mu}=(E,p\sin\theta_{p}\cos\phi_{p},p\sin\theta_{p}\sin\phi_{p},p\cos\theta_{p}),$
(52)

where $\theta_{p}$ and $\phi_{p}$ denote the polar and azimuthal angle,
respectively, in the momentum space. Spherical symmetry can be employed to
restrict our considerations to $\theta_{p}=0$, which leads to

$\displaystyle u\cdot p$
$\displaystyle=\gamma(\zeta)\,E-\gamma(\zeta)\,v_{r}(\zeta)\,p\cos\theta,$
(53) $\displaystyle p\cdot d\Sigma$
$\displaystyle=\left(E\,\frac{dr}{d\zeta}-p\cos\theta\,\frac{dt}{d\zeta}\right)r^{2}(\zeta)\sin\theta\,d\zeta\,d\theta\,d\phi.$
(54)

Substituting the above expressions in Eq. (47), along with the classical
Maxwell-Jüttner form of the phase-space distribution function, we get

$E\frac{dN}{d^{3}p}=\frac{g}{(2\pi)^{2}}\int_{0}^{1}d\zeta\
r^{2}\int_{0}^{\pi}d\theta\,\sin\theta\left(E\,\frac{dr}{d\zeta}-p\cos\theta\,\frac{dt}{d\zeta}\right)e^{-\beta
u\cdot p}.$ (55)

We notice that the $\theta$ integral in the above equation can also be done
analytically.

After performing the $\theta$ integral in Eq. (55), the expression for the
momentum distribution is given by

$E\frac{dN}{d^{3}p}=\frac{g}{(2\pi)^{2}}\int_{0}^{1}e^{-\beta\gamma
E}\left[E\frac{dr}{d\zeta}\frac{\sinh(\alpha)}{\alpha}+T\frac{dt}{d\zeta}\left(\frac{\sinh(\alpha)-\alpha\cosh(\alpha)}{\alpha\gamma
v_{r}}\right)\right]r^{2}\,d\zeta,$ (56)

where $\alpha\equiv\beta\gamma v_{r}p$. To proceed further, we make use of the
condition for constant temperature freezeout hypersurface resulting in
$\tau^{2}_{3f}=\text{const.}=t^{2}-r^{2}$. This choice of freezeout
hypersurface leads to

$\frac{dt}{d\zeta}=\frac{r}{t}\,\frac{dr}{d\zeta}.$ (57)

Using the above relation, we obtain the final form of the particle momentum
distribution to be

$E\frac{dN}{d^{3}p}=\frac{gR^{3}}{2\pi^{2}}\int_{0}^{1}e^{-\beta
E\sqrt{\nu^{2}+\chi^{2}}/\nu}\left[E\,\frac{\sinh(a\chi)}{a\chi}+T\nu\left(\frac{\sinh(a\chi)-a\chi\cosh(a\chi)}{a\chi\sqrt{\nu^{2}+\chi^{2}}}\right)\right]\chi^{2}\
d\chi,$ (58)

where $\chi\equiv r/R$, $\nu\equiv\tau_{3\text{f}}/R$, and $a\equiv p/(T\nu)$.
At this point, we would like to emphasize that the above expression for hadron
spectrum is exact analytical hydrodynamic result for spherically expanding
fireball produced in relativistic heavy-ion collisions.

It is instructive to calculate the average radial velocity of the fireball in
terms of the parameters on the freezeout hypersurface. The average radial
velocity can be defined as

$\langle v_{r}\rangle=\frac{\int
v_{r}\sqrt{d\Sigma_{\mu}d\Sigma^{\mu}}}{\int\sqrt{d\Sigma_{\mu}d\Sigma^{\mu}}}.$
(59)

With $v_{r}=r/t$ and employing Eqs. (50) and (57), we obtain the final
expression for the average radial velocity as

$\langle
v_{r}\rangle=\frac{1+\nu^{2}\log\left(\frac{\nu^{2}}{1+\nu^{2}}\right)}{\sqrt{1+\nu^{2}}-\nu^{2}\log\left[\frac{\nu}{\sqrt{1+\nu^{2}}-1}\right]}.$
(60)

It is straightforward to verify that the maximum value of $\langle
v_{r}\rangle$ is $1$ for $\nu=0$.

###  III.2 Cylindrical geometry

Points on a longitudinally boost-invariant and cylindrically symmetric
hypersurface can be parametrized as [60]

$x^{\mu}=\left[\tau(\zeta)\cosh\eta_{s},\,\rho(\zeta)\cos\phi,\,\rho(\zeta)\sin\phi,\,\tau(\zeta)\sinh\eta_{s}\right]$
(61)

where $\rho=\sqrt{x^{2}+y^{2}}$ and $\tau=\sqrt{t^{2}-z^{2}}$. Here again,
$\zeta$ is the parameter which defines the cylindrical hypersurface with
$0<\zeta<1$ such that $\rho(0)=0$ and $\rho(1)=R$, with $R$ being the
transverse radius of the fireball at freezeout 222We use the same notation for
transverse freezeout radius as well as spherical radius at freezeout as
introduced in the previous section, for reasons that will become apparent
later.. Note that space-time rapidity $\eta_{s}$ is independent of the
parameter $\zeta$ owing to the boost-invariance symmetry considered here. The
surface element of the hypersurface can be obtained as

$d\Sigma^{\mu}=\left(\frac{d\rho}{d\zeta}\cosh\eta_{s},\,\frac{d\tau}{d\zeta}\cos\phi,\,\frac{d\tau}{d\zeta}\sin\phi,\,\frac{d\rho}{d\zeta}\sinh\eta_{s}\right)\tau(\zeta)\rho(\zeta)d\zeta
d\phi d\eta_{s}.$ (62)

Defining longitudinal fluid flow rapidity, $\eta_{f}$, and transverse flow
rapidity, $\eta_{T}$ as

$\displaystyle\eta_{f}$
$\displaystyle=\frac{1}{2}\log\left(\frac{1+v_{z}}{1-v_{z}}\right),$ (63)
$\displaystyle\eta_{T}$
$\displaystyle=\tanh^{-1}\left(\frac{v_{\rho}}{\sqrt{1-v_{z}^{2}}}\right),$
(64)

the fluid four-velocity can be written as

$u^{\mu}=\left(\cosh\eta_{T}\cosh\eta_{f},\,\sinh\eta_{T}\cos\phi,\,\sinh\eta_{T}\sin\phi,\,\cosh\eta_{T}\sinh\eta_{f}\right).$
(66)

In terms of the transverse momentum and rapidity variables, the four-momentum
of a particle, $p^{\mu}=(E,\,p^{x},\,p^{y},\,p^{z})$, is given by

$p^{\mu}=\left(m_{T}\cosh
y,\,p_{T}\cos\phi_{p},\,p_{T}\sin\phi_{p},\,m_{T}\sinh y\right),$ (67)

where $y$ is the longitudinal rapidity of the particle and $\phi_{p}$ is the
azimuthal angle in momentum space. Here $m_{T}\equiv\sqrt{m^{2}+p_{T}^{2}}$ is
the transverse mass with $p_{T}$ being the transverse momentum of the
particles.

The longitudinal boost-invariance for fluid velocity profile is guaranteed for
$v_{z}=z/t$ which leads to $\eta_{f}=\eta_{s}$. Imposing this boost-invariance
condition, we have

$\displaystyle u\cdot p$
$\displaystyle=m_{T}\cosh(\eta_{T})\cosh(y-\eta_{s})-p_{T}\sinh(\eta_{T})\cos(\phi-\phi_{p}),$
(68) $\displaystyle p\cdot d\Sigma$
$\displaystyle=\left[m_{T}\cosh(y-\eta_{s})\frac{d\rho}{d\zeta}-p_{T}\cos(\phi-\phi_{p})\frac{d\tau}{d\zeta}\right]\tau(\zeta)\,\rho(\zeta)\,d\zeta\,d\phi\,d\eta_{s}.$
(69)

Substituting the above expressions in Eq. (47), along with the classical
Maxwell-Jüttner form of the phase-space distribution function, and performing
the $\phi$ and $\eta_{s}$ integrations, we get

$E\,\frac{dN}{d^{3}p}=\frac{g}{2\pi^{2}}\\!\int_{0}^{1}\\!\left[m_{T}\,I_{0}(\beta
p_{T}\sinh\eta_{T})\,K_{1}(\beta
m_{T}\cosh\eta_{T})\frac{d\rho}{d\zeta}-p_{T}\,I_{1}(\beta
p_{T}\sinh\eta_{T})\,K_{0}(\beta
m_{T}\cosh\eta_{T})\frac{d\tau}{d\zeta}\right]\rho\,\tau\,d\zeta,$ (70)

where, $I_{n}(x)$ and $K_{n}(x)$ are modified Bessel functions of first kind
and second kind, respectively.

Figure 1: Average radial and transverse velocity as a function of the
freezeout parameter $\nu$. These average velocities are calculated using Eqs.
(60) and (74)for spherical and cylindrical expansion geometries, respectively.

To proceed further, we make use of the condition for constant temperature
freezeout hypersurface resulting in
$\tau^{2}_{3f}=\text{const.}=\tau^{2}-\rho^{2}$. This choice of freezeout
hypersurface leads to

$\frac{d\tau}{d\zeta}=\frac{\rho}{\tau}\,\frac{d\rho}{d\zeta}.$ (71)

Considering Hubble-like flow in the transverse plane, i.e., $v_{\rho}=\rho/t$
leads to $\rho=\tau_{3f}\sinh\eta_{T},\,\tau=\tau_{3f}\cosh\eta_{T}$ on the
freezeout hypersurface. Using Eq. (71) and these relations, we obtain the
final form of the particle momentum distribution to be

$E\,\frac{dN}{d^{3}p}=\frac{gR^{3}}{2\pi^{2}}\\!\int_{0}^{1}\\!\left[m_{T}\sqrt{\nu^{2}+\chi^{2}}I_{0}\\!\left(\\!\frac{\beta
p_{T}\chi}{\nu}\\!\right)K_{1}\\!\left(\\!\beta
m_{T}\sqrt{1+\frac{\chi^{2}}{\nu^{2}}}\right)-\chi
p_{T}I_{1}\\!\left(\\!\frac{\beta
p_{T}\chi}{\nu}\\!\right)K_{0}\\!\left(\\!\beta
m_{T}\sqrt{1+\frac{\chi^{2}}{\nu^{2}}}\right)\right]\chi\,d\chi,$ (72)

where $\chi\equiv\rho/R$ and $\nu\equiv\tau_{3f}/R$. Comparing the above
equation with Eq. (58), we find that the form of the overall normalization
factor, appearing outside the integration on r.h.s, is identical. At this
juncture, we would again like to emphasize that the above expression for the
hadron spectrum is exact analytical hydrodynamic result for a cylindrically
expanding fireball produced in relativistic heavy-ion collisions.

Further, it is informative to calculate the average transverse velocity of the
fireball in terms of the parameters on the freezeout hypersurface. The average
transverse velocity can be defined as

$\langle v_{T}\rangle=\frac{\int
v_{\rho}\sqrt{d\Sigma_{\mu}d\Sigma^{\mu}}}{\int\sqrt{d\Sigma_{\mu}d\Sigma^{\mu}}}.$
(73)

With $v_{\rho}=\rho/t$ and employing Eqs. (62) and (71), we obtain the final
expression for the average transverse velocity as

$\langle
v_{T}\rangle=\sqrt{1+\nu^{2}}-\nu^{2}\log\left[\frac{\nu}{\sqrt{1+\nu^{2}}-1}\right].$
(74)

It is again straightforward to verify that the maximum value of $\langle
v_{T}\rangle$ is $1$ for $\nu=0$.

It is important to understand the physical meaning of the freezeout parameter
$\nu$. To achieve this, we have calculated the average radial and transverse
velocities in Eqs. (60) and (74). In Fig. 1, we show the $\nu$-dependence of
average radial velocity for spherical expansion, given in Eq. (60), as well as
average transverse velocity for cylindrical expansion, given in Eq. (74). We
find a monotonic decrease of the average velocities as a function of $\nu$;
recall that $\nu=\tau_{3f}/R$ has a dimension of inverse velocity. Therefore
the observed monotonic decrease is expected. Moreover, we see that for the
same value of $\nu$, the average radial velocity is slightly larger than
average transverse velocity. In the next section, we shall use the expressions
for the hadron momentum distribution, Eqs. (58) and (72), to fit the
experimental data for the transverse momentum spectrum of hadrons. The fit
parameters, i.e., temperature $T$ and $\nu$ which represent average expansion
velocity, would be interpreted in the context of experimental results.

##  IV Numerical results and discussions

Figure 2:  (Left) Mid-rapidity $p_{T}$ spectra for $\pi^{-}$(black), $k^{+}$
(red) and $k^{-}$ (blue) at SPS collision energy of 30 AGeV. Square boxes
denote the data from Ref. [61] and the line represents the model estimations.
(Right) $p_{T}$ integrated rapidity spectra $dN/dy$ of the same species of
particles calculated using parameters from fit to $p_{T}$ spectra.

In this section, we apply our exact analytical results for hadron momentum
spectra to fit experimental data. In low-energy heavy-ion collisions, one
expects the expansion geometry of the fireball to be spherical. To this end,
we apply our analytical results for spherical geometry to Pb-Pb collisions at
$30$ AGeV in Super Proton Synchrotron (SPS). We consider $y=0$ in Eq. (58) to
obtain the transverse momentum, $p_{T}$, distributions corresponding to the
mid-rapidity range. We then employ chi-square minimization to obtain the best-
fit result by analyzing the SPS experimental data for mid-rapidity $p_{T}$
spectra. In this case, apart from the overall normalization, there are two fit
parameters, $T$ and $\nu$ (or equivalently $\langle v_{r}\rangle$). In the
left panel of Fig. 2, we have successfully characterized the charged pion and
kaon spectra using our calculated spectra, with a freeze-out temperature of
$81$ MeV and average radial velocity $\langle v_{r}\rangle=0.65$. These fitted
results are in agreement with earlier analyses [56, 57].

With a suitable fitting of the $p_{T}$ spectra with a reduced $\chi^{2}$ of
$5.6$, it would be interesting to check the momentum-integrated rapidity
spectra. The deviation from boost-invariance necessitates that the rapidity
spectra exhibit a Gaussian distribution, initially proposed by Landau [1].
Recently, in Ref. [10], a more generalized form of rapidity spectra has been
proposed, accommodating a non-conformal equation of state. In the right column
of Fig. 2, we present the rapidity spectra, computed using the parameters
derived from the fitting of the mid-rapidity differential $p_{T}$ data, as
discussed earlier. While the Gaussian-like structure is apparent in the
resulting rapidity distribution, there exists a noticeable quantitative
discrepancy with the data at higher rapidities. This discrepancy can be
attributed to the fact that, for fitting the $p_{T}$ spectra, mid-rapidity
data was employed, whereas the same parameter set was utilized for computing
this differential rapidity variation across the entire rapidity range.
Moreover, the assumption of a completely spherical geometry for this SPS
collision energy may not be valid. For a suitable agreement between the data
and model estimation of the rapidity spectra, one has to perform fits to the
$p_{T}$ spectra for different rapidity windows, as discussed in Ref.[62],
which is beyond the scope of the present study.

Figure 3: (Left) Mid-rapidity $p_{T}$ spectra for $\pi^{-}$(black), $k^{-}$
(blue) and $p$ (red) at LHC collision energy of 2.76 ATeV. Square boxes denote
the data from Ref. [63] and the line represents the model estimations.

At ultra-relativistic high energies, the system exhibits longitudinal boost-
invariance, a characteristic we have incorporated into our analytical
hydrodynamics framework by adopting cylindrical geometry with boost-invariance
along the $z$ direction. This solution has been translated into the form of
transverse momentum spectra as given in Eq. (72). We have applied this result
to fit the transverse momentum spectra observed in LHC collisions at
$\sqrt{s_{NN}}=2.76$ TeV energy. Utilizing mid-rapidity $p_{T}$ spectra data
for pions, kaons, and protons, we conducted a chi-square minimization process.
The fitted spectra are plotted along with the data in Fig. 3. The resulting
temperature is found to be $112$ MeV, which is higher compared to the value
observed at SPS ($81$ MeV). These reference values align with our general
expectation of a higher freeze-out temperature at LHC due to the increased
collision energies. Additionally, we obtain $\langle v_{T}\rangle=0.72$ at LHC
which aligns with the conventional blast wave parameterization detailed in
Ref. [63]. In contrast, the average radial velocity ($\langle v_{r}\rangle$)
at SPS energies is lower, indicating a higher value of $\nu$. This observation
suggests a relatively later freeze-out time at lower collision energies.

##  V Summary and conclusions

Analytical hydrodynamic solutions play a crucial role in establishing a
benchmark for understanding the evolution of QCD matter created in heavy-ion
collisions. While providing a general solution remains challenging, one can
focus on specific geometries and flow profiles that facilitate formulation. In
this study, we obtained analytical solutions for relativistic hydrodynamics in
a system characterized by cylindrical symmetry with boost-invariant
longitudinal expansion and Hubble-like transverse expansion. Additionally, we
explored an analytical solution for a Hubble-like spherically expanding
system. These straightforward profiles allowed us to derive a simple
analytical solution for the temperature evolution relative to the proper time,
facilitating a connection with the freeze-out hypersurface. For both
cylindrical and spherical flow profiles, we computed the analytical
expressions for the transverse momentum spectra of hadrons on a constant
temperature freeze-out hypersurface, utilizing the Cooper-Frye prescription.

We emphasize that the expressions obtained for hadron momentum distributions
are exact analytical hydrodynamic results. We compared these results for
transverse momentum spectra with experimental data from SPS and LHC energies.
Adhering to the symmetry of the expanding system, we have employed the
spherical solution for lower-energy collisions at SPS while opting for the
cylindrical description at the LHC. We used $\chi^{2}$ minimization to analyze
available $p_{T}$ spectra data for SPS collisions at $30$ AGeV and LHC
collisions at $2.76$ TeV. The successful description of the data with our
analytical results provides a validation to this solution. However, we note
that the inclusion of equation of state in the analytical solutions to
hydrodynamic equations are rather restrictive and does not allow for
incorporating realistic equation of state such as lattice QCD results. On the
other hand, it may be possible to extract an average speed-of-sound from
analysis of experimental observables. A more detailed phenomenology of heavy-
ion collision experiments, within the presented framework, is left for future
work.

## References

  * [1] L. D. Landau, On the multiparticle production in high-energy collisions, Izv. Akad. Nauk Ser. Fiz. 17 (1953) 51–64.
  * [2] C.-Y. Wong, Landau Hydrodynamics Revisited, Phys. Rev. C78 (2008) 054902.  arXiv:0808.1294, doi:10.1103/PhysRevC.78.054902.
  * [3] R. C. Hwa, Statistical Description of Hadron Constituents as a Basis for the Fluid Model of High-Energy Collisions, Phys. Rev. D10 (1974) 2260.  doi:10.1103/PhysRevD.10.2260.
  * [4] J. D. Bjorken, Highly Relativistic Nucleus-Nucleus Collisions: The Central Rapidity Region, Phys. Rev. D 27 (1983) 140–151.  doi:10.1103/PhysRevD.27.140.
  * [5] M. J. Murray, Scanning the phases of QCD with BRAHMS, J. Phys. G30 (2004) S667–S674.  arXiv:nucl-ex/0404007, doi:10.1088/0954-3899/30/8/004.
  * [6] I. G. Bearden, et al., Charged meson rapidity distributions in central Au+Au collisions at s(NN)**(1/2) = 200-GeV, Phys. Rev. Lett. 94 (2005) 162301.  arXiv:nucl-ex/0403050, doi:10.1103/PhysRevLett.94.162301.
  * [7] M. Murray, Flavor dynamics, J. Phys. G35 (2008) 044015.  arXiv:0710.4576, doi:10.1088/0954-3899/35/4/044015.
  * [8] P. Steinberg, Bulk dynamics in heavy ion collisions, Nucl. Phys. A752 (2005) 423–432.  arXiv:nucl-ex/0412009, doi:10.1016/j.nuclphysa.2005.02.139.
  * [9] P. Steinberg, Entropy Production at High Energy and mu(B), PoS CPOD2006 (2006) 036.  arXiv:nucl-ex/0702019, doi:10.22323/1.029.0036.
  * [10] D. Biswas, K. Deka, A. Jaiswal, S. Roy, Viscosity, nonconformal equation of state, and sound velocity in Landau hydrodynamics, Phys. Rev. C 102 (1) (2020) 014912.  arXiv:1910.13368, doi:10.1103/PhysRevC.102.014912.
  * [11] S. V. Akkelin, T. Csorgo, B. Lukacs, Y. M. Sinyukov, M. Weiner, Simple solutions of fireball hydrodynamics for selfsimilar elliptic flows, Phys. Lett. B 505 (2001) 64–70.  arXiv:hep-ph/0012127, doi:10.1016/S0370-2693(01)00351-3.
  * [12] T. Csörgő, F. Grassi, Y. Hama, T. Kodama, Simple solutions of relativistic hydrodynamics for longitudinally and cylindrically expanding systems, Phys. Lett. B 565 (2003) 107–115.  arXiv:nucl-th/0305059, doi:10.1016/S0370-2693(03)00747-0.
  * [13] T. Csorgo, M. I. Nagy, M. Csanad, A New family of simple solutions of perfect fluid hydrodynamics, Phys. Lett. B 663 (2008) 306–311.  arXiv:nucl-th/0605070, doi:10.1016/j.physletb.2008.04.038.
  * [14] S. S. Gubser, Symmetry constraints on generalizations of Bjorken flow, Phys. Rev. D 82 (2010) 085027.  arXiv:1006.0006, doi:10.1103/PhysRevD.82.085027.
  * [15] T. Csörgő, M. I. Nagy, New family of exact and rotating solutions of fireball hydrodynamics, Phys. Rev. C 89 (4) (2014) 044901.  arXiv:1309.4390, doi:10.1103/PhysRevC.89.044901.
  * [16] H. Bantilan, T. Ishii, P. Romatschke, Holographic Heavy-Ion Collisions: Analytic Solutions with Longitudinal Flow, Elliptic Flow and Vorticity, Phys. Lett. B 785 (2018) 201–206.  arXiv:1803.10774, doi:10.1016/j.physletb.2018.08.038.
  * [17] T. Hirano, K. Tsuda, Collective flow and two pion correlations from a relativistic hydrodynamic model with early chemical freezeout, Phys. Rev. C 66 (2002) 054905.  arXiv:nucl-th/0205043, doi:10.1103/PhysRevC.66.054905.
  * [18] P. F. Kolb, U. W. Heinz, Hydrodynamic description of ultrarelativistic heavy ion collisions (2003) 634–714arXiv:nucl-th/0305084.
  * [19] P. Romatschke, U. Romatschke, Viscosity Information from Relativistic Nuclear Collisions: How Perfect is the Fluid Observed at RHIC?, Phys. Rev. Lett. 99 (2007) 172301.  arXiv:0706.1522, doi:10.1103/PhysRevLett.99.172301.
  * [20] H. Song, U. W. Heinz, Causal viscous hydrodynamics in 2+1 dimensions for relativistic heavy-ion collisions, Phys. Rev. C 77 (2008) 064901.  arXiv:0712.3715, doi:10.1103/PhysRevC.77.064901.
  * [21] M. Luzum, P. Romatschke, Conformal Relativistic Viscous Hydrodynamics: Applications to RHIC results at s(NN)**(1/2) = 200-GeV, Phys. Rev. C78 (2008) 034915, [Erratum: Phys. Rev.C79,039903(2009)].  arXiv:0804.4015, doi:10.1103/PhysRevC.78.034915,10.1103/PhysRevC.79.039903.
  * [22] M. Luzum, P. Romatschke, Viscous Hydrodynamic Predictions for Nuclear Collisions at the LHC, Phys. Rev. Lett. 103 (2009) 262302.  arXiv:0901.4588, doi:10.1103/PhysRevLett.103.262302.
  * [23] H. Song, S. A. Bass, U. Heinz, T. Hirano, C. Shen, 200 A GeV Au+Au collisions serve a nearly perfect quark-gluon liquid, Phys. Rev. Lett. 106 (2011) 192301, [Erratum: Phys. Rev. Lett.109,139904(2012)].  arXiv:1011.2783, doi:10.1103/PhysRevLett.106.192301,%****␣main.bbl␣Line␣175␣****10.1103/PhysRevLett.109.139904.
  * [24] B. Schenke, S. Jeon, C. Gale, Elliptic and triangular flow in event-by-event (3+1)D viscous hydrodynamics, Phys. Rev. Lett. 106 (2011) 042301.  arXiv:1009.3244, doi:10.1103/PhysRevLett.106.042301.
  * [25] M. Luzum, Elliptic flow at energies available at the CERN Large Hadron Collider: Comparing heavy-ion data to viscous hydrodynamic predictions, Phys. Rev. C83 (2011) 044911.  arXiv:1011.5173, doi:10.1103/PhysRevC.83.044911.
  * [26] B. Schenke, S. Jeon, C. Gale, Anisotropic flow in $\sqrt{s}=2.76$ TeV Pb+Pb collisions at the LHC, Phys. Lett. B702 (2011) 59–63.  arXiv:1102.0575, doi:10.1016/j.physletb.2011.06.065.
  * [27] C. Gale, S. Jeon, B. Schenke, P. Tribedy, R. Venugopalan, Event-by-event anisotropic flow in heavy-ion collisions from combined Yang-Mills and viscous fluid dynamics, Phys. Rev. Lett. 110 (1) (2013) 012302.  arXiv:1209.6330, doi:10.1103/PhysRevLett.110.012302.
  * [28] U. Heinz, R. Snellings, Collective flow and viscosity in relativistic heavy-ion collisions, Ann. Rev. Nucl. Part. Sci. 63 (2013) 123–151.  arXiv:1301.2826, doi:10.1146/annurev-nucl-102212-170540.
  * [29] C. Gale, S. Jeon, B. Schenke, Hydrodynamic Modeling of Heavy-Ion Collisions, Int. J. Mod. Phys. A 28 (2013) 1340011.  arXiv:1301.5893, doi:10.1142/S0217751X13400113.
  * [30] R. S. Bhalerao, A. Jaiswal, S. Pal, Collective flow in event-by-event partonic transport plus hydrodynamics hybrid approach, Phys. Rev. C92 (1) (2015) 014903\.  arXiv:1503.03862, doi:10.1103/PhysRevC.92.014903.
  * [31] A. Jaiswal, V. Roy, Relativistic hydrodynamics in heavy-ion collisions: general aspects and recent developments, Adv. High Energy Phys. 2016 (2016) 9623034\.  arXiv:1605.08694, doi:10.1155/2016/9623034.
  * [32] F. Cooper, G. Frye, Comment on the Single Particle Distribution in the Hydrodynamic and Statistical Thermodynamic Models of Multiparticle Production, Phys. Rev. D 10 (1974) 186.  doi:10.1103/PhysRevD.10.186.
  * [33] P. J. Siemens, J. O. Rasmussen, Evidence for a blast wave from compress nuclear matter, Phys. Rev. Lett. 42 (1979) 880–887.  doi:10.1103/PhysRevLett.42.880.
  * [34] K. S. Lee, U. W. Heinz, Collective Flow Model for the Pion Transverse Momentum Spectra From Relativistic Nuclear Collisions at CERN, Z. Phys. C 43 (1989) 425\.  doi:10.1007/BF01506538.
  * [35] K. S. Lee, U. W. Heinz, E. Schnedermann, Search for Collective Transverse Flow Using Particle Transverse Momentum Spectra in Relativistic Heavy Ion Collisions, Z. Phys. C 48 (1990) 525–541.  doi:10.1007/BF01572035.
  * [36] E. Schnedermann, J. Sollfrank, U. W. Heinz, Thermal phenomenology of hadrons from 200-A/GeV S+S collisions, Phys. Rev. C 48 (1993) 2462–2475.  arXiv:nucl-th/9307020, doi:10.1103/PhysRevC.48.2462.
  * [37] W. Florkowski, W. Broniowski, Hydro-inspired parameterizations of freeze-out in relativistic heavy-ion collisions, Acta Phys. Polon. B 35 (2004) 2895–2910.  arXiv:nucl-th/0410081.
  * [38] W. Florkowski, Particle spectra and hydro-inspired models, Nucl. Phys. A 774 (2006) 179–188.  arXiv:nucl-th/0509039, doi:10.1016/j.nuclphysa.2006.06.039.
  * [39] F. Retiere, M. A. Lisa, Observable implications of geometrical and dynamical aspects of freeze out in heavy ion collisions, Phys. Rev. C 70 (2004) 044907\.  arXiv:nucl-th/0312024, doi:10.1103/PhysRevC.70.044907.
  * [40] T. Csorgo, B. Lorstad, Bose-Einstein correlations for three-dimensionally expanding, cylindrically symmetric, finite systems, Phys. Rev. C 54 (1996) 1390–1403.  arXiv:hep-ph/9509213, doi:10.1103/PhysRevC.54.1390.
  * [41] T. Csorgo, Particle interferometry from 40-MeV to 40-TeV, Acta Phys. Hung. A 15 (2002) 1–80.  arXiv:hep-ph/0001233, doi:10.1556/APH.15.2002.1-2.1.
  * [42] H. Dobler, J. Sollfrank, U. W. Heinz, Kinetic freezeout and radial flow in 11.6-GeV/A Au + Au collisions, Phys. Lett. B 457 (1999) 353–358.  arXiv:nucl-th/9904018, doi:10.1016/S0370-2693(99)00551-1.
  * [43] P. Huovinen, P. F. Kolb, U. W. Heinz, P. V. Ruuskanen, S. A. Voloshin, Radial and elliptic flow at RHIC: Further predictions, Phys. Lett. B 503 (2001) 58–64.  arXiv:hep-ph/0101136, doi:10.1016/S0370-2693(01)00219-2.
  * [44] C. Adler, et al., Identified particle elliptic flow in Au + Au collisions at s(NN)**(1/2) = 130-GeV, Phys. Rev. Lett. 87 (2001) 182301.  arXiv:nucl-ex/0107003, doi:10.1103/PhysRevLett.87.182301.
  * [45] M. Csanad, T. Csorgo, B. Lorstad, Buda-Lund hydro model for ellipsoidally symmetric fireballs and the elliptic flow at RHIC, Nucl. Phys. A 742 (2004) 80–94.  arXiv:nucl-th/0310040, doi:10.1016/j.nuclphysa.2004.06.023.
  * [46] D. Teaney, The Effects of viscosity on spectra, elliptic flow, and HBT radii, Phys. Rev. C 68 (2003) 034913.  arXiv:nucl-th/0301099, doi:10.1103/PhysRevC.68.034913.
  * [47] J. G. Cramer, G. A. Miller, J. M. S. Wu, J.-H. Yoon, Quantum opacity, the RHIC HBT puzzle, and the chiral phase transition, Phys. Rev. Lett. 94 (2005) 102302, [Erratum: Phys.Rev.Lett. 95, 139901 (2005)].  arXiv:nucl-th/0411031, doi:10.1103/PhysRevLett.94.102302.
  * [48] G. A. Miller, J. G. Cramer, Polishing the lens: I. Pionic final state interactions and HBT correlations: Distorted wave emission function (DWEF) formalism and examples, J. Phys. G 34 (2007) 703–740.  arXiv:nucl-th/0507004, doi:10.1088/0954-3899/34/4/009.
  * [49] T. Renk, Model for the spacetime evolution of 200-a-gev au-au collisions, Phys. Rev. C 70 (2004) 021903.  arXiv:hep-ph/0404140, doi:10.1103/PhysRevC.70.021903.
  * [50] T. Renk, A Comprehensive description of multiple observables in heavy ion collisions at SPS, J. Phys. G 30 (2004) 1495.  arXiv:hep-ph/0403239, doi:10.1088/0954-3899/30/10/015.
  * [51] W. Broniowski, W. Florkowski, Explanation of the RHIC p(T) spectra in a thermal model with expansion, Phys. Rev. Lett. 87 (2001) 272302.  arXiv:nucl-th/0106050, doi:10.1103/PhysRevLett.87.272302.
  * [52] W. Broniowski, W. Florkowski, Strange particle production at RHIC in a single freezeout model, Phys. Rev. C 65 (2002) 064905.  arXiv:nucl-th/0112043, doi:10.1103/PhysRevC.65.064905.
  * [53] W. Broniowski, A. Baran, W. Florkowski, Thermal approach to RHIC, Acta Phys. Polon. B 33 (2002) 4235–4258.  arXiv:hep-ph/0209286.
  * [54] P. Ghosh, S. Muhuri, J. K. Nayak, R. Varma, Indication of transverse radial flow in high-multiplicity proton–proton collisions at the Large Hadron Collider, J. Phys. G 41 (2014) 035106.  arXiv:1402.6813, doi:10.1088/0954-3899/41/3/035106.
  * [55] A. Jaiswal, V. Koch, A viscous blast-wave model for relativistic heavy-ion collisions (8 2015).  arXiv:1508.05878.
  * [56] S. P. Rode, P. P. Bhaduri, A. Jaiswal, A. Roy, Kinetic freeze-out conditions in nuclear collisions with $2A$ \- $158A$ GeV beam energy within a non-boost-invariant blast-wave model, Phys. Rev. C 98 (2) (2018) 024907.  arXiv:1805.11463, doi:10.1103/PhysRevC.98.024907.
  * [57] S. P. Rode, P. P. Bhaduri, A. Jaiswal, Flow fluctuations and kinetic freeze-out of identified hadrons at energies available at the CERN Super Proton Synchrotron, Phys. Rev. C 108 (1) (2023) 014906.  arXiv:2303.10947, doi:10.1103/PhysRevC.108.014906.
  * [58] J. P. Bondorf, S. I. A. Garpman, J. Zimanyi, A Simple Analytic Hydrodynamic Model for Expanding Fireballs, Nucl. Phys. A 296 (1978) 320–332.  doi:10.1016/0375-9474(78)90076-3.
  * [59] G. Baym, B. L. Friman, J. P. Blaizot, M. Soyeur, W. Czyz, Hydrodynamics of Ultrarelativistic Heavy Ion Collisions, Nucl. Phys. A 407 (1983) 541–570.  doi:10.1016/0375-9474(83)90666-8.
  * [60] W. Florkowski, Phenomenology of Ultra-Relativistic Heavy-Ion Collisions, World Scientific Publishing Company, 2010.
  * [61] C. Alt, et al., Pion and kaon production in central Pb + Pb collisions at 20-A and 30-A-GeV: Evidence for the onset of deconfinement, Phys. Rev. C 77 (2008) 024903.  arXiv:0710.0118, doi:10.1103/PhysRevC.77.024903.
  * [62] M. Waqas, H.-M. Chen, G.-X. Peng, A. A. K. H. Ismail, M. Ajaz, Z. Wazir, R. Shehzadi, S. Jamal, A. AbdelKader, Study of Kinetic Freeze-Out Parameters as a Function of Rapidity in pp Collisions at CERN SPS Energies, Entropy 23 (10) (2021) 1363.  arXiv:2109.05831, doi:10.3390/e23101363.
  * [63] B. Abelev, et al., Centrality dependence of $\pi$, K, p production in Pb-Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV, Phys. Rev. C 88 (2013) 044910.  arXiv:1303.0737, doi:10.1103/PhysRevC.88.044910.


# Cell theories for the chiral crystal phase of hard equilateral triangles

Yuri Martínez-Ratón<EMAIL_ADDRESS>Grupo Interdisciplinar de Sistemas
Complejos (GISC), Departamento de Matemáticas, Escuela Politécnica Superior,
Universidad Carlos III de Madrid, Avenida de la Universidad 30, E-28911,
Leganés, Madrid, Spain  Enrique Velasco<EMAIL_ADDRESS>Departamento
de Física Teórica de la Materia Condensada, Instituto de Física de la Materia
Condensada (IFIMAC) and Instituto de Ciencia de Materiales Nicolás Cabrera,
Universidad Autónoma de Madrid, E-28049, Madrid, Spain

###### Abstract

We derive several versions of the cell theory for a crystal phase of hard
equilateral triangles. To that purpose we analytically calculated the free
area of a frozen oriented or freely rotating particle inside the cavity formed
by its neighbours in a chiral configuration of their orientations. From the
most successful versions of the theory we predict an equation of state which,
despite being derived from a crystal configuration of particles, describes
very reasonably the equation of state of the 6-atic liquid-crystal phase at
packing fractions not very close from the isotropic-6-atic bifurcation. Also,
the same equation of state performs well when compared to that from MC
simulations for the stable crystal phase. The agreement can even be improved
by selecting adequate values for the angle of chirality. Despite the success
of two of the versions of the theory for the pressure, we show that the free-
energy is an increasing function of the angle of chirality, implying that the
most stable phase is the achiral phase. Furthermore, we show that possible
clustering effects, such as the formation of perfect chiral hexagonal
clusters, which in turn crystallize into an hexagonal lattice, cannot explain
the presence of the chirality observed in simulations.

##  I Introduction

The last two decades have witnessed the growth of theoretical and experimental
works on the equilibrium liquid-crystal phase behavior of two-dimensional (2D)
hard anisotropic particles. This effort is mainly motivated by the design,
using micro-lithographic techniques, of prism-like particles with different
cross-sectional polygonal geometries such as squares or rectangles, rhombuses,
triangles and kites Zhao1 ; Zhao2 ; Zhao3 ; Zhao4 . These particles, after
being adsorbed into monololayers, form effectively 2D systems of Brownian
particles which approximately interact through entropic forces, making these
systems ideal experimental platforms to explore the effect of particle shape
on the liquid-crystal symmetries in stable phases at equilibrium, and
constitute a test bed for the vast number of theoretical works on 2D liquid
crystals in the last years Schlaken ; Frenkel1 ; Dijkstra1 ; Frenkel2 ; MR1 ;
Donev ; MR2 ; Escobedo ; Varga1 ; Varga2 ; Odriozola ; Heras ; Wittmann .

Of particular interest are the symmetries of the liquid-crystal phases that
can stabilized due to the specific particle shapes. For example, hard
rectangles (and squares) of relatively small aspect ratio stabilize, before
crystallization, a 4-atic (tetratic) phase with two equivalent, perpendicular
directors. The orientational distribution function of this phase has a
fourfold symmetry, $h(\phi)=h(\phi+\pi/2)$ Schlaken ; Frenkel2 ; MR1 ; MR2 .
For larger aspect ratios, rectangles undergo the standard isotropic
(I)-nematic (N) transition before crystallization. Equilateral triangles
stabilize a 6-atic phase with sixfold symmetry, $h(\phi)=h(\phi+\pi/3)$
Dijkstra2 ; MR3 . Hard-right triangles exhibit stable 4-atic liquid-crystal
and crystal phases, but when the fluid is compressed, strong 8-atic
orientational correlations are observed Dijkstra2 ; MR4 . We recently proposed
a version of density-functional theory that captures the formation of clusters
obtained from different monomeral triangular units, that in turn explains the
peculiar phase behavior of hard-right triangles MR5 .

It has been shown recently that some 2D nonchiral particle shapes tend to
generate chiral liquid-crystal or crystalline phases Dijkstra2 ; Zhao6 ; Zhao7
; Zhao8 ; Carmichael . For instance, triangular shapes with rounded corners
have been shown experimentally to generate a chiral 6-atic liquid-crystalline
phase Zhao6 . The presence of roundness and curvature in the edges and corners
of the particles can promote certain configurations of dimers of particles
that explain the presence of chirality Zhao6 ; Carmichael . Note however that
interactions other than entropic, such as electrostatic or polymer-mediated
interactions, could induce the observed chirality. Monte Carlo (MC)
simulations on perfect hard equilateral triangles (HET) also stabilize chiral
crystal phases: in equilibrium two coexisting domains, each one populated by
one of two possible _enantiomers_ , contain equally oriented particles with
respect to the main crystalline axes but with a nonzero chiral angle Dijkstra2
. The coexisting domains are separated by a non-chiral dividing interface
Dijkstra2 .

A physical explanation of chirality in the crystal phase of the purely
entropic HET system is still lacking. A numerically implemented cell theory
was recently used to explain chirality in the crystal phase of hard rhombuses
Zhao8 , and it was shown that the entropic nature of interactions, together
with the particular symmetry of the cell, can explain this phenomenon.

In this work we develop different versions of cell theories for the
crystalline phase of HET. We analytically calculate the free area of a given
particle inside the cell formed by its neighbours in a particular chiral
configurations. While the orientation of the neighbours is fixed, the _caged_
particle can freely rotate. We show that some of these versions predict an
equation of state, not only for the crystalline, but also for the 6-atic
liquid-crystal phase, in reasonable agreement with MC simulations. This fact
shows that the liquid-crystal phase has a local spatial order similar to that
of the crystal. By varying the degree of chirality (making it similar to that
found in simulations) the agreement between theory and simulation can be
improved. A very simple equation of state showing good agreement with MC
results is obtained by assuming a frozen orientation for the caged particle,
but with an extra global factor of $3/2$. We theoretically explain the origin
of this factor.

By contrast, as far as the free energy is concerned, we show that the cell-
theory free-energy at a fixed packing fraction is an increasing function of
chirality angle, implying that the nonchiral phase is the most stable phase.
This fact points to the presence of many-body correlations between particles,
not included in the cell theory, that could explain the presence of chirality.
The entropic self-assembly or clustering of 2D hard particles Glotzer1 ;
Glotzer2 ; Glotzer3 ; Wittmann2  is a physical mechanism to explain the
stability of some liquid-crystal symmetries, not obtained from the particle
shape of the ‘monomeral’ units of clusters MR5 . Motivated by this result, we
also explore the effect of clustering on the stability of the chiral phase by
calculating the free area of hexagonal chiral cluster in the cell formed by
its neighbours located in a hexagonal cell. We show that the free area is
always lower than that of the perfect nonchiral hexagonal clusters, so that
this particular arrangement of particles is not behind the observed chirality.

The article is organized as follows. Section II presents the model and the
analytical expressions for the free energies and equations of state in their
different versions, and make a comparison with MC simulations. Details on the
analytical calculations are relegated to Sec. A. In Sec. III some conclusions
are drawn.

##  II Theoretical Model and Results

Figure 1: Fragments of chiral lattices of HET with different values of
chirality $\theta$ as indicated. Small hexagons from which depart the
triangular vertexes are shown. The large hexagon (dashed line) of size
$\sigma$ contain six triangles from which the packing fraction is obtained.
The angle $\phi$ between the triangular sides and that of the large hexagon is
also indicated.

The model for the crystal phase of HET is a simple 2D lattice where the
centers of mass of triangles of side-length $l$ are located at the nodes of a
triangular lattice. However we extend this model to include chirality, i.e.
the rotation of the main particle axes with respect to the principal lattice
directions. A chiral triangular lattice can be viewed as follows (a schematic
is shown in Fig. 1). The vertexes of neighbour triangles are in contact with
the vertexes of small hexagons. These hexagons are in turn rotated with
respect to the triangle axes (which are fixed and parallel to each other) by
an angle $\theta\in[0,\pi/2]$ (the chiral angle). This procedure is equivalent
to keeping the hexagons fixed and rotating the triangles by the same angle.
Note that the angle $\theta$ is not the one used to characterize chirality in
MC simulations Dijkstra2 . Joining the centers of small hexagons we obtain a
large hexagon (shown with dashed lines) inside which six triangles are
located. The relative angle between the sides of triangles and those of the
large hexagon (see Fig. 1) is called $\phi$; this is the angle used in Ref.
Dijkstra2 .

Figure 2: (a) Fragment of the non-chiral ($\theta=0$) triangular lattice of
HET where the size $s$ of the inscribed into small hexagons equilateral
triangle is indicated. It is also shown the free-area (of size $2s$) of the
caged triangle in the cell formed by its neighbours labeled as 1, 2 and 3.
(b): Fragment of the chiral ($\theta\neq 0$) triangular lattice of HET where
the free area is also shown together with its size $2s\cos\theta$.

First we formulate a very simple model: a cell theory with fixed orientation
for all triangles. We calculate the free area of the center of mass of one
triangle with the same (fixed) orientation as its neighbours inside the cavity
formed by them. The neighbours are located at the nodes of the chiral lattice
(as defined above). See Fig. 2 for the sketch of the free-area, which has the
geometry of an equilateral triangle of side-length $2s\cos\theta$, where $s$
is the length of the (also equilateral) triangle inscribed into the small
hexagons (see Fig. 2). Thus,

$\displaystyle A_{\rm free}=\sqrt{3}s^{2}\cos^{2}\theta.$ (1)

The side length $\sigma$ or the large hexagon as a function of $s$ and
$\theta$ can also be calculated:

$\displaystyle\sigma=\sqrt{l^{2}+s^{2}+2ls\cos\theta},$ (2)

and from here the area of the large hexagon, $\displaystyle{a_{\rm
hex}=\frac{3\sqrt{3}}{2}\sigma^{2}}$. The packing fraction is

$\displaystyle\eta=\frac{6a_{0}}{a_{\rm
hex}}=\left(\frac{l}{\sigma}\right)^{2},$ (3)

which is six times the area of one triangular particle,
$\displaystyle{a_{0}=\frac{\sqrt{3}}{4}l^{2}}$, divided by the area of the
large hexagon. Eqns. (2) and (3) allow to calculate $s$ in units of $l$ as a
function of $\eta$ and $\theta$:

$\displaystyle\tau\equiv\frac{s}{l}=\left[\sqrt{\frac{1}{\eta}-\sin^{2}\theta}-\cos\theta\right],$
(4)

and hence the free area $A_{\rm free}$ from Eqn. (1).

The angle $\phi$ between the edge-length of the large hexagon and the
triangular sides can be calculated as a function of $\theta$ and $\eta$ as

$\displaystyle\cos\phi=\sqrt{\eta}\sin^{2}\theta+\cos\theta\sqrt{1-\eta\sin^{2}\theta},\quad
0\leq\phi\leq\phi_{\rm max}=\arccos\left(\sqrt{\eta}\right).$ (5)

This angle will be used later to plot different quantities and to compare with
MC results. Note that $\phi$ is a function of $\theta$ and $\eta$, and that
$\phi\to 0$ when $\eta\to 0$ regardless of the value of $\theta$. This in turn
means that even when chirality is large ($\theta\sim\pi/3$), $\phi$ could be
rather small if $\eta\sim 1$. Fig. 3 is a 3D plot of the function
$\phi(\theta,\eta)$. Note that, along particular curves $\theta(\eta)$, the
angle $\phi$ could be a nonmonotonic function of $\eta$, exhibiting a maximum.

Figure 3: $\phi$ as a function of $\theta$ and $\eta$. An example is shown of
the path $\theta(\eta)=180\sqrt{\eta-0.75}$, where $\phi$ has a nonmonotonic
behavior with $\eta$.

Finally the cell-theory free energy per particle, using the parallel
orientation approximation, is

$\displaystyle\varphi=-\ln\left(\frac{A_{\rm
free}}{\Lambda^{2}}\right)=-2\log\left(\frac{\sqrt{3}l}{\Lambda}\right)-2\ln\left[\cos\theta\left(\sqrt{\frac{1}{\eta}-\sin^{2}\theta}-\cos\theta\right)\right],$
(6)

with $\Lambda$ the thermal length. From the free energy, the pressure $p^{*}$
as a function of $\theta$ can be calculated. In reduced units,

$\displaystyle p^{*}\equiv\beta
pa_{0}=\eta^{2}\frac{\partial\varphi}{\partial\eta}=\left[\sqrt{\frac{1}{\eta}-\sin^{2}\theta}\left(\sqrt{\frac{1}{\eta}-\sin^{2}\theta}-\cos\theta\right)\right]^{-1}.$
(7)

A much simpler expression results for $p^{*}$ as a function of $\phi$:

$\displaystyle p^{*}=\frac{\eta}{1-\sqrt{\eta}\cos\phi}$ (8)

As shown below, this equation of state performs poorly when compared with MC
simulation results. Moreover the derivatives of $\varphi$ and $p^{*}$ with
respect to $\theta$ are

$\displaystyle\frac{\partial\varphi}{\partial\theta}=2\sin\theta\left(\frac{1}{\cos\theta}-\frac{1}{\sqrt{\eta^{-1}-\sin^{2}\theta}}\right)\geq
0,$ (9) $\displaystyle\frac{\partial
p^{*}}{\partial\theta}=-\frac{\sin\theta}{\left(\eta^{-1}-\sin^{2}\theta\right)^{3/2}}\leq
0.$ (10)

Eqn. (9) in turn implies that chirality increases the free energy. Therefore,
at this level of approximation, cell theory cannot explain the observed
chirality in MC simulations.

Figure 4: Free area $A_{\rm free}/a_{0}$ (in units of the particle area
$a_{0}$) as a function of $\theta$ (a) or $\phi$ (b) (measured in degrees)
from the freely rotating approximation of the caged triangle inside the cell
formed by its neighbours in the chiral lattice (of chirality $\theta$ or
$\phi$). Results for three different values of packing fraction $\eta$ (as
they are labeled) are shown. In the inset of panel (b) we show the derivative
of $A_{\rm free}/a_{0}$ with respect to $\phi$.

The next step is to remove the constraint of a frozen orientation for the
_caged_ triangle (that which wanders inside the cell formed by its
neighbours). The calculation of the free area is now cumbersome, and we
relegate the details to the appendix A. The procedure to calculate the free
area is the following: (i) Fix the angle $\gamma$ between the axis of the
caged triangle and the rest of its neighbours, and move the center of mass of
this particle, keeping contact with the neighbours. This step will define the
envelope of the free area, ${\cal A}(\gamma)$, with a geometry in general more
complicated than equilateral: For certain angular intervals of $\gamma$, the
shape of the free area is made of an equilateral core plus three identical
scalene triangles in contact with the equilateral core; while for other
intervals the are reduces to a single equilateral core. The areas of all of
these triangles can be computed analytically as a function of $l$, $s$,
$\theta$ and $\gamma$ (note that $s$ from Eqn. (4) is a function of $\theta$
and $\eta$). (ii) Analytically calculate the maximum allowed values of the
angle $\gamma$, i.e. $\gamma_{\pm}$, following clockwise ($+$) and
counterclockwise ($-$) rotations, for which the free area collapses to zero,
${\cal A}(\gamma_{\pm})=0$. (iii) The total free area can be computed via
integration, $A_{\rm free}=3\int_{-\gamma_{-}}^{\gamma_{+}}d\gamma{\cal
A}(\gamma)$, which provides an analytical expression. The factor $3$ is
related to the fact that, due to the symmetry of the equilateral triangle, the
condition $\gamma\in[-\gamma_{-},\gamma_{+}]$ reproduces only one third of the
total free area. In Fig. 4 the total free area $A_{\rm free}$ as a function of
$\theta$ or $\phi$, and for different values of the packing fraction $\eta$,
is plotted. Clearly $A_{\rm free}(\theta)$ is a decreasing function of
$\theta$. Fig. 4(b) shows the derivative $\displaystyle{\frac{dA_{\rm
free}(\phi)}{d\phi}}$ vs. $\phi$. It can be seen that the derivative behaves
linearly for angles $\phi\sim 0$, and consequently the function $A_{\rm
free}(\phi)$ decreases quadratically for $\phi\sim 0$. The monotonically
decreasing behavior of $A_{\rm free}$ means that the cell-theory free energy
per particle, $\varphi=-\log\left(A_{\rm free}/\Lambda^{2}\right)$, is an
increasing function of $\theta$ or $\phi$, which confirms that the new version
of cell theory does not explain the stabilization of the chiral crystal.

The total free area $A_{\rm free}$ (see Sec. A) as a function of the small
parameter $\tau$ can be Taylor expanded, which results in

$\displaystyle\frac{A_{\rm
free}}{a_{0}}=\frac{2}{\sqrt{3}}\left(2\tau\cos\theta\right)^{3}+o\left(\tau^{3}\right).$
(11)

Neglecting all terms beyond order $\tau^{3}$, and calculating the pressure
using this approximation, we obtain the expression

$\displaystyle p^{*}=\frac{3}{2}\cdot\frac{\eta}{1-\sqrt{\eta}\cos\phi},$ (12)

which is $3/2$ higher than the pressure obtained from the parallel particle
approximation (see Eqn. (8)), and consequently performs much better when
compared with MC simulations.

In addition, an extension of cell theory to also include the average of the
free area with respect to the angular distribution function
$f_{\rm\theta}(\gamma)$ of the caged particle inside a chiral cell with
chirality angle $\theta$ has been constructed. As usual, the total free energy
has ideal and excess part contributions,

$\displaystyle\varphi\left[f_{\theta}\right]=\int_{-\gamma_{-}}^{\gamma_{+}}d\gamma
f_{\theta}(\gamma)\ln f_{\theta}(\gamma)-1-\ln\frac{\langle{\cal
A}\rangle}{a_{0}},\quad\langle{\cal
A}\rangle\equiv\int_{-\gamma_{-}}^{\gamma_{+}}d\theta f_{\theta}(\gamma){\cal
A}(\gamma),$ (13)

where we drop the trivial term $-\ln\left(a_{0}/\Lambda^{2}\right)$ and
$\gamma_{\pm}>0$ denote the maximum rotation angles in the clockwise ($+$) and
counterclockwise ($-$) directions. Functional minimization of
$\varphi\left[f_{\theta}\right]$ with respect to $f_{\theta}(\gamma)$, with
the constraint $\int_{-\gamma_{-}}^{\gamma_{+}}d\gamma f_{\theta}(\gamma)=1$,
gives

$\displaystyle f_{\theta}(\gamma)=\frac{e^{{\cal A}(\gamma)/\langle{\cal
A}\rangle}}{\int_{-\gamma_{-}}^{\gamma_{+}}d\gamma^{\prime}e^{{\cal
A}(\gamma^{\prime})/\langle{\cal A}\rangle}},$ (14)

while the value of $\langle{\cal A}\rangle$ can be calculated self-
consistently from its definition:

$\displaystyle\langle{\cal
A}\rangle=\frac{\int_{-\gamma_{-}}^{\gamma_{+}}d\gamma{\cal A}(\gamma)e^{{\cal
A}(\gamma)/\langle{\cal
A}\rangle}}{\int_{-\gamma_{-}}^{\gamma^{+}}d\gamma^{\prime}e^{{\cal
A}(\gamma^{\prime})/\langle{\cal A}\rangle}}.$ (15)

At equilibrium we find the free energy

$\displaystyle\varphi\left[f_{\theta}^{(\rm
eq)}\right]=-\ln\left(\int_{-\gamma_{-}}^{\gamma_{+}}d\gamma e^{{\cal
A}(\gamma)/\langle{\cal A}\rangle}\right)-\ln\frac{\langle{\cal
A}\rangle}{a_{0}},$ (16)

and the pressure in reduced units

$\displaystyle p^{*}=-\frac{\eta^{2}}{\langle{\cal
A}\rangle}\left\langle\frac{\partial{\cal A}}{\partial\eta}\right\rangle.$
(17)

Figure 5: (a) Averaged, with respect to $f_{\theta}(\phi)$, free area as a
function of $\phi$ (in degrees) for three different values of packing
fractions (as labeled). (b) Free energy per particle $\varphi$ vs. $\theta$ at
equilibrium for the same values of $\eta$. A constant term
$-\ln\left(3a_{0}/\Lambda^{2}\right)$ was subtracted.

The averaged free area $\langle{\cal A}\rangle$ and free energy per particle
$\varphi$, obtained from Eqn. (16), are plotted in Fig. 5 as a function of
$\phi$ and $\theta$, respectively. The main difference between the values of
$\langle{\cal A}\rangle$ and $A_{\rm free}$ (nonaveraged free area) which can
be seen in Figs. 4 and 5, is related to the fact that $A_{\rm free}$ is not
weighted by the uniform distribution
$\displaystyle{f_{\theta}(\gamma)=\frac{1}{\gamma_{+}+\gamma_{-}}}$. Although
the free energy per particle is an increasing function of $\theta$, and
consequently chiral configurations are not entropically favored within the
present theory, we note that the slope $d\varphi/d\theta$ is very small,
especially close to $\theta=0$. This suggests that any hypothetical mechanism
explaining chirality should involve small entropy gains overcoming the small
entropy loss due to a smaller free area. As an example, many body correlations
could favor the presence of hexagonal clusters formed by six triangles in a
chiral configuration, which crystallize into a hexagonal crystal phase (see
Fig. 6 for a sketch). We explored this avenue and conducted MC integrations to
calculate the free area of chiral hexagonal clusters and their nonchiral
versions (perfect hexagons), crystallizing into hexagonal lattices of the same
packing fraction. We proved that nonchiral perfect hexagons always exhibit
higher free areas than their chiral counterparts. As a consequence, we can
conclude at this point that particle clustering does not seem to be a viable
mechanism to explain the chirality of the crystal phase of HET. We note that
clustering structures other than hexagonal could be explored. However, MC
simulations do not provide any hint as to alternative structures, as only
pictures at very high density are provided in Ref. Dijkstra2 .

Figure 6: Sketch of the hexagonal chiral clusters on a hexagonal lattice.

Figure 7: EOS of the chiral crystal phase of the fluid of HET from cell
theories with chirality fixed to $\theta=0^{\circ}$ (a), $30^{\circ}$ (b),
$45^{\circ}$ (c) and $60^{\circ}$ (d). With symbols we show MC simulation
results. The dotted line indicates the simplest cell theory model with frozen
orientation of the caged particle while with solid and dashed lines (both are
very close to each other) we show the results from the cell theory of freely
rotating caged triangle with uniform and minimized distribution functions
$f_{\theta}(\gamma)$ respectively. Finally with dot-dashed line we show the
EOS from Eq. (8). For the values of packing fraction in the range [0.70,0.89]
MC simulations predict an stable 6-atic liquid-crystal phase with the
transition to the crystal phase at $\eta=0.89$ (labeled with arrow) Dijkstra2
. Thus our cell theory describe very well also the EOS of the liquid-crystal
phase for $\eta>0.75$ (not too close from the I–6-atic transition).

Despite this failure, the two theories presented (namely freely rotating caged
particle with uniform orientational distribution or else obtained via
functional minimization) enjoy an important success: their fair performance in
describing the equation of state (EOS) obtained from MC simulations. Fig. 7
compares both theories, together with that derived from the frozen orientation
approximation (see Eqn. (6)). The free-rotation approximation accurately
describes the EOS in an interval of packing fraction from 0.75 to 0.9, while
the frozen-orientation approach provides a poor representation. In fact is is
possible to optimize the value of $\theta$ (or $\phi$) for which the
theoretical description is better. Fig. 7 allows to discard values
$\theta>45^{\circ}$, while values $\theta\sim 30^{\circ}$ (or $\phi\sim 3.44$,
2.55 and 1.68 for $\eta=0.8$, 0.85 and 0.9 respectively) seem to optimize the
agreement with simulations Dijkstra2 . Actually, as simulations show, in the
liquid-crystal interval of packing fractions, the chirality of the fluid is
completely absent. Chirality emerges beyond a particular packing fraction in
the crystal phase. From Fig. 7 we can infer that, from $0.75$ up to $0.87$,
the best fitting from cell theory is obtained for $\theta=0$, while beyond
$0.87$ a chirality angle between $30^{\circ}$ and $45^{\circ}$ improves the
description of the EOS in the crystal phase. Unfortunately, the dependence of
chirality on packing fraction cannot be obtained in the context of our theory,
due to the monotonically increasing behavior of the free energy as a function
of $\theta$.

##  III Discussion and conclusions

We have proposed three versions of cell theory for the chiral crystal phase of
HET. The simplest theory involves the constraining of the caged particle
inside a cell defined by the fixed positions of its neighbours, located at the
nodes of the chiral triangular lattice, and with their axes parallel to each
other. Moreover the orientation of the caged particle is identical to those of
its neighbours. We have shown that this version fails to predict the EOS of
the liquid crystal and crystal phases obtained from MC simulations. However,
when multiplied by a factor of $3/2$, whose origin comes from the leading
order in a small lattice-parameter expansion of the free area for a freely-
rotating caged particle, the performance in the description of the EOS greatly
improves.

We derived a second version of the theory by removing the constraint of frozen
orientation of the caged particle (keeping the orientation of its neighbours
frozen) and allowing it to rotate inside the cell with a uniform angular
distribution function. We have calculated analytically the integral of the
free area over the angular degrees of freedom and the associated free energy.
The ensuing pressure considerably improves the EOS when compared with MC
results pertaining to the EOS of the 6-atic liquid-crystal and crystalline
phases. By optimizing the chiral angle one can improve the agreement with MC
simulations, the conclusion being that the chiral angle cannot be larger than
$45^{\circ}$, also in agreement with simulations.

The third and last version of the theory involves the averaging of the free
area over the angular distribution function. This function is obtained from
minimization of a proposed free-energy density functional consisting of ideal
angular term plus excess part, approximated as minus the logarithm of the
averaged free area. The latter is obtained self-consistently using a numerical
procedure. The EOS obtained from this version is practically indistinguishable
from that obtained with the previous theory.

We have shown that the free energy per particle from the three theories are
always increasing functions of the chiral angle, implying that the most stable
configuration of the particles are achiral. However, a more detailed study of
the integrated free area with respect to clockwise and counterclockwise
rotations (assuming a clockwise chirality), shows that the former exhibits a
maximum at a very small chiral angle. This fact, and the small value of the
slope of the free-energy as a function of the chiral angle, makes it plausible
the existence of another mechanism, based most likely on many-body particle
correlations, with an associated increase in entropy, sufficient to overcome
the entropy loss involved in a chiral rotation. Exploring this avenue will
have to wait until a density-functional theory for inhomogeneous densities,
incorporating exact correlations beyond two particles, can be developed.
Looking for a possible mechanism to explain chirality we proved that, in the
context of cell theory, there is no free-area gain when particles are
clustered into perfect chiral hexamers, as compared with their nonchiral
counterparts: this particular clustering effect is not behind the chirality of
the crystal phase observed in MC simulations.

##  Appendix A Free area of the freely rotating caged triangle

Figure 8: (a) Configuration where the caged triangle is in contact with two of
its neighbours (with their closest vertexes labeled with filled circles). The
known quantities are the angles $\alpha=2\pi/3+\theta$ and $\gamma$ and the
lengths $l$ and $s$ while the unknowns are $x^{\prime}$, $\beta$, $x$ and $y$.
(b) The geometry of the free area for relatively small values of $\gamma$: an
equilateral triangle of size $t$ and three equal scalene triangles of sizes
$x$, $y$ and $t$ with the angle $\pi/3-\gamma$ being the interior angle
between the sides of lengths $x$ and $y$.

In this section we show how to obtain the analytic expressions for the free
area of a freely rotating caged triangle inside the cell formed by its
neighbours in the chiral triangular lattice configuration. In Fig. 8(a) we
show the caged triangle rotated counterclockwise in a configuration where it
is in contact with two of its neighbours (two of their vertexes labeled with
filled circles). The known quantities are the lengths $l$, $s$, and the angles
$\gamma>0$ (the rotation angle) and $\alpha=2\pi/3+\theta$. The unknown
quantities, necessary to find the free area, are $x$ and $y$. From the sinus-
law of the triangle construction, shown in Fig. 8(a), we can find $x$ and $y$
as

$\displaystyle
l-x=\frac{\sqrt{3}l/2-s\sin\left(\pi/3+\gamma+\theta\right)}{\sin\left(\pi/3-\gamma\right)},$
(18) $\displaystyle
y=\frac{s\sin(\pi/3-\theta)-l\sin\gamma}{\sin\left(\pi/3-\gamma\right)}.$ (19)

In Fig. 8(b) a sketch is shown of the free area for relatively small angles
$\gamma$. It consists of an equilateral triangular core of length $t$ and
three identical scalene triangles of sizes $x$, $y$ and $t$, with the interior
angle between the sides of sizes $x$ and $y$ being equal to $\pi/3-\gamma$.
The free area for this case can be computed as

$\displaystyle{\cal A}^{({\rm
L},1)}(\gamma)=\frac{\sqrt{3}}{4}\left(x^{2}+y^{2}+4xy\cos(\pi/3+\gamma)\right),$
(20)

which, after substituting Eqns. (18) and (19), becomes (in units of $a_{0}$)

$\displaystyle\frac{{\cal A}^{({\rm
L},1)}(\gamma)}{a_{0}}=\frac{1}{2}+2\tau\left[\cos(\psi)+\sqrt{3}\sin(\psi)\right]+\tau^{2}\left[1+\sqrt{3}\sin(2\psi)\right]$
$\displaystyle-4\left[1+\tau\cos(\psi)\right]\cos(\gamma^{\prime})+4\tau\sin(\psi)\sin(\gamma^{\prime})-\frac{6\tau\sin(\psi)}{\sin(\gamma^{\prime})}$
$\displaystyle+\sqrt{3}\left[\frac{3}{2}+2\tau^{2}\sin^{2}(\psi)\right]\cot(\gamma^{\prime}),\quad
0\leq\gamma\leq\gamma_{-}^{(1)}$ (21)

where we have used the superindexes L and 1 to indicate that it is the first
(1) part of the free area when the caged triangle is rotated to the left (L),
i.e. counterclockwise. In (21) we have used the notations
$\tau=s/l=\sqrt{\eta^{-1}-\sin^{2}\theta}-\cos\theta$, $\psi=\pi/3-\theta$,
and $\gamma^{\prime}=\pi/3-\gamma$, while
$\gamma_{-}^{(1)}=\arcsin(\tau\sin\psi)$ is the angle for which the length $y$
becomes zero and the free area collapses to an equilateral triangle. Note that
$\gamma_{-}^{(1)}=0$ for $\theta=\pi/3$. Integrating (21) with respect to
$\gamma$ from 0 to $\gamma_{-}^{(1)}$ we obtain

$\displaystyle\frac{A_{\rm free}^{({\rm
L},1)}}{a_{0}}\equiv\frac{\int_{0}^{\gamma_{-}^{(1)}}d\gamma{\cal A}^{({\rm
L},1)}(\gamma)}{a_{0}}$
$\displaystyle=\left[\frac{1}{2}+2\tau\left(\cos\psi+\sqrt{3}\sin\psi\right)+\tau^{2}\left(1+\sqrt{3}\sin(2\psi)\right)\right]\gamma_{-}^{(1)}$
$\displaystyle+4(1+\tau\cos\psi)\left(\sin\delta_{-}-\frac{\sqrt{3}}{2}\right)+4\tau\sin\psi\left(\cos\delta_{-}-\frac{1}{2}\right)$
$\displaystyle+6\tau\sin\psi\ln\left[\sqrt{3}\tan(\delta_{-}/2)\right]+\sqrt{3}\left(\frac{3}{2}+2\tau^{2}\sin^{2}\psi\right)\ln\left[\frac{\sqrt{3}}{2}\csc\delta_{-}\right],$
(22)

with $\delta_{-}=\pi/3-\gamma_{-}^{(1)}$.

Figure 9: The geometric construction to find the variables $x$ and $y$,
necessary to calculate the free area for $\gamma>\gamma_{-}^{(1)}$. Filled
circles indicate the vertexes of the neighbours to the caged triangles (not
shown). The known quantities are $l$, $s$, $\alpha=2\pi/3+\theta$ and
$\gamma$, while the unknowns are $\beta$, $x^{\prime}$, $x$ and $y$.

For $\gamma>\gamma_{-}^{(1)}$ the main geometric construction to find the
variables $x$ and $y$ is sketched in Fig. 9. This time the free area is an
equilateral triangle of size $x-y$, so the free area is ${\cal A}^{({\rm
L},2)}(\gamma)=\sqrt{3}\left(x-y\right)^{2}/4$. From Fig. 9 we find:

$\displaystyle
l-x=\frac{2}{\sqrt{3}}\left[l\sin\left(\pi/3+\gamma\right)-s\sin\left(\pi/3+\theta\right)\right],$
(23) $\displaystyle
y=\frac{2}{\sqrt{3}}\left[l\sin\gamma-s\sin(\pi/3-\theta)\right],$ (24)

Thus we have

$\displaystyle\frac{{\cal A}^{({\rm
L},2)}(\gamma)}{a_{0}}=3\tau+4\tau\cos\theta+4\tau^{2}\cos^{2}\theta-4(1+2\tau\cos\theta)\cos(\pi/3-\gamma)-2\cos(\pi/3+2\gamma),$
(25)

valid for $\gamma_{-}^{(1)}\leq\gamma\leq\gamma_{-}^{(2)}$, where
$\gamma_{-}^{(2)}$ is the angle for which $x-y=0$:

$\displaystyle\gamma_{-}^{(2)}=\frac{\pi}{3}-\arccos\left(\frac{1}{2}+\tau\cos\theta\right).$
(26)

The integration of (25) with respect to $\gamma$ from $\tilde{\gamma}_{-}$ to
$\gamma_{-}^{(2)}$, where $\tilde{\gamma}_{-}=\gamma_{-}^{(1)}$ for
$\theta\leq\pi/3$ while $\tilde{\gamma}_{-}=0$ for $\theta>\pi/3$, gives

$\displaystyle\frac{A_{\rm free}^{({\rm
L},2)}}{a_{0}}=\frac{\int_{\tilde{\gamma}_{-}}^{\gamma_{-}^{(2)}}d\gamma{\cal
A}^{({\rm
L},2)}(\gamma)}{a_{0}}=\left(3+4\tau\cos\theta+4\tau^{2}\cos^{2}\theta\right)\left(\gamma_{-}^{(2)}-\tilde{\gamma}_{-}\right)$
$\displaystyle+4(1+2\tau\cos\theta)\left(\sin\left(\frac{\pi}{3}-\gamma_{-}^{(2)}\right)-\sin\left(\frac{\pi}{3}-\tilde{\gamma}_{-}\right)\right)+\sin\left(\frac{\pi}{3}+2\tilde{\gamma}_{-}\right)-\sin\left(\frac{\pi}{3}+2\gamma_{-}^{(2)}\right).$

Thus the total integrated free area for counterclockwise rotation can be found
from

$\displaystyle A_{\rm free}^{({\rm L})}=\left\\{\begin{matrix}3\left(A_{{\rm
free}}^{({\rm L},1)}+A_{\rm free}^{({\rm L},2)}\right),&\text{if}\
\displaystyle{0\leq\theta\leq\frac{\pi}{3}},\\\ 3A_{\rm free}^{({\rm
L},2)},&\text{if}\
\displaystyle{\frac{\pi}{3}\leq\theta\leq\frac{\pi}{2}.}\end{matrix}\right.$
(28)

where the factor 3 is due to the presence of three sectors of equal areas due
to the symmetry of the cell geometry.

Figure 10: The left free area $A_{\rm free}^{({\rm L})}/a_{0}$ as a function
of $\phi$ measured in degrees for three different values of $\eta$ as they are
labeled.

Fig. 10 plots the function $A_{\rm free}^{({\rm L})}(\theta)$ for three
different values of $\eta$. Clearly it is a decreasing function of $\theta$.

When we rotate the caged triangle clockwise (in the same direction of
chirality), the free areas are different. For small enough angle $\gamma$, the
triangular constructions to find the magnitudes $x$ and $y$ are similar to
that given in Fig. 8(a) and Fig. 9, but with an important difference. Now we
have $\alpha=2\pi/3-\theta$ (instead of $2\pi/3+\theta$), so the magnitudes
$x$ and $y$ are given by the same Eqns. (18)-(19) and (23)-(24) except for the
substitution $\theta\to-\theta$. Thus we have that the free areas have the
same expressions as in Eqns. (21) and (25) with this substitution, i.e.

$\displaystyle{\cal A}^{({\rm R},1)}(\gamma;\theta)={\cal A}^{({\rm
L},1)}(\gamma;-\theta),\ \text{for}\ 0\leq\gamma\leq\gamma_{+}^{(1)}$ (29)
$\displaystyle{\cal A}^{({\rm R},2)}(\gamma;\theta)={\cal A}^{({\rm
L},2)}(\gamma;-\theta),\ \text{for}\
\gamma_{+}^{(1)}\leq\gamma\leq\gamma_{+}^{(2)},$ (30)

where we have

$\displaystyle\gamma_{+}^{(i)}(\theta)=\gamma_{-}^{(i)}(-\theta),\quad i=1,2.$
(31)

Figure 11: Triangle constructions from which the magnitudes $y_{1}$ and
$y_{2}$, and consequently the excluded area, can be found for
$\theta>\theta^{*}$. The lengths $l$ and $s$ and the angles $\alpha$ and
$\gamma$ are known while the lengths $x^{\prime}$, $y_{1}$ and $y_{2}$ and the
angle $\delta$ are unknown.

But now the length $x$ can be equal to zero (which in turn implies that the
free area becomes a equilateral triangle of side $y$), and this takes place at
an angle

$\displaystyle\gamma_{+}^{(3)}=\arcsin\left[\frac{\sqrt{3}/2+\tau\sin(\pi/3-\theta)}{{\cal
L}}\right]-\arcsin\left[\frac{\sqrt{3}/2}{{\cal L}}\right].$ (32)

where we have defined the quantity

$\displaystyle{\cal L}=\sqrt{1+2\tau\cos(\pi/3+\theta)+\tau^{2}}.$ (33)

The equality $\gamma_{+}^{(1)}=\gamma_{+}^{(2)}=\gamma_{+}^{(3)}$ takes place
at
$\displaystyle{\theta=\theta^{*}=\arcsin\left(\frac{1}{2\sqrt{\eta}}\right)}$.
Thus, for $\theta\leq\theta^{*}$, the integrated free areas $A_{\rm
free}^{({\rm R},i)}$ ($i=1,2$) has the same expressions as those from Eqns.
(22) and (A), except for the substitutions $\theta\to-\theta$ and
$\gamma_{-}^{(i)}\to\gamma_{+}^{(i)}$. For $\theta>\theta^{*}$ and
$0\leq\gamma\leq\gamma_{+}^{(3)}$, we have again the equality (29) with the
substitution $\gamma_{+}^{(1)}\to\gamma_{+}^{(3)}$, necessary to calculate the
integrated area $A_{\rm free}^{({\rm R},1)}$ in this interval of angles.
However, for $\gamma>\gamma_{+}^{(3)}$, we obtain that the free area is an
equilateral triangle of length $y_{1}-y_{2}$, where the magnitudes $y_{i}$ can
be calculated from the triangular construction shown in Fig. 11, with the
result

$\displaystyle
y_{1}=l\left[1-\frac{2}{\sqrt{3}}\sin(\pi/3-\gamma)\right]-\frac{2}{\sqrt{3}}s\sin(\pi/3+\gamma-\theta),$
(34) $\displaystyle
y_{2}=\frac{2}{\sqrt{3}}\left[s\sin(\pi/3+\theta-\gamma)-l\sin\gamma\right].$
(35)

Thus we obtain the free area

$\displaystyle\frac{{\cal A}^{({\rm
R},2)}(\gamma)}{a_{0}}=3+4\tau\cos(\pi/3+\theta)+2\tau^{2}$
$\displaystyle+\left[4\tau\cos(\pi/3-\theta)+2\tau^{2}\cos(2\theta)-1\right]\cos(2\gamma)$
$\displaystyle-2(1+2\tau\cos\theta)\cos\gamma+2\left(\sqrt{3}-2\tau\sin\theta\right)\sin\gamma$
$\displaystyle-2\left[\frac{\sqrt{3}}{2}+2\tau\sin(\pi/3-\theta)-\tau^{2}\sin(2\theta)\right]\sin(2\gamma),$
(36)

which has the primitive

$\displaystyle{\cal S}(\gamma)=\int d\gamma\frac{{\cal A}^{({\rm
R},2)}(\gamma)}{a_{0}}=3\gamma+\sin(\pi/3-2\gamma)-4\sin(\pi/3+\gamma)$
$\displaystyle+2\tau\left[2\gamma\cos(\pi/3+\theta)+\sin(\pi/3-\theta+2\gamma)+2\sin(\theta-\gamma)\right]$
$\displaystyle+\tau^{2}\left[2\gamma-\sin(2(\theta-\gamma))\right].$ (37)

The excluded area (36) becomes zero when $y_{1}=y_{2}$, a condition that takes
place at the angle

$\displaystyle\gamma_{+}^{(4)}=\arccos\left(\frac{1/2}{{\cal
L}}\right)-\arccos\left(\frac{1/2+\tau\cos\theta}{{\cal L}}\right).$ (38)

Note that $\gamma_{+}^{(3)}=0$ for $\theta=\pi/3$ while $\gamma_{+}^{(4)}=0$
for $\theta=\pi/2$. Then we have that the second integrated area for
$\theta>\theta^{*}$ can be calculated from

$\displaystyle\frac{A_{\rm free}^{({\rm
R},2)}}{a_{0}}=\left\\{\begin{matrix}{\cal
S}\left(\gamma_{+}^{(4)}\right)-{\cal
S}\left(\gamma_{+}^{(3)}\right),&\displaystyle{\theta^{*}\leq\theta\leq\frac{\pi}{3}},\\\
{\cal S}\left(\gamma_{+}^{(4)}\right)-{\cal
S}(0),&\displaystyle{\frac{\pi}{3}\leq\theta\leq\frac{\pi}{2}}.\end{matrix}\right.$
(39)

Figure 12: The right free area as a function of $\phi$ (in degrees) for three
different values of packing fractions as they are labeled.

Finally, the total area for clockwise rotation is given again by the same Eqn.
(28), and it is plotted in Fig. 12 as a function of $\phi$ for three different
values of $\eta$. Note that they exhibit small maxima at angles close, but not
equal, to $0$. However the total free area, $A_{\rm free}=A_{\rm free}^{({\rm
L})}+A_{\rm free}^{({\rm R})}$, exhibits a monotonic behavior, as seen in Fig.
4. Therefore, the presence of these maxima are not enough to generate a non-
monotonic behavior of the total free-area (which would explain the presence of
chirality in MC simulations).

Figure 13: (a) The functions ${\cal A}^{({\rm L})}(\gamma)$ (solid) and ${\cal
A}^{({\rm R})}(\gamma)$ (dashed) for different values of $\theta$ as they
labeled and for a fixed value of $\eta=0.8$. With different symbols we show
that maximum values of $\gamma$ for each curve. (b) The equilibrium angular
distribution functions $f_{\theta}(\gamma)$ for the same values of $\theta$
and $\eta$.

The reason for the nonmonotonic behavior of $A_{\rm free}^{({\rm R})}$ as a
function of $\theta$ or $\phi$ can be easily explained by examining the
functions ${\cal A}^{(\nu)}(\gamma)={\cal A}^{(\nu,1)}(\gamma)+{\cal
A}^{(\nu,2)}(\gamma)$ ($\nu=\left\\{L,R\right\\}$) for different values of
$\theta$, see Fig. 13 (a). It can be seen that, for certain values of
$\theta_{2}>\theta_{1}$, the inequality ${\cal A}^{({\rm
R})}_{\theta_{2}}(\gamma)>\left\\{{\cal A}^{({\rm
L})}_{\theta_{1}}(\gamma),{\cal A}^{({\rm R})}_{\theta_{1}}(\gamma)\right\\}$
is fulfilled for high enough values of $\gamma$. Thus, there exist values of
$\theta_{1},\theta_{2}\sim 0$ for which the inequality between the integrals
$A_{{\rm free},\theta_{2}}^{({\rm R})}=3\int d\gamma{\cal A}^{({\rm
R})}_{\theta_{2}}(\gamma)>A_{{\rm free},\theta_{1}}^{({\rm R})}=3\int
d\gamma{\cal A}^{({\rm R})}_{\theta_{1}}(\gamma)$ is also fulfilled which, in
turn, explains the presence of maxima of $A_{\rm free}^{({\rm R})}$ as a
function of $\theta$ or $\phi$ at a vanishing angle. Fig. 13(b) is a plot of
the equilibrium angular distribution function $f_{\theta}(\gamma)$ calculated
from Eqn. (14) with ${\cal A}(\gamma)={\cal A}^{({\rm L})}(|\gamma|)$ if
$\gamma<0$, and ${\cal A}(\gamma)={\cal A}^{({\rm R})}(|\gamma|)$ if
$\gamma\geq 0$. It can be seen that it is not symmetric with respect to
$\gamma=0$. This in turn implies that the value of
$\langle\gamma\rangle\equiv\int_{-\gamma_{-}}^{\gamma_{+}}d\gamma\gamma
f_{\theta}(\gamma)\neq 0$. We obtain $\langle\gamma\rangle=0.7130^{\circ}$,
$0.3815^{\circ}$, $0.1561^{\circ}$ and $0.0405^{\circ}$ for
$\theta=75^{\circ}$, $60^{\circ}$, $45^{\circ}$ and $30^{\circ}$,
respectively.

###### Acknowledgements.

Financial support from Grant No. PID2021-126307NB-C21 (MCIU/AEI/FEDER,UE) is
acknowledged.

## References

  * (1) K. Zhao, C. Harrison, D. Huse, W. B. Russel, and P. M. Chaikin, Phys. Rev. E 76, 040401(R) (2007).
  * (2) K. Zhao, R. Bruinsma, and T. G. Mason, Proc. Natl. Acad. Sci. USA 108, 2684 (2011).
  * (3) Z. Hou, K. Zhao, Y. Zong, and T. G. Mason, Phys. Rev. Materials 3, 015601 (2019).
  * (4) Z. Hou, Y. Zong, Z. Sun, F. Ye, T. G. Mason, and K. Zhao, Nat. Commun. 11, 2064 (2020).
  * (5) H. Schlacken, H.-J. Mogel, and P. Schiller, Mol. Phys. 93, 777 (1998).
  * (6) M. A. Bates and D. Frenkel, J. Chem. Phys. 112, 10034 (2000).
  * (7) M. C. Lagomarsino, M. Dogterom, and M. Dijkstra, J. Chem. Phys. 119, 3535 (2003).
  * (8) _Computational Methods in Science and Technology_ , ed. by K. W. Wojciechowski and D. Frenkel (Science Publishers OWN, Poznan, Poland, 2004), 10, 235.
  * (9) Y. Martínez-Ratón, E. Velasco, and L. Mederos, J. Chem. Phys. 122, 064903 (2005).
  * (10) A. Donev, J. Burton, F. H. Stillinger, and S. Torquato, Phys. Rev. B 7J. Chem. Phys. 140, 204502 (2014)3, 054109 (2006).
  * (11) Y. Martínez-Ratón, E. Velasco, and L. Mederos, J. Chem. Phys. 125, 014501 (2006).
  * (12) C. Avendaño and F. A. Escobedo, Soft Matter 8, 4675 (2012).
  * (13) J. A. Martínez-González, S. Varga, P. Gurin and J. Quintana-H, J. Mol. Liq. 185, 26 (2013).
  * (14) S. Mizani, P. Gurin, R. Aliabadi, H. Salehi, and S. Varga, J. Chem. Phys. 153, 034501 (2020).
  * (15) G. Bautista-Carbajal and G. Odriozola, J. hem. Phys. 140, 204502 (2014).
  * (16) T. Geigenfeind and D. de las Heras, J. Chem. Phys. 150, 184906 (2019).
  * (17) R. Wittmann, C. E. Sitta, and F. Smallenburg, and H. Löwen, J. Chem. Phys. 147, 134908 (2017).
  * (18) A. P. Gantapara, W. Qi, and M. Dijkstra, Soft Matter 11, 8684 (2015).
  * (19) Y. Martínez-Ratón, A. Díaz-De Armas, and E. Velasco, Phys. Rev. E 97, 052703 (2018).
  * (20) Y. Martínez-Ratón and E. Velasco, Phys. Fluids 34, 037110 (2022).
  * (21) Y. Martínez-Ratón and E. Velasco, J. Mol. Liq. 397, 124044 (2024).
  * (22) K. Zhao, R. Bruisma, and T. G. Mason, Nat. Commun. 3, 801 (2012).
  * (23) K. Zhao and T. G. Mason, J. Phys.: Condens. Matter 26 152101 (2014).
  * (24) K. Mayoral and T. G. Mason, Soft Matter 10, 4471 (2014).
  * (25) S. P. Carmichael and M. S. Shell, J. Chem. Phys 139, 164705 (2013).
  * (26) T. Vo and S. C. Glotzer, Proc. Natl. Acad. Sci. USA 119, 211641119 (2022).
  * (27) S. Lee, T. Vo, and S. C. Glotzer, Nat. Chem. 15, 905 (2023).
  * (28) T. Dwyer, T. C. Moore, J. A. Anderson, and S. C. Glotzer, Soft Matter 36, 711 (2023).
  * (29) R. Wittmann, P.A. Monderkamp, and H. Löwen, Phys. Rev. E 107, 064602 (2023).