# Phase Behaviour and Dynamics of Three-Dimensional Active Dumbbell Systems

C. B. Caporusso1∗, G. Negro1111These authors contributed equally to this
work., A. Suma1, P. Digregorio2,3, L. N. Carenza4, G. Gonnella1, L. F.
Cugliandolo5,6 1Dipartimento Interateneo di Fisica, Universitá degli Studi di
Bari and INFN, Sezione di Bari, via Amendola 173, Bari, I-70126, Italy,
2Departement de Fisica de la Materia Condensada, Facultat de Fisica,
Universitat de Barcelona, Martí i Franquès 1, E08028 Barcelona, Spain,
3UBICS University of Barcelona Institute of Complex Systems, Martí i Franquès
1, E08028 Barcelona, Spain,
4Instituut-Lorentz, Universiteit Leiden, P.O. Box 9506, 2300 RA Leiden,
Netherlands,
5Sorbonne Université, Laboratoire de Physique Théorique et Hautes Energies,
CNRS UMR 7589, 4 Place Jussieu, 75252 Paris Cedex 05, France
6Institut Universitaire de France, 1 rue Descartes, 75231 Paris Cedex 05,
France

###### Abstract

We present a comprehensive numerical study of the phase behavior and dynamics
of a three-dimensional active dumbbell system with attractive interactions. We
demonstrate that attraction is essential for the system to exhibit nontrivial
phases. We construct a detailed phase diagram by exploring the effects of the
system’s activity, density, and attraction strength. We identify several
distinct phases, including a disordered, a gel, and a completely phase-
separated phase. Additionally, we discover a novel dynamical phase, that we
name percolating network, which is characterized by the presence of a spanning
network of connected dumbbells. In the phase-separated phase we characterize
numerically and describe analytically the helical motion of the dense cluster.

##  I Introduction

Active systems are a striking class of soft matter which employs internal
energy stored in some kind of reservoir and converts it into motion
marchetti2013 ; Bechinger2016 ; Shankar2022 ; carenza2019 . Active systems
greatly vary in their extension, ranging from macroscopic to microscopic
assemblies of active constituents, and a multitude of different realizations.
Focusing on the microscopic world, active systems can either have synthetic
origin (_e.g._ Janus particles Walther2013 ; Das2015 ; ebbens2016 ) or
biological origin. In this latter case, cytoskeletal suspensions sanchez2012
and cellular cultures bookhoward ; saw2017 ; Guillamat2022 ; Armengol2023  are
examples of _in vitro_ systems which are often employed to study active
behaviors in a controlled environment, but more complex living organisms, such
as developing embryos or living tissues Haeger2019 ; MaroudasSacks2021 ;
Brauns2023  have recently been considered for their promising implications in
understanding a variety of different biological phenomena from morphogenesis
to cancer progression and spreading of infections in living organisms.
Regardless of their particular realization, active systems can exploit energy
to interact with the surrounding environment and perform autonomous motion,
leading to a plethora of collective behaviors.

In the past two decades much effort has been spent to identify the physical
rules underlying the behavior of living and active systems marchetti2013 ;
Hatwalne2004 ; Kruse2004 ; cates2015 , with important repercussions on our
understanding of their non-equilibrium dynamics. Among others, bacteria have
gained the attention of the physics community as they represent a simple
natural realization of self-propelled particles, by virtue of their biological
simplicity and limited ability in interacting with both the environment and
other units. Therefore, bacterial systems provide an elementary yet relevant
example of how the dynamics of active constituents may lead to self-assembly
Bratanov2015 ; Bi2018 ; MassanaCid2022 . Indeed bacteria exhibit enhanced
chemotactic properties which allow them to respond to external stimuli such as
variations of temperature, nutrient availability, oxygen concentration, _etc._
to develop different kinds of colonies, ranging from biofilms to fluidic
suspensions, often characterized by the emergence of elaborate patterns, even
in absence of guidance and only because of the uncoordinated evolution of the
separate units Yaman2019 .

The chemotactic properties of bacteria are often exploited in experiments to
trigger a particular response. For instance, most bacteria preferentially
reproduce at room temperature –process during which they mostly remain
motionless Basaran2022a ; Basaran2022b – while at higher temperatures they can
either develop biofilms if the growing substrate is dry or become motile by
growing flagella –a proteic protrusion which is used by some bacterial species
to swarm and swim in humid environments. In this latter case, oxygen
availability or light intensity determines the dynamical response of the
system and can be used in experiments as a control parameter to tune the
typical speed of migration. When this is large enough, dense suspensions of
flagellate bacteria in a fluidic environment arrange in a liquid crystalline
fashion Yaman2019  and develop a chaotic flowing state characterized by
whirling patterns which resemble those observed in an isotropic fluid flowing
at large Reynolds numbers –a feature which earned this dynamical state the
title of bacterial (or active) turbulence Dunkel2013 ; giomi2015 ; carenza2020
; carenza2020_bif ; Aranson2022 ; Alert2022 .

Inspired by this huge variety of possible states, many models have been
advanced to capture the dynamics of bacterial suspensions. They range from
particle (or molecular) models marchetti2013 ; Locatelli2015 ; Solon2015
–where each unit in the system evolves according to a given equation of motion
and interactions– to continuous models cates2022 ; wittkowski2014 ;
Arciprete2018 ; negro2019b ; Giordano2021 ; Carenza2020b ; Favuzzi2021 ;
Carenza2019d  –where the status of the system is described in terms of a few
continuous fields capturing only the slowly-varying hydrodynamic features of
the system (local group velocity, polarization, density, _etc._).

In the following we will take the former _molecular_ approach and describe
every constituent as a dumbbell –namely two joint spheres linked with a rigid
spring – autonomously moving in a viscous medium and subject to thermal
fluctuations. This particular realization falls in the wide class of the so-
called _self-propelled active Brownian_ constituents. Many recent studies have
focused on different elements ranging from disks digregorio2019 ;
caporusso2020 ; Caporusso-2023 ; cagnetta2017 ; cates2015  to rods Bar2020 ,
and contributed to explicate the effects of non-equilibrium in phenomena like
motility induced phase separation (MIPS) cates2015  and arrested phase
separation Sciortino2008b , providing significant insights on the physics of
collective effects of bacterial suspensions and other similar active systems.
Interestingly enough, it was repeatedly suggested that active systems could
hinge on a minimal toolkit of physical mechanisms, independently of the
particular realization, to achieve a specific dynamical response. However,
most research has focused on two-dimensional systems, while much less is known
on three-dimensional realizations, relevant for a full characterization of
real biological systems, as these are most often subject to a full $3D$
dynamics.

To better outline the scope of our study, we find it useful to summarize some
important results obtained in systems of active Brownian particles (ABP) both
in $2D$ and $3D$, before proceeding with the presentation of our findings. In
Sec. I.1, we provide a brief summary of previous studies, regarding the
properties of MIPS at varying the system’s dimensionality and the particles’
interaction. Far from being a complete review on ABP, we suggest the
interested reader to refer to Ref. cates2015 . Next, in Sec. I.2 we motivate
our study and summarize our results.

###  I.1 Brief summary of previous studies

Most literature on active Brownian systems focused on two-dimensional
realizations, both for their biological relevance as for the important
implications on the melting mechanism in $2D$ Halperin78 . In recent years,
there has been significant progress in understanding the non-equilibrium
behavior of active Brownian particles, with much attention on disks and
dumbbells. Both types of systems exhibit a common characteristic, namely the
segregation into dense and dilute phases caused solely by repulsive two-body
forces. This process, known in the literature as motility-induced phase
separation (or shortly MIPS), consists of two steps. First, scattering of
active particles results into the slowdown of the particles themselves;
second, slower particles tend to aggregate into dense clusters which
themselves displace and also aggregate caporusso2022dynamics ; Caporusso-2023
. This gives rise to a positive feedback which ultimately produces a
macroscopic change of state, from a gas-like to fluid-like phase, mediated by
the non-equilibrium injection of energy at the level of the individual
constituents. More specifically, for dumbbells a coexistence region between
dense (hexatic) and dilute (liquid) phases, found in the absence of activity,
extends continuously to high values of activity cugliandolo2017 . For disks,
instead, between the hexatic-liquid coexistence, and the MIPS, occurring at
low and high activity, respectively, there is a region at intermediate
activity where no coexistence between different phases is found digregorio2018
. Other special features of the active disk systems include micro-separation
of hexatic domains or the formation of cavitation bubbles in the interior of
the dense aggregates digregorio2018 ; cugliandolo2017 ; caporusso2020  which
take different form in active diatomic systems sum_gonn_mar ; Suma2014 ;
Suma2014pre ; petrelli2018 ; Caporusso-2023 .

While the panorama is pretty much understood for two-dimensional realizations
of active systems, much less is known about the phase behaviour of these same
systems in three-dimensional space Stenhammar2014 ; wysocki2014 ; vandamme2019
; nie2020 ; omar2021 ; Turci21 ; Venkatareddy23 . Stenhammar _et al._
Stenhammar2014  first observed that the MIPS transition is significantly
suppressed in 3D for active Brownian disks in absence of any attractive
interactions. More recently, van Damme _et al._ vandamme2019  pointed out that
MIPS is also suppressed when the active constituents are spherocylinder-
shaped. These authors argued that MIPS suppression is ultimately due to the
presence of passive torques, causing a reduction in the particle-particle
collision time, and showed that MIPS eventually disappears altogether for
sufficiently elongated particles. These studies suggest that in 3D activity
alone is not enough to drive the system towards a phase-separated state. These
observations were further investigated by Omar _et al._ omar2021  who have
recently conducted a comprehensive study of the phase diagram of active
Brownian spheres in 3D focusing also on the process of crystallization, _i.e._
the coexistence between liquid and solid phases. They found that phase
coexistence not only occurs out of equilibrium but it is also the most stable
state in 3D and dominates over the coexistence between liquid and gas typical
of MIPS.

