# Orbital order in a bosonic $p$-band triangular lattice

Hua Chen  Electronic address<EMAIL_ADDRESS>he/him/his Department of
Physics, Zhejiang Normal University, Jinhua 321004, China  X. C. Xie
International Center for Quantum Materials, School of Physics, Peking
University, Beijing 100871, China  Beijing Academy of Quantum Information
Sciences, Beijing 100193, China  CAS Center for Excellence in Topological
Quantum Computation, University of Chinese Academy of Sciences, Beijing
100190, China

###### Abstract

We present a detailed study of the Bose-Hubbard model in a $p$-band triangular
lattice by focusing on the evolution of orbital order across the superfluid-
Mott insulator transition. Two distinct phases are found in the superfluid
regime. One of these phases adiabatically connects the weak interacting limit.
This phase is characterized by the intertwining of axial $p_{\pm}=p_{x}\pm
ip_{y}$ and in-plane $p_{\theta}=\cos\theta p_{x}+\sin\theta p_{y}$ orbital
orders, which break the time-reversal symmetry and lattice symmetries
simultaneously. In addition, the calculated Bogoliubov excitation spectrum
gaps the original Dirac points in the single-particle spectrum but exhibits
emergent Dirac points. The other superfluid phase in close proximity to the
Mott insulator with unit boson filling shows a detwined in-plane ferro-orbital
order. Finally, an orbital exchange model is constructed for the Mott
insulator phase. Its classical ground state has an emergent SO$(2)$ rotational
symmetry in the in-plane orbital space and therefore enjoys an infinite
degeneracy, which is ultimately lifted by the orbital fluctuation via the
order by disorder mechanism. Our systematic analysis suggests that the in-
plane ferro-orbital order in the Mott insulator phase agrees with and likely
evolves from the latter superfluid phase.

##  I Introduction

Orbital order is a long-standing issue tracing back to the transition metal
oxides [1, 2, 3]. The precise mechanism driving orbital order remains largely
unknown due to the intricate interplay among spin, orbital, charge, and
lattice degrees of freedom in host crystals. In particular, a recent example
of relevance is the nematic phase in iron-based superconductors, which
entwines with spin Ising order, orbital order, and lattice structural
distortion as dictated by symmetry [4]. Among diverse theoretical proposals in
addressing the origin of nematicity [5, 6], one interesting finding is that
the orbital order in the nematic phase manifests its essential role in the
metal-insulator transition and promotes an intermediate phase, i.e., the
orbital-selective Mott phase [7]. This phase is characterized by the orbital-
dependent Mott localization and interpolates the itinerant and Mott localized
limits, validating the incipient Mott picture [8, 9]. By contrast, a natural
question may raise for bosonic systems: how the orbital order evolves in the
superfluid-Mott insulator (SF-MI) transition.

Yet, much efforts have been denoted to the understanding of orbital order in
electronic materials. While, the studies in bosonic systems are rare [10, 11].
Experimentally, artificial systems, such as ultracold atomic [12, 13, 14, 15,
16] and photonic [17, 18, 19] systems, have been shown the exciting
possibility of stimulating the crystals with $p$-orbital bosons in the first
excited band. For instance, the Dirac points in the $p$-band hexagonal lattice
are theoretically predicted by the early study [20] and experimentally
observed in photonic systems [17]. Later, the orbital edge state, which is
extensively studied in graphene [21, 22, 23], is confirmed in the subsequent
photonic experiment [18]. More recently, the evidence of nematic superfluid
(SF) phase in a hexagonal lattice, which is attributed to the orbital order,
is also reported in ultracold atomic systems [24].

The main purpose of our study is to give a comprehensive understanding of the
orbital order in the SF-MI transition. The single-particle spectrum of the
$p$-band triangular lattice exhibits a pair of Dirac points at the corners of
hexagonal Brillouin zone (HBZ), resembling the low energy physics of graphene
[21, 22, 23]. The evolution of orbital order across the SF-MI transition is
then studied based on the Bose-Hubbard model. In the weak-interacting limit,
the $p$-band triangular lattice is frustrated due to the inability to
simultaneously minimize both the kinetic and interacting energies. This weak-
coupling SF phase is characterized by the intertwining of the axial
$p_{\pm}=p_{x}\pm ip_{y}$ and in-plane $p_{\theta}=\cos\theta p_{x}+\sin\theta
p_{y}$ orbital orders. Interestingly, the Bogoliubov excitation spectrum in
this SF phase gaps the original Dirac points in the single-particle spectrum
but exhibits emergent Dirac points. In the strong-interacting limit, the
orbital order is also studied based on the orbital exchange model. We show
that the classical ground state is of ferro-orbital type and enjoys an
emergent SO$(2)$ rotational symmetry, which ensures an infinite degeneracy.
The orbital fluctuation ultimately lifts the degeneracy and selects discrete
quantum ground states through the order by disorder mechanism. Moreover, the
phase diagram established by Gutzwiller approach interpolates these two
limits. Besides these two phases, we find an intermediate SF phase with the
ferro-orbital order. This intermediate phase survives in a wide range of low
boson filling and gradually increase the occupation in the preferable in-plane
orbital when approaching the Mott insulator (MI) with unit boson filling
$n=1$. Our study provides strong clues that the ferro-orbital order in the MI
$n=1$ phase likely evolves from the intermediate SF phase, facilitating the
understanding on the role of orbital order in the SF-MI transition.

The reminder of this paper is organized as follow. In Sec. II, we introduce
the $p$-band tight-binding model in the triangular lattice as well as the
Bose-Hubbard model. We establish the ground-state phase diagram by utilizing
Gutzwiller approach in Sec. III. The orbital order in the weak-interacting
limit is further studied with Bogoliubov approximation by treating the Bose-
Hubbard interaction perturbatively in Sec IV. In Sec. V the orbital exchange
model is constructed to study the orbital order by treating the hopping
processes as perturbations. Finally, we summarise and discuss the results in
Sec. VI.

Figure 1:  (a) The structure of triangular lattice and the hexagonal Brillouin
zone. The blue lines mark the Wigner-Seitz cell and the reduced Brillouin zone
due to the umklapp scattering between the band minima at $M$ points. (b) The
band structure of the tight-binding model in Eq. (1) with the $\sigma$-bonding
$t_{\sigma}=1$ being the energy unit. The pseudovector fields
$\bm{d}\equiv\left(d_{z},d_{x}\right)$ near the Dirac point at $K$ (c) and
$K^{\prime}$ (d) resemble the vortex in $XY$ systems with winding number
$W=1$.

##  II Minimal model

We begin with the tight-binding model that describes the hopping processes of
bosons in the $p$-band triangular lattice depicted in Fig. 1(a). Introducing
an orbital pseudospin representation, the momentum-space Hamiltonian in the
basis $p_{\bm{k}}=\left[p_{x{\bm{k}}},p_{y{\bm{k}}}\right]^{\text{T}}$ reads

$\mathcal{H}_{\bm{k}}=d_{0}\left({\bm{k}}\right)\tau_{0}+d_{x}\left({\bm{k}}\right)\tau_{x}+d_{z}\left({\bm{k}}\right)\tau_{z},$
(1)

where $\tau_{0}$ and $\bm{\tau}$ are the identity matrix and Pauli matrices
respectively, and the coefficients
$d_{0}\left(\bm{k}\right)=\left(t_{\sigma}+t_{\pi}\right)\sum_{i}\cos k_{i}$,
and
$\\{d_{x}\left(\bm{k}\right),d_{z}\left(\bm{k}\right)\\}=\left(t_{\sigma}-t_{\pi}\right)/2\times\\{\sqrt{3}\left(\cos
k_{3}-\cos k_{2}\right),\cos k_{1}+\sum_{i}\cos k_{i}\\}$. Here, the crystal
momenta $\\{k_{1},k_{2},k_{3}\\}$ are measured along reciprocal lattice
vectors $\\{\bm{b}_{1},\bm{b}_{2},\bm{b}_{3}\equiv-\bm{b}_{1}-\bm{b}_{2}\\}$,
and the hopping integrals $t_{\sigma}$ and $t_{\pi}$ denote the $\sigma$ and
$\pi$ bonding of $p$ orbitals, respectively. For the $\pi$ bonding, the bond
vector lies in the nodal plane of $p$ orbitals. As a result, the strength of
$\pi$ bonding is typically much weaker than that of $\sigma$ bonding. The band
structure of the tight-binding model in Eq. (1) is plotted in Fig. 1(b).
Notably, two bands cross at the Dirac points located at $K$ and $K^{\prime}$
points of HBZ. To describe the corresponding low-energy behavior around $K$
and $K^{\prime}$ points, we derive the effective $k\cdot p$ model

$\displaystyle\mathcal{H}_{K/K^{\prime}}\left({\bm{q}}\right)=d_{0}\tau_{0}+d_{x}\tau_{x}+d_{z}\tau_{z}+\mathcal{O}\left(q^{2}\right)$
(2)

with the coefficients

$\displaystyle
d_{0}=-\frac{3}{2}\left(t_{\sigma}+t_{\pi}\right),\\{d_{x},d_{z}\\}=\pm\frac{3}{4}\sqrt{3}\left(t_{\sigma}-t_{\pi}\right)\\{q_{x},-q_{y}\\}.$

Diagonalizing $\mathcal{H}_{K/K^{\prime}}\left({\bm{q}}\right)$ gives two non-
interacting bands
$E^{\pm}_{K/K^{\prime}}\left({\bm{q}}\right)=d_{0}\pm\sqrt{d_{x}^{2}+d_{z}^{2}}$,
resulting in a linear dispersed Dirac point with the velocity
$v=3\sqrt{3}/4(t_{\sigma}-t_{\pi})$. The pseudovector fields
$\bm{d}\equiv\left(d_{z},d_{x}\right)$ around $K$ and $K^{\prime}$ points,
shown in Figs. 1(c) and 1(d) respectively, have a $p$-wave symmetry. The
topological charge of Dirac point is given by the winding number of
pseudovector field:
$W=\frac{1}{2\pi}\oint_{\mathcal{C}}\nabla\theta\left({\bm{q}}\right)\cdot
d{\bm{q}}=1$, where $\theta\equiv\text{arctan}\left(d_{x}/d_{z}\right)$ and
$\mathcal{C}$ is a contour enclosing the singular $K/K^{\prime}$ point,
indicating that the Dirac point carries a $\pi$ Berry flux. The band minima
are located at three inequivalent centres $M$ of HBZ edges, promoting a
finite-momentum Bose-Einstein condensate for weakly interacting bosons. For
non-interacting bosons, an infinite degenerate manifold of the single-particle
ground state can be constructed by the linear superposition of the Bloch
functions at these band minima. The umklapp scattering between the band minima
transfers a lattice phonon which carries the momentum of multiple primitive
reciprocal vectors. This process folds three $M$ points to $\Gamma$ point and
underlies the reduced Brillouin zone (RBZ) and the enlarged Wigner-Seitz cell,
as illustrated by the blue lines in Fig. 1(a).

Figure 2:  Gutzwiller approach. (a) The ground-state phase diagram in $\mu/U$
vs $U/t_{\sigma}$ plane accommodates three distinct phases, including (1) ISF,
the superfluid phase intertwining axial and in-plane orbital orders; (2) DSF,
the superfluid phase with the detwinned in-plane orbital order, and (3) MI
$n=1$, the Mott insulator phase with unit filling $n=1$. The transition lines
separating these three phases merge at a triple point
$\left(U/t_{\sigma},\mu/U\right)\approx\left(12.2,0.66\right)$. (b) Evolution
of amplitude of condensate order parameter
$\bm{\phi}=\left(\phi_{x},\phi_{y}\right)$, compressibility $\kappa=\partial
n/\partial\mu$, and orbital pseudospin $\bm{\tau}$ at fixed $\mu/U=0.5$. The
black and blue dots in (a) mark the the critical points with the vanishing
order parameters $|\bm{\phi}|$ and $|\tau_{y}|$, respectively. In the
numerical calculations, the truncation of local Fock space $N_{\pm}=10$ for
the maximum occupation in axial orbitals $p_{\pm}=p_{x}\pm ip_{y}$ and
$t_{\pi}=0$ are used.

Having established the single-particle physics, we are then in a position to
study the effects of many-particle interactions. The interacting Hamiltonian
can be generally constructed in terms of Haldane pseudopotentials by
projecting a pair of particles into relative angular momenta, respecting the
quantum statistics [25, 26]. The Bose-Hubbard interaction which is
mathematically described by zero relative angular momentum takes the form

$H_{\text{I}}=\frac{3}{2}U\sum_{i}\left[\hat{n}_{i}\left(\hat{n}_{i}-\frac{2}{3}\right)-\frac{1}{3}\hat{L}_{zi}^{2}\right],$
(3)

