Transverse Single Spin Asymmetries of Heavy Flavor Electrons and Charged Pions
in 200 GeV $p+p^{\uparrow}$ Collisions at Midrapidity

Dillon S. Fitzgerald1${\dagger}\star$

1 University of Michigan, Ann Arbor, USA

${\dagger}$ On behalf of the PHENIX collaboration

$\star$<EMAIL_ADDRESS>
Proceedings for the XXVIII International Workshop
on Deep-Inelastic Scattering and Related Subjects,
Stony Brook University, New York, USA, 12-16 April 2021
10.21468/SciPostPhysProc.?

## Abstract

Transverse single spin asymmetries of particles produced in $p+p^{\uparrow}$
collisions provide insight on the partonic spin and momentum structure of
hadrons; heavy flavor electrons provide access to initial state spin-momentum
correlations of gluons in the proton, while charged pions provide access to
initial and final state transverse spin effects. Charged particles are
measured at midrapidity at PHENIX using the silicon vertex detector and
central arm spectrometer, made of an electromagnetic calorimeter, a ring-
imaging Cherenkov detector, and drift and pad chambers. Recent results for
both electron and charged pion measurements from the 2015 running period will
be presented.

###### Contents

  1. 1 Introduction
  2. 2 Experimental Setup
    1. 2.1 Charged Particle Identification at PHENIX
    2. 2.2 Analysis Procedure
  3. 3 Heavy Flavor Electron TSSA
  4. 4 Charged Pion TSSA
  5. 5 Conclusion

##  1 Introduction

Our understanding of proton structure has evolved greatly over the past few
decades – from a naive picture where the valence quarks carry all of the
proton’s spin, to a much more complex picture involving spin and orbital
angular momentum of quarks and gluons. The complex structure of hadrons gives
rise to emergent properties of bound states such as spin-momentum and spin-
spin correlations within hadronic systems, similar to spin-orbit and spin-spin
couplings in atomic systems. These correlations are manifest in experimentally
observable ways, such as azimuthal modulations of particle yields in
collisions involving polarized hadrons. One such observable, the transverse
single-spin asymmetry (TSSA), involves collisions with an unpolarized beam
particle and a transversely polarized hadron in the initial state (e.g.
$p+p^{\uparrow}$ at the Relativistic Heavy Ion Collider (RHIC)), and is
defined as follows,

$A_{N}(\phi)=\frac{\sigma^{\uparrow}(\phi)-\sigma^{\downarrow}(\phi)}{\sigma^{\uparrow}(\phi)+\sigma^{\downarrow}(\phi)}=A_{N}\cos\phi.$
(1)

$\sigma^{\uparrow,\downarrow}(\phi)$ are cross sections measured as a function
of $\phi$ for collision systems with $p^{\uparrow}$ and $p^{\downarrow}$
initial state protons respectively [1]. Perturbative Quatnum Chromodynamics
(pQCD) calculations predict negligible contributions to TSSAs from partonic
hard scattering effects (<1%) [2]. However, large TSSAs have been measured as
a function of $x_{F}=2p_{z}/\sqrt{s}$ and $p_{T}$ of produced particles, even
persisting at $\sqrt{s}=200$ GeV and $p_{T}$ $\approx 7$ GeV/c, well into the
regime where pQCD applies[3, 4]. This implies the existence of nonperturbative
spin-momentum and spin-spin correlations within hadrons, as contributions to
large TSSAs must be coming from the nonperturbative elements of factorized
cross sections.

Two theoretical frameworks have emerged that have been successful in
describing the large observed TSSAs, one with transverse momentum dependent
(TMD) nonperturbative functions that depend explicitly on the transverse
momentum of partons within the hadron, and one with higher twist effects,
where higher order power suppressed terms in the collinear factorization
expansion must be considered. Both have the notable feature that one must look
beyond standard twist-2 collinear factorization – TMD functions require access
to both a hard and soft scale with sufficient scale separation (e.g.
$Q>>k_{T}$), while the higher twist framework preserves the collinear
factorization scheme, allowing for access from observables sensitive to only
the hard scale, which is the case for the measurements presented in this note.
Looking to higher twist leads to terms in the polarized cross section
including the convolution of correlation functions dependent on momentum
fractions of 2 partons in a given hadron with standard parton distribution
functions (PDFs) and fragmentation functions (FFs). This leads to the
following proportionality for TSSAs in the collinear twist-3 framework [5],

$\displaystyle
A_{N}\propto\sum_{a,b,c}\phi^{(3)}_{a/A}(x_{1},x_{2},\vec{s_{\perp}})\otimes\phi_{b/B}(x^{\prime})\otimes\hat{\sigma}\otimes
D_{c\rightarrow C}(z)$ (2) $\displaystyle+\sum_{a,b,c}\delta
q_{a/A}(x,\vec{s_{\perp}})\otimes\phi^{(3)}_{b/B}(x^{\prime}_{1},x^{\prime}_{2})\otimes\hat{\sigma^{\prime}}\otimes
D_{c\rightarrow C}(z)$ $\displaystyle+\sum_{a,b,c}\delta
q_{a/A}(x,\vec{s_{\perp}})\otimes\phi_{b/B}(x^{\prime})\otimes\hat{\sigma^{\prime\prime}}\otimes
D^{(3)}_{c\rightarrow C}(z_{1},z_{2}).$

Here $\delta q_{a/A}$ is the transversity distribution of parton $\it{a}$ in
polarized hadron $\it{A}$, and any term with a superscript $(3)$ is a twist-3
correlator, also known as a parton correlation function (PCF). Unification of
these two theoretical frameworks has been demonstrated, with twist 3
correlators related to moments in transverse momentum of the corresponding
parton TMD [6]. Hence one can learn something interesting about different
spin-momentum and spin-spin correlations from each term in Equation 2. The
first term is related to the Sivers TMD, while the second is related to
transversity $\otimes$ Boer-Mulders TMD and the third to transversity
$\otimes$ Collins TMD FF [7]. By measuring $A_{N}$ for particular final state
particles, different elements of Equation 2 can be isolated and constrained.
Heavy flavor electron production is dominated by gluon-gluon fusion at
$\sqrt{s}=200$ GeV at midrapidity, and the transversity distribution for
gluons $\delta_{g/X}=0$, providing access to the twist 3 correlator in the
polarized proton $\phi^{(3)}_{g/A}(x_{1},x_{2},\vec{s_{\perp}})$ while the two
terms that include convolution with transversity distributions are zero [5].
On the other hand, charged pion production is dominated by quark-gluon
scattering in the same kinematic region, providing broader access to the
different terms in Equation 2 and sensitivity to quark flavor.

##  2 Experimental Setup

Data for the analyses presented here were taken from PHENIX at RHIC. In
particular, the measurements presented are from the 2015 running year with
$\sqrt{s}=200$ GeV $p+p^{\uparrow}$ initial state, and inclusive final state
particles measured at midrapidity $\mid\eta\mid<0.35$.

###  2.1 Charged Particle Identification at PHENIX

The PHENIX detector contains two central arm spectrometers with acceptance
$\mid\eta\mid<0.35$, and $\Delta\phi=0.5$ per arm. Each consist of a drift
chamber and pad chambers for tracking and momentum measurements, ring-imaging
Cherenkov (RICH) detectors for particle identification, and electromagnetic
calorimeters (EMCal) for measuring energy deposition and triggering on charged
particles. In addition, for the open heavy flavor electron analysis, the hit
pattern in the silicon vertex detector close to the beam line is measured in
order to veto conversion electrons. All of these components allow for curation
of samples of charge particles at PHENIX, and additional information from the
RICH and EMCal subsystems can be used to separate samples of charged hadrons
from electrons. The RICH detector has a Cherenkov threshold of $\gamma=35$,
providing $e$, $\pi$ separation from $20$ MeV/c $<p_{T}<5$ GeV/c. Furthermore,
electrons deposit all or most of their energy in the EMCal, while charged
hadrons only deposit a fraction of their energy, allowing one to use the $E/p$
distribution of charged tracks to improve electron purity.

###  2.2 Analysis Procedure

The TSSA observable $A_{N}$ is calculated with the relative luminosity
formula, integrating over the $\phi$ ranges of the east and west spectrometer
arm due to the limited detector acceptance,

$A_{N}=\frac{1}{\langle\mid\cos\phi\mid\rangle}\frac{1}{P}\frac{N_{L}^{\uparrow}-RN_{L}^{\downarrow}}{N_{L}^{\uparrow}+RN_{L}^{\downarrow}}.$
(3)

The $\uparrow,\downarrow$ superscripts represent spin up vs spin down bunch
crossings, the $L$ subscript represents the the left spectrometer arm, $P$ is
the average polarization fraction of the beam,
$\langle\mid\cos\phi\mid\rangle$ is a correction factor applied for
integrating over $\phi$, and
$R=\mathcal{L}^{\uparrow}/\mathcal{L}^{\downarrow}$ is the relative
luminosity, defined as the ratio of luminosity for spin up bunch crossings vs
spin down bunch crossings. Note that there is an equivalent formula to
Equation 3 for the opposite spectrometer arm that comes with an additional
factor of $-1$. As mentioned in Section 2.1, the central arm spectrometers
measure charged particle yields. In addition to this, the luminosity
$\mathcal{L}$ is measured by the beam beam counter (BBC), and beam
polarization $P$ with the RHIC polarimeters, giving all the necessary
ingredients for Equation 3. Once $A^{S+B}_{N}$ is calculated for an inclusive
sample, it is necessary to quantify the signal to background ratio in each bin
as well as the background asymmetries $A_{N}^{B_{i}}$ in order to extract the
signal asymmetry $A^{S}_{N}$. For the heavy flavor electron analysis, the
backgrounds consist of electrons from other sources
($\pi^{0},\eta,\gamma,J/\psi,K^{0}_{S}/K^{\pm}$), and charged hadrons
(primarily $\pi^{\pm}$), while for the $\pi^{\pm}$ analysis, the backgrounds
consist of mostly electrons, with additional contributions from heavier
hadrons ($p/\bar{p},K^{\pm}$). However, these heavier hadrons contribute
negligibly in the measured range of the $\pi^{\pm}$ results
($5$GeV/c$<p_{T}<15$GeV/c). An electron cocktail method is used to estimate
background fractions of various electronic decay channels in simulation, and
data driven methods can be applied using information recorded in the EMCal and
RICH detector subsystems to quantify the $e^{\pm}/h^{\pm}$ fraction. The EMCal
and RICH selection requirements can be tuned to balance electron purity vs.
selection efficiency. The background asymmetries for the heavy flavor analysis
were taken from [8, 9] for the photonic sources ($\pi^{0},\eta,\gamma$), from
[10] for the hadron contamination $h^{\pm}$, and [11] for the $J/\psi$, while
$K^{0}_{S}/K^{\pm}$ was measured to be negligible in the measured $p_{T}$
range, and was not considered in the background correction. The background
asymmetry for the $J/\psi$ is a large source of statistical uncertainty for
the final heavy flavor electron results. For this reason, a non photonic
electron asymmetry was also measured that only takes
$(\pi^{0},\eta,\gamma,h^{\pm})$ into account in the background correction
procedure.