Three-dimensional systems of active spherical particles Mognetti2013  and
dumbbells schwarz2012  –both in presence of attractive interactions– have been
studied numerically with the purpose of addressing the experimental
realization of bacterial systems in the presence of non-absorbing polymers
Melaugh2019 ; schwarz2012 . In this case, demixing into coexisting phases is
tamed by the aggregation strength, ultimately leading to the formation of
highly dynamical structures as rotating clusters schwarz2012 . However,
activity tends to suppress aggregation –consistently with experimental
observations of suspensions of motile bacteria. The non-equilibrium nature of
the phase coexistence is also illustrated by the formation of new dynamical
states, such as percolating networks caused by the interplay between
attraction and motility prymidis2015 .

Furthermore, the phase behavior of 3D passive attractive atomic and molecular
systems is very rich in itself and is not fully understood yet Zaccarelli2007
. This is because the way in which many-body systems interact is highly
influenced by the geometry of the elementary constituents and different
macroscopic behaviors are observed in systems with different shapes. The
interest in these systems has exploded in the last 15 years or so, thanks to
the possibility of engineering colloids with different shapes –sphere’s
diameter and separation length– and tuning particle-particle interactions
Kraft12 , paving the way towards the design of colloidal crystals with useful
optical properties. Indeed, dumbbells model diatomic molecules and show
interesting crystal structures at equilibrium Vega92 ; Marechal08 ; Kowalik08
; Dennison12 . In the context of self-assembly studies Jack15 , (patchy)
dumbbells have been the focus of much attention Whitelam10 ; Munao13 ; Munao14
; Munao15 ; Avvisati15 ; OToole17 . In the limit of weak activity, the 3D
model that we will consider in this article, would be relevant to the
description and characterization of these very interesting passive systems.

###  I.2 Motivation and summary of results

In this introduction, we discussed how clustering of bacterial suspensions can
be explained in terms of simple physical mechanisms when modelled as systems
of active Brownian constituents. Importantly, most of these results and
observations only hold in two-dimensions and cannot be invoked to explain the
coarsening of bacterial aggregates in three-dimensions. However, collective
effects in bacterial suspensions are also observed in 3D.

Our goal in this paper is to clarify the relevance of motility in the dynamics
of three-dimensional suspensions by performing a systematic analysis of the
phase behavior of a system of propelled attractive dumbbells.

The article is structured as follows. In Section II we will present the model
for active Brownian dumbbells. The discrete mesoscopic approach here
implemented allows us to retain only the main features of bacterial suspension
(local orientation, attractive interactions) and to wipe out more system-
dependent features. Section III is devoted to the analysis of the phase
behavior of the system in a three-dimensional unconfined geometry. We will
identify four phases (gel, disordered, phase-separated and percolating
network) while varying two control parameters, namely the Péclet number ${\rm
Pe}$ measuring the strength the activity with respect to the thermal
fluctuations and the total packing fraction $\phi$. The effect of the
attraction strength will also be considered. Next, in Sec. IV we will focus on
the description of the motion of a single isolated cluster in the phase-
separated phase. We close the paper with some conclusions in Sec. V.

##  II Model and numerical methods

We simulate a system of $N$ dumbbells immersed in a three dimensional space.
Each dumbbell is a dimer composed by two spherical beads of mass $m$ and
diameter $\sigma$. The beads are rigidly connected and their center-to-center
distance is kept fixed at distance $\sigma$. Each bead $l$ evolves in time
according to the Langevin equation

$m\ddot{\bm{r}}_{l}=-\gamma\dot{\bm{r}_{l}}+{\bm{f}_{\text{act}}}_{l}-\bm{\nabla}_{l}U+\bm{\xi}_{l}\;,$
(1)

where $l=1,..,2N$. $\gamma$ is the damping coefficient and $U$ is the total
potential energy. The term $\bm{\xi}_{l}$ is a Gaussian white noise, with zero
average and variance fixed by the fluctuation-dissipation theorem as

$\displaystyle\langle\xi_{l\alpha}(t)\rangle$ $\displaystyle=$ $\displaystyle
0\;,$ (2) $\displaystyle\langle\xi_{l\alpha}(t_{1})\xi_{p\beta}(t_{2})\rangle$
$\displaystyle=$ $\displaystyle
2k_{B}T\gamma\delta_{lp}\delta_{\alpha\beta}\delta(t_{1}-t_{2})\;,$ (3)

where $\alpha,\beta=1,2,3$ label the spatial coordinates, $T$ is the
temperature of the equilibrium bath and $k_{B}$ the Boltzmann constant.

Figure 1: Simulation set up. (a) Mie potential with truncated attractive tail
at distance $r_{c}=1$ (blue curve) and three other attractive energy,
$\epsilon$, and cut-off distance, $r_{c}$, choices. In all cases $n=32$. (b) A
snapshot of an active dumbbell system with pure repulsive interbead
interactions in a 3D periodic box. $\phi=0.01$ and Pe $=150$. Without
attraction between the dumbbells MIPS is absent even at higher densities and
activity.

The self-propulsion is represented by the active force
${\bm{f}_{\text{act}}}_{l}$. This force has the same fixed modulus $f_{\rm
act}$ for all disks. It points from the tail to the head monomer in the
dumbbell and acts along the line that joins the centers of the two beads. The
tail and head identity of the beads in a molecule are attributed randomly at
the beginning of the simulation and they are kept fixed along the numerical
experiment.

Beads belonging to different dumbbells interact through the Mie short-range
potential:

$U_{\text{Mie}}(r)=4\epsilon\left[\left(\frac{\sigma_{\text{Mie}}}{r}\right)^{2n}-\left(\frac{\sigma_{\text{Mie}}}{r}\right)^{n}\right]\theta(r_{c}-r)\;,$
(4)

with $n=32$, $\sigma_{\text{Mie}}=2^{-1/n}\sigma$, $r$ the distance between
the concerned beads, $r_{c}$ the potential cutoff and $\theta(r)$ the
Heaviside step function. The variation of the potential $U_{\rm Mie}$ with the
parameter $\epsilon$ controlling the potential energy scale is shown in Fig.
1(a). This potential has a pronounced minimum at $r_{\rm
min}=2^{1/n}\sigma_{\rm Mie}=\sigma$ and it allows colloids to behave
similarly to hard disks at distances shorter than $\sigma$ digregorio2018 ;
cugliandolo2017 , while it acts attractively at longer distances. Unless
otherwise stated, we use $r_{c}=1.5\sigma$. The active dumbbells move in a
cubic box with fully periodic boundary conditions, as shown in Fig. 1(b).

Before entering into the heart of our study, we confirmed that there is no
stable MIPS Turci21  for purely repulsive interactions between the active
dumbbells. The condition for purely repulsive dumbbells is achieved by setting
the cut-off lengthscale $r_{c}=\sigma$, so that only the repulsive core in the
Mie potential is preserved (blue curve in Fig. 1(a)). A typical configuration
of the steady state of repulsive active dumbbells interacting via a truncated
potential with $n=32$ in the cubic box is shown in Fig. 1(b), clearly showing
the absence of particle aggregation.

The system can be characterized through the following dimensionless numbers.
The relevance of active injection with respect to thermal fluctuations is
captured by the Péclet number

$\text{Pe}=\frac{2f_{\text{act}}\sigma}{k_{B}T}\;,$ (5)

which is the dimensionless ratio between the advective transport rate and the
diffusive transport rate. Moreover, the aggregation capability of the system
can be captured by measuring the ratio between the strength of the attractive
force and the active force Mognetti2013

$P_{\rm agg}=\frac{\epsilon}{f_{\rm act}\sigma}\;.$ (6)

Finally, the global density or packing fraction is quantified by the volume
fraction

$\phi=N\,\frac{\pi\sigma^{3}}{3V}\;,$ (7)

with $V$ the volume of the box where the dumbbells are placed and $N$ the
total number of dumbbells, that is, twice the number of spherical beads
222Although in the following we will mostly consider relatively low packing
fractions, $\phi\leq 0.3$, it is worth mentioning that the close-packing value
$\phi_{\rm cp}\simeq 0.74$ and the random close packing value $\phi_{\rm
rcp}\simeq 0.64$ for hard spheres in 3D. Close packing is achieved by the face
center cubic and the hexagonal close packing arrangements. In both cases each
sphere touches 12 neighboring spheres Paddy23 ..

