# The Thermodynamic Origins of Chiral Twist in Monolayer Assemblies of Hard
Rod-like Colloids

Yawei Liu<EMAIL_ADDRESS>ARC Centre of Excellence in Exciton
Science, School of Chemistry, University of Sydney, Sydney, New South Wales
2006, Australia  Beijing Key Laboratory of Ionic Liquids Clean Process, CAS
Key Laboratory of Green Process and Engineering, State Key Laboratory of
Multiphase Complex Systems, Institute of Process Engineering, Chinese Academy
of Sciences, Beijing, 100190, China  Jared A. Wood  ARC Centre of Excellence
in Exciton Science, School of Chemistry, University of Sydney, Sydney, New
South Wales 2006, Australia  The University of Sydney Nano Institute,
University of Sydney, Sydney, New South Wales 2006, Australia  Achille
Giacometti  Dipartimento di Scienze Molecolari e Nanosistemi, Università Ca’
Foscari di Venezia Campus Scientifico, Edificio Alfa, via Torino 155,30170
Venezia Mestre, Italy  European Centre for Living Technology (ECLT) Ca’
Bottacin, 3911 Dorsoduro Calle Crosera, 30123 Venice, Italy  Asaph Widmer-
Cooper<EMAIL_ADDRESS>ARC Centre of Excellence in Exciton
Science, School of Chemistry, University of Sydney, Sydney, New South Wales
2006, Australia  The University of Sydney Nano Institute, University of
Sydney, Sydney, New South Wales 2006, Australia

###### Abstract

The propagation of chirality across scales is a common but poorly understood
phenomenon in soft matter. In this work, we use computer simulations to study
chiral monolayer assemblies formed by hard rod-like colloidal particles in the
presence of non-adsorbing polymer and characterize the thermodynamic driving
forces responsible for the twisting. Simulations show that straight (achiral)
rods assemble into monolayers with a spontaneous twist that is either left- or
right-handed, while helical (chiral) rods lead to assemblies with preferential
chiral features that depend on their handedness and curliness. The onset of
chirality in these monolayers can be traced back to small clusters formed at
the initial stage of the self-assembly. In these microscopic monolayers,
entropy drives twisting in ways that differ from the assumptions on which
existing continuum theory is built. Depending on the geometry of the
constituent rods, the preferred chiral twist can be driven by entropy gain of
the polymers, or of the rods, or both. In addition, the variation of the
polymer entropy with twist depends on changes in both the surface area and the
volume of the monolayer. Rod fluctuations perpendicular to the monolayer also
play an important role in stabilising the twisting.

##  I Introduction

Colloidal suspensions composed of anisotropic particles can undergo self-
assembly that involves the propagation of chirality from the single-particle
level to the macroscopic level, and so have emerged as a versatile platform
for understanding this common phenomenon in soft matter Barry  _et al._
(2006); Tombolato  _et al._ (2006); Greco and Ferrarini (2015); Dussi and
Dijkstra (2016); Sharma  _et al._ (2020). Cholesteric liquid crystals are a
well-known example, but there are many others. This includes the behaviour of
colloidal suspensions of DNA, viruses, peptides, polysaccharides and various
synthetic nanoparticles Siavashpouri  _et al._ (2017); Miller  _et al._
(2019); Nyström  _et al._ (2018a, b); Aggeli  _et al._ (2001); Lv  _et al._
(2022); Fernández-Rico  _et al._ (2020). While discussion continues about the
physical levers that can be used to control the phenomenon Harris  _et al._
(1999); Sharma  _et al._ (2009); Morrow  _et al._ (2017); Wang  _et al._
(2017); Yeom  _et al._ (2015); Sun  _et al._ (2018); Chiappini  _et al._
(2019), its potential application in areas including optics, catalysis, and
sensing is already being explored Hentschel  _et al._ (2017); Li  _et al._
(2019); Hao  _et al._ (2020), and will likely accelerate in light of advances
in the synthesis of anisotropic and chiral nanoparticles Glotzer and Solomon
(2007); Lee  _et al._ (2018); González-Rubio  _et al._ (2020); Fernández-Rico
_et al._ (2020). It is therefore important to have a better understanding of
the forces which control the propagation of chirality in these systems, and
which can even drive changes in surface topology Khanra  _et al._ (2022).

One typical such colloidal suspension is a mixture consisting of rod-like
particles and non-adsorbing polymers in a good solvent. In these rod-polymer
mixtures, the polymers, which behave as random coils with a radius of gyration
$r_{g}$, can provide an effective attraction between the rods via depletion
forces Asakura and Oosawa (1954), and thus drive the rods to assemble into
diverse ordered structures Lekkerkerker and Stroobants (1994); Siavashpouri
_et al._ (2019). For example, two-dimensional colloidal membranes can form in
a suspension of filamentous viruses and dextran Barry and Dogic (2010). These
colloidal membranes are liquid-like monolayer assemblies, and often have a
round-shaped edge in which the constituent viruses are twisted and exhibit a
chiral distribution of their orientations Gibaud  _et al._ (2012, 2017).

This chiral twist is characteristic of these nearly two-dimensional systems
and is very different from the more common cholesteric twist observed in bulk
(i.e., three dimensional) chiral assembly Straley (1976). The former one is
commonly known as “double twist” to distinguish it from the cholesteric
(single) twist. The double twist cannot be spatially uniform in the bulk and
always occurs with other deformations, with typical examples being the twist-
bend and splay-twist textures Selinger (2018). While the driving mechanism for
the cholesteric twist is relatively well understood, it remains elusive for
the double twist in colloidal membranes and has so far defied a complete
explanation, notwithstanding recent attempts. For instance, an entropically-
motivated continuum theory has been developed to explain the experimental
behaviour of these colloidal membranes Kang  _et al._ (2016a). Briefly, the
entropy, manifested through the viruses as Frank elastic energy for the twist
distortion and through the polymers as an effective surface tension for the
excluded volume, drives the chiral twist of the membranes. This description
further assumes that the membranes are incompressible in the continuum limit.

To test the generality of this theory, and to serve as an important
complementary tool to interpret experimental results, it would be useful to
study such membranes using a particle-based simulation approach. This would be
especially useful for analyzing small clusters formed at the onset of the
self-assembly process where continuum descriptions often break down. To our
knowledge, however, existing simulation studies of twisted membranes have been
limited to the case of achiral rods which lack intrinsic chirality Gibaud  _et
al._ (2012). In this work, we therefore study membranes formed by both achiral
rods and chiral rod-like helices Frezza  _et al._ (2013) using Langevin
dynamics (LD) simulations and characterise the thermodynamic driving forces
responsible for the propagation of chirality in these systems.

##  II Model and Method

###  II.1 Models for rod-polymer suspensions

The rod-polymer suspensions were described using a continuous potential model
that approximates the well-known Asakura-Oosawa-Vrij (AO) model Asakura and
Oosawa (1954); Liu and Widmer-Cooper (2019): (achiral) straight rods,
described as hard spherocylinders, were represented by a rigid linear chain of
length $L$ consisting of overlapping hard spheres of diameter $D$ (Fig. 1a);
chiral rods, described as hard helices, were modeled as a set of hard spheres
having diameter $D$ evenly arranged along a helical line of contour length
$L$, pitch $p$ and radius $r$ (Fig. 1b); and the non-adsorbing polymers were
modelled as spheres with diameter $d=2r_{g}$ that are freely interpenetrable
to each other but experience a hard repulsion from the rod spheres. For
simplicity, we set the diameter of polymer spheres $d=D$. In our simulations,
the hard-core potential between rod-rod (rr) and between rod-polymer (rp)
sphere pairs was replaced by a continuous pseudo-hard-core potential, i.e.,
$U^{\alpha\beta}(r)=50(50/49)^{49}\epsilon[(\sigma/r)^{50}-(\sigma/r)^{49}]$
($\alpha\beta\in\\{\text{rr, rp}\\}$) truncated and shifted at
$r^{\alpha\beta}_{cut}=(50/49)\sigma$, where $r$ is the centre-to-centre
distance between the spheres, $\epsilon$ is the energy parameter, and $\sigma$
is the distance parameter with $\sigma=D$. Besides, for all rods used in this
work, the distance between consecutive spheres is $0.5D$, which is sufficient
to remove side effects associated with surface roughness (see Appendix. A).

While an implicit polymer model for (achiral) straight rods such as that in
Refs. Savenko and Dijkstra (2006); Cherstvy (2008); Patti and Dijkstra (2009)
can allow us to simulate large systems, the corresponding model for helical
rods is lacking and developing an accurate implicit polymer model for hard
helices, especially in the case of large polymers, could be quite challenging
Wood  _et al._ (2021).

###  II.2 Langevin dynamics simulation details

All LD simulations were carried out using LAMMPS Plimpton (1995) at a
dimensionless temperature $k_{B}T/\epsilon=1$ (where $k_{B}$ is the Boltzmann
constant and $T$ is the temperature). In the simulations, rod and polymer
spheres are subjected to three forces: the conservative force $f^{C}$ computed
via the pairwise interactions (i.e., the pseudo-hard-core potential); the
friction force $f^{F}=-(m/\gamma)v$ with $m$ the mass, $\gamma$ the damping
factor, and $v$ the velocity of the sphere; and the random force
$f^{R}\propto\sqrt{k_{B}Tm/(\Delta t\gamma)}$ with $\Delta t$ the time step.
All simulations were performed in a box with periodic boundary conditions. The
velocity-Verlet algorithm was used to integrate the equations of motion with a
time step $\Delta t=0.001\tau$ where $\tau=D\sqrt{m/(k_{B}T)}$, and the
damping factor was set to be $\gamma=1\tau$. In all simulations, we set the
masses of one polymer and one rod $m_{p}=m_{r}=m=1$. Simulations of large
monolayers were performed in an isothermal-isobaric ($NPT$) ensemble. A
Berendsen barostat with a time constant of $1\tau$ was applied. Most
simulations were initialised with $N_{r}=480$ rods in a single hexagonally
packed layer surrounded by $N_{p}=40000$ polymer spheres in a box with initial
dimensions $44\times 44\times 21D^{3}$. Initial configurations with different
chiral twists were also used to confirm that only one handedness was stable.
At least $10$ independent simulations with different initial configurations
were performed for each rod shape, and all simulations were run for at least
$5\times 10^{6}$ steps to collect enough configurations at the equilibrium
state.

###  II.3 Free energy calculations

Simulations used for measuring changes in the free energy
($\Delta\Omega_{total}$) as a function of the twist ($\langle\psi_{i}\rangle$,
see its definition in next section) were performed in a semi-grand canonical
($\mu_{p}VT$) ensemble with $N_{r}=2-61$ rods. During the simulations, $1000$
GCMC insertion and deletion moves were performed every $1000$ LD steps to
maintain the chemical potential of the polymers ($\mu_{p}$). Simulations were
initialised with $N_{r}$ rods in a single hexagonally packed layer surrounded
by $\sim 2000$ polymer spheres in a box with dimensions $15\times 15\times
15D^{3}$. The values of $\Delta\Omega_{total}$ as a function of
$\langle\psi_{i}\rangle$ were evaluated by means of the umbrella sampling (US)
method Torrie and Valleau (1977). We imposed a harmonic spring biasing
potential given by
$U=0.5k[\langle\psi_{i}\rangle-\langle\psi_{i}\rangle_{0}]^{2}$ on the system
using the Colvars package Fiorin  _et al._ (2013). Here, $k$ is the spring
constant, $\langle\psi_{i}\rangle_{0}$ is the desired twist, and
$\langle\psi_{i}\rangle$ is the actual twist in the monolayer. Under the
biasing potential, the monolayer is forced to stay in a pseudo-equilibrium
state with $\langle\psi_{i}\rangle$ fluctuating around
$\langle\psi_{i}\rangle_{0}$. Different twisted states can be described by a
series of values with
$\langle\psi_{i}\rangle_{0}\in(\langle\psi_{i}\rangle_{min},\langle\psi_{i}\rangle_{max})$.
In our simulations, $k=1k_{B}T/deg^{2}$,
$\langle\psi_{i}\rangle_{min}=-24^{\circ}$,
$\langle\psi_{i}\rangle_{max}=24^{\circ}$, and the increment of
$\langle\psi_{i}\rangle_{0}$ was $1^{\circ}$ or $2^{\circ}$. For each given
$\langle\psi_{i}\rangle_{0}$, the system was equilibrated for $2\times 10^{6}$
steps followed by another $2\times 10^{6}$ steps production run in which data
was accumulated every $1000$ LD steps. Finally, the WHAM algorithm Grossfield
was used to calculate the free energy change $\Delta\Omega_{total}$ as a
function of $\langle\psi_{i}\rangle$ . For each monolayer, $10$ independent
simulations were carried out to obtain good statistics. Meanwhile, to prevent
the disassembly of small monolayers, additional spring forces were imposed on
the rod to move it back when the distance from the centre of the rod to the
centre of the monolayer is larger than a critical value $r_{c}$, where
$r_{c}=0.75D$ for $N_{r}=2$, $r_{c}=1.5D$ for $N_{r}=7$, $r_{c}=3.0D$ for
$N_{r}=19$, $r_{c}=4.5D$ for $N_{r}=37$ and $r_{c}=6.0D$ for $N_{r}=61$. These
critical values are larger than the equilibrium radii of the respective stable
monolayers.