where $\hat{n}_{i}=\sum_{\alpha=x,y}p_{\alpha i}^{\dagger}p_{\alpha i}$ is the
occupation operator and
$\hat{L}_{zi}=-i\sum_{\alpha,\beta=x,y}\epsilon_{z\alpha\beta}p_{\alpha
i}^{\dagger}p_{\beta i}$ is the $z$-component orbital angular momentum at
$i$-th site [27, 28, 29]. Here $\epsilon_{\alpha\beta\gamma}$ is the Levi-
Civita antisymmetric tensor. The interaction can be experimentally realized
through the Feshbach resonance [30] and optical nonlinearities [31] for
ultracold atomic and photonic systems, respectively.

##  III Gutzwiller approach

To gain an overall understanding on the ground-state phase diagram, Gutzwiller
approach [32, 33, 34, 35] has its advance in capturing the physics in the
intermediate regime of Hubbard interaction $U$, and straddles the limits of
the weakly interacting SF and strongly interacting MI phases. It has been
utilized to establish the phase diagram of $p$-band Bose-Hubbard model with
different lattice geometries in the early study [27]. This approach starts
from the factorized local Fock state

$\displaystyle\left|\Psi_{\text{GW}}\right\rangle=$
$\displaystyle\prod_{i}\sum_{\text{F}}\eta_{\text{F}}^{i}\left|\text{F}\right\rangle_{i},$
$\displaystyle\left|\text{F}\right\rangle_{i}=\frac{1}{\sqrt{n_{+}^{\text{F}}!n_{-}^{\text{F}}!}}$
$\displaystyle\left(p_{+i}^{\dagger}\right)^{n_{+}^{\text{F}}}\left(p_{-i}^{\dagger}\right)^{n_{-}^{\text{F}}}\left|0\right\rangle$

where $n_{\pm}^{\text{F}}$ is the occupation of bosons in the axial orbitals
$p_{\pm}=p_{x}\pm ip_{y}$ and $\eta_{\text{F}}^{i}$ is the probability
weighting factor determined variationally. It takes into account that the
multioccupation of bosons in the local orbitals are energetically costly. In
numerical calculations, a truncation of the local Fock space is imposed and
the filling of bosons is dictated by chemical potential $\mu$ for the grand
canonical ensemble. We decompose the hopping terms in the tight-binding
Hamiltonian as $p^{\dagger}_{\alpha i}p_{\beta j}\approx p^{\dagger}_{\alpha
i}\phi_{\beta j}+\phi^{*}_{\alpha i}p_{\beta j}-\phi^{*}_{\alpha i}\phi_{\beta
j}$ with the condensate order parameter $\phi_{\alpha
i}=\sum_{\text{F}\text{F}^{\prime}}\eta_{\text{F}}^{i*}\eta_{\text{F}^{\prime}}^{i}{}_{i}\left\langle\text{F}|p_{\alpha
i}|\text{F}^{\prime}\right\rangle_{i}$ entangling the local Fock states. The
self-consistent solution of the order parameters $\phi_{\alpha i}$ requires an
iterative minimization of the energy functional over the Wigner-Seitz cell in
Fig. 1(a). The calculated phase diagram for $t_{\pi}=0$ shown in Fig. 2 (a)
accommodates three different phases including two distinct SF phases and the
MI phase. We have also verified that the phase diagram remains qualitatively
robust against the perturbative $\pi$ bonding $t_{\pi}=-0.1t_{\sigma}$. To
characterize the orbital order, we numerically evaluate the ground-state
expectation of orbital pseudospin
$\bm{\hat{\tau}}=\sum_{\alpha\beta}p^{\dagger}_{\alpha}\bm{\tau}_{\alpha\beta}p_{\beta}$.
The axial orbital order $p_{\pm}=p_{x}\pm ip_{y}$ is characterized by the
orbital polarization $\tau_{y}$ in $y$ axis, while the in-plane orbital order
$p_{\theta}=\cos\theta p_{x}+\sin\theta p_{y}$ directing at angle $\theta$
with $x$ axis corresponds to the orbital polarization
$\left(\tau_{z},\tau_{x}\right)=\tau\left(\cos\left[2\theta\right],\sin\left[2\theta\right]\right)$
in $zx$ plane. Figure 2(b) shows the detailed evolution of order parameters at
fixed $\mu/U=0.5$, which determines the phases across the SF-MI transition.
The stability of each phase is further checked with various sets of supercell
sizes up to $8\bm{a}_{1}\times 8\bm{a}_{2}$. As shown in Fig. 2(b), two
distinct SF phases share a non-vanishing uniform order parameter
$|\bm{\phi}_{i}|$ and can be distinguished by the ground-state expectation
value of orbital pseudospin $\hat{\bm{\tau}}$. Initially, the ground state at
weak Hubbard interaction $U$ develops an intertwined order by entangling both
the axial and in-plane orbital orders. The former is characterized by the
alternating signs in adjacent rows but an identical amplitude of
$\tau_{y}^{i}$ therefore suggesting the ordering of antiferro-orbital angular
momentum, while the latter is indicated by the uniform
$\left(\tau^{i}_{z},\tau^{i}_{x}\right)$ implying the ferro-orbital order.
Therefore, this superfluid phase intertwining axial and in-plane orbital
orders is denoted as ISF. The detailed pattern of orbital orders will be
further discussed in Sec. IV. At the critical Hubbard interaction $U_{c1}$,
the orbital pseudospin is completely aligned in the $zx$ plane
$\left(\tau_{z},\tau_{x}\right)$, showing a detwinned ferro-orbital order.
This superfluid phase is thus denoted as DSF. Since the inter-site hopping
process is treated at the mean-filed level in Gutzwiller approach, we will
show in Sec. V that the orientation of in-plane orbital is solely determined
by the quantum fluctuation due to the orbital anisotropy. With further
increasing Hubbard interaction $U>U_{c2}$, the MI phase obtained within
Gutzwiller approximation is simply a product of local Fock states and is thus
featureless. As shown in Fig. 2(b), the ISF-DSF-MI transition driven by he
Hubbard interaction $U$ is well detected by the discontinuous jumps of the
compressibility $\kappa=\partial n/\partial\mu$. Therefore, the phase
transition discussed here may be experimentally probed by measuring the boson
filling $n$. Below we shall justify the orbital orders above from two extrema
limits.

Figure 3:  Bogoliubov approach. (a) The condensed ground state intertwines the
ferro-orbital order with the in-plane orbital being parallel to the bond, and
the antiferro-orbital angular momentum with the spatial pattern indicated by
red and blue lines, supporting staggered fluxes. (c) The Bogoliubov excitation
spectrum along the high symmetry line indicated in (b). (d) The detailed
Bogoliubov excitation spectrum indicated by the pink box in (c) exhibits Dirac
bosons. (e) The pseudovector field $\bm{d}$ of Dirac boson in the lowest two
excitation branches around $M$ point has winding number $W=1$. The parameters
used in numerics are $\left(t_{\sigma},t_{\pi},nU\right)=\left(1,0,1\right)$.

##  IV Weak-coupling approach

In the weakly interacting limit, the Hubbard interaction $U$ is treated
perturbatively. The operators can be decomposed in terms of quantum
fluctuations $\tilde{p}_{\alpha\bm{k}\ell}$ around the condensed ground-state
wave function $\phi_{\alpha\ell}$

$p_{\alpha\bm{k}\ell}=\phi_{\alpha\ell}\delta\left(\bm{k}\right)+\tilde{p}_{\alpha\bm{k}\ell}\text{,
}\alpha=x,y$ (4)

with $\ell$ specifying the sublattice in the Wigner-Seitz cell. In the spirit
of Bogoliubov approximation [36, 37], the Hamiltonian is expanded in powers of
quantum fluctuations and is truncated up to the quadratic order. Detailed
derivations are presented in Appendix A. The zeroth-order terms in this
expansion determine the energy functional
$\varepsilon\left(\bm{\phi}^{*},\bm{\phi}\right)$ of the condensate at
$\Gamma$ point in RBZ. The time-dependent Gross-Pitaevskii equation can be
readily derived via the Euler-Lagrange equation

$\frac{\partial\mathcal{L}}{\partial\phi_{\alpha\ell}^{*}}-\frac{d}{dt}\left(\frac{\partial\mathcal{L}}{\partial\dot{\phi}_{\alpha\ell}^{*}}\right)=0,$
(5)

where the Lagrangian
$\mathcal{L}=\sum_{\alpha\ell}i\hbar\left(\phi_{\alpha\ell}^{*}\dot{\phi}_{\alpha\ell}-\phi_{\alpha\ell}\dot{\phi}_{\alpha\ell}^{*}\right)-\varepsilon\left(\bm{\phi}^{*},\bm{\phi}\right)$
[38]. The ground state can be numerically solved through the imaginary-time
evolution of Gross-Pitaevskii equation by propagating an initial trial state
[39]. Mathematically, this procedure is equivalent to the minimization of the
energy functional $\varepsilon\left(\bm{\phi}^{*},\bm{\phi}\right)$, which
causes the linear order terms in $\tilde{p}_{\alpha\bm{k}\ell}$ to vanish. The
calculated ground-state condensate develops an intertwined order, confirming
the results from Gutzwiller approach. The ferro-orbital order is characterized
by orientating the in-plane orbital parallel to the bond, which breaks the
lattice rotational symmetry. The antiferro-orbital angular momentum is
characterized by the alternating sign of $\tau_{y}$ along the direction
perpendicular to the in-plane orbital, and breaks the time-reversal symmetry
as well as the lattice translational symmetry. Interestingly, the symmetry
breaking of this weak-coupling phase, which has been studied in details, is
shown to be universal in the strong-coupling regime with boson filling $n\geq
2$ [40]. As schematically depicted in Fig. 3(a), the staggered flux pattern of
the ground-state condensate arising from the time-reversal symmetry breaking
is characterized by the bond current
$J_{ij}=-i\sum_{\alpha\beta}t_{\gamma}^{\alpha\beta}(\langle
p^{\dagger}_{\alpha i}p_{\beta
j}\rangle-\text{c.c.})\delta_{j,i\pm\bm{a}_{\gamma}}$ where the hopping matrix
$t_{\gamma}=[(t_{\sigma}+t_{\pi})\tau_{0}+(t_{\sigma}-t_{\pi})(\cos[2\theta_{\gamma}]\tau_{z}+\sin[2\theta_{\gamma}]\tau_{x})]/2$
and $\theta_{\gamma}$ is the azimuthal angle of $\bm{a}_{\gamma}$. The early
studies find pure axial orbital orders with different lattice geometries,
which support bond currents as a natural consequence [28, 41, 42]. In
contrast, the intertwined orbital order in the present study, due to the
inability to simultaneously minimize both kinetic and interacting energies,
originates from the geometric frustration of the triangular lattice. Having
settled the ground state, we then proceed with the quadratic order

$\mathcal{H}^{(2)}_{\bm{k}}=\frac{1}{2}\left[\tilde{\bm{p}}^{\dagger}_{\bm{k}},\tilde{\bm{p}}_{-\bm{k}}\right]\left[\begin{matrix}X_{\bm{k}}&Y\\\
Y^{\dagger}&X_{-\bm{k}}\end{matrix}\right]\left[\begin{matrix}\tilde{\bm{p}}_{\bm{k}}\\\
\tilde{\bm{p}}^{\dagger}_{-\bm{k}}\end{matrix}\right],$ (6)

which describes the Bogoliubov excitation on top of the ground-state
condensate. The diagonal terms $X_{\pm\bm{k}}$ in Eq. (6) receive
contributions from both the hopping processes and the self-energy correction
of Hubbard interaction $U$, while the off-diagonal terms $Y$ ($Y^{\dagger}$)
describe the anomalous processes in which a pair of bosons scatter with each
other into the excited states (condensates). As shown in Fig. 3(c), the
Bogoliubov excitation spectrum has a gapless Goldstone mode around $\Gamma$
point in RBZ as the signature of $U(1)$ symmetry breaking. Interestingly, we
find that the Bogoliubov spectrum along the high symmetry line
$K$-$K^{\prime}$ gaps out the original Dirac points at $K$ and $K^{\prime}$
points in the band structure but exhibits emergent Dirac bosons at $M$ point.
As a representative example, the Dirac point in the lowest two branches is
identified by numerically calculating the winding number shown in Fig. 3(e).
It is worth mentioning that the Dirac bosons only exist on one edge of RBZ
parallel to the bond that is selected by the in-plane orbital order, serving
as a fingerprint of lattice rotational symmetry breaking. We have also checked
that the Dirac related physics is robust for $t_{\pi}=-0.1t_{\sigma}$. While
the early study focuses on the Dirac fermion in the band structure with
staggered fluxes [43], the present study investigates the Dirac bosons in the
elementary excitation on top of the Bose-Einstein condensate with staggered
fluxes instead.

Figure 4:  Strong-coupling approach. (a) The eigenvalues of orbital
interaction matrix $\Lambda_{\bm{k}}$ within the classical approximation. (b)
The emergent SO$(2)$ rotational symmetry in the easy plane
$(\tau_{z},\tau_{x})$ of orbital pseudospin space. (c) The zero-point energy
$E_{\text{ZP}}\left(\theta\right)$ in Eq. (8) from the correction of orbital
fluctuations has minima at $\theta=\pi/6+\mathbb{Z}\pi/3$. The energies in (a)
and (c) are in units of $t_{\sigma}^{2}/16U$.

