# A confined rod: mean field theory for hard rod-like particles

Jamie M. Taylor 1, Thomas G. Fai 2, Epifanio G. Virga3,
Xiaoyu Zheng 4, Peter Palffy-Muhoray 4,5
_1 CUNEF Universidad, Madrid, Spain_
_2 Department of Mathematics and Volen Center for Complex Systems,_
_Brandeis University Waltham, MA, USA_
_3 Dipartimento di Matematica, Università di Pavia, Pavia, Italy _
_4 Department of Mathematical Sciences, Kent State University, Kent, OH, USA _
_5 Advanced Materials and Liquid Crystal Institute, Kent State University,
Kent, OH, USA_ 0000-0002-5423-828X<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper, we model the configurations of a system of hard rods by viewing
each rod in a cell formed by its neighbors. By minimizing the free energy in
the model and performing molecular dynamics, where, in both cases, the shape
of the cell is a free parameter, we obtain the equilibrium orientational order
parameter, free energy and pressure of the system. Our model enables the
calculation of anisotropic stresses exerted on the walls of the cell due to
shape change of the rod in photoisomerization. These results are a key step
towards understanding molecular shape change effects in photomechanical
systems under illumination.

##  1 Introduction

Our motivation for this work is to gain insights into photomechanical
materials which convert light energy directly into mechanical work. Azo-dye
containing nematic liquid crystal elastomers (LCEs) are one example of such
materials [1, 2, 3, 4]. Under illumination, dye molecules can undergo
photoisomerization, changing their shape, resulting in the macroscopic
deformation of samples, which can exert forces and do mechanical work.

Here we aim to understand the mechanics of bulk shape change of such materials
from a microscopic perspective. Specifically, we consider a system of hard
rod-like particles which interact via hard-core steric interactions,
conserving linear and angular momentum and energy. At low number densities, we
expect the system to be isotropic, but as the number density is increased, a
transition to an orientationally ordered phase is expected. Unlike in Onsager
theory [5], where the pair-excluded volume plays an essential role, we regard
each rod as being effectively confined in a local cell – a rectangular box or
an ellipsoid – created by neighboring particles, where the free volume of each
particle in its cell determines the phase behavior.

There are various cell models for hard spheres [6, 7]. In cell theories, the
free volume for non-interacting rigid spheres is defined as that volume of the
cell in which the center of mass of a particular molecule is able to move when
all of the other molecules in the liquid or dense gas are held fixed at their
mean lattice positions. The single occupancy model, with one particle per
cell, works well at high densities, but less well at low densities. To rectify
the situation, the concept of ‘communal entropy’ was introduced to describe
the solid-fluid phase transition [8].

As a first attempt for rod-like particles, in this paper, we shall only
consider the single occupancy model where each rod is confined to its private
cell. Rather than considering the precise geometry of the boundary of the cell
described by the neighboring particles, we assume that the cell can be
approximated via a geometric mean field, formed from a simple geometry
(cuboidal or ellipsoidal) of given volume. We anticipate our results to
produce quantitatively accurate thermodynamic properties, such as
orientational order, pressure and energy in the high density regime.

The paper is organized as follows. In section 2, we lay out our mean-field
theory, where each particle is regarded as being confined to a cell defined by
its neighbors. The free volume of a particle in a given cell depends on its
orientation; the free energy depends on the average free volume over all
possible particle orientations available to the center of the particle in its
cell. At a given volume fraction, the cell adopts a shape which maximizes the
average free volume, or, equivalently, minimizes the free energy. In section
3, we apply the theory to the case of an ellipsoid in a rectangular cell, and
compare the results with those from molecular dynamics simulations. We also
consider the limiting case of a needle in an ellipsoidal cell. In section 4,
we discuss the effects of shape change in both rectangular and ellipsoidal
cells. The paper is summarized in the conclusion section 5.

##  2 Mean-field theory

The Helmholtz free energy for a system of $N$ indistinguishable particles in
volume $V$ at temperature $T$ is given by [7]

$F=-kT\ln Z_{tot},$ (1)

where $k$ is Boltzmann’s constant and $Z_{tot}$ is the total partition
function. After integrating over momenta, this is given by

$Z_{tot}=\frac{1}{\Lambda_{dB}^{3N}N!}\int
e^{-\frac{H(\mathbf{q}_{1},...,\mathbf{q}_{N})}{kT}}d\mathbf{q}_{1}...d\mathbf{q}_{N},$
(2)

where $\Lambda_{dB}$ is the de Broglie wavelength,
$\mathbf{q}_{i}=(\mathbf{r}_{i},\mathbf{\hat{m}}_{i})$ is the generalized
coordinate of the $i-$th particle, where $\mathbf{r}_{i}$ is the position of
center of mass and $\mathbf{\hat{m}}_{i}$ is a unit vector along the symmetry
axis of the particle.

The essence of mean field theory is to approximate the potential energy of the
system by the sum of identical single-particle potentials, which describe the
interaction of each particle with an effective mean field [9]. Single-particle
potentials are not unique, but they must satisfy certain consistency
conditions [10]. In this work, we assume that there are no long-range,
attractive interactions and the particles are kept in a prescribed volume $V$
by an adjustable external pressure. Since the interactions are steric, the
particles cannot overlap. As a first attempt, then, following the single-
particle potential construction procedure in [10], we obtain the single
particle potential,