###  II.4 Excluded volume calculations

During the production stage of US simulations, we sampled configurations every
$5000$ steps and computed the excluded volume (i.e., $\Delta V_{exc}$) for
polymer spheres due to rods and the corresponding contributions from the
volume and the surface area of the monolayer (i.e., $\Delta V_{exc}^{bulk}$
and $\Delta V_{exc}^{surf}$). For a given configuration, all rod spheres that
had at least one polymer neighbour within $1.5D$ from their centre were
classified as surface rod-spheres, and the rest were classified as bulk rod-
spheres. To compute the excluded volumes, the whole system was divided into
many small cubic bins with an edge length of $l=0.5D$. We confirmed that using
a smaller value of $l$ (e.g., $l=0.25D$, see Appendix. B) gave similar
results. A bin was occupied by rods if there was at least one rod-sphere whose
centre was less than $1.0D$ (corresponding to the polymer diameter) from the
bin’s centre, and the volume of the bin contributed to $\Delta V_{exc}^{surf}$
if all rod-spheres occupying this bin were surface rod-spheres, otherwise it
contributed to $\Delta V_{exc}^{bulk}$. The final value of the excluded volume
at a given $\langle\psi_{i}\rangle$ was averaged over all configurations
collected at the corresponding $\langle\psi_{i}\rangle_{0}$. For these
calculations, the Freud Python package Ramasubramani  _et al._ (2020) was used
to analyse the simulation data.

###  II.5 Suppressing perpendicular fluctuations

In the simulations for monolayers without rod fluctuations perpendicular to
the monolayer, the centres of mass of all rods were constrained on a common
plane via harmonic spring forces using a spring constant of
$1000k_{B}T/D^{2}$.

##  III Results and discussion

Figure 1: Spontaneous twist in monolayers of rods. (a) An achiral straight rod
of length $L$ and diameter $D$. (b) A chiral rod (left-handed helix) with
contour length $L$, diameter $D$, pitch $p$ and radius $r$. $\theta$
quantifies the inclination angle of the rod Frezza  _et al._ (2014). (c-d)
Side-view and top-view of (c) left-handed ($\mathcal{L}$) and (d) right-handed
($\mathcal{R}$) monolayers composed of straight rods with the color indicating
the normalised tilted angle ($\psi_{i}/|\psi_{i}^{max}|$) between the rod axis
$\hat{\textbf{u}}_{i}$ and the nematic director $\hat{\textbf{n}}$.

###  III.1 Spontaneous twist in monolayers

We first considered monolayers composed of (achiral) straight rods with length
of $L=10D$ (Fig. 1a). Figure 1c, d show equilibrium configurations obtained
from simulations with $N_{r}=480$ rods surrounded by $N_{p}=40000$ polymer
spheres at the pressure $P=1.2k_{B}T/D^{3}$. The rods are parallel to the
normal axis (i.e., the nematic director $\hat{\textbf{n}}$) at the centre, but
tilt with increasing magnitude around the radial axis away from the centre. In
multiple independent simulations started with an untwisted configuration, the
monolayer shows nearly equal probability to end up showing left-handed
($\mathcal{L}$) or right-handed ($\mathcal{R}$) twist. Such monolayers, with
roughly square edge profiles, are also predicted by the continuum theory and
have been observed in experiments for small colloidal membranes Kang  _et al._
(2016a).

To evaluate the degree of twist in the monolayers, we used the average tilt
angle of rods with respect to the nematic director, defined as

$\langle\psi_{i}\rangle=\left\langle\dfrac{\hat{\textbf{r}}_{i}\cdot(\hat{\textbf{n}}\times\hat{\textbf{u}}_{i})}{|\hat{\textbf{r}}_{i}\cdot(\hat{\textbf{n}}\times\hat{\textbf{u}}_{i})|}\cos^{-1}(|\hat{\textbf{n}}\cdot\hat{\textbf{u}}_{i}|)\right\rangle,$
(1)

where $\hat{\textbf{r}}_{i}$ is unit vector connecting the center-of-mass of
the monolayer and the center-of-mass of rod $i$, $\hat{\textbf{u}}_{i}$ is the
a unit vector along the long axis of the rod and $\langle\dots\rangle$
indicates an average over all rods in the monolayer and all configurations
collected at the equilibrium state. $\langle\psi_{i}\rangle$ is negative for
the $\mathcal{L}$ twist and positive for the $\mathcal{R}$ twist (Fig. 1c, d).

###  III.2 Phase diagram of chirality in monolayers

We then studied monolayers of left-handed hard helices with $L=10$ and varying
$r$ and $p$ (Fig. 1b). A summary of the results obtained from simulations with
$N_{r}=480$ at $P=1.2k_{B}T/D^{3}$ is reported in Fig. 2. In the phase diagram
(Fig. 2a), we can identify the values of $r$ and $p$ that give rise to a
chiral twist, whose handedness with respect to that of the constituent rods is
(i) the same (e.g., $r=0.1$, $p=2$), (ii) the opposite (e.g., $r=0.1$,
$p=12$), or (iii) mixed with either $\mathcal{R}$ or $\mathcal{L}$ (e.g.,
$r=0.1$, $p=22$) (Fig. 2b).

Figure 2: Phase diagram of chirality in monolayers. (a) The phase diagram
shows the chirality of monolayer assemblies of hard rods with $L=10$ and
varying $r$ and $p$ at $P=1.2$. All symbols are colored according to
$\langle\psi_{i}\rangle$ (Eq. 1) of the corresponding monolayer. The
handedness of the monolayer with respect to that of the constituent hard rods
can be the same (blue square), opposite (red square), or mixed with either
$\mathcal{R}$ or $\mathcal{L}$ (bi-colored square). Note that all symbols with
$r=0$ or $p=\infty$ represent the case of straight rods. Lines (i) and (ii)
indicate the approximate phase boundaries for monolayers. Lines (iii-v),
provided for comparison, indicate the phase boundaries between the same and
opposite regions for bulk cholesteric phases. Line (iii) is given by the
critical inclination angle $\theta=45^{\circ}$ Tombolato  _et al._ (2006);
Cherstvy (2008), while Lines (iv-v) were obtained using density functional
theory at (iv) low and (v) high volume fractions, respectively Belli  _et al._
(2014). (b) Typical snapshots of stable monolayers obtained from simulations
using different left-handed helices.

In Fig. 2a, Line (v) is the phase boundary between same and opposite regimes
for the corresponding cholesteric phases at high volume fractions obtained
using density functional theory Belli  _et al._ (2014). We can see that, in
comparison, the corresponding phase boundary for monolayers [i.e., Line (i)]
is shifted toward larger values of $p$ at $r>0.1$. Such shifting is also
observed in cholesteric phases when the packing density of helices increases
[compare Lines (iii)/(iv) to Line (v) in Fig. 2a]. Thus the difference between
Line (i) and Line (v) is likely due to the higher rod packing fraction in the
monolayers (0.6-0.7) compared to in the bulk cholesteric phases (0.35-0.5)
Belli  _et al._ (2014).

For most helices, however, their monolayer assemblies and cholesteric phases
have the same handedness, supporting the experimental observation of
consistent chirality between the two for rod-shaped viruses Gibaud  _et al._
(2017). For weakly curled helices in the mixed regime, the monolayers can be
either $\mathcal{R}$ or $\mathcal{L}$, while the cholesteric phases may only
exhibit weak opposite handedness Belli  _et al._ (2014). As will be elaborated
later, the driving mechanism of this chiral monolayer assembly is very
different from the bulk cholesteric chiral assembly that was originally
predicted by Straley Straley (1976) and recently confirmed by density
functional theory Frezza  _et al._ (2014); Belli  _et al._ (2014) and
numerical simulations Cinacchi  _et al._ (2017).

From the phase diagram, we also can see that the degree of twist (i.e.,
$\langle\psi_{i}\rangle$) is a non-monotonic function of the intrinsic pitch
of the rods (i.e., $p$), which is consistent with the behaviour of bulk
cholesteric phases formed by hard helices Frezza  _et al._ (2014); Cinacchi
_et al._ (2017). Starting at $p=\infty$ (i.e., straight rods), the magnitude
of $\langle\psi_{i}\rangle$ increases as $p$ decreases, and reaches a maximum
for moderately curled helices (e.g., $r=0.1$, $p=12$ and $r=0.3$, $p=16$) in
the opposite regime, before decreasing to $0$ for helices at the phase
boundary between same and opposite regimes (e.g., $r=0.1$, $p=4$) (Fig. 2b).
In the same regime, $|\langle\psi_{i}\rangle|$ is small, but our results at
$r=0.1$ show that here again $|\langle\psi_{i}\rangle|$ first increases and
then deceases as $p$ decreases (Fig. 2a).

###  III.3 Thermodynamic origins of chiral twist

To study the thermodynamic origins of chiral twist in these monolayer
assemblies, we considered a monolayer of $N_{r}$ rods in a sea of polymer
spheres at fixed volume $V$ and temperature $T$. This system was kept in
osmotic equilibrium with a large reservoir containing the pure polymer
solution at fixed fugacity $z_{p}=\exp(\mu_{p}/k_{B}T)$ where $\mu_{p}$ is the
polymer chemical potential. The grand potential of the system can be written
as $\Omega_{\text{total}}(N_{r},V,T,\mu_{p})=F_{r}-z_{p}(V-V_{exc})k_{B}T$,
where $F_{r}$ is the Helmholtz energy of the rods and $V_{exc}$ is the volume
excluded to the polymers by the hard rods Lekkerkerker and Stroobants (1994);
Bolhuis  _et al._ (1997). The second term on the right is the free energy of
the polymers $\Omega_{p}$. $V_{exc}$ can be further divided into a bulk term
and a surface term associated to the volume and surface area of the monolayer,
respectively. Thus, we obtain the change in free energy expressed as

$\displaystyle\Delta\Omega_{total}$ $\displaystyle=\Delta
F_{r}+\Delta\Omega_{p}$ (2) $\displaystyle=\Delta F_{r}+z_{p}k_{B}T\Delta
V_{exc}$ $\displaystyle=\Delta F_{r}+z_{p}k_{B}T(\Delta V_{exc}^{bulk}+\Delta
V_{exc}^{surf}).$

Both $\Delta F_{r}$ and $\Delta V_{exc}$ depend on the twisting state of rods
in the monolayer. We measured $\Delta\Omega_{total}$ as a function of
$\langle\psi_{i}\rangle$ in a semi-grand canonical ($\mu_{p}VT$) ensemble with
fixed $N_{r}$ at $z_{p}=1.2$ (corresponding to $P=1.2k_{B}T/D^{3}$ in the
previous simulations). The twist was constrained using the US approach, while
$\Delta V_{exc}^{bulk}$ and $\Delta V_{exc}^{surf}$ were numerically
calculated.

We performed a series of US simulations to calculate the changes in free
energy as a function of the twist for monolayers formed by $N_{r}=2-61$ rods
with varying $p$ (see Appendix C). The average twist of stable monolayers
monotonically increases as the monolayer size increases, which is
qualitatively consistent with the theoretical description for small colloidal
membranes in which the tile angles of rods at the edge have yet to reach the
limiting value of $90^{\circ}$ Kang  _et al._ (2016a). More importantly, we
found that simulating tens of rods is sufficient to capture the chiral
behaviour exhibited by the large monolayers shown in Fig. 2, indicating that
the chirality of these monolayers is determined already during the onset of
the self-assembly process.

Figure 3a shows $\Delta\Omega_{total}$ vs. $\langle\psi_{i}\rangle$ for three
typical monolayers made up of $37$ straight rods or left-handed helices. For
the monolayer of straight rods (i.e., $r=0$), the two identical minima at
$\langle\psi_{i}\rangle>0$ and $\langle\psi_{i}\rangle<0$ in the curve of
$\Delta\Omega_{total}$ indicate the stable twist is equally likely to be
$\mathcal{R}$ or $\mathcal{L}$. For the monolayer of helices, only one local
minimum appears in the curve of $\Delta\Omega_{total}$ and is located at
$\langle\psi_{i}\rangle>0$ for the left-handed moderately curled helices
(i.e., $r=0.1$, $p=12$) and at $\langle\psi_{i}\rangle<0$ for the left-handed
highly curled helices (i.e., $r=0.1$, $p=2$), consistent with the behaviour of
the larger monolayers summarised in Fig. 2a.

The decomposition of $\Delta\Omega_{total}$ in Eq. 2 reveals that chiral twist
in these monolayers is stabilised by different driving forces depending on the
rod shape. As shown in Fig. 3a (i), the twist in monolayers of straight rods
is driven by the entropy gain of the polymers with respect to the untwisted
state (i.e., the decease in $\Delta\Omega_{p}$), but further twisting beyond
the equilibrium state is also prevented by the rapidly increasing entropy loss
of the polymers at larger $|\langle\psi_{i}\rangle|$. The rod entropy in this
case shows an almost opposite dependence on $\langle\psi_{i}\rangle$, but the
polymer entropy dominates and stabilises the twist in the monolayer. For the
monolayer of moderately curled helices [Fig. 3a (ii)], the polymer entropy
also dominates and leads to a single stable twist, but now entropy gain from
the rods also contributes. In sharp contrast, for the monolayer of highly
curled helices [Fig. 3a (iii)], the single stable twist is entirely driven by
the rod entropy, competing against the entropy loss of the polymers.