##  V Strong-coupling approach

Finally, we turn to the strong-coupling limit in which the virtual hopping
processes are treated perturbatively. Since the charge excitation in MI phase
is suppressed by the charge gap proportional to Hubbard interaction $U$, the
orbital fluctuation is the remaining low energy degree of freedom. Following
the standard second-order perturbation theory, the effective low-energy
physics for the MI $n=1$ phase is captured by the following orbital exchange
model

$\displaystyle H_{\text{OE}}=J\sum_{\langle
ij\rangle\parallel\bm{a}_{\gamma}}\tau^{i}_{\gamma}\tau^{j}_{\gamma}+J^{\prime}\sum_{\langle
ij\rangle}\left(\bm{\tau}^{i}\cdot\bm{\tau}^{j}+2\tau^{i}_{y}\tau^{j}_{y}\right)$
(7)

with

$\displaystyle\tau_{\gamma}=\tau_{z}\cos\left[2\theta_{\gamma}\right]+\tau_{x}\sin\left[2\theta_{\gamma}\right].$

Detailed derivations are presented in Appendix B. The ferro-orbital exchange
$J=-(t_{\sigma}-t_{\pi})^{2}/16U$ in Eq. (7) is inherently anisotropic
originating from the anisotropic shape of $p$ orbitals. In contrast, the
exchange $J^{\prime}=-t_{\sigma}t_{\pi}/8U$ is antiferro-orbital due to the
opposite sign of $t_{\sigma}$ and $t_{\pi}$. The first term in Eq. (7)
involving interacting orbital degrees of freedom is coined as the compass
model [44, 45]. The ground state is first studied by treating the orbital
pseudospin $\bm{\tau}$ as a classical vector. The orbital interaction can be
minimized via the diagonalization of the orbital exchange Hamiltonian in
momentum space
$H_{\text{OE}}=\sum_{\bm{k}}\bm{\tau}^{-\bm{k}}\Lambda_{\bm{k}}\bm{\tau}^{\bm{k}}$.
As shown in Fig. 4(a), the lowest eigenvalue of $\Lambda_{\bm{k}}$ is found at
$\Gamma$ point in RBZ and has a twofold degeneracy, suggesting that the
classical ground state is ferro-orbital ordering. The degenerated eigenvalue
has important implications on the structure of orbital order. A close
inspection on the Hamiltonian in Eq. (7) reveals that the $y$ component of
orbital pseudospin $\tau_{y}$ is decoupled from the other two components
$\tau_{z,x}$, in which the lowest degenerate eigenvalue arises. It indicates
that the ordering of orbital pseudospin $\bm{\tau}$ lies in the $zx$ plane.
More importantly, this degeneracy renders a continuous $\text{SO}(2)$
rotational symmetry of orbital pseudospin
$(\tau_{z},\tau_{x})=\tau(\cos\theta,\sin\theta)$ as shown in Fig. 4(b). Note
that this symmetry restricted to the classical ground state is emergent and is
not an exact symmetry of the orbital exchange model, which is only invariant
under finite point group rotations. The orbital order of the classical ground
state evolves in the $zx$ plane without any energy cost, which makes the
system particularly susceptible to quantum fluctuations. Following Holstein-
Primakoff spin wave theory [46], the zero-point energy arises from the
correction of quantum fluctuations, and is studied as a function of the
rotation $\theta$ about the $y$ axis of orbital pseudospin. To the leading
order, the zero-point energy takes the form

$E_{\text{ZP}}(\theta)=\frac{1}{2N}\sum_{\bm{k}}\omega_{\bm{k}}\left(\theta\right)+6J+12J^{\prime},$
(8)

where $N$ is the number of sites, the orbital excitation
$\omega_{\bm{k}}(\theta)=2\sqrt{[\varphi_{\bm{k}}(\theta)+2\varphi^{\prime}_{\bm{k}}-6J-12J^{\prime}]^{2}-[\varphi_{\bm{k}}(\theta)+\varphi_{\bm{k}}^{\prime}]^{2}}$,
and the auxiliary functions
$\\{\varphi_{\bm{k}}(\theta),\varphi_{\bm{k}}^{\prime}\\}=\\{2J\sum_{\gamma}\sin^{2}[2\theta_{\gamma}+\theta],4J^{\prime}\\}\cos[\bm{k}\cdot\bm{a}_{\gamma}]$.
Detailed derivations are presented in Appendix C. The numerical evaluation of
zero-point energy is shown in Fig. 4(c). The orbital fluctuation lifts the
degeneracy protected by the continuous $\text{SO}(2)$ rotational symmetry, and
selects the quantum ground state at $\theta=\pi/6+\mathbb{Z}\pi/3$
($\mathbb{Z}$ is an integer). This mechanism is known as order by disorder in
frustrated spin systems [47, 48, 49, 50, 51].

##  VI Conclusion and discussion

To summarize, we have studied the evolution of orbital ordering across the SF-
MI transition in the $p$-band triangular lattice. The ground-state phase
diagram is first established by Gutzwiller approach, which interpolates
continuously between two extreme limits, deep in SF phase and deep in MI
phase. The orbital orders in these two limits are further examined by the
perturbation approaches. With systematic analyses, we identify an intermediate
SF phase with the detwined in-plane ferro-orbital order, which correctly
reproduce the orbital order in the MI $n=1$ phase. It is worth remarking
several directions for further studies. The quantum fluctuations, which can be
partially restored with the cluster Gutzwiller approach [52, 53], may deserve
to be studied for its role in selecting the orbital order in the vicinity of
SF-MI transition. Alternatively, it is also interesting to investigate the
details of SF-MI transition within a single unified method, e.g. the quantum
Monte-Carlo simulation. Finally, we close by briefly discussing the
dissipation. The experimental realization of the Bose-Hubbard model in
photonic systems involves light-matter interactions, which may be better
described as an open system. Therefore, another direction to generalize our
work is to study the effect of dissipation.

## Acknowlegdgement

We thank Congjun Wu and W. Vincent Liu for helpful discussions. This work is
supported by the National Natural Science Foundation of China under Grants No.
11704338, No. 11534001, No. 11504008, and the National Basic Research Program
of China under Grant No. 2015CB921102.

##  Appendix A Details of Bogoliubov approach

In the weak-coupling limit, the band minima are occupied by a macroscopic
number of bosons $N_{0}$ at zero crystal momentum $\bm{k}=0$ in RBZ, which can
be expressed in terms of the ground-state wave function
$N_{0}=\sum_{\alpha\ell}|\phi_{\alpha\ell}|^{2}$. Here $\alpha$ and $\ell$
specify the orbital and sublattice degree of freedom, respectively. The
quantum effects arising from the commutation relation
$[p_{\alpha\ell\bm{k}=0},p^{\dagger}_{\alpha^{\prime}\ell^{\prime}\bm{k}=0}]=\delta_{\alpha,\alpha^{\prime}}\delta_{\ell,\ell^{\prime}}$
are suppressed by the macroscopic occupation $N_{0}$. The occupation of the
excited states at non-zero momentum $\bm{k}\neq 0$ with the corresponding
operator $p_{\alpha\ell\bm{k}\neq 0}$ is generally small and is treated as
quantum fluctuations $\tilde{p}_{\alpha\ell\bm{k}}$. The Hamiltonian can be
reexpressed by splitting the operators into the ground-state wave function
$\phi_{\alpha\ell}$ and the quantum fluctuations
$\tilde{p}_{\alpha\ell\bm{k}}$. To the lowest order in
$\tilde{p}_{\alpha\ell\bm{k}}$, one neglect all fluctuations and obtain the
energy functional per unit cell

$\displaystyle\varepsilon\left(\bm{\phi}^{*},\bm{\phi}\right)$
$\displaystyle=$
$\displaystyle\mathcal{T}_{\Gamma}\left(\bm{\phi}^{*},\bm{\phi}\right)+\mathcal{I}\left(\bm{\phi}^{*},\bm{\phi}\right)$
(9)

with the energy functional

$\displaystyle\mathcal{T}_{\bm{k}}\left(\bm{\phi}^{*},\bm{\phi}\right)=\frac{1}{N_{\text{uc}}}\sum_{\alpha\alpha^{\prime}}\sum_{\ell\ell^{\prime}}t^{\alpha\alpha^{\prime}}_{\ell\ell^{\prime}}(\bm{k})\phi_{\alpha\ell}^{*}\phi_{\alpha^{\prime}\ell^{\prime}}$
(10)

and

$\displaystyle\mathcal{I}\left(\bm{\phi}^{*},\bm{\phi}\right)$
$\displaystyle=$
$\displaystyle\frac{3U}{2N_{\text{uc}}^{2}}\sum_{\ell}\left[\right.|\phi_{x\ell}|^{4}+|\phi_{y\ell}|^{4}+\frac{4}{3}|\phi_{x\ell}|^{2}|\phi_{y\ell}|^{2}$
(11) $\displaystyle+$
$\displaystyle\frac{1}{3}\phi_{x\ell}^{*2}\phi_{y\ell}^{2}+\frac{1}{3}\phi_{x\ell}^{2}\phi_{y\ell}^{*2}\left.\right]$

arising from the hopping processes and the Hubbard interaction, respectively.
Here $N_{\text{uc}}$ is the number of unit cells and
$t^{\alpha\alpha^{\prime}}_{\ell\ell^{\prime}}(\bm{k})$ is the matrix elements
describing the process that the bosons with crystal momentum $\bm{k}$ hop
among the sublattices in the Wigner-Seitz cell. The ground-state wave function
is obtained by minimizing the energy functional in Eq. (9), which ensures the
linear order terms in $\tilde{\bm{p}}_{\alpha\ell\bm{k}}$ to vanish. The
quadratic order terms can be written as

$\displaystyle\mathcal{H}^{(2)}_{\bm{k}}=\frac{1}{2}\left[\tilde{\bm{p}}^{\dagger}_{\bm{k}},\tilde{\bm{p}}_{-\bm{k}}\right]\left[\begin{matrix}X_{\bm{k}}&Y\\\
Y^{\dagger}&X_{-\bm{k}}\end{matrix}\right]\left[\begin{matrix}\tilde{\bm{p}}_{\bm{k}}\\\
\tilde{\bm{p}}^{\dagger}_{-\bm{k}}\end{matrix}\right].$ (12)

The elements of diagonal matrix in the orbital-sublattice basis
$\tilde{\bm{p}}_{\bm{k}}=\left[\tilde{p}_{x1\bm{k}},\tilde{p}_{x2\bm{k}},\tilde{p}_{x3\bm{k}},\tilde{p}_{x4\bm{k}},\tilde{p}_{y1\bm{k}},\tilde{p}_{y2\bm{k}},\tilde{p}_{y3\bm{k}},\tilde{p}_{y4\bm{k}}\right]$
are

$\displaystyle X_{\bm{k}}^{\alpha\ell,\alpha^{\prime}\ell^{\prime}}$
$\displaystyle=$
$\displaystyle-(\mathcal{T}_{\Gamma}+2\mathcal{I})\delta_{\alpha,\alpha^{\prime}}\delta_{\ell,\ell^{\prime}}+t^{\alpha\alpha^{\prime}}_{\ell\ell^{\prime}}(\bm{k})$
(13) $\displaystyle+$
$\displaystyle\frac{2U}{N_{\text{uc}}}\delta_{\alpha,\alpha^{\prime}}\delta_{\ell,\ell^{\prime}}[|\phi_{x\ell}|^{2}+|\phi_{y\ell}|^{2}+2|\phi_{\alpha\ell}|^{2}]$
$\displaystyle+$
$\displaystyle\frac{2U}{N_{\text{uc}}}(1-\delta_{\alpha,\alpha^{\prime}})\delta_{\ell,\ell^{\prime}}[\phi_{x\ell}^{*}\phi_{y\ell}+\phi_{y\ell}^{*}\phi_{x\ell}].$

It is worth mentioning that the first term in Eq. (13) arises from the
conservation of bosons. The total bosons consist of condensed bosons at zero
momentum and excited bosons in the fluctuating fields
$N=N_{0}+\sum_{\alpha\ell\bm{k}}\tilde{p}_{\alpha\ell\bm{k}}^{\dagger}\tilde{p}_{\alpha\ell\bm{k}}$.
The elements of off-diagonal matrix are as follow

$\displaystyle Y^{\alpha\ell,\alpha^{\prime}\ell^{\prime}}$ $\displaystyle=$
$\displaystyle\frac{U}{N_{\text{uc}}}\delta_{\alpha,\alpha^{\prime}}\delta_{\ell,\ell^{\prime}}[\phi_{x\ell}^{2}+\phi_{y\ell}^{2}+2\phi_{\alpha\ell}^{2}]$
(14) $\displaystyle+$
$\displaystyle\frac{2U}{N_{\text{uc}}}(1-\delta_{\alpha,\alpha^{\prime}})\delta_{\ell,\ell^{\prime}}\phi_{x\ell}\phi_{y\ell}.$