##  3 Heavy Flavor Electron TSSA

Heavy flavor production at midrapidity is dominated by the gluon-gluon fusion
channel at $\sqrt{s}=200$ GeV [12], receiving only a small contribution from
$q\bar{q}$ annihilation in the measured range $1.0$ GeV/c $<p_{T}<5.0$ GeV/c
(note that the RICH provides good $e/\pi$ separation in that range). Open
charm in the final state is also dominant at these kinematics. The gluon-gluon
fusion channel is of particular interest because the twist 3 trigluon (ggg)
correlator is not well constrained from previous measurements, while the
Efremov-Teryaev-Qiu-Sterman (qgq) correlator is somewhat constrained from
previous measurements. In addition, since the gluon is a massless gauge boson,
its transversity distribution is zero, isolating the term involving the
trigluon correlator in Equation 2. It should be noted that the qgq correlator
is used as an input to extract information on the ggg correlator. Figure 1
shows the results of the charge combined analysis, with heavy flavor $e^{\pm}$
$A_{N}(p_{T})$ in blue and nonphotonic $e^{\pm}$ $A_{N}(p_{T})$ in green, both
are the most statistically precise measurements of this observable at
midrapidity. It can be seen that both heavy flavor and nonphotonic $e^{\pm}$
$A_{N}(p_{T})$ results are consistent with zero in the measured range.

Figure 1: Heavy flavor $e^{\pm}$ $A_{N}(p_{T})$ (blue), nonphotonic $e^{\pm}$
$A_{N}(p_{T})$ (green)

Figure 2 shows the results of the charge separated analysis, with heavy flavor
$e^{+}$ $A_{N}(p_{T})$ in blue and nonphotonic $e^{+}$ $A_{N}(p_{T})$ in green
in the left panel and, and similarly for $e^{-}$ $A_{N}(p_{T})$ in the right
panel. It can again be observed that both heavy flavor and nonphotonic
$e^{+,-}$ $A_{N}(p_{T})$ results are consistent with zero in the measured
range.

Figure 2: Left panel: heavy flavor $e^{+}$ $A_{N}(p_{T})$ (blue), nonphotonic
$e^{+}$ $A_{N}(p_{T})$ (green). Right panel: Heavy flavor $e^{-}$
$A_{N}(p_{T})$ (blue), nonphotonic $e^{-}$ $A_{N}(p_{T})$ (green)

Most interestingly, Figure 3 shows the results of the charge separated
analysis again for only the heavy flavor $e^{+,-}$ samples, this time plotted
alongside theoretical results from [12]. $A_{N}^{D^{0},\bar{D}^{0}}(p_{T})$
were calculated at midrapidity in [12] and PYTHIA simulations were conducted
to convert these asymmetries to $A_{N}^{D^{0}\rightarrow
e^{+},\bar{D}^{0}\rightarrow e^{-}}(p_{T})$. In particular, the left panel
shows $e^{+}$ $A_{N}(p_{T})$ in blue plotted alongside theoretical predictions
of $A_{N}^{D^{0}\rightarrow e^{+}}(p_{T})$ from [12] \+ PYTHIA simulations for
various parameter combinations $(\lambda_{f},\lambda_{d})$, while the right
panel shows $e^{-}$ $A_{N}(p_{T})$ in blue plotted alongside theoretical
predictions of $A_{N}^{\bar{D}^{0}\rightarrow e^{-}}(p_{T})$ from [12] \+
PYTHIA simulations those same parameter combinations
$(\lambda_{f},\lambda_{d})$. The $\lambda$ parameters correspond to
normalization factors of trigluon correlators
$\phi^{(3),(f,d)}_{g/A}(x,x,\vec{s_{\perp}})$ at a particular fixed point in
x, to unpolarized gluon PDFs $g(x)$, with
$\phi^{(3),f}_{g/A}(x,x,\vec{s_{\perp}})$ = $\lambda_{f}g(x)$ and
$\phi^{(3),d}_{g/A}(x,x,\vec{s_{\perp}})$ = $\lambda_{d}g(x)$. The theoretical
curves are ordered differently for the different charges, as shown in Figure
3, allowing for some sensitivity to constrain $\lambda_{f}$ and $\lambda_{d}$
by looking at the charges separately.

Figure 3: Left panel: heavy flavor $e^{+}$ $A_{N}(p_{T})$ (blue), theoretical
predictions of $A_{N}^{D^{0}\rightarrow e^{+}}(p_{T})$ from [12] \+ PYTHIA
simulations (red). Right panel: heavy flavor $e^{-}$ $A_{N}(p_{T})$ (blue),
theoretical predictions of $A_{N}^{\bar{D}^{0}\rightarrow e^{-}}(p_{T})$ from
[12] \+ PYTHIA simulations (red).

##  4 Charged Pion TSSA

Charged pion production at midrapidity is dominated by quark-gluon scattering
at $\sqrt{s}=200$ GeV in the measured range $5$ GeV/c $<p_{T}<15$ GeV/c,
potentially allowing for sensitivity to light quark flavors when looking at
$\pi^{\pm}$ separately. In addition, $(\pi^{\pm},\pi^{0})$ is an isospin
triplet; comparing $A_{N}$ for these different particles is a good test for
theoretical models. The $\pi^{0}$ $A_{N}(p_{T})$ was previously measured at
$\sqrt{s}=200$ GeV at midrapidity [8]. This is compared with new results for
$\pi^{\pm}$ $A_{N}(p_{T})$ in Figure 4, with the red points representing
$\pi^{-}$ $A_{N}(p_{T})$, blue points representing $\pi^{+}$ $A_{N}(p_{T})$,
and gray points representing $\pi^{0}$ $A_{N}(p_{T})$ from [8]. Both the
$\pi^{0}$ $A_{N}(p_{T})$ and $\pi^{\pm}$ $A_{N}(p_{T})$ results are consistent
with zero within the measured range, but the $\pi^{\pm}$ $A_{N}(p_{T})$
results indicate that $\pi^{+}$ and $\pi^{-}$ may behave differently,
potentially exhibiting a dependence on quark flavor.

Figure 4: $\pi^{-}$ $A_{N}(p_{T})$ (red), $\pi^{+}$ $A_{N}(p_{T})$ (blue),
$\pi^{-}$ $A_{N}(p_{T})$ from [8] (gray)

##  5 Conclusion

TSSAs provide access to non-perturbative spin-momentum and spin-spin
correlations within the proton. While both TMD factorization and twist 3
collinear factorization can explain large measured TSSAs, twist 3 collinear
factorization requires only a single scale (the hard scale), for which the
measured particle’s $p_{T}$ can be taken as a proxy. Measurements presented in
this note can therefore be used to constrain twist 3 correlators. Midrapidity
heavy flavor $e^{\pm}$ and $\pi^{\pm}$ TSSAs from PHENIX were presented for
$p+p^{\uparrow}$ collisions at $\sqrt{s}=200$ GeV. The heavy flavor $e^{\pm}$
TSSA is the most precise to date, and is consistent with zero in the measured
range for both the charge combined and separated samples. The charge separated
heavy flavor $e^{+,-}$ TSSAs were compared with theoretical predictions from
[12], and will serve useful in providing constraints to the trigluon
correlator. The $\pi^{\pm}$ TSSAs are also consistent with zero in the
measured range, but the different charges seem to display different trends,
which may point to flavor sensitivity of the separate charges. The $\pi^{\pm}$
TSSAs were also compared with previously measured $\pi^{0}$ TSSAs from [8].
Both presented results are in preparation for publication, and provide
additional insight on gluon spin-momentum correlations inside protons

## Acknowledgements

Thank you to Zhongbo Kang for providing the parametrizations for the trigluon
correlator curves from [12] shown in Figure 3.

#### Funding information

The funding for this project was provided by the Department of Energy, grant
number DE-SC0013393.

## References

  * [1] PHENIX Collaboration (Aidala, C. et al.), Cross section and transverse single-spin asymmetry of muons from open heavy-flavor decays in polarized p+p collisions at $\sqrt{s}=200$ GeV, Phys. Rev. D 95, 112001 (2017), 10.1103/PhysRevD.95.112001.
  * [2] G. L. Kane, J. Pumplin, and W. Repko, Transverse Quark Polarization in Large$-p_{T}$ Reactions, $e^{+}e^{-}$ Jets, and Leptoproduction: A Test of Quantum Chromodynamics, Phys. Rev. Lett. 41, 1689 (1978), 10.1103/PhysRevLett.41.1689.
  * [3] C.A. Aidala, S.D. Bass, D. Hasch, and G.K. Mallot, The spin structure of the nucleon, Rev. Mod. Phys. 85, 655 (2013), 10.1103/RevModPhys.85.655.
  * [4] C. Aidala, Spin-Momentum Correlations, Aharonov-Bohm, and Color Entanglement in Quantum Chromodynamics, Colloquium at KIT (2019) https://www.phenix.bnl.gov/WWW/publish/caidala/Colloquia/ColloqKIT2019.pdf.
  * [5] PHENIX Collaboration (Aidala, C. et al.), Single-spin asymmetry of $J/\psi$ production in $p+p$, $p+Al$, and $p+Au$ collisions with transversely polarized proton beams at $\sqrt{s_{NN}}=200$ GeV, Phys. Rev. D 98, 012006 (2018), 10.1103/PhysRevD.98.012006.
  * [6] D. Pitonyak, Transverse spin observables in hard-scattering hadronic processes within collinear factorization, Int. J. Mod. Phys. A Vol. 31 No. 32 1630049 (2016), 10.1142/S0217751X16300490
  * [7] Metz, A., Pitonyak, D., Schäfer, A. et al., Transverse Single-Spin Asymmetries: Challenges and Recent Progress, Few-Body Syst 56, 331-336 (2015), 10.1007/s00601-014-0929-9
  * [8] PHENIX Collaboration, Transverse single-spin asymmetries of midrapidity $\pi^{0}$ and $\eta$ mesons in polarized p+p collisions at $\sqrt{s}=200$ GeV, Phys. Rev. D 103, 052009 (2021), 10.1103/PhysRevD.103.052009.
  * [9] PHENIX Collaboration, Transverse Single-Spin Asymmetries of Midrapidity Direct Photons, Neutral Pions, and Eta Mesons in 200 GeV Polarized Proton-Proton Collisions at PHENIX, https://arxiv.org/abs/2008.04283.
  * [10] PHENIX Collaboration (S. S. Adler et al.), Measurement of Transverse Single-Spin Asymmetries for Midrapidity Production of Neutral Pions and Charged Hadrons in Polarized $p+p$ Collisions at $\sqrt{s}=200$ GeV, Phys. Rev. Lett. 95 , 202001 (2005), 10.1103/PhysRevLett.95.202001
  * [11] PHENIX Collaboration (Adare, A. et al.), Measurement of transverse single-spin asymmetries for $J/\psi$ production in polarized $p+p$ collisions at $\sqrt{s}=200$GeV, Phys. Rev. D 82 112008 (2010), Erratum: Phys. Rev. D 86 (2012), 10.1103/PhysRevD.82.112008.
  * [12] Z.B Kang, J.W. Qiu, W. Vogelsang, and F. Yuan, Accessing trigluon correlations in the nucleon via the single spin asymmetry in open charm production, Phys. Rev. D 78, 114013 (2008), 10.1103/PhysRevD.78.114013.