The polymer entropy is related to changes in the volume excluded to polymers
($\Delta V_{exc}$), which is determined by both the volume and the surface
area of the monolayer. Figure 3b shows the change in $\Delta V_{exc}$ and its
volume/surface components when twisting the three monolayers discussed in the
previous paragraph. This reveals that not only the surface area but also the
volume of the monolayers changes significantly during the twisting process.
Especially for monolayers of straight rods [Fig. 3b (i)], the decrease in the
volume acts as the major driving force for twisting.

Figure 3: The thermodynamic origins of chiral twist in monolayers of different
rods. (a) The changes in free energy ($\Delta\Omega_{total}$, $\Delta F_{r}$
and $\Delta\Omega_{p}$) as a function of the twist ($\langle\psi_{i}\rangle$)
for monolayers formed by $N_{r}=37$ (i) straight rods ($r=0.0$, $p=\infty$),
(ii) left-handed moderately curled helices ($r=0.1$, $p=12$), and (iii) left-
handed highly curled helices ($r=0.1$, $p=2$). (b) The corresponding changes
of excluded volume ($\Delta V_{exc}$, $\Delta V_{exc}^{bulk}$ and $\Delta
V_{exc}^{surf}$).

###  III.4 Comparison with continuum theory

Having shown that different entropy components can drive rod monolayers to
twist, we now compare our results with the continuum theory developed to
describe such colloidal membranes Kang  _et al._ (2016a). The continuum theory
is based on a relatively simple physical picture (see Appendix D): that the
twist is driven mainly by the entropy gained by the polymers when the membrane
surface area is minimised at constant membrane volume. In this model, the
polymer entropy is invariant under chirality inversion and does not contribute
to the preference of the handedness, regardless of the chirality of the rods.
The preferred handedness is instead attributed to an entropy term in the Frank
elastic energy of the rods, whose magnitude depends on the preferred twist
wavenumber that implicitly contains the chiral features of the rods.

In contrast, our simulation results reveal the existence of more complex
thermodynamic behaviour. First, while the polymer entropy often drives
twisting, it can also oppose twisting, with the rod entropy instead driving
twisting in those cases [e.g., Fig 3a (iii)]. Second, the polymer entropy is
asymmetric under chirality inversion for monolayers of helical particles and
contributes to the preference of the handedness in these cases [e.g., Fig 3a
(ii)]. This indicates that, at best, the Frank elastic energy in the continuum
theory can depend on polymer concentration. Third, the constant-volume
assumption in the continuum theory clearly breaks down, at least for the small
assemblies considered here, indicating that the variation of the polymer
entropy involves contributions from not only the surface area but also the
volume of the monolayer (Fig. 3b).

Recent theoretical and experimental work Chaturvedi and Kamien (2020); Kamien
and Machon (2020); Miller  _et al._ (2020) has also concluded that the volume
change upon twist plays a crucial role in determining the geometry and
stability of colloidal membranes of rod-like particles. The geometric
frustration between double-twist and splay causes the twisted monolayer to
have a hyperbolic edge (i.e., "splay-twist" texture Chaturvedi and Kamien
(2020)), and the splay of rods away from the monolayer midplane leads to a
local volume expansion, which is most significant at the top and the bottom of
the monolayer edge. Based on this geometric argument, using a combination of
experiments and theory in which the variation of rod density is considered,
Miller et al. Miller  _et al._ (2020) demonstrated that for the colloidal
rafts in the membranes composed of rigid rods of different lengths, the splay
deformation causes expansion and compression of the inner and outer raft
edges,respectively, and their competition results in spontaneous twist even
for achiral systems and non-monotonic dependence of the stable twist as a
function of the raft size. As for our simulated monolayers which are assembled
from monodisperse rods in non-adsorbing polymers, the diffuse interfacial
region exhibits a clear decline of the rod density away from the midplane,
especially at the edges (see Fig. E.1a in Appendix E), but we did not observe
an obvious hyperboloid-like shape. Only when the perpendicular fluctuations of
the rods are suppressed and their centers are confined at the 2D midplane
(which is the same as the theoretical model in Refs. Miller  _et al._ (2020)),
does our monolayer exhibit a hyperboloid-like shape (see Fig. E.1b in Appendix
E). Moreover, in our systems, the total volumes of small monolayers (which are
made up of $37$ rods in our US simulations) are always decreasing during the
initial twisting process (see Fig. E.2 Appendix E), suggesting that the volume
expansion due to the splay deformation does not dominate, at least in these
small monolayers. All these results suggest that the volume change due to
twist is important to the stable texture in various colloidal membranes of
rod-like particles and thus should not be ignored in theoretical models.

Figure 4: The role of perpendicular fluctuations. (a) The changes in free
energy ($\Delta\Omega_{total}$, $\Delta F_{r}$ and $\Delta\Omega_{p}$) as a
function of the twist ($\langle\psi_{i}\rangle$) for monolayers formed by
$N_{r}=37$ (i) straight rods ($r=0.0$, $p=\infty$), (ii) left-handed
moderately curled helices ($r=0.1$, $p=12$), and (iii) left-handed highly
curled helices ($r=0.1$, $p=2$) when the fluctuations perpendicular to the
plane of the monolayer are artificially suppressed. (b) Typical snapshots of
large equilibrated monolayers formed by $N_{r}=480$ rods (corresponding to
those in a), obtained from simulations without rod fluctuations perpendicular
to the monolayer.

###  III.5 The role of perpendicular fluctuations

Finally, we considered a special entropy contribution from rods related to
their fluctuations perpendicular to the monolayer. Figure 4c shows
$\Delta\Omega_{total}$ vs. $\langle\psi_{i}\rangle$ for the three example
monolayers when the centres-of-mass of all rods are constrained to the
midplane of the monolayer. The rod fluctuations out of the plane are expected
to produce surface roughness and so increase the volume excluded to the
polymers. We found that suppressing the fluctuations resulted in more entropy
gain for the polymers upon twisting for monolayers of straight rods and
moderately curled helices [see larger changes of $\Delta\Omega_{p}$ in Fig. 4a
(i) and (ii) compared to that in Fig. 3a (i) and (ii)]. This stabilises the
twisted states for these monolayers, and even adds a new metastable twisted
state for the monolayer of moderately curled helices [Fig. 4a (ii)]. In
contrast, for the monolayer of highly curled helices, the polymer entropy
increases dramatically upon twisting when the fluctuations are suppressed,
causing the original weakly-stable twisted state to disappear [Fig. 4a (iii)].
These results are consistent with unconstrained simulations of large
monolayers (Fig. 4b), and clearly show that rod fluctuations perpendicular to
the monolayer have important effects on the stability of the chiral twists
that depend on the shape of the individual rods.

We note that such contributions from rod fluctuations perpendicular to the
monolayer are either ignored in continuum models of colloidal membranes or
only taken into account in a simplistic manner (which is not curliness-
dependent) Kang  _et al._ (2016a) (see Appendix D). Our simulation results,
however, indicate that these fluctuations can play a crucial role in the
stability of the chiral twist, and thus may need to be accurately described in
order to predict the stable chiral twist.

##  IV Conclusions

In summary, we have used a simple model to characterise spontaneous chiral
twist in monolayers assembled from either achiral or chiral rods in non-
adsorbing polymer solutions, and thus to reveal the thermodynamic driving
forces responsible for the chirality propagation from single particles to
their assemblies. Note that the chiral twist discussed in this work is the
double twist, which is distinct from the cholesteric (single) twist Selinger
(2018).

Depending on the geometry of the constituent rods, their monolayer assemblies
exhibit a broad range of chiral behaviour, including variations in handedness
and twist magnitude. Compared to the constituent rods, the (achiral) straight
rods and weakly curled helices form monolayers with either $\mathcal{R}$ or
$\mathcal{L}$ twist (i.e., the mixed regime), moderately curled helices form
monolayers with the opposite handedness (i.e., the opposite regime), and
highly curled helices form monolayers with the same handedness (i.e., the same
regime). Moreover, the degree of twist in the monolayers is a non-monotonic
function of the intrinsic pitch of the helices, with the most twisted
monolayers forming from moderately curled helices [Fig. 2a].

The thermodynamic forces responsible for spontaneous chiral twist also vary
dramatically between different particle shapes. In the mixed and the opposite
regimes, the twist in monolayers is mainly driven by the polymer entropy [Fig.
3a (i)]. As the rods becomes more curled, the rod entropy also contributes to
the twist, and only the twisted state with the opposite handedness remains
stable [Fig. 3a (ii)]. For even more curled helices, only one weakly twisted
state with the same handedness is stable, and is entirely driven by the rod
entropy [Fig. 3a (iii)].

Our simulation results also indicate important contributions from the volume
change upon twist and the rod fluctuations perpendicular to the monolayer that
have so far been ignored in continuum theories. Our preliminary results,
obtained from Monte Carlo (MC) simulations, for rods held together by explicit
attraction rather than polymer depletion indicate a similar complexity (see
Appendix F). Overall, we find increasing deviations from current continuum
theory as the attractive forces holding the rods together become weaker,
regardless of whether they are due to direct energetic or indirect entropic
effects (see Appendix F). All these results contribute to our understanding of
chirality transmission across scales when chiral objects assemble into larger
aggregates.

While our simulations were based on a simplified model for rod-like colloids
(i.e., hard spherocylinders and helices), they clearly show that twisted
colloidal membranes can also be formed by helical rods, which could be an
interesting behaviour to investigate in future experiments by using similar
natural and synthetic particles Barry  _et al._ (2006); Zhu  _et al._ (2014);
Feng  _et al._ (2017); Tao  _et al._ (2019). Our current work also offers a
helpful reference for understanding the behaviour in more complex systems. It
would be very useful to consider models which are closer to the chiral rods
(e.g. fd-virus and DNA origami rods) used in experiments of colloidal
membranes. For example, using a “straight and helically-decorated” model
Tombolato  _et al._ (2006); Wensink and Ferreiro-Córdova (2017) would allow us
to compare the computational and experimental results more directly.
Meanwhile, recent experiments have shown that the shape fluctuation of chiral
rods also dramatically affects their assembled structures Tortora  _et al._
(2020), thus it would be interesting to consider the flexibility of rod-like
particles in future work.

###### Acknowledgements.

This work was supported by the Australian Research Council under Grant
CE170100026 and by the MIUR PRIN-COFIN2017 Soft Adaptive Networks grant
2017Z55KCW. Computational resources were provided by the Sydney Informatics
Hub, a Core Research Facility of the University of Sydney.