The Bogoliubov excitation is determined by the eigenvalues of the Bogoliubov
dispersion matrix [54]

$\displaystyle\bm{\sigma}_{z}\left[\begin{matrix}X_{\bm{k}}&Y\\\
Y^{\dagger}&X_{-\bm{k}}\end{matrix}\right]=\left[\begin{matrix}X_{\bm{k}}&Y\\\
-Y^{\dagger}&-X_{-\bm{k}}\end{matrix}\right],$ (15)

where $\bm{\sigma}_{z}=\text{diag}\left(+\bm{1},-\bm{1}\right)$ is the block
Pauli matrix.

##  Appendix B Derivation of orbital exchange Hamiltonian

Table 1:  The $i$-th eigenstate $\Gamma_{n}^{i}$ of the Hamiltonian $H_{\text{I}}$ in Eq. (16) with the corresponding eigenenergy $E_{\Gamma_{n}^{i}}$ for $p^{n=0,1,2}$ configurations. |  $p^{0}$ configuration |  $p^{1}$ configuration |  $p^{2}$ configuration
---|---|---|---
$i$ | $1$ | $1$ | $2$ | $1$ | $2$ | $3$
$E_{\Gamma_{n}^{i}}$ | $0$ | $0$ | $0$ | $2U$ | $4U$ | $2U$
$\left|n_{+}n_{-}\right\rangle$ | $\left|0,0\right\rangle$ | $\left|1,0\right\rangle$ | $\left|0,1\right\rangle$ | $\left|2,0\right\rangle$ | $\left|1,1\right\rangle$ | $\left|0,2\right\rangle$

To derive the effective low-energy Hamiltonian, we shall first diagonalize the
on-site Hubbard interaction

$\displaystyle H_{\text{I}}$ $\displaystyle=$
$\displaystyle\frac{3}{2}U\left[\hat{n}\left(\hat{n}-\frac{2}{3}\right)-\frac{1}{3}\hat{L}_{z}^{2}\right].$
(16)

It is easy to show that the Hamiltonian in Eq. (16) commutes with both the
total occupation operator $\hat{n}$ and the $z$-component angular momentum
$\hat{L}_{z}$. Moreover, the matrix representation of the operators
$\hat{n}=\hat{n}_{+}+\hat{n}_{-}$ and $\hat{L}_{z}=\hat{n}_{+}-\hat{n}_{-}$ is
diagonal in the basis of axial orbitals $p_{\pm}=p_{x}\pm ip_{y}$. Here,
$\hat{n}_{\pm}$ correspond to the occupation operators of axial orbitals
$p_{\pm}$. Therefore, the eigenstate of the Hamiltonian in Eq. (16) can be
labelled by the corresponding quantum number

$\displaystyle\left|n_{+}n_{-}\right\rangle=\frac{1}{\sqrt{n_{+}!n_{-}!}}\left(p_{+}^{\dagger}\right)^{n_{+}}\left(p_{-}^{\dagger}\right)^{n_{-}}\left|0,0\right\rangle.$
(17)

The orbital exchange model in the MI $n=1$ phase involves the $p^{n=0,1,2}$
configurations. The eigenstates $\Gamma_{n}^{i}$ with the corresponding
eigenenergies $E_{\Gamma_{n}^{i}}$ are listed in Table. 1. The $p^{1}$
configuration with zero energy is an orbital doublet with one boson occupying
either $p_{+}$ or $p_{-}$ orbital. Note that the charge excitation
$\left(p^{1}\right)_{i}\left(p^{1}\right)_{j}\rightleftharpoons\left(p^{2}\right)_{i}\left(p^{0}\right)_{j}$
through the hopping processes $t_{\alpha\beta}d^{\dagger}_{i\alpha}d_{j\beta}$
has an energy gap that is proportional to the Hubbard interaction $U$. In the
large-$U$ limit, the effective low-energy model is described by the second-
order hopping process with both the initial and final states in $p^{1}$
configuration, which involves no charge gap. Let us first derive the orbital
exchange interaction along the $\bm{a}_{1}=\hat{\bm{x}}$ bonds. Employing the
second-order perturbation theory [55], the matrix form of orbital exchange
interaction is given by

$\displaystyle\mathbb{J}_{kl,k^{\prime}l^{\prime}}=$ $\displaystyle-$
$\displaystyle\sum_{mn}\frac{1}{E_{\Gamma_{2}^{m}}+E_{\Gamma_{0}^{n}}}$ (18)
$\displaystyle\times$ $\displaystyle\left\langle\underset{i\text{-th
site}}{\Gamma_{1}^{k}}\underset{j\text{-th
site}}{\Gamma_{1}^{l}}\right|\sum_{\alpha\beta}t_{\alpha\beta}^{*}p^{\dagger}_{j\beta}p_{i\alpha}\left|\underset{i\text{-th
site}}{\Gamma_{2}^{m}}\underset{j\text{-th
site}}{\Gamma_{0}^{n}}\right\rangle$ $\displaystyle\times$
$\displaystyle\left\langle\underset{i\text{-th
site}}{\Gamma_{2}^{m}}\underset{j\text{-th
site}}{\Gamma_{0}^{n}}\right|\sum_{\delta\gamma}t_{\delta\gamma}p_{i\alpha}^{\dagger}p_{j\gamma}\left|\underset{i\text{-th
site}}{\Gamma_{1}^{k^{\prime}}}\underset{j\text{-th
site}}{\Gamma_{1}^{l^{\prime}}}\right\rangle$ $\displaystyle+$
$\displaystyle{i\leftrightarrow j}.$

To describe the orbital exchange Hamiltonian, we introduce the orbital
pseudospin operators
$\\{\tau_{+},\tau_{-}\\}\equiv\\{p^{\dagger}_{x}p_{y},p^{\dagger}_{y}p_{x}\\}$,
which flip the states of the orbital doublet in $p^{1}$ configuration. The $z$
component of pseudospin $\bm{\tau}$-vector follows through the spin-$1/2$
angular momentum algebra $\tau_{z}=\left[\tau_{+},\tau_{-}\right]$. A lengthy
but straightforward algebra leads the orbital exchange Hamiltonian along the
$\bm{a}_{1}$ bonds

$H_{\text{OE}}^{1}=\sum_{\langle
ij\rangle}\left[J\tau_{z}^{i}\tau_{z}^{j}+J^{\prime}\left(\bm{\tau}^{i}\cdot\bm{\tau}^{j}+2\tau^{i}_{y}\tau^{j}_{y}\right)\right]$
(19)

with

$\displaystyle\\{J,J^{\prime}\\}$ $\displaystyle=$
$\displaystyle\\{-\frac{\left(t_{\sigma}-t_{\pi}\right)^{2}}{16U},\frac{t_{\sigma}t_{\pi}}{8U}\\}.$
(20)

Having derived the orbital exchange interaction $H_{\text{OE}}^{1}$ along bond
vector $\bm{a}_{1}$, the interaction $H_{\text{OE}}^{2,3}$ has exactly the
same form with $H_{\text{OE}}^{1}$ if the orbital pseudospin operators
$\bm{\tau}$ are defined in the local coordinate. Thus, the connection between
the local and global coordinates (the global $x$ axis along $\bm{a}_{1}$ bond
vector) is linked by a rotation of $\theta=\frac{2\pi}{3},\frac{4\pi}{3}$
about $z$ axis, corresponding to the $\bm{a}_{2},\bm{a}_{3}$ bonds,
respectively. Under this rotation, the $p$ orbital wave functions transform as

$\displaystyle p_{x}$ $\displaystyle\to$ $\displaystyle\cos\theta
p_{x}-\sin\theta p_{y},$ (21a) $\displaystyle p_{y}$ $\displaystyle\to$
$\displaystyle\sin\theta p_{x}+\cos\theta p_{y}.$ (21b)

Accordingly, the pseudospin operators $\bm{\tau}$ transform as follow

$\displaystyle\tau_{z}$ $\displaystyle\to$
$\displaystyle\sin\left[2\theta\right]\tau_{x}+\cos\left[2\theta\right]\tau_{z},$
(22a) $\displaystyle\tau_{x}$ $\displaystyle\to$
$\displaystyle\cos\left[2\theta\right]\tau_{x}-\sin\left[2\theta\right]\tau_{z},$
(22b) $\displaystyle\tau_{y}$ $\displaystyle\to$ $\displaystyle\tau_{y}.$
(22c)

The pseudospin vector $\bm{\tau}$ is rotated by $2\theta$ about its $y$ axis
in the pseudospin space. It is now straightforward to obtain the Hamiltonian
$H_{\text{OE}}^{2,3}$ by replacing the pseudospin $\bm{\tau}$ in
$H_{\text{OE}}^{1}$. Finally, the total orbital exchange Hamiltonian takes the
form

$\displaystyle H_{\text{OE}}=J\sum_{\langle
ij\rangle\parallel\bm{a}_{\gamma}}\tau^{i}_{\gamma}\tau^{j}_{\gamma}+J^{\prime}\sum_{\langle
ij\rangle}\left(\bm{\tau}^{i}\cdot\bm{\tau}^{j}+2\tau^{i}_{y}\tau^{j}_{y}\right)$
(23)

with

$\displaystyle\tau_{\gamma}=\tau_{z}\cos\left[2\theta_{\gamma}\right]+\tau_{x}\sin\left[2\theta_{\gamma}\right],\\{\theta_{1},\theta_{2},\theta_{3}\\}=\\{0,\frac{2}{3}\pi,\frac{4}{3}\pi\\}.$

Unlike spin systems, the orbital exchange is anisotropic. It roots in the fact
that the hopping processes are bond dependent due to the spatial orientation
of $p$ orbitals.

##  Appendix C Linear orbital wave

Having established the classical ground state of the orbital exchange
Hamiltonian $H_{\text{OE}}$ in Eq. (23), we then proceed to derive the
Hamiltonian that describes the orbital excitation. The orbital pseudospin
operators, obeying the angular momentum algebra of spin $T=1/2$, can be
expressed in terms of Holstein-Primakoff bosons [46]

$\displaystyle\tau_{x}$ $\displaystyle=$
$\displaystyle\sqrt{2T}\left(a^{\dagger}+a\right),$ (25a)
$\displaystyle\tau_{y}$ $\displaystyle=$ $\displaystyle
i\sqrt{2T}\left(a^{\dagger}-a\right),$ (25b) $\displaystyle\tau_{z}$
$\displaystyle=$ $\displaystyle 2\left(T-a^{\dagger}a\right).$ (25c)

Note that the classical ground state enjoys a continuous SO$(2)$ rotational
symmetry, which transforms the orbital pseudospin as

$\displaystyle\tau_{z}$ $\displaystyle\to$
$\displaystyle\cos\theta\tau_{z}+\sin\theta\tau_{x},$ (26a)
$\displaystyle\tau_{x}$ $\displaystyle\to$
$\displaystyle-\sin\theta\tau_{z}+\cos\theta\tau_{x}.$ (26b)

The pseudospin in Eq. (23) is first replaced by the above transformation.
Expanding in powers of $T$ followed by Fourier transformation then leads to
the following $\theta$-dependent Hamiltonian

$\displaystyle
H_{\text{OE}}(\theta)=4T^{2}NE_{\text{c}}+2T\sum_{\bm{k}}\hat{\mathcal{H}}_{\text{LOW}}\left(\bm{k},\theta\right)+\mathcal{O}(\sqrt{T})\hskip
22.76219pt$ (27)

The first term $E_{\text{c}}=3J+6J^{\prime}$ in Eq. (27) recovers the
classical ground energy per site. The second term in Eq. (27) describes the
linear orbital wave Hamiltonian

$\displaystyle\hat{\mathcal{H}}_{\text{LOW}}\left(\bm{k},\theta\right)$
$\displaystyle=$
$\displaystyle-\left(12J+24J^{\prime}\right)a^{\dagger}_{\bm{k}}a_{\bm{k}}-3\varphi_{\bm{k}}^{\prime}\left(a_{\bm{k}}a_{-\bm{k}}+a_{-\bm{k}}^{\dagger}a_{\bm{k}}^{\dagger}\right)+\left[\varphi_{\bm{k}}\left(\theta\right)+2\varphi_{\bm{k}}^{\prime}\right]\left(a_{-\bm{k}}^{\dagger}+a_{\bm{k}}\right)\left(a_{\bm{k}}^{\dagger}+a_{-\bm{k}}\right)$
(28) $\displaystyle=$
$\displaystyle\left[a^{\dagger}_{\bm{k}},a_{-\bm{k}}\right]\left[\begin{matrix}\varphi_{\bm{k}}(\theta)+2\varphi^{\prime}_{\bm{k}}-6J-12J^{\prime}&\varphi_{\bm{k}}(\theta)-\varphi^{\prime}_{\bm{k}}\\\
\varphi_{\bm{k}}(\theta)-\varphi^{\prime}_{\bm{k}}&\varphi_{\bm{k}}(\theta)+2\varphi^{\prime}_{\bm{k}}-6J-12J^{\prime}\end{matrix}\right]\left[\begin{matrix}a_{\bm{k}}\\\
a^{\dagger}_{-\bm{k}}\end{matrix}\right]+N\left(6J+12J^{\prime}\right)$ (29)