# Network topology of interlocked chiral particles

Paul A. Monderkamp<EMAIL_ADDRESS>Institut für Theoretische Physik
II: Weiche Materie, Heinrich-Heine-Universität Düsseldorf, 40225 Düsseldorf,
Germany  Rika S. Windisch  Institut für Theoretische Physik II: Weiche
Materie, Heinrich-Heine-Universität Düsseldorf, 40225 Düsseldorf, Germany
René Wittmann  Institut für Theoretische Physik II: Weiche Materie, Heinrich-
Heine-Universität Düsseldorf, 40225 Düsseldorf, Germany  Hartmut Löwen
Institut für Theoretische Physik II: Weiche Materie, Heinrich-Heine-
Universität Düsseldorf, 40225 Düsseldorf, Germany

###### Abstract

Self-assembly of chiral particles with an L-shape is explored by Monte-Carlo
computer simulations in two spatial dimensions. For sufficiently high packing
densities in confinement, a carpet-like texture emerges due to the
interlocking of L-shaped particles, resembling a distorted smectic liquid
crystalline layer pattern. From the positions of either of the two axes of the
particles, two different types of layers can be extracted, which form distinct
but complementary entangled networks. These coarse-grained network structures
are then analyzed from a topological point of view. We propose a global charge
conservation law by using an analogy to uniaxial smectics and show that the
individual network topology can be steered by both confinement and particle
geometry. Our topological analysis provides a general classification framework
for applications to other intertwined dual networks.

##  I Introduction

The response of any liquid crystal to external constraints, such as
confinement, intricately depends on the density but crucially also on the
geometry of the particles. In fluids of uniaxial rod-like particles, for
instance, the particles tend to align at intermediate densities, forming a so-
called nematic phase. At the same time, the rods favor certain alignments with
confining walls, such that the material becomes continuously deformed, to
balance the competition between these two factors Yao et al. (2018); Yao and
Chen (2020); Lewis et al. (2014); Geigenfeind et al. (2015). Upon increasing
packing fraction, the liquid crystal tends to transition into a smectic phase,
where the particles additionally stack into layers. Consequently, confinement
typically leads to a fragmentation into several domains, separated by grain
boundaries Wittmann et al. (2021); Basurto et al. (2020); Monderkamp et al.
(2021, 2022).

Figure 1:  Circularly confined L-particles on different levels of coarse-
graining according to the yellow arrows. (a): Simulation snapshot. The long
axis of the L is colored in red, the short axis in green. (b): The coarse-
grained layer structures of both axes are visualized (long axis layers:
magenta, short axis layers: green). (c): Further coarse-grained structures
with their vertices (dots). We assign the vertices a topological charge $q$,
according to the number of adjacent edges (blue: $q>0$, yellow $q<0$). The
number of these charges $q$ serves the quantification of disorder in the
confined systems.

Advances in the synthesis of molecular and colloidal liquid crystals, enables
the study of systems with increasing geometrical complexity of the particles.
Non-convex particles such as polygonal rings, banana-shaped particles and
colloidal alphabet soups allow for geometric interlocking, increasing the
rigidity of the material Avendaño et al. (2016); Hernandez and Mason (2007);
Niori et al. (1996); Heppke et al. (2000); Dingemans and Samulski (2000); Ros
et al. (2005). In particular, particles with a characteristic chiral shape can
exhibit interlocking at high densities providing an ideal playground for a
wealth of interesting ordered structures Kamien and Selinger (2001); Harris et
al. (1997, 1999); Pollard et al. (2019); Meyer (1977); Dierking (2014); Hoell
and Löwen (2016).

In this work, we study the topology of smectic-like layers for interlocked
chiral particles under the influence of confinement (see Fig. 1). To this end,
with the aid of Monte-Carlo simulations, we generate systems of particles with
the chiral shape of the letter L (see Fig. 2), which are confined to circular
and annular cavities. We observe the emergence of highly complex structures
through the interlocking of the chiral, non-convex particles (see Fig. 1.(a)).
In particular, both particle axes display a tendency for layering, while the
competition between these two rigidly connected building blocks prevents the
formation of regular layered patterns found in ordinary smectics. Each layer
associated with one axis of the L-shaped particles can be interpreted as a
network such that the whole confined liquid crystal can be understood as an
interwoven structure of two coexisting species of networks (see Fig. 1.(b)).
We analyze the topology of the systems by only considering these coarse-
grained networks and assigning an index to the vertices, depending on the
number of adjacent edges (see Fig. 1.(c)).

The algorithm for the extraction of these layers is designed to create an
output analogous to the topological picture of conventional uniaxial smectics.
There, topological charge conservation is guaranteed by explicitly considering
smectic layers (density peaks) and so-called half-layers in between (density
minima) as topological entities, that carry topological charge Machon et al.
(2019); Aharoni et al. (2017); Beller et al. (2014); Hocking et al. (2022). In
that way, the indices of the vertices become topological, i.e., adhere to
topological charge conservation in analogy to conventional electrodynamics,
where the total charge, consisting of inside and boundary charge, adjusts to
the topology of the confining container. In that way, the circle and the
annulus yield different total topological charges, due to their different
connectivity.

Moreover, we show that the structure of the individual networks can be largely
steered by the particle shape. The vertical axis of the letter L is denoted as
leg, while the horizontal axis is denoted as foot. In particular, the spatial
distribution as well as the amount of inside topological charges depends
delicately on the foot to leg ratio
$p=L_{\text{\text{f}}}/L_{\text{\text{l}}}$ (see Fig. 2) of the L-shaped
particle. Specifically, the total amount of charge within the interior depends
non-monotonically on $p$.

Finally, We expect that our general topological treatment of the convoluted
network structures is also of practical value as a classification framework
for other systems, where intertwined dual network structures can be found,
such as gyroid phases in celestial nuclear matter Schütrumpf et al. (2015);
Nakazato et al. (2009), technical applications Crossland et al. (2009), blood
vessels in living organisms or traffic networks Barthélemy and Flammini
(2006); Papageorgiou (1990).

The paper is arranged as follows. In Sec. II.1, we present our simulation
protocol, while we elaborate on the graph-theoretical foundations of the
network topology in Sec. II.2. We present our results in Sec. III, before we
conclude in Sec. IV.

##  II Methods

###  II.1 Simulations

We perform canonical Monte-Carlo simulations for particles that have the shape
of the letter L (see Fig. 2). We confine the particles to two-dimensional
spherical and annular cavities. The particles are modeled as a pair of rigidly
connected almost hard discorectangles with equal diameters $D$, and core
lengths $L_{\text{\text{f}}}$, $L_{\text{\text{l}}}$, expressed in units of
$D$. Throughout this work, we vary
$p=L_{\text{\text{f}}}/L_{\text{\text{l}}}$, while
$L_{\text{\text{f}}}+L_{\text{\text{l}}}=28D$ stays constant. The interactions
of the L-shaped particles are conveniently specified through those of the
individual rod-like building blocks.

For any two rods, $i$ and $j$, not constituting to the same L-particle, with
positions $\mathbf{r}_{i}$, $\mathbf{r}_{j}$ and orientations
$\mathbf{\hat{u}}_{i}$ and $\mathbf{\hat{u}}_{j}$, the pair potential is
defined as harmonic repulsion

