# Conversations with Flaschka: Kac-Moody groups and Verblunsky coefficients

Mohammad Javad Latifi<EMAIL_ADDRESS>Department
of Mathematics, Dartmouth College, Hanover, NH 03755  and  Doug Pickrell
<EMAIL_ADDRESS>Mathematics Department, University of Arizona,
Tucson, AZ 85721

###### Abstract.

In this paper we discuss two items which in one way or another originated from
conversations with Hermann Flaschka and his students. The first is an
application of the Toda lattice to the question of whether there exists a
complex Lie group (with certain properties) associated to an indefinite type
Kac-Moody Lie algebra. The second concerns a new example of the Verblunsky
correspondence.

##  0\. Preface

This paper is dedicated to the memory of Hermann Flaschka. Hermann made
seminal contributions to completely integrable systems and many related areas
of mathematics, reflecting his broad interests. He was beloved as a colleague
and friend for his humility, wisdom and wit. His intensity in pursuit of facts
and solutions, whether related to mathematical problems or raises for adjuncts
or climate change, was legendary.

Many decades ago, after a sanitized introduction to abstract mathematics, I
(the second author) took a course from Hermann which emphasized explicit
formulas for tau functions and solitons. This was a formative experience. A
decade ago Hermann, Bole Yang (Hermann’s last PhD student) and I spent a
semester learning about the Verblunsky correspondence, CMV matrices and the
discrete nonlinear Schroedinger equation, using the analogy with Jacobi
matrices and the Toda lattice as a guide.

##  1\. Introduction

The goal of this paper is to illustrate some ideas which in one way or another
originated from some conversations with Hermann Flaschka and his students.

###  1.1. Kac-Moody Groups

Kac-Moody algebras fall into three classes: finite dimensional, e.g.
$sl(n,\mathbb{C})$, which is the setting for Flaschka’s original reformulation
of the nonperiodic Toda equations; affine, e.g. $sl(n,\mathbb{C}[z,z^{-1}]))$
(matrices with finite Laurent series coefficients), which for example are
relevant to the periodic Toda equations; and indefinite. Kac and Peterson
associated an ‘algebraic group’ to a general Kac-Moody algebra. For the affine
and indefinite classes, this algebraic group lacks an exponential map, hence
it is natural to ask whether there exists a Lie group (with an exponential
map) which is associated to a completion of the algebra, e.g. $C^{0}$ or
$C^{1}$ or $C^{\infty}(S^{1},SL(n,\mathbb{C}))$ in the affine case. Toda
equations can be associated to the generalized Cartan matrix for any Kac-Moody
algebra, and in conventional cases solving the equations involves
exponentiation and LDU type factorization. Our main point is that properties
of the Toda equation solutions are intimately related to the question of
whether there exists a complex Lie group with exponential map associated to a
given Kac-Moody algebra. The dissertation related work of Hermann’s student
Maorong Zou suggests that there does not exist a complex Lie group, having
certain properties, in the indefinite cases, but to our knowledge there does
not exist a definitive theorem.

######  Remark 1.

In his work on moonshine conjectures and the Monster finite simple group,
Borcherds enlarged an indefinite type Lie algebra to the so called Monster Lie
algebra (see [4] for a short introduction). For recent work on the question of
how to associate groups to this algebra and others, see [3] and the
references.

###  1.2. The Verblunsky Correspondence

Let $Prob(S^{1})$ denote the metrizable compact convex set of probability
measures on the unit circle $S^{1}$, let $Prob^{\prime}(S^{1})$ denote the
subset consisting of measures with infinite support, and let
$\Delta:=\\{z\in\mathbb{C}:|z|<1\\}$. The Verblunsky correspondence refers to
a miraculous homeomorphism

(1.1)
$Prob^{\prime}(S^{1})\leftrightarrow\prod_{n=1}^{\infty}\Delta:\mu\leftrightarrow\alpha$

or more generally to a homeomorphism of $Prob(S^{1})$ and a compactification
of the product space. This correspondence can be described in several
different ways. One formulation is as follows. Given $\mu\in
Prob^{\prime}(S^{1})$, if $p_{0}=z^{0}$, $p_{1}(z)$, $p_{2}(z)$,… are the
monic orthogonal polynomials corresponding to $\mu$, then
$\alpha_{n}=p_{n}(0)^{*}$ (where $(\cdot)^{*}$ is complex conjugation).
Conversely, given $\alpha$, $\mu$ is the weak* limit of probability measures

(1.2)
$d\mu=\lim_{N\to\infty}\frac{\prod_{n=1}^{N}(1-|\alpha_{n}|^{2})}{|p_{N}(z)|^{2}}\frac{d\theta}{2\pi}$

where $p_{0}=1$,

(1.3) $p_{n}(z)=zp_{n-1}(z)+\alpha_{n}^{*}z^{n-1}p_{n-1}^{*}(z),\quad n>0,$

and $p_{n}^{*}(z)=p_{n}(\frac{1}{z^{*}})^{*}$.

######  Remarks.

(a) This correspondence is truly ‘miraculous’ because it yields a
parametrization of (nontrivial) measures on the circle by parameters which are
functionally independent - this has consequences, see e.g. [8].

(b) For a masterful survey of the history and ramifications of the Verblunsky
theory, see Simon’s two volume work, [13], or for a short survey, see [14].
For more recent connections with completely integrable systems, see [10] and
[15].

There are a number of known explicit examples of this correspondence, see
section 1.6 of Volume I of [13], including the case $\beta=1$ of the following

######  Theorem 1.1.

Suppose that $\beta>0$. The probability measure
$\mu=\frac{1}{\mathfrak{Z}}(1-cos(\theta))^{1/\beta}\frac{d\theta}{2\pi}$
corresponds to the sequence $\alpha_{n}=\frac{1}{1+n\beta}$, $n=1,2,...$.

This example illuminates the relation between decay properties of the
coefficients and regularity properties of the corresponding measures, which we
will highlight in the text.

It is a simple matter to use the definition of the Verblunsky coefficients to
compute a few of the coefficients; presumably we could have discovered this
example using these calculations. However for us this arose from mulling the
following formula,

(1.4)
$\sum_{i,j}\prod_{u=1}^{L}\frac{1}{(i(u)\beta+1)(j(u)\beta+1)}=\prod_{k=1}^{n}(\frac{1}{k}\beta^{-1}+\frac{k-1}{k})$

where the sum is over all integral sequences $i(1)>j(1)>...>i(L)>j(L)\geq 0$
satisfying $\sum_{u=1}^{L}(i(u)-j(u))=n$, for some $1\leq L\leq n$ (see [8]).
This leads to ‘super-telescoping formula’ (see [9])

(1.5)
$\sum_{n=0}^{\infty}\left(\sum_{i,j}\prod_{u=1}^{L}\frac{1}{(i(u)\beta+1)(j(u)\beta+1)}\right)z^{n}=\frac{1}{(1-z)^{1/\beta}}$

It turns out that inverting the Verblunsky map is related to ‘root subgroup
factorization’ (for an introduction to root subgroup factorization, see [1]).
Theorem 1.1 is essentially a corollary of the following

######  Theorem 1.2.

Suppose that $\alpha_{n}=1/(1+n\beta)$ as before and let $m=1/\beta$. Then

(1.6) $\lim_{n\to\infty}\left(\begin{matrix}1&\alpha_{n}z^{-n}\\\
\bar{\alpha}_{n}z^{n}&1\end{matrix}\right)...\left(\begin{matrix}1&\alpha_{1}z^{-1}\\\
\bar{\alpha}_{1}z&1\end{matrix}\right)=\left(\begin{matrix}\frac{1+F_{m}^{*}(z)}{2(1-z^{*})^{m}}&\frac{1-F_{m}^{*}(z)}{2(1-z^{*})^{m}}\\\
\frac{1-F_{m}(z)}{2(1-z)^{m}}&\frac{1+F_{m}(z)}{2(1-z)^{m}}\end{matrix}\right)$

where

$F_{m}(z)=1+\sqrt{\pi}2^{1-m}\frac{\Gamma(m+1)}{\Gamma(m+\frac{1}{2})}\sum_{n=1}^{\infty}\sum_{k=0}^{\infty}\left(\begin{matrix}m\\\
n+2k\end{matrix}\right)\left(\begin{matrix}n+2k\\\
k\end{matrix}\right)(\frac{-1}{2})^{n+2k}z^{n}$

(If $m$ is an integer, this is a polynomial; in general the binomial
coefficient is defined in terms of gamma functions).

One interesting thing to observe is that this infinite product is not bounded,
reflecting the non-summability of the coefficients. This has consequences.

###  1.3. Plan of the Paper

The prerequisites for Section 2 are steep. Our main goals are to given an
accessible and motivated introduction to questions about groups associated to
Kac-Moody algebras, and to explain why properties of generalized Toda
equations are inextricably linked to questions about Kac-Moody groups.

Section 2 is more straightforward. Our main purpose is to prove Theorems 1.1
and 1.2, and some related calculations.

###  1.4. Notation

If $f(z)=\sum f_{n}z^{n}$, then we will write

$f=f_{-}+f_{0}+f_{+}$

where $f_{-}(z)=\sum_{n<0}f_{n}z^{n}$ and $f_{+}(z)=\sum_{n>0}f_{n}z^{n}$,
$(f)_{0+}=f_{0}+f_{+}$, and $f^{*}(z)=\sum(f_{-n})^{*}z^{n}$, where
$w^{*}=\bar{w}$ is the complex conjugate of the complex number $w$. If the
Fourier series is convergent at a point $z\in S^{1}$, then $f^{*}(z)$ is the
conjugate of the complex number $f(z)$. If $f\in H^{0}(\Delta)$, then
$f^{*}\in H^{0}(\Delta^{*})$, where $\Delta$ is the open unit disk,
$\Delta^{*}$ is the open unit disk at $\infty$, and $H^{0}(U)$ denotes the
space of holomorphic functions for a domain $U\subset\mathbb{C}$.

Our conventions regarding Verblunsky coefficients differ from those in [13],
where the nontrivial coefficients are indexed by $n=0,1,...$. We use the
convention that $\alpha_{0}=1$ and the nontrivial Verblunsky coefficients are
$\alpha_{1},\alpha_{2},...$.

##  2\. An Application of Toda Lattice Theory to Kac-Moody Groups

Chevalley and Serre famously discovered that a finite dimensional simple
complex Lie algebra can be defined in terms of generators and relations using
a Cartan matrix $A$, i.e. an integral $n\times n$ matrix having positive
principal minors and satisfying $A_{ii}=2$, $A_{ij}\leq 0$ and
$A_{ij}=0\implies A_{ji}=0$ for $i\neq j$. Kac and Moody independently deleted
the condition that $A$ has positive principal minors, and discovered that the
theory remained robust (one can also allow $n$ to be infinite, which is
essential in Borcherds work on the Monster). The corresponding Kac-Moody
algebra $\mathfrak{g}=\mathfrak{g}(A)$ is generated by $3n$ letters $e_{i},$
$f_{i},$ and $h_{i}$ $(1\leq i<n+1)$ satisfying the Chevalley-Serre relations

$[h_{i},h_{j}]=0,$ $[e_{i},f_{i}$ $]=h_{i}$, $[e_{i}$ $,f_{j}$ $]=0,$ $i\neq
j$

$[h_{i}$ $,e_{j}]=A_{ij}e_{j}$, $[h_{i},f_{j}]=-A_{ij}f_{j},$

$(ad(e_{i}))^{1-A_{ij}}e_{j}=0,$ $(ad(f_{i}))^{1-A_{ij}}f_{j}=0,$ $i\neq j$

The main take away is this: for each $i=1,...,n$, there is an isomorphism
$sl(2,\mathbb{C})\leftrightarrow span(f_{i},h_{i},e_{i})$,

(2.1) $\left(\begin{matrix}0&0\\\ 1&0\end{matrix}\right)\leftrightarrow
f_{i},\quad\left(\begin{matrix}1&0\\\ 0&-1\end{matrix}\right)\leftrightarrow
h_{i},\quad\left(\begin{matrix}0&1\\\ 0&0\end{matrix}\right)\leftrightarrow
e_{i}$

and the other relations detail how these $n$ copies of $sl(2,\mathbb{C})$ (the
‘root subalgebras’ corresponding to simple roots) interact. A technical aside:
We will assume that $A$ is irreducible and symmetrizable, meaning that there
exists an essentially unique nondegenerate symmetric bilinear form on
$\mathfrak{g}$, denoted $\langle\cdot,\cdot\rangle$.

As in the finite dimensional case, the Lie algebra $\mathfrak{g}$ has a
triangular decomposition

$\mathfrak{g}=\mathfrak{n}^{-}\oplus\mathfrak{h}\oplus\mathfrak{n}^{+}$

where $\mathfrak{n}^{-}$ is generated by $\\{f_{j}\\}$ (think lower triangular
matrices, or creation operators), $\mathfrak{h}$ by $\\{h_{j}\\}$ (think
diagonal matrices), and $\mathfrak{n}^{+}$ by $\\{e_{j}\\}$ (think upper
triangular matrices, or annihilation operators).