with the auxiliary functions

$\displaystyle\varphi_{\bm{k}}\left(\theta\right)$ $\displaystyle=$
$\displaystyle
2J\sum_{\gamma}\sin^{2}[2\theta_{\gamma}+\theta]\cos[\bm{k}\cdot\bm{a}_{\gamma}],$
(30a) $\displaystyle\varphi_{\bm{k}}^{\prime}$ $\displaystyle=$ $\displaystyle
4J^{\prime}\cos[\bm{k}\cdot\bm{a}_{\gamma}].$ (30b)

It can be diagonalized via the Bogoliubov transformation [54]

$\displaystyle\left[\begin{matrix}a_{\bm{k}}\\\
a^{\dagger}_{-\bm{k}}\end{matrix}\right]=T_{\bm{k}}\left[\begin{matrix}b_{\bm{k}}\\\
b^{\dagger}_{-\bm{k}}\end{matrix}\right],$ (31)

which relates the Holstein-Primakoff bosons with orbital excitation modes.
Notably, the bosonic statistics require that $T_{\bm{k}}$ satisfies the para-
unitary condition

$\displaystyle T^{\dagger}_{\bm{k}}\sigma_{z}T_{\bm{k}}=\sigma_{z}.$ (32)

Therefore, the linear orbital wave Hamiltonian can be written in terms of
diagonalized bosons

$\displaystyle
H_{\text{LOW}}\left(\theta\right)=\sum_{\bm{k}}\omega_{\bm{k}}(\theta)b^{\dagger}_{\bm{k}}b_{\bm{k}}+NE_{\text{ZP}}\left(\theta\right)$
(33)

where the orbital wave excitation is

$\displaystyle\omega_{\bm{k}}(\theta)=2\sqrt{[\varphi_{\bm{k}}(\theta)+2\varphi^{\prime}_{\bm{k}}-6J-12J^{\prime}]^{2}-[\varphi_{\bm{k}}(\theta)+\varphi_{\bm{k}}^{\prime}]^{2}}$

and the energy from zero-point motion per site is

$\displaystyle
E_{\text{ZP}}\left(\theta\right)=\frac{1}{2N}\sum_{\bm{k}}\omega_{\bm{k}}(\theta)+6J+12J^{\prime}.$
(34)

Finally, the zero-point energy can be evaluated numerically.

## References

  * Tokura and Nagaosa [2000] Y. Tokura and N. Nagaosa, Science 288, 462 (2000).
  * Maekawa  _et al._ [2004] S. Maekawa, T. Tohyama, S. E. Barnes, S. Ishihara, W. Koshibae, and G. Khaliullin,  _Physics of Transition Metal Oxides_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2004).
  * Khomskii [2014] D. I. Khomskii,  _Transition Metal Compounds_ (Cambridge University Press, 2014).
  * Fernandes  _et al._ [2014] R. M. Fernandes, A. V. Chubukov, and J. Schmalian, Nature Physics 10, 97 (2014).
  * Fradkin  _et al._ [2010] E. Fradkin, S. A. Kivelson, M. J. Lawler, J. P. Eisenstein, and A. P. Mackenzie, Annual Review of Condensed Matter Physics 1, 153 (2010).
  * Fernandes  _et al._ [2019] R. M. Fernandes, P. P. Orth, and J. Schmalian, Annual Review of Condensed Matter Physics 10, 133 (2019).
  * Yu  _et al._ [2018] R. Yu, J.-X. Zhu, and Q. Si, Phys. Rev. Lett. 121, 227003 (2018).
  * Si  _et al._ [2009] Q. Si, E. Abrahams, J. Dai, and J.-X. Zhu, New Journal of Physics 11, 045001 (2009).
  * Si [2009] Q. Si, Nature Physics 5, 629 (2009).
  * Wu [2009] C. Wu, Modern Physics Letters B 23, 1 (2009).
  * Li and Liu [2016] X. Li and W. V. Liu, Reports on Progress in Physics 79, 116401 (2016).
  * Müller  _et al._ [2007] T. Müller, S. Fölling, A. Widera, and I. Bloch, Phys. Rev. Lett. 99, 200405 (2007).
  * Wirth  _et al._ [2011] G. Wirth, M. Ölschläger, and A. Hemmerich, Nature Physics 7, 147 (2011).
  * Soltan-Panahi  _et al._ [2012] P. Soltan-Panahi, D.-S. Luhmann, J. Struck, P. Windpassinger, and K. Sengstock, Nature Physics 8, 71 (2012).
  * Kock  _et al._ [2016] T. Kock, C. Hippler, A. Ewerbeck, and A. Hemmerich, Journal of Physics B 49, 042001 (2016).
  * Niu  _et al._ [2018] L. Niu, S. Jin, X. Chen, X. Li, and X. Zhou, Phys. Rev. Lett. 121, 265301 (2018).
  * Jacqmin  _et al._ [2014] T. Jacqmin, I. Carusotto, I. Sagnes, M. Abbarchi, D. D. Solnyshkov, G. Malpuech, E. Galopin, A. Lemaître, J. Bloch, and A. Amo, Phys. Rev. Lett. 112, 116402 (2014).
  * Milicevic  _et al._ [2017] M. Milicevic, T. Ozawa, G. Montambaux, I. Carusotto, E. Galopin, A. Lemaître, L. Le Gratiet, I. Sagnes, J. Bloch, and A. Amo, Phys. Rev. Lett. 118, 107403 (2017).
  * Milicevic  _et al._ [2019] M. Milicevic, G. Montambaux, T. Ozawa, O. Jamadi, B. Real, I. Sagnes, A. Lemaître, L. Le Gratiet, A. Harouri, J. Bloch, and A. Amo, Phys. Rev. X 9, 031010 (2019).
  * Wu and Das Sarma [2008] C. Wu and S. Das Sarma, Phys. Rev. B 77, 235107 (2008).
  * Castro Neto  _et al._ [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * Novoselov [2011] K. S. Novoselov, Rev. Mod. Phys. 83, 837 (2011).
  * Geim [2011] A. K. Geim, Rev. Mod. Phys. 83, 851 (2011).
  * Jin  _et al._ [2021] S. Jin, W. Zhang, X. Guo, X. Chen, X. Zhou, and X. Li, Phys. Rev. Lett. 126, 035301 (2021).
  * Haldane [1990] F. D. M. Haldane, The hierarchy of fractional states and numerical studies, in  _The Quantum Hall Effect_, edited by R. E. Prange and S. M. Girvin (Springer New York, New York, NY, 1990) pp. 303–352.
  * Chen and Xie [2018] H. Chen and X. C. Xie, Phys. Rev. A 98, 053611 (2018).
  * Isacsson and Girvin [2005] A. Isacsson and S. M. Girvin, Phys. Rev. A 72, 053604 (2005).
  * Liu and Wu [2006] W. V. Liu and C. Wu, Phys. Rev. A 74, 013607 (2006).
  * Umucalılar and Oktel [2008] R. O. Umucalılar and M. O. Oktel, Phys. Rev. A 78, 033602 (2008).
  * Chin  _et al._ [2010] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
  * Carusotto and Ciuti [2013] I. Carusotto and C. Ciuti, Rev. Mod. Phys. 85, 299 (2013).
  * Gutzwiller [1963] M. C. Gutzwiller, Phys. Rev. Lett. 10, 159 (1963).
  * Gutzwiller [1965] M. C. Gutzwiller, Phys. Rev. 137, A1726 (1965).
  * Rokhsar and Kotliar [1991] D. S. Rokhsar and B. G. Kotliar, Phys. Rev. B 44, 10328 (1991).
  * Krauth  _et al._ [1992] W. Krauth, M. Caffarel, and J.-P. Bouchaud, Phys. Rev. B 45, 3137 (1992).
  * Bogolyubov [1947] N. Bogolyubov, J. Phys. (USSR) 11, 23 (1947).
  * Abrikosov  _et al._ [1963] A. A. Abrikosov, L. P. Gorkov, and I. E. Dzyaloshinski,  _Methods of Quantum Field Theory in Statistical Physics_ (Dover, 1963).
  * Pethick and Smith [2008] C. J. Pethick and H. Smith,  _Bose-Einstein Condensation in Dilute Gases_, 2nd ed. (Cambridge University Press, 2008).
  * Dalfovo  _et al._ [1999] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 71, 463 (1999).
  * Wu  _et al._ [2006] C. Wu, W. V. Liu, J. Moore, and S. Das Sarma, Phys. Rev. Lett. 97, 190406 (2006).
  * Xu  _et al._ [2016] Z.-F. Xu, L. You, A. Hemmerich, and W. V. Liu, Phys. Rev. Lett. 117, 085301 (2016).
  * Di Liberto  _et al._ [2016] M. Di Liberto, A. Hemmerich, and C. Morais Smith, Phys. Rev. Lett. 117, 163001 (2016).
  * Zhou  _et al._ [2017] S. Zhou, K. Jiang, H. Chen, and Z. Wang, Phys. Rev. X 7, 041018 (2017).
  * Kugel and Khomskii [1982] K. I. Kugel and D. I. Khomskii, Soviet Physics Uspekhi 25, 231 (1982).
  * Nussinov and van den Brink [2015] Z. Nussinov and J. van den Brink, Rev. Mod. Phys. 87, 1 (2015).
  * Holstein and Primakoff [1940] T. Holstein and H. Primakoff, Phys. Rev. 58, 1098 (1940).
  * Villain  _et al._ [1980] J. Villain, R. Bidaux1, J.-P. Carton, and R. Conte, J. Phys. France 41, 1263 (1980).
  * Henley [1989] C. L. Henley, Phys. Rev. Lett. 62, 2056 (1989).
  * Lacroix  _et al._ [2011] C. Lacroix, P. Mendels, and F. Mila,  _Introduction to Frustrated Magnetism_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2011).
  * Diep [2013] H. T. Diep,  _Frustrated Spin Systems_, 2nd ed. (WORLD SCIENTIFIC, 2013).
  * Green  _et al._ [2018] A. G. Green, G. Conduit, and F. Krüger, Annual Review of Condensed Matter Physics 9, 59 (2018).
  * Lühmann [2013] D.-S. Lühmann, Phys. Rev. A 87, 043619 (2013).
  * Bai  _et al._ [2018] R. Bai, S. Bandyopadhyay, S. Pal, K. Suthar, and D. Angom, Phys. Rev. A 98, 023606 (2018).
  * Blaizot and Ripka [1986] J. P. Blaizot and G. Ripka,  _Quantum Theory of Finite Systems_ (MIT Press, 1986).
  * Kuklov and Svistunov [2003] A. B. Kuklov and B. V. Svistunov, Phys. Rev. Lett. 90, 100401 (2003).


.tgapng.pngconvert #1 .tga

††thanks: Lead Contact

# Aggregation and structural phase transitions of semiflexible polymer
bundles: a braided circuit topology approach

Jonas Berx  Medical Systems Biophysics and Bioengineering, Leiden Academic
Centre for Drug Research, Faculty of Science, Leiden University, 2333CC
Leiden, The Netherlands  Alireza Mashaghi
<EMAIL_ADDRESS>Medical Systems Biophysics and
Bioengineering, Leiden Academic Centre for Drug Research, Leiden University,
2333CC Leiden, The Netherlands

###### Abstract

We present a braided circuit topology framework for investigating topology and
structural phase transitions in aggregates of semiflexible polymers. In the
conventional approach to circuit topology, which specifically applies to
single isolated folded linear chains, the number and arrangement of contacts
within the circuitry of a folded chain give rise to increasingly complex fold
topologies. Another avenue for achieving complexity is through the interaction
and entanglement of two or more folded linear chains. The braided circuit
topology approach describes the topology of such multiple-chain systems and
offers topological measures such as writhe, complexity, braid length, and
isotopy class. This extension of circuit topology to multichains reveals the
interplay between collapse, aggregation, and entanglement. In this work, we
show that circuit topological motif fractions are ideally suited order
parameters to characterise structural phase transitions in entangled systems
that can detect structural re-ordering other measures cannot.

Circuit topology, conformational phase transitions, entanglement, braiding

## Introduction

It has been shown that, for a single semiflexible chain, the polymer stiffness
leads to a structural “phase diagram”, displaying a multitude of of
conformations such as globular, hairpin, knotted or extended structures, both
for off-lattice [1, 2, 3] and on-lattice systems [4, 5]. By keeping a constant
monomer-monomer interaction energy scale and decreasing the temperature, the
polymer minimizes its energy by either collapsing or stiffening locally.
Adjusting the stiffness $\kappa$ allows one to study the competition between
these two effects, which leads to fundamentally distinct structural motifs.