## References

  * Barry  _et al._ (2006) E. Barry, Z. Hensel, Z. Dogic, M. Shribak, and R. Oldenbourg, Phys. Rev. Lett. 96, 018305 (2006).
  * Tombolato  _et al._ (2006) F. Tombolato, A. Ferrarini, and E. Grelet, Phys. Rev. Lett. 96, 258302 (2006).
  * Greco and Ferrarini (2015) C. Greco and A. Ferrarini, Phys. Rev. Lett. 115, 147801 (2015).
  * Dussi and Dijkstra (2016) S. Dussi and M. Dijkstra, Nat. Commun. 7, 11175 (2016).
  * Sharma  _et al._ (2020) A. Sharma, J. P. Wojciechowski, Y. Liu, T. Pelras, C. M. Wallace, M. Müllner, A. Widmer-Cooper, P. Thordarson, and G. Lakhwani, Cell Reports Phys. Sci. 1, 100148 (2020).
  * Siavashpouri  _et al._ (2017) M. Siavashpouri, C. H. Wachauf, M. J. Zakhary, F. Praetorius, H. Dietz, and Z. Dogic, Nat. Mater. 16, 849 (2017).
  * Miller  _et al._ (2019) J. M. Miller, C. Joshi, P. Sharma, A. Baskaran, A. Baskaran, G. M. Grason, M. F. Hagan, and Z. Dogic, Proc. Natl. Acad. Sci. 116, 15792 (2019).
  * Nyström  _et al._ (2018a) G. Nyström, M. Arcari, and R. Mezzenga, Nat. Nanotechnol. 13, 330 (2018a).
  * Nyström  _et al._ (2018b) G. Nyström, M. Arcari, J. Adamcik, I. Usov, and R. Mezzenga, ACS Nano 12, 5141 (2018b).
  * Aggeli  _et al._ (2001) A. Aggeli, I. A. Nyrkova, M. Bell, R. Harding, L. Carrick, T. C. B. McLeish, A. N. Semenov, and N. Boden, Proc. Natl. Acad. Sci. 98, 11857 (2001).
  * Lv  _et al._ (2022) J. Lv, X. Gao, B. Han, Y. Zhu, K. Hou, and Z. Tang, Nat. Rev. Chem. 6, 125 (2022).
  * Fernández-Rico  _et al._ (2020) C. Fernández-Rico, M. Chiappini, T. Yanagishima, H. de Sousa, D. G. A. L. Aarts, M. Dijkstra, and R. P. A. Dullens, Science 369, 950 (2020).
  * Harris  _et al._ (1999) A. B. Harris, R. D. Kamien, and T. C. Lubensky, Rev. Mod. Phys. 71, 1745 (1999).
  * Sharma  _et al._ (2009) V. Sharma, M. Crne, J. O. Park, and M. Srinivasarao, Science 325, 449 (2009).
  * Morrow  _et al._ (2017) S. M. Morrow, A. J. Bissette, and S. P. Fletcher, Nat. Nanotechnol. 12, 410 (2017).
  * Wang  _et al._ (2017) P.-P. Wang, S.-J. Yu, A. O. Govorov, and M. Ouyang, Nat. Commun. 8, 14312 (2017).
  * Yeom  _et al._ (2015) J. Yeom, B. Yeom, H. Chan, K. W. Smith, S. Dominguez-Medina, J. H. Bahng, G. Zhao, W. S. Chang, S. J. Chang, A. Chuvilin, D. Melnikau, A. L. Rogach, P. Zhang, S. Link, P. Král, and N. A. Kotov, Nat. Mater. 14, 66 (2015).
  * Sun  _et al._ (2018) J. Sun, Y. Li, F. Yan, C. Liu, Y. Sang, F. Tian, Q. Feng, P. Duan, L. Zhang, X. Shi, B. Ding, and M. Liu, Nat. Commun. 9, 2599 (2018).
  * Chiappini  _et al._ (2019) M. Chiappini, T. Drwenski, R. van Roij, and M. Dijkstra, Phys. Rev. Lett. 123, 068001 (2019).
  * Hentschel  _et al._ (2017) M. Hentschel, M. Schäferling, X. Duan, H. Giessen, and N. Liu, Sci. Adv. 3, 1 (2017).
  * Li  _et al._ (2019) S. Li, J. Liu, N. S. Ramesar, H. Heinz, L. Xu, C. Xu, and N. A. Kotov, Nat. Commun. 10, 4826 (2019).
  * Hao  _et al._ (2020) C. Hao, L. Xu, H. Kuang, and C. Xu, Adv. Mater. 32, 1 (2020).
  * Glotzer and Solomon (2007) S. C. Glotzer and M. J. Solomon, Nat. Mater. 6, 557 (2007).
  * Lee  _et al._ (2018) H. E. Lee, H. Y. Ahn, J. Mun, Y. Y. Lee, M. Kim, N. H. Cho, K. Chang, W. S. Kim, J. Rho, and K. T. Nam, Nature 556, 360 (2018).
  * González-Rubio  _et al._ (2020) G. González-Rubio, J. Mosquera, V. Kumar, A. Pedrazo-Tardajos, P. Llombart, D. M. Solís, I. Lobato, E. G. Noya, A. Guerrero-Martínez, J. M. Taboada, F. Obelleiro, L. G. MacDowell, S. Bals, and L. M. Liz-Marzán, Science 368, 1472 (2020).
  * Khanra  _et al._ (2022) A. Khanra, L. L. Jia, N. P. Mitchell, A. Balchunas, R. A. Pelcovits, T. R. Powers, Z. Dogic, and P. Sharma, Proc. Natl. Acad. Sci. 119, e2204453119 (2022).
  * Asakura and Oosawa (1954) S. Asakura and F. Oosawa, J. Chem. Phys. 22, 1255 (1954).
  * Lekkerkerker and Stroobants (1994) H. N. W. Lekkerkerker and A. Stroobants, Nuovo Cim. D 16, 949 (1994).
  * Siavashpouri  _et al._ (2019) M. Siavashpouri, P. Sharma, J. Fung, M. F. Hagan, and Z. Dogic, Soft Matter 15, 7033 (2019).
  * Barry and Dogic (2010) E. Barry and Z. Dogic, Proc. Natl. Acad. Sci. 107, 10348 (2010).
  * Gibaud  _et al._ (2012) T. Gibaud, E. Barry, M. J. Zakhary, M. Henglin, A. Ward, Y. Yang, C. Berciu, R. Oldenbourg, M. F. Hagan, D. Nicastro, R. B. Meyer, and Z. Dogic, Nature 481, 348 (2012).
  * Gibaud  _et al._ (2017) T. Gibaud, C. N. Kaplan, P. Sharma, M. J. Zakhary, A. Ward, R. Oldenbourg, R. B. Meyer, R. D. Kamien, T. R. Powers, and Z. Dogic, Proc. Natl. Acad. Sci. 114, E3376 (2017).
  * Straley (1976) J. P. Straley, Phys. Rev. A 14, 1835 (1976).
  * Selinger (2018) J. V. Selinger, Liq. Cryst. Rev. 6, 129 (2018).
  * Kang  _et al._ (2016a) L. Kang, T. Gibaud, Z. Dogic, and T. C. Lubensky, Soft Matter 12, 386 (2016a).
  * Frezza  _et al._ (2013) E. Frezza, A. Ferrarini, H. B. Kolli, A. Giacometti, and G. Cinacchi, J. Chem. Phys. 138, 164906 (2013).
  * Liu and Widmer-Cooper (2019) Y. Liu and A. Widmer-Cooper, J. Chem. Phys. 150, 244508 (2019).
  * Savenko and Dijkstra (2006) S. V. Savenko and M. Dijkstra, J. Chem. Phys. 124, 234902 (2006).
  * Cherstvy (2008) A. G. Cherstvy, J. Phys. Chem. B 112, 12585 (2008).
  * Patti and Dijkstra (2009) A. Patti and M. Dijkstra, Phys. Rev. Lett. 102, 128301 (2009).
  * Wood  _et al._ (2021) J. A. Wood, Y. Liu, and A. Widmer-Cooper, J. Chem. Phys. 154, 244505 (2021).
  * Plimpton (1995) S. Plimpton, J. Comput. Phys. 117, 1 (1995).
  * Torrie and Valleau (1977) G. Torrie and J. Valleau, J. Comput. Phys. 23, 187 (1977).
  * Fiorin  _et al._ (2013) G. Fiorin, M. L. Klein, and J. Hénin, Mol. Phys. 111, 3345 (2013).
  * (45) A. Grossfield, “Wham: the weighted histogram analysis method,” Version 2.0.9.1.
  * Ramasubramani  _et al._ (2020) V. Ramasubramani, B. D. Dice, E. S. Harper, M. P. Spellings, J. A. Anderson, and S. C. Glotzer, Comput. Phys. Commun. 254, 107275 (2020).
  * Frezza  _et al._ (2014) E. Frezza, A. Ferrarini, H. Bindu Kolli, A. Giacometti, and G. Cinacchi, Phys. Chem. Chem. Phys. 16, 16225 (2014).
  * Belli  _et al._ (2014) S. Belli, S. Dussi, M. Dijkstra, and R. van Roij, Phys. Rev. E 90, 020503(R) (2014).
  * Cinacchi  _et al._ (2017) G. Cinacchi, A. Ferrarini, A. Giacometti, and H. B. Kolli, J. Chem. Phys. 147, 224903 (2017).
  * Bolhuis  _et al._ (1997) P. G. Bolhuis, A. Stroobants, D. Frenkel, and H. N. W. Lekkerkerker, J. Chem. Phys. 107, 1551 (1997).
  * Chaturvedi and Kamien (2020) N. Chaturvedi and R. D. Kamien, Proc. R. Soc. A Math. Phys. Eng. Sci. 476, 20190824 (2020).
  * Kamien and Machon (2020) R. D. Kamien and T. Machon, Proc. Natl. Acad. Sci. 117, 24102 (2020).
  * Miller  _et al._ (2020) J. M. Miller, D. Hall, J. Robaszewski, P. Sharma, M. F. Hagan, G. M. Grason, and Z. Dogic, Sci. Adv. 6, eaba2331 (2020).
  * Zhu  _et al._ (2014) Y. Zhu, J. He, C. Shang, X. Miao, J. Huang, Z. Liu, H. Chen, and Y. Han, J. Am. Chem. Soc. 136, 12746 (2014).
  * Feng  _et al._ (2017) W. Feng, J.-Y. Kim, X. Wang, H. A. Calcaterra, Z. Qu, L. Meshi, and N. A. Kotov, Sci. Adv. 3, e1601159 (2017).
  * Tao  _et al._ (2019) X. Tao, H. Li, B. Yu, X. Wu, Y. Lu, Y. Wang, and H. Chen, Nanoscale 11, 19729 (2019).
  * Wensink and Ferreiro-Córdova (2017) H. H. Wensink and C. Ferreiro-Córdova, Soft Matter 13, 3885 (2017).
  * Tortora  _et al._ (2020) M. M. C. Tortora, G. Mishra, D. Prešern, and J. P. K. Doye, Sci. Adv. 6, eaaw8331 (2020).
  * Meijer and Frenkel (1994) E. J. Meijer and D. Frenkel, J. Chem. Phys. 100, 6873 (1994).
  * Shirts and Chodera (2008) M. R. Shirts and J. D. Chodera, J. Chem. Phys. 129, 124105 (2008), 0801.1426 .
  * Kang  _et al._ (2016b) L. Kang, T. Gibaud, Z. Dogic, and T. C. Lubensky, Soft Matter 12, 386 (2016b).

##  Appendix A Effect of the number of spheres in one rod

For all results reported in the main text, the contour length $L=10$ and there
are $21$ spheres evenly distributed in each rod (i.e. the distance between
consecutive spheres is $0.5D$). We confirmed that this model is smooth enough
to remove side effects associated with grooves between overlapping spheres.
Fig. A.1 shows that using more spheres in each rod yields similar results.

Figure A.1: The total free energy ($\Delta\Omega_{total}$) and its
decomposition ($\Delta F_{r}$ and $\Delta\Omega_{p}$)as a function of the
degree of the twist ($\langle\psi_{i}\rangle$) for a monolayer formed by
$N_{r}=37$ straight rods with (a) $21$ and (b) $41$ spheres in each rod.

##  Appendix B Effect of the bin size to determine the excluded volume

To compute the excluded volumes during the umbrella sampling simulations, the
whole system was divided into many small cubic bins with an edge length of
$l=0.5D$. We confirmed that using a smaller value of $l$ (e.g., $l=0.25D$,
Fig. B.1) gave similar results.

Figure B.1: The change in the total excluded volume ($\Delta V_{exc}$) as a
function of $\langle\psi_{i}\rangle$ and its decomposition into bulk and
surface contributions ($\Delta V_{exc}^{bulk}$ and $\Delta V_{exc}^{surf}$)
for monolayers of $37$ straight rods. In the calculations, the edge length of
the cubic volumetric bins was $0.5D$ in (a) and $0.25D$ in (b).

##  Appendix C Free energy changes for monolayers composed of different
numbers of rods

Figure C.1: The total free energy change ($\Delta\Omega_{total}$) as a
function of the twist ($\langle\psi_{i}\rangle$) for monolayers formed by
$N_{r}=2-61$ rods with $L=10$, $r=0.1$ and varying $p$.

We performed a series of umbrella sampling simulations to calculate the change
in free energy as a function of the twist for monolayers formed by
$N_{r}=2-61$ rods with $L=10$, $r=0.1$ and varying $p$. The obtained results
shown in Fig. C.1 suggest that when $N_{r}>=19$, the free energy minima in
these curves are consistent with the phase behaviours reported in Fig. 2 (a)
in the main text, as discussed in the main text.

Figure C.2: The change in Gibbs free energy ($\Delta G$) as a function of the
twist ($\langle\psi_{i}\rangle$) for monolayers formed by $N_{r}=61$ straight
rods (Blue), $N_{r}=120$ (Orange), $N_{r}=240$ (Green), and $N_{r}=480$ (Red).
Shaded areas indicate one standard error of the mean.

We also examined the free energy of twisting for larger monolayers. For
computational efficiency, this work was performed using Monte Carlo (MC)
simulations of the standard hard AO interaction modelWood  _et al._ (2021).
The same rod length and polymer size was used as in the rest of the paper,
with natural twisting observed at a polymer packing fraction of $0.467$ and
pressure of $1.2$. Staged umbrella samplingMeijer and Frenkel (1994) with the
same bias as in the main text was then used to determine the change in Gibbs
free energy as a function of the average tilt angle for several different
monolayer sizes, with the free energy profiles assembled using the MBAR
techniqueShirts and Chodera (2008). The results are shown in Fig. C.2 and are
consistent with the results presented in Fig. C.1 (at bottom right) for
$N_{r}=61$.

There are two notable changes as the size of the monolayer increases. First,
the magnitude of the preferred average twist increases with the size of the
cluster, despite the fraction of rods at the edge of the cluster shrinking
with increasing cluster size. This indicates that the local twist at the edge
has yet to reach its limiting value of 90°. Once that occurs, we expect the
average twist to reach a maximum value and then to gradually decrease. Second,
the depth of the free energy minimum, and consequently the free energy barrier
to reversing twist direction, increases with size. This increase is roughly
linear with the number of rods and monolayer area.