The evolution equations were integrated using the velocity-Verlet algorithm
with the open-source code LAMMPS LAMMPS . Both beads of each dumbbell are
constrained at fixed distance via the SHAKE algorithm, which adjusts the
positions of the beads by applying an additional force so that to fix the bond
length at the next iteration step. The required force vector at each iteration
is found by solving a system of linear equations refSHAKE . We apply the
DBSCAN algorithm for cluster identification dbscan ; digregorio2019 ;
caporusso2020 . DBSCAN is a density-based clustering algorithm that can
identify clusters of arbitrary shape and size, grouping together points that
are closely packed in some space, and marking as outliers points that lie
alone in low-density regions. DBSCAN works with two parameters: $R$ and
$n_{\text{min}}$. $R$ specifies the radius of a neighborhood around a point,
and $n_{\text{min}}$ specifies the minimum number of points required to form a
dense region. A point is classified as a core point if there are at least
$n_{\text{min}}$ other points within $R$ distance from it (including itself).
A point is classified as a border point if there are less than
$n_{\text{min}}$ points within $R$ distance from it, but is within $R$
distance of a core point. A point is classified as a noise point if it is
neither a core nor a border point. In our case, we choose $R=1.5\sigma$, where
$\sigma$ is the particle diameter, and $n_{\text{min}}=12$, based on the
preferred structures at close packing. Then, DBSCAN assigns each core point to
a cluster, along with all the points that are density-reachable from it (i.e.,
there is a path of core points connecting them). Border points are assigned to
the cluster of their nearest core point, and noise points are left unassigned.

All the quantities are made dimensionless in terms of the bead mass $m$, its
diameter $\sigma$ and the unit of energy $\epsilon$. We chose $k_{B}T=0.05$
and $\gamma=10.0$, which ensure that the dynamics is over-damped. We use
$N=4096$ dumbbells and vary the box size to select the volume fraction. In the
following we present the phase diagram resulting by varying $\phi$ from
$0.001$ to $0.3$ and the Péclet number Pe between $0$ and $200$. We use a Mie
potential with $n=32$ and the default value $\epsilon=0.3$ for the energy
scale, though we have also explored cases with $\epsilon=0.2$ and
$\epsilon=0.4$.

##  III Phase Behaviour

We start by presenting the phase diagram obtained by scanning various values
of the packing fraction $\phi$ and the Péclet number $\rm{Pe}$. In Fig. 2,
snapshots are presented for different combinations of these two parameters. We
find that the system exhibits a variety of diverse phases. We start by
discussing them at a pictorial level and we later give a quantitative
characterization by studying the distribution of local densities, the
distribution of cluster sizes, and the dynamics.

###  III.1 The structures

Figure 2: Phase behavior of attractive active dumbbells. Snapshots of the
dumbbells’ system for different global densities $\phi$ and Péclet numbers Pe.
The attraction $\epsilon=0.3$ for all cases presented here. Particles are
colored according to the local density. From top to bottom, the rows
correspond to Pe $=0,\,50,\,100,\,150$. From left to right, columns correspond
to $\phi=0.05,\,0.1,\,0.2,\,0.3$. At Pe = $0$, panels (a)-(d), the system is
in the gel phase. Increasing activity for the same densities, panels (e)-(h)
for ${\rm Pe}=50$, or for even higher ${\rm Pe}=100$ and small densities
(i)-(j), aggregation is favoured and there is complete phase separation, in
the sense described in the main text. At Pe = 150 and low density, panels
(m)-(n), the structure breaks down into small pieces. At higher densities,
panels (o)-(p), the percolating network subsists. The bar on top of the plot
shows the scale for the local density with which the beads are colored, the
scale we use to distinguish the phases, see Fig. 3 and its discussion,
together with the dynamic properties to be presented in Sec. III.5.

We commence our analysis by setting Pe = 0 and examining the structures formed
at long times for the global densities $\phi=0.05,\,0.1,\,0.2,\,0.3$, as shown
in the first row of Fig. 2. In the passive limit the dumbbells aggregate in a
gel phase, in which particles form a rather static percolating network. This
behavior, called gelation Poon02 ; Coniglio04 ; Zaccarelli2007 ;
Sciortino2008a ; Sciortino2008b , is due to a rapid cooling of an initially
disordered conformation which in conjunction with the dumbbells’ anisotropy
and the short-range attraction induce the formation of such a metastable
state. The stable state would be instead phase-separated – a configuration
which, nonetheless, cannot be reached since the thermal energy is insufficient
for the dumbbells to rearrange into a single compact cluster.

In the second row of Fig. 2 we consider a higher activity Pe $=50$, while
maintaining the same densities as in the first row. The dumbbells now
aggregate to form a single cluster at any density. We will refer to this state
as phase-separated. The resulting cluster exhibits both translational and
rotational persistent motion, which we will characterize in detail in Sec. IV.
Remarkably, a relatively small increase in ${\rm Pe}$ with respect to the
passive case, induced the dumbbells to reorganize destroying the gel network
and undergoing phase separation, thanks to the activity that acts as an
additional source of noise.

When the activity is further increased to ${\rm Pe}=100$ (third row), an
interesting effect is observed upon varying the density. Specifically, at
densities $\phi=0.05$ (i) and $\phi=0.1$ (j), we still find phase separation.
However, at higher densities $\phi=0.2$ (k) and $\phi=0.3$ (l), the system
becomes more disordered, with small clusters that do not grow in size over
time. These clusters connect to each other, forming a phase that is similar to
the gel observed at vanishing activity but is more dynamic in nature, with the
network connections constantly forming and breaking throughout the simulation
(see SM Movie 1) movies . This phase will be referred to as percolating phase,
following the nomenclature of Ref. prymidis2015  where a similar behavior was
found for a system of self-propelled Lennard-Jones disks.