This line of reasoning can be extended to aggregates of semiflexible polymers.
It has been shown [6, 2] that for small systems consisting of $M=2,4,8$
polymers of length $N=13$ at a high temperature $T$, the system is fragmented
and individual polymers can be considered isolated. The structural properties
in this regime follow the single-chain results in good approximation. For a
decreasing temperature, however, flexible polymers aggregate in an amorphous
globular configuration, while stiffer polymers form (twisted) bundles. The
concomitant disentanglement of long polymers during cooling also depends
fundamentally on the stiffness [7], the study of which offers insights into
effects such as the asymmetry between crystallisation and melting, or glass
formation [8]. Moreover, for polymer melts lowering the temperature can lead
to a transition from a disordered melt phase to an ordered nematic and
crystalline phase, as a consequence of increasing stiffness [9, 10].

The collapse and aggregation transitions are not separate processes. If the
energy reduction associated with the formation of more interchain contacts
(aggregation) is more favourable than the formation of intrachain contacts
(collapse or folding), multichain aggregation may be expected to undo single-
chain collapse. Hence, by studying the inter-, and intrachain contacts, we can
deduce information about the structural properties and phase transitions of
the system. Both types of contacts can be neatly studied within the framework
of circuit topology (CT), which was originally introduced to describe the
topology of folded proteins and folding reactions [11, 12], and was
subsequently extended to include entanglement (i.e., ‘soft contacts’) [13, 14,
15] and multiple chains [16]. By combining circuit topology with a topological
description on the level of entanglement, we get a more complete picture of
aggregation processes with a small number of chains.

The outline of this paper is as follows: in section I we set the stage for our
analysis by reiterating the basic concepts of circuit topology for multichain
systems. Subsequently, in section II, we discuss the theoretical framework
necessary for our concomitant braid-theoretic description. Section III
investigates the structural phase transitions in a system of $M=4$ chains of
varying monomer number $N\in\\{10,\,30\\}$ and connects this with the circuit
topology and braiding analysis. Finally, in section IV, we present conclusions
and a future outlook.

##  I Multichain Circuit Topology

We start by revisiting multichain circuit topology. To distinguish between the
aforementioned ‘soft contacts’, which are units of single-chain entanglement,
we will henceforth only consider ‘hard contacts’, which we will hereafter
simply refer to as ‘contacts’. To model the formation of intra- and inter-
chain bonds in molecular systems and the resulting topology, we can examine
the mutual topological relation of binary contacts that represent these bonds
in the circuit topology framework. Two different contacts, that we name $A$
and $B$, each consisting of two contact ‘sites’ bearing the same name -
totalling four sites - will be assigned to a topological relation. Note that
such ‘bonds’ pertain to any types of atomic or molecular interactions that
contribute to folding and assembly, including covalent and non-covalent ones
(e.g., hydrogen bonds, disulfide bonds, beta-beta interactions in
biomolecules). Depending on the application, one may focus on a specific type
of bonds, and model them as binary contacts.

On a single strand, only three possible arrangements of $A$ and $B$ are
possible: $AABB$, $ABBA$ and $ABAB$, corresponding respectively to a series
($S$), parallel ($P$) or cross ($X$) motif, where we take renaming
$(A\leftrightarrow B)$ into account.

Similar to the case of single chains, a motif that comprises pairs of contact
sites is referred to as a “circuit”. A circuit can be braided, as described in
subsequent sections. When considering a system consisting of multiple distinct
open-ended chains, a number of motifs needs to be added to the above set of
circuits. For $n=2$ chains, we can discern three topologically distinct
motifs: independent ($I_{2}$), loop ($L_{2}$) and tandem ($T_{2}$). Within
these motifs, degeneracies are possible. The independent relation is unique up
to renaming of the contact sites, the loop relation can form either a
“parallel” or a “cross” loop (not to be confused with the single-chain CT
motifs), and the tandem relation occurs in either an “umbrella” or an “arc”
motif (Fig. 1).

For $n=3$, only independent ($I_{3}$) and tandem ($T_{3}$) circuits can be
formed, which are both non-degenerate. Lastly, for $n=4$, only the independent
relation ($I_{4}$) can occur. Four is the maximum number of chains that can
participate in creating circuits with four contact sites. In Fig. 1, we list
all possible multichain motifs.

Extending now the string notation for the single-chain CT to the multichain
framework requires the introduction of “ghost contacts”, which we will denote
with $\mathcal{O}$ in the string. To construct the string notation
corresponding to a particular motif, we proceed as follows. The two chains are
placed parallel and aligned vertically. We now read from top to bottom and
from left to right, adding a ghost contact wherever a real contact site does
not possess partners on the same level on other chains.

A simple example will clarify our formalism. Consider the $T_{2}$ umbrella CT
motif, i.e., the tandem relation on two chains. We assume that the first chain
counted from the left possesses the ordered contact sites $A\,,B\,,A$, and
that the second chain only possesses a single $B$ contact site. We supplement
the right chain with two ghost nodes. The string describing $T_{2}$ is then
$\mathcal{S}_{2}=ABB\mathcal{O}A\mathcal{O}$, where the subscript indicates
the number of chains that participate in the motif. Note that the same motif
is obtained independent of the exact position of the ghost contacts on the
second chain.

The string notations for all multichain motifs can be found in Fig. 1. Note
that all motifs within a degenerate circuit class (e.g., $T_{2}$) can be
obtained by cyclically permuting the characters corresponding to a single
chain.

Figure 1: The set of circuit topology motifs for $n=1,2,3,4$ strands. Contacts
$A$ and $B$ are indicated by red and turquoise filled circles, respectively,
while ghost contacts $\mathcal{O}$ are indicated by filled green circles.
Different chains are indicated by different colours. The coloured chain ends
indicate the orientation of the strands; we choose to orient every strand from
the yellow to the violet end. Corresponding string notations are given below
each motif, where the string is read in $n$-tuples.

##  II Topological entanglement and braiding measures

We now assume that the different chains can cross one another in the embedding
space, and henceforth call the number of chains the _braid index_ $n$, where
we will limit ourselves to $n\leq 4$. For $n=2$, a twist of the chains is not
topologically protected, i.e., it can be undone by rotating one of the planes
in which the endpoints are fixed. we will henceforth assume that the resulting
braid projection is confined to the plane spanned by the parallel lines
connecting the first and last contact-contact or contact-ghost pair. We label
the chains by $\nu=1\,,2\,,3\,,...$ from the left and introduce notation that
is common in the theory of braiding, by following the Artin representation
[17]. The operators $\sigma_{i}$, with $i=1,2,...,n-1$, indicate that the
$i$th chain, as counted from the left passes _over_ the $(i+1)$th chain. The
inverse operators $\sigma_{1}^{-1}\,,\sigma_{2}^{-1}\,,...\,,\sigma_{i}^{-1}$
indicate that the $i$th chain passes _underneath_ the $(i+1)$th chain. Note
that the index $i$ is independent of the labelling $\nu$ of the strands; $i$
represents only the order of the strands, as read from left to right, while
$\nu$ represents the strands themselves. This way, we can represent a braid
with index $n$ by a string of $\sigma-$operators. For example, the string
$\beta=\sigma_{2}\sigma_{1}^{-1}\sigma_{3}^{-1}\sigma_{1}\sigma_{1}\sigma_{2}\sigma_{2}^{-1}$
corresponds to the braid shown in Fig. 2(a).

Figure 2: (a) A braid with index $n=4$ with string notation
$\mathcal{S}_{4}=\sigma_{2}\sigma_{1}^{-1}\sigma_{3}^{-1}\sigma_{1}\sigma_{1}\sigma_{2}\sigma_{2}^{-1}$.
(b) A braid containing contacts $A,\,B$ indicated by red and turquoise points,
respectively. The string notation is
$\mathcal{S}_{3}=A\mathcal{O}A\,\sigma_{2}^{-1}\sigma_{1}^{-1}\sigma_{2}\,BB\mathcal{O}\,\sigma_{2}$.

We remark here that our string notation for CT and the Artin braid notation
can quite naturally be combined into a generic framework by inserting the
braid operators $\sigma_{i}$ in the string after every $n$-tuple of contacts.
An example is given in Fig. 2(b). Hard contacts can then be modelled within
this framework as separate operators, with their own set of Reidemeister moves
and bond moves [18]. We defer to ongoing theoretical research on this topic
and will not use the string notation any further in this manuscript.

Let us formulate some useful properties of braids before continuing. Similar
to knot theory, braids belong to the same equivalence class if they are
related by a set of moves that are related to the so-called Reidemeister
moves. The three fundamental moves can be formulated as follows:

  1. 1.

$\sigma_{i}\sigma_{j}=\sigma_{j}\sigma_{i}$ for $|i-j|\geq 2$ (disjoint strand
relation)

  2. 2.

$\sigma_{i}\sigma_{i+1}\sigma_{i}=\sigma_{i+1}\sigma_{i}\sigma_{i+1}$ (Skein
relation)

  3. 3.

$\sigma_{i}\sigma_{i}^{-1}=e$ (annihilation relation)

The element $e$ is the unit element of the braid group; if we encounter such
an operator pair we can eliminate it from the string. The braid group $B_{n}$
itself is then defined from these moves.

A natural, yet nontrivial question now arises: which topological quantities
can we compute from this braid representation to characterise the entanglement
of our system of semiflexible polymers? We briefly discuss four possible
measures that can quantify certain topological properties of a braid: the
writhe, braid length, complexity index and isotopy class.

###  II.1 Writhe

The writhe $\mathrm{W}(\beta)$ of a braid $\beta$ is defined as the sum of the
exponents of the braid operators, i.e., for
$\beta=\sigma_{i_{1}}^{a_{1}}\sigma_{i_{2}}^{a_{2}}...\sigma_{i_{m}}^{a_{m}}$
it is $\mathrm{W}(\beta)=\sum_{k}a_{k}$. It characterises the net twist of a
braid and can be used to detect whether underlying (chemical) chirality of the
monomers or interactions between them influences the mesoscopic twist of the
system [19]. When there are no interactions that prefer one handedness over
the other, the writhe can be expected to be symmetrically distributed with
mean zero. Since two equivalent braids have the same writhe, it is a
topological braid invariant. In, e.g., Fig. 2(a) the writhe is $\mathrm{W}=1$,
while for (b) is is $\mathrm{W}=0$.

Since in the next section we will use the Kremer-Grest model, which does not
intrinsically include any chiral interactions, to numerically simulate the
system, we compute the writhe only as a sanity check.

###  II.2 Braid length

The braid length $L(\beta)$ is a topological measure, not a geometrical one.
It counts the number of operators in the braid word $\beta$. Since we can
always add the trivial combination $\sigma_{i}\sigma_{i}^{-1}$ to the braid
word, the length is not an invariant. However, the minimal braid length
$L_{c}(\beta)$, obtained by reducing the braid word to its normal form, is a
braid invariant. Henceforth, we will only consider the minimal length and will
just refer to it as “the length $L$”, without ambiguity. In Fig. 2(a), the
minimum braid length is $L=3$, since the normal form of the braid is
$\beta_{c}=\sigma_{2}\sigma_{3}^{-1}\sigma_{1}$.

To some extent, the length can be used as a measure for the complexity of a
braid, i.e., the greater the length, the more complex the associated braid.
Since it is a feature of the normal form, trivial entanglement such as, e.g.,
projection-induced crossings $\sigma_{i}\sigma_{i}^{-1}$ are eliminated and
only actual crossings are counted. Naturally, since we only consider systems
with a finite number of monomers, we can assume that the maximum braid length
is bounded and that it heavily depends on the number of chains $M$, number of
monomers per chain $N$ and the stiffness $\kappa$.

###  II.3 Complexity index