$U(\mathbf{q)=}\left\\{\begin{array}[]{cc}0,&\text{if particle is wholly
within the cell, and}\\\ \infty&\text{if it is not,}\end{array}\right.$ (3)

which satisfies the consistency requirements. Each particle is thus confined
to a cell $\mathcal{C}$ defined by the surrounding particles. Although in a
real system the cells have non-uniform distributions and fluctuate in time,
here, in order to make the model tractable, we assume that the cells are
frozen in time and all of equal size.

The total partition function is

$Z_{tot}=z^{N}=\left(\frac{1}{\Lambda_{dB}^{3}}\int
e^{-\frac{U(\mathbf{q}_{1})}{kT}}d\mathbf{q}_{1}\right)^{N},$ (4)

where $z$ is the single-particle partition function, and $U(\mathbf{q}_{1})$
is the single-particle potential in Eq. (3). Being now in individual cells,
all particles are effectively distinguishable and there is no need for the
Boltzmann factor $1/N!$ as in Eq. (2).

The integration over position $\mathbf{r}_{1}$ in Eq. (4) gives the free
volume $V_{f}(\mathbf{\hat{m},}\mathcal{C})$, the volume accessible to the
center of the particle in the cell $\mathcal{C}$ for a given particle
orientation $\mathbf{\hat{m}}$. That is,

$z=\frac{1}{\Lambda_{dB}^{3}}\int
e^{-\frac{U(\mathbf{q}_{1})}{kT}}d\mathbf{q}_{1}=\frac{1}{\Lambda_{dB}^{3}}\int_{S^{2}}\int_{\mathcal{C}}e^{-\frac{U(\mathbf{q}_{1})}{kT}}d\mathbf{r}_{1}d\mathbf{m}_{1}=\frac{1}{\Lambda_{dB}^{3}}\int_{S^{2}}V_{f}(\mathbf{\hat{m}}_{1},\mathcal{C})d\mathbf{\hat{m}}_{1}.$
(5)

The total partition function in Eq. (4) can then be rewritten as

$Z_{tot}(\mathcal{C)}=\left(\frac{1}{\Lambda_{dB}^{3}}\int_{S^{2}}V_{f}(\mathbf{\hat{m},}\mathcal{C})d\mathbf{\hat{m}}\right)^{N}.$
(6)

The free energy in Eq. (1) can be approximated, to within an inessential
additive constant, as

$F=-NkT\ln\frac{1}{\Lambda_{dB}^{3}}\int_{S^{2}}V_{f}(\mathbf{\hat{m},}\mathcal{C})d\mathbf{\hat{m}}.$
(7)

The free energy density $\mathcal{F}=F/(NV_{cell})$ may be written in the more
familiar form

$\mathcal{F}=kT\rho_{0}\ln\rho_{0}-kT\rho_{0}\ln\frac{1}{\Lambda_{dB}^{3}}\frac{1}{V_{cell}}\int_{S^{2}}V_{f}(\mathbf{\hat{m}},\mathcal{C})d\mathbf{\hat{m}},$
(8)

where $\rho_{0}=1/V_{cell}$ is the number density, but for our purposes Eq.
(7) suffices.

The orientational probability density function $\rho(\mathbf{\hat{m})}$ is

$\rho(\mathbf{\hat{m}})=\frac{V_{f}(\mathbf{\hat{m}},\mathcal{C})}{\int_{S^{2}}V_{f}(\mathbf{\hat{m}},\mathcal{C})d\mathbf{\hat{m}}},$
(9)

and the orientational order parameter tensor $\mathbf{Q}$ is

$\mathbf{Q}=\frac{\int_{S^{2}}\frac{1}{2}(3\mathbf{\hat{m}\hat{m}}-\mathbb{I})V_{f}(\mathbf{\hat{m},}\mathcal{C})d\mathbf{\hat{m}}}{\int_{S^{2}}V_{f}(\mathbf{\hat{m},}\mathcal{C})d\mathbf{\hat{m}}},$
(10)

where is $\mathbb{I}$ the identity tensor. The scalar order parameters are
projections of $\mathbf{Q}$ onto the principal directions of the cell
$\mathbf{\hat{e}}_{i}\mathbf{\hat{e}}_{i}$, and

$S_{2i}=\left\langle\mathcal{P}_{2}(\mathbf{\hat{m}\cdot\hat{e}}_{i})\right\rangle,$
(11)

where $\mathcal{P}_{2}(x)$ is the second Legendre polynomial. For simiplicity,
we drop the subscript $i$, and understand that $\mathbf{\hat{e}}$ refers to
the distinguished direction of the cell. The order parameter tensor is the
normalized traceless second moment of the orientational distributions
function; the first moment vanishes due to the quadrupolar symmetry of the
rodlike particles. We focus primarily on $\mathbf{Q}$ in this work, however,
higher order moments are also present, and shall be discussed below.

In general, self-consistency requires that free parameters of the mean-field
pseudopotential be chosen so as to minimize the free energy. In the present
setting, it is the cell shape itself that needs to be so chosen. We next
address this issue in general for a given cell volume $V_{cell}$ before
considering more specific instances.

###  2.1 Cell shape and pressure

We first show that, independent of the shape of the cell, the Cauchy stress
tensor $\bm{\sigma}$, associated at equilibrium with the free energy $F$, is
isotropic, and so its effect on the walls of the cell reduces to uniform
pressure. Here we consider the minimiser $\mathcal{C}$ over a family of cells
described as a linear, volume preserving transformation of a reference cell
$\mathcal{C}_{0}$. We assume that the cell $\mathcal{C}_{0}$ is of some
arbitrary shape (say, a parallelepiped, for definiteness), and let
$\mathbf{F}$ denote a corresponding deformation gradient. Thus $F$ can be
regarded as a function of $\mathbf{F}$, whose explicit form need not be
specified at this stage, provided that it is sufficiently smooth. The Cauchy
stress $\bm{\sigma}$ exerted on the walls of $\mathcal{C}$ is [11]

$\bm{\sigma}=\frac{1}{\det{\mathbf{F}}}\frac{\partial
F}{\partial{\mathbf{F}}}\cdot{\mathbf{F}}^{T}.$ (12)

For $\mathcal{C}$ and $\mathcal{C}_{0}$ to have the same volume $V_{cell}$,
$\mathbf{F}$ must obey the constraint

$\det\mathbf{F}\equiv 1.$ (13)

Stationarity of $F$ over this constraint requires that

$\frac{\partial
F}{\partial\mathbf{F}}=\lambda\frac{\partial}{\partial\mathbf{F}}\det\mathbf{F},$
(14)

where $\lambda$ is a Lagrange multiplier. Since, for an invertible tensor
$\mathbf{F}$,

$\frac{\partial}{\partial\mathbf{F}}\det{\mathbf{F}}=(\det\mathbf{F})\mathbf{F}^{-T},$
(15)

it follows from Eqs. (12), (14), and (15) that

$\bm{\sigma}=\frac{1}{\det{\mathbf{F}}}\frac{\partial
F}{\partial{\mathbf{F}}}\cdot{\mathbf{F}}^{T}=\lambda\mathbb{I}.$ (16)

Eq. (16) indicates that the pressure is isotropic in equilibrium, and the
Lagrange multiplier reduces to the (isotropic) pressure.

We remark that a stress also acts on the particle inside the cell, and if the
pressure is isotropic on the cell walls, it will not, in general, be isotropic
on the particle - and vice versa.

##  3 Results

In this section, we first present results for a uniaxial ellipsoid in a
rectangular cell where the free volume is calculated by numerical integration,
and then compare these results with those from molecular dynamics. We then
present the case of a needle in an ellipsoidal cell, where an analytical
solution is possible.

###  3.1  Uniaxial ellipsoid in a rectangular cell

Here we consider the case when the cell is a parallelepiped and the particle
is a uniaxial ellipsoid. We conjecture, as suggested by numerics, that the
orientational averaged free volume of an ellipsoid in a parallelepiped cannot
attain its maximum and free energy its minimum when the cell faces are not
orthogonal. Furthermore, if the particle has uniaxial symmetry, the maximum
free volume occurs when the cell is uniaxial. Therefore here we consider
uniaxial rectangular cells, and apply the mean field theory of Section 2 to
the case of a rigid uniaxial ellipsoid in a uniaxial rectangular cell.

Let $a$ and $b$ be the lengths of the semiaxes of the uniaxial ellipsoid in
the directions along and perpendicular to the symmetry axis. The equation of
the ellipsoid is

$\mathbf{rAr}=1,$ (17)

where

$\mathbf{A}=\frac{1}{b^{2}}\mathbb{I}+(\frac{1}{a^{2}}-\frac{1}{b^{2}})\mathbf{\hat{m}\hat{m},}$
(18)

and $\mathbf{\hat{m}}$ is along the symmetry axis. The aspect ratio, the ratio
of the length of the distinguished semi-axis to that of one of the others, is
$\ \eta=a/b$. Alternately, we can consider the uniaxial affine deformation of
a sphere of radius $R$ to form this ellipsoid; the principal stretches are
$\lambda_{1}=a/R$, $\lambda_{2}=\lambda_{3}=b/R$. If the volume is conserved,
$\lambda_{2}=\lambda_{3}=1/\sqrt{\lambda_{1}}$, and the aspect ratio

$\eta=\lambda^{3/2},$ (19)

where $\lambda$ indicates stretch along the distinguished principal direction.
The aspect ratio $\eta$ and stretch $\lambda$ will be used interchangeably to
describe particle and cell shapes.

The particle volume is $v_{p}=4\pi ab^{2}/3$ and the occupied volume fraction
is $\phi=v_{p}/V_{cell}=v_{p}\rho_{0}$.

The extent of the ellipsoid with orientation $\mathbf{\hat{m}}$ in direction
$\mathbf{\hat{N}}$ is

$h(\mathbf{\hat{N})}=\sqrt{a^{2}(\mathbf{\hat{m}\cdot\hat{N}})^{2}+b^{2}(1-(\mathbf{\hat{m}\cdot\hat{N}})^{2})},$
(20)

and in the principal directions of the cell this is

$h_{i}=\sqrt{a^{2}m_{i}^{2}+b^{2}(1-m_{i}^{2})},i=1,2,3.$ (21)

The free volume available to the center of the ellipsoid is then simply

$V_{f}(\mathbf{\hat{m}},\mathcal{C})=\prod_{i=1}^{3}\max\\{d_{i}-2h_{i},0\\},$
(22)

where $d_{i}$’s are lengths of the cell along its principal directions. A 2D
illustration of the free volume of an ellipse in a rectangular cell is given
in Fig. 1.

Figure 1: Illustration of the free volume available to a convex particle in a
polygonal cavity.

The free energy of one cell, from Eq. (7), is

$F_{c}=-kT\ln\int_{S^{2}}\prod_{i=1}^{3}\max(d_{i}-2h_{i}(\mathbf{\hat{m}}),0)d\mathbf{\hat{m}},$
(23)

where we have omitted $-3\ln\Lambda_{dB}$. For simplicity, we will drop the
subscript, and with the understanding that the free energy $F$ is for a single
cell. The free energy $F$ can therefore be calculated by straightforward
numerical integration of Eq. (23) using, say, the quasi-Monte Carlo method.
Local minima/maxima can be identified via standard gradient descent/ascent
methods.

Let $d_{2}=d_{3}$ be the lengths of the cell along two principal directions,
and $d_{1}$ be the length along the third, the distinguished principal
direction. The aspect ratio of the cell is $\eta_{cell}=d_{1}/d_{2}$; the
stretch is $\lambda_{cell}=\eta_{cell}^{2/3}$. The cell is considered prolate
if $\eta_{cell}>1$ and oblate if $\eta_{cell}<1$. In calculating the
equilibrium free energy, at each volume fraction, we use $\eta_{cell}$, which
minimizes the free energy, to satisfy self-consistency.

We note that if the maximum (linear) dimension of the ellipsoid exceeds the
minimum dimension of the cell, then the cell will exclude some orientations of
the ellipsoid. In the case of a cubic domain, with $\eta_{cell}=1$, this
occurs at a critical volume fraction; the corresponding critical volume
fractions are

$\phi_{cp}=\frac{\frac{4\pi}{3}ab^{2}}{(2a)^{3}}=\frac{\pi}{6}\left(\frac{b}{a}\right)^{2}$
(24)

for prolate and

$\phi_{co}=\frac{\frac{4\pi}{3}ab^{2}}{(2b)^{3}}=\frac{\pi}{6}\left(\frac{a}{b}\right)$
(25)

for oblate ellipsoids.

(a)

(b)

(c)

(d)

Figure 2: Results for uniaxial ellipsoids with aspect ratio
$\eta_{particle}=2$ (left column) and $\eta_{particle}=0.5$ (right column) in
a rectangular box. Top row: pressure $P$ (left axes) and cell aspect ratio
$\eta_{cell}$ (right axes) vs. volume fraction $\phi$. Bottom row: free energy
$F$ (left axes) and scalar order parameter $S_{2}$ (right axes) vs. volume
fraction $\phi$. Only locally stable solutions are shown. Dotted curves are
for isotropic phases where the cell is a cube, the dashed curves are for
oblate phases with $S_{2}<0$ and solid curves are for prolate phase with
$S_{2}>0$.

(a) $(\phi,\eta_{cell})=(0.13,1)$

(b) $(\phi,\eta_{cell})=(0.153,1.17)$

(c) $(\phi,\eta_{cell})=(0.23,1.57)$

(d) $(\phi,\eta_{cell})=(0.15,0.78)$

(e) $(\phi,\eta_{cell})=(0.29,0.65)$

Figure 3: Orientational distribution function for a uniaxial prolate ellipsoid
with aspect ratio $\eta_{particle}=2$ at various volume fractions $\phi$ in a
cube $\eta_{cell}=1$(a), a prolate cell with $\eta_{cell}>1$ (b,c) and an
oblate cell with $\eta_{cell}<1$ (d,e). Each plot is scaled so that its
maximum value is 1, and the regions in black correspond to where the
probability is zero.

Figure 2 shows the results for two representative ellipsoids with aspect
ratios $2$ and $1/2$, respectively, as function of occupied volume fraction.
Volume fraction was varied by changing particle volumes at fixed particle
aspect ratio. Anisotropy at large volume fractions is apparent both in the
cell shape, and in the orientational distribution function. Four scenarios
have been observed: prolate or oblate particles in prolate or oblate cells. If
the particle and cell shape are similar, then $S_{2}>0$. If the particle and
cell shapes are dissimilar, then $S_{2}<0$.

For the prolate ellipsoid with $\eta_{particle}=2$, the system is dilute with
$\phi<\phi_{cp}$, the whole orientation state space is accessible, the cubic
cell with $\eta_{cell}=1$ is the only minimiser of $F$, and the orientational
order parameter $S_{2}=0$, the system is isotropic. At $\phi=\phi_{cp}$, the
isotropic state loses its stability and two orientationally ordered solutions
bifurcate from the isotropic branch. The prolate branch with
$S_{2}>1,\eta_{cell}>1$ has lower energy and pressure than those of the oblate
branch with $S_{2}<0,\eta_{cell}>1$ for $\phi>\phi_{cp}$.

For the oblate ellipsoid with $\eta_{particle}=0.5$, although the bifurcation
is at $\phi_{c}$, orientational order can be observed both below and above
$\phi_{c}$. In fact, there are two first order transitions from isotropic to
both prolate nematic phase and oblate nematic phase at two different volume
fractions slightly below $\phi_{c}$. The prolate branch with
$S_{2}>0,\eta_{cell}<1$ possesses lower energy and pressure than those of the
oblate branch with $S_{2}>0,\eta_{cell}>1$.

For sufficiently large volume fractions, certain orientations are disallowed
by the shape of the cell, and the orientational distribution function has
nontrivial compact support. This is indicated in Fig. 3(b-e), where
representative equilibrium orientational probability density functions for a
prolate ellipsoid with $\eta_{particle}=2$ in equilibrium prolate
$\eta_{cell}>1$ (Figs. 3(b) and (c)) and oblate $\eta_{cell}<1$ (Figs. 3(d)
and (e)) cells are shown. One can readily see that the distribution function
lacks axial symmetry, even in the very dilute limit where $S_{2}=0$, as shown
in Fig. 3(a). This is because an uniaxial ellipsoid has more free volume when
it aligns along the diagonals of a rectangular cell compared to other
orientations.

As the volume fraction is increased, a topological transition takes place: the
connected accessible orientational space becomes disconnected forming a number
of disconnected regions, where the orientation of the ellipsoid is locked in a
small region, as can be observed in going from Figs. 3(b) to (c) and from (d)
to (e). These topological changes also manifest themselves as kinks on the
bifurcation diagrams as in Fig. 2 and Fig. 4, although the topological change
on the prolate branch is not as apparent as the one on the oblate branch.

Since the orientational distribution function is known, in addition to
$S_{2}$, higher moments
$S_{n}=\langle\mathcal{P}_{n}(\mathbf{\hat{m}\cdot\hat{e}})\rangle$ can be
readily calculated. Here $\mathcal{P}_{n}$ is the Legendre polynomial of
degree $n$. Although in the isotropic phase $S_{2}=0$, interestingly, higher
moments may exist as nonzero functions of the occupied volume fraction $\phi$,
as shown in Fig. 4(b). Given the edges and corners of the box, it is not
surprising that the distribution function has higher moments which are nonzero
at all volume fractions.

The free volume of an arbitrarily shaped particle in a rectangular cell can be
calculated using a support function based approach, where the support function
plays the role of the function $h$ in Eq. (22). We also remark that as long as
the cell is convex, the free volume for a concave particle will be identical
to that of its convex hull. Applying the support function approach to extend
our model to particles with more complex shapes has already been initiated by
one of us (JMT).

We realize that cell shapes other than rectangular boxes may allow more free
volume to ellipsoidal particles; we expect that to every particle shape there
corresponds an optimal cavity shape with the lowest free energy in
equilibrium. The rectangular box is just our first attempt; an ellipsoidal
cell is considered in Section 3.3.

###  3.2 Molecular dynamics

Molecular dynamics (MD) is a useful tool for studying properties of systems of
interacting particles. MD simulations keep track of the motion of individual
particles relying on classical mechanics. The results of MD simulations can be
used to test the validity of statistical mechanical models. Here, we perform
MD simulations of a hard rigid ellipsoidal particle moving in a rectangular
cell, undergoing collisions with the cell walls. We imagine that this
approximates the real physical situation where thermally excited molecules
collide with nearest neighbors. As we show in this subsection, MD gives the
same predictions for the equilibrium shape of the cell, for the average
orientation of the ellipsoid, and for the pressure as our statistical
mechanical model.

The dynamics consists of free flights and collisions with cell walls. When the
particle does not undergo a collision, the position of its center of mass is
updated as

$\mathbf{r(}t_{i+1})=\mathbf{r}(t_{i})+\mathbf{v}_{cm}\Delta t_{i}.$ (26)

However, the angular velocity of the particle changes in time even without
external torque. The angular acceleration in an inertial reference frame is
given by Euler’s equations for rigid body dynamics,

$\dot{\bm{\omega}}(t_{i})=\mathbf{-I}^{-1}\cdot(\bm{\omega}\times(\mathbf{I}\cdot\bm{\omega})),$
(27)

where $\mathbf{I}$ is the tensor of inertia of the ellipsoid, given by

$\mathbf{I}=(\frac{2}{5}mb^{2}-\frac{1}{5}m(a^{2}+b^{2}))\mathbf{\hat{m}\hat{m}}+\frac{1}{5}m(a^{2}+b^{2})\mathbb{I}\mathbf{,}$
(28)

where $m$ is its uniformly distributed mass. We then update the angular
velocity and the orientation of the particle as follows,

$\displaystyle\bm{\omega}(t_{i})$ $\displaystyle=$
$\displaystyle\bm{\omega}(t_{i-1})+\bm{\dot{\omega}}\Delta t_{i},$ (29)
$\displaystyle\mathbf{\hat{m}}(t_{i+1})$ $\displaystyle=$
$\displaystyle\frac{\mathbf{\hat{m}}(t_{i})+\bm{\omega}(t_{i})\times\mathbf{\hat{m}}(t_{i})\Delta
t_{i}}{|\mathbf{\hat{m}}(t_{i})+\bm{\
\omega}(t_{i})\times\mathbf{\hat{m}}(t_{i})\Delta t_{i}|},$ (30)

where we use the default value of $\Delta t_{i}=10^{-3}$. If the updated
configuration causes any part of the particle to leave the cell, we return to
the previous step, and halve $\Delta t_{i}$. The process is repeated until the
new configuration of the particle is completely within the cell. When the time
step size reaches a lower threshold (e.g., $\Delta t_{i}=10^{-14}$), collision
with the wall occurs.

We assume that both the particle and cell walls are frictionless, thus the
impulses are along the normals to the surfaces in contact. In a collision, the
conservation of linear momentum, angular momentum and kinetic energy, gives
the postcollision center of mass velocity $\mathbf{v}_{f}$ and angular
velocity $\bm{\omega}_{f}$ by ([12])

$\mathbf{v}_{f}=\mathbf{v}_{i}-\frac{J}{m}\mathbf{\hat{N}},$ (31)
$\bm{\omega}_{f}=\bm{\omega}_{i}-J\mathbf{I}^{-1}\cdot(\mathbf{p\times\hat{N}}),$
(32)

where $\mathbf{v}_{i}$ and $\bm{\omega}_{i}$ are precollision velocity and
angular velocity, $\mathbf{\hat{N}}$ is the outward normal of the ellipsoid
pointing towards the cell wall. The vector $\mathbf{p}$ is a body-fixed vector
joining the center of mass of the ellipsoid to the point of contact with the
cell wall at the instant of collision,

$\mathbf{p}=\frac{\mathbf{A}^{-1}\cdot\mathbf{\hat{N}}}{\sqrt{\mathbf{\hat{N}}\cdot\mathbf{A}^{-1}\cdot\mathbf{\hat{N}}}}$
(33)

and the positive scalar $J$ is the magnitude of the impulse transferred to the
ellipsoid at the collision,

$J=2\frac{(\mathbf{\mathbf{v}}_{i}+\mathbf{\mathbf{\omega}}_{i}\times\mathbf{p)\cdot\hat{N}}}{1/m-\mathbf{\hat{N}}\cdot\mathbf{p}\times\mathbf{I}^{-1}\times\mathbf{p}\cdot\mathbf{\hat{N}}}.$
(34)

The pressure on each wall is the average impulse per area on the wall. If the
pressure happens not to be isotropic, we then update the relevant length of
the cell, denoted $d_{l}$, according to

$d_{l}\rightarrow d_{l}+\gamma\frac{P_{l}-P_{t}}{P_{l}},$ (35)

where $P_{l}$ and $P_{t}$ represent the averaged pressures on two square walls
and four rectangular walls, respectively, and $\gamma=0.01$ plays the role of
compressibility. The updates for the other two dimensions are determined by
maintaining constant volume and a uniaxial shape.

A very great deal of work has been done on simulations of systems of hard
ellipsoids[13],[14],[15],[16]. Our system is different, however, due to the
mean field aspect of our model, where only a single particle is present,
enclosed by a cell representing the effects of the other particles.
Confinement induced ordering of hard ellipsoids has also been studied between
two hard walls [17] however our confinement is in 3D; furthermore the shape of
the confining shape is adaptive.

The scalar order parameter is calculated according to
$S_{n}=\langle\mathcal{P}_{n}(\mathbf{\hat{m}\cdot\hat{e}})\rangle$ where the
average is over all instantaneous orientations of the symmetry axis
$\mathbf{\hat{m}}$ of the ellipsoid between collisions. All simulation results
were averaged over $10^{3}$ random initial conditions.

Figure 4 shows the agreement of the results of free volume calculation in
Sect. 3.1 and MD simulations for an ellipsoid with aspect ratio $4$. The
transition from the isotropic phase to the prolate phase is continuous,
whereas the transition from the isotropic phase to the oblate phase is
discontinuous. This is different from the case of $\eta_{particle}=2$ where
both transitions are continuous.

(a)

(b)

Figure 4:  Results for ellipsoids with aspect ratio $\eta_{particle}=4$. (a)
Pressure $P$ (left axes) and cell aspect ratio $\eta_{cell}$ (right axis) vs.
volume fraction $\phi$. (b) Scalar order parameters $S_{2}$ (left axes) and
$S_{4}$ (right axes) vs. volume fraction $\phi$. Only locally stable phases
are shown. Solid lines are prolate phase, and dashed lines are for oblate
phase, and dotted lines for isotropic phase. The MD results are shown as
discrete points with error bars.

###  3.3 Needle in an ellipsoidal cell

In this subsection, we consider a needle within an ellipsoidal cell. This
geometry allows for an analytic free volume expression, and reveals a topology
common to the case of very slender particles.

(a)

(b)

Figure 5: The undeformed spherical globule $\mathcal{C}_{0}$ enclosing a rigid
rod (a) and its ellipsoidal image $\mathcal{C}$ under affine deformation (b).

Consider $\mathcal{C}_{0}$ to be a spherical cavity $\mathscr{B}_{R}$ of
radius $R$ enclosing a rigid rod of length $2a$ (see Fig. 5(a)) and subject
the sphere to a homogeneous, isochoric deformation $\bm{f}$ with gradient
$\mathbf{F}$, which, modulo a rotation, changes the sphere into the ellipsoid
$\mathcal{C}$ depicted in Fig. 5(b).

Now, the unit vector $\mathbf{\hat{m}}$ denotes the orientation of the rod,
such that if the center of the rod is at the point $\bm{y}$, its ends occupy
the points $\bm{y}_{\pm}=\bm{y}\pm a\mathbf{\hat{m}}$, respectively. We want
to determine the volume of the free region
$\mathscr{R}_{\mathrm{f}}(\mathbf{\hat{m}})$ accessible to $\bm{y}$ within the
ellipsoid $\mathcal{C}$, for a given orientation $\mathbf{\hat{m}}$ of the
rod. Writing $\bm{y}_{\pm}=\mathbf{F}\bm{x}_{\pm}$ and
$\bm{y}=\mathbf{F}\bm{x}$, so that $\bm{x}_{\pm}$ are the preimages of
$\bm{y}_{\pm}$, and $\bm{x}$ the preimage of $\bm{y}$, the region
$\bm{f}^{-1}(\mathscr{R}_{\mathrm{f}})$ comprises all points $\bm{x}$ such
that

$\mathbf{F}\bm{x}\pm a\mathbf{\hat{m}}=\mathbf{F}\bm{x}_{\pm}\quad\text{for
both}\quad\bm{x}_{\pm}\in\mathscr{B}_{R},$ (36)

from which, since $\mathbf{F}$ is invertible, we have that

$\bm{x}\pm
a\mathbf{F}^{-1}\mathbf{\hat{m}}=\bm{x}_{\pm},\quad\bm{x}_{\pm}\in\mathscr{B}_{R},$
(37)

and so

$\bm{f}^{-1}(\mathscr{R}_{\mathrm{f}})=\\{\bm{x}:|\bm{x}+a\mathbf{F}^{-1}\mathbf{\hat{m}}|\leq
R\\}\cap\\{\bm{x}:|\bm{x}-a\mathbf{F}^{-1}\mathbf{\hat{m}}|\leq R\\}.$ (38)

Geometrically, this set is the region comprised of two spheres of radius $R$
whose centers are $2b$ apart, where

$b:=a|\mathbf{F}^{-1}\mathbf{\hat{m}}|.$ (39)

Since $\det\mathbf{F}=1$, the free volume $V_{\mathrm{f}}$ is the same as the
volume of $\bm{f}^{-1}(\mathscr{R}_{\mathrm{f}})$, that is, twice the volume
of a spherical frustum of radius $r:=\sqrt{R^{2}-b^{2}}$ and height $h:=R-b$,

$V_{\mathrm{f}}=\frac{\pi
h}{3}(3r^{2}+h^{2})=V_{cell}\left(1-\frac{3}{2}\alpha\beta+\frac{1}{2}\alpha^{3}\beta^{3}\right),$
(40)

where $V_{cell}$ is the volume of $\mathscr{B}_{R}$ and we have set

$\alpha=\frac{a}{R}\leq 1\quad\text{and}\quad\frac{b}{R}=\alpha\beta\leq 1.$
(41)

It readily follows from (38) and (40) that

$\beta=\sqrt{\mathbf{\hat{m}}\cdot\mathbf{B}^{-1}\mathbf{\hat{m}}},$ (42)

where $\mathbf{B}=\mathbf{F}\mathbf{F}^{\mathsf{T}}$ is the left Cauchy-Green
tensor, which in the frame
$(\mathbf{\hat{e}}_{1},\mathbf{\hat{e}}_{2},\mathbf{\hat{e}}_{3})$ of the
principal directions of stretching is represented as

$\mathbf{B}=\lambda_{1}^{2}\mathbf{\hat{e}}_{1}\mathbf{\hat{e}}_{1}+\lambda_{2}^{2}\mathbf{\hat{e}}_{2}\mathbf{\hat{e}}_{2}+\lambda_{3}^{2}\mathbf{\hat{e}}_{3}\mathbf{\hat{e}}_{3},$
(43)

where $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$ are the principal
stretches subject to $\lambda_{1}\lambda_{2}\lambda_{3}=1$. Letting $m_{i}$ be
the components of $\mathbf{\hat{m}}$ in the frame
$(\mathbf{\hat{e}}_{1},\mathbf{\hat{e}}_{2},\mathbf{\hat{e}}_{3})$, they can
be written as

$m_{1}=\sin\theta\cos\varphi,\quad m_{2}=\sin\theta\sin\varphi,\quad
m_{3}=\cos\theta,$ (44)

for $\theta\in[0,\pi]$ and $\varphi\in[0,2\pi]$. By use of (44) and (43) in
(42) , we arrive at

$\beta=\sqrt{\frac{1}{\lambda_{1}^{2}}\sin^{2}\theta\cos^{2}\varphi+\frac{1}{\lambda_{2}^{2}}\sin^{2}\theta\sin^{2}\varphi+\frac{1}{\lambda_{3}^{2}}\cos^{2}\theta}.$
(45)

For simplicity, we shall assume that $\mathbf{B}$ has the uniaxial symmetry
around $\hat{\mathbf{e}}_{3}$, so that
$\lambda_{1}=\lambda_{2}=1/\sqrt{\lambda_{3}}$ and $\beta$ reduces to

$\beta=\sqrt{{\lambda}\sin^{2}\theta+\frac{1}{\lambda^{2}}\cos^{2}\theta},\text{
with }\lambda=\lambda_{3}.$ (46)

With this choice, the ellipsoid is prolate along $\hat{\mathbf{e}}_{3}$ if
$\lambda>1$ and oblate if $\lambda<1$.

For given $\lambda$, the domain $I_{\lambda}$ of admissible values of $\theta$
is restricted by the second inequality in (41). A direct inspection shows that
where

$I_{\lambda}=\left\\{\begin{array}[]{cc}[\chi,\pi-\chi],&0\leq\lambda\leq\alpha\\\
[0,\pi],&\alpha\leq\lambda\leq\frac{1}{\alpha^{2}}\\\
[0,\chi]\cup[\pi-\chi,\pi],&\lambda\geq\frac{1}{\alpha^{2}}\end{array}\right.$
(47)

where

$\chi=\arccos\left(\frac{\lambda}{\alpha}\sqrt{\frac{1-\lambda\alpha^{2}}{1-\lambda^{3}}}\right).$
(48)

The free energy of the rod (in units $kT$) is then given by

$F(\alpha,\lambda)=-\ln\left(\pi\int_{I_{\lambda}}(2-3\alpha\beta+\alpha^{3}\beta^{3})\sin\theta
d\theta\right)+\ln[2\pi(\alpha-1)^{2}(\alpha+2)],$ (49)

where the added constant has been gauged so as to ensure that $F(\alpha,1)=0$.
The integral in (49) can be expressed in terms of elementary functions, not
given explicitly here in the interest of space. Plots of $F$ as function of
$\lambda$ are shown in Fig. 6 for several values of $\alpha$.

(a) $\alpha=\alpha_{c}^{(1)}\dot{=}0.737$

(b) $\alpha=\alpha_{t}\dot{=}0.757$

(c) $\alpha=\alpha_{c}^{(2)}\dot{=}0.805$

(d) $\alpha=\sqrt{2/3}\dot{=}0.816$

(e) $\alpha=\alpha_{c}^{(3)}\dot{=}0.826$

(f) $\alpha=0.850$

Figure 6: The dimensionless free energy $f$ as in (49) plotted against
$\lambda$ for different values of $\alpha$.

For $\alpha<\alpha_{\mathrm{c}}^{(1)}\doteq 0.737$, $F$ has a single critical
point at $\lambda=1$, its absolute minimum representing the isotropic state
where the spherical globule remains undeformed. For
$\alpha>\alpha_{\mathrm{c}}^{(1)}$, a local minimum and a local maximum
develop in the oblate branch (where $\lambda<1$), while the isotropic state
remains the absolute minimum, until $\alpha$ reaches the transition value
$\alpha_{\mathrm{t}}\doteq 0.757$, where both the isotropic state and the
oblate minimum have the same energy. A first-order transition takes place at
$\alpha=\alpha_{\mathrm{t}}$, with the oblate minimum at $\lambda\doteq
0.518$. At $\alpha=\alpha_{c}^{(2)}\dot{=}0.805$, an inflection point develops
on the prolate branch, from which a local maximum and a local minimum emanate,
the former approaching the equilibrium isotropic state at $\lambda=1$. At
$\alpha=\sqrt{\frac{2}{3}}\doteq 0.816$, the isotropic state becomes unstable
and two branches emanate from it, a local prolate minimum and a local oblate
maximum: the former connects with the prolate maximum at
$\alpha=\alpha_{c}^{(3)}\dot{=}0.826$, while the latter connects with the
oblate minimum at $\alpha=\alpha_{c}^{(1)}$. The locally stable prolate branch
converges to $\lambda\doteq 2.627$ as $\alpha\rightarrow 1$, whereas the
globally stable oblate branch converges to $\lambda\dot{=}0.254$ as $\alpha\to
1$.

Figure 7 summarizes the above results, which are also compared with the
bifurcation diagrams for a needle in a rectangular cell. We note that the
topologies are mildly distinct. Whilst the oblate branch has a similar
topology in each case, two unstable branches, one prolate and one oblate,
bifurcate from the isotropic state at the critical value $\alpha=1$ in the
cuboidal case, with the former being a very small region that is difficult to
see in the figure.

(a)

(b)

Figure 7: Bifurcation diagrams. (a) Needle in an ellipsoid cell. (b) Needle in
a rectangular cell. In the former case $\alpha=2a/\sqrt[3]{3V_{cell}/(4\pi)}$,
while in the latter $\alpha=2a/\sqrt[3]{V_{cell}}$. Solid red lines are local
minima, dotted blue lines are local maxima. Red circles mark the transition
from the isotropic state to the absolute oblate minimiser. Blue circles
represent the critical points discussed in the text.

##  4 Work by shape change

Our mean field model gives some insights about the effects of shape change by
individual molecules on their neighbors. A number of photomechanical materials
consist of a nematic liquid crystal elastomer, containing a small fraction
($\sim 2\%$) of photoisomerizable azo dyes. These, in their ground state
trans- configuration are elongated and in essence very similar to the liquid
crystal molecules of the host. When illuminated, however, they undergo a
trans-cis photoisomerization, and their shape changes from the elongated
trans- form to the more compact – cis- form. Based on our model, we can
estimate the free energy increase and stress on the host due to this shape
change.

The aspect ratio of the trans- isomer is (roughly) known; we can also estimate
the volume fraction. We know that in equilibrium the free energy is a minimum;
so we minimize the free energy with respect to the cell aspect ratio. A
possible starting configuration is shown in Fig. 8(a); the box shape is
optimal for the given volume fraction $\phi=0.322$ and ellipsoid aspect ratio
of $2$; the sample is in equilibrium and the pressure on the box is isotropic.
The sample is then illuminated with UV light, and there is photoisomerization
where the dye molecule changes its shape from the elongated trans- to the more
spherical cis- shape. As the ellipsoid suddenly becomes more spherical, the
system loses equilibrium, the free energy dramatically increases and the
pressure becomes anisotropic as shown in Fig. 8(b). This is the initial free
energy for the photowork, shown in Fig. 8(b), together with pressure, which is
now strongly anisotropic.

If allowed, the cell shape then changes to aspect ratio close to $1$ to
minimize the free energy for the new ellipsoid shape. The photomechanical work
is carried out by the stress on the cell walls. The free energy corresponding
to box aspect ratio $1$ is the final energy of the system after the photowork.
The free energy difference between the initial and final states is the energy
available for photowork.

In an isotropic medium, dissolved azo dyes behave like standard two-level
systems, governed by a nearly constant potential energy difference between the
tran- and cis- states. In a liquid crystal host, however, potential energy of
each isomer111internal energy plus the work required to insert the isomer into
the standard cell in the system would sensitively depend on the anisotropy of
the cell. One would expect therefore that the temperature dependence of isomer
populations in liquid crystals to significantly differ from that of a simple
two-level system.

Energy for the photoexcitation comes from the absorbed photon; what is not
used in increasing the free energy is released as heat. Matching the free
energy increase to the absorbed photon energy may be a viable strategy to
increase the efficiency of photomechanical materials.

(a)

(b)

Figure 8: The particle shape change (decrease in aspect ratio
$\eta_{particle}$) at $\phi=0.322$ causes a large ($\sim$ factor of 2)
increase in the free energy $F$, and a large increase ($\sim$ factor of 10) in
the pressure perpendicular to the symmetry axis.

##  5 Conclusion

In this work, we have considered a mean field model of $N$ indistinguishable
particles with hard core steric interactions in a volume $V$ at temperature
$T$. The region is divided into $N$ identical cells, each with one particle in
each cell. Each cell has a volume $V/N$. At nonzero temperatures, the
particles undergo collisions with the cell walls which approximate the effects
of neighboring particles. In equilibrium, satisfying self-consistency, the
cell adopts a shape which minimizes the free energy, and leads to an isotropic
pressure tensor.

We have applied this mean field cell model to the case of hard uniaxial
ellipsoids in rectangular and ellipsoidal cells.

At low occupied volume fractions, the shortest dimension of the cell is longer
than the longest dimension of the particles. Hence all orientations of the
particles are possible in the cells; here the particles are orientationally
disordered, and the orientational order parameter - the second moment of the
orientational distribution - is zero. As the occupied volume fraction
increases, there is a critical volume fraction where the shortest dimension of
the cell is equal to the longest dimension of the particles. Above this
critical volume fraction, the orientational distribution function has compact
support; there is a bifurcation and the system becomes orientationally
ordered. One of the orientationally ordered solutions is prolate, while the
other is oblate; the stability is determined by the shape of the cell. The
results of mean-field theory are fully compatible with the the results of
molecular dynamics simulations.

The above results give insights towards understanding the shape change
exhibited by photomechanical materials. In such azo-dye doped liquid crystal
elastomer materials, illumination causes photoisomerization of the dye. Due to
the shape change, the system loses equilibrium, the free energy increases and
anistoropic stress appears in the cell, which changes the shape of the cell
and of the bulk host, and can do mechanical work. Energy for the work comes
from the absorbed light. The model not only describes the mechanism whereby
particle shape change results in local increased free energy and anisotropic
internal stress, but allows quantitative estimates of these.

This paper has reported our first attempts to implement a simple mean field
theory for rod-like hard particles. Work to extend and improve the model is
clearly needed. Allowing the confined particles to respond to the anisotropic
internal pressures acting on them seems promising; as does consideration of
cell and particle shapes more general than those considered here. Efforts in
these directions are under way.

#### Acknowledgements

P.P-M. acknowledges support from the Office of Naval Research through the MURI
on Photomechanical Material Systems (ONR N00014-18-1-2624).

## References

  * [1] Camacho-Lopez M., Finkelmann, H.,Palffy-Muhoray, P. and Shelley, M. Fast liquid crystal elastomer swims into the dark, Nature Materials, 3, 307-310 (2004)
  * [2] White T.J., Broer, D.J. Programmable and adaptive mechanics with liquid crystal polymer networks and elastomers, Nature Materials, 14, 1087-1098 (2015)
  * [3] M. G. Kuzyk and N.J. Dawson, Photomechanical materials and applications: a tutorial, Advances in Optics and Photonics 12, 847 (2020)
  * [4] Guo, T., Svanidze,A., Zheng, X. Palffy-Muhoray, P. Regimes in the Response of Photomechanical Materials, Appl. Sci. 12, 7723, 2022
  * [5] Onsager L. The effects of shape on the interaction of colloidal particles. Annals of the New York Academy of Sciences. 1949 May;51(4):627-59.
  * [6] Lennard-Jones JE, Devonshire AF. Critical phenomena in gases-I. Proceedings of the Royal Society of London. Series A-Mathematical and Physical Sciences. 1937 Nov 5;163(912):53-70.
  * [7] Kincaid JF, Eyring H, Stearn AE. The Theory of Absolute Reaction Rates and its Application to Viscosity and Diffusion in the Liquid State. Chemical Reviews. 1941 Apr;28(2):301-65.
  * [8] Kirkwood JG. Critique of the free volume theory of the liquid state. The Journal of Chemical Physics. 1950 Mar;18(3):380-2.
  * [9] A.L. Kuzemsky, Variational Principle of Bogoliubov and Generalized Mean Fields in Many-Particle Interacting Systems, Int.J. Mod. Phys. B 29, (2015)
  * [10] Palffy-Muhoray, P., The single particle potential in mean field theory, Am. J. Phys. 70, 422-437 (2002)
  * [11] Spencer, A.J.M., 2004. Continuum mechanics. Courier Corporation.
  * [12] Palffy-Muhoray, P., Virga, E. G., Wilkinson, M., Zheng, X., On a paradox in the impact dynamics of smooth rigid bodies. Mathematics and Mechanics of Solids, 24(3), 573-597 (2019).
  * [13] D. Frenkel and B. M. Mulder, Mol. Phys. 55, 1171 (1985)
  * [14] J.W.Perram, M.S.Wertheim, J.L. Lebowitz ,and G.O.Williams,Chem. Phys.Lett.105,277(1984).
  * [15] G. Bautista-Carbajal, A. Moncho-Jorda and G. Odriozola, Further details on the phase diagram of hard ellipsoids o frevolution, J. Chem. Phys. 138, 064501 (2013)
  * [16] B.J.Berne and P.Pechukas,J.Chem.Phys.56,4213(1972)
  * [17] H. Miao and H. Ma, Confinement Induced Ordering in Fluid of Hard Ellipsoids, Chin. J. Chem. Phys. 29, 212 (2016)