As we noted in the introduction, Kac-Moody Lie algebras fall into three
classes: (1) $\mathfrak{g}$ is finite dimensional, in which case there exists
a corresponding complex Lie group $G(A)$. (2) $\mathfrak{g}$ is of affine
type, in which case there exists a finite dimensional simple complex Lie
algebra $\dot{\mathfrak{g}}$ such that $\mathfrak{g}$ is a central extension
of the loop algebra of functions $S^{1}\to\dot{\mathfrak{g}}$ having finite
Fourier series (or a twisted version of this); in this case there exist
infinite dimensional Lie groups corresponding to completions of
$\mathfrak{g}$, e.g. (a central extension of) the loop group
$C^{\infty}(S^{1},\dot{G}(A))$ (see [11]). And finally (3) $\mathfrak{g}$ is
of indefinite type.

######  Question.

In the indefinite case, does there exist a corresponding complex Lie group?

This is an interesting question, because for example in conventional
situations one uses the global group structure to solve associated
differential equations, such as the Toda equations.

Kac and Peterson showed in general how to construct an ‘algebraic group’
$G(A)$ which is associated to $\mathfrak{g}(A)$ (see [5] and [6]). The basic
idea is simple: $G(A)$ is generated by $n$ ‘root subgroups’, i.e. subgroups
isomorphic to $SL(2,\mathbb{C})$ corresponding to the $n$ copies of
$sl(2,\mathbb{C})$ in (2.1)). One imposes all the relations that arise from
realizing these groups in a highest weight representation of the Lie algebra
(meaning that the $f_{i}$ are creation operators and the $e_{i}$ are
annihilation operators). The paper [6] is a masterful analysis of the
essential relations, which in the words of the authors are ’analytic
continuations of the relations defining the Weyl group’, which are Coxeter
group relations.

$G(A)$ is not a Lie group, except in the finite case. $G(A)$ is useful in
formulating global algebraic notions, such as LDU factorization: for $g\in
G(A)$, $g=ldu$, where $l\in N^{-}$, $d\in H$, $u\in N^{+}$ (the groups
corresponding to $\mathfrak{n}^{-},\mathfrak{h},\mathfrak{n}^{+}$,
respectively) if and only if $\sigma_{j}(g)\neq 0$, where the $\sigma_{j}$ are
the fundamental matrix coefficients, $\sigma_{j}(g)=\langle g\cdot v,v\rangle$
where $v$ is the vacuum for the $i$th fundamental highest weight
representation (For $g\in SL(n,\mathbb{C})$ this is the statement that $g$ has
a LDU factorization iff the principal minors are nonsingular, and for $g\in
SL(n,\mathbb{C}[z,z^{-1}])$ this is essentially the statement that g has a
Riemann-Hilbert factorization (with factors having finite Laurent expansions,
a nontrivial condition) iff the associated Toeplitz operator is invertible; in
all cases $d=\prod_{j=1}^{n}\sigma_{j}(g)^{h_{j}}$, and - what is essential
for us - the $l$ and $u$ factors are rational in any given representation).
But the ordinary exponential function is transcendental, and $G(A)$ does not
have an exponential map, except in the finite case. We would like to complete
$G(A)$ so that it has an exponential map, while preserving its algebraic
properties.

To formalize this, assume that there exists a complex Lie group $\overline{G}$
which has a Lie algebra $\overline{\mathfrak{g}}$ containing $\mathfrak{g}$ as
a dense subalgebra. Without further comment we will denote the completion of
$\mathfrak{n}^{-}$ by $\overline{\mathfrak{n}^{-}}$ and so on. We will assume
that $\overline{G}$ has the following properties. First, the form
$\langle\cdot,\cdot\rangle$ extends continuously in each variable to
$\overline{\mathfrak{g}}$; secondly, the exponential map
$exp:\overline{\mathfrak{g}}\to\overline{G}$ is well-defined and holomorphic;
and thirdly, $\overline{G}$ has a Birkhoff stratification parameterized by the
elements of the Weyl group $W=N(H)/H$, compatible with, and having the same
basic properties as, the Birkhoff stratification for $G(A)$. We will not spell
out the third condition in full detail, but simply note $g$ will be in the top
stratum $\Sigma_{1}$ (a dense open set) iff the holomorphic functions
$\sigma_{j}(g)\neq 0$ iff $g$ has a unique LDU factorization, $g=ldu$,
$l\in\overline{N^{-}}$, $d=\prod_{j}\sigma_{j}(g)^{h_{j}}\in H$,
$u\in\overline{N^{+}}$, and analogous to the $G(A)$ case, $l$ and $u$ are
meromorphic in any given representation.

We now mimic (what is usually referred to as) the Kostant formulation of the
Toda lattice equations (see [7], although many other sources could be cited).
Define $\mathfrak{b}^{\pm}=\mathfrak{h}+\mathfrak{n}^{\pm}$. Since the form
$\langle\cdot,\cdot\rangle$ extends continuously to $\overline{\mathfrak{g}}$
there is an induced map with dense image,

$\overline{\mathfrak{b}^{-}}\to(\overline{\mathfrak{b}^{+}})^{*}:x\to\langle
x,\cdot\rangle$

This induces a Poisson structure on $\overline{\mathfrak{b}^{-}}$. Define
$\epsilon=\sum_{i=1}^{n}e_{i}$ and $M=\epsilon+\overline{\mathfrak{b}^{-}}$.
Via the identification

$\overline{\mathfrak{b}^{-}}\to M:x\to\epsilon+x$

there is an induced Poisson structure on $M$.

Suppose that $x(t)$ is a curve in $M$. The complex Kostant version of the Toda
equations is given by

(2.2) $\frac{d}{dt}x(t)=[proj_{\mathfrak{n}^{-}}(x(t)),x(t)]$

Since $x=\epsilon+x_{0-}$, $x_{0-}\in\overline{\mathfrak{b}^{-}}$, the right
hand side of (2.2) is in $\overline{\mathfrak{b}^{-}}$, the tangent space to
$M$. The basic facts about these equations for finite type Cartan matrices
carry over directly to the more general setting, namely

######  Proposition 1.

(a) The Toda equations are Hamiltonian for $M$. More precisely, define

$H:M\to\mathbb{C}:x\to\frac{1}{2}\langle x,x\rangle$

Then for any smooth function $f:M\to\mathbb{C}$,

$\frac{d}{dt}f(x(t))=\\{H,f\\}|_{x(t)}$

for any solution $x(t)$.

(b) Given an initial condition $x_{0}$ with $e^{x_{0}}\in\Sigma_{1}$,
$x(t)=Ad(l(t)^{-1})(x_{0})$ is a solution of (2.2), up to the time when
$e^{tx_{0}}$ exits $\Sigma_{1}$, where $e^{tx_{0}}$ has LDU factorization
$e^{tx_{0}}=l(t)d(t)u(t)$. $x(t)$ is a meromorphic function of $t$.

(c) Suppose that $x_{0}=\epsilon+\sum_{j=1}^{n}(a_{j}h_{j}+b_{j}f_{j})$ (note
in the $sl(n,\mathbb{C})$ case this is tridiagonal). Then the Toda equations
are equivalent to

$\frac{d}{dt}a_{j}=b_{j}\qquad\frac{d}{dt}b_{j}=\sum_{i=1}^{n}a_{i}A_{ij}b_{j},\qquad
j=1,...,n$

We do not need (a), although we will mention its relevance below.

###### Proof.

Define $x(t)=l(t)^{-1}x_{0}l(t)$ (in some matrix representation), as in (b).
Directly from this definition one calculates that

(2.3) $\frac{d}{dt}x=[l^{-1}\frac{d}{dt}l,x]$

Since $ldu=e^{tx_{0}}$, it follows that

$\frac{d}{dt}ldu+l\frac{d}{dt}{du}=e^{tx_{0}}x_{0}$

This implies

(2.4)
$l^{-1}\frac{d}{dt}l+\frac{d}{dt}(du)(du)^{-1}=l^{-1}e^{tx_{0}}x_{0}(du)^{-1}=l^{-1}x_{0}e^{tx_{0}}(du)^{-1}$

The last two expressions simplify to

(2.5) $=(du)x_{0}(du)^{-1}=l^{-1}x_{0}l$

and this is equal to $x(t)$. The equality of $x$ to the first expression in
2.4 implies that

$proj_{\mathfrak{n}^{-}}(x)=l^{-1}\frac{d}{dt}l$

Together with (2.3), this proves part (b).

Now consider part (c). It is convenient to write
$x_{0}=\epsilon+h\vec{a}(0)+f\vec{b}(0)$, where $h=(h_{1},...,h_{n})$,
$f=(f_{1},...,f_{n})$ and $a$ and $b$ are columns. We also need a bit more
structure. The Lie algebra $\mathfrak{g}$ is graded by height, where
$f_{j},h_{j},e_{j}$ have heights $-1,0,+1$ respectively. For example
$\epsilon$ is homogeneous of height $+1$, and $x_{0}$ is ‘tridiagonal’. Since
$l\in\overline{N}^{-}$ and $du\in\overline{B}^{+}$, the first formula for
$x(t)$ in (2.5) implies the homogeneous components of $x(t)$ have heights
$\geq-1$, and the second formula in (2.5) implies the homogeneous components
of $x(t)$ have heights $\leq+1$. Thus $x(t)$ is ‘tridiagonal’, and since
$x(t)\in M$, $x(t)=\epsilon+h\vec{a}(t)+f\vec{b}(t)$. The right hand side of
(2.2) is

(2.6)
$=[f\vec{b},\epsilon+h\vec{a}+f\vec{b}]=h\vec{b}+[f\vec{b},h\vec{a}]=h\vec{b}+f\vec{(\sum_{i}a_{i}A_{ij}b_{j})}$

where the second equality uses the Chevalley-Serre relations. This implies
(c).

∎

######  Corollary 1.

If the solutions in (c) are not meromorphic, then there does not exist a Lie
completion of $G(A)$ having the properties that we ascribed to $\overline{G}$
above.

The Toda equations in (c) were studied by Zou in the simplest rank 2 cases,
i.e. when the Cartan matrix $A=\left(\begin{matrix}2&-k\\\
-l&2\end{matrix}\right)$, see chapter 5 of [16], who was mainly interested in
integrability criteria from the Hamiltonian point of view.

######  Remark 2.

To put Zou’s equations (5.15)-(5.18) into our form, viewing $a$ and $b$ as
columns, make the linear change

$y=\left(\begin{matrix}\frac{2}{A}&\frac{2B}{B^{2}+C^{2}}\\\
0&\frac{2C}{B^{2}+C^{2}}\end{matrix}\right)a,\quad
x=-\left(\begin{matrix}\frac{2}{A^{2}}&0\\\
0&\frac{2}{B^{2}+C^{2}}\end{matrix}\right)b$

The ‘Cartan matrix’ is then $\left(\begin{matrix}2&\frac{2AB}{B^{2}+C^{2}}\\\
\frac{2B}{A}&2\end{matrix}\right)$. The off diagonal entries are not
necessarily nonnegative integers, so Zou is considering more general
equations.

In connection with this work, Zou found that the solutions are not generally
meromorphic, although in retrospect a definitive theorem seems to be lacking.
This appears to be a practical empirical criterion in particular cases, but to
our knowledge there are not any general theorems.

######  Remarks.

(a) The Virasoro algebra is a simple Lie algebra which is very similar to a
Kac-Moody Lie algebra and which does not have a global complexification (see
Proposition 3.3.2 of [11]).

(b) Physicists have been speculating for a long time that hyperbolic type Lie
algebras might be concretely related to supergravity. A list of hyperbolic
algebras, and some references to the physics literature, can be found in [2].

##  3\. An Example of the Verblunsky Correspondence

In the introduction we recalled one way to define the Verblunsky map
$Prob^{\prime}(S^{1})\to\prod_{n=1}^{\infty}\Delta$, using orthogonal
polynomials and the Szego recursion. This definition may seem to lack
motivation, and truthfully for us the definition seems a bit ad hoc. As it
turns out, there are other definitions, very different in appearance, such as
the Schur algorithm, see the remarkable theorems of Geronimus in section 3
[14]). Hence the definition is fully justified, even if it does not seem
entirely natural.

Here, in order to prove Theorems 1.1 and 1.2 of the introduction, we will
briefly recall from Section 2 of [8] our perspective on how to invert the
Verblunsky map, which emphasizes the connection with so called root subgroup
factorization for the loop group $LSU(1,1)$ (see [1] for how this can be
motivated).

Suppose that $\alpha\in\prod_{n=1}^{\infty}\Delta$, and suppose initially that
$\sum_{n=1}^{\infty}|\alpha_{n}|<\infty$. Then