Just how entangled is a braid? This is not an easy question to answer, since
complexity is not an objective observable. Therefore, we make the choice that
henceforth entanglement can be characterised by the complexity index
$C(\beta)$ [20]. The braid complexity index, together with the braid length,
constitutes our main measure to characterise the topology and entanglement
properties of a braid. It is defined as the number of intersections (which we
denote by #) of a curve spanning a punctured disk of $n$ punctures with the
real axis, after application of the braiding operations. This is illustrated
in Fig. 3. Normalising by the number of initial intersections $\\#E$ with the
real axis, which is $n-1=3$ in our case, and subsequently taking the natural
logarithm, the complexity can be written as $C(\beta)=\log\left(\frac{\\#\beta
E}{\\#E}\right)$.

Figure 3: The punctured disk representation of the braid complexity for four
different braids. Punctures corresponding to strands are indicated by coloured
dots. The operators $\sigma_{i},\,\sigma_{i}^{-1}$ change punctures $i$ and
$i+1$ CW or CCW, respectively. Intersections with the central axis are
indicated by open circles.

So are for example the braids $\beta_{1}=\sigma_{1}\sigma_{2}$ and
$\beta_{2}=\sigma_{1}\sigma_{2}^{-1}$ very similar, differing only in the sign
of the second crossing, but their complexities are $C_{1}=1.099$ and
$C_{2}=1.386$, respectively, as can be seen in Fig. 3(c)-(d). This indicates
that braid $\beta_{2}$ is more complex (i.e., more entangled) than braid
$\beta_{1}$, which was intuitively clear. In particular, since we only
consider systems consisting of $M=4$ chains, we will use the base-three
instead of the natural logarithm, i.e., $C=\log_{3}(\frac{\\#\beta E}{\\#E})$.
In this formulation, the braid $\beta=\sigma_{1}\sigma_{2}^{-1}$ has a
complexity equal to unity.

###  II.4 Isotopy class

Every braid can be classified according to three types of isotopy classes
given by the Thurston–Nielsen (TN) classification: reducible (RE), finite-
order (FO) or pseudo-Anosov (PA) [21, 22, 23]. The isotopy class of a braid
carries global information about the topology. For instance, if a braid is
finite-order or, alternatively, periodic if it can be isotoped to the identity
after $n$ iterations. Reducible braids can be decomposed in distinct
subbraids, where we can imagine “tubes” encompassing the subbraids, which then
themselves form a braid. Finally, braids that are not finite-order or
reducible are called pseudo-Anosov, and they represent well-entangled braids
that cannot trivially be decomposed or untied. We will use the TN isotopy
class to differentiate between twisted and braided polymer bundles.

##  III Aggregation and structural phase transitions

To study the braiding properties of a multichain system, we simulated a
coarse-grained bead-spring model system of $M=4$ chains, each with $N=30$ or
$N=10$ monomers with mass $m$, using the Kremer-Grest model [24, 25] by means
of the LAMMPS software. The total number of monomers in the system is then
$\mathcal{N}=M\cdot N$. The interaction potential between monomers is of
Lennard-Jones (LJ) type, whose van der Waals length and cohesive strength
scales are $\sigma$ and $\epsilon k_{B}T$, respectively. The hard-sphere
monomer beads are connected by strong nonlinear springs, characterised by a
finite-extensible-nonlinear spring (FENE) potential. The bending energy
between two successive bonds spanning an angle $\theta$ is equal to $U_{\rm
bend}(\theta)=\kappa k_{B}T(1-\cos(\theta))$, where $\kappa$ is the bending
parameter or elastic constant. In all simulations, we set
$k_{B}=\epsilon=\sigma=1$, such that the stiffness parameter and temperature
tune the bending energy of the chains in units of $\epsilon$. The
characteristic timescale is $\tau=\sigma\sqrt{m/k_{B}T}$ and the average bond
length is $\ell_{b}=0.965\sigma$ [26]. The time step for all simulations is
set to $\Delta t=0.01\tau$. After equilibration, the circuit topological
content is recorded. For all configurations studied, we perform the
simulations in the NVT ensemble, where we fix the temperature $T$ to be low
enough as to be in the aggregation regime. Since we do not focus on the
thermodynamics of the system in this work, we will not study thermodynamic
phase transitions and will only vary the temperature as a means to control the
aggregation of the system before varying the stiffness, ensuring that all
polymers participate in the resulting amorphous or bundled aggregate.

To determine the mutual entanglement of collections of fluctuating chains, a
primitive path analysis (PPA) is subsequently carried out to reduce the system
to a collection of straight segments that are interlocked [27]. The PPA
algorithm essentially contracts the contours of the chains while keeping the
ends fixed and without allowing strand passages, see Fig. 4. In this manner,
we combine circuit topology with entanglement and make subsequent braid
analysis easier to perform.

Figure 4: (a) The reduction of a fluctuating polymer (black line) within its
confining tube (gray) to a primitive path (orange line). Hard contacts are
indicated by colored points. (b) Topological arrangement of the individual
hard contact points (colored) within the polymer’s primitive path.

After the PPA is performed, the spatial coordinates of the monomers are
recorded and connected by means of a linear interpolation, which is justified
due to the volume exclusion we use in the simulations and the PPA algorithm;
no other monomer can occupy the space between two bonded monomers in the
polymer backbone. The resulting “trajectories” are subsequently analysed by
means of the Matlab package BraidLab [28, 23] and the writhe, braid length,
complexity index and isotopy class are computed. More details on the
computational aspects can be found in [29].

We now study the braided circuit in two aspects: the circuit topological
content and the braiding measures. For every configuration, the writhe,
complexity, braid length and isotopy class are recorded. To track the
structure of the polymer system, we will use the multichain radius of gyration
$R_{g}$, which is defined as

$R_{g}^{2}=\frac{1}{\mathcal{M}}\sum\limits_{i=1}^{\mathcal{N}}m_{i}|r_{i}-r_{CM}|^{2}\,,$
(1)

where $\mathcal{M}=\sum_{i}m_{i}$ is the total mass of all monomers, $r_{i}$
and $m_{i}$ are respectively the position and mass of the $i$th monomer, and
$r_{CM}$ is the center of mass of the full system. We subsequently make
$R_{g}$ dimensionless by rescaling by the single-chain radius of gyration,
i.e, $\widetilde{R}_{g}=R_{g}/R_{g,\theta}$, where
$R_{g,\theta}=b\,\sqrt{(N-1)/6}$ is the entropically governed radius of
gyration of a three-dimensional random walk with equilibrium bond length $b$
and stiffness $\kappa=0$.

Additionally, we define another end-to-end correlation parameter $C_{R}$ as

$C_{R}=\frac{2}{M(M-1)}\sum\limits_{i<j}\left({\bf\hat{R}_{i}}\cdot{\bf\hat{R}_{j}}\right)^{2}\,,$
(2)

where ${\bf\hat{R}_{i}}$ is the unit end-to-end vector of the $i$th polymer
[6]. This parameter is $C_{R}=1$ for completely aligned bundles, and
$C_{R}=1/3$ for uncorrelated polymer systems, i.e., in the amorphous regime.
It plays a role similar to a nematic order parameter.

Due to long equilibration times for longer polymers, we analyse the smaller
$N=10$ system in more detail, averaged over a larger ensemble of
configurations, while for the longer $N=30$ polymers we present more rough
data.

Figure 5: Density plots of the topological observables for $N=10$ as a
function of the stiffness $\kappa$. It can be easily seen that there is a
transition from an amorphous to an aligned aggregate around $\kappa\approx 7$.
The colours indicate the probability of an observable having a value given on
the $y$-axis. For high stiffness values, all observables stabilise around a
fixed distribution.

Let us first consider the case $N=10$ at a temperature of $T=1$ and density
$\rho=0.01$, where we vary the parameter $\kappa$ from $\kappa=0$ to
$\kappa=12$ in increments of $\Delta\kappa=0.2$, covering the entire stiffness
range. We average over 1000 realisations of the system per value of $\kappa$,
after equilibrating the system for a time of $20000\tau$, such that
$\widetilde{R}_{g},\,C_{R}$ have saturated and only fluctuate around their
steady-state value.

The end-to-end correlation parameter $C_{R}$ varies smoothly from
$C_{R}\approx 1/3$ to $C_{R}\approx 0.84$, not completely reaching $C_{R}=1$.
This indicates that the four chains are aligned along a common axis and that
they have transitioned from an uncorrelated, amorphous state to a more ordered
state, where the chains form an oriented bundle. These bundles can twist and
form (braided) bundles in order to minimise their energy, or as a means for
kink stabilisation as a result of defects [30]. This can be observed by
studying the braid length and complexity, which are increasing functions of
$\kappa$, and which saturate at the transition at $\kappa\approx 7$. Note also
that the average writhe fluctuates around zero, as should be the case. A more
detailed breakdown of the distributions of the aforementioned quantities as
functions of $\kappa$ is shown in Fig. 5

The isotopy class also reveals information on the topology of the system.
Generally, there is high fraction of finite-order braids, for all values of
$\kappa$. A natural explanation for this behaviour is that since for low
values of $\kappa$, the chains first collapse onto themselves and only then
aggregate. As a result, when the PPA reveals the resulting braid $\beta$, it
is unbraided, i.e., equal to the identity $\beta=e$. Since the latter is
finite-order, we see an overabundance of this isotopy class for low values of
$\kappa$. The increasing fraction of reducible and pseudo-Anosov braids for
higher values of $\kappa$ is then a natural extension of this reasoning; when
the chains become stiffer, they have a higher probability of aggregating and
forming more complex braids. For the reducible isotopy class, a subset of the
chains aggregate first and form separate braids, e.g.,
$\beta=\sigma_{1}\sigma_{3}$, and subsequently aggregate into a single
structure. The pseudo-Anosov class then results from a collective aggregation
and braiding. These statements about the isotopy class are not exhaustive; a
pseudo-Anosov braid can as well result from taking a different projection of
the system. As an example, consider the braid $\beta=\sigma_{1}\sigma_{2}$, in
which three strands cross and the fourth is left unbraided. This is a
reducible braid, with “tubes” around the first three, and around the fourth.
By rotating our projection angle slightly, we could make the third strand
undercross the fourth, making it a pseudo-Anosov braid. Since there is no
preferred direction or external force that influences the system, the
projection angle depends on the specifics of the polymer coordinates. We
simulate the system multiple times, and thus we can assume that such artefacts
can be neglected to have an influence.

Figure 6: Observables in the polymer system with $M=4$, $N=10$ at $T=1$,
$\rho=0.01$. (a) Writhe $W$, complexity $C$ and length $L$. (b) The rescaled
radius of gyration $\widetilde{R}_{g}$ and correlation parameter $C_{R}$. (c)
Circuit topology fractions of the different motifs; error bars are smaller
than symbol size. (d) Isotopy class given by the TN types (RE, FO, PA), and
representative structural conformations of the system at different $\kappa$.
All results are averaged over 1000 runs.

Let us now study the circuit topology in more detail. The ensemble averaged
circuit topological motif fractions are shown in Fig. 6(c); a more detailed
breakdown of the distribution of the topological motif fractions as a function
of the stiffness is given in the Supplemental Material [29]. From these
figures, it can be seen that already for $\kappa=0$, the fraction of single-
chain motifs, i.e., $S,\,P,\,X$ is very low. The reason is that these motifs
require the formation of two loops on a single chain, which is highly unlikely
for the short polymers we study here. As a result, the motifs that either
involve two loops on two separate chains (i.e., $I_{2}$), or one loop on a
single chain ($T_{2}$ and $I_{3}$), are somewhat more represented. Since in
our setup we study the aggregated phase, the motifs that consist of only
interchain connections, ($L_{2}$, $T_{3}$, $I_{4}$) are dominant for
$\kappa\gtrsim 2$. Within the latter motifs, the relative abundance of $T_{3}$
is a consequence of the low number of polymers in the system; a disordered
chain can more easily form two distinct bonds with two other chains, in
contrast with forming two distinct bonds with the same chain, which tends to
align both chains. Moreover, the $I_{4}$ motif requires that two bonds are
distinctly _not_ sharing a polymer. Since the aggregates are closely packed,
this situation is also more unlikely than either $T_{3}$ or $L_{2}$.
Therefore, we see that, on average, $[T_{3}]>[L_{2}]>[I_{4}]$ for all values
of $\kappa$. For increasing stiffness $\kappa$ the aggregate aligns and the
motifs that contain loops vanish. It can then be seen that the motifs
containing only interchain contacts, i.e, $L_{2}$, $T_{3}$ and $I_{4}$ tend to
dominate the system. For $\kappa\gtrsim 7$, the sum of these motifs is equal
to one, indicating that all intrachain contacts (i.e., loops) have
disappeared. Hence, the total fraction of $L_{2}$, $T_{3}$ and $I_{4}$ motifs
can be used as an order parameter to measure the degree of alignment of an
aggregate.

Let us now consider longer chains where the number of monomers is now $N=30$.
At a temperature of $T=0.1$ and density $\rho=0.01$, we equilibrate the system
for a time of $350000\tau$ until $\widetilde{R}_{g}$ has, on average, reached
its equilibrium value. Taking stiffness steps of $\Delta\kappa=1$, we again
compute the braiding quantities $W$, $C$ and $L$, along with
$\widetilde{R}_{g}$, $C_{R}$ and the CT fractions. The results are shown in
Figs. 7 and 8.

Figure 7: Observables in the polymer system with $M=4$, $N=30$ at $T=0.1$,
$\rho=0.01$. (a) Writhe $W$, complexity $C$ and length $L$. (b) The rescaled
radius of gyration $\widetilde{R}_{g}$ and correlation parameter $C_{R}$. (c)
Circuit topology fractions of the different motifs; error bars are smaller
than symbol size. All results are averaged over 1000 runs. Figure 8: Density
plots of the topological observables for $N=30$ as a function of the stiffness
$\kappa$. It can be easily seen that there is a transition from an amorphous
to an aligned aggregate around $\kappa\approx 8$. The colours indicate the
probability of an observable having a value given on the $y$-axis. For high
stiffness values, all observables stabilise around a fixed distribution.

In contrast with the $N=10$ system, we see that the topological quantities are
not monotonous anymore, but exhibit a region of increased braid length and
complexity for $\kappa\lesssim 8$, while for $\kappa\gtrsim 8$ these
quantities saturate at the values $L\approx 3$ and $C\approx 1$. From this
behaviour, it becomes clear that for very flexible chains individual collapse
occurs before aggregation, decreasing the probability of finding heavily
entangled aggregates. When the stiffness increases, however, there is a
competition between aggregation and collapse; the individual chains can now
aggregate and subsequently collectively collapse to form more intricate
entangled structures. As a consequence, the length and complexity increase
while $\widetilde{R}_{g}$ and $C_{R}$ indicate that the system is still not
aligned and are unable to capture this behaviour. For a stiffness of
$\kappa\approx 8$, the individual chains aggregate but do not collapse,
leading to lower average values for the braid length and complexity.

Note that $L$ and $C$ converge to the values $L=3$ and $C=1$ or, equivalently,
the ratio $C/L$ converges to $1/3$, which corresponds to the braid
$\beta=\sigma_{1}^{\pm 1}\sigma_{2}^{\pm 1}\sigma_{3}^{\pm 1}$ (and isotopic
equivalent braids). For four strands, this indicates that a single strand
crosses over (under) all others; none of the strands in such a braid are
interwoven. Consequently, the aligned polymer bundle is not twisted nor
braided, since this situation is a result of the choice of projection angle.
We can see that this also holds for the $N=10$ system, where $C/L\approx
0.346$.

Considering now the CT fractions, shown in Fig. 7(c) and in the Supplemental
Material [29], we see a distinct difference with respect to the $N=10$ system.
While the overall monotonous behaviour of the fractions remains the same as
before, a plateau arises for values $4\lesssim\kappa\lesssim 7$, when the
braiding quantities reach their maximal value. In this region, one can see
that there is coexistence between $T_{2}$ and $L_{2}$, and between $I_{3}$ and
$I_{4}$ motifs. This is a consequence of the energy balance between the
formation of intra- and interchain contacts, i.e., forming loops versus
forming interchain bonds. Since loops can enclose or trap other chains, they
generally lead to a higher degree of braid entanglement. Hence, in the
coexistence region, the braiding measures saturate at their maximal value as a
consequence of loops being broken.

For larger stiffness, aggregation wins out over the collapse and the CT
fractions associated with interchain interactions, i.e., $L_{2}$, $T_{3}$ and
$I_{4}$ again dominate the motif space, since the polymers are stretched along
a single direction when forming the aggregate. While $\widetilde{R}_{g}$ and
$C_{R}$ might be able to detect the structural phase transition between the
amorphous and aligned configurations, circuit topology provides a more nuanced
and rich approach to the topology of entangled systems.

Considering the $N=30$ system, we can consider the absence of the plateau
phase for the $N=10$ situation. It is indeed possible that such a plateau
phase exists for the latter situation, but since the chains are too short to
form complex braided aggregates, the stiffness range for which this occurs is
too small to see in simulations.

##  IV Discussion and conclusions

In summary, our investigation of aggregates of semiflexible polymers has
provided insights into their structural properties and phase transitions. We
observed that varying the polymer stiffness leads to a structural "phase
diagram" where different conformations such as amorphous or extended
structures are favoured, depending on the stiffness. By adjusting the latter,
we were able to study the competition between collapse and stiffening effects.
The inter- and intrachain contacts provided valuable information about the
system’s structural properties and phase transitions, which could be analyzed
within the framework of braided multichain circuit topology (CT). The
combination of circuit topology and entanglement analysis allowed for a more
comprehensive understanding of aggregation processes in systems with a small
number of chains.

The analysis of braided circuits revealed the importance of topological
measures such as writhe, complexity, braid length, and isotopy class. The
isotopy class analysis showed a higher fraction of finite-order braids for
lower stiffness values, while higher stiffness values favoured more complex
pseudo-Anosov and reducible braids. The circuit topological analysis reveals
the dominance of interchain connections, particularly motifs $T_{3}$, $L_{2}$,
and $I_{4}$, which reflects the closely packed nature of the aggregates.

When considering longer chains, we observed a region of increased braid length
and complexity, indicating a competition between aggregation and collapse. For
very flexible chains, individual collapse occurred before aggregation,
reducing the probability of finding heavily entangled aggregates. However, as
stiffness increased, the chains aggregated and collectively collapsed, forming
intricate entangled structures. The length and complexity of the braids
increased in this regime, while the radius of gyration and end-to-end
correlation parameter indicated that the system was not fully aligned. The
circuit topology fractions revealed differences compared to the smaller
system, with a plateau phase emerging for intermediate stiffness values when
the braiding quantities reached their maximum. For larger stiffness, the
dominance of interchain motifs indicated the stretching of polymers along a
single direction during aggregate formation.

Our study provided a comprehensive understanding of the structural properties,
phase transitions, and braiding characteristics of multichain systems. The
combination of circuit topology analysis and braiding measures allowed for a
deeper exploration of the system’s behavior, shedding light on the interplay
between collapse, aggregation, and entanglement. These findings contribute to
the broader understanding of semiflexible polymer systems and pave the way for
further investigations in this field. One future application we envision is
the study of the coil-globule transition for water-soluble hydrophobic polymer
chains in aqueous solutions [31]. In such systems, the explicit presence of a
polar solvent weakens the attractive interactions among monomers at
temperatures close to room temperature and the polymers collapse upon heating,
which is in contrast to polymers in inorganic solvents. Circuit topology can
then be used to study the conformational changes of such systems at the
critical temperature. Finally, the general framework of braided circuit
topology is expected to be applicable to the broader field of materials
science, including the design and analysis of nanotube, nanofiber, and
nanowire assemblies.

##  V Limitations of the study

Since we only study short identical polymers ($N=10$ or $N=30$ monomers), the
current results might not be directly applicable to longer chains or to
systems with a non-uniform length distribution; crossings can occur more often
when the chain length increases and endpoint defects can induce a higher
degree of braiding. Moreover, we here considered systems of only $M=4$ chains,
while it is known that $M$ strongly influences the structural collapse and
aggregation transitions [6]. For systems where $M\geq 4$, we expect the CT
motif fractions to heavily favour the $I_{4}$ motif when stiffness increases.
These limitations, however, can be readily addressed in future research with
the same methodology and minor software adjustments.

###### Acknowledgements.

We are grateful for inspiring discussions with I. Diamantis on the topic of
braids, and to K. Koga for suggesting the water-soluble polymer chains as an
avenue for future research.

###  V.1 Author contributions

J.B.: Conceptualization, Formal Analysis, Methodology, Software, Writing-
Original draft preparation, Visualization.
A.M.: Conceptualization, Supervision, Methodology, Writing- Original draft
preparation, Project administration, Funding acquisition, Resources.

###  V.2 Declaration of interests

The authors declare no competing interests.

## References

  * Seaton  _et al._ [2013] D. T. Seaton, S. Schnabel, D. P. Landau, and M. Bachmann, From flexible to stiff: Systematic analysis of structural phases for single semiflexible polymers, Phys. Rev. Lett. 110, 028103 (2013).
  * Zierenberg  _et al._ [2016] J. Zierenberg, M. Marenz, and W. Janke, Dilute semiflexible polymers with attraction: Collapse, folding and aggregation, Polymers 8 (2016).
  * Marenz and Janke [2016] M. Marenz and W. Janke, Knots as a topological order parameter for semiflexible polymers, Phys. Rev. Lett. 116, 128301 (2016).
  * Doniach  _et al._ [1996] S. Doniach, T. Garel, and H. Orland, Phase diagram of a semiflexible polymer chain in a $\theta$ solvent: Application to protein folding, J. Chem. Phys. 105, 1601 (1996).
  * Bastolla and Grassberger [1997] U. Bastolla and P. Grassberger, Phase transitions of single semistiff polymer chains, J. Stat. Phys. 89, 1061 (1997).
  * Zierenberg and Janke [2015] J. Zierenberg and W. Janke, From amorphous aggregates to polymer bundles: The role of stiffness on structural phases in polymer aggregation, EPL 109, 28002 (2015).
  * Luo and Sommer [2012] C. Luo and J.-U. Sommer, Disentanglement of linear polymer chains toward unentangled crystals, ACS Macro Lett. 2, 31–34 (2012).
  * Nguyen  _et al._ [2015] H. T. Nguyen, T. B. Smith, R. S. Hoy, and N. C. Karayiannis, Effect of chain stiffness on the competition between crystallization and glass-formation in model unentangled polymers, J. Chem. Phys. 143, 10.1063/1.4932193 (2015).
  * Kawak  _et al._ [2021] P. Kawak, D. S. Banks, and D. R. Tree, Semiflexible oligomers crystallize via a cooperative phase transition, J. Chem. Phys. 155, 10.1063/5.0067788 (2021).
  * Martínez-Fernández  _et al._ [2023] D. Martínez-Fernández, M. Herranz, K. Foteinopoulou, N. C. Karayiannis, and M. Laso, Local and global order in dense packings of semi-flexible polymers of hard spheres, Polymers 15, 551 (2023).
  * Mashaghi  _et al._ [2014] A. Mashaghi, R. van Wijk, and S. Tans, Circuit topology of proteins and nucleic acids, Structure 22, 1227 (2014).
  * Heidari  _et al._ [2020] M. Heidari, H. Schiessel, and A. Mashaghi, Circuit topology analysis of polymer folding reactions, ACS Cent. Sci. 6, 839 (2020).
  * Golovnev and Mashaghi [2020] A. Golovnev and A. Mashaghi, Generalized circuit topology of folded linear chains, iScience 23, 10.1016/j.isci.2020.101492 (2020).
  * Golovnev and Mashaghi [2022] A. Golovnev and A. Mashaghi, Topological analysis of folded linear molecular chains, in  _Topological Polymer Chemistry: Concepts and Practices_, edited by Y. Tezuka and T. Deguchi (Springer Nature Singapore, Singapore, 2022) pp. 105–114.
  * Berx and Mashaghi [2023] J. Berx and A. Mashaghi, Decoding chirality in circuit topology of a self entangled chain through braiding, Soft Matter 19, 5888 (2023).
  * Heidari  _et al._ [2022] M. Heidari, D. Moes, O. Schullian, B. Scalvini, and A. Mashaghi, A topology framework for macromolecular complexes and condensates, Nano Res. 15, 9809 (2022).
  * Birman and Brendle [2005] J. S. Birman and T. E. Brendle, Chapter 2 - braids: A survey, in  _Handbook of Knot Theory_, edited by W. Menasco and M. Thistlethwaite (Elsevier Science, Amsterdam, 2005) pp. 19–103.
  * Gügümcü  _et al._ [2022] N. Gügümcü, B. Gabrovsek, and L. H. Kauffman, Invariants of bonded knotoids and applications to protein folding, Symmetry 14, 10.3390/sym14081724 (2022).
  * Lee  _et al._ [2013] D. J. Lee, R. Cortini, A. P. Korte, E. L. Starostin, G. H. M. van der Heijden, and A. A. Kornyshev, Chiral effects in dual-DNA braiding, Soft Matter 9, 9833 (2013).
  * Dynnikov and Wiest [2007] I. A. Dynnikov and B. Wiest, On the complexity of braids, J. Eur. Math. Soc. 9, 801–840 (2007).
  * Nielsen [1944] J. Nielsen,  _Surface Transformation Classes of Algebraically Finite Type_ , Danske Vid. Selsk. Math.-Phys. Medd., (Munksgaard, 1944).
  * Thurston [1988] W. P. Thurston, On the geometry and dynamics of diffeomorphisms of surfaces, Bull. Amer. Math. Soc. 19, 417 (1988).
  * Thiffeault [2022] J.-L. Thiffeault,  _Braids and Dynamics_ (Springer International Publishing, 2022).
  * Grest and Kremer [1986] G. S. Grest and K. Kremer, Molecular dynamics simulation for polymers in the presence of a heat bath, Phys. Rev. A 33, 3628 (1986).
  * Kremer and Grest [1990] K. Kremer and G. S. Grest, Dynamics of entangled linear polymer melts: A molecular-dynamics simulation, J. Chem. Phys. 92, 5057 (1990).
  * Hunt and Todd [2009] T. A. Hunt and B. Todd, A comparison of model linear chain molecules with constrained and flexible bond lengths under planar couette and extensional flows, Mol. Simul. 35, 1153–1167 (2009).
  * Caraglio  _et al._ [2017] M. Caraglio, C. Micheletti, and E. Orlandini, Physical links: defining and detecting inter-chain entanglement, Sci. Rep. 7, 1156 (2017).
  * Thiffeault and Budisic [2014] J.-L. Thiffeault and M. Budisic, Braidlab: A software package for braids and loops (2014).
  * [29] Link_to_Supplemental_Material.
  * Slepukhin  _et al._ [2021] V. M. Slepukhin, M. J. Grill, Q. Hu, E. L. Botvinick, W. A. Wall, and A. J. Levine, Topological defects produce kinks in biopolymer filament bundles, PNAS 118, e2024362118 (2021).
  * Hatano  _et al._ [2016] I. Hatano, K. Mochizuki, T. Sumi, and K. Koga, Hydrophobic polymer chain in water that undergoes a coil-to-globule transition near room temperature, J. Phys. Chem. B 120, 12127 (2016).