##  Appendix D Brief review of the continuum theory for colloidal membranes

In the continuum theory Kang  _et al._ (2016b), the colloidal membrane is
treated as a continuum medium composed of rods at constant density, and the
membrane has a fixed number of rods and a constant volume. In the membrane, a
rod at x (i.e. a given position in the membrane plane) is tilted by
$\theta(\textbf{x})$ with respect to the normal director. The membrane half-
thickness can be written as
$h(\textbf{x})=t\cos\theta(\textbf{x})+b(\textbf{x})$ where $t$ is the half-
length of the rod and $b(\textbf{x})$ is the height fluctuation amplitude of
the rod. Assuming a circularly-symmetric membrane of radius $R$ and using
cylindrical coordinates, one can have $h(r)$, $b(r)$, and $\theta(r)$ that
only depend on the radial coordinate.

The free energy associated with the twist distortion of the rods is described
by the Frank elastic free energy. Using the one-constant approximation, the
free energy is given by

$F_{\text{Frank}}=2\pi K\int_{0}^{R}drh\left[r(\partial_{r}\theta)^{2}+\sin
2\theta\partial_{r}\theta+\frac{\sin^{2}\theta}{r}-2qr\partial_{r}\theta-q\sin
2\theta\right],$ (3)

where $K$ is the Frank elastic constant and $q$ is the preferred twist
wavenumber associated with the intrinsic chirality of the constituent rods.

The free energy of the polymers is related to the volume excluded to them by
the rods. For polymers small compared to the dimensions of the membrane, this
excluded volume is approximately $V_{0}+aA$, where $V_{0}$ is the volume of
the membrane, $A$ is the surface area of the membrane, and $a$ is the polymer
characteristic radius. $V_{0}$ is assumed to be constant, so the corresponding
free energy is similar to an effective surface tension, which is expressed as

$F_{\text{polymer}}=4\pi
nak_{B}T\left[\int_{0}^{R}drr\sqrt{1+(\partial_{r}h)^{2}}+Rh(R)\right],$ (4)

where $n$ is the polymer concentration.

If the rod fluctuations perpendicularly to the membrane plane are ignored,
then $h\approx t\cos\theta$, and the profile of the membrane is can be
obtained by minimizing the total free energy over $h(r)$ using a volume-
conserving Lagrange multiplier $\lambda$, i.e.,

$F=F_{\text{Frank}}+F_{\text{polymer}}+\lambda\left[V_{0}-4\pi\int_{0}^{R}drrh\right],$
(5)

In this continuum theory, the key points for the twist are:

  * •

The preferred handedness in membranes composed of chiral rods is entirely
determined by the non-zero $q$ through the free energy of the rods. For
example, when $q>0$, twisted membranes with $\partial_{r}\theta>0$ have lower
energy than those with $\partial_{r}\theta<0$.

  * •

The polymer entropy drives the twist via minimizing the excluded volume, but
it is invariant under the chirality inversion $\theta\rightarrow-\theta$, and
thus does not contribute to the preferred handedness regardless of the
intrinsic chirality of the rods (i.e. $q$).

  * •

The excluded volume only depends on the surface area of the membrane since the
volume of the membrane is assumed to be constant.

The rod fluctuations perpendicular to the membrane plane have complicated,
non-linear effects on the free energy, and so are very difficult to accurately
account for in the theory. The continuum theory in Ref. Kang  _et al._ (2016b)
only considered fluctuations of single rods and thus ignored their
interactions and correlated motion. The corresponding free energy was
calculated in the small rod angle and small fluctuation amplitude limit, and
the rods were assumed to be packed hexagonally and to maintain a constant
perpendicular distance $\xi$ between nearest-neighbors. Under these
assumptions, the rod fluctuation free energy is given by

$F_{\text{fluctuation}}=\frac{8\pi^{2}nak_{B}T}{\sqrt{3}\xi^{2}}\int_{0}^{R}drr\cos\theta\left[h-t\cos\theta-(2\pi
na)^{-1/2}\right].$ (6)

##  Appendix E Density profiles and volume changes of assembled monolayers

Figure E.1: The snapshots and the density profile for (a) a free monolayer
suspended in nonadsorbing polymers (b) a confined monolayer with centers of
all rods fixed at the 2D midplane. The dashed lines are density contours at
$\rho=0.5$, $1.0$, and $1.5$ $D^{-3}$. The slightly concave contour lines in
(b) are the results of the hyperboloid-like shape. All monolayers are made up
of (achiral) straight rods. Figure E.2: The changes of excluded volume for the
two monolayers made up of $37$ straight rods. (a) A free monolayer suspended
in nonadsorbing polymers. (b) A confined monolayer with centers of all rods
fixed at the 2D midplane.

##  Appendix F Monolayers held together by explicit attraction

We observed similar twisting behavior in the absence of polymer depletion when
the rods interacted with each other via an explicit attractive potential. In
this model, the rods consisted of $40$ spheres over a length $L=10$ with the
spheres in different rods interacting with each other via a square well
potential of width $w=1$ and depth $1/20^{2}=0.0025\epsilon_{0}$, and the
temperature scaled to $T=k_{B}T/\epsilon_{0}$. There was no polymer in the
model, with rod-attraction coming only from the square well interactions. All
work was performed using monolayers of 61 rods. Using the same MC simulation
and umbrella sampling approach described in section C, this time in the $NVT$
ensemble, we calculated the Helmholtz free energy change ($\Delta F$) along
with the change in the mean potential energy ($\Delta U$) as a function of the
average twist. Using this, we also determined the entropy contribution to
$\Delta F$ using $-T\Delta S=\Delta F-\Delta U$.

Figure F.1: The change in free energy $\Delta F$ (Blue) and its decomposition
into $-T\Delta S$ (Orange) and $\Delta U$ (Green) as a function of twist for
monolayers formed by $N_{r}=61$ square-well rods at temperatures of a)
$T=0.1$, b) $T=0.12$, and c) $T=0.16$. d) Shows $\Delta U$ normalized by the
number of rods for all rods in the monolayer (Green), for rods located at the
edge of the monolayer (Red), and for rods located in the interior (Purple).
Shaded areas indicate one standard error of the mean.

The energy changes obtained as a function of the twist angle at three
different temperatures are shown in Fig. F.1 (a-c). At the lower temperatures
($T=0.1$ and $T=0.12$), the initial twisting is driven entirely by the
potential energy of the rods, with the rod entropy only becoming a significant
driver at higher angles. This is similar to what we observed for straight rods
when the twisting was driven by depletion-attraction [Fig. 3a (i) in the main
text]. In that case, the rod free energy $\Delta F_{r}$, which is analogous to
$-T\Delta S$ in this model, opposes twisting at small angles and only drives
it at larger angles. Similarly, the force holding the monolayer together in
both models ($\Delta\Omega_{p}$ and $\Delta U$) favors small twists but
opposes larger ones.

Breaking down the potential energy into contributions from (i) rods in the
interior of the monolayer and (ii) those at the edge, provides insight into
whether the initial twisting is driven by minimization of the surface energy
of the monolayer. This breakdown is shown for $T=0.12$ in Fig. F.1 (d).
Surprisingly, we find very little difference in the potential energy change
per rod between the different subgroups until past the optimum twist angle,
counter-intuitively indicating that the edge of the monolayer, where the
twisting is the greatest and there are fewer interactions per rod, is not more
energetically advantaged or disadvantaged by the twisting than the rest of the
cluster.

At higher temperature [$T=0.16$, Fig. F.1 (c)], the driving forces change,
with the rod entropy now driving the twist at all angles and the potential
energy always opposing the twist. This indicates that twisting becomes
increasingly driven by the rod entropy as the monolayer density decreases.
Consistent with this, we observe a similar change for depletion-driven
twisting as the polymer fugacity decreases, i.e. with the rod entropy $F_{r}$
increasingly favoring twisting at small twists (e.g., $\mu_{p}=0.8$) in the
depletion-driven case of straight rods [see Fig. F.2 compared to Fig. 3a in
the main text].

Figure F.2: The change in total free energy ($\Delta\Omega_{total}$) and its
decomposition ($\Delta F_{r}$ and $\Delta\Omega_{p}$) as a function of the
twist ($\langle\psi_{i}\rangle$) for a monolayer formed by $N_{r}=37$ straight
rods at a polymer fugacity of (a) $z_{p}=1.0$ and (b) $z_{p}=0.8$.


# Quantum phase transition in skewed ladders: an entanglement entropy
and fidelity study

Sambunath Das  Solid State and Structural Chemistry Unit, Indian Institute of
Science, Bangalore 560012, India  S. N. Bose National Centre for Basic
Sciences, Block - JD, Sector - III, Salt Lake, Kolkata - 700106, India
Dayasindhu Dey  Solid State and Structural Chemistry Unit, Indian Institute of
Science, Bangalore 560012, India  UGC-DAE Consortium for Scientific Research,
University Campus, Khandwa Road, Indore 452001, India  S. Ramasesha  Solid
State and Structural Chemistry Unit, Indian Institute of Science, Bangalore
560012, India  Manoranjan Kumar  S. N. Bose National Centre for Basic
Sciences, Block - JD, Sector - III, Salt Lake, Kolkata - 700106, India

###### Abstract

Entanglement entropy (EE) of a state is a measure of correlation or
entanglement between two parts of a composite system and it may show
appreciable change when the ground state (GS) undergoes a qualitative change
in a quantum phase transition (QPT). Therefore, the EE has been extensively
used to characterise the QPT in various correlated Hamiltonians. Similarly
fidelity also shows sharp changes at a QPT. We characterized the QPT of
frustrated antiferromagnetic Heisenberg spin-1/2 systems on 3/4, 3/5 and 5/7
skewed ladders using the EE and fidelity analysis. It is noted that all the
non-magnetic to magnetic QPT boundary in these systems can be accurately
determined using the EE and fidelity, and the EE exhibits a discontinuous
change, whereas fidelity shows a sharp dip at the transition points. It is
also noted that in case of the degenerate GS, the unsymmetrized calculations
show wild fluctuations in the EE and fidelity even without actual phase
transition, however, this problem is resolved by calculating the EE and the
fidelity in the lowest energy state of the symmetry subspaces, to which the
degenerate states belong.

##  I Introduction

Quantum phase transition (QPT) involves a qualitative change in the nature of
the ground state (GS) in the phase space of the Hamiltonian parameters. The
prerequisite for a Hamiltonian to exhibit QPT is that it should consist of
noncommuting terms. In this case the quantum fluctuations drive the QPT. In
1-D and 2-D Hamiltonians the fluctuations are dominant, making them viable
candidates to exhibit QPT. Many interesting and exotic quantum phases observed
include the dimer phase Majumdar and Ghosh (1969); Chubukov (1991); Chitra
_et al._ (1995); White and Affleck (1996); Itoi and Qin (2001); Mahdavifar
(2008); Sirker (2010); Kumar  _et al._ (2015); Soos  _et al._ (2016); Kumar
_et al._ (2010), spin liquid phase Hamada  _et al._ (1988); White and Affleck
(1996); Itoi and Qin (2001); Mahdavifar (2008); Kumar  _et al._ (2015); Soos
_et al._ (2016), charge density waves Kumar  _et al._ (2011, 2009); Hirsch
(1984, 1985), spin density waves Dey  _et al._ (2016); White and Huse (1993);
So/rensen and Affleck (1994), vector chiral phase Chubukov (1991); Hikihara
_et al._ (2008); Sudan  _et al._ (2009); Parvej and Kumar (2016), the valence
bond solid Affleck  _et al._ (1987, 1988); Schollwöck  _et al._ (1996) and
topological phases Gu and Wen (2009); Pollmann  _et al._ (2012); Haldane
(1983, 1988), to name a few. The QPT is different from the classical or
thermal phase transition in which phase transition is driven by the
competition between energy and thermal entropy at finite temperature. Quantum
fluctuations are dominant in confined systems like one, quasi-one or two
dimensional systems and lead to many interesting quantum phases. In a quantum
system the EE which indirectly measures the correlation between different
parts of the system is expected to change in QPT. The extent of entanglement
depends on Hamiltonian parameters such as system geometry, competing
interaction parameters, magnetic and electric fields.

Effective low dimensional systems such as spin chains and spin ladders are of
enduring interest for the past many decades. Fundamentally different spectral
gaps in spin-1/2 and spin-1 chains predicted by Haldane, spin-Peierls’
distortion in spin-1/2 chains predicted from theory, and occurrence of
magnetization plataeus in spin ladders have all been observed experimentally
Kikuchi  _et al._ (2005, 2006); Gu and Su (2006); Hase  _et al._ (2006);
Maignan  _et al._ (2004); Hardy  _et al._ (2004, 2006); Ishiwata  _et al._
(2005); Wang  _et al._ (2011); Yao (2012); Lenertz  _et al._ (2011); He  _et
al._ (2009); Shiramura  _et al._ (1998). Thus, most theoretical predictions of
spin chains have been vindicated by experiments, mostly in molecular systems
and magnetically low-dimensional systems characterized by strong exchange in
one direction and much weaker exchange in orthogonal directions Dagotto
(2005). In fact, experimentally, there is evidence of systems Castilla  _et
al._ (1995) which belong to different regions of the quantum phase diagram in
the well studied $J_{1}$ – $J_{2}$ – $\delta$ model of spin chains, where
$J_{1}$ is the first neighbour exchange interaction, $J_{2}$ is the second
neighbour exchange interaction and $\delta$ is the dimerization in the nearest
neighbour exchange strength.