$g_{2}(z):=\prod_{k=1}^{\stackrel{{\scriptstyle\infty}}{{\leftarrow}}}(1-|\alpha_{k}|^{2})^{-1/2}\left(\begin{matrix}1&\alpha_{k}^{*}z^{-k}\\\
\alpha_{k}z^{k}&1\end{matrix}\right)$

is a continuous function $g_{2}:S^{1}\to SU(1,1)$ (This is an example of a
(special) root subgroup factorization; the name arises because the factors
correspond to root subgroups of the loop group). This $SU(1,1)$ loop has the
very special form

(3.1) $g_{2}(z)=\left(\begin{matrix}d_{2}^{*}(z)&c_{2}^{*}(z)\\\
c_{2}(z)&d_{2}(z)\end{matrix}\right)$

where $c_{2}$ and $d_{2}$ are (boundary values of) holomorphic functions in
$\Delta$ satisfying $c_{2}(0)=0$ and $d_{2}(0)=1$. This loop also has a
triangular factorization (essentially a Riemann-Hilbert factorization) of the
especially simple form

(3.2) $g_{2}(z)=\left(\begin{matrix}1&\bf x^{*}(z)\\\
0&1\end{matrix}\right)\left(\begin{matrix}\mathbf{a}_{2}&0\\\
0&\mathbf{a}_{2}^{-1}\end{matrix}\right)\left(\begin{matrix}\alpha_{2}(z)&\beta_{2}(z)\\\
\gamma_{2}(z)&\delta_{2}(z)\end{matrix}\right)$

where the first matrix is holomorphic in $\Delta^{*}$ and $=1$ at $\infty$,
$\mathbf{a}_{2}$ is a positive constant, and the third matrix is a holomorphic
map $\overline{\Delta}\to SL(2,\mathbb{C}))$ and unipotent upper triangular at
$z=0$ (see [1]). It is important to note that this depends on the continuity
of $g_{2}$ (so that the Toeplitz operator corresponding to this symbol is
Fredholm), in addition to its special form, hence on the assumption that
$\alpha$ is absolutely summable.

######  Remark 3.

As an aside regarding completely integrable systems, in the Segal-Wilson
framework for KdV, loops into $SL(2,\mathbb{C})$ map to subspaces of a
Grassmannian, which in turn correspond to solutions of the KdV equation. In
the notation of [12], if $W$ denotes the subspace corresponding to the loop
$g_{2}(\alpha)$ then the (initial condition for) the Baker function is

$\psi_{W}(z)=1+z(x^{*}(z^{2}))=1+x_{1}^{*}z^{-1}+x_{2}^{*}z^{-3}+x_{3}^{*}z^{-5}+...$

Thus $x^{*}(z)$ is essentially the initial condition for the Baker function
(which itself has dependence on time). Whether root subgroup factorization (as
in [1]) has some use in the Segal-Wilson framework is unknown to us.

The following is explained in section 2.1 of [8]:

######  Proposition 2.

Via the Verblunsky correspondence, the measure $\mu\in Prob^{\prime}(S^{1})$
that corresponds to $\alpha\in\prod_{n=1}^{\infty}\Delta$, assuming
$(\alpha_{n})$ is absolutely summable, is given by

$d\mu=\frac{\prod_{k=1}^{\infty}(1-|\alpha_{k}|^{2})}{|\gamma_{2}+\delta_{2}|^{2}}\frac{d\theta}{2\pi}$

The absolute summability of $\alpha$ guarantees that the denominator for the
density of $d\mu$ is nonvanishing around the circle, hence the density is
continuous and invertible. Theorem 1.1 will show this fails for the marginally
non-summable sequence $\alpha_{n}=\frac{1}{1+n\beta}$. In general given
$\alpha\in\prod\Delta$, $\alpha$ and $\delta$ are holomorphic in $\Delta$, and
the sum $\alpha+\delta$ is nonvanishing in $\Delta$ (see the Caratheodory
function in the next paragraph). It is the boundary behavior which is in
question.

The measure $\mu$ is determined by a number of ancillary objects which are
explained in Section 2 of [14]. From our perspective on the inverse of the
Verblunsky map, these are expressed in the following way. Because $g_{2}$ has
values in $SU(1,1)$, the so called Schur function

$\mathbf{f}(z):=-\frac{c_{2}(z)}{zd_{2}(z)}=-\frac{\gamma_{2}(z)}{z\delta_{2}(z)}$

is a holomorphic map $\Delta\to\Delta$ (this follows from
$|d_{2}|^{2}-|c_{2}|^{2}=1$ on $S^{1}$, $c_{2}(0)=0$, and the maximum
principle). The Caratheodory function, implicitly defined by
$F(z)=\frac{1+z\mathbf{f}(z)}{1-z\mathbf{f}(z)}$ is a holomorphic function
which maps the disk to the right half plane. In our case, using the expression
for the Schur function,

(3.3) $F=\frac{\delta-\gamma}{\delta+\gamma}$

Note that on $S^{1}$

$F=\frac{(\delta_{2}-\gamma_{2})(\delta_{2}+\gamma_{2})^{*}}{(\delta_{2}+\gamma_{2})(\delta_{2}+\gamma_{2})^{*}}=\frac{\delta_{2}\delta_{2}^{*}-\gamma_{2}\gamma_{2}^{*}+Imaginary}{|\delta_{2}+\gamma_{2}|^{2}}$

implying

$Re(F)\frac{d\theta}{2\pi}=\frac{\prod_{k=1}^{\infty}(1-|\alpha_{k}|^{2})}{|\gamma_{2}+\delta_{2}|^{2}}\frac{d\theta}{2\pi}$

as in the Proposition.

The formula for the Caratheodory formula implies the following

######  Lemma 1.

$\delta=\frac{1}{2}(1+F)(\delta+\gamma)$

This implies that the ‘multiplicative regular part of $g_{2}$’, the
holomorphic function

$\left(\begin{matrix}\alpha_{2}&\beta_{2}\\\
\gamma_{2}&\delta_{2}\end{matrix}\right)=\left(\begin{matrix}\alpha_{2}&\beta_{2}\\\
\frac{1}{2}(1-F)(\delta_{2}+\gamma_{2})&\frac{1}{2}(1+F)(\delta_{2}+\gamma_{2})\end{matrix}\right)$

To clarify this lemma, to guarantee that $g_{2}$ has a triangular
factorization as in (3.2), we need to assume that $\alpha$ is absolutely
summable - for our example $\alpha_{n}=\frac{1}{1+n\beta}$, $g_{2}$ does not
have a triangular factorization. However the formula for $\delta$ is valid in
general, as an equality in $\Delta$.

###  3.1. Proof of Theorems 1.1 and 1.2

To explicitly understand the Verblunsky correspondence, one needs to calculate
$F$ and $\delta_{2}+\gamma_{2}$. In the remainder of this section we assume
$\beta>0$ and set $m=1/\beta$.

######  Theorem 3.1.

Suppose that $\alpha_{n}=\frac{1}{n\beta+1}$. Then