$U_{i,j}(\mathbf{r}_{i},\mathbf{r}_{j},\mathbf{\hat{u}}_{i},\mathbf{\hat{u}}_{j})=\left\\{\begin{array}[]{ll}U_{0}(\tau)(D-d_{ij})^{2}&\textnormal{
}d_{ij}<D\\\ 0&\textnormal{ }d_{ij}\geq D,\\\ \end{array}\right.$ (1)
$d_{i,j}=\min_{\left|\alpha,\beta\right|<\frac{L}{2}}\left\|\mathbf{r}_{i}+\alpha\mathbf{\hat{u}}_{i}-(\mathbf{r}_{j}+\beta\mathbf{\hat{u}}_{j})\right\|,$
(2)

with $\alpha,\beta\in[-L/2,L/2]$, with $L$ in
$\\{L_{\text{\text{f}}},L_{\text{\text{l}}}\\}$ defines the overlap Vega and
Lago (1994) and $U_{0}(\tau)$ is the interaction coefficient, which is
linearly increasing as a function of the simulation progress $\tau\in[0,1]$,
i.e., fraction of completed total of $10^{6}$ MC-cycles. $U_{0}(\tau=1)$ is
chosen as $10^{3}k_{\text{B}}T$ to mimic almost hard repulsion in the
equilibrated system.

Figure 2:  Schematic of the particles with the shape of the letter L as
considered in this work. The lengths of the two axes are referred to as
$L_{\text{\text{f}}}$ (foot of the letter L, green) and $L_{\text{\text{l}}}$
(leg of the letter L, red), respectively. The ratio of the axes is denoted by
$p=L_{\text{\text{f}}}/L_{\text{\text{l}}}$.

We model the wall interaction by representing the rods as three virtual point-
like particles sitting at the ends, and middle of the $k$-th rod
$\mathbf{r}_{\lambda}=\mathbf{r}_{k}+\lambda\mathbf{\hat{u}}_{k}$, with
$\lambda\in\\{-L/2,0,L/2\\}$. The interaction potential

$V(x)=\left\\{\begin{array}[]{ll}V_{0}\left(\tau\right)\,x^{2}&\textnormal{for
}\mathbf{r}_{\lambda}\textnormal{ outside,}\\\ 0&\textnormal{for
}\mathbf{r}_{\lambda}\textnormal{ inside the cavity }\\\ \end{array}\right.$
(3)

with the walls is once again harmonic, where $x$ denotes the closest distance
of $\mathbf{r}_{\lambda}$ to the inside of the respective cavity. Similarly to
$U_{0}$, $V_{0}$ is linearly increased as a function of the completed Monte-
Carlo cycles $\tau\in[0,1]$ up to as $10^{3}k_{\text{B}}T$ to mimic hard
walls. We simulate the systems at packing fractions
$\eta_{\text{1}}=Na_{\text{L}}/a_{\text{cav}}=0.4$ and $0.6$, with the area of
a single L-shaped particle $a_{\text{L}}$ and the area of the cavity
$a_{\text{cav}}$. The radii of the confinements are typically in the range
$4L\lesssim R\lesssim 50L$. Correspondingly, typical particle numbers $N$ are
between several hundreds and several thousands. In the annular confinement, we
keep the width of the annular ring constant ($\Delta R\approx
1.83L_{\text{\text{l}}}$), while we vary the inner radius ($R_{\text{in}}$)
and outer radius ($R_{\text{in}}+\Delta R$).

We follow a compression protocol, where we randomly initialize the system at a
low volume fraction $\eta_{\text{0}}=0.1\eta_{\text{1}}$. Each MC-cycle
consists of a trial displacement or rotation of each particle. The acceptance
probability $P=\min(1,\exp(-\Delta U/k_{\text{B}}T))$, for any trial move, is
given by the Metropolis criterion from the difference $\Delta U$ of the
energies (see Eqs. (1) and (3)) in the system before and after trial
Metropolis et al. (1953). Over the course of the simulation, we gradually
rescale the positions of the walls and particles to increase $\eta$ according
to

$\eta(\tau)=\eta(\tau)=(\eta_{\text{1}}-\eta_{\text{0}}){\tau}^{\frac{1}{3}}+\eta_{\text{0}}$
(4)

until the volume fraction $\eta_{\text{1}}=0.4$ is reached. This type of
decelerating compression aids the equilibration speed, since the system is
quickly compressed in the dilute stage, while being allowed to undergo a
larger fraction of MC-cycles in the regime, where self-assembly of the ordered
structures is expected (for more details on the protocol, see Appendix A).

###  II.2  Network topological charge analysis

As illustrated in Fig. 2, we denote the vertical (red) axis by leg which we
distinguish from the horizontal foot (green) axis by presuming a parallel wall
alignment of the leg (see Appendix B for more details). As visible in Fig. 1,
each ensemble of smectic-like layers associated with either foot or leg
(perceived as individual rods) forms convoluted networks and correspondingly
can be analyzed with the help of mathematical graph theory.

Each network consists of a set of vertices . Those vertices are connected via
a set of edge lines Euler (1758); Alama (2008).

Figure 3:  Schematic of the network topological model (see Eq. (5)). (a):
Serving charge conservation, our algorithm is designed such that there can not
be empty loops. As such, any loop contains a network of the respective other
species. The total network topological charge in the system is conserved under
this condition, since the charges of a loop and a simply connected network
cancel. (b): The total charge $q$, within any area of a single network, is
calculated from the number of in-/outgoing edges through the boundary,
inserted into Eq. (5). Dangling ends, such as within the dashed circle, are
charge neutral. Vertices with two adjacent edges, carry no charge, and are
therefore not explicitly labeled. Non charge-neutral operations are (i) adding
loops ($\Delta q=-1$) or (ii) new networks ($\Delta q=+1$) adds to the total
charge within the system. (c): A vertex with four outgoing edges
($\mathpzc{d}=4$, $q=-1$ see Eq. 5) can be understood as two infinitesimally
close vertices with $\mathpzc{d}=3$, $q=-1/2$, respectively.

As known from the treatment of the topology of layers in conventional smectics
Machon et al. (2019); Aharoni et al. (2017); Beller et al. (2014); Hocking et
al. (2022), a charge conservation follows, if the species of networks
alternate. In other words, between any two smectic layers has to be a density
minimum, i.e., a half-layer. Accordingly, our graph generation is designed,
such that empty loops are contracted into a single vertex (see Appendix B).
Therefore, we only observe the occurrence of loops (see Fig. 1), that each
envelop a simply connected graph (without loops) of the other species. This is
schematically visualized in Fig. 3.(a). More seldom, multiple loops of
alternating species are nestled into each other, with a simply connected graph
in the middle. To characterize this behavior in general, we assign a
topological network charge to any vertex in the network as

$q=1-\frac{\mathpzc{d}}{2}$ (5)

with the degree $\mathpzc{d}$, i.e., the adjacent number of edges. Note, that
this definition of the network charges is analogous to the layer topological
charges, typically considered in smectic liquid crystals, where the edge lines
represent smectic layers (see Appendix C). As such, the total charge of a
network species reads as

$Q_{a}=\sum_{\begin{subarray}{c}\text{vertices}\\\ \text{in
$a$}\end{subarray}}q$ (6)

with $a\in\\{\text{f},\text{l}\\}$, standing for foot and leg. where
$Q_{\text{\text{f}}}+Q_{\text{\text{l}}}$ is a conserved quantity.

Figure 4:  Results of three independent simulation runs for L-shaped particles
confined to spherical cavities with packing fraction $\eta=0.4$ and (a):
$p=0.4$, (b): $p=2.5$ and the hard-rod limit (c): $p\approx 0$. The upper row
(1) displays particle snapshots superimposed with the layer networks of both
species (see sec. II.2). The lower row (2) shows the isolated networks. The
negative network charges are indicated in brown, while the positive are shown
in blue. In all cases, the boundary charges are assigned to the magenta leg-
network. Therefore the inside charge of the foot-species, matches the
respective total charge
($Q_{\text{\text{f}}}^{\text{(in)}}=Q_{\text{\text{f}}}$), which is not true
for the leg-species ($Q_{\text{\text{l}}}^{\text{(in)}}\neq
Q_{\text{\text{l}}}$). Accordingly (a) and (b) correspond to a color swap,
only in the inside of the cavity. The total network charge in the system
$Q_{\text{\text{f}}}+Q_{\text{\text{l}}}$ matches the Euler characteristic of
the confinement $\chi=1$.

The total charge of a network is only determined by its connectivity and not
by the total number of vertices. As visible from Fig. 3.(b), adding a dangling
end to a previously existing network, is a charge neutral operation.
Consequently, any simply connected network, carries the charge of a single
isolated vertex $q=1$. Only the addition of loops, i.e., adding an edge
between two existing vertices, decreases the net charge by $1$. This network
topological charge can only be conserved, if every loop coincides with a
simply connected (loop-less) network, e.g., at its center.

This definition of the charge is consistent, such that the charge within any
area, can be calculated from the number of in-/outgoing edges, similar to
Gauss’s law in classical electrodynamics (see Fig. 3.(c)). Any vertex of
higher degree can be viewed as a composition of merged $q=-1/2$ charges.
Vertices with $\mathpzc{d}=2$ carry no charge, and can be therefore, together
with the edges, viewed as constituting the layers, in between the rest of the
charges. Respecting these properties of the network charge, our algorithm
which generates the final networks (see Appendix B) is designed to
systematically eliminate vertices, such that the final network structure, as
shown in Fig. 1.(c), both illustrates the network connectivity and allows us
to properly determine the total charges $Q_{\text{\text{f}}}$ and
$Q_{\text{\text{l}}}$.

Typically, considering the conservation of topological charges in confined
geometries, e.g., orientational topological defects in nematically ordered
fluids Kim et al. (2013); Monderkamp et al. (2021), requires the
identification of boundary defects on the system walls if no global alignment
condition is obeyed. Due to the invariance of the total charge within the
cavity $Q=Q_{\text{\text{f}}}+Q_{\text{\text{l}}}$ (determined by the topology
of the confining domain), we have the liberty to assign the outer walls of the
confinement to any of the two network species. In the following, without loss
of generality, we choose to assign the boundary defects to the layer network
of the leg, generating the boundary charge $Q_{\text{l}}^{\text{(b)}}$. We
will denote the inside network charges, i.e., charges without explicit
inclusion of boundary charge, as $Q_{\text{\text{f}}}^{\text{(in)}}$ and
$Q_{\text{\text{l}}}^{\text{(in)}}$. As detailed in Sec. III, the total sum

$Q_{\text{\text{f}}}^{\text{(in)}}+Q_{\text{\text{l}}}^{\text{(in)}}=Q-Q_{\text{l}}^{\text{(b)}}$
(7)

is not conserved. Still, these quantities contain structural information about
the confined state. Through our choice of assigning the boundary charges to
the leg-network, clearly
$Q_{\text{\text{f}}}^{\text{(in)}}=Q_{\text{\text{f}}}$ but
$Q_{\text{\text{l}}}^{\text{(in)}}\neq
Q_{\text{\text{l}}}=Q_{\text{\text{l}}}^{\text{(in)}}+Q_{\text{l}}^{\text{(b)}}$.

##  III Results

###  III.1 Circular confinement

For circular cavities, we show simulation results for systems of L-shaped
particles, with three different ratios of axes lengths $p$ in Fig. 4. at
packing fraction $\eta=0.4$. In Fig. 4.(a1) a typical snapshot, for $p=2.5$,
is depicted. The particles locally interlock, while the longer foot-axes
display a strong tendency for alignment, leading to elongated clusters. Global
orientational ordering, however, is not visible, such as would be expected
from, e.g., a conventional smectic liquid crystal. Through the interlocking,
the shorter leg-axes fill the spaces in between. Below, in Fig. 4.(a2), we
show the corresponding graph network, resulting from our analysis of the
layers. The elongated foot-clusters are represented in Fig. 4.(a2) by the
green network. The leg-network is depicted in magenta. The relative higher
stiffness of the wider foot-layers results in a favoring of longer simply-
connected networks, each contained in a loop of a single large magenta leg-
network. This results in the presence of majorly negative defects (indicated
in brown), due to the loops. Conversely, every separated component of the
foot-network contains a charge $\sum q=1$. As elaborated in Sec. II.2, this
results in a total network charge within the system of $Q=1$. Similarly, Fig.
4.(b) shows a simulation snapshot for the inverse length ratio $p=0.4$.
Accordingly, the behavior of the networks within the interior of the
confinement, visible in Fig. 4.(b2) is analogous to the former case with
inverted species. The assignment of the boundary charges remains with the leg-
network, as in Fig. 4.(a). As the total charge is invariant of our choice of
this alignment condition, the inclusion of the boundary charges still retains
total network charge within the confined system as $Q=1$. Finally, Fig. 4.(c)
shows a simulation snapshot, where $p\approx 0$, i.e., hard rods. The confined
system resides in a state, where global orientational ordering is present.
Through the lack of interlocking, no strong indications of layers are visible,
as expected from the unconfined nematic bulk phase which is stable at the
chosen packing fraction $\eta=0.4$ Bates and Frenkel (2000); Bolhuis and
Frenkel (1997). Nevertheless, the network analysis can still be used to
quantify the global topology. Naturally, the total charge is still conserved
as $Q=1$.

Figure 5:  Simulation results of liquid crystals composed of L-shaped
particles for a range of different length ratios $p$ of both particle axes,
illustrating a topological transition in the networks at $p=1$, where the
chirality of the particles flips. (a): Inside charges
$Q_{\text{\text{f}}}^{\text{(in)}}$ (green),
$Q_{\text{\text{l}}}^{\text{(in)}}$ (magenta) and sum
$Q_{\text{\text{f}}}^{\text{(in)}}+Q_{\text{\text{l}}}^{\text{(in)}}$ (dashed
black) (b): Global orientational order parameter $S$. We observe, that at each
point of the horizontal axes, the respective longer axis network has a
positive charge, and the shorter axis network a negative charge. Above
approximately $0.2$, the systems reside in interlocked layered states without
global orientational order. Below $0.2$, the systems converge against
conventional nematic order. The results are shown at packing fraction
$\eta=0.4$ .The lines correspond to averages over $75$ simulations per $p$.

In Fig. 5.(a) we show the inside charges $Q_{\text{\text{f}}}^{\text{(in)}}$
and $Q_{\text{\text{l}}}^{\text{(in)}}$ of both network species as well as
their sum, i.e., the total charge without the inclusion of boundary charges,
as a function of $p$. Additionally, along the same horizontal axis, in Fig.
5.(b), we show the global orientational order parameter
$S=|\left<\exp(i2\phi)\right>|$, where $\left<...\right>$ denotes an average
over all particles. The data shown in both figures can be roughly divided into
three characteristic regions along the horizontal axis, namely $p>1$ (cf. Fig.
4.(a)), $0.2\lesssim p<1$ (cf. Fig. 4.(b), as well as $p\lesssim 0.2$ (cf.
Fig. 4.(c)). In the case $p>1$, the simply-connected, but isolated foot-layers
(green) result in positive inside charge. At the same time, the single large
leg-network with loops results in negative charges (magenta). Analogously, the
same holds true in the regime $0.2\lesssim p<1$, only with inverted network
species. Therefore, the signs of $Q_{\text{\text{f}}}^{\text{(in)}}$ and
$Q_{\text{\text{l}}}^{\text{(in)}}$ are flipped. For $p\lesssim 0.2$, the
shapes of the particles are approaching the hard-rod limit. The absence of an
entropic contribution from the interlocking mechanism of the L-shaped
particles, allows the liquid crystal to undergo a transition into a nematic
state with global orientational ordering. Correspondingly, Fig. 5.(b) shows an
increased orientational order parameter $S$. In the transition regime,
complicated packing effects dominate the system, causing non-trivial behavior
in $Q_{\text{\text{f}}}^{\text{(in)}}$ and
$Q_{\text{\text{l}}}^{\text{(in)}}$. We observe that the first peak in
$Q_{\text{\text{l}}}^{\text{(in)}}$ decreases for larger systems, where the
interior of the confinement is less influenced by the system boundaries (not
shown). We therefore infer, that the behavior in Fig. 5 stems from the extreme
influence exerted by the confinement. Due to the symmetry of the particles,
all the observables are symmetric around $1$, i.e., $S(p)=S(1/p)$, as well as
$Q_{\text{\text{l}}}^{\text{(in)}}(p)=Q_{\text{\text{f}}}^{\text{(in)}}(1/p)$
(see Appendix D). Naturally, the total charges $Q_{\text{\text{f}}}$ and
$Q_{\text{\text{l}}}$ do not adhere to this symmetry, which is broken by the
assignment of the boundary charges to the leg-network.

In general, the total inside charge
$Q_{\text{\text{f}}}^{\text{(in)}}+Q_{\text{\text{l}}}^{\text{(in)}}$ within
the cavity, i.e., the sum of all network charges in Fig. 5 is constantly
greater than zero. This is consistent with our observation, of the presence of
isolated networks within the inside, without the existence of empty loops.
Moreover, including the boundary to the leg-network does by construction only
add negative defects, i.e., $Q_{\text{l}}^{\text{(b)}}<0$, compare, e.g., Fig.
4.(2). Taking a closer look at the behavior of
$Q_{\text{\text{f}}}^{\text{(in)}}+Q_{\text{\text{l}}}^{\text{(in)}}$ in Fig.
5.(a), we notice two local maxima. The first one at $p\approx 0.1$ coincides
with the onset of global orientational order, i.e., close to the transition
between confined nematic and interlocked layer states. The second one is
located around $p=1$, i.e., where neither network dominates the structure,
such that the mutual interruption of the layering of the two building blocks
is most pronounced. We thus conclude that the total inside charge
$Q_{\text{\text{f}}}^{\text{(in)}}+Q_{\text{\text{l}}}^{\text{(in)}}$ provides
a useful measure for disorder in interlocked or frustrated systems.

###  III.2 Annular confinement

Figure 6:  Simulation results for a system of L-shaped particles in annular
confinement. (a): The system resides in a smectic-$C$-like state, where due to
the interlocking, the layer director is at an angle to the local orientation.
(b): The network structure roughly shows one foot-network per smectic layer.
The boundary charges are assigned to the leg-network, therefore the network
charge is conserved. Through the empty loop on the inner wall, the total
network charge matches the Euler characteristic $\chi=0$ of the confinement.
Packing fraction and axes length ratio are $\eta=0.6$ and $p=0.1$. Figure 7:
Network topological charge of the foot-species
$Q_{\text{\text{f}}}=-Q_{\text{\text{l}}}$ within annular confinement as a
function of the inner radius $R_{\text{in}}$. Packing fraction and particle
dimensions match Fig. 6. We show results of $20$ simulations per considered
$R_{\text{in}}$ and a linear fit. The slope of the line $m$ divided by $2\pi$
is approximately one, confirming that there is slightly less than one charge
per particle length around the perimeter. This represents, as expected, one
smectic-like layer per particle length (cf. Fig. 6).

One of the important characteristics of a topological conservation law within
a confined liquid crystal system is the relevance of the topology of the
confining container Fulton (1995); Alexander et al. (2012). To further explore
this, we introduce a confining domain with annular shape, that possesses Euler
characteristic $\chi=0$. In Fig. 6.(a), a corresponding particle snapshot is
shown ($R_{\text{in}}=4L_{\text{\text{l}}}$, $\Delta R\approx
1.83L_{\text{\text{l}}}$, $\eta=0.6$, $p=0.1$). Through the high packing
fraction and the relatively long leg-axes (red), the particles show a visible
tendency to reside in locally parallel layers, where the long axes of the
particles are parallel align with the outer walls. At the same time, the
protruding foot-axes (green) cause a relative shift of neighboring particles,
resulting in a characteristic smectic-$C$-like shape, where the particles are
tilted with respect to the layers. The associated network structure is visible
in Fig. 6.(b). As in the previous section, the boundary charges are assigned
to the leg-network. Each smectic block of particles results in an isolated
simply connected foot-network, wrapped into a loop of the large leg-network.
The inner confinement walls cause an additional empty loop. As a result, all
charges sum up to the Euler characteristic of an annulus $\chi=Q=0$.

Figure 8:  Simulation results in annuli with large inner radii
$R_{\text{in}}=40L_{\text{\text{l}}}$. The other simulation parameters match
Fig. 6. (a): Excerpt from a simulation snapshot. The tilt angle of the layers
with respect to the walls of the annular ring is clearly visible. (b): Annular
pair distribution function $g_{\circ}(\mathbf{r})$ between the positions of
the foot-particles (green) in the $xy$-plane up to a radial distance of
$L_{\text{\text{l}}}$ (cf. Eq. (8)). The slanted peaks represent the slanted
foot-layers in (a). The angle relative to the horizontal is equal to $\approx
34.4^{\circ}$. (c): schematic of L-particles in close packing. For
$D_{\text{foot}}=D_{\text{leg}}=D$, the optimal tilt angles of the layers is
equal to $45^{\circ}$ (vertical and horizontal sides of the right triangle in
the figure have to be equal). The deviation in the measured angle in (b)
possibly stems from complex entropic interactions of the layer species at
lower packing fraction.

The previous results are for systems with relatively small particle numbers
$N=320$. We make this particular choice to illustrate the network charges on
the particle-resolved level. Naturally, all of the above holds true also for
larger systems. In Fig. 7, we present the total topological charge
$Q_{\text{\text{f}}}=Q_{\text{\text{f}}}^{\text{(in)}}=-Q_{\text{\text{l}}}$
of the foot-species as a function of the inner radius $R_{\text{in}}$ of the
annular confinement, while keeping the width of the annulus $\Delta
R=1.83L_{\text{\text{l}}}$ constant. Since the boundary charges are assigned
to the leg-network, the total charge $Q_{\text{\text{f}}}$ coincides with the
inside charge $Q_{\text{\text{f}}}^{\text{(in)}}$. Each depicted point
represents an individual simulation. As visible in the plot, the data points
lie fairly accurately on a straight line, indicated by a linear fit with slope
$m=0.912\times 2\pi$. More specifically, there is slightly less than one
positive charge per particle length along the inner circumference of the
annulus. This is consistent with the observations made in Fig. 6, where each
particle layer forms a new network. Since any particle layer has a width,
which is slightly larger than $L_{\text{\text{l}}}$, $m/(2\pi)$ is slightly
smaller than one.

To further investigate the origin of the smectic-$C$-like tilted layers for
these L-shaped particles, additional results are presented in Fig. 8. In
conventional hard-rod smectics, confined to annular cavities, the layers
typically align with the outer walls, while the direction of the layer is
typically in radial direction of the confinement, i.e., orthogonal to the
walls Wittmann et al. (2021); Yao and Chen (2020). Here, we observe layers
which are tilted with respect to the radial direction. This is nicely visible
in Fig. 6 and Fig. 8.(a). The latter shows an excerpt from a simulation in an
annulus with inner radius $R_{\text{in}}=40L_{\text{\text{l}}}$. We see
similar tilt angles in both snapshots, irrespective of the large difference in
curvature of the respective confinement walls. In order to show the local
structure within the annular ring, we introduce the two-dimensional annular
pair distribution function

$\displaystyle
g_{\circ}(\mathbf{r})=\frac{1}{N\rho}\left<\sum_{\begin{subarray}{c}i,j=1\\\
i\neq
j\end{subarray}}^{N}\delta(\mathbf{r}-\overline{\overline{\mathbf{R}}}\cdot\left(\mathbf{r}_{i}-\mathbf{r}_{j}\right))\right>.$
(8)

This is obtained via the matrix

$\displaystyle\overline{\overline{\mathbf{R}}}=\overline{\overline{\mathbf{R}}}(\phi_{i})=\left(\begin{array}[]{rr}\cos(\pi/2-\phi_{i})&-\sin(\pi/2-\phi_{i})\\\
\sin(\pi/2-\phi_{i})&\cos(\pi/2-\phi_{i})\\\ \end{array}\right),$ (11)

where $\phi_{i}$ is the polar angle of the position $\mathbf{r}_{i}$ of the
$i$-th particle with respect to the annulus center. Here, $N$ and $\rho$ are
the global particle numbers and densities, respectively, and
$\delta(\mathbf{r})$ denotes the delta-distribution. The physical
interpretation of $g_{\circ}(\mathbf{r})$ is a distribution function of
rotated relative vectors. The data to obtain $g_{\circ}(\mathbf{r})$ are
sampled from $20$ independent simulation runs. We compute
$g_{\circ}(\mathbf{r})$ for the positions of the green foot-particles, as
depicted in Fig. 8.(b). In the depiction, the horizontal axis denotes the
tangential direction of the annular walls. The center of the diagram shows an
anisotropic depletion zone around the particles, resulting from the almost
hard repulsion. Furthermore, slanted density peaks are visible, representing
the tilted layers within the annular ring. These peaks form an approximate
angle of $\alpha=55.6^{\circ}$ with the radial direction. In Fig. 8.(c), it is
illustrated that for identical width of the particle axes
($D_{\text{foot}}=D_{\text{leg}}=D$), an angle of $\alpha=45^{\circ}$ leads to
efficient packing at extreme packing fractions. At lower packing fractions,
however, this would lead to a drastic decrease of the free length out of the
layer, parallel to the longer leg-particles. More specifically, the measured
angle can be understood as a result of two competing entropic factors: (i) the
aforementioned free length along the leg orientation, which tends to increase
the tilt angle relative to the radial direction, and (ii) the systems tendency
to reside in a smectic-like layered structure perpendicular to the leg
orientation, due to the high density, which tends to decrease the tilt angle
(as presumably, the layered structure breaks apart for
$\alpha\rightarrow\pi/2$).

##  IV Conclusion

In this work, we present a formalism for the analysis of the topology of
liquid crystals with chiral particle shapes, which give rise to irregular
defect structures in confinement. This approach is based on the consideration
of the entangled network structure formed by the respective axes positions. In
order to analyze the structures with two layer species, we generate networks
utilizing Delaunay triangulations and identify topological charges, from the
degrees (number of adjacent edges) of the network vertices. This procedure
leans on the concept of layers and half-layers that characterize the topology
in conventional smectics Machon et al. (2019); Aharoni et al. (2017), as,
e.g., formed by uniaxial hard rods. Like these previous applications, our
study relies on coexisting but disjoint network species, which as a whole
adhere to a common topological conservation law, where the sum of the
respective topological charges in a confined system equals the Euler
characteristic $\chi$ of the finite container. However, there are two crucial
differences. (i) In our study, both network species explicitly refer to an
axis of the particles. In contrast, in conventional smectics, one layer
species represents the density minima in between the particle layers. (ii) Our
approach is based on the analysis of the layer network structure on the
particle-resolved level, while the observation of conventional smectic point
defects (see Appendix C) typically happens on much larger length scales and
was hitherto majorly employed in approaches that describe the coarse-grained
smectic layer structure Machon et al. (2019); Aharoni et al. (2017); Beller et
al. (2014); Hocking et al. (2022); Stannarius and Harth (2016); Harth and
Stannarius (2020). Moreover, we exemplify that the bare inside charges serve
as a useful quantifier of disorder.

We use this framework of network topology to analyze the structure of liquid
crystals composed of particles with the shape of the letter L, confined to
circular ($\chi=1$) and annular ($\chi=0$) cavities. These are generated via
the use of canonical Monte-Carlo simulation. We find that the global inside
charges of the two network species intricately depend on both the local
particle symmetries and the global degree of order. We observe a variety of
remarkable states at different packing fractions and L-axes lengths: at
intermediate densities, the particles prefer interlocking with suppression of
global orientational order, when both axes of the L-shape have comparable
sizes. Otherwise the system tends towards a conventional nematic state Yao et
al. (2018, 2021); Bates and Frenkel (2000). At high packing fractions in the
annulus, the particles arrange in more rigid layers, which resembles a
smectic-$C$ structure. In regards to the latter, we present additional
simulation for annuli with large radii elucidating the origin of the tilt
angle of the layers.

Based on this insight, we expect that the formalism used throughout this work
can positively contribute to a variety of future endeavors. In particular we
expect, that it will be useful in the interpretation of future computational,
theoretical and experimental studies of systems with complicated particle
shapes Martínez-Ratón and Velasco (2022); Zhao et al. (2007, 2012); Avendaño
and Escobedo (2017); Barmes et al. (2003); Kraft et al. (2013); Yuan et al.
(2018); Lapointe et al. (2009); Zerrouki et al. (2008). Even though we
introduce this method as a tool for the investigation of relatively complex
chiral particles, we also suspect that it will be insightful to apply it to
hard-rod fluids, since the network analysis can be readily applied to
visualize the half-layers in conventional smectics. We expect it to provide
also insight, applied in the analysis of gyroid phases in systems of particles
without chirality Schönhöfer et al. (2017, 2018). This will possibly be a
powerful tool in understanding the emergent topological structures Bramble et
al. (2007); Kim et al. (2009); Liarte et al. (2016). Recent years have also
seen an increased interest in non-equilibrium systems, where the individual
particles consume energy from their surrounding to propel themselves Tan et
al. (2019); Loewe and Shendruk (2021); Giomi (2015); Bär et al. (2020). It has
been shown, that in these systems, the orientational defects have dynamical
properties, which in turn depend on the respective topological charges DeCamp
et al. (2015); Tóth et al. (2002); Huang et al. (2022). It seems therefore
very reasonable to assume, that network topological charges will carry
interesting dynamical properties in collectively moving L-shaped or other
chiral swimmers Hernández et al. (2020); Kümmel et al. (2013); Caprini and
Marconi (2019); Löwen (2016).

## Acknowledgements

The authors would like to thank Michael te Vrugt and Raphael Wittkowski for
stimulating discussions. This work is funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) – LO 418/20-2.

## Author declarations

### Conflict of Interest

The authors have no conflicts to disclose.

### Author contributions

Paul A. Monderkamp: Conceptualization (equal); Data Curation (lead); Formal
Analysis (lead); Methodology (equal); Project Administration (equal); Software
(equal); Supervision (equal); Validation (lead); Visualization (lead);
Writing/Original Draft Preparation (lead); Writing/Review & Editing (equal);
Rika S. Windisch: Conceptualization (equal); Data Curation (supporting);
Formal Analysis (supporting); Methodology (equal); Project Administration
(supporting); Software (equal); René Wittmann: Conceptualization (equal);
Formal Analysis (supporting); Funding Acquisition (supporting); Methodology
(equal); Project Administration (equal); Supervision (equal); Writing/Review &
Editing (equal); Hartmut Löwen: Conceptualization (equal); Funding Acquisition
(lead); Methodology (equal); Project Administration (equal); Resources (lead);
Supervision (equal); Writing/Review & Editing (equal);

### Data availability

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

Figure 9:  Delaunay triangulation of the union of all coordinate positions of
both axes of all L-shaped particles, used as a base point for the generation
of the graphs networks, considered in this manuscript. The corresponding
particle snapshot is shown in Fig. 1. The edges, connecting the foot-positions
are colored in green, while the edges, connecting the leg-axes are colored
magenta. Grey edges connect foot\- and leg-vertices. In order to obtain lines,
representing the layer structure of the respective axes, the gray edges are
deleted, and the remaining colored networks are transformed into simple lines.

##  Appendix A Details on the equilibration protocol

In the main manuscript, the simulation results are generated with the help of
a canonical Monte-Carlo simulation, as described in Sec. II.1. The fundamental
goal of the procedure is to find a configuration at a relatively high packing
fraction which reflects the equilibrium state. Since it is practically
impossible to obtain a randomized configuration at the target packing fraction
(since guessing a configuration which fits into the cavity corresponds to
guessing the final simulation result), the systems are initialized several
orders of magnitude below the target packing fraction. They are subsequently
compressed with a decelerating compression rate (cf. Eq. (4)). Additionally,
the interaction constant $U_{0}$ (cf. Eq. (1)) between the particles and the
wall interaction constant $V_{0}$ (cf. Eq. (3)) are linearly increased as a
function of the fraction of completed Monte-Carlo cycles $\tau\in[0,1]$. The
initial softness of the interactions helps the particles heal-out overlaps,
which may occur in the beginning of the simulation due to random
initialization. At the end of the simulation, $U_{0}$ and $V_{0}$ are
$10^{3}k_{\text{B}}T$. We find that the $\tau$-dependence of these constants
becomes less relevant, as the simulation progresses, since the particles
effectively feature hard repulsion at a certain point. This positively
contributes to the equilibration.

To obtain configurations, which reflect the equilibrium configuration, without
explicitly evaluating free energies, the equilibration is performed slowly
enough, such that two assumptions can be made about the transient
configurations over the course of the simulation. (i): The system is ergodic
(samples the whole configuration space), in the stage of self-assembly, such
that the system configuration space at higher densities is also sampled fairly
over independent different simulation runs. (ii): The compression of the
system and increase of the interaction strengths occurs slowly enough, that
the system can equilibrate faster than the parameters change, such that the
final configuration reflects a state close to equilibrium.

Below, we show additional simulation results, which aim to support both
assumptions. The simulations are performed as described in Sec. II.1 with axis
length ratio $p=1$ and fraction $\eta_{1}=0.4$ inserted into the formula for
compression (Eq. (4)). But instead of equilibrating up to $\eta_{1}=0.4$, the
compression and increase of the interaction constants are stopped at the half-
point of the simulation $\tau=1/2$. This value coincides with $\eta\approx
0.32$ (cf. Eq. (4)) Subsequently, the simulation runs at constant parameters
for the second half to illustrate the change of the system without compression
at intermediate densities.

Figure 10:  Trajectory of a single particle within a Monte-Carlo simulation,
where the compression and change of interaction constants is stopped at half
the number of total Monte-Carlo cycles $\tau=1/2$. The simulation starts as
described in Sec. II.1, but runs with constant parameters after $\tau=1/2$.
The displayed particle positions are rescaled with the radius of the
confinement, i.e. the depicted trajectory shows $\mathbf{r}(\tau)/R(\tau)$.
The green arrow denotes the direction of motion within the trajectory. Note
that the trajectory reflects Monte-Carlo displacement (contrary to a standard
equation of motion). The relatively unconstrained motion suggests, that the
system is in principle able to occupy all configurations. The particle also
samples the whole orientation space (not shown).

###  A.1 Ergodicity

Here, we support the assumption, that the simulation protocol samples the
configuration space ergodically at intermediate densities. This is to
guarantee that the configuration space at higher densities is sampled fairly
across many simulations. To this end, we show the trajectory of the leg-axis
of a single particle over the course of the simulation in Fig. 10. The
displayed positions are rescaled with radius of the confinement $R(\tau)$
($R(\tau)^{2}\propto 1/\eta(\tau)^{1/2}$, cf. Eq. (4)). As such, the figure
shows $\mathbf{r}(\tau)/R(\tau)$. The fact that the particle moves throughout
the whole cavity strongly suggests, that the entire configuration space is
sampled fairly, since the particles can freely rearrange at intermediate
densities, and therefore the system is not hindered from occupying specific
configurations. Note that not every single particle has to be able to traverse
the whole cavity, since they can be considered physically indistinguishable in
terms of the states. We consider this numerical indication, that the system
arrests into a high-density state, which is close to the equilibrium and
independent of the initialization.

Figure 11:  Total network charge of the foot-network as a function of
simulation progress $\tau$, i.e., the fraction of completed Monte-Carlo
cycles. The compression and adjustment of interaction constants is stopped at
the half-way point, after which the simulation runs at constant parameters. We
pick six instances of simulation progress at which we average
$Q_{\text{\text{f}}}$. Those averages are denoted $\overline{\mathcal{Q}}_{f}$
(see red box). We perform a hypothesis test for equality of all six
$\overline{\mathcal{Q}}_{f}$, via a one-way $F$-test. We find an $F$-value of
approximately $0.70$ and a $p$-value of approximately $0.62$. We view this as
indication, that the system at $\tau=1/2$ already reflects an equilibrium
configuration well.

###  A.2 Slow compression

In a simulation protocol, such as described in Sec. II.1, one needs to
guarantee, that the final configuration reflects/is close to a true
equilibrium configuration. In practice, this corresponds to the execution of
the protocol slowly enough, such that the system does not get quenched into an
unfavorable kinetically arrested state. In particular, we execute the protocol
slowly enough, that all configurations throughout the simulation (beyond the
initial fast compression given by Eq. (4)) reflect a state close to
equilibrium for any instance of parameters across the simulation.

To reinforce this claim, we show simulation results in this section, where we
stop the compression and increase of the interaction constants, $U_{0}$ and
$V_{0}$, after half the number of Monte-Carlo cycles is completed
($\tau=1/2$). More specifically, we start the simulation with the same
equilibration protocol as described in Sec. II.1. We abruptly terminate the
change of parameters at $\tau=1/2$. We continue the Monte-Carlo simulation at
constant parameters and show that the foot-network charge
$Q_{\text{\text{f}}}$ does not significantly change afterwards. We do this, to
illustrate, that the system is close to equilibrium throughout the simulation.
In Fig. 11, we show distributions of $Q_{\text{\text{f}}}$ over $20$
independent simulations as a function of simulation progress $\tau$ (completed
number of Monte-Carlo cycles). We denote the average value of
$Q_{\text{\text{f}}}$ for a constant $\tau$ by $\overline{\mathcal{Q}}_{f}$.
We test whether the $\overline{\mathcal{Q}}_{f}$ vary significantly after
stopping the compression at $\tau=1/2$. To this end, we set up a hypothesis
test for the equality of all values with $\tau\geq 1/2$ Miller Jr (1997) (the
conditions for the hypothesis test: equal variance, approximate normality,
independence, were checked). The null hypothesis is given by $\mathrm{H}_{0}$:
“All measured $\overline{\mathcal{Q}}_{f}(\tau\geq 1/2)$ do not differ
significantly” i.e., the true average $Q_{\text{\text{f}}}$ from all possible
configurations, which can theoretically be obtained in simulation, is
constant. The alternate hypothesis is $\mathrm{H}_{a}$: “At least one pair is
not equal.”. If the alternate hypothesis is true, we have to assume that the
$\overline{\mathcal{Q}}_{f}(\tau\geq 1/2)$ vary significantly, and so the
simulation has not reached near-equilibrium. Through a one-way $F$-test
Virtanen et al. (2020), we find $F\approx 0.70$ and a $p$-value of $0.62$. The
value $F$ denotes the variance between $\overline{\mathcal{Q}}_{f}$ for
different $\tau$ divided by the variance for a constant $\tau$. The $p$-value
indicates that a distribution of $\overline{\mathcal{Q}}_{f}$ with such a
variance, or more, is expected to occur with probability $p=0.62$. We
conclude, that given this data, the null hypothesis $\mathrm{H}_{0}$ can not
be rejected with any reasonable significance level. We furthermore infer that,
in combination with the absence of a clear up- or downwards trend after
$\tau=1/2$, we have reason to believe, that the observables do not change
after stopping compression. This can be considered indication, that at
$\tau=1/2$, the system resides in a state which reflects equilibrium well.

##  Appendix B Network generation

The networks, considered in this work, are generated by considering the
Delaunay triangulation Aurenhammer et al. (2013) of the union of the
geometrical centers of both axes of the L-shaped particles (see Fig. 9,
corresponding to the particle snapshot in Fig. 1.(a)). The emerging network is
disconnected into the two network species, by deleting edges between opposing
foot and leg vertices (Fig. 9, gray). Furthermore, to obtain a physical
picture of these two intertwined networks as a pair of smectic-like layers
with a topological charge conservation, we require two further systematic
modifications. Firstly, we assign the boundary to the leg-network. This
corresponds to assuming that the leg of the particles align preferably
parallel with the wall, while other local configurations are interpreted as a
defect. Such a presumed uniform alignment rule stands at the basis of any
topological conservation law in confinement. Secondly, to ensure that the
final networks represent global layer structures, we transform the triangular
meshes into simple lines through merging vertices that form empty triangles
(without changing the hierarchy between the two networks). Thereby we delete
any empty loops which are not compatible with the concept of alternating
layers required for topological charge conservation. This is in accordance
with the layer and half-layer picture of conventional smectics and can be
readily applied to uniaxial rods by considering the limit
$L_{\text{\text{l}}}=0$. The final network, after applying this protocol to
Fig. 9, is shown in Fig. 1.(c).

Figure 12:  Schematic of an orientational topological defect in a smectic
liquid crystal. The defect is a point, where the local orientation
$\mathbf{\hat{n}}(\mathbf{r})$, typically orthogonal to the smectic layers,
indicated as black lines, is ill-defined. The rotation $\Delta\phi_{i}$ from
$\mathbf{A}$ to $\mathbf{B}$ is equal to $\alpha_{i}-\pi$. Therefore the total
rotation $\Phi$ is equal to $\sum_{i}^{n}\Delta\phi_{i}=2\pi(1-n/2)$, with the
number of outgoing layers $n$. The topological charge is equal to
$\Phi/2\pi=1-n/2$, in analogy to Eq. (5).

##  Appendix C Analogy to topological charges in smectics

The network topological charges, introduced in Sec. II.2 share a close
relation with the orientational topological charges, typically considered in
smectic liquid crystals Kamien and Mosna (2016); Chen et al. (2009); Kurik and
Lavrentovich (1988). In this manuscript, we assign a network topological
charge $q=1-\mathpzc{d}/2$ to any vertex in the observed network via its
degree $\mathpzc{d}$ (see Eq. (5)). Similarly, topological charges of defects
in smectic systems can be understood in terms of adjacent layers Machon et al.
(2019); Aharoni et al. (2017).

Figure 13:  Simulation results of liquid crystals composed of L-shaped
particles for a range of different length ratios $p$ of both particle axes as
Fig. 5 in the main text. Through the logarithmic scaling of the horizontal
axis, the symmetry around $p=1$ becomes apparent.

Orientational topological defects in liquid crystals, that display local
alignment of the molecules, can be understood as singular locations in space,
where the local orientation $\mathbf{\hat{n}}(\mathbf{r})$ undergoes a
discontinuous jump and is therefore ill-defined. In smectic systems, where the
particles additionally arrange in layers, this can happen across grain
boundaries in two and three dimensions, or across points in two dimensions. A
schematic of a two-dimensional point defect is depicted in Fig. 12. Smectic
layers are represented by black lines. In this particular example, five layers
meet in a singular point. Around this point, the layers, and thus
$\mathbf{\hat{n}}(\mathbf{r})$, typically at a constant angle to the layers,
are continuously bent. The strength of the defect is characterized by the
total rotation of $\mathbf{\hat{n}}(\mathbf{r})$ traversing the defect in
counterclockwise direction. Consider one wedge of the rotation
$\mathbf{A}\rightarrow\mathbf{B}$ between two layers at an angle $\alpha_{i}$:
The rotation angle is equal to $\Delta\phi_{i}=\alpha_{i}-\pi$. Consequently,
the total rotation around the defect is equal to

$\Phi=\sum_{i}^{n}\Delta\phi_{i}=2\pi(1-n/2),$ (12)

with the number of outgoing layers $n$. The topological charge of the defect
is defined by $\Phi/2\pi$ resulting in Eq. (5), with the vertex degree
$\mathpzc{d}$ identified with $n$.

##  Appendix D Larger axes length ratios

We denote the conventionally short (length $L_{\text{\text{f}}}$) horizontal
axis of the letter L by foot and the conventionally long (length
$L_{\text{\text{l}}}$) vertical axis by leg. The length ratio $p$ of the axes
can in practice vary between $0$ and $\infty$. If the symmetry between the two
axes is not broken, e.g., by assignment of boundary charges to any of the two
corresponding network species, the physical observables should generally be
symmetric around $p=1$. This is confirmed in Fig. 13, where (a) the inside
network charges and (b) the global orientational order parameter $S$ are
depicted. The figure features the simulation data depicted in Fig. 5
supplemented with simulation data for larger $p$. In particular, we verify
that $S(p)=S(1/p)$, as well as
$Q_{\text{\text{l}}}^{\text{(in)}}(p)=Q_{\text{\text{f}}}^{\text{(in)}}(1/p)$,
i.e., the symmetry around $p=1$ is clearly visible.

## References

  * Yao et al. (2018) X. Yao, H. Zhang, and J. Z. Chen, Phys. Rev. E 97, 052707 (2018).
  * Yao and Chen (2020) X. Yao and J. Z. Y. Chen, Phys. Rev. E 101, 062706 (2020).
  * Lewis et al. (2014) A. H. Lewis, I. Garlea, J. Alvarado, O. J. Dammone, P. D. Howell, A. Majumdar, B. M. Mulder, M. Lettinga, G. H. Koenderink, and D. G. A. L. Aarts, Soft Matter 10, 7865 (2014).
  * Geigenfeind et al. (2015) T. Geigenfeind, S. Rosenzweig, M. Schmidt, and D. de las Heras, J. Chem. Phys. 142, 174701 (2015).
  * Wittmann et al. (2021) R. Wittmann, L. B. G. Cortes, H. Löwen, and D. G. A. L. Aarts, Nat. Commun. 12, 623 (2021).
  * Basurto et al. (2020) E. Basurto, P. Gurin, S. Varga, and G. Odriozola, Phys. Rev. Research 2, 013356 (2020).
  * Monderkamp et al. (2021) P. A. Monderkamp, R. Wittmann, L. B. G. Cortes, D. G. A. L. Aarts, F. Smallenburg, and H. Löwen, Phys. Rev. Lett. 127, 198001 (2021).
  * Monderkamp et al. (2022) P. A. Monderkamp, R. Wittmann, M. Te Vrugt, A. Voigt, R. Wittkowski, and H. Löwen, PCCP 24, 15691 (2022).
  * Avendaño et al. (2016) C. Avendaño, G. Jackson, E. A. Müller, and F. A. Escobedo, PNAS 113, 9699 (2016).
  * Hernandez and Mason (2007) C. J. Hernandez and T. G. Mason, J. Phys. Chem. C 111, 4477 (2007).
  * Niori et al. (1996) T. Niori, T. Sekine, J. Watanabe, T. Furukawa, and H. Takezoe, J. Mater. Chem. 6, 1231 (1996).
  * Heppke et al. (2000) G. Heppke, D. Parghi, and H. Sawade, Liq. Cryst. 27, 313 (2000).
  * Dingemans and Samulski (2000) T. J. Dingemans and E. T. Samulski, Liq. Cryst. 27, 131 (2000).
  * Ros et al. (2005) M. B. Ros, J. L. Serrano, M. R. de La Fuente, and C. L. Folcia, J. Mater. Chem. 15, 5093 (2005).
  * Kamien and Selinger (2001) R. D. Kamien and J. V. Selinger, Journal of Physics: Condensed Matter 13, R1 (2001).
  * Harris et al. (1997) A. B. Harris, R. D. Kamien, and T. C. Lubensky, Phys. Rev. Lett. 78, 1476 (1997).
  * Harris et al. (1999) A. B. Harris, R. D. Kamien, and T. C. Lubensky, Rev. Mod. Phys. 71, 1745 (1999).
  * Pollard et al. (2019) J. Pollard, G. Posnjak, S. Čopar, I. Muševič, and G. P. Alexander, Phys. Rev. X 9, 021004 (2019).
  * Meyer (1977) R. B. Meyer, Mol. Cryst. Liq. 40, 33 (1977).
  * Dierking (2014) I. Dierking, Sym. 6, 444 (2014).
  * Hoell and Löwen (2016) C. Hoell and H. Löwen, J. Chem. Phys. 144, 174901 (2016).
  * Machon et al. (2019) T. Machon, H. Aharoni, Y. Hu, and R. D. Kamien, Commun. Math. Phys. 372, 525 (2019).
  * Aharoni et al. (2017) H. Aharoni, T. Machon, and R. D. Kamien, Phys. Rev. Lett. 118, 257801 (2017).
  * Beller et al. (2014) D. A. Beller, T. Machon, S. Čopar, D. M. Sussman, G. P. Alexander, R. D. Kamien, and R. A. Mosna, Phys. Rev. X 4, 031050 (2014).
  * Hocking et al. (2022) B. J. Hocking, H. S. Ansell, R. D. Kamien, and T. Machon, Proc. R. Soc. 478, 20210725 (2022).
  * Schütrumpf et al. (2015) B. Schütrumpf, M. A. Klatt, K. Iida, G. E. Schröder-Turk, J. A. Maruhn, K. Mecke, and P.-G. Reinhard, Phys. Rev. C 91, 025801 (2015).
  * Nakazato et al. (2009) K. Nakazato, K. Oyamatsu, and S. Yamada, Phys. Rev. Lett. 103, 132501 (2009).
  * Crossland et al. (2009) E. J. Crossland, M. Kamperman, M. Nedelcu, C. Ducati, U. Wiesner, D.-M. Smilgies, G. E. Toombes, M. A. Hillmyer, S. Ludwigs, U. Steiner, et al., Nano Lett. 9, 2807 (2009).
  * Barthélemy and Flammini (2006) M. Barthélemy and A. Flammini, JSTAT 2006, L07002 (2006).
  * Papageorgiou (1990) M. Papageorgiou, TRANSPORT RES B-METH 24, 471 (1990).
  * Vega and Lago (1994) C. Vega and S. Lago, Comput. Chem. 18, 55 (1994).
  * Metropolis et al. (1953) N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953).
  * Euler (1758) L. Euler, Novi Comment. Acad. Sci. I. Petropolitanae 4, 109 (1758).
  * Alama (2008) J. Alama, Formaliz. Math. 16, 7 (2008).
  * Kim et al. (2013) Y.-K. Kim, S. V. Shiyanovskii, and O. D. Lavrentovich, J. Condens. Matter Phys. 25, 404202 (2013).
  * Bates and Frenkel (2000) M. A. Bates and D. Frenkel, J. Chem. Phys. 112, 10034 (2000).
  * Bolhuis and Frenkel (1997) P. Bolhuis and D. Frenkel, J. Chem. Phys. 106, 666 (1997).
  * Fulton (1995) W. Fulton, _Algebraic Topology - A First Course_ (Springer-Verlag, 1995).
  * Alexander et al. (2012) G. P. Alexander, B. G.-g. Chen, E. A. Matsumoto, and R. D. Kamien, Rev. Mod. Phys. 84, 497 (2012).
  * Stannarius and Harth (2016) R. Stannarius and K. Harth, Phys. Rev. Lett. 117, 157801 (2016).
  * Harth and Stannarius (2020) K. Harth and R. Stannarius, Front. Phys. 8, 112 (2020).
  * Yao et al. (2021) X. Yao, L. Zhang, and J. Z. Y. Chen, arXiv:2112.07889 (2021).
  * Martínez-Ratón and Velasco (2022) Y. Martínez-Ratón and E. Velasco, Phys. Fluids 34, 037110 (2022).
  * Zhao et al. (2007) K. Zhao, C. Harrison, D. Huse, W. Russel, and P. Chaikin, Phys. Rev. E 76, 040401 (2007).
  * Zhao et al. (2012) K. Zhao, R. Bruinsma, and T. G. Mason, Nat. Commun. 3, 801 (2012).
  * Avendaño and Escobedo (2017) C. Avendaño and F. A. Escobedo, COCIS 30, 62 (2017).
  * Barmes et al. (2003) F. Barmes, M. Ricci, C. Zannoni, and D. Cleaver, Phys. Rev. E 68, 021708 (2003).
  * Kraft et al. (2013) D. J. Kraft, R. Wittkowski, B. Ten Hagen, K. V. Edmond, D. J. Pine, and H. Löwen, Phys. Rev. E 88, 050301 (2013).
  * Yuan et al. (2018) Y. Yuan, A. Martinez, B. Senyuk, M. Tasinkevych, and I. I. Smalyukh, Nat. Mater 17, 71 (2018).
  * Lapointe et al. (2009) C. P. Lapointe, T. G. Mason, and I. I. Smalyukh, Science 326, 1083 (2009).
  * Zerrouki et al. (2008) D. Zerrouki, J. Baudry, D. Pine, P. Chaikin, and J. Bibette, Nature 455, 380 (2008).
  * Schönhöfer et al. (2017) P. W. Schönhöfer, L. J. Ellison, M. Marechal, D. J. Cleaver, and G. E. Schröder-Turk, Interface Focus 7, 20160161 (2017).
  * Schönhöfer et al. (2018) P. W. Schönhöfer, D. J. Cleaver, and G. E. Schröder-Turk, J. Phys. D 51, 464003 (2018).
  * Bramble et al. (2007) J. P. Bramble, S. D. Evans, J. R. Henderson, T. J. Atherton, and N. J. Smith, Liq. Cryst. 34, 1137 (2007).
  * Kim et al. (2009) Y. H. Kim, D. K. Yoon, M.-C. Choi, H. S. Jeong, M. W. Kim, O. D. Lavrentovich, and H.-T. Jung, Langmuir 25, 1685 (2009).
  * Liarte et al. (2016) D. B. Liarte, M. Bierbaum, R. A. Mosna, R. D. Kamien, and J. P. Sethna, Phys. Rev. Lett. 116, 147802 (2016).
  * Tan et al. (2019) A. J. Tan, E. Roberts, S. A. Smith, U. A. Olvera, J. Arteaga, S. Fortini, K. A. Mitchell, and L. S. Hirst, Nat. Phys. 15, 1033 (2019).
  * Loewe and Shendruk (2021) B. Loewe and T. N. Shendruk, arXiv:2111.07364 (2021).
  * Giomi (2015) L. Giomi, Phys. Rev. X 5, 031003 (2015).
  * Bär et al. (2020) M. Bär, R. Großmann, S. Heidenreich, and F. Peruani, Annu. Rev. Condens. Matter Phys. 11, 441 (2020).
  * DeCamp et al. (2015) S. J. DeCamp, G. S. Redner, A. Baskaran, M. F. Hagan, and Z. Dogic, Nat. Mater. 14, 1110 (2015).
  * Tóth et al. (2002) G. Tóth, C. Denniston, and J. M. Yeomans, Phys. Rev. Lett. 88, 105504 (2002).
  * Huang et al. (2022) Z.-F. Huang, H. Löwen, and A. Voigt, Commun. Phys. 5, 294 (2022).
  * Hernández et al. (2020) R. J. Hernández, F. J. Sevilla, A. Mazzulla, P. Pagliusi, N. Pellizzi, and G. Cipparrone, Soft Matter 16, 7704 (2020).
  * Kümmel et al. (2013) F. Kümmel, B. Ten Hagen, R. Wittkowski, I. Buttinoni, R. Eichhorn, G. Volpe, H. Löwen, and C. Bechinger, Phys. Rev. Lett. 110, 198302 (2013).
  * Caprini and Marconi (2019) L. Caprini and U. M. B. Marconi, Soft Matter 15, 2627 (2019).
  * Löwen (2016) H. Löwen, EPJ ST 225, 2319 (2016).
  * Miller Jr (1997) R. G. Miller Jr, _Beyond ANOVA: basics of applied statistics_ (CRC press, 1997).
  * Virtanen et al. (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, et al., Nat. Methods 17, 261 (2020).
  * Aurenhammer et al. (2013) F. Aurenhammer, R. Klein, and D.-T. Lee, _Voronoi diagrams and Delaunay triangulations_ (World Scientific Publishing Company, 2013).
  * Kamien and Mosna (2016) R. D. Kamien and R. A. Mosna, New J. Phys. 18, 053012 (2016).
  * Chen et al. (2009) B. G.-g. Chen, G. P. Alexander, and R. D. Kamien, Proc. Natl. Acad. Sci. U.S.A. 106, 15577 (2009).
  * Kurik and Lavrentovich (1988) M. V. Kurik and O. Lavrentovich, Phys.-Uspekhi 31, 196 (1988).