Belonging to the class of spin ladders is a class of skewed ladders where one
or more rung bonds in the unit cell is slanted. This ladder system can also be
mapped into zigzag ladders with periodically missing bonds as shown in Fig. 1.
The skewed ladder is classified as n/m ladder where n and m are the ring sizes
of adjacent rings in the ladder. These systems, in our opinion, are realizable
in concatenated transition metal complexes. This class of systems show a very
interesting quantum phase diagram as a function of the ratio of the rung
exchange ($J_{1}$) to the nearest neighbour leg exchange ($J_{2}$). Depending
upon this ratio n/m and ${J_{1}}/{J_{2}}$ and system size, the skewed ladders
show remarkable switching in the GS spin, exhibit bond order wave, spin
density wave and chiral vector phases Giri  _et al._ (2017). The Heisenberg
antiferromagnetic (HAF) spin-1/2 system on 5/7, 3/4, 3/5 and 5/5 skewed
ladders shows interesting magnetic and non-magnetic QPT in the parameter space
of the rung and the leg spin exchange constants Giri  _et al._ (2017); Dey
_et al._ (2020). Recent work on HAF spin-1 on 5/7 Das  _et al._ (2021a), 3/4
and 3/5 Das  _et al._ (2021b) skewed ladders has also shown existence of
similar magnetic and non-magnetic phase transitions. The quantum phases of HAF
spin-1/2 system on 3/4, 5/5 and 3/5 skewed ladders studied at finite external
axial magnetic field $B$ are found to show various magnetic plateau phases
which follow the Oshikawa, Yamanaka, and Affleck (OYA) Oshikawa  _et al._
(1997) rule, in the $B-J_{1}-J_{2}$ parameter space.

In general correlation functions, energy crossovers, GS symmetry changes and
local order parameters are studied to characterize the QPT. However, hidden
order like topological properties can not be detected by conventional
procedures. Therefore, to study the long range correlations in the system and
to characterize the phase boundaries sharply, we have studied the entanglement
properties of this system. Where entanglement does not show sharp features,
fidelity of the GS is a useful tool for determining the phase boundary. The
QPT of frustrated $J_{1}-J_{2}$ model has been characterised using EE in
earlier studies Li  _et al._ (2016). Simillarly the QPT of XXZ model has been
studied using the fidelity as a property obtained from exact diagonalization
and density matrix renormalization group methods Luo and Xu (2015). The
fidelity approach to characterize QPT is reviewed in GU (2010).

We find that fidelity is very sensitive to changes in the GS and shows sharp
changes even when the spin gap vs $J_{1}/J_{2}$ plots show only a change in
the slope or when there is a cross over in the two low-lying excited states.
When there is symmetry breaking, the fidelity in one of the subspaces shows a
sharp change while in another subspace there is no change. This allows us to
identify the symmetry before and after a quantum phase transition. In this
work, we follow the changes in EE and fidelity of different skewed spin-1/2
ladders. We mainly focus on exact diagonalization studies of these ladders. We
also employ symmetries such as the conservation of the z-component of total
spin, $S_{z}$, and the reflection symmetry at degeneracies to follow changes
in EE and fidelity. We find that in the 3/4 skewed ladder, EE shows sharp
changes at the critical value of $J_{1}$. Fidelity also shows a similar sharp
changes and we can identify the critical $J_{1}$ value very accurately. In the
skewed ladders 3/5 and 5/7, both EE and fidelity show sharp changes and the
critical $J_{1}$ values for QPT can be determined highly accurately. In what
follows, we discuss the results of our studies on the 3/4, 3/5 and 5/7 ladder
systems and show how both EE and fidelity can be used to accurately find
critical $J_{1}$ values at the transition points.

This paper is divided into five sections. In section II we provide a brief
introduction to entanglement entropy and fidelity. We then present the model
Hamiltonian and numerical methods in section III. The results are discussed in
section IV under three subsections. The summary of our results are presented
in section V.

##  II Entanglement entropy and Fidelity

Entanglement entropy (EE) has been studied extensively from a quantum
information perspective due to the promise that quantum computing holds. It
has also been shown in many studies that EE is a useful tool for studying
phase transition in quantum systems. Given two subsystems A and B and spanning
Hilbert spaces $H_{A}$ and $H_{B}$ with full system Hilbert space
$H=H_{A}\times H_{B}$, we can always find an orthonormal bases
$\\{\phi_{A}\\}$ and $\\{\phi_{B}\\}$ such that any state
$|{\psi}_{AB}\rangle$ in the Hilbert space belonging to $H$, can be expressed
as

$|\psi_{AB}\rangle=\sum_{i,j}\alpha_{ij}|\phi_{i}\rangle_{A}|\phi_{j}\rangle_{B}.$
(1)

The reduced density matrix of the subsystem B is given by

$\rho_{B}=Tr_{A}\left(|\psi_{AB}\rangle\langle\psi_{AB}|\right),$ (2)

and the $jj^{{}^{\prime}th}$ matrix element of the reduced density matrix
$\rho_{B}$ is,

${(\rho_{B})}_{jj^{{}^{\prime}}}=\sum_{i}\alpha_{ij}\alpha{{}^{*}}_{ij{{}^{{}^{\prime}}}}.$
(3)

If the eigenvalues of $\rho_{B}$ are $\lambda_{i}$ then we can define EE in
many different ways. If we restrict ourselves to extensivity of entropy, i.e.
$S_{AB}=S_{A}+S_{B}$, we can define entropy of a state in two ways. In
condensed matter physics, the widely used entropy is the von Neumann entropy
which is given by Petz (2001)

$S=-\sum_{i}\lambda_{i}\log_{2}\lambda_{i},$ (4)

where $\lambda_{i}^{\prime}$s are the eigenvalues of the reduced density
matrix. Von Neumann entropy is the generalization of the classical Shannon
entropy. Another entropy function which is Rényi entropy $S_{A,\alpha}^{R}$ is
defined as Rényi (1961)

$S_{A,\alpha}^{R}=\frac{1}{1-\alpha}\textrm{log}\left[\sum_{i}(\rho_{i}^{A})^{\alpha}\right]$
(5)

where $\alpha$ can take values between 0 and $\infty$. Rényi entropy reduces
to von Neumann entropy in the limit $\alpha\to 1$ and it is maximum when
$\alpha=0$ and minimum when $\alpha\to\infty$.

In this work, we employ von Neumann entropy to analyze the QPT as von Neumann
entropy is the most widely used entropy to study the QPT. We also restrict
ourselves to equal sizes of the subsystems A and B as von Neumann entropy is
maximum in this case. In frustrated one-dimensional spin chains, which are
known as the $J_{1}-J_{2}$ chains, the GS EE cannot identify the critical
$J_{2}$ value ($J_{2c}$) below which the system is gapless and above which the
system is gapped. However, the GS EE can identify the Majumdar-Ghosh point
which is not surprising since the GS corresponds to nearest neighbor singlets
Goli  _et al._ (2013). The $J_{2c}$ is correctly identified from EE when the
$J_{2}$ corresponding to crossing of the EE of the first excited triplet and
the second singlet is extrapolated from finite size calculations to the
thermodynamic limit Goli  _et al._ (2013). EE has also been studied in spin
ladders and it has been shown that the area law of entanglement is valid upto
seven legs Kallin  _et al._ (2009).

A more sensitive property to follow QPT is fidelity and fidelity
susceptibility. Fidelity ($F(\omega)$) at a parameter $\omega$ of the
Hamiltonian is given by the overlap of the desired state of the Hamiltonian at
$\omega$ with that in its neighborhood $\omega+\delta\omega$ i.e.

$F(\omega)=\langle\psi(\omega)|\psi(\omega+\delta\omega)\rangle.$ (6)

We can also define Fidelity susceptibility $\chi(\omega)$ as

$\chi(\omega)=\frac{2(1-F(\omega))}{(\delta\omega)^{2}}$ (7)

which follows from a Taylor series expansion of $F(\omega)$. It has been shown
that fidelity of the GS in the $J_{1}-J_{2}$ model can give accurate value of
$J_{2c}$ Thesberg and Sørensen (2011). In the $J_{1}-J_{2}$ model, the
critical value of $J_{2}$ has also been determined by following the fidelity
of the first excited state of short Heisenberg chains to determine $J_{2}(N)$
and extrapolating it to the thermodynamic limit $N\to\infty$ Chen  _et al._
(2007). In the $J_{1}-J_{2}-\delta$ model, using EE, it is not possible to
obtain the quantum phase diagram accurately; the EE contours show phase
changes although phase boundries remain fuzzy. We have obtained the EE and
fidelity of the 3/4, 3/5, and 5/7 skewed ladders; these systems show quantum
phase transition for fixed next nearest neighbor exchange $J_{2}$ and varying
nearest neighbor $J_{1}$ exchange.

##  III Model Hamiltonian and Numerical method

In our study, we consider a spin-1/2 model on 3/4, 3/5 and 5/7 skewed ladders
as shown in Fig. 1(a), (b) and (c). The rung bond interaction is $J_{1}$ and
the leg bond interaction is $J_{2}$. All the exchange interactions in the
Hamiltonians we consider are antiferromagnetic in nature. The sites are
numbered in such a way that the even numbered sites are on the top of the leg
and odd numbers are on the bottom of the leg. The leg bond interaction $J_{2}$
is set to 1 and it defines the energy scale. The model Hamiltonian of 3/4
skewed ladder is written as

$\displaystyle H_{3/4}$ $\displaystyle=$ $\displaystyle
J_{1}\sum_{i}\left[\left(\vec{S}_{i,1}+\vec{S}_{i,3}\right)\cdot\vec{S}_{i,2}+\left(\vec{S}_{i,4}+\vec{S}_{i,6}\right)\cdot\vec{S}_{i,5}\right]$
(8)
$\displaystyle+J_{2}\sum_{i}\bigg{(}\vec{S}_{i,5}\cdot\vec{S}_{i+1,1}+\vec{S}_{i,6}\cdot\vec{S}_{i+1,2}$
$\displaystyle\qquad\qquad+\sum_{k=1}^{4}\vec{S}_{i,k}\cdot\vec{S}_{i,k+2}\bigg{)}$

where $i$ labels the unit cell and $k$ the spins within the unit cell. The
first term denotes the rung bond interaction and the second term denotes the
leg bond interaction. Similarly, the model Hamiltonian for the 3/5 and 5/7
systems can be written as

Figure 1: Schematic diagram of (a) 3/4 skewed ladder, (b) 3/5 skewed ladder,
and (c) 5/7 skewed ladder. The nearest neighbor or rung interaction is $J_{1}$
and the next-nearest neighbor (according to our numbering scheme) interaction
is $J_{2}$. The sites on the top leg are even numbered and on the bottom leg
are odd numbered. The unit cell for each ladder system is indicated by a
rectangle with dashed red lines. The mirror plane for each ladder system is
represented by a perpendicular blue line.

$\displaystyle H_{3/5}$ $\displaystyle=$ $\displaystyle
J_{1}\sum_{i}\left(\vec{S}_{i,1}+\vec{S}_{i,3}\right)\cdot\vec{S}_{i,2}+J_{2}\sum_{i}\bigg{(}\vec{S}_{i,3}\cdot\vec{S}_{i+1,1}$
(9)
$\displaystyle+\vec{S}_{i,4}\cdot\vec{S}_{i+1,2}+\sum_{k=1}^{2}\vec{S}_{i,k}\cdot\vec{S}_{i,k+2}\bigg{)},$

and

$\displaystyle H_{5/7}$ $\displaystyle=$ $\displaystyle
J_{1}\sum_{i}\left(\vec{S}_{i,1}\cdot\vec{S}_{i,2}+\vec{S}_{i,4}\cdot\vec{S}_{i,5}\right)$
(10)
$\displaystyle+J_{2}\sum_{i}\bigg{(}\vec{S}_{i,7}\cdot\vec{S}_{i+1,1}+\vec{S}_{i,8}\cdot\vec{S}_{i+1,2}$
$\displaystyle\qquad\qquad+\sum_{k=1}^{6}\vec{S}_{i,k}\cdot\vec{S}_{i,k+2}\bigg{)}.$

We have carried out the computations using exact diagonalization technique on
a system with 24 spin-1/2 sites, in all cases. We have employed cyclic
boundary conditions and the EE computed is the von Neumann entropy. The
environment and the system, both have equal number of sites, i.e. 12 sites.
The reduced density matrix of the left half of the system is obtained by
tracing over the right half of the system,

$\rho_{L{L^{{}^{\prime}}}}=\sum_{R}{C_{LR}}C_{L^{{}^{\prime}}R}$ (11)