$\mu=\sqrt{\pi}2^{-m}\frac{\Gamma(m+1)}{\Gamma(m+\frac{1}{2})}(1-cos(\theta)^{m}\frac{d\theta}{2\pi}$

and hence

$Re(F)=\sqrt{\pi}2^{-m}\frac{\Gamma(m+1)}{\Gamma(m+\frac{1}{2})}(1-cos(\theta)^{m}$

###### Proof.

For the infinite product

$\left(\begin{matrix}\delta_{2}^{*}(z)&\gamma_{2}^{*}(z)\\\
\gamma_{2}(z)&\delta_{2}(z)\end{matrix}\right)=\prod_{k=1}^{\stackrel{{\scriptstyle\infty}}{{\leftarrow}}}\left(\begin{matrix}1&\alpha_{k}^{*}z^{-k}\\\
\alpha_{k}z^{k}&1\end{matrix}\right)$

it is straightforward to calculate that
$\gamma_{2}(z)=\sum_{n=1}^{\infty}\gamma_{2,n}z^{n}$ where

$\gamma_{2,n}=\sum\bar{\alpha}_{i_{1}}\alpha_{j_{1}}...\bar{\alpha}_{i_{r}}\alpha_{j_{r}}\bar{\alpha}_{i_{r+1}}$

the sum over multiindices satisfying

$0<i_{1}<j_{1}<..<j_{r}<i_{r+1},\quad\sum i_{*}-\sum j_{*}=n,$

and $\delta_{2}(z)=1+\sum_{n=1}^{\infty}\delta_{2,n}z^{n}$ where

$\delta_{2,n}=\sum\alpha_{i_{1}}\bar{\alpha}_{j_{1}}...\alpha_{i_{r}}\bar{\alpha}_{j_{r}}$

the sum over multiindices satisfying

$0<i_{1}<j_{1}<..<j_{r},\quad\sum(j_{*}-i_{*})=n.$

This implies that the $n$th coefficient of $\delta_{2}+\gamma_{2}$ is the left
hand side of the ‘super-telescoping formula’ (1.4). The right hand side of
(1.4) can be written as the ‘binomial coefficient’
$\left(\begin{matrix}m+n-1\\\ n\end{matrix}\right)$ (which is implicitly
defined in terms of $\Gamma$ functions). The negative binomial formula then
implies that

(3.4) $\delta_{2}+\gamma_{2}=\frac{1}{(1-z)^{m}}$

Proposition 2 now implies that

$\mu=\mathrm{Re}(F)\frac{d\theta}{2\pi}=\prod_{n=1}^{\infty}(1-|\alpha_{n}|^{2})|1-z|^{2m}\frac{d\theta}{2\pi}$
$=\prod_{n=1}^{\infty}(1-|\alpha_{n}|^{2})(2(1-cos(\theta))^{m}\frac{d\theta}{2\pi}$

Observe

$\int_{\theta=0}^{2\pi}(1-cos(\theta))^{m}d\theta=\int_{\theta=0}^{2\pi}(2sin^{2}(\frac{\theta}{2}))^{m}d\theta=2\int_{\phi=0}^{\pi}(2sin^{2}(\phi))^{m}d\phi$
$=2^{1+m}2\int_{\phi=0}^{\pi/2}sin(\phi)^{2m}d\phi=2^{1+m}\sqrt{\pi}\frac{\Gamma(m+\frac{1}{2})}{\Gamma(m+1)}$

Therefore, since $\mu$ is a probability measure

$\prod_{n=1}^{\infty}(1-|\alpha_{n}|^{2})=\sqrt{\pi}2^{-2m}\frac{\Gamma(m+1)}{\Gamma(m+\frac{1}{2})}$

hence

$\mu=\sqrt{\pi}2^{-2m}\frac{\Gamma(m+1)}{\Gamma(m+\frac{1}{2})}(2(1-cos(\theta))^{m}\frac{d\theta}{2\pi}$
$=\sqrt{\pi}2^{-m}\frac{\Gamma(m+1)}{\Gamma(m+\frac{1}{2})}(1-cos(\theta)^{m}\frac{d\theta}{2\pi}$

and this immediately implies the formula for $Re(F)$.

∎

To prove Theorem 1.2, we need to prove the following

######  Theorem 3.2.

In general

$F_{m}(z)=1+\sqrt{\pi}2^{1-m}\frac{\Gamma(m+1)}{\Gamma(m+\frac{1}{2})}\sum_{n=1}^{\infty}\sum_{k=0}^{\infty}\left(\begin{matrix}m\\\
n+2k\end{matrix}\right)\left(\begin{matrix}n+2k\\\
k\end{matrix}\right)(\frac{-1}{2})^{n+2k}z^{n}$

If $m$ is integral, then

$F_{m}(z)=1+2\frac{(m!)^{2}}{(2m-1)!!}\sum_{n=1}^{m}\sum_{k=0}^{[(m-n)/2]}\frac{1}{(m-n-2k)!k!(n+k)!}(\frac{-1}{2})^{n+2k}z^{n}$

######  Remark 4.

This $g_{2}$ is not continuous, reflecting the fact that $\alpha_{n}$ is not
summable. At least in simple cases, e.g. $m=1$, it can be checked that $g_{2}$
does not have a triangular factorization. Hence we cannot associate a solution
to KdV, as in Remark 3, to this kind of loop.

###### Proof.

This uses facts about Chebyshev expansions. Let
$\widetilde{T}_{0}=\frac{1}{2}T_{0}=\frac{1}{2}$ and $\widetilde{T}_{n}=T_{n}$
for $n>0$. If $r\geq 1$, then

$x^{r}=2^{1-r}\sum_{k=0}^{[r/2]}\left(\begin{matrix}r\\\
k\end{matrix}\right)\widetilde{T}_{r-2k}(x)$

The binomial formula implies

$(1-x)^{m}=2^{1-r}\sum_{r\geq 0}\left(\begin{matrix}m\\\
r\end{matrix}\right)(-1)^{r}x^{r}$ $=\sum_{r=0}^{m}\left(\begin{matrix}m\\\
r\end{matrix}\right)(-1)^{r}2^{1-r}\sum_{k=0}^{[r/2]}\left(\begin{matrix}r\\\
k\end{matrix}\right)\widetilde{T}_{r-2k}(x)$

Upon substituting $x=cos(\theta)$, this implies

$(1-cos(\theta))^{m}=2\sum_{r=0}^{m}\sum_{k=0}^{[r/2]}\left(\begin{matrix}m\\\
r\end{matrix}\right)\left(\begin{matrix}r\\\
k\end{matrix}\right)(\frac{-1}{2})^{r}\widetilde{T}_{r-2k}(cos(\theta))$

Since $T_{r-2k}(cos(\theta))=cos((r-2k)\theta)$, $r-2k\neq 0$,

$(1-cos(\theta))^{m}=2Re\sum_{r=0}^{m}\sum_{k=0}^{[r/2]^{\prime}}\left(\begin{matrix}m\\\
r\end{matrix}\right)\left(\begin{matrix}r\\\
k\end{matrix}\right)(\frac{-1}{2})^{r}z^{r-2k}$

where the prime indicates that we have to halve the sum when $r-2k=0$.
Therefore (incorporating the normalization factor for $F_{m}$)

$F_{m}(z)=1+\sqrt{\pi}2^{1-m}\frac{\Gamma(m+1)}{\Gamma(m+\frac{1}{2})}\sum\left(\begin{matrix}m\\\
r\end{matrix}\right)\left(\begin{matrix}r\\\
k\end{matrix}\right)(\frac{-1}{2})^{r}z^{r-2k}$

where the sum is over pairs $r,k\geq 0$, $r-2k>0$.

Now suppose we make the change of variables $n=r-2k$.

$F_{m}(z)=1+\sqrt{\pi}2^{1-m}\frac{\Gamma(m+1)}{\Gamma(m+\frac{1}{2})}\sum_{n=1}^{\infty}\sum_{k=0}^{\infty}\left(\begin{matrix}m\\\
n+2k\end{matrix}\right)\left(\begin{matrix}n+2k\\\
k\end{matrix}\right)(\frac{-1}{2})^{n+2k}z^{n}$

∎

######  Examples.

In case someone recognizes something!:

$F_{1}(z)=1-z$ $F_{2}(z)1-\frac{4}{3}z+\frac{1}{3}z^{2}$
$F_{3}=1-\frac{15}{10}z+\frac{6}{10}z^{2}-\frac{1}{10}z^{3}$
$F_{4}=1-\frac{56}{35}z+\frac{48}{35}z^{2}-\frac{8}{35}z^{3}+\frac{1}{36}z^{4}$

###  3.2. Moments of $\mu$

To compute the moments (Fourier coefficients) of
$d\mu=\frac{1}{\mathfrak{Z}}(1-\cos(\theta))^{m}d\theta/2\pi$ write

$d\mu=\frac{1}{\mathfrak{Z}}2^{-m}\left(2-e^{i\theta}-e^{-i\theta}\right)^{m}\frac{d\theta}{2\pi}$

Then

(3.5) $c_{n}=\frac{1}{\mathfrak{Z}}2^{-m}\int
e^{-in\theta}\left(2-e^{i\theta}-e^{-i\theta}\right)^{m}\frac{d\theta}{2\pi}$

Expand

$\left(2-e^{i\theta}-e^{-i\theta}\right)^{m}=\sum_{k=0}^{m}{m\choose
k}2^{m-k}(-1)^{k}\left(e^{i\theta}+e^{-i\theta}\right)^{k}=\sum_{k=0}^{m}\sum_{l=0}^{k}{m\choose
k}{k\choose l}2^{m-k}(-1)^{k}e^{i(k-2l)\theta}$

and rewrite the sum in terms of $r=k-2l$

(3.6)
$=\left(2-e^{i\theta}-e^{-i\theta}\right)^{m}=\sum_{r=-m}^{m}\sum_{l=0}^{\lfloor\frac{m-r}{2}\rfloor}{m\choose
r+2l}{r+2l\choose l}2^{m-r-2l}(-1)^{r+2l}e^{ir\theta}$

This implies

######  Theorem 3.3.

The Fourier coefficients of the measure
$\frac{1}{\mathfrak{Z}}(1-\cos(\theta))^{m}d\theta/2\pi$ on $S^{1}$ are given
by

(3.7)
$c_{n}=\sqrt{\pi}2^{-m}\frac{\Gamma(m+1)}{\Gamma(m+\frac{1}{2})}\sum_{l=0}^{\lfloor\frac{m-n}{2}\rfloor}{m\choose
n+2l}{n+2l\choose l}(\frac{-1}{2})^{n+2l}$

$n=1,2,...$

In comparison with the other explicit examples in [13], we do lack transparent
expressions for the orthogonal polynomials for $\mu$.

## References

  * [1] E. Basor and D. Pickrell, Loops in $SL(2,\mathbb{C})$ and root subgroup factorization, Random Matrices and Applications 7, no. 3 (2018).
  * [2] L. Carbone, S. Chung, L. Cobbs, R. Mcrae, D. Namdi, Y. Naqvi, and D. Penta, Classification of hyperbolic Dynkin diagrams, root lengths, and Weyl group orbits, J. Phys. A 43, no. 15 (2010).
  * [3] L. Carbone, E. Jurisich, S. Murray, Constructing a Lie group analog for the Monster Lie algebra Lett. Math. Phys. 112 (2022), no. 3, Paper No. 43 (2022).
  * [4] P. Goddard, The work of Richard Borcherds, Proceedings of the International Congress of Mathematicians, Berlin, (1998) 99-108.
  * [5] V. Kac, Constructing groups from infinite dimensional Lie algebras, Infinite Dimensional Groups with Applications, edited by V Kac, MSRI publication, Springer-Verlag (1985) 167-216.
  * [6] V. Kac and D. Peterson, Defining relations of certain infinite dimensional groups, Elie Cartan et les mathematiques d’aujourd’hui - Lyon, 25-29 juin 1984, Asterisque, no. S131 (1985), 44 p.
  * [7] B. Kostant, The solution to a generalized Toda lattice and representation theory. Adv. in Math. 34 (1979), no. 3, 195-338.
  * [8] M. J. Latifi Jebelli and D. Pickrell, Exponential of the $S^{1}$ Trace of the Free Field and Verblunsky Coefficients, Rocky Mountain J. Math. 52 (2022), no. 3, 899-924.
  * [9] M.J. Latifi Jebelli, Lattice models and super telescoping formula, arXiv:2108.07770, (2021)
  * [10] I. Nenciu, Lax Pairs for the Ablowitz-Ladik System via Orthogonal Polynomials on the Unit Circle, dissertation, Cal Tech (2005)
  * [11] A. Pressley and G. Segal, Loop groups, Oxford Mathematical Monographs, Oxford Science Publications, Oxford University Press, New York (1986).
  * [12] G. Segal and G. Wilson, Loop groups and equations of KdV type. Inst. Hautes Etudes Sci. Publ. Math. No. 61 (1985), 5-65.
  * [13] B. Simon, Orthogonal Polynomials on the Unit Circle, Part I Classical Theory, Part 2. Spectral theory. AMS Colloquium Series, Amer. Math. Soc., Providence, RI (2005)
  * [14] B. Simon, OPUC on one foot, Bull. Amer.Math. Soc., Vol. 42, No 4 (2005) 431-460.
  * [15] B. Yang, Algebraic Aspects of the Dispersionless Limit of the Discrete Nonlinear Schrodinger Equation, dissertation, University of Arizona (2013).
  * [16] M. Zou, The Geometry of Two Degree of Freedom Integrable Hamiltonian Systems, dissertation, University of Arizona (1992).


# Two-stage assembly of patchy ellipses: From bent-core particles to liquid
crystal analogs

Anuj Kumar Singh1 Arunkumar Bupathy1,2 Jenis Thongam3 Emanuela Bianchi3,4
Gerhard Kahl3 Varsha Banerjee1 1Department of Physics, Indian Institute of
Technology Delhi, New Delhi 110016, India.  2School of Chemistry, University
of Birmingham, Birmingham B15 2TT, UK.  3Institut für Theoretische Physik, TU
Wien, Wiedner Hauptstraße 8-10, A-1040 Wien, Austria  4CNR-ISC, Uos Sapienza,
Piazzale A. Moro 2, 00185 Roma, Italy.

###### Abstract

We investigate the two-dimensional behavior of colloidal patchy ellipsoids
specifically designed to follow a two-step assembly process from the monomer
state to mesoscopic liquid-crystal phases, via the formation of so-called
bent-core units at the intermediate stage. Our model comprises a binary
mixture of ellipses interacting via the Gay-Berne potential and decorated by
surface patches, with the binary components being mirror-image variants of
each other - referred to as left-handed and right-handed ellipses according to
the position of their patches. The surface patches are designed so as in the
first stage of the assembly the monomers form bent-cores units, i.e. V-shaped
dimers with a specific bent angle. The Gay-Berne interactions, which act
between the ellipses, drive the dimers to subsequently form the characteristic
phase observed in bent-core liquid crystals. We numerically investigate – by
means of both Molecular Dynamics and Monte Carlo simulations – the described
two-step process: we first optimize a target bent-core unit and we then fully
characterize its state diagram in temperature and density, defining the
regions where the different liquid crystalline phases dominate.

##  I Introduction

During the past years bent-core liquid crystals (BCLCs) have attracted a
rapidly increasing interest both in experimental and theoretical
investigations Lubensky and Radzihovsky (2002); Ros  _et al._ (2005); Takezoe
and Takanishi (2006); Jákli  _et al._ (2018). The constituent units of these
systems are bent-core (BC) colloidal entities (typically ranging in size from
nm to $\mu$m) that have a (sharp or smooth) bent shape. This characteristic
shape can be realized in different ways, for instance by merging two
ellipsoids (or sphero-cylinders) via their tips at a (possibly flexible) angle
Johnston  _et al._ (2002); Bates and Luckhurst (2005); Bates (2006) or by
creating bead-resolved, bent strings of spherical particles (often referred to
as “banana-shaped” colloids) Dewar and Camp (2004, 2005); Chiappini  _et al._
(2019). The reason for the considerable interest in such BC particles is the
large variety of (ordered) phases that these entities can form: fascinating
chiral nematic, smectic, splay, and twist-bend phases (despite the achiral
nature of the BC constituents) are just a few examples.

Here we propose a model system where the BC units result from the spontaneous
assembly of suitably designed particles characterized by shape and bond
anisotropy, namely we design anisotropic, patchy colloids such that they form
V-shaped dimers. To introduce shape anisotropy, we choose – in our two-
dimensional setup – elliptic particles, which interact via the well-
established Gay-Berne (GB) potential Berardi  _et al._ (1995). The relevant
parameters that characterize this interaction are the aspect ratio of the
ellipses and their interaction range. The other ingredient to our model is the
patchiness of the colloids Pawar and Kretzschmar (2010); Bianchi  _et al._
(2011); Ravaine and Duguet (2017); Bianchi  _et al._ (2017): patchy colloids
have surfaces decorated with well defined regions (both in terms of position
and of spatial extent) that support a different interaction behaviour with
respect to the bare surface. Thus, depending on their mutual interaction
properties these regions can either foster or suppress patch-patch directional
bonding. Patchiness has become a widely used standard tool (both in
experiments and in theoretical/simulation-based investigations) to trigger in
bottom-up processes the self-assembly of decorated particles into highly
complex bonding patterns, ranging from liquid-crystalline mesophases, over
complex ordered structures and quasi-crystalline lattices to structural
glasses Smallenburg  _et al._ (2013); Chen  _et al._ (2011); Wang  _et al._
(2012); Morphew  _et al._ (2018); He  _et al._ (2020); Posnjak  _et al._
(2024); Liu  _et al._ (2024); Karner and Bianchi (2024).

In the present contribution we describe in detail the process that starts from
the formation of BCs via the self-assembly of elliptic patchy colloids and
ends up in the self-organization of these entities into BCLCs. Our design of
the elliptic patchy particles envisages a binary mixture of mirror-symmetric
units: while both species have the same elliptic shape, the two types of
particles are mirror symmetric with respect to their patch decoration (leading
thus to left- and right-handed particles). In fact, the two patches (per
particle) are located close to the tip of the colloid. The form of the patch-
patch interaction is based on the work of Russo et al. Russo  _et al._ (2009),
while the choice of the position, interaction range and interaction strength
of the patches influences the features of the BC units. In a first step we
demonstrate that our particles are indeed able to self-assemble via the
related patch-patch bonding mechanism into stable BC entities; we demonstrate
that our model is characterized by a high versatility, as the bonding angle
can be tailored via the patch positions and/or the patch-patch interaction
strength. The impact of all these parameters on the actual shape, stability
and yield of these units is studied in detail. In a second step we consider an
ensemble of such BCs and study their self-assembly into mesophases. Focusing
for the moment on a particular set of model parameters we are able to trace
out in a qualitative manner a diagram of states (in the temperature-density
plane), which shows the emergence of isotropic and smectic phases. All
investigations have been carried out via extensive computer-simulations: we
have applied both molecular dynamics (MD) and – in a complementary manner –
Monte Carlo (MC) simulations, using different types of thermodynamic
ensembles. As the proposed BC model can be easily extended to support BCs with
different sets of features, it reveals itself as a powerful tool to gain a
deeper insight into this fascinating class of functional materials with
potential applications in optics, photonics, and sensing. The rest of our
paper is organized as follows. Sec. II provides the theoretical framework
which includes the model and the methodology to characterize the bent-cores.
The computational details of the MD and MC procedures are provided in Sec.
III. The results and discussion follow in Sec. IV. Finally, Sec. V provides
the summary and conclusion.

##  II Theoretical framework

###  II.1 Model

We want to design a simple model in which colloidal monomers form a bent-core
(BC) entity first, which then self-assemble into mesoscopic liquid crystalline
phases. Such a targeted design process would allow for exploring a large
number of phases which emerge by changing the shape as well as the surface
chemistry of the colloidal monomers.

To simulate such a system we use a simple Gay-Berne (GB) potential Berardi
_et al._ (1995) for the elliptic colloidal monomers which are decorated on
their surface by suitably placed patches (index ‘P’). The form of the patch
potential is adopted from the work of Russo et al. Russo  _et al._ (2009).
Thus the interaction potential between a pair of particles $i$ and $j$ is
given by (where we have suppressed for the moment the arguments)

$V_{ij}=V_{\rm GB}+V_{\rm P},$ (1)

where $V_{\rm GB}$ is the Gay-Berne potential acting between the centers of
mass of the two particles, and $V_{\rm P}$ is the interaction potential due to
all pairs of patches located on the surfaces of particles $i$ and $j$.

The GB potential $V_{\rm GB}$ characterizes the shape-anisotropy of the
elliptical particles and is defined as follows:

$V_{GB}(\hat{u}_{i},\hat{u}_{j},\boldsymbol{r}_{ij})=4\epsilon(\hat{u}_{i},\hat{u}_{j},\hat{r}_{ij})\left[\left(\frac{\sigma_{0}}{r_{ij}-\sigma(\hat{u}_{i},\hat{u}_{j},\hat{r}_{ij})+\sigma_{0}}\right)^{12}\right.\\\
\left.-\left(\frac{\sigma_{0}}{r_{ij}-\sigma(\hat{u}_{i},\hat{u}_{j},\hat{r}_{ij})+\sigma_{0}}\right)^{6}\right],$
(2)

where $\sigma_{0}$ is the unit distance of the interactions.
$\boldsymbol{r}_{ij}$ is the vector connecting the centers of mass of the
particles and $\hat{u}_{i}$ and $\hat{u}_{j}$ are the unit vectors
representing the orientations of the particles as shown in Fig. 1(a).
$\epsilon=\epsilon_{0}\epsilon_{1}^{\mu}\epsilon_{2}^{\nu}$ is the orientation
dependent interaction strength, where $\epsilon_{0}$ represents the energy
unit of the interactions; the parameters $\mu$ and $\nu$ modify the well depth
of the potential; $\sigma$ is the anisotropic contact distance. Further
$\epsilon_{1}$, $\epsilon_{2}$ and $\sigma$ are defined as:

$\displaystyle\epsilon_{1}(\hat{u}_{i},\hat{u}_{j},\hat{r}_{ij})$
$\displaystyle=1-\frac{\chi^{\prime}}{2}\left\\{\frac{(\hat{r}_{ij}.\hat{u}_{i}+\hat{r}_{ij}.\hat{u}_{j})^{2}}{1+\chi^{\prime}\hat{u}_{i}.\hat{u}_{j}}+\frac{(\hat{r}_{ij}.\hat{u}_{i}-\hat{r}_{ij}.\hat{u}_{j})^{2}}{1-\chi^{\prime}\hat{u}_{i}.\hat{u}_{j}}\right\\},$
(3) $\displaystyle\epsilon_{2}(\hat{u}_{i},\hat{u}_{j})$
$\displaystyle=\left[1-\chi^{2}(\hat{u}_{i}.\hat{u}_{j})^{2}\right]^{-\frac{1}{2}},$
(4) $\displaystyle\sigma(\hat{u}_{i},\hat{u}_{j},\hat{r}_{ij})$
$\displaystyle=\sigma_{0}\left[1-\frac{\chi}{2}\left\\{\frac{(\hat{r}_{ij}.\hat{u}_{i}+\hat{r}_{ij}.\hat{u}_{j})^{2}}{1+\chi\hat{u}_{i}.\hat{u}_{j}}+\frac{(\hat{r}_{ij}.\hat{u}_{i}-\hat{r}_{ij}.\hat{u}_{j})^{2}}{1-\chi\hat{u}_{i}.\hat{u}_{j}}\right\\}\right]^{-\frac{1}{2}}.$
(5)

The energy anisotropy parameter $\chi^{\prime}$ and the shape anisotropy
parameter $\chi$ are defined as:

$\displaystyle\chi^{\prime}$ $\displaystyle=(\kappa^{\prime
1/\mu}-1)/(\kappa^{\prime 1/\mu}+1),$ (6) $\displaystyle\chi$
$\displaystyle=(\kappa^{2}-1)/(\kappa^{2}+1),$ (7)

where $\kappa^{\prime}=\epsilon_{s}/\epsilon_{e}$ is the energy anisotropy
with $\epsilon_{s}$ and $\epsilon_{e}$ being the well depths of the side-to-
side and end-to-end configurations. $\kappa=\sigma_{\parallel}/\sigma_{\perp}$
is the aspect ratio of the elliptic particles with $\sigma_{\parallel}$ and
$\sigma_{\perp}$ being the major and minor axes of the particles. The GB
potential therefore contains four key parameters to describe interacting
elliptic particles: $\kappa$, $\kappa^{\prime}$, $\mu$, and $\nu$. There exist
a number of GB homologues in the context of nematic LCs depending on the
values chosen for the four parameters Gay and Berne (1981); Berardi and RJ
(1993); Bates and Luckhurst (1999); Zannoni (2001). We will discuss the
appropriate choice for obtaining the BCs and their LC-like self-assembly
shortly.

The patch-patch interactions $V_{\rm P}$ are described by a sum of inverted
Gaussian potentials Russo  _et al._ (2009):

$V_{\rm P}(d)=-\sum_{k\in P_{i}}\sum_{l\in
P_{j}}\epsilon_{kl}\exp\left[-\frac{1}{2}\left(\frac{d_{kl}}{\alpha}\right)^{2}\right],$
(8)

where $\epsilon_{kl}$ is the strength of the interaction between patches $k$
and $l$, $P_{i}$ and $P_{j}$ are the set of patches on the particles $i$ and
$j$, respectively, and $d_{kl}$ is the distance between the centers of patches
$k$ and $l$, which decorate the surface of the elliptical particles, as shown
in the schematic of Fig. 1(a). The width of the inverted Gaussian, $\alpha$,
defines the range of the patch-patch interactions. In addition to being able
to compute gradients required for molecular dynamics (MD) simulations, this
choice of the patch potential also allows for oblique patch placements, i.e.,
at an angle from the tip of the ellipses. Such placements promote the assembly
of stable BCs.

Let us consider the schematic in Fig. 1(a) that shows ellipsoidal particles
with two patches placed obliquely; the reason for a two-patch model will be
clear soon. In our model there are two kinds of patch placements on the
ellipsoids: one to the “left” (red particles) and the other to the “right”
(green particles). The green particles have patches $G_{1}$ and $G_{2}$ at
angles $-\theta_{1}$ and $-\theta_{2}$ from the major axis. Similarly, the red
particles have patches $R_{1}$ and $R_{2}$ at angles $\theta_{1}$ and
$\theta_{2}$ from the major axis. The interaction strengths
$\epsilon_{G_{1}R_{1}}$ and $\epsilon_{G_{2}R_{2}}$ are set equal to
$\epsilon_{p}$, while all other interactions are set to 0: these include
notably $\epsilon_{G_{1}G_{1}}$, $\epsilon_{G_{1}G_{2}}$,
$\epsilon_{R_{1}R_{1}}$, $\epsilon_{R_{1}R_{2}}$, $\epsilon_{G_{1}R_{2}}$ and
$\epsilon_{G_{2}R_{1}}$. We consider a pair of left-handed and right-handed
particles to be a BC unit if the patch interaction strength $V_{p}$ is
sufficiently large to prevent dissociation of this arrangement by thermal
fluctuations. In our simulations, we find that a choice of $V_{p}<-3k_{\rm
B}T$ ($T$ being the temperature and $k_{\rm B}$ being the Boltzmann constant)
ensures that the BCs once formed remain bonded.

Fig. 1(b) shows a schematic of a BC colloid corresponding to the above model
description. This entity is characterized by the director ${\hat{n}}$, the
polarization $\hat{p}$, and the bend angle $\gamma$ (as indicated). In our
investigations we have found that single patch particles are characterized by
a broad distribution of the bend angles, as shown in Fig. 1(c). Reducing the
width of the angle distribution would require lowering the value of $\alpha$,
which, as we shall see later, negatively affects the kinetics of BC self-
assembly. Further, multiple single-patch ellipses can adhere in clusters (or
asters). Having two patches on the other hand results in smaller magnitudes of
the bend angle and lesser flexibility. Thus the two-patch ellipsoids represent
an optimal choice to obtain BCs at a high yield. The model parameters
$\theta_{i}$, $\epsilon_{p}$, and $\alpha$ have a significant impact on the BC
yield and phases. We will study their influence in the rest of the paper.

###  II.2 Methods

####  II.2.1 Bent-core characteristics

Let us refer to the schematic in Fig. 1 that depicts a pair of left-handed and
right-handed patchy ellipsoids. Let their center of mass, after the formation
of the BC, be represented by $\boldsymbol{r}_{\rm cm}^{\rm l}$ and
$\boldsymbol{r}_{\rm cm}^{\rm r}$, respectively. The center of mass of the BC
unit is given by the vector average:

$\boldsymbol{R}_{\rm cm}=\frac{\boldsymbol{r}_{\rm cm}^{\rm
l}+\boldsymbol{r}_{\rm cm}^{\rm r}}{2}.$ (9)

Further, the local director $\hat{n}$ and the polarization $\hat{p}$ can be
expressed in terms of the orientations $\boldsymbol{u}_{l}$ and
$\boldsymbol{u}_{r}$ as follows:

$\hat{n}=\frac{\boldsymbol{u}_{\rm l}-\boldsymbol{u}_{\rm
r}}{|\boldsymbol{u}_{\rm l}-\boldsymbol{u}_{\rm
r}|},\quad\hat{p}=\frac{\boldsymbol{u}_{\rm l}+\boldsymbol{u}_{\rm
r}}{|\boldsymbol{u}_{\rm l}+\boldsymbol{u}_{\rm r}|},$ (10)

Similarly, the bend angle $\gamma$ can also be defined in terms of the
orientation vectors as follows:

$\gamma=\cos^{-1}\left(\frac{\boldsymbol{u}_{\rm l}\cdot\boldsymbol{u}_{\rm
r}}{|\boldsymbol{u}_{\rm l}||\boldsymbol{u}_{\rm r}|}\right).$ (11)

As discussed in the context of Fig. 1(b), the BC unit is characterized by
$\hat{n}$, $\hat{p}$, and $\gamma$. We will use Eqs. (10) and (11) for the
evaluation of these quantities in all our simulations.

####  II.2.2 Q-Tensor

To characterize the global orientational order of the mesoscopic phases formed
by the BC units, we compute the $Q$-tensor as an appropriate order parameter
for quantifying both the nematic and the polar order De Gennes and Prost
(1993); Wei and Patey (1992); Allen and Tildesley (2017). Its components are
defined as follows:

$Q_{\alpha\beta}=\frac{1}{N_{\rm b}}\sum_{i=1}^{N_{\rm
b}}\frac{1}{2}\left(3n_{i,\alpha}n_{i,\beta}-\delta_{\alpha\beta}\right),$
(12)

where $N_{\rm b}$ is the number of BCs in the ensemble, $\alpha$,
$\beta\in\\{x,y\\}$ represent the Cartesian components of $\hat{n}_{i}$
corresponding to the $i^{th}$ BC. The largest eigenvalue of $Q$ gives the
nematic order parameter $\mathcal{S}$ and the corresponding eigenvector
represents the global director ${\bf n}$. The nematic order parameter can also
be interpreted as $\mathcal{S}=\overline{P_{2}(\hat{n}_{i}\cdot\mathbf{n})}$,
where $P_{2}(x)=(3x^{2}-1)/2$ is the Legendre polynomial of degree two and the
over-bar indicates an average over all the molecules. The isotropic phase
corresponds to ${\cal S}=0$ and a fully aligned nematic phase has ${\cal
S}=1$. The polar order can be computed as
$\mathcal{P}=\left|\overline{\hat{p}_{i}\cdot\mathbf{p}}\right|$, where
$\mathbf{p}$ is normal to $\mathbf{n}$. A polarization $\mathcal{P}=1$ implies
that the polar vectors $\hat{p}_{i}$ of the BCs are all pointing in the same
direction as $\mathbf{p}$. In this case, aggregates of BCs exhibit nematic as
well as polar order. $\mathcal{P}=0$ indicates random or staggered (e.g.,
anti-ferroelectric) orientation of $\hat{p}_{i}$ of the BCs. Note that in this
case, the aggregate could still exhibit nematic order.

####  II.2.3 Radial distribution function

The radial distribution function $g(r)$ (RDF) which measures correlations
between pairs of particles can provide information about the presence of
density modulations in the system and is a useful tool for identifying the
structure of the aggregates or domains. This function is defined as Frenkel
and Smit (2023); Billeter and Pelcovits (1998)

$g(r)=\frac{1}{\rho_{\rm b}N_{\rm b}}\left\langle\sum_{\stackrel{{\scriptstyle
i,j}}{{i\neq j}}}^{N_{\rm b}}\delta(r-r_{ij})\right\rangle,$ (13)

where $r_{ij}=|\boldsymbol{r}_{ij}|$ is the separation distance between the
centers of mass of BCs $i$ and $j$, $N_{\rm b}$ is the number of BCs, and
$\rho_{\rm b}=N_{\rm b}/L^{2}$ is their number density. In practice, the
$\delta$-function is replaced by a counting function which is normalized
appropriately by the area of the finite but small region in which the counting
is done. The division by $\rho_{\rm b}$ ensures that $g(r)=1$ when the
particles are completely uncorrelated, as this is, for instance realized in an
ideal gas. Sharp peaks indicate the presence of structure and oscillatory
behavior marks the presence of translational order in the system. We can also
define a characteristic lengthscale (or correlation length) over which the
particle positions are correlated. Typically, it is the distance $r$ over
which the peaks in $g(r)$ decay to – say – 0.2 times the maximum value of the
RDF Puri (2004).

In order to identify layered phases such as the smectic one, we resolve the
RDF along axes in the local frame of reference of the BC units and average the
results. The RDF resolved parallel to the local director $\hat{n}_{i}$ of the
BCs, $g_{\parallel}(r_{\parallel})$, is defined as follows Weis and Levesque
(1993); Cañeda-Guzmán  _et al._ (2014):

$g_{\parallel}(r_{\parallel})=\frac{1}{\rho_{b}N_{b}}\left\langle\sum_{\stackrel{{\scriptstyle
i,j}}{{i\neq
j}}}^{N_{b}}\frac{\delta(r_{\parallel}-r_{ij,\parallel})\Theta(h/2-r_{ij,\perp})}{h}\right\rangle,$
(14)

where $r_{ij,\parallel}=|\boldsymbol{r_{ij}}\cdot\hat{n}_{i}|$ and
$r_{ij,\perp}=|\boldsymbol{r_{ij}}\times\hat{n}_{i}|$ are the parallel and
perpendicular components of the separation $r_{ij}$ with respect to
$\hat{n}_{i}$; further, $h$ is the average height of the BCs (i.e., average
extension of the BCs in the direction perpendicular to $\hat{n}_{i}$, see Fig.
1b) and $\Theta(x)$ is the Heaviside step function, which ensures that only
BCs within the same layer as the reference BC (with index $i$) are included in
the computation. Similarly, $g_{\perp}(r_{\perp})$ which resolves RDF in the
direction perpendicular to the local director is defined as Weis and Levesque
(1993); Cañeda-Guzmán  _et al._ (2014)

$g_{\perp}(r_{\perp})=\frac{1}{\rho_{\rm b}N_{\rm
b}}\left\langle\sum_{\stackrel{{\scriptstyle i,j}}{{i\neq j}}}^{N_{\rm
b}}\frac{\delta(r_{\perp}-r_{ij,\perp})\Theta(w/2-r_{ij,\parallel})}{w}\right\rangle,$
(15)

where $w$ is the average width of the BCs (see Fig. 1b).

##  III Computational Details

Our starting point to obtain BC colloids and their phases is a collection of
an equal number of right-handed and left-handed two-patch ellipsoids, and
specified patch-patch interactions. We study the self-assembly using extensive
molecular dynamics (MD) and Monte Carlo (MC) simulations. As discussed in
detail in Sec. II, the ellipsoids interact via the GB potential which is
characterized by four essential parameters: $\kappa$, $\kappa^{\prime}$, $\mu$
and $\nu$. There are a large variety of GB homologues which differ from each
other in terms of the values chosen for the four parameters. The primary focus
in selecting their values has been (i) to obtain the nematic and smectic LC
phases and (ii) to obtain convergence with experimental data Gay and Berne
(1981); Berardi and RJ (1993); Bates and Luckhurst (1999); Zannoni (2001). The
choice of $\kappa=3.0$ is elementary as for real LC systems, the length-to-
breadth ratio of the constituent molecules is typically equal to or greater
than 3:1. A phase diagram of the GB model in the $\kappa^{\prime}-T$ indicates
that the ordered phase as $\kappa^{\prime}\rightarrow 1$ is primarily nematic
Birdi  _et al._ (2022). A value of $\kappa^{\prime}>1$ introduces energy
anisotropy, favoring side-to-side alignment of ellipses. These values promote
nematic as well as smectic order characteristics of LCs. The parameters $\mu$
and $\nu$ modify the well depths of the potential, and hence their impact on
the nematicity and smecticity is very subtle Gay and Berne (1981); Berardi and
RJ (1993); Bates and Luckhurst (1999); Zannoni (2001). A popular choice of
parameters is $\kappa=3$, $\kappa^{\prime}=5$, $\mu=1$ and $\nu=3$ due to
Berardi et al. Berardi and RJ (1993). They ensure the formation of LC phases
and agreement with experimental data in many real systems Berardi and RJ
(1993). In the present work, our focus is to first create BC units, and then
study their assembly in ordered phases. Consequently, we choose
$\kappa^{\prime}=1$ to promote patch-patch aggregation for the formation of
BCs in the first stage, which then assemble to form ordered phases due to the
GB interactions in the second stage.

###  III.1 Details of MD simulations

MD simulations are performed in NVT and NPT ensembles using the LAMMPS
software LAM ; Plimpton (1995), employing the Nosé-Hoover thermostat and
barostat for temperature and pressure coupling, as needed Frenkel and Smit
(2023); Hoover (1985, 1999). Periodic boundaries with minimum image convention
have been used for the simulations. All calculations were done in reduced LJ
units: $T^{*}=k_{B}T/\epsilon_{0}$, $\Delta t^{*}=\Delta
t/\sqrt{m\sigma_{0}^{2}/\epsilon_{0}}$, where $k_{B}$ is the Boltzmann
constant. The star is dropped in subsequent discussions. The equations of
motion are integrated using the velocity-Verlet algorithm with simulation time
step $\Delta t=0.001$ Andersen (1983).

We initialize the system with $N=4000$ patchy ellipsoids in a square box of
length $L$ depending on the desired density
$\rho=(\pi\sigma_{\parallel}\sigma_{\perp})N/L^{2}$, where
$\sigma_{\parallel}$ and $\sigma_{\perp}$ are the major and minor axes of the
ellipsoids. For all studied densities, we initialize the particles on a
regular lattice and allow the system to equilibrate at $T=3$ for $10^{5}$ MD
steps. The resulting homogeneous state is then quenched to the desired
temperature. Simulations were performed up to $10^{7}$ steps and relevant
quantities measured at intervals of $5\times 10^{3}$ steps. All presented data
have been averaged over 10 independent runs.

###  III.2 Details of MC simulations

In a complementary fashion we have used Monte Carlo (MC) Frenkel and Smit
(2023); Allen and Tildesley (2017) simulations to investigate the structural
properties of our BC entities. To be more specific we have used so-called
“virtual-move” MC (VMMC) simulations Whitelam and Geissler (2007); Whitelam
_et al._ (2009); Hedges , an algorithm that proposes simultaneous
translational and rotational moves of clusters of particles according to
gradients of interaction energies of the particles. In short, a seed particle
is subjected to a “virtual” trial move; based on the emerging bond energies
before and after the trial move, neighbouring particles are added to a virtual
cluster. This process is iterated for each particle in the cluster and the
positions of all particles involved are updated simultaneously (for more
details we refer to Refs. Whitelam and Geissler (2007); Whitelam  _et al._
(2009)).

The VMMC simulations were performed in the NVT ensemble, applying standard
periodic boundary conditions. Due to the fact that the degree of
parallelization of the VMMC algorithm is strongly limited, we had to restrict
ourselves to considerably smaller ensembles (as compared to the MD
simulations): a typical ensemble size of $N=900$ particles was used. Starting
in a square box (whose size is imposed by the desired density $\rho$) half of
the particles are left- and right-handed, respectively. At a given temperature
the system is equilibrated during 4 $\times$ 107 VMMC steps, while the
production part of the simulations extended over 2 $\times$ 107 VMMC steps.
For each state point we considered ten independent MC runs, each of them
starting from different initial conditions. For the highest densities we
started from ordered arrays of parallel unbounded ellipses, while for the
lowest densities we started from random configurations. All quantities were
averaged over 4 $\times$ 103 independent configurations per single run, adding
to 4 $\times$ 104 per state point.

##  IV Results and Discussion

###  IV.1 First step assembly stage

Before studying the mesoscopic phases of the system, it is useful to
understand the effect of the model parameters on the formation of BCs. The
presented patchy colloidal analog of BCLCs has three model parameters and two
system parameters to control the properties of the BCs: patch interaction
strength $\epsilon_{p}$ (defined in units of $\epsilon_{0}$), patch
interaction range $\alpha$, patch positions $\theta_{i}$, density $\rho$, and
temperature $T$. We perform simulations in the $NVT$ ensemble with density
$\rho=0.4$ and temperature $T=1.2$, corresponding to an equilibrium isotropic
fluid state of the underlying Gay-Berne particles. Unless stated otherwise,
the patch positions are fixed at $\theta_{1}=5$ and $\theta_{2}=10$.

We start by demonstrating the effect of patch-patch interaction strength
$\epsilon_{p}$ on the formation of BCs, keeping the patch range constant at
$\alpha=0.1$. The important metrics to characterize the BCs are: (i) the
fraction of particles $f$ that have self-assembled into BCs, (ii) the average
bend angle $\bar{\gamma}$ of the BCs and (iii) the flexibility of the BCs.
Figs. 2(a)-(c) show prototypical equilibrated snapshots at time $t=10^{4}$ for
$\epsilon_{p}$ = 7.5, 15 and 26.25. The red and green patchy ellipsoids which
yield a BC unit are coloured purple. The asters, formed when more than two
ellipsoidal particles adhere, are coloured cyan. (See insets in Figs.
2(a)-(c).) The BC yield $f$ for different values of $\epsilon_{p}$ is shown in
Fig. 2(d). The inset shows the patch-patch interaction $V_{p}(d)$ (inverted
Gaussian) for the different values of $\epsilon_{p}$. The primary effect of
varying $\epsilon_{p}$ is to vary the well-depth. Fig. 2(e) shows that the
corresponding distribution of bend angles while Fig. 2(f) shows $f$ and
$\bar{\gamma}$ as a function of $\epsilon_{p}$. For $\epsilon_{p}=7.5$, $f$
remains quite low even after long simulation times. This is because $V_{p}$ in
this case is weaker than thermal energy $k_{B}T$ and is not sufficient to keep
the left-handed and right-handed particles bonded, see the morphology in Fig.
2(a). Fig. 2(e) shows that the corresponding distribution of bend angles is
broad. For an intermediate value of $\epsilon_{p}=15$, $V_{p}$ is sufficiently
strong with respect to $k_{B}T$, thus keeping the individual building blocks
bonded, as shown in Fig. 2(b). As a result, we achieve a BC yield close to
$1$. For $\epsilon_{p}=30$, $V_{p}$ is even stronger, making single bonds also
stable against thermal fluctuations. In addition to the BCs, this leads to the
creation of asters formed due to the patch at $\theta_{1}=5$. Single bonds
with more than one monomer of the opposite handedness are also stable in this
case as seen in the morphology of Fig. 2(c). Consequently, the BC yield is
also reduced, as observed in Fig. 2(d). The corresponding bend angle
distribution shown in Fig. 2(e), which excludes the asters, is narrower than
that for $\epsilon_{p}=7.5$, indicating the reduced flexibility of the BCs, as
is expected with stronger patch-patch interactions. Interestingly, for the
range of $\epsilon_{p}$ values over which the BC formation is efficient, we do
not notice a significant variation in the flexibility of the bend angle as
revealed by Fig. 2(f).

The other crucial parameter for controlling the self-assembly of BCs is the
patch interaction range $\alpha$. Graphically, it represents the well width of
inverted Gaussian potential. Guided by the observations in Fig. 2, we set the
patch strength $\epsilon_{p}=15$ to explore the effect of $\alpha$. Figs.
3(a)-(c) show the snapshots for $\alpha$ = 0.03, 0.1 and 0.15. For the small
value of 0.03, the effective region within which two complementary patches can
make a stable bond is significantly reduced. This not only requires the
bonding monomers to be close to each other but also that they are in near-
perfect alignment in order to make a double bond. As a result, the probability
of opposite-handed particles forming a BC is quite small. This negatively
affects the kinetics of the self-assembly leading to a low rate of BC
production, as seen in Fig. 3(d). The resultant BCs are also less flexible and
prone to dissociate into monomers upon the smallest fluctuations. As $\alpha$
increases to 0.1 in Fig. 3(b), the monomers have a higher probability of
making bonds and the rate of BC self-assembly increases. The flexibility of
BCs also increases as the increased patch range allows for a larger window of
angles over which the monomers can remain bonded, as shown in Fig 3(e) and
(f). Both these factors contribute to a larger BC yield. For larger patch
ranges $\alpha\geq 0.15$ in Figs. 3(f), the rate of self-assembly of BCs is
high. However, the yield is significantly reduced due to the formation of
asters. They coexist with the BCs as shown in Fig 3(c), and the BCs display a
large fluctuation in their bend angle for larger values of $\alpha$. From
Figs. 2(f) and 3(f), it can be seen that for fixed values of the patch angles,
there is a sizable region in the ($\epsilon_{p}-\alpha$) parameter space which
yields an efficient formation of BCs. The distribution of the bend angle
$\gamma$ is predominantly dictated by the patch range $\alpha$ and to a lesser
extent by the patch strength $\epsilon_{p}$.

It is also educative to understand the role of the patch placements
$[\theta_{1}$, $\theta_{2}]$ on the BCs. Choosing the $\epsilon_{p}=15$ and
$\alpha=0.1$, Fig. 4 shows typical morphologies for four different choices of
patch positions: (a) [0, 2]; (b) [0, 5]; (c) [3, 6]; (d) [5, 10]. The
corresponding bend angle distributions $P(\gamma)$ vs. $\gamma$ and the BC
yield $f$ are also shown in Figs. 4(e) and (f). Their choice affects the bend
angle. For example [0, 5] yields $\bar{\gamma}\simeq 145$, but shifted
placements [5, 10] results in $\bar{\gamma}\simeq 82$. Note that this angle is
dependent on the shape anisotropy $\kappa$ as well. For instance, the mean
bend angle would be larger for less eccentric ellipses, i.e., for $\kappa<3$.
BC yields are also affected as in cases such as [0, 2], there is a larger
possibility of non-BC aggregations such as the asters. Further, we can also
expect distinct phases and their textures as the patch positions are varied:
[0, 2] nudges the arrangement in a nematic phase while [0, 5] introduces a
splay in the texture. Our preferred choice has been [5, 10] as it yields BC
formation $\simeq 1$. The interplay of $\alpha$, $\theta_{1}$ and $\theta_{2}$
can be expected to yield a range of bend angles and novel phases. This is an
exhaustive investigation that deserves attention. However, first it is
essential to understand the methodology to create model BCs from patchy
ellipsoids and subsequently study their self-assembly into states with long-
range liquid crystalline order. We prioritize establishing these protocols in
the present work and will address the $\alpha$-$\theta_{1}$-$\theta_{2}$ state
diagram in an ensuing study.

###  IV.2 Second step assembly stage

We now turn our attention to the ordering assemblies from the initial
monomeric state. Figs. 5(a)-(d) show the time evolution of the morphology. The
corresponding BC yields $f$ are shown in Fig. 5(e). At time $t=0$, when the
system is quenched from $T=2$ to $1.2$, the system is in a homogeneous
isotropic state where the monomers have not formed any BCs [Fig. 5(a)]. In the
absence of any patches, this would remain an isotropic fluid, given our choice
of the GB parameters. Due to the presence of attractive patches, opposite-
handed monomers start adhering. The BC yield $f$ rapidly approaches 1 as time
increases. An unusual observation from Figs. 5(c) and (d) is that the evolving
ordered state has anti-ferroelectric smectic order. We also evaluate the
nematic and polar order parameters ($\mathcal{S}$ and $\mathcal{P}$) for the
evolving morphology in Fig. 5(f). Note that the polar direction $\hat{p}$ of
the neighbouring BCs in the ordered phase are antiparallel (or staggered).
Consequently, we evaluate the polarization by considering only the magnitude
of $\hat{p}$ along the direction of polarization ${\bf\boldsymbol{p}}$. In
other words, we take the absolute value of $(\hat{p}_{i}\cdot\mathbf{p})$
before averaging over all BCs. As seen in Fig. 5(f), $\mathcal{S}$ is close to
zero at early times as the system is in an isotropic state. As time evolves,
there is the emergence of BCs and subsequent liquid crystalline order due to
the GB interactions between the BC units. The BCs start forming small locally
aligned regions, which then coalesce into larger domains with liquid
crystalline order [Fig. 5(c) and (d)]. Consequently, there is an increase in
the nematic and polar order parameters, as seen in Fig. 5(f). These
observations, along with Figs. 5(e), emphasize the two-stage self-assembly in
our model: the BC formation followed by long-range liquid crystalline order in
the smectic anti-ferroelectric phase.

The full set of BC states (phases) for a specific choice of model parameters
($\epsilon_{p}=15$, $\alpha=0.1$) in the optimal window is shown in the
$\rho-T$ state diagram of Fig. 6(a). The different symbols indicate the
different states observed in our simulations: isotropic state (indigo
circles), anti-ferroelectric smectic (magenta triangles), arrested or frozen
state (green squares) and disordered state (blue diamonds). The background
color represents the BC yield. We see that there is a sizeable region in the
$\rho-T$ space where the BC yield is greater than $0.9$. Panels (b) to (e)
show typical morphologies obtained from MD (left) and MC (right) simulations
for the four states obtained at indicated parameter values of the
corresponding filled symbols in the state diagram. In (b), there is formation
of BCs, but the system is in an isotropic state as the Gay-Berne attraction
between them is not sufficient to overcome the thermal fluctuations. The
region between the dashed lines is where we obtain the smectic anti-
ferroelectric phase shown in the snapshot (c). For these parameter values, the
BCs are fluid enough to self-assemble into large ordered domains. At low
temperatures ($T<0.9$), there are locally aligned regions of BCs, but there is
a lack of long-range order as these clusters are frozen. Such a kinetically
arrested or frozen state due to low temperatures is shown in panel (d). At
very high densities $\rho>0.6$, the particles lack the free volume to move and
reorient, leading to a dense disordered BC state, as shown in panel (e).
Similar morphologies are also observed in the NPT ensemble with pressure-tuned
to achieve corresponding densities, see Fig. 7.

In order to characterize the structure of the phases, we evaluate the
directional radial distribution functions (RDF) resolved along axes in the
local frame of reference of the BCs as described in section II.2.3. This is
schematically illustrated in Fig. 8(a) where for each BC $i$, we compute the
RDF with other BCs along parallel or perpendicular directions to its local
alignment $\hat{n}_{i}$. Only those BCs that are within a finite region
perpendicular to the direction of the correlation evaluation are included. The
directional RDFs at $\rho=0.4$ and $T=1.5$, corresponding to the isotropic
phase, are shown in Fig. 8(b). We observe only a few peaks in
$g_{\perp}(r_{\perp})$ and hardly any in $g_{\parallel}(r_{\parallel})$,
indicating that at most, only two or three BCs stack along their short axis.
In the case of the smectic anti-ferroelectric phase ($\rho=0.4$ and $T=1.2$),
we observe multiple equidistant peaks along either direction, indicating the
arrangement of BCs into multiple layers, as shown in Fig. 8(c). Interestingly,
the correlations decay much more slowly in the direction of the long axis of
the BCs ($g_{\perp}(r_{\perp})$). In the frozen state at $\rho=0.4$ and
$T=0.9$ [Fig. 8(e)], we find equidistant peaks that decay rapidly,
corresponding to short-range correlations along either direction. This is due
to the presence of multiple small domains of locally aligned BCs. In the dense
disordered state ($\rho=0.7$ and $T=1.2$), we only find short-range
correlations along the short axis of the BCs, as seen from Fig. 8(f). In all
cases, our MD and MC results are in agreement with each other. The variation
in the correlation lengths along the two directions ($l_{n}$ and $l_{p}$) with
temperature $T$ is shown in Fig. 8(d), for two densities ($\rho=0.4$ and
$0.7$). At the higher density, the system is in a dense disordered state, and
hence the correlation lengths do not vary with temperature. At the lower
density, we see a gradual increase in the correlation lengths with a peak
centered around $T=1.2$, corresponding to the increased smectic anti-
ferroelectric alignment, beyond which the correlations decay rapidly as the
system enters the high-temperature isotropic phase. The plots also show that
the size of the LC domains along the long axis of the BCs is larger than that
in the orthogonal direction. The current investigation of parameters belonging
to patch-patch interactions provides only isotropic and smectic anti-
ferroelectric phases. This is because we have not fully explored the effect of
patch positions on the formation of mesoscopic phases. Our preliminary
explorations of the $\alpha$-$\theta_{1}$-$\theta_{2}$ space reveal a variety
of phases, e.g., nematic and splay-bends, through appropriate choice of the
model parameters.

##  V Summary and Conclusion

We have developed a patchy ellipsoid model capable of generating mesoscopic
liquid crystalline phases in a hierarchical two-stage self-assembly process,
via the formation of BC dimers in the first stage. This is achieved through a
binary mixture of ellipses having a pair of complementary and obliquely placed
patches with specific combination rules to guide the formation of BCs. The
patches are modeled as points on the surface of the ellipses interacting via
an inverted Gaussian potential. The second stage of the assembly into ordered
liquid phases is governed by the Gay-Berne interactions acting between the
ellipses themselves. Our model has three crucial parameters: (i) the patch
interaction strength $\epsilon_{p}$, (ii) the patch interaction range $\alpha$
and (iii) the patch positions $\theta_{1}$ and $\theta_{2}$, that control the
formation of the BCs, their mean bend angle and flexibility.

We perform extensive MD and MC simulations to elucidate the model and the
effect of the model parameters. Our main observations are summarized below:
(a) For fixed patch positions $\theta_{1}$ and $\theta_{2}$ and patch range
$\alpha$, we observe a window of $\epsilon_{p}$ values for which the BC yield
$f$ is close to $1$. For small values of $\epsilon_{p}$, the particles do not
make stable bonds, whereas for large values of $\epsilon_{p}$, assemblies of
more than two complementary monomers (asters) are formed due to increased
stability of single bonds.
(b) Similarly, with all other parameters fixed, we obtain high BC yields for a
wide range of values of $\alpha$. Although smaller values of $\alpha$ also can
produce BCs, we observe low yields in our simulations due to slow kinetics.
Larger values of $\alpha$ promote the formation of asters as the monomers can
make erroneous bonds due to increased range.
(c) For a choice of parameters in the optimal window, we study the evolution
of the BC yield and the nematic and polar orders, which show that the assembly
into ordered phases occurs in two-stages with the BC formation preceding the
long-range ordering.
(d) We obtain the state diagram of the model in the $\rho-T$ plane, which
exhibits isotropic, smectic anti-ferroelectric and disordered arrested states
(due to low temperatures or/and high densities) for the given choice of model
parameters. We compute correlation functions and order parameters of these
states and show that there is consistency between MD and MC simulations.
Further, there is consistency in the morphologies between NVT and NPT
ensembles as well, at appropriately chosen parameter values.
(e) The patch positions $\theta_{1}$ and $\theta_{2}$ can also be varied to
generate BCs with different configurations leading to a myriad of interesting
liquid crystalline phases, which we shall explore in subsequent works.

The assembly scenarios of these systems are very rich and, so far, the
available models for BC entities from patches have not been able to predict
all the observed BCLC phases. Our work sets the stage for these future
investigations. In contrast to the molecular BCs often studied in experiments,
our model provides a colloidal route to create BCLCs. Generating them at the
colloidal scale in experiments is challenging but could have huge potential
applications. We believe that our study offers one possible solution for the
experimental realization of these interesting fluids.

## Acknowledgements

AKS gratefully acknowledges the HPC computing facility of IIT Delhi. JT, EB
and GK gratefully acknowledge the generous share of CPU time, offered by the
Vienna Scientific Cluster (VSC) under project number fs71898. JT acknowledges
financial support from the EU-funded doctorate college ENROL. EB acknowledges
support from the Austrian Science Fund (FWF) under Proj. No. Y-1163-N27. VB
acknowledges SERB (India) for research grants CRG/2022/004871 and
CRG/2019/001836. AB acknowledges partial financial support from
CRG/2019/001836 during a visit to IIT Delhi that initiated this work.

Figure 1: (a) Schematic representation of the patchy colloidal model used in
this contribution with right-handed (green) and left-handed (red) elliptic
particles with positions $\boldsymbol{r}_{i}$, $\boldsymbol{r}_{j}$ and
orientations $\hat{u}_{i}$, $\hat{u}_{j}$. The right-handed particle has
patches $G_{1}$ and $G_{2}$ at angles $-\theta_{1}$ and $-\theta_{2}$ and the
left-handed particle have patches $R_{1}$ and $R_{2}$ at angles $-\theta_{1}$
and $-\theta_{2}$ (as indicated).
$r_{ij}=|\boldsymbol{r}_{j}-\boldsymbol{r}_{i}|$ represents the separation
between the center of masses of the two particles, and $d_{G_{1}R_{1}}$ the
separation between patches $G_{1},R_{1}$. (b) Schematic depiction of a BC
particle having nematic $\hat{n}$ and polar $\hat{p}$ directions and bend
angle $\gamma$. (c) Distribution of bend angles $P(\gamma)$ vs. $\gamma$
obtained using single patch particles (dashed line) and two patch particles
(bold line) for $\epsilon_{p}=15$, $\alpha=0.1$, $\theta_{1}=5$,
$\theta_{2}=10$, at $\rho=0.4$ and $T=1.2$. Figure 2: Prototypical
morphologies at $t=10^{4}$, $\alpha=0.1$, $T=1.2$ and $\rho=0.4$ for (a)
$\epsilon_{p}=7.5$, (b) $\epsilon_{p}=15$, and (c) $\epsilon_{p}=26.25$. In
the snapshots, BCs are shown in purple, and asters in cyan. (d) Corresponding
plots of BC yield $f$ vs. $t$. The inset depicts the inverted Gaussian
potential $V_{p}$ for $\epsilon_{p}$ = 7.5, 15, and 26.25 at fixed
$\alpha=0.1$. (e) Distribution of bend angles $P(\gamma)$ vs $\gamma$ with the
same color coding as in (d). The average bend angle $\bar{\gamma}$ and
standard deviation $\sigma_{\gamma}$ are mentioned within the plot. (f) Plots
of $f$ and $\bar{\gamma}$ as a function of patch interaction strength
$\epsilon_{p}$ at the end of the simulations. The error bars show the standard
deviation of the bend angle. Figure 3: Prototypical morphologies at
$t=10^{4}$, $\epsilon_{p}=15$, $T=1.2$ and $\rho=0.4$ for (a) $\alpha=0.03$,
(b) $\alpha=0.1$, and (c) $\alpha=0.15$. (d) Corresponding plots of BC yield
$f$ vs. $t$. The inset depicts the inverted Gaussian potential $V_{p}$ for
$\alpha$ = 0.03, 0.1, and 0.15 at fixed $\epsilon_{p}=15$. (e) Distribution of
bend angles $P(\gamma)$ vs $\gamma$ following the same color coding as panel
(d), with the corresponding $\bar{\gamma}$ and $\sigma_{\gamma}$ indicated
within the plot. (f) Plots of $f$ and $\bar{\gamma}$ as a function of $\alpha$
at the end of the simulations. Figure 4: Prototypical morphologies at
$\epsilon_{p}=15$, $\alpha=0.1$, $T=1.2$ and $\rho=0.4$ for different patch
positions: (a) $\theta_{1}=0$, $\theta_{2}=2$, (b) $\theta_{1}=0$,
$\theta_{2}=5$, (c) $\theta_{1}=3$, $\theta_{2}=6$, and (d) $\theta_{1}=5$,
$\theta_{2}=10$. (e) Corresponding bend angle distributions $P(\gamma)$ vs
$\gamma$. (f) Corresponding BC yields $f$ vs. $t$. Figure 5: Evolution
morphologies for $\epsilon_{p}=15$, $\alpha=0.1$, at $T=1.2$ and $\rho=0.4$
from a typical run at (a) $t=0$, (b) $t=250$, (c) $t=2000$, and (d) $t=10000$.
(e) Corresponding BC yields $f$ vs. $t$. The times corresponding to snapshots
(a)-(d) are also indicated (open circles). (f) Evolution of the nematic
($\mathcal{S}$) and polar ($\mathcal{P}$) order parameters with simulation
time $t$. Figure 6: (a) $\rho-T$ state diagram for the choice of parameters
$\epsilon_{p}=15$, $\alpha=0.1$, in the NVT ensemble. The symbols correspond
to different phases: isotropic (I, indigo circle), anti-ferroelectric smectic
(magenta triangle), arrested state due to low-temperature (AS1, green square),
and arrested state due to high density (AS2, blue diamond). The background
color represents the BC yield $f$ at the given state point. The region between
dashed lines corresponds to smectic anti-ferroelectric. Snapshots of typical
morphologies obtained from MD (left) and MC (right) simulations at state
points marked with filled symbols in panel (a): (b) $T=1.5$, $\rho=0.4$ (I),
(c) $T=1.2$, $\rho=0.4$ (smectic anti-ferroelectric), (d) $T=0.9$, $\rho=0.4$
(AS1), and (e) $T=1.2$, $\rho=0.7$ (AS2). Figure 7: The plot of BC yield $f$
vs. $t$ for fixed values of $\epsilon_{p}=15$ and $\alpha=0.1$ in NPT ensemble
for different pressures $P$ and temperatures $T$. Corresponding typical final
morphologies at $t=10^{4}$: (b) $P=0.01$, $T=1.5$; (c) $P=0.01$, $T=1.2$; (d)
$P=0.01$, $T=0.9$; (e) $P=0.8$, $T=1.2$. (f) Corresponding density $\rho$
variation with simulation time $t$. The morphologies resemble those in Fig. 6
for corresponding values of $\rho$ and $T$. Figure 8: (a) Schematic
illustration of the evaluation of the directional RDFs along the polar and
nematic alignment of the central BC, with correlation lengths $l_{n}$ along
nematic and $l_{p}$ along polar directions. Plots of the directional RDFs
$g_{\parallel}(r_{\parallel})$ (in red), $g_{\perp}(r_{\perp})$ (in blue) vs.
$r_{\parallel},r_{\perp}$ in the NVT ensemble for various values of $T$ and
$\rho$: (b) $T=1.5$, $\rho=0.4$, (c) $T=1.2$, $\rho=0.4$, (e) $T=0.9$,
$\rho=0.4$, and (f) $T=1.2$, $\rho=0.7$. All data correspond to the choice of
model parameters: $\epsilon_{p}=15$, $\alpha=0.1$, $\theta_{1}=5$, and
$\theta_{2}=10$. (d) Correlation lengths $l_{n}$ and $l_{p}$ as a function of
temperature $T$ at densities $\rho$ = 0.4 and $\rho$ = 0.7.

## References

  * Lubensky and Radzihovsky (2002) T. Lubensky and L. Radzihovsky, Phys. Rev. E 66, 031704 (2002).
  * 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).
  * Takezoe and Takanishi (2006) H. Takezoe and Y. Takanishi, Jpn. J. Appl. Phys. 45, 597 (2006).
  * Jákli  _et al._ (2018) A. Jákli, O. D. Lavrentovich, and J. V. Selinger, Rev. Mod. Phys. 90, 045004 (2018).
  * Johnston  _et al._ (2002) S. J. Johnston, R. J. Low, and M. P. Neal, Phys. Rev. E 65, 051706 (2002).
  * Bates and Luckhurst (2005) M. A. Bates and G. R. Luckhurst, Phys. Rev. E 72, 051702 (2005).
  * Bates (2006) M. A. Bates, Phys. Rev. E 74, 061702 (2006).
  * Dewar and Camp (2004) A. Dewar and P. J. Camp, Phys. Rev. E 70, 011704 (2004).
  * Dewar and Camp (2005) A. Dewar and P. J. Camp, J. Chem. Phys. 123, 174907 (2005).
  * Chiappini  _et al._ (2019) M. Chiappini, T. Drwenski, R. van Roij, and M. Dijkstra, Phys. Rev. Lett. 123, 068001 (2019).
  * Berardi  _et al._ (1995) R. Berardi, C. Fava, and C. Zannoni, Chem. Phys. Lett. 236, 462 (1995).
  * Pawar and Kretzschmar (2010) A. Pawar and I. Kretzschmar, Macromol. Rapid Commun. 31, 150 (2010).
  * Bianchi  _et al._ (2011) E. Bianchi, R. Blaak, and C. N. Likos, Phys. Chem. Chem. Phys. 13, 6397 (2011).
  * Ravaine and Duguet (2017) S. Ravaine and E. Duguet, Current opinion in colloid & interface science 30, 45 (2017).
  * Bianchi  _et al._ (2017) E. Bianchi, B. Capone, I. Coluzza, L. Rovigatti, and P. D. J. van Oostrum, Phys. Chem. Chem. Phys. 19, 19847 (2017).
  * Smallenburg  _et al._ (2013) F. Smallenburg, L. Leibler, and F. Sciortino, Phys. Rev. Lett. 111, 188002 (2013).
  * Chen  _et al._ (2011) Q. Chen, E. Diesel, J. K. Whitmer, S. C. Bae, E. Luijten, and S. Granick, J. Am. Chem. Soc. 133, 7725 (2011).
  * Wang  _et al._ (2012) Y. Wang, Y. Wang, D. R. Breed, V. N. Manoharan, L. Feng, A. D. Hollingsworth, M. Weck, and D. J. Pine, Nature 491, 51 (2012).
  * Morphew  _et al._ (2018) D. Morphew, J. Shaw, C. Avins, and D. Chakrabarti, ACS Nano 12, 2355 (2018).
  * He  _et al._ (2020) M. He, J. P. Gales, É. Ducrot, Z. Gong, G.-R. Yi, S. Sacanna, and D. J. Pine, Nature 585, 524 (2020).
  * Posnjak  _et al._ (2024) G. Posnjak, X. Yin, P. Butler, O. Bienek, M. Dass, I. D. Sharp, and T. Liedl, Science 384, 781 (2024).
  * Liu  _et al._ (2024) H. Liu, M. Matthies, J. Russo, L. Rovigatti, R. P. Narayanan, T. Diep, D. McKeen, O. Gang, N. Stephanopoulos, F. Sciortino, H. Yan, F. Romano, and P. Šulc, Science 384, 776 (2024).
  * Karner and Bianchi (2024) C. Karner and E. Bianchi, Nanoscale Advances 6, 443 (2024).
  * Russo  _et al._ (2009) J. Russo, P. Tartaglia, and F. Sciortino, J. Chem. Phys. 131 (2009).
  * Gay and Berne (1981) J. Gay and B. Berne, J. Chem. Phys. 74, 3316 (1981).
  * Berardi and RJ (1993) R. Berardi and A. RJ, J. Chem. Soc. Farady Trans 89, 4069 (1993).
  * Bates and Luckhurst (1999) M. Bates and G. Luckhurst, J. Chem. Phys. 110, 7087 (1999).
  * Zannoni (2001) C. Zannoni, J. Mater. Chem. 11, 2637 (2001).
  * De Gennes and Prost (1993) P.-G. De Gennes and J. Prost,  _The physics of liquid crystals_ , 83 (Oxford university press, 1993).
  * Wei and Patey (1992) D. Wei and G. Patey, Phys. Rev. Lett. 68, 2043 (1992).
  * Allen and Tildesley (2017) M. P. Allen and D. J. Tildesley,  _Computer simulation of liquids_ (Oxford university press, 2017).
  * Frenkel and Smit (2023) D. Frenkel and B. Smit,  _Understanding molecular simulation: from algorithms to applications_ (Elsevier, 2023).
  * Billeter and Pelcovits (1998) J. Billeter and R. Pelcovits, Computer in Physics 12, 440 (1998).
  * Puri (2004) S. Puri, Phase Transitions 77, 407 (2004).
  * Weis and Levesque (1993) J. J. Weis and D. Levesque, Phys. Rev. E 48, 3728 (1993).
  * Cañeda-Guzmán  _et al._ (2014) E. Cañeda-Guzmán, J. Moreno-Razo, E. Díaz-Herrera, and E. Sambriski, Mol. Phys. 112, 1149 (2014).
  * Birdi  _et al._ (2022) N. Birdi, T. L. Underwood, N. B. Wilding, S. Puri, and V. Banerjee, Phys. Rev. E 105, 024706 (2022).
  * (38) https://lammps.sandia.gov.
  * Plimpton (1995) S. Plimpton, J. Comput. Phys. 117, 1 (1995).
  * Hoover (1985) W. G. Hoover, Phys. Rev. A 31, 1695 (1985).
  * Hoover (1999) W. G. Hoover,  _Time Reversibility, Computer Simulation, and Chaos_ (World Scientific, Singapore, 1999).
  * Andersen (1983) H. C. Andersen, J. Comput. Phys. 52, 24 (1983).
  * Whitelam and Geissler (2007) S. Whitelam and P. L. Geissler, J. Chem. Phys. 127, 154101 (2007).
  * Whitelam  _et al._ (2009) S. Whitelam, E. H. Feng, M. F. Hagan, and P. L. Geissler, Soft Matter 5, 1251 (2009).
  * (45) L. O. Hedges, “Libvmmc,” .