Finally, at ${\rm Pe}=150$ (fourth row), we notice that at low densities the
activity causes the system to become fully disordered (Fig. 2(m,n). In this
case, the activity is strong enough to overcome the aggregation force, leading
to the cluster’s breakup. Conversely, at higher densities, (o) and (p), we
observe again the formation of a percolating network.

Despite the presence of attraction, we have observed four phases occurring in
distinct ways. Such findings suggest the existence of a complex phase diagram,
which we will now proceed to characterize.

###  III.2 The local densities

First of all, we reckon that there are no isolated dumbbells once the steady
state is reached. In order to identify the different phases, we developed a
method involving the construction of a Voronoi tessellation Okabe2000 . Such a
tessellation partitions the simulation box into a set of regions –the Voronoi
cells– surrounding each bead such that all points within each cell are closer
to the pertaining bead than any other beads. By computing the local packing
fraction within each region, $\phi_{i}$, defined as the ratio between a single
bead’s volume and its Voronoi cell volume, we obtain valuable information
about the local packing and structure of the system. In fact, this quantity
expresses how much free space the bead has and thus whether it is free, or
located in a cluster.

Figure 3: Local Voronoi density. Distribution of the local packing fraction of
Voronoi cells for different global density and Péclet numbers in the four
phases: $\phi=0.10$ and $\text{Pe}=0$ in the gel, $\phi=0.10$ and
$\text{Pe}=100$ in the phase-separated state, $\phi=0.10$ and $\text{Pe}=200$
in the disordered phase, and $\phi=0.30$ and $\text{Pe}=150$ in the
percolating network. The high density peak in the gel and phase-separated
phase are at the close packing value, shown with a dotted vertical line, while
the next one close to it is in between $\phi_{\rm rcp}\sim 0.64$ and
$\phi_{\rm cp}\sim 0.74$, shown with a dashed vertical line, respectively.

In Fig. 3, we present the probability distributions of the local density
$\phi_{i}$ in the four distinct phases identified in the study (the global
density and ${\rm Pe}$ values are specified in the caption). In each phase the
distribution has unique characteristics which can be used as a fingerprint of
each phase to classify different configurations.

In the gel phase, the distribution has three well-defined peaks, one at low
and two at high density. The low density peak is due to the beads placed on
the border of the gel; these beads are located inside very large Voronoi cells
and hence return a very low $\phi_{i}$. Instead, the high density peaks are
due to beads inside the dense regions of the gel. One is centered around the
close packing value $\phi_{\rm cp}\sim 0.74$, while the other is at a value in
between the random close packing $\phi_{\rm rcp}\sim 0.64$ and $\phi_{\rm
cp}\sim 0.74$. We note that for these parameters the low density peak carries
more weight than each of the two high density ones. This feature also holds
for other values of the parameters corresponding to a gel network.

Figure 4:  Structural analysis of active clusters. (a) Snapshot of the system
for Pe $=50$ and $\phi=0.10$ (phase-separated phase), with particles colored
accordingly to the number of neighbouring particles inside a cutoff distance
of $1.5\sigma$. (b) Snapshot of the same configurations with particles colored
accordingly to the common neighbor analysis FAKEN1994279 . Particles with FCC
order are in blue, those with HCP order in red, while particles with no
crystalline order are in grey. (c) Distribution of the local packing fraction
of Voronoi cells for a configuration in the phase-separated state. The dotted
colored lines represent the contribution to the distribution of the beads
having a selected Coordination Number (CN), reported in the legend.

In the phase-separated systems, the low density peak also corresponds to disks
which are on the single cluster surface. At high local densities, there is a
double peak structure with one of the peaks sitting on the close packing
density and the other one at a slightly lower value, very similar to what we
found in the gel. There is, though, an important difference with the gel data,
which resides on the height of these peaks: in the phase-separated phase a
much larger portion of particles have the large densities of these two peaks.

To elucidate the origin of the two peaks in the local density distributions of
the gel and phase-separated phase, we analyze the relative frequency of beads
with different coordination numbers. The coordination number of a bead is
defined as the number of its nearest neighbours in the Voronoi sense, that is,
two beads are considered to be nearest neighbours if they share a common face
of the tesselation, and their centers are located within a distance of
$1.5\,\sigma$, Fig. 4(a). The data in Fig. 4 demonstrate that beads with
coordination number 12 contribute to the peak at the close packing value. More
precisely, these beads belong to either a face-centered cubic (fcc) or a
hexagonal close-packed (hcp) local crystalline structure in the interior of
the cluster. As shown in Fig. 4 (b), the two structures coexist and both have
local density $\phi\simeq\phi_{\rm cp}$, consistently with the location of the
purple peak in Fig. 4 (c). The secondary peak is associated with beads that
have a coordination number equal to 11, and are topological defects that
locally disrupt the regular lattice and lower the local packing fraction. The
crystalline structure in the bulk of the cluster resembles very closely the
high-density ground state of a system of passive dumbbells, which is a so-
called aperiodic crystal, with all beads arranged on a close-packed ordered
structure and disordered bonds between beads Dijkstra_3d-dumbbells . The only
difference here is that we see a mix of fcc and hcp structures. This
arrangement is very common in close-packed spheres, and the interfacial planes
are called stacking faults. We note that our configurations are not ground
state and the system is active, so it is reasonable to expect other close-
packed structures to be stable. Particles with a lower coordination number
contribute to the rest of the distribution and, in particular, to the lower
density peak which, as mentioned above, is due to the beads on the surface of
the cluster, thus with a smaller number of neighbors compared to the core
ones, see Fig. 4(a).

In order to distinguish between the gel and phase-separated phases, it is
necessary to check the system’s motility: in the gel, dumbbells are almost
frozen and cannot diffuse, conversely in the phase-separated case diffusion is
observed, see Fig. 8. More detailed characterizations of the motility of the
full system, and of the single clusters of the phase-separated phase, are
given in Sec. III.5 and Sec. IV, respectively.

In the disordered phase, the beads are more homogeneously distributed,
resulting in an equally partitioned space and the development of only a single
peak at small density, around $\phi_{i}\sim 0.1$. The percolating network has
a significantly different $\phi_{i}$ distribution. It is broad, with the
appearance of all intermediate densities, with almost equal probability. As
noted in the fourth row of Fig. 2, the system smoothly transitions between the
disordered and the percolating phases.

###  III.3 The cluster size distribution

In order to identify the global density $\phi$ at which the percolating
network appears, one can proceed as in the usual analysis of percolation
Stauffer  and gel Sciortino2008a ; Zaccarelli2007 ; Sciortino2008b
transitions. A common observable to detect a percolating structure is the
cluster size distribution, displayed in Fig. 5 for ${\rm Pe}=150$ and two
packing fractions, $\phi=0.1,\,0.2$. After clustering each conformation using
DBSCAN (see the method description in Sec. II), we compute the distribution of
cluster sizes. Near the density where clusters start to connect to form a
percolating network, the cluster size distribution takes a power-law form with
an exponential cutoff Stauffer  ($\phi=0.1$ in Fig. 5). Notably, the Fisher
exponent $\tau$ in the algebraic law $N_{c}^{-\tau}$ matches the one of
standard random percolation for a three-dimensional system, $\tau\sim 2.18$,
associated to a fractal dimension $d_{f}\sim 2.53$ by a hyperscaling relation.
Increasing density, percolation starts to occur and the system becomes
interconnected, with a large cluster and a few small ones. In this case, the
cluster size distribution appears to have a peak at a cluster size of the
order of the number of beads in the system ($\phi=0.2$ in Fig. 5).

Figure 5:  Cluster size distribution in the gel, disordered and percolating
network cases. (a) Probability distribution of cluster sizes $N_{c}$ at Pe
$=5$, for two values of the global packing fraction $\phi=0.01$ (grey curve
corresponding to the disordered phase) and $\phi=0.05$ (green curve
correspoding to the gel phase). Panel (b) shows the cluster size distribution
for $Pe=150$, and $\phi=0.1$ (black curve corresponding to the disordered
phase) and $\phi=0.2$ (yellow curve for the percolating network phase). The
dotted line corresponds to an algebraic decay $N_{c}^{-\tau}$ with Fisher
exponent $\tau=2.18$. The cluster size distribution was obtained by counting
the number of coarse-grained clusters identified with the DBSCAN algorithm
dbscan ; digregorio2019 ; caporusso2020 ; negro2022hydrodynamic . For panel
(a) 10 different simulations where performed, and for each of them $10^{3}$
configurations where used for sampling. Panel (b) has been obtained using a
single run and $5\times 10^{3}$ configurations, exploiting the shuffling
induced by activity.

###  III.4 The phase diagram

The aforementioned observations allowed us to identify at each Pe and $\phi$
the corresponding phase, and build a complete phase diagram, which is
displayed in Fig. 6 for the particular case $\epsilon=0.3$. The analysis of
the dynamics which will be presented in Sec. III.5 will further support this
classification.

Let us fix $\phi$ and discuss the various transitions found upon increasing Pe
in Fig. 7. At $\phi=10^{-3}$ (very dilute system), and for all ${\rm Pe}$, the
system appears to be always disordered. Upon considering $\phi=0.05$, we
observe a small region where the system is in a gel phase, followed by the
appearance at ${\rm Pe}\sim 10$, or $P^{-1}_{\rm agg}\sim 1$, of phase
separation 333We stress that the numerical values for the dimensionless
numbers that we give here are approximated.. This is due to the activity
acting as an additional noise, which enables particles to rearrange and d to
not get stuck in a metastable conformation. At ${\rm Pe}=150$, or $P^{-1}_{\rm
agg}\sim 12$, activity is high enough so that the phase-separated
configurations start to break, and the disordered phase takes over. At a
higher packing fraction, e.g. $\phi=0.15$, we also observe a transition
between gel and phase-separation, now at a slightly higher ${\rm Pe}$ or
$P^{-1}_{\rm agg}$. We attribute this increase with respect to what was found
for lower $\phi$ to the fact that at higher densities the gel becomes thicker,
and thus more difficult to break. At ${\rm Pe}\sim 120$, or $P^{-1}_{\rm
agg}\sim 10$, we observe a transition between a phase-separated phase and a
percolating network. Again, activity is strong enough to break the single
dense phase, but the density is high enough so that the small clusters are
interconnected, forming the network. At still higher $\phi$, the behavior
remains the same as for $\phi=0.15$, with the ${\rm Pe}$ threshold between gel
and phase separation slightly increasing with density, and the one between
phase separation and percolating network weakly decreasing.

Figure 6: Bulk phase behaviour. Phase diagram in the Pe-$\phi$ plane, at fixed
$\epsilon=0.3$.

We can now explore the effects of the beads attraction strength $\epsilon$.
For both $\epsilon=0.2$ and $\epsilon=0.4$ we find the four phases described
in Fig. 3, and a subdivision similar to Fig. 6 of the regions pertaining to
each phase. However, the critical Pe between these phases are shifted.

In more detail, for $\epsilon=0.2$ the transition between gel and phase-
separated at low $\phi$ and the one between phase-separated and percolating
network at sufficiently high densities take place at smaller Pe than for
$\epsilon=0.3$. This can be attributed to the fact that both activity and
temperature have a higher role in degrading the clustered system as the ratio
Pe$/\epsilon$ is increased. At the extreme case $\epsilon=0$ (not shown) we do
not see ordering nor structure in the system, at any given Pe or $\phi$
considered. The phase diagram here presents only a disordered phase. This
confirms that attraction plays an essential role in the formation of the
aggregates. Conversely, for $\epsilon=0.4$, the aggregation is stronger than
in Fig. 6, and it takes a higher activity to break the clusters. Both gel and
phase-separated phases extend towards slightly higher values of Pe at the same
packing fraction.

Figure 7: Effect of the attraction strength. Phase diagram in the
$1/\phi-\epsilon$ plane for Pe $=100$. The color code is the same as in Fig. 6
with phase-separated in red, disordered in gray and percolating network in
yellow. At this high Pe there is no gel.

The effect of $\epsilon$ can also be visualized in a new phase diagram, on the
$(\phi,1/\epsilon)$ plane, built at fixed $k_{B}T$ and Pe and shown in Fig. 8.
We note that the topology of this phase diagram is very similar to the one in
Fig. 2 of Ref. Zaccarelli2007  derived for a passive system.

###  III.5 Dumbbells dynamics

We now complement the analysis of the structural properties that led us to the
phase diagram in Fig. 6 and Fig. 7 with a characterization of the dynamics. In
particular, we measure the global dumbbells’ mean square displacement (MSD) in
each of the four phases, and the motion of the single clusters formed via
phase separation when a non-vanishing active force is applied.

The total mean square displacement (MSD) is the result of the sum of the MSDs
of the $2N$ colloids composing the $N$ dumbbells. It is defined as

$\Delta^{2}(t,t_{0})=\frac{1}{2N}\sum_{l=1}^{2N}\,\langle(x_{l}(t)-x_{l}(t_{0}))^{2}\rangle\;.$
(8)

Fig. 8 shows $\Delta^{2}$ for four combinations of parameters $\phi$ and
$\rm{Pe}$ representative of the four phases we described beforehand.

In the gel phase, the MSD is sub-diffusive (green curve in Fig. 8(a)). This
can be understood in terms of particles being confined within the gel,
resulting in limited movement. In contrast, the disordered phase and
percolating network (grey and yellow curves in Fig. 8(a)) show diffusive
motion after a brief ballistic regime.

Figure 8: Mean-square displacement (MSD) in the different phases. Left scale:
MSD in the gel (green, Pe = 0 and $\phi=0.1$), disordered phase (gray, Pe
$=200$ and $\phi=0.1$) and percolating network (yellow, Pe $=150$ and
$\phi=0.3$). The percolating network and disordered phases exhibit diffusive
behavior at late time delays, as evidenced by the linear fit (dashed line),
while the dynamics of the gel is much slower and may also tend to freeze.
Right scale: MSD in the phase-separated case (red curve) with Pe = 50 and
$\phi=0.1$. The dotted line represents a quadratic law, suggesting that the
motion of the dumbbell clusters is, on average, ballistic in this case.

When considering the phase-separated system (red curve in Fig. 8(b)), there
is, instead, a ballistic behaviour at long time-scales, characterised by a
quadratic behavior $\Delta^{2}\sim t^{2}$, associated to the motion of a
single aggregate, as we do not see any particles in the dilute phase. This
indicates that the cluster exhibits persistent motion over time, and suggests
to analyze its dynamics in more detail.

##  IV Single cluster motion

We now characterize the ballistic motion of a single cluster in the phase
separated regime. We use a kinematic and dynamical approach to decipher the
interplay between frictional dissipation and active forcing which leads to the
persistent motion of the active cluster.

Figure 9: Approximately helical trajectory of a typical cluster. (a) In red
the trajectory of the center of mass of a cluster with $N_{c}=3728$ beads,
with no surrounding gas, formed by aggregation at $f_{\rm{act}}=1$ and
evolving at $T=0$. The size of the cluster has been reduced by a factor of
$2.5$ to ease the visualization of the trajectory. The scales of the Cartesian
axes are measured in units of the bead’s diameter $\sigma$. The aggregate is,
approximately, an ellipsoid of revolution, or spheroid. The approximate
lengths of the semi-axes of the associated inertia tensor are
$(10.11,\,9.71,\,13.34)$. The principal eigenvector is roughly parallel to
$\hat{\bm{k}}$ – the direction along which the net force momentum is null –
while the remaining eigenvalues associated to the other two eigenvectors are
roughly degenerate, taking values $(17.7,\,16.2,\,37.5)$. In dashed gray the
axis of the helix. The $\hat{\bm{i}},\,\hat{\bm{j}},\,\hat{\bm{k}}$ Cartesian
and inertial coordinate system, superimposed to the cluster, moves along the
helical axis with constant velocity. The unit vectors do not change their
orientation with respect to the laboratory reference frame, with
$\hat{\bm{k}}$ oriented in the direction of the helical axis and
$\hat{{\bm{i}}},\,\hat{{\bm{j}}}$ and perpendicular to it. (b) One period of
the trajectories of the center of mass (dashed line, right scale), and of one
bead in the cluster (blue line, left scale), on the
$\hat{\bm{k}},\,\hat{\bm{i}}$ plane. Two movies in the SM highlight these two
trajectories in the full three-dimensional space.

Figure 9(a) shows a typical trajectory of the center of mass (CoM) of a single
cluster (red curve). This trajectory was obtained by extracting one typical
cluster from the system and placing it in a box with no gas dumbbells around
(in the figure the cluster has been resized to make the trajectory visible).
More details of the motion can be appreciated in the video SM Movie 2. The
cluster moves along a trajectory resembling a helix. Indeed, its motion can be
seen as the composition of a linear motion at constant velocity along the
direction defined by the axis of the helix (dashed grey straight line), and an
additional rotation around this axis444Notice that the radius of the helix is
generally comparable with the beads size, therefore rotational motion around
the persistence axis does not significantly contribute to the MSD, as it will
be rationalized in the following.. Importantly, once the cluster is formed,
the beads do not change their arrangement, so that the cluster effectively
behaves as a rigid body.

We introduce a new reference frame, with its center placed on a fixed point
within the cluster and moving along the axis of the helix with the same
persistence velocity of the cluster itself. The axes of this new coordinate
system are defined by the unit vectors $\hat{\bm{i}}$, $\hat{\bm{j}}$,
$\hat{\bm{k}}$ and are chosen in such a way that $\hat{\bm{k}}$ is the
direction of the helical axis and $\hat{\bm{i}}$, $\hat{\bm{j}}$ are mutually
perpendicular so that the triad forms an orthonormal basis for the three-
dimensional space. This is shown in Fig. 9(a). We stress that such a reference
frame is inertial, as it only translates (and does not rotate) with respect to
the laboratory frame with constant velocity.

In order to rationalize the dynamics of the cluster, we plot in Fig. 9(b) the
CoM trajectory during a single helical period on the
$\hat{\bm{k}},\hat{\bm{i}}$ plane (dashed black line), along with the motion
of a generic bead in the cluster (solid blue curve). Two movies in the SM give
a three-dimensional vision of these motions. By comparing the two, one finds
that their motion is periodic with the same period. However, while the CoM
trajectory has a sinusoidal profile (with a rather small amplitude),
consistently with a circular motion on the $\hat{\bm{i}},\hat{\bm{j}}$ plane
normal to the helix axis, the bead’s motion describes an epicycloid. This is
the fingerprint of the composite dynamics of the beads which rotate around the
CoM while the latter rotates around the helix axis. Such dynamics is
reminiscent of that of the Moon orbiting the Earth and, in turn, rotating on
itself with the same period. Analogously, the cluster rotates around the helix
axis showing always the same face (see SM Movie 3).

We stress that the motion displayed by the cluster in Fig. 9 is not a
particular case; it is consistently observed for any typical cluster in the
phase-separated region of the phase diagram. We have checked that the analysis
that we develop below describes the motion of such clusters with number of
beads ranging from 100 to 4000. Other clusters with more complex forms and
dynamics can also exist (as the result, for example, of the aggregation of two
colliding ones) but we do not discuss them here.

Figure 10: Schematic representation of the forces on the cluster and its
helical and processional motion. (a) Top view, showing the components of
${\bm{F}}_{\mathrm{act}}$ and ${\bm{F}}_{\mathrm{drag}}$ on the
$\hat{\bm{i}},\hat{\bm{j}}$ plane, as well as the radial force ${\bm{F}}_{r}$.
(b) Side view, showing the processional motion of the angular momentum
$\bm{L}$ around the axis $\hat{\bm{k}}$. (c) Time evolution of the $k$ and
modulus of the $ij$ component of the active force $\bm{F}_{\mathrm{act}}$, and
$F_{r}$, compared to the theoretical predictions (dashed lines) from the
equations described in the text. (d) Time evolution of the active torque
$\bm{\mathcal{T}}_{\mathrm{act}}$ and angular momentum $\bm{L}$, also compared
to the theoretical predictions (dashed lines).

To go beyond the kinematic description of the cluster’s motion, we now solve
its dynamics building on the assumption that the cluster is, for all practical
effects, a rigid body. Therefore, the motion can be described by separating
positional and rotational degrees of freedom.

The CoM motion is governed by the effects of the total active force and
friction. Newton’s equation for the CoM position ${\bm{R}}_{\rm cm}$ reads

$M\ddot{\bm{R}}_{\rm cm}=\bm{F}_{\rm drag}+\bm{F}_{\text{act}}\;,$ (9)

with $M=mN_{c}$ the total mass of the cluster, $N_{c}$ the number of particles
in the cluster, ${\bm{R}}_{\rm cm}=N_{c}^{-1}\sum_{l=1}^{N_{c}}{\bm{r}}_{l}$
the position of the CoM,
$\bm{F}_{\text{act}}=\sum_{l=1}^{N_{c}}{\bm{f}_{\text{act}}}_{l}$ the total
active force, and $\bm{F}_{\rm drag}=-(M/m)\gamma\dot{\bm{R}}_{\rm cm}$ the
total drag acting on the center of mass. We stress that internal forces,
arising from the Mie potential and bond constraints, do not contribute to the
dynamics of the center of mass.

The equations describing the rotational dynamics are given by

$\dot{\bm{L}}=\bm{\mathcal{T}}_{\rm drag}+\bm{\mathcal{T}}_{\rm act}\;.$ (10)

Here
${\bm{L}}=\sum_{l=1}^{N_{c}}\bm{r}^{\prime}_{l}\times\dot{\bm{r}}^{\prime}_{l}$
is the angular momentum of the cluster computed choosing the CoM as pole, with
$\bm{r}^{\prime}_{l}=\bm{r}_{l}-\bm{R}_{\rm cm}$ the position of the $l-$th
bead with respect to the CoM. The right-hand-side is the total force momentum,
decomposed in a contribution from the active force $\bm{\mathcal{T}}_{\rm
act}=\sum_{l=1}^{N_{c}}\bm{r}^{\prime}_{l}\times{\bm{f}_{\rm act}}_{l}$ and
another one from the drag force $\bm{\mathcal{T}}_{\rm drag}=-\gamma{\bm{L}}$.
Once more, torques originating from internal forces do not contribute to the
overall rotational dynamics.

Before proceeding, we stress that a peculiarity of this system is that the
active forces $\bm{F}_{\text{act}}$ and momenta $\bm{\mathcal{T}}_{\rm act}$
are solidal to the body itself, so that they undergo the same rotational and
translational motion as the rigid body, as beads do not reposition during the
trajectory.

With the equations of motion at hand, we are now ready to derive an expression
relating forces and momenta with kinematic quantities. First, we observe that
the CoM acceleration is null along the $\hat{\bm{k}}$ direction, therefore the
$k$-th component of the active force is counterbalanced by the drag:

$F_{\rm{act}}^{k}=\dfrac{M}{m}\gamma V^{k}\;,$ (11)

where $\bm{V}=\dot{\bm{R}}_{\rm cm}$, so that, by construction, $V^{k}$ is the
persistence velocity of the cluster along the helix axis, measured in the
laboratory reference frame.

Conversely, in the $\hat{\bm{i}},\,\hat{\bm{j}}$ plane, the CoM velocity
rotates with constant angular velocity $\bm{\omega}=\omega^{k}\hat{\bm{k}}$,
see Fig. 10(a). In the overdamped regime, we expect that the active force is
approximately counterbalanced by the frictional force leading to
$\bm{F}_{\rm{act}}^{ij}=M\gamma\bm{V}^{ij}/m$, except from a small
contribution, directed radially with respect to the circle, computed as
$\bm{F}_{r}=\bm{F}^{ij}_{\rm drag}+\bm{F}^{ij}_{\text{act}}$, which is
ultimately responsible for the rotation of the cluster around the $k-$axis
(and is equal in fact to $M\ddot{\bm{R}}_{\rm cm}$). Moreover, because the
circular motion occurs with uniform velocity $V^{ij}=\omega^{k}R$, with $R$
the radius of the circle, one can write these forces in terms of the
rotational velocity, as follows:

$\displaystyle F_{\rm{act}}^{ij}$
$\displaystyle=\dfrac{M}{m}R\gamma\omega^{k}\;,$ (12) $\displaystyle F_{r}$
$\displaystyle=MR{(\omega^{k})}^{2}\;.$ (13)

(Note that in Eq. (12) we have dropped the inertial contribution, consistently
with the observation that $\omega^{k}$ does not change in time.)

We now test whether the formulæ (11)-(13) are in agreement with the simulation
data. In Fig. 10(c), we show the forces considered (continuous lines) as
functions of time, and we compare them to the corresponding quantities
appearing on the right-hand-sides of Eqs. (11)-(13) (dashed lines). We find
good agreement in all cases. In particular, we notice that the $\hat{\bm{k}}$
direction roughly coincides with that defined by the principal eigenvector of
the inertia tensor $\bm{I}$ while, instead, the radial force $F_{r}$ is
significantly smaller than the other ones. Moreover, the forces are
hierarchically organized, with $F_{\rm act}^{k}\sim 100$, $F_{\rm
act}^{ij}\sim 20$, and $F_{\rm r}\sim 10^{-2}$, in units of $f_{\rm act}$.

We now proceed with the analysis of the rotational dynamics by assuming that
the angular momentum component $L^{k}$ in the $\hat{\bm{k}}$ direction remains
constant over time – or equivalently drag and active force momenta compensate
along $\hat{\bm{k}}$, i.e. $\gamma L^{k}=\mathcal{T}_{\rm act}^{k}$. The
normal components on the $\hat{\bm{i}},\hat{\bm{j}}$ plane are constant in
modulus, and rotate with angular velocity $\bm{\omega}$ around the helix axis
under the action of the total force momentum $\bm{\mathcal{T}}_{\rm
tot}=\bm{\mathcal{T}}_{\rm drag}+\bm{\mathcal{T}}_{\rm act}$. Under these
assumptions, it is straightforward to rewrite the time derivative in Eq. (10)
so that

$\bm{\omega}\times\bm{L}=\bm{\mathcal{T}}_{\rm tot}\;.$ (14)

This expression suggests that $\bm{L}$ performs a precessional motion around
the helical axis under the action of transversal force momenta, see also Fig.
10(a). Moreover, one can also relate the angular momentum along $\hat{\bm{k}}$
with the $kk$ component of the inertia tensor so that

$L^{k}=I^{kk}\omega^{k}\;,$ (15)

being $I^{kk}$ the helical axis component of the inertia tensor $\bm{I}$ of
the cluster.

These equations are tested against simulations in Fig. 10(d), and we find
again good agreement. This confirms that the angular momentum performs a
precession around the $\hat{\bm{k}}$ axis with precession angular velocity
equal to that of the cluster, confirming _a posteriori_ our initial
assumption.

As further tests, we checked the dependence of the period on the magnitude of
the active force $f_{\rm act}$ derived from Eq. (12),

$\tau=\frac{2\pi M\gamma R}{\kappa m}\,\frac{1}{f_{\rm{act}}}\;,$ (16)

where $F_{\rm{act}}^{ij}=\kappa f_{\rm{act}}$ with $\kappa$ a proportionality
coefficient measured from simulation data that turns out to be independent of
the active force $f_{\rm act}$. To this purpose, we considered the same
cluster and we varied $f_{\rm act}$. The results of this analysis show that
the period and the active force are indeed inversely proportional, see Fig.
11(a), with the proportionality coefficient $\kappa$ depending on the
structural properties of each individual cluster considered, and being
consistent with the parameter dependencies in Eq. (12)

Figure 11: Motion of a cluster with mass $M=3728$ at temperature $T=0$. (a)
Rotation period $\tau$ as a function of the active force $f_{\rm act}$ on a
single bead. The slashed straight line represents the theoretical prediction
$\tau\propto f_{\rm act}^{-1}$ in Eq. (16), with $\kappa=48.71$ and $R=0.204$
measured independently. (b) Linear velocity of the CoM, $V^{k}$, as a function
of the single bead active force $f_{\rm act}$. The dashed line is the
theoretical prediction inverting Eq. (11), $V^{k}=m\ell/(M\gamma)f_{\rm act}$,
using $\ell=F^{k}_{\rm act}/f_{\rm act}\approx 96.34$ (measured independently)
and $\gamma=10$.

At the same time, we find that the radius $R$ remains constant while we vary
$f_{\rm act}$ (not shown). In fact, inverting Eq. (12), one finds that
$R=F_{\rm{act}}^{ij}m/(M\gamma\omega^{k})$, with both $F_{\rm{act}}^{ij}$ and
$\omega^{k}$ proportional to $f_{\rm{act}}$, leaving no dependence on the
active force.

Finally, we put Eq. (11) to the numerical test. We show in Fig. 11(b) that
$V^{k}$ is proportional to $f_{\rm act}$ with a proportionality constant which
agrees with the parameter dependence of the prefactor in Eq. (11).

We also checked that the direction of motion ${\bf\hat{k}}$ does not change in
time.

We stress that the description of the cluster’s dynamics is based on the
assumption that the rotational velocity $\bm{\omega}$ is oriented parallel to
the helical axis, with negligible transversal components. The latter could
lead in principle to more complex effects – such as nutation for instance –
which, nevertheless, were not observed in simulations. This also resonates
with the fact that in absence of net torques, a rigid body with two degenerate
eigenvalues in its inertia tensor would rotate by keeping constant the
component of the angular momentum parallel to the non-degenerate eigenvector
Goldstein . Interestingly enough, in our case the particle aggregates indeed
exhibit elongated shapes with the principle eigenvector roughly parallel to
$\hat{\bm{k}}$ and the two dimensions being roughly similar, see the values
given in the caption of Fig. 9. Nevertheless, the motion of the cluster
considered here is in general significantly more complicated as, in this case,
the net force and torque are not null, leading to a more complex dynamics,
with the radial force sustaining the rotational dynamics of the CoM of the
cluster, ultimately leading to the intriguing helical motion described above.

##  V Conclusions

To conclude, we investigated the nature of the phase behavior and the dynamics
of an active dumbbell system with attractive interactions in three dimensions.
We characterized the phase diagram in the (Pe, $\phi\leq 0.3$) plane at fixed
$\epsilon=0.3$ and we elucidated the effect of the strength of the attraction
$\epsilon$ by studying the stationary state reached at Pe = 100 and parameters
in the ($\phi\leq 0.3,\,1/\epsilon\leq 5$) plane. In this way we showed that
four cases are realized: a disordered state, a percolating network, phase
separation and an active gel.

Next, we focused on characterising the motion of a typical dense cluster in
the phase-separated phase. First, we found that these clusters typically take
a spheroid form, and for $f_{\rm act}<10$ displace with constant velocity in a
direction which is very close to its main axis of symmetry while performing a
rotational motion of very small radius in the transverse plane. All in all,
the motion is very close to helical. With some simple arguments, explained in
Sec. IV, we then related the linear and angular velocities to the strength of
the active force acting on the single molecules finding very good agreement
with the numerical measurements. The dynamics of formation of these clusters,
a full characterization of their morphology and statistics, and many other
details are very interesting but fall beyond the scope of this work.

## Acknowledgements

This numerical work was carried out on the Dutch national e-infrastructure
with the support of SURF through Grant 2021.028 for computational time (L. N.
C. and G. N.) and on ReCas HPC-Cluster in Bari (Italy). We acknowledge funding
from MIUR Project No. PRIN 2020/PFCXPE and ANR-20-CE30-0031. We thank P.
Royall and F. Turci for very useful discussions.

## References

  * [1] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha.  Hydrodynamics of soft active matter.  Rev. Mod. Phys., 85:1143, 2013.
  * [2] C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe.  Active particles in complex and crowded environments.  Rev. Mod. Phys., 88:045006, 2016.
  * [3] S. Shankar, A. Souslov, M. J. Bowick, M. C. Marchetti, and V. Vitelli.  Topological active matter.  Nature Reviews Physics, 4:380–398, 2022.
  * [4] L. N. Carenza, G. Gonnella, A. Lamura, G. Negro, and A. Tiribocchi.  Lattice boltzmann methods and active fluids.  Eur. Phys. J. E, 42:81, 2019.
  * [5] A. Walther and A. H. E. Müller.  Janus Particles: Synthesis, Self-Assembly, Physical Properties, and Applications.  Chem. Rev., 113:5194–5261, 07 2013.
  * [6] S. Das, A. Garg, A. I. Campbell, J. Howse, A. Sen, D. Velegol, R. Golestanian, and S. J. Ebbens.  Boundaries can steer active janus spheres.  Nature Comm., 6:8999, 2015.
  * [7] S. J. Ebbens.  Active colloids: Progress and challenges towards realising autonomous applications.  Curr. Op. Coll. Int. Sci., 21:14, 2016.
  * [8] T. Sanchez, D. T. N. Chen, S. J. DeCamp, M. Heymann, and Z. Dogic.  Spontaneous motion in hierarchically assembled active matter.  Nature, 491:431, 2012.
  * [9] J. Howard.  Mechanics of Motor Proteins and the Cytoskeleton.  Sinauer Associates, Publishers, 2001.
  * [10] T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans, and B. Ladoux.  Topological defects in epithelia govern cell death and extrusion.  Nature, 544:212–216, 2017.
  * [11] P. Guillamat, C. Blanch-Mercader, G. Pernollet, K. Kruse, and A. Roux.  Integer topological defects organize stresses driving tissue morphogenesis.  Nat. Mat., 21:588, 2022.
  * [12] J.-M. Armengol-Collado, L. N. Carenza, J. Eckert, D. Krommydas, and L. Giomi.  Epithelia are multiscale active liquid crystals.  Nat. Phys., page (in press), 2023.
  * [13] A. Haeger, S. Alexander, M. Vullings, F. M.P. Kaiser, C. Veelken, U. Flucke, G. E. Koehl, M. Hirschberg, M. Flentje, R. M. Hoffman, E. K. Geissler, S. Kissler, and P. Friedl.  Collective cancer invasion forms an integrin-dependent radioresistant niche.  Journal of Experimental Medicine, 217, 10 2019.  e20181184.
  * [14] Y. Maroudas-Sacks, L. Garion, L. Shani-Zerbib, A. Livshits, E. Braun, and K. Keren.  Topological defects in the nematic order of actin fibers as organization centers of Hydra morphogenesis.  Nat. Phys., 17:251, 2021.
  * [15] F. Brauns, N. H. Claussen, E. F. Wieschaus, and B. I. Shraiman.  Epithelial flow by controlled transformation of internal force-balance geometry.  bioRxiv, 2023.
  * [16] Y. Hatwalne, S. Ramaswamy, M. Rao, and R. A. Simha.  Rheology of active-particle suspensions.  Phys. Rev. Lett., 92:118101, 2004.
  * [17] K. Kruse, J. F. Joanny, F. Jülicher, J. Prost, and K. Sekimoto.  Asters, vortices, and rotating spirals in active gels of polar filaments.  Phys. Rev. Lett., 92:078101, Feb 2004.
  * [18] M. E. Cates and J. Tailleur.  Motility-induced phase separation.  Annu. Rev. Condens. Matter Phys., 6:219, 2015.
  * [19] V. Bratanov, F. Jenko, and E. Frey.  New class of turbulence in active fluids.  Proc. Nat. Acad. Sci., 112:15048, 2015.
  * [20] S. Bi, F. Jin, and V. Sourjik.  Inverted signaling by bacterial chemotaxis receptors.  Nature Comm., 9:2927, 2018.
  * [21] H. Massana-Cid, C. Maggi, G. Frangipane, and R. Di Leonardo.  Rectification and confinement of photokinetic bacteria in an optical feedback loop.  Nature Comm., 13:2740, 2022.
  * [22] Y. I. Yaman, E. Demir, R. Vetter, and A. Kocabas.  Emergence of active nematics in chaining bacterial biofilms.  Nature Comm., 10:2285, 2019.
  * [23] M. Basaran, Y. I. Yaman, T. C. Yüce, R. Vetter, and A. Kocabas.  Large-scale orientational order in bacterial colonies during inward growth.  eLife, 11:e72187, mar 2022.
  * [24] M. Basaran, T.C. Yuce, Y. I. Yaman, R. Vetter, and A. Kocabas.  Bacterial biofilms use chiral branches to escape crowded environments by tracking oxygen gradient.  arXiv:2208.09730, 2022.
  * [25] J. Dunkel, S. Heidenreich, K. Drescher, H. H. Wensink, M. Bär, and R. E. Goldstein.  Fluid dynamics of bacterial turbulence.  Phys. Rev. Lett., 110:228102, 2013.
  * [26] L. Giomi.  Geometry and topology of turbulence in active nematics.  Phys. Rev. X, 5, 2015.
  * [27] L. N. Carenza, L. Biferale, and G. Gonnella.  Multiscale control of active emulsion dynamics.  Phys. Rev. Fluids, 5:011302, 2020.
  * [28] L. N. Carenza, L. Biferale, and G. Gonnella.  Cascade or not cascade? Energy Transfer and Elastic Effects in Active Nematics.  EPL, 132:44003, 2020.
  * [29] I. S. Aranson.  Bacterial active matter.  Reports on Progress in Physics, 85:076601, 2022.
  * [30] Ricard Alert, Jaume Casademunt, and Jean-François Joanny.  Active turbulence.  Annu. Rev. Condens. Matter Phys., 13(1):143–170, 2022.
  * [31] E. Locatelli, F. Baldovin, E. Orlandini, and M. Pierno.  Active brownian particles escaping a channel in single file.  Phys. Rev. E, 91:022109, Feb 2015.
  * [32] A. P. Solon, M. E. Cates, and J. Tailleur.  Active brownian particles and run-and-tumble particles: A comparative study.  The European Physical Journal Special Topics, 224(7):1231–1262, 2015.
  * [33] M. E. Cates.  Active Field Theories.  In Active Matter and Nonequilibrium Statistical Physics. Oxford University Press, 11 2022.
  * [34] R. Wittkowski, A. Tiribocchi, J. Stenhammar, R. J. Allen, D. Marenduzzo, and M.E. Cates.  Scalar $\phi^{4}$ field theory for active-particle phase separation.  Nat. Comm., 5, 2014.
  * [35] D. Dell’Arciprete, M. L. Blow, A. T. Brown, F. D. C. Farrell, J. S. Lintuvuori, A. F. McVey, D. Marenduzzo, and W. C. K. Poon.  A growing bacterial colony in two dimensions as an active nematic.  Nature Communications, 9(1):4190, 2018.
  * [36] G. Negro, A. Lamura, G. Gonnella, and D. Marenduzzo.  Hydrodynamics of contraction-based motility in a compressible active fluid.  EPL, page 58001, 2018.
  * [37] M. G. Giordano, F. Bonelli, L. N. Carenza, G. Gonnella, and G. Negro.  Activity-induced isotropic-polar transition in active liquid crystals.  EPL, 133(5):58004, 2021.
  * [38] L. N. Carenza, G. Gonnella, A. Lamura, D. Marenduzzo, G. Negro, and A. Tiribocchi.  Soft channel formation and symmetry breaking in exotic active emulsions.  Scientific Reports, 10(1):15936, 2020.
  * [39] I. Favuzzi, L.N. Carenza, F. Corberi, G. Gonnella, A. Lamura, and G. Negro.  Rheology of active emulsions with negative effective viscosity.  Soft Materials, 19(3):334–345, 2021.
  * [40] L.N. Carenza, G. Gonnella, A. Lamura, and G. Negro.  Dynamically asymmetric and bicontinuous morphologies in active emulsions.  Int. J. Mod. Phys. C, 30(10):1941002, 2019.
  * [41] P. Digregorio, D. Levis, L. F. Cugliandolo, G. Gonnella, and I. Pagonabarraga.  Clustering of topological defects in two-dimensional melting of active and passive disks.  Soft Matter, 18:566, 2022.
  * [42] C. B. Caporusso, P. Digregorio, D. Levis, L. F. Cugliandolo, and G. Gonnella.  Motility-induced microphase and macrophase separation in a two-dimensional active Brownian particle system.  Phys. Rev. Lett., 125:178004, 2020.
  * [43] C. Caporusso et al., 2023.  in preparation.
  * [44] F. Cagnetta, F. Corberi, G. Gonnella, and A. Suma.  Large fluctuations and dynamic phase transition in a system of self-propelled particles.  Phys. Rev. Lett., 119:158002, 2017.
  * [45] M. Bär, R. Großmann, S. Heidenreich, and F. Peruani.  Self-propelled rods: Insights and perspectives for active matter.  Annual Rev. of Condens. Matt. Phys., 11:441, 2020.
  * [46] E. Zaccarelli, P. J. Lu, F. Ciulla, D. A. Weitz, and F. Sciortino.  Gelation as arrested phase separation in short-ranged attractive colloid-polymer mixtures.  J. Phys.: Condens. Matter, 20:494242, 2008.
  * [47] B. I. Halperin and David R. Nelson.  Theory of two-dimensional melting.  Phys. Rev. Lett., 41:121–124, Jul 1978.
  * [48] C. B. Caporusso, L. F. Cugliandolo, P. Digregorio, G. Gonnella, D. Levis, and A. Suma.  Dynamics of motility-induced clusters: coarsening beyond ostwald ripening.  arXiv:2211.12361, 2022.
  * [49] L. F. Cugliandolo, P. Digregorio, G. Gonnella, and A. Suma.  Phase coexistence in two-dimensional passive and active dumbbell systems.  Phys. Rev. Lett., 119:268002, 2017.
  * [50] P. Digregorio, D. Levis, A. Suma, L. F. Cugliandolo, G. Gonnella, and I. Pagonabarraga.  Full phase diagram of active Brownian disks: From melting to motility-induced phase separation.  Phys. Rev. Lett., 121:098003, 2018.
  * [51] G. Gonnella, D. Marenduzzo, A. Suma, and A. Tiribocchi.  Motility-induced phase separation and coarsening in active matter.  Comptes Rendus Physique, 16:316, 2015.
  * [52] A. Suma, G. Gonnella, D. Marenduzzo, and E. Orlandini.  Motility-induced phase separation in an active dumbbell fluid.  EPL, 108:56004, 2014.
  * [53] A. Suma, G. Gonnella, G. Laghezza, A. Lamura, A. Mossa, and L. F. Cugliandolo.  Dynamics of a homogeneous active dumbbell system.  Phys. Rev. E, 90:052130, 2014.
  * [54] I. Petrelli, P. Digregorio, L. F. Cugliandolo, G. Gonnella, and A. Suma.  Active dumbbells: Dynamics and morphology in the coexisting region.  Eur. Phys. J. E Soft Matter, 41:128, 2018.
  * [55] J. Stenhammar, D. Marenduzzo, R. J. Allen, and M. E. Cates.  Phase behaviour of active Brownian particles: the role of dimensionality.  Soft Matter, 10:1489, 2014.
  * [56] A. Wysocki, R. G. Winkler, and G. Gompper.  Cooperative motion of active Brownian spheres in three-dimensional dense suspensions.  EPL, 105:48004, 2014.
  * [57] R. van Damme, J. Rodenburg, R. Van Roij, and M. Dijkstra.  Interparticle torques suppress motility-induced phase separation for rodlike particles.  J. Chem. Phys., 150:164501, 2019.
  * [58] P. Nie, J. Chattoraj, A. Piscitelli, P. Doyle, R. Ni, and M. Pica Ciamarra.  Stability phase diagram of active Brownian particles.  Phys. Rev. Research, 2:023010, 2020.
  * [59] A. K. Omar, K. Klymko, T. GrandPre, and P. L. Geissler.  Phase diagram of active Brownian spheres: Crystallization and the metastability of motility-induced phase separation.  Phys. Rev. Lett., 126:188002, 2021.
  * [60] F. Turci and N. B. Wilding.  Phase separation and multibody effects in three-dimensional active Brownian particles.  Phys. Rev. Lett., 126:038002, 2021.
  * [61] N. Venkatareddy, S.-T. Lin, and K. Prabal.  Phase behavior of active and passive dumbbells.  Phys. Rev. E, 107:034607, 2023.
  * [62] B. M. Mognetti, A. Šarić, S. Angioletti-Uberti, A. Cacciuto, C. Valeriani, and D. Frenkel.  Living clusters and crystals from low-density suspensions of active colloids.  Phys. Rev. Lett., 111:245702, 2013.
  * [63] J. Schwarz-Linek, C. Valeriani, A. Cacciuto, M. E. Cates, D. Marenduzzo, A. N. Morozov, and W. C. K. Poon.  Phase separation and rotor self-assembly in active particle suspensions.  Proc. Nat. Acad. Sc., 109:4052, 2012.
  * [64] G. Melaugh, D. Marenduzzo, A. Morozov, and R. J. Allen.  A simulation study of aggregation mediated by production of cohesive molecules.  Soft Matter, 15:9120, 2019.
  * [65] V. Prymidis, H. Sielcken, and L. Filion.  Self-assembly of active attractive spheres.  Soft Matter, 11:4158, 2015.
  * [66] E. Zaccarelli.  Colloidal gels: equilibrium and non-equilibrium routes.  J. Phys.: Condens. Matter, 19:323101, 2007.
  * [67] D. J. Kraft, R. Ni, F. Smallenburg, D. A. Weitz M. Hermes, K. Yoon, A. van Blaaderen, J. Gronewold, M. Dijkstra, and W. K. Kegel.  Surface roughness directed self-assembly of patchy particles into colloidal micelles.  Proc. Nat. Acad. Sc., 109:10787, 2012.
  * [68] C. Vega, E. P. A. Paras, and P. A. Monson.  Solid–fluid equilibria for hard dumbbells via monte carlo simulation.  J. Chem. Phys., 96:9060, 1992.
  * [69] M. Marechal and M. Dijkstra.  Stability of orientationally disordered crystal structures of colloidal hard dumbbells.  Phys. Rev. E, 77:061405, 2008.
  * [70] M. Kowalik and K. W. Wojciechowski.  The free energy of hard dimer solids revisited.  J. Non-Cryst. Solids, 354:4354, 2008.
  * [71] M. Dennison, K. Milinković, and M. Dijkstra.  Phase diagram of hard snowman-shaped particles.  J. Chem. Phys., 137:044507, 2012.
  * [72] S. Whitelam and R. L. Jack.  The statistical mechanics of dynamic pathways to self-assembly.  Annu. Rev. Phys. Chem., 66:143, 2015.
  * [73] S. Whitelam and S. A. F. Bon.  Self-assembly of amphiphilic peanut-shaped nanoparticles.  J. Chem. Phys., 132:074901, 2010.
  * [74] G. Munaò, D. Costa, A. Giacometti, C. Caccamo, and F. Sciortino.  Structure and phase behavior of colloidal dumbbells with tunable attractive interactions.  Phys. Chem. Chem. Phys., 15:20590, 2013.
  * [75] G. Munaò, P. O’Toole, T. S. Hudson, D. Costa, C. Caccamo, A. Giacometti, and F. Sciortino.  Phase separation and self-assembly of colloidal dimers with tunable attractive strength: from symmetrical square-wells to janus dumbbells.  Soft Matter, 10:5269, 2014.
  * [76] G. Munaò, P. O’Toole, T. S. Hudson, D. Costa, C. Caccamo, A. Giacometti, and F. Sciortino.  Cluster formation and phase separation in heteronuclear janus dumbbells.  J. Phys. Condens. Matt., 27:234101, 2015.
  * [77] G. Avvisati, T. Vissers, and M. Dijkstra.  Self-assembly of patchy colloidal dumbbells.  J. Chem. Phys., 142:084905, 2015.
  * [78] P. O’Toole, A. Giacometti, and T. Hudson.  Phase diagram of heteronuclear janus dumbbells.  Soft Matter, 13:803, 2017.
  * [79] C. P. Royall, P. Charbonneau, M. Dijkstra, J. Russo, F. Smallenburg, T. Speck, and C. Valeriani.  Colloidal hard spheres: Triumphs, challenges and mysteries.  arXiv:2305.02452.
  * [80] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton.  LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales.  Comp. Phys. Comm., 271:108171, 2022.
  * [81] J.-P. Ryckaert, G. Ciccotti, and H. J. C. Berendsen.  Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes.  J. Comp. Phys., 23:327–341, 1977.
  * [82] M. Ester, H.-P. Kriegel, J. Sander, and X. Xu.  A density-based algorithm for discovering clusters in large spatial databases with noise.  page 226, 1996.
  * [83] W. C. K. Poon.  The physics of a model colloid-polymer mixture.  J. Phys.: Condens. Matter, 14:R859, 2002.
  * [84] A. Coniglio, L. De Arcangelis, E. Del Gado, A. Fierro, and N. Sator.  Percolation, gelation and dynamical behaviour in colloids.  J. Phys.: Condens. Matter, 16:S4831, 2004.
  * [85] P. J. Lu, E. Zaccarelli, F. Ciulla, A. B. Schofield, F. Sciortino, and D. A. Weitz.  Gelation of particles with short-range attraction.  Nature, 453:499, 2008.
  * [86] Supplemental movies can be found at https://www.dropbox.com/scl/fo/julh85npfkhv7rw3tke5r/h?rlkey=botjlcrqs9wpxgtiz2czfmtia&dl=0.
  * [87] A. Okabe, B. Boots, K. Sugihara, and S. N. Chiu.  Spatial Tessellations – Concepts and Applications of Voronoi Diagrams (2nd ed.).  John Wiley, 2000.
  * [88] D. Faken and H. Jónsson.  Systematic analysis of local atomic structure combined with 3d computer graphics.  Computational Materials Science, 2:279, 1994.
  * [89] M. Marechal and M. Dijkstra.  Stability of orientationally disordered crystal structures of colloidal hard dumbbells.  Phys. Rev. E, 77:061405, 2008.
  * [90] D. Stauffer and A. Aharony.  Introduction To Percolation Theory: Second Edition.  Taylor & Francis, London, UK, 2017.
  * [91] G. Negro, C.B. Caporusso, P. Digregorio, G. Gonnella, A. Lamura, and A. Suma.  Hydrodynamic effects on the liquid-hexatic transition of active colloids.  Eur Phys J E Soft Matter, 45(9):75, 2022.
  * [92] H. Goldstein.  Classical Mechanics.  Addison-Wesley, 1959.