where the GS wavefunction $\psi_{g}$ is expressed in the direct product basis
of the Fock space of the system ($L$) and the surroundings ($R$). The
dimensionality of the Fock space of the system block, and hence the order of
the density matrix is $2^{12}\times 2^{12}$ or $4096\times 4096$. After
obtaining all the eigenvalues $\\{\lambda_{i}\\}$ of the reduced density
matrix, the von Neumann entropy is obtained from Eqn. 4. We have computed the
fidelity $F(J_{1})$ of the state using the normalized ground states and Eqn.
6. We have varied $J_{1}$ in very small steps of $\delta(=0.001)$ near the
quantum phase transitions obtained from the spin gap data. The spin gap
$\Gamma_{l}$ is defined as the energy gap between the lowest eigenstates
$S_{z}=l$ and $S_{z}=0$ manifolds, where $S_{z}$ is the z-component of the
total spin.

$\Gamma_{l}=E_{0}(S_{z}=l)-E_{0}(S_{z}=0),$ (12)

where $l$ is an integer. The GS spin $(S_{G})$ $=l$, for $l$ that satisfies
$\Gamma_{l}=0$ and $\Gamma_{l+1}>0$. The excitation gap $\Gamma_{\sigma}$ at a
fixed $S_{G}$ between the lowest energy states in different reflection
symmetry subspaces $\sigma(-)$ and $\sigma(+)$ is defined as

$\Gamma_{\sigma}=E_{0}(S_{G},\sigma(-))-E_{0}(S_{G},\sigma(+)).$ (13)

The GS is in odd subspace ($\sigma(-)$) when $\Gamma_{\sigma}<0$, even
subspace ($\sigma(+)$) when $\Gamma_{\sigma}>0$ and doubly degenerate when
$\Gamma_{\sigma}=0$.

##  IV Results and discussions

###  IV.1 3/4 skewed ladder

This is the simplest of the skewed ladders to show a transition in the spin of
the GS. However, there is a transition in the GS near $J_{1}=1.3$ which does
not involve change in the GS spin.

Figure 2:  Spin gaps $\Gamma_{l}$ for a 3/4 skewed ladder of 24 spins with PBC
shown as a function of $J_{1}$. The GS spin $S_{G}=0$ for $J_{1}<1.581$ and
$S_{G}=4$ for $J_{1}\geq 1.581$.

We notice in the $\Gamma_{l}$ vs $J_{1}$ plot (Fig. 2) that there is a kink
near $J_{1}=1.3$, but it is not possible to pin point the $J_{1}$ value at
which this kink appears from the $\Gamma_{l}$ plot. However, the EE shows a
sharp change at $J_{1}=1.340$ and it is also seen more clearly in the
fidelity, where it vanishes at $J_{1}=1.340$ but has a value of 1.0 for other
$J_{1}$ values (Fig. 3). Closer examination of the full eigenvalue spectrum of
the Hamiltonian shows that at this value of $J_{1}$, the singlet GS is doubly
degenerate and there is a cross over from one singlet GS to another singlet
GS. The transition to the highest spin states occurs sharply, when the GS spin
changes to $S_{G}=4$.

Figure 3:  The behaviour of EE (S) and fidelity (F) in the range $1.337\leq
J_{1}\leq 1.343$ of a 3/4 skewed ladder with N = 24 spins are shown in
($a_{1}$) and ($a_{2}$) respectively. The EE changes abruptly and a sharp drop
in fidelity occurs at $J_{1}=1.340$. At other values of $J_{1}$ fidelity is
constant (=1). Simillary the behaviour of EE and fidelity in the range
$1.578\leq J_{1}\leq 1.584$ are shown in ($b_{1}$) and ($b_{2}$). At
$J_{1}=1.581$, the entropy changes and fidelity shows a sharp drop.

The transition occurs at $J_{1}=1.581$ and there is a sharp drop in EE at this
point. Fidelity also goes to zero at $J_{1}=1.581$ but assumes a value of $1$,
at the neighboring points. For larger system sizes, the level crossing method
is used to calculate the transition points, and it has been reported that the
transition points do not change significantly with system size (ref. Giri  _et
al._ (2017)). In this system, we have been able to determine the parameter at
which the QPT occurs very accurately and we also find that there is a
transition at $J_{1}=1.340$ which is not evident from the plots of the spin
gaps.

###  IV.2 3/5 skewed ladder

The GS of 3/5 skewed ladder is nonmagnetic for $J_{1}<2.026$ and the spin of
the GS progressively increases from $S_{G}=0$ to $S_{G}=3$ form $J_{1}=2.026$
to $2.297$ (Fig. 4(a)).

Figure 4: (a) The spin gaps $\Gamma_{1}$, $\Gamma_{2}$ and $\Gamma_{3}$ are
shown as a function of $J_{1}$ for N = 24 sites in a 3/5 skewed ladder with
PBC. For $J_{1}$ $<$ 2.026 and $J_{1}$ $>$ 6.859, $\Gamma_{1}$ becomes 0 and
the system shows nonmagnetic behaviour. In the region $2.026\leq J_{1}\leq
2.297$ the GS spin gradually changes from 0 to 3. (b) shows the EE and (c)
shows the fidelity for the unsymmetrized GS as a function of $J_{1}$. The EE
exhibits a discontinuous change, while fidelity shows a sharp drop at the
transition points. The thick line near the vicinity of $J_{1}=2$ indicates
many transitions both in EE and fidelity.

The spin of the GS remains 3 as $J_{1}$ is increased further until
$J_{1}=6.859$ and for $J_{1}>6.859$ the state reenters the nonmagnetic phase.
EE and fidelity are shown in (Fig. 4(b),(c)) over the entire range of our
study, namely $1\leq J_{1}\leq 8$. In this ladder system, both EE and fidelity
give the critical value of $J_{1}$ to be 1.217 (Fig. 5) and at this value of
$J_{1}$ both the lowest energy state and first excited state are same upto
third decimal space in $S_{z}=0$ sector. Here again it is likely due to broken
spatial symmetry in the GS. At $J_{1}=6.859$, we find a sharp transition in
both EE and fidelity to the reentrant nonmagnetic phase.

Figure 5:  The behaviour of EE (S) and fidelity (F) in the range of $J_{1}$
($1.214\leq J_{1}\leq 1.22$) of a 3/5 skewed ladder with N = 24 spins are
shown in ($a_{1}$) and ($a_{2}$) respectively. The change in entropy and sharp
drop in fidelity occurs at $J_{1}=1.217$. For other values of $J_{1}$ the
fidelity is constant (=1). Simillary behaviour of the entropy and fidelity in
the range of $J_{1}$ ($6.856\leq J_{1}\leq 6.862$) are shown in ($b_{1}$) and
($b_{2}$) respectively. At $J_{1}=6.859$, the entropy changes and fidelity
shows a sharp drop.

For $J_{1}$ in the neighborhood of 2, fidelity and EE seem to indicate many
transitions. The system has a reflection symmetry, as shown in Fig. 1(b). The
lowest energy states for $2.026\leq J_{1}\leq 2.046$ and $2.290\leq J_{1}\leq
2.296$, in the $\sigma(+)$ and $\sigma(-)$ manifold are degenerate. A linear
combination of the degenerate states can lead to a broken symmetry state (see
ref. Giri  _et al._ (2017) for details). The intervals of $J_{1}$ for which
the GS is doubly degenerate are shown in third column of Table 1.

Table 1:  The interval of $J_{1}$ for the GS in different subspaces of a 3/5
skewed ladder with $N=24$ spins. The first column represents the range of
$J_{1}$ for which the GS is in $\sigma(+)$ subspace, the second column
represents the GS in $\sigma(-)$ subspace. The third column gives the range of
$J_{1}$ for which the difference between the lowest energies of two different
subspaces is zero.

$\sigma(+)$ |  $\sigma(-)$ |  $|\Gamma_{\sigma}|=0$
---|---|---
$1<J_{1}\leq 1.217$ | $1.217\leq J_{1}<2.025$ | $2.026\leq J_{1}\leq 2.046$
and |  and |  and
$2.046<J_{1}<2.290$ |  $J_{1}>2.296$ | $2.290\leq J_{1}\leq 2.296$

The EE of the eigenstate with symmetry shows a kink at $J_{1}=2.025$ below
which the GS is in the $\sigma(-)$ space ($\sigma$ is the reflection symmetry)
and the EE in the $\sigma(+)$ space shows a kink at $2.046$ (Fig. 6(a)). The
fidelity of the GS in $\sigma(-)$ space show a drop at $2.025$ and that in the
$\sigma(+)$ space shows a drop at 2.046 (Fig. 6(b)).

Figure 6:  The EE (S) and fidelity (F) of lowest energy states of a 3/5 skewed
ladder for 24 sites with periodic boundary condition in different reflection
symmetry subspaces (($\sigma(+)$) and ($\sigma(-)$). Both the EE and fidelity
for the unsymmetrized GS is represented by the black curve whereas for the
symmetrised GS these are represented by the red and blue curves.

In the region between these values, the GS is doubly degenerate and any linear
combination of these states is also an eigenstate and hence the EE and
fidelity, calculated from unsymmetrized GS show wild fluctuations. For
$2.290\leq J_{1}\leq 2.296$ the GS is degenerate. So the unsymmetrized
calculation for both EE and fidelity shows wild fluctuations in this range of
$J_{1}$. In symmetrized calculation, both the EE and fidelity change at
$J_{1}=2.289$ and 2.296 (Fig. 6(c) and 6(d)). The $|\Gamma_{\sigma}|$ value
vanishes in the interval $2.290\leq J_{1}\leq 2.296$. The GS below 2.290 has
$\sigma(+)$ symmetry whereas above 2.296 the GS has $\sigma(-)$ symmetry.
Therefore fidelity for the $\sigma(+)$ subspace shows a sharp dip at
$J_{1}=2.289$ and that of state $\sigma(-)$ shows a dip at $J_{1}=2.296$.

###  IV.3 5/7 skewed ladder

Among the skewed ladders we have studied, the 5/7 skewed ladder has by far the
most interesting quantum phase transitions. For fixed $J_{2}=J_{1}=1$, the GS
spin of the system, $S_{G}$, increases with system size, a feature we have not
observed in the other systems studied here. In this work, we focus on the EE
and fidelity for a 24 spin 5/7 ladder with periodic boundary condition as a
function of $J_{1}$. From the plot of spin gaps vs. $J_{1}$(Fig. 7(a)) and

Figure 7: (a) The spin gaps $\Gamma_{1}$, $\Gamma_{2}$ and $\Gamma_{3}$ are
shown as a function of $J_{1}$ for 24 site spins in a 5/7 skewed ladder with
PBC. For $J_{1}<1.427$, $\Gamma_{1}$ becomes nonzero and the system shows
nonmagnetic behaviour. The system enters a reentrant nonmagnetic phase for
$1.734<J_{1}<1.872$ where $\Gamma_{1}$ is nonzero. The EE for the
unsymmetrized GS is shown in (b), and the fidelity is shown in (c) as a
function of $J_{1}$ for the unsymmetrized GS. The EE exhibits a discontinuous
change, while fidelity shows a sharp drop at the transition points. Wild
fluctuations (many transitions) in both the EE and fidelity around $J_{1}=1.8$
and 2.3 are due to the degeneracy in the GS energy.

earlier studies it is established that the GS switches from a singlet to a
triplet at $J_{1}=1.427$. For $1.734<J_{1}<1.872$ the GS of the system again
becomes a singlet and for $J_{1}>1.872$ the GS changes spin, finally attaining
a $S_{G}=3$ for $J_{1}\geq 2.355$. In (Fig. 7 (b) and (c)) we show how the EE
and fidelity varies as a function of $J_{1}$. At the first transition, both
these quantities show a sharp change. In the region of the reentrant phase,
the variations are violent. In the region $1.872<J_{1}<2.176$, EE shows smooth
variation and the fidelity also does not change. In this region the GS spin,
$S_{G}=1$. Although the $S_{G}$ value continues to be 1 in the region
$2.176<J_{1}<2.355$, the EE as well as fidelity show very sharp changes.
Finally beyond $J_{1}=2.355$, EE is almost constant and fidelity stays at 1.

The first transition from a singlet to a triplet is simple and the critical
$J_{1}$ value can be pinned down both from EE and fidelity. This occurs at
${J_{1}}=1.427$. To understand the behavior at other values of $J_{1}$, we
recognize that the system has a reflection symmetry perpendicular to the legs.
So, the states of the system can be classified as $\sigma{(+)}$ or
$\sigma{(-)}$ depending on the symmetry of the space in which the GS is found
Giri  _et al._ (2017). The range of $J_{1}$ for the lowest energy state in
different subspaces is shown in Table 2.

Table 2:  The range of $J_{1}$ over which the GS have different symmetries in
a 5/7 skewed ladder with $N=24$ spins. The first column represents the range
of $J_{1}$ for which the GS symmetry is $\sigma(+)$ and similarly the second
column represents the GS is $\sigma(-)$ subspace. The third column gives the
intervals of $J_{1}$ for which the GS is doubly degenerate.

$\sigma(+)$ |  $\sigma(-)$ |  $|\Gamma_{\sigma}|=0$
---|---|---
$1<J_{1}\leq 1.427$ | $1.427\leq J_{1}\leq 1.753$ | $1.753<J_{1}\leq 1.871$
and |  and |  and
$1.871<J_{1}\leq 2.176$ |  $J_{1}>2.355$ | $2.176<J_{1}\leq 2.355$

We see that the change in the spin of the GS from $S_{G}=0$ to $S_{G}=1$ at
$J_{1}=1.427$ is also accompanied by vanishing of $|\Gamma_{\sigma}|$. EE
shows a jump and fidelity shows a sharp dip at $J_{1}=1.615$ and 1.725. The
change in EE and fidelity at 1.615 is associated with the maximum in the
$|\Gamma_{\sigma}|$ gap. Beyond $J_{1}=1.615$, the $|\Gamma_{\sigma}|$ gap
begins to close and at $J_{1}=1.754$ the gap vanishes. In region between 1.753
and 1.871, $|\Gamma_{\sigma}|$ vanishes. In the unsymmetrized calculations
both fidelity and entropy show sharp fluctuations, which are suppressed when
the degenerate ground states are obtained in the $\sigma{(+)}$ and
$\sigma{(-)}$ subspaces (Fig. 8). In the $J_{1}$ region between 2.176 to 2.355
the GS is a spin triplet and $|\Gamma_{\sigma}|$ also vanishes. Fidelity of
the $\sigma{(+)}$ and $\sigma{(-)}$ show a sharp dip at $J_{1}=2.176$ and
$2.355$ respectively; the GS below 2.176 has $\sigma{(+)}$ symmetry and
becomes degenerate with the lowest state in $\sigma{(-)}$ symmetry in the
region of $J_{1}$ ($2.176<J_{1}\leq 2.355$) and for $J_{1}$ value above 2.355
the GS switches back to $\sigma{(-)}$ subspace. There are no sharp changes in
either entropy or fidelity beyond $J_{1}=2.355$.

Figure 8:  For a 5/7 skewed ladder with 24 sites in PBC, (a) EE and (b)
fidelity are plotted as a function of $J_{1}$ for both the reflection symmetry
subspaces ($\sigma(+)$ and $\sigma(-)$). Black curve corresponds to the
$\sigma(+)$ subspace, whereas the red curve represents the $\sigma(-)$
subspace.

##  V Summary

We have studied EE and fidelity to characterize the quantum phase transitions
in skewed ladder systems. While correlation function also give a measure of
the extent of interactions between sites, EE embodies long range correlation
between different parts of the system. In the regions of $J_{1}$ where the
system has a nondegenerate GS without change in spin or spatial symmetry the
entropy change is gradual and fidelity remains at one. At the transition
points, the entropy shows a discontinuous change and fidelity shows a sharp
dip. The transition points are accurately determined from these characteristic
changes. In regions where the GS is degenerate, the unsymmetrized calculations
show sharp fluctuations in EE and fidelity. However, in the symmetrized
calculations they show abrupt changes only at the transition. We also find
that for regions where $|\Gamma_{\sigma}|$ vanishes, the fidelity does not
change in either subspace. But at the beginning and at the end of the
transition, fidelity shows sharp dip in the lowest eigenstate of one of the
symmetry subspaces.

## acknowledgments

S. Ramasesha acknowledges the Indian National Science Academy and DST-SERB for
supporting this work. Manoranjan Kumar acknowledges the SERB for financial
support through Project File No. CRG/2020/000754.

## Author contribution statement

The project was conceived by S. Ramasesha and Manoranjan Kumar. Sambunath Das
and Dayasindhu Dey have performed the numerical calculations. All the authors
were involved in designing the calculations and interpretation of results. The
work was jointly written up by all the authors.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Majumdar and Ghosh (1969) C. K. Majumdar and D. K. Ghosh, Journal of Mathematical Physics 10, 1399 (1969).
  * Chubukov (1991) A. V. Chubukov, Phys. Rev. B 44, 4693 (1991).
  * Chitra  _et al._ (1995) R. Chitra, S. Pati, H. R. Krishnamurthy, D. Sen, and S. Ramasesha, Phys. Rev. B 52, 6581 (1995).
  * White and Affleck (1996) S. R. White and I. Affleck, Phys. Rev. B 54, 9862 (1996).
  * Itoi and Qin (2001) C. Itoi and S. Qin, Phys. Rev. B 63, 224423 (2001).
  * Mahdavifar (2008) S. Mahdavifar, Journal of Physics: Condensed Matter 20, 335230 (2008).
  * Sirker (2010) J. Sirker, Phys. Rev. B 81, 014419 (2010).
  * Kumar  _et al._ (2015) M. Kumar, A. Parvej, and Z. G. Soos, Journal of Physics: Condensed Matter 27, 316001 (2015).
  * Soos  _et al._ (2016) Z. G. Soos, A. Parvej, and M. Kumar, Journal of Physics: Condensed Matter 28, 175603 (2016).
  * Kumar  _et al._ (2010) M. Kumar, S. Ramasesha, and Z. G. Soos, Phys. Rev. B 81, 054413 (2010).
  * Hamada  _et al._ (1988) T. Hamada, J.-i. Kane, S.-i. Nakagawa, and Y. Natsume, Journal of the Physical Society of Japan 57, 1891 (1988).
  * Kumar  _et al._ (2011) M. Kumar, B. J. Topham, R. Yu, Q. B. D. Ha, and Z. G. Soos, The Journal of Chemical Physics 134, 234304 (2011).
  * Kumar  _et al._ (2009) M. Kumar, S. Ramasesha, and Z. G. Soos, Phys. Rev. B 79, 035102 (2009).
  * Hirsch (1984) J. E. Hirsch, Phys. Rev. Lett. 53, 2327 (1984).
  * Hirsch (1985) J. E. Hirsch, Phys. Rev. B 31, 6022 (1985).
  * Dey  _et al._ (2016) D. Dey, M. Kumar, and Z. G. Soos, Phys. Rev. B 94, 144417 (2016).
  * White and Huse (1993) S. R. White and D. A. Huse, Phys. Rev. B 48, 3844 (1993).
  * So/rensen and Affleck (1994) E. S. So/rensen and I. Affleck, Phys. Rev. B 49, 15771 (1994).
  * Hikihara  _et al._ (2008) T. Hikihara, L. Kecke, T. Momoi, and A. Furusaki, Phys. Rev. B 78, 144404 (2008).
  * Sudan  _et al._ (2009) J. Sudan, A. Lüscher, and A. M. Läuchli, Phys. Rev. B 80, 140402 (2009).
  * Parvej and Kumar (2016) A. Parvej and M. Kumar, Journal of Magnetism and Magnetic Materials 401, 96 (2016).
  * Affleck  _et al._ (1987) I. Affleck, T. Kennedy, E. H. Lieb, and H. Tasaki, Phys. Rev. Lett. 59, 799 (1987).
  * Affleck  _et al._ (1988) I. Affleck, T. Kennedy, E. H. Lieb, and H. Tasaki, Commun.Math. Phys. 115, 477–528 (1988).
  * Schollwöck  _et al._ (1996) U. Schollwöck, O. Golinelli, and T. Jolicœur, Phys. Rev. B 54, 4038 (1996).
  * Gu and Wen (2009) Z.-C. Gu and X.-G. Wen, Phys. Rev. B 80, 155131 (2009).
  * Pollmann  _et al._ (2012) F. Pollmann, E. Berg, A. M. Turner, and M. Oshikawa, Phys. Rev. B 85, 075125 (2012).
  * Haldane (1983) F. Haldane, Physics Letters A 93, 464 (1983).
  * Haldane (1988) F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).
  * Kikuchi  _et al._ (2005) H. Kikuchi, Y. Fujii, M. Chiba, S. Mitsudo, T. Idehara, T. Tonegawa, K. Okamoto, T. Sakai, T. Kuwai, and H. Ohta, Phys. Rev. Lett. 94, 227201 (2005).
  * Kikuchi  _et al._ (2006) H. Kikuchi, Y. Fujii, M. Chiba, S. Mitsudo, T. Idehara, T. Tonegawa, K. Okamoto, T. Sakai, T. Kuwai, and H. Ohta, Phys. Rev. Lett. 97, 089702 (2006).
  * Gu and Su (2006) B. Gu and G. Su, Phys. Rev. Lett. 97, 089701 (2006).
  * Hase  _et al._ (2006) M. Hase, M. Kohno, H. Kitazawa, N. Tsujii, O. Suzuki, K. Ozawa, G. Kido, M. Imai, and X. Hu, Phys. Rev. B 73, 104419 (2006).
  * Maignan  _et al._ (2004) A. Maignan, V. Hardy, S. Hébert, M. Drillon, M. R. Lees, O. Petrenko, D. M. K. Paul, and D. Khomskii, J. Mater. Chem. 14, 1231 (2004).
  * Hardy  _et al._ (2004) V. Hardy, D. Flahaut, M. R. Lees, and O. A. Petrenko, Phys. Rev. B 70, 214439 (2004).
  * Hardy  _et al._ (2006) V. Hardy, C. Martin, G. Martinet, and G. André, Phys. Rev. B 74, 064413 (2006).
  * Ishiwata  _et al._ (2005) S. Ishiwata, D. Wang, T. Saito, and M. Takano, Chemistry of Materials 17, 2789 (2005).
  * Wang  _et al._ (2011) X. X. Wang, J. J. Li, Y. G. Shi, Y. Tsujimoto, Y. F. Guo, S. B. Zhang, Y. Matsushita, M. Tanaka, Y. Katsuya, K. Kobayashi, K. Yamaura, and E. Takayama-Muromachi, Phys. Rev. B 83, 100410 (2011).
  * Yao (2012) X. Yao, The Journal of Physical Chemistry A 116, 2278 (2012), pMID: 22364513.
  * Lenertz  _et al._ (2011) M. Lenertz, J. Alaria, D. Stoeffler, S. Colis, and A. Dinia, The Journal of Physical Chemistry C 115, 17190 (2011).
  * He  _et al._ (2009) Z. He, J.-I. Yamaura, Y. Ueda, and W. Cheng, Journal of the American Chemical Society 131, 7554 (2009), pMID: 19489641\.
  * Shiramura  _et al._ (1998) W. Shiramura, K.-i. Takatsu, B. Kurniawan, H. Tanaka, H. Uekusa, Y. Ohashi, K. Takizawa, H. Mitamura, and T. Goto, Journal of the Physical Society of Japan 67, 1548 (1998).
  * Dagotto (2005) E. Dagotto, Science 309, 257 (2005).
  * Castilla  _et al._ (1995) G. Castilla, S. Chakravarty, and V. J. Emery, Phys. Rev. Lett. 75, 1823 (1995).
  * Giri  _et al._ (2017) G. Giri, D. Dey, M. Kumar, S. Ramasesha, and Z. G. Soos, Phys. Rev. B 95, 224408 (2017).
  * Dey  _et al._ (2020) D. Dey, S. Das, M. Kumar, and S. Ramasesha, Phys. Rev. B 101, 195110 (2020).
  * Das  _et al._ (2021a) S. Das, D. Dey, M. Kumar, and S. Ramasesha, Phys. Rev. B 104, 125138 (2021a).
  * Das  _et al._ (2021b) S. Das, D. Dey, S. Ramasesha, and M. Kumar, Journal of Applied Physics 129, 223902 (2021b).
  * Oshikawa  _et al._ (1997) M. Oshikawa, M. Yamanaka, and I. Affleck, Phys. Rev. Lett. 78, 1984 (1997).
  * Li  _et al._ (2016) Y.-C. Li, Y.-H. Zhu, and Z.-G. Yuan, Physics Letters A 380, 1066 (2016).
  * Luo and Xu (2015) D.-W. Luo and J.-B. Xu, Annals of Physics 354, 298 (2015).
  * GU (2010) S.-J. GU, International Journal of Modern Physics B 24, 4371 (2010).
  * Petz (2001) D. Petz, “Entropy, von neumann and the von neumann entropy,” in  _John von Neumann and the Foundations of Quantum Physics_, edited by M. Rédei and M. Stöltzner (Springer Netherlands, Dordrecht, 2001) pp. 83–96.
  * Rényi (1961) A. Rényi, In: Proceedings of the 4th Berkeley Symposium on Mathematics, Statistics and Probability 1, 547 (1961).
  * Goli  _et al._ (2013) V. M. L. D. P. Goli, S. Sahoo, S. Ramasesha, and D. Sen, Journal of Physics: Condensed Matter 25, 125603 (2013).
  * Kallin  _et al._ (2009) A. B. Kallin, I. González, M. B. Hastings, and R. G. Melko, Phys. Rev. Lett. 103, 117203 (2009).
  * Thesberg and Sørensen (2011) M. Thesberg and E. S. Sørensen, Phys. Rev. B 84, 224435 (2011).
  * Chen  _et al._ (2007) S. Chen, L. Wang, S.-J. Gu, and Y. Wang, Phys. Rev. E 76, 061108 (2007).