11institutetext:  Dipartimento di Fisica, Università degli Studi di Bari and
INFN, Sezione di Bari, via Amendola 173, Bari, I-70126, Italy 22institutetext:
Centre Européen de Calcul Atomique et Moléculaire (CECAM), Ecole Polytechnique
Fédérale de Lausanne (EPFL), Batochimie, Avenue Forel 2, 1015 Lausanne,
Switzerland 33institutetext: Istituto Applicazioni Calcolo, CNR, Via Amendola
122/D, I-70126 Bari, Italy

We study numerically the role of hydrodynamics in the liquid-hexatic
transition of active colloids at intermediate activity, where motility induced
phase separation (MIPS) does not occur. We show that in the case of active
Brownian particles (ABP), the critical density of the transition decreases
upon increasing the particle’s mass, enhancing ordering, while self-propulsion
has the opposite effect in the activity regime considered. Active hydrodynamic
particles (AHP), instead, undergo the liquid-hexatic transition at higher
values of packing fraction $\phi$ than the corresponding ABP, suggesting that
hydrodynamics have the net effect of disordering the system. At increasing
densities, close to the hexatic-liquid transition, we found in the case of AHP
the appearance of self-sustained organized motion with clusters of particles
moving coherently.

# Hydrodynamic effects on the liquid-hexatic transition of active colloids

G. Negroaddr1  C. B. Caporussoaddr1  P. Digregorioaddr2  G. Gonnellaaddr1  A.
Lamuraaddr3  A. Sumaaddr1

(Received: date / Accepted: date)

††journal: Eur. Phys. J. E

##  1 Introduction

Self-propelled particles (SPP) are the fundamental units of a broad class of
theoretical models for active matter. In the context of SPP models, injected
energy from the environment fuels a persistent motion of the single
constituents, driving the system out of thermal equilibrium. Simplified models
of SPP filyABPs ; romanczuk2012 ; speckCRYSTALSPP ; bechinger-revmodphys  are
of crucial importance, because they offer a minimal setup to explore some of
the large variety of collective behaviours observed in nature for systems of
motile living bodies at different length scales, from flocking of birds and
fish bialekFLOCKS , to swarming in bacterial colonies beerBACTERIA  and
dynamics in cells’ cytoskeleton dogicCELLS .

Active Brownian Particles (ABP) models are very popular among SPP models
nostroDISCHI ; claudioDISCHI ; filyABPs . Active colloids are usually
spherical particles undergoing directed motion due to an active force, while
both translational and rotational degrees of freedom are in contact with a
stochastic thermal bath. Although the model is very simple, ABP show
paradigmatic collective phenomena like motility-induced phase separation
(MIPS) catesMIPS ; filyABPs ; gonnellaMIPS ; rednerMIPS ; buttSPP  and are
therefore very interesting in order to characterize the fundamental principles
governing active matter systems. Moreover, ABP are of primary use for
comparisons with experimental systems of synthetic micro-swimmers buttSPP ;
demianSPP , opening the perspectives for a systematic control of active
systems and collective motion, with the purpose to exploit some of their
unique features for technological uses, for instance in robotics robotswarm ;
smarticles ; robotlocomotion , realisation of biological machines
kapral_nanomotors , or understanding of flocking intelligence PNAS_schooling ;
biderale_swimming-strategies .

Of particular interest is the characterization of ABP in the dense regime, see
e.g. spontaneous flow yeomansDENSE  or glassy behaviour berthierDENSE ;
mandalDENSE  in biological tissues, biofilms, cell mono-layers henkesDENSE ;
biACTGLASS , and can be considered a target for the development of new
materials smart_materials_1 . In two dimensions (2D), ABP present ordering
phase transitions when the density of the system is increased krauthACTIVE ;
nostroDISCHI ; claudioDISCHI ; nostroDIFETTI , which are connected to those
encountered for passive hard colloids krauthPASSIVE1 ; krauthPASSIVE2 ;
thorneyworkACTIVEMELTING . At intermediate values of the self-propelling
force, a liquid-hexatic critical transition is followed by a hexatic-solid
transition, where the solid phase has quasi-long-range (QLR) positional and
long-range (LR) orientational order, the hexatic phase has short-range (SR)
positional and QLR orientational order, while the liquid phase is homogeneous
and has SR positional and orientational order. This scenario is very similar
to the theoretical Kosterlitz, Thouless, Halperin, Nelson, and Young (KTHNY)
two-step scenarioKT_1 ; KT_2 ; KT_3 . If activity is high enough, instead,
MIPS takes place, as a phase separation between a dense phase and a gaseous
one nostroDISCHI .

The aforementioned features of the ABP phase diagram have been well
established in the context of over-damped motion and without an explicit
underlying thermo-hydrodynamic bath. At the same time, there are other
interesting questions that remain to be considered. The first question
concerns the role of particles mass, and in particular the interplay between
inertial and active diffusion timescales, which can be varied independently
pica2020 ; capriniINERTIA . It has been pointed out in omar2021tuning  that in
three-dimensional active systems, inertia should attenuate the destabilizing
effect of activity on the ordered phase. The presence of large inertia has
also been shown to strongly affect the kinetic energy of the particles into
the highly dense phase of MIPS isabellaTEFF , and to highly inhibit phase
segregation lowenINERTIA . However, the role of inertia in the context of
dense ABP, and in particular how the particle’s mass affects the hexatic
phase, has not yet been characterized. The second question concerns the role
of hydrodynamic interactions in the dense phase. Regarding the influence of
hydrodynamics in MIPS, it is found that in 2D MIPS is suppressed theersSRD ;
navarro2014 ; gompperROADMAP , as hydrodynamics favour reorientation of
particles’ self propulsion direction, while in quasi-2D systems MIPS has been
observed for low-density fluidsstark  and not when the fluid was made
incompressible thee2018 ; gompperROADMAP . For elongated colloids, steric
alignment and hydrodynamics show highly non-trivial interplay, such that MIPS
is enhanced for pullers and suppressed for pushers theersSRD .

As a first step in the direction of answering these two questions, we
characterize how the critical density for the liquid-hexatic transition of
active particles is modified, in an intermediate activity regime where MIPS
does not occur for ABP, by i) the inertial effects due to mass changes, and
ii) the presence of non-isotropic interactions between colloids introduced by
hydrodynamics. Hydrodynamics has been implemented by using the multi-particle
collision method male1999 ; gommper_review_2009 , which seamlessly integrate
with the dynamics of active Brownian particles plipton_colloids_lammps . In
particular, we implement thermal slip boundary conditions, decoupling colloids
rotational diffusion from the solvent and test the consistency of this
implementation with known benchmark tests. We focus here only on 2D systems
where the rotational diffusion follows the same equations as for ABP. This
allows us to have an active hydrodynamic particle (AHP) model with the same
friction, temperature and rotational diffusion as the ABP model, providing a
way to quantitatively compare them.

We find that changing the colloids mass and introducing hydrodynamic
interactions affect the critical density at which the liquid-hexatic
transition occurs. In particular, mass changes lower this density with respect
to over-damped ABP, while hydrodynamics increases the critical density. We
also find that the system with hydrodynamics undergoes a transition from a
disorganized to a self-sustained flow regime upon increasing the density, with
particles moving on the same direction at high densities.

The work is organized in the following way. In Sec. II we discuss the
numerical methods and parameter choice for the ABP model and for the AHP
model, with Sec. II C providing several tests for implementation of the latter
model. In Sec. III A we discuss how the liquid-hexatic scenario changes by
varying the active colloids mass, while in Sec. III B we discuss the effects
due to hydrodynamics interactions. Finally we draw some conclusions discussing
the main findings.

##  2 Numerical methods

In this Section we describe the numerical models. We will start with the ABP
model, which follows a Langevin equation and does not include hydrodynamic
interactions. We will then describe the AHP model, where hydrodynamic is
accounted explicitly, and provide some numerical tests of the implementation.

###  2.1 Active Brownian particles (ABP)

We consider a two-dimensional system with $N_{c}$ disks of mass $m_{c}$ and
diameter $\sigma_{c}$ in a square box size of side $L$. Each disk $i$ has also
an associated axis $\mathbf{n}_{i}=(\cos{\theta_{i}(t)},\sin{\theta_{i}(t)})$,
where $\theta_{i}$ is the angle between the axis and the x-axis and which
evolves over time. $\mathbf{n}_{i}$ represents the direction in which the
self-propulsion occurs.

The particles interact with each other via a short-ranged repulsive potential:

$U(r)=4\epsilon\Big{[}\Big{(}\frac{\sigma}{r}\Big{)}^{64}-\Big{(}\frac{\sigma}{r}\Big{)}^{32}+\frac{1}{4}\Big{]}\Theta(\sigma_{c}-r)$
(1)

where $r$ is the inter-particle distance between the center of masses of each
colloid, $\Theta(r)$ is the Heaviside function ($\Theta(r)=0$ for $r<0$ and
$\Theta(r)=1$ for $r\geq 0$), and $\sigma=2^{-1/32}\sigma_{c}$.

The evolution of the centre of mass of disks is described by a Langevin
equation, with activity modelled as a force $F_{\rm act}$ of constant
magnitude acting along the particle axis $\mathbf{n}_{i}$, while the
propulsion axis changes its direction in time through a diffusion equation:

$\displaystyle
m_{c}\ddot{\bm{r}_{i}}=-\gamma\dot{\bm{r}_{i}}+F_{\text{act}}\mathbf{n}_{i}-\bm{\nabla_{i}}\sum_{j\neq
i}U(r_{ij})+\bm{\xi}_{i},$ (2) $\displaystyle\dot{\theta}_{i}=\eta_{i}\ ,$ (3)

where $i=1,..,N_{c}$, $r_{ij}=|\mathbf{r}_{i}-\mathbf{r}_{j}|$ and $\gamma$ is
the damping coefficient. The terms $\bm{\xi_{i}}$ and $\eta_{i}$ are Gaussian
white noises that mimic the interaction with a thermal bath, with average zero
and variance fixed by the fluctuation-dissipation theorem:

$\displaystyle\langle\xi_{i\alpha}(t)\rangle=0,\;\langle\eta_{i}(t)\rangle=0,$
(4)
$\displaystyle\langle\xi_{i\alpha}(t_{1})\xi_{j\beta}(t_{2})\rangle=2k_{B}T\gamma\delta_{ij}\delta_{\alpha\beta}\delta(t_{1}-t_{2}),$
(5)
$\displaystyle\langle\eta_{i}(t_{1})\eta_{j}(t_{2})\rangle=2D_{\theta}\delta_{ij}\delta(t_{1}-t_{2}),$
(6)

where $\alpha,\beta=1,2$ are the indices of the spatial coordinates, $T$ the
temperature of the system, $k_{B}$ the Boltzmann constant and $D_{\theta}$ the
rotational diffusion coefficient. We express all the quantities in units of
mass, length and energy ($\tilde{m}$, $\tilde{\sigma}$ and $\varepsilon$,
respectively), with the time unit expressed as
$\tau=(\tilde{m}{\tilde{\sigma}}^{2}/\epsilon)^{1/2}$. Note that we fix
$\sigma_{c}=1\tilde{\sigma}$, while $m_{c}$ is varied with respect to the mass
unit $\tilde{m}$. From now on we will drop the units for simplicity.

The density of the system is expressed in terms of the packing fraction
$\phi=\pi{\sigma^{2}_{c}}N_{c}/(4L^{2})$, ratio between the surface occupied
by the colloids and the total system surface $L^{2}$. An important
adimensional number, which measures the ratio between the active work required
to move a particle by $\sigma_{c}$ and the typical thermal energy $k_{B}T$, is
the Péclet number Pe = $F_{\rm act}{\sigma_{c}}/(k_{B}T)$. Another useful
adimensional number is the active Reynolds number, which measures the ratio
between inertial and viscous forces acting on the colloids,
$Re_{act}=\frac{m_{c}F_{a}}{\sigma_{c}\gamma^{2}}$ PhysRevE.90.052130 .

The typical time scales for a single ABP are the inertial time
$t_{I}=m_{c}/\gamma$ and the persistence time $t_{p}=1/D_{\theta}$, with the
latter signaling the crossover to the final diffusive regime and that depends
only on the rate of rotational diffusion and not on the activity parameter. We
can define a useful adimensional number as the ratio between $t_{I}$ and
$t_{p}$, to which we refer to hereafter as the persistence number
$pn=t_{I}/t_{p}$.

We fix in our numerical simulations $\gamma=10$, as previously done for ABP
nostroDISCHI  where the choice $m_{c}=1$ was adopted, which corresponds to
limit inertial effects at small times $t_{I}=0.1$. In the following we keep
fixed $\gamma$ and vary the disk mass $m_{c}$ to consider different inertial
contributions, $k_{B}T=0.05$ and
$D_{\theta}=3k_{B}T/(\sigma^{2}_{c}\gamma)=0.015$. We fix $N_{c}=16384$ and
vary $L$ in order to obtain the correct packing fraction $\phi$. We use
LAMMPSplim1995  to integrate numerically the equations of motion, using a
timestep $\Delta t_{c}=0.001$ and periodic boundary conditions. We fix the
Péclet number to $\textit{Pe}=$ $5$, $10$ and $20$, and vary the packing
fraction $\phi$ between $0.60$ and $0.88$. Within the range of chosen
parameters, $Re_{act}$ is always smaller than one. For each set of parameters
a single realization was considered which was run between $10^{4}$ and
$10^{5}$ simulation time units after steady state was reached. In this time
frame averaged quantities were computed.

###  2.2 Active hydrodynamics particles (AHP)

The ABP model described beforehand does not account explicitly for the
solvent. In order to add this effect, we choose as model a mesoscopic method
known as multi-particle collision (MPC) dynamics, first introduced in male1999
. After briefly describing the MPC model, we will introduce two possible ways
to couple solvent and disks, their dynamics and the specific parameters used
for simulations. Tests of this implementation are presented in section 2.3.

####  2.2.1 Solvent dynamics

The solvent consists of $N_{s}$ identical point-like particles of mass $m_{s}$
embedded in a two-dimensional square box of size $L$. Each particle $i$ is
characterized by a position ${\mathbf{r}}_{i}$ and velocity
${\mathbf{v}}_{i}$, both of which are continuous variables. In this algorithm,
the time is discretized in units $\Delta t_{s}$, and the evolution of the
system is composed by two steps, propagation and collision, which are applied
consecutively for each $\Delta t_{s}$.

In the propagation step, particles are freely streamed according to their
velocities as

${\mathbf{r}}_{i}(t+\Delta
t_{s})={\mathbf{r}}_{i}(t)+{\mathbf{v}}_{i}(t)\Delta t_{s}.$ (7)

In order to perform the collision step, the system is partitioned into cells
of a square lattice with mesh size $\sigma_{s}$. Each cell is the scattering
area where a MPC occurs, which updates particles velocities according to the
rule male1999 ; male2000

${\mathbf{v}}_{i}(t+\Delta
t_{s})={\mathbf{u}}(t)+\Omega[{\mathbf{v}}_{i}(t)-{\mathbf{u}}(t)],$ (8)

where ${\mathbf{u}}=(\sum_{i=1}^{m}{\mathbf{v}}_{i})/m$ is the mean velocity
of the $m$ colliding particles in the cell, also assumed to be the macroscopic
velocity of the fluid. $\Omega$ is a rotation matrix with angle $\pm\alpha$
($0<\alpha<\pi$). The angle $\alpha$ is fixed at the beginning of the
simulation while its sign is assigned with equal probability to every cell at
each time step. In each cell all the $m$ relative velocities are rotated with
the same angle. Linear momentum and kinetic energy are conserved under this
dynamics.

The transport coefficients of this model can be analytically derived. In
particular, for our purposes the kinematic viscosity $\nu_{s}$ and the self-
diffusion coefficient $D_{s}$ will be useful. In 2D the viscosity is equal
tokiku2003 ; ihle2004 :

$\begin{split}\nu_{s}&=\frac{\sigma_{s}^{2}}{2\Delta
t_{s}}\Big{[}\left(\frac{\lambda}{\sigma_{s}}\right)^{2}\Big{(}\frac{n_{s}}{(n_{s}-1+\exp{(-n_{s})})\sin^{2}(\alpha)}-1\Big{)}\\\
&+\frac{(n_{s}-1+\exp{(-n_{s})})(1-\cos(\alpha))}{6n_{s}}\Big{]},\end{split}$
(9)

while the coefficient $D_{s}$ is tuze2006 :

$D_{s}=\frac{\lambda^{2}}{2\Delta
t_{s}}\Big{(}\frac{2n_{s}}{(n_{s}-1+\exp{(-n_{s})})(1-\cos(\alpha))}\Big{)},$
(10)

where $n_{s}=N_{s}\sigma_{s}^{2}/L^{2}$ is the average number of particles per
cell and $\lambda=\Delta t_{s}\sqrt{k_{B}T/m_{s}}$ is the mean-free path.

####  2.2.2 Solvent-colloids coupling

The next step would be to integrate the solvent particles with the colloids,
which means that we need to decide how to couple colloids and solvent
dynamics. Different strategies are possible and a review for MPC with passive
colloids can be found in gommper_review_2009 ; here we adopted the one
implemented in the LAMMPS softwareplipton_colloids_lammps .

In this implementation colloids are evolved for $n$ timesteps $\Delta t_{c}$,
following the equation of motion (2) without the force terms $\bm{\xi}_{i}$
and $\gamma\dot{\bm{r}_{i}}$, which accounted implicitly for the thermal bath
in the ABP model, and are substituted here by the MPC bath. Afterwards solvent
particles are propagated for a timestep equal to $\Delta t_{s}=n\Delta t_{c}$.
Note that both $\Delta t_{c}$ and $\Delta t_{s}$ are expressed in the same
time unit as in the ABP model. Before computing the collision (8), the
algorithm checks if solvent particles are overlapping with disks having
diameter $\sigma_{c}$ and mass $m_{c}$, that is if the position of point-like
solvent particles is inside the disks area. In this case, an exchange of
momentum occurs, followed by a change in the position of solvent particles to
place them out of the colloids, and, finally, the collision step for solvent
particles is applied.

The exchange of momentum is decided by the proper colloid–solvent boundary
condition (BC) adopted, which can be either no-slip or slip. No-slip BC means
that both linear and angular momentum are exchanged between colloid and
solvent particlesbocq1994 , while for slip BC only linear momentum is
transferred as in the case with radial interactions male2000 . Several
implementations of the BC are available, such as the so called thermal BC
inou2002  and the bounce-back collision rulelamu2001 . The latter, used in the
case of no-slip BC, requires the use of phantom particles inside the colloid
while the former does not. Here we choose the thermal BC method, described
below for the slip and no-slip cases, as it is in general useful under forced
flow conditions, like the case of active particles, and is particularly suited
when the solvent mean free path is much smaller than the disk radius padd2005
; padd2006 .

In the no-slip thermal BC, when a solvent particle of velocity ${\mathbf{v}}$
overlaps with a disk, it is moved back to the disk surface along the shortest
vector ${\mathbf{r}}_{d}$ and then streamed for a distance
${\mathbf{v}}^{\prime}\Delta t_{s}\varepsilon$, where ${\mathbf{v}}^{\prime}$
is the updated velocity and $\varepsilon$ is a uniformly distributed random
number in the interval $[0,1]$ hech2005 . The new velocity
${\mathbf{v}}^{\prime}$ is divided in the normal $v_{N}$ and tangential
$v_{T}$ velocity components with respect to the particle-colloid distance, and
chosen according to the stochastic distributions

$\displaystyle p_{N}(v_{N})$ $\displaystyle=$
$\displaystyle(m_{s}v_{N}/k_{B}T)\exp{(-m_{s}v_{N}^{2}/2k_{B}T)},v_{N}>0$ (11)
$\displaystyle p_{T}(v_{T})$ $\displaystyle=$ $\displaystyle\sqrt{m_{s}/2\pi
k_{B}T}\exp{(-m_{s}v_{T}^{2}/2k_{B}T)},$ (12)

centred around the local velocity ${\mathbf{v}}_{d}$ of the colloid surface,
where
${\mathbf{v}}_{d}={\mathbf{V}}+{\mathbf{\omega}}\times({\mathbf{r}}_{d}-{\mathbf{R}})$,
with ${\mathbf{R}}$ being the position of the colloid centre, ${\mathbf{V}}$
and ${\mathbf{\omega}}$ the linear and angular velocities of the colloid.
Regarding the change in momentum for the colloid after the collision, all the
linear and angular momenta variations of the overlapping solvent particles are
summed up as
$\Delta{\mathbf{P}}=\sum_{s}m_{s}({\mathbf{v}}-{\mathbf{v}}^{\prime})$ and
$\Delta{\mathbf{L}}=\sum_{s}m_{s}({\mathbf{r}}_{d}-{\mathbf{R}})\times({\mathbf{v}}-{\mathbf{v}}^{\prime})$,
and the linear and angular velocities of the colloid are updated as:
${\mathbf{V}}^{\prime}={\mathbf{V}}+\Delta{\mathbf{P}}/m_{c}$ and
${\mathbf{\omega}}^{\prime}={\mathbf{\omega}}+\Delta{\mathbf{L}}/I$ where
$I=m_{c}\sigma_{c}^{2}/8$ is the moment of inertia of a disk. In case of high
packing fraction of colloids, it may happen that a single solvent particle can
scatter with several disks in the same timestep $\Delta t_{s}$. Ignoring such
multiple collisions would cause an attractive depletion-like force between
disks hech2005 . This effect can be kept under control allowing a maximum
number $N_{M}$ of multiple collisions. It was found empirically that
$N_{M}\simeq 10$ is the best choice to optimize computational speed and
accuracy.

Figure 1:  VACF for different values of the number of solvent particles per
cell $n_{s}$. Panels (a-c) show the VACF for three different values of
$n_{s}$, namely $n_{s}=20$ (blue curves) $n_{s}=10$ (purple curves), and
$n_{s}=5$ (orange curves). For short times (a), the autocorrelation function
shows a clear exponential decay, which overlaps well with the theoretical
predictions of the the Enskog time, $t_{E}$, shown as a dashed line for each
case. At late times (b) simulations show a long time tail $t^{-1}$ (grey
dashed lines in panels (b) and (c), and dotted purple line in panel (a)). All
the data collapse to the same curve if time is rescaled by $t_{\nu}$ (c).

Figure 2:  VACF for different values of the temperature $k_{B}T$. Panels (a-c)
show the VACF for the values $k_{B}T=0.05$ (red curves) and $k_{B}T=0.1$
(yellow curves), for the same number of solvent particles per cell $n_{s}=10$.
For short times (a), the autocorrelation function shows a clear exponential
decay, which overlaps well with the theoretical prediction of the Enskog time
$t_{E}$ shown as dashed line for each case. At late times (b) simulations show
a long time tail $t^{-1}$ (dotted line in panel (b)). All the data collapse to
the same curve if time is rescaled by $t_{\nu}$ (c). (d) Time evolution, in
semi-logarithmic scale, of the diffusion coefficient computed from the
integral of the VACF, for the same parameters of the yellow curves of panel
(a). The dotted line has the slope $\frac{k_{B}T}{8\pi\rho_{s}\nu_{s}}$.

In the case of slip thermal BC, the tangential component of the fluid particle
velocity is preserved during the scattering with disks; thus no torque is
imparted to colloids. The normal component $v_{N}$ of the solvent particle new
velocity ${\mathbf{v}}^{\prime}$ is sampled from a Gaussian distribution
according to the distribution of Eq. (11) which is centered around the disk
velocity ${\mathbf{V}}$ (the angular velocity is irrelevant since collisions
are now treated as central) plipton_colloids_lammps .

The choice between no-slip and slip BC is directly connected to the way the
axis of colloids ${\mathbf{n}}_{i}$ is evolved. In the first case, the
solvent-disk interaction determines directly through torque exchange how
colloids diffuse rotationally. In the second case, the rotational diffusion is
accounted independently using Eq.(3). In this paper we choose the slip thermal
BC for two reasons. The first one is that in this way we can choose the value
of $D_{\theta}$ independently and match it with the one used in the ABP model.
The second reason is that the integration of slip conditions is much faster
than no-slip ones, since there is no need of considering the integration of
disks angular velocities.

Since we will mostly deal with non-equilibrium simulations, solvent particles
must be coupled to a thermostat to maintain constant temperature. We use the
method of locally rescaling fluid particles velocities $\mathbf{v}_{i}$
relative to the centre of mass velocity $\mathbf{u}$ for each cell by a proper
factor that enforces the correct temperature thermostat . We do not expect
that this approach may alter flow profiles since, as later shown, we will
adopt a very small cell size $\sigma_{s}$ compared to variations in flow
patterns and a very large value of $n_{s}$, the average number of solvent
particles per cell. Note that this implementation ensures only local linear
momentum conservation, while angular and total linear momenta are not
conserved, as typically ensured in simulations of swimmers thee2018 .

Figure 3: Hexatic order parameter color map in the ABP model. (a)-(c) Color
maps of the projection of the local hexatic order parameter of each particle,
$\psi_{6,j}$, onto the direction of the system’s global average,
$\bm{\Psi}=1/N\sum_{j}\psi_{6,j}$, at fixed $\textrm{Pe}=10$ and $m_{c}=44$
for $\phi={0.710,0.730,0.760}$ respectively, for a system of size
$L=256\sigma_{c}$.

####  2.2.3 Parameter choice

In the case of the MPC fluid, an additional set of simulation parameters has
to be set – $n_{s}$, $m_{s}$, $\sigma_{s}$, $\alpha$, $\Delta t_{s}$ – which
will be expressed in terms of the colloids units – $\tilde{m}$,
$\tilde{\sigma}$, $\epsilon$. In order to decide the MPC parameter values, a
set of criteria, listed below, has to be satisfied.

The first criterion is that the solvent has to behave as a fluid (we remind
that MPC particles satisfy an ideal gas equation of state); for such purpose
we need to have a Schmidt number $Sc\simeq 10^{2}-10^{3}$, typical of liquids
ripo2004 . The Schmidt number represents in fact the ratio between the rate of
momentum diffusion and the rate of mass transfer, and for large values of $Sc$
the dynamics resembles the one of a liquid ripo2005 . $Sc$ is defined as
$Sc=\nu_{s}/D_{s}$. Values $Sc\sim O(10)$ can be obtained by requiring small
values of $\lambda$ and large rotation angle ripo2005 . Note that the choice
$\lambda<\sigma_{s}$ is known to break the Galilean invariance ihle2001 ,
although this problem is cured by implementing the random shift procedure
ihle2001  which is here implemented. By using the expressions of $\nu_{s}$ and
$D_{s}$ in the limit of $\lambda/\sigma_{s}\ll 1$, we find that the Schmidt
number depends only on the mean-free path and takes the simple formripo2004 :

$Sc\simeq\frac{1}{12(\lambda/\sigma_{s})^{2}}\ ,$ (13)

where the dependence on $n_{s}$ and $\alpha$ has been omitted since the
dominant contribution is with $\lambda$.

The second criterion is that we want to have the same value of the friction
$\gamma$ as in ABP simulations, where $\gamma$ has the same role as in the
Langevin equation. For the MPC dynamics, this formula is:

$\gamma=C_{2D}\pi\nu_{s}\rho_{s}(\sigma_{c}/2),$ (14)

where $\rho_{s}=n_{s}m_{s}/\sigma_{s}^{2}$ is the solvent density. The
coefficient $C_{2D}$ depends on dimensionality merda2d  and the MPC model
consideredtheersSRD . We performed simulations measuring the velocity of a
colloid dragged by a constant force along a direction in 2D and we fitted a
value of $C_{2D}=1.84\pm 0.1$, using six different values of forces and
averaging over ten realizations. It is evident that also the choice of
$\gamma$ depends directly only on $\lambda$, when all the other parameters are
fixed.

Regarding the active force and the rotational diffusion of the colloids axis,
we do not need any change in the parameters chosen for the ABP, as the active
force and the rotational diffusion are the same as the ones described in the
equation of motion of the ABP model (equation (3)). Thus, the Pe number
depends only on the colloids parameters, and is already set.

The last criterion that we need to follow is to have a very low
compressibility in presence of the active force, in order for the fluid to
remain homogeneous during the time evolution. This criterion was discussed in
padd2006 ; theersSRD . The correct parameters to look at are the Mach number
and the Pumping number. The Mach number $Ma$ is given by the ratio between the
average fluid velocity $\textrm{v}_{s}$ due to the external forces (in our
case due to activity) and the sound velocity
$\textrm{v}_{\textrm{sound}}=\sqrt{2k_{B}T/m_{s}}$ inside the fluid:

$Ma=\frac{\mathrm{v_{s}}}{\mathrm{v}_{sound}}.$ (15)

Its value depends directly on flow velocity. In order to reduce
compressibility effects of the MPC fluid it should be $Ma<0.2$ lamu2001_1 ;
guyon . The Pumping number $Pu$, instead, is the ratio between the active
stationary colloid velocity $F_{act}/\gamma$ and the fluid self-diffusion:

$Pu=\frac{\sigma_{c}F_{act}}{6\gamma D_{s}},$ (16)

and should be less than $1$ theersSRD  in order for the fluid-particle
diffusion to be faster than activity-induced advection, thus avoiding strong
density inhomogeneities in the fluid.

Following these criteria, we chose the cell size to be
$\sigma_{s}=0.2\sigma_{c}$. This guarantees that there is a sufficiently large
number of cells covering a colloid hech2005 . We fix $\alpha=\pi/2$,
$m_{s}=0.15$ and $n_{s}=15$ for the fluid. Typically the colloids and solvent
mass density should match in order for the colloids to be buoyant, so we set
$m_{c}=44.15$ such that
$n_{s}m_{s}/\sigma_{s}^{2}=4m_{c}/(\pi\sigma_{c}^{2})$. This choice provides a
good compromise between avoiding compressibility effectsthee2018 , which for
example arises if we choose lower $n_{s}$, and computational cost, which
arises with higher values of $n_{s}$. We use as $\Delta t_{c}=10^{-4}$ and
$\Delta t_{s}=410\Delta t_{c}$. The temperature $T$ for the solvent and the
other parameters relative to the active force and rotational diffusion remain
the same as the one used for ABP. These parameters lead to the required values
of $\gamma=10.04$ ($\nu_{s}=0.061$ and $\rho_{s}=56.24$), $Sc=99.48$, $Ma=0.1$
and $Pu=0.9$ for the highest $Pe=20$ value considered. We note that the
Reynolds number of the fluid is given by

$Re=u_{0}\sigma_{c}/\nu_{s}=F_{act}\sigma_{c}/\gamma\nu_{s}=6PuD_{s}/\nu_{s}=6Pu/Sc,$
(17)

which is always much less than one for our choice of the parameters. Thus we
are in the low Reynolds number regime.

We start from a close-packed initial configuration of particles positioned in
a triangular lattice, forming a slab, and with the orientation of the self-
propelled force uniformly distributed. The initial velocities of all particles
(fluid particles and colloids) were extracted from a Gaussian distribution
with zero mean and variance $k_{B}T/m_{s}$ and $k_{B}T/m_{c}$ for solvent
particles and colloids, respectively. Given that all the MPC and MD parameters
are the same we are able to consider the same exact active colloids system,
except for the presence of long range hydrodynamic interactions. We fix the
Péclet number to $\textit{Pe}=10$ and $20$, and vary the packing fraction
$\phi$ between $0.60$ and $0.88$, where the hexatic-liquid transition was
found to be critical for $m_{c}=1$ nostroDISCHI . For each set of parameters a
single realization was considered, run between $10^{4}$ and $10^{5}$
simulation time units after steady state was reached, and averaged quantities
where performed during this time frame. To limit the computational cost for
MPC simulations we always fix the box side to $L=128\sigma_{c}$, unless
otherwise specified.

###  2.3 Validation of slip boundary conditions

We focus here on the behaviour of passive colloids embedded in a solvent to
test the accuracy of the previously described slip boundary conditions with
respect to known results for 2D hydrodynamics. Following Ref. sane2009 , we
measure the velocity auto-correlation function (VACF) and the diffusion
coefficient $D_{c}$ of colloids. The parameters chosen for the simulation are
the same as described in the previous section, except that we also varied
either the average number of solvent particles per cell $n_{s}$ or the
temperature $k_{B}T$, always keeping the Schmidt number $Sc\simeq 100$. We
considered large systems, with $L=900\sigma_{c}$ to reduce periodic boundary
effects.

Figure 4: Hexatic order correlation functions for the ABP model. Hexatic order
correlation functions $g_{6}(r)$ for $m_{c}=44$ at $\textrm{Pe}=10$ (left) and
$\textrm{Pe}=20$ (right) for different global packing fractions given in the
keys. Figure 5:  Effects of colloids mass on the Liquid-Hexatic transition. On
the left panel, orientational correlation functions, $g_{6}(r)$, at fixed $\rm
Pe=10$ and $\phi=0.74$ for different values of the mass of the particles given
in the keys. On the right panel, the liquid-Hexatic critical density,
$\phi_{c}$, at fixed $\rm Pe=10$, as a function of the mass of the colloids
$m_{c}$. The solid line is a fit of the data using the function
$\phi_{c}(m_{c})=a+be^{-m_{c}/c}$, with parameters $a=0.71,\,b=0.05,\,c=10.37$
for $\rm Pe=5$; $a=0.72,\,b=0.08,\,c=12.92$ for $\rm Pe=10$; and
$a=0.76,\,b=0.08,\,c=14.30$ for $\rm Pe=20$. The error bars correspond to the
gap $\Delta\phi$ between the densities scanned in the simulations for each
value of the mass $m_{c}$.

At very short times, when hydrodynamics effects can be neglected, the main
contribution to the overall diffusion comes from the local random collisions
between colloid and solvent particles. The VACF is given by

$C_{u}(t)=<u(t)u(0)>=\frac{k_{B}T}{m_{c}}exp(-t/t_{E})\ \ ,$ (18)

where $u$ is a Cartesian components (either $x$ or $y$) of colloids velocity,
$t_{E}=m_{c}/\xi$ is the Enskog time, that is the typical velocity
decorrelation time, and $\xi$ the Enskog friction coefficient given in two
spatial dimensions by sane2009

$\xi=\frac{3\sqrt{2}}{4}\sigma_{c}n_{s}\pi^{3/2}\left(k_{B}T\frac{m_{c}m_{s}}{m_{c}+m_{s}}\right)^{1/2}\
\ .$ (19)

The integral of the VACF is related to the diffusion coefficient $D_{c}$
through the Green-Kubo relation,

$D_{c}=\int_{0}^{\infty}<u(t)u(0)>dt=\frac{k_{B}T}{\xi}\ \ .$ (20)

However, as well known alder1970 ; Ernst1970  fluid dynamic interactions have
an important effect on the long-time behaviour of the VACF. Indeed, due to
momentum conservation, the asymptotic form of the VACF shows an algebraic
decay of the form

$C_{u}(t)=\left(\frac{1}{2\rho_{s}}\right)\frac{k_{B}T}{[4\pi(D_{c}+\nu_{s})t]}\
\ ,$ (21)

for slip boundary conditions in two dimensions. The VACF has a $t^{-1}$ tail,
meaning that the diffusion coefficient $D_{c}$ diverges logarithmically with
time. The long time tail can be expected to appear on the kinematic time scale
$t_{\nu}=\sigma_{c}^{2}/\nu_{s}$, that is the time required by the kinematic
viscosity $\nu_{s}$ to diffuse over the colloid radius. We validate the slip
coupling method introduced in Sec. 2.2.2 between solvent and passive colloids
by testing these predictions.

Since the kinematic viscosity (9) depends only very weakly on $n_{s}$, for
large values of $n_{s}$, and given that

$C_{u}(0)=\frac{k_{B}T}{m_{c}}\\\ ,$ (22)

from equipartition, the long-time tails should all scale onto the same curve
if time is rescaled by $t_{\nu}$.

Figure 1 shows the VACF for three different values of $n_{s}$ and
$k_{B}T=0.1$. As shown in panels (a), for short times, the autocorrelation
function shows clear exponential decay, while at late times (panels (b)-(c))
simulations show a long time tail $t^{-1}$. When plotted as functions of the
reduced time $t/t_{\nu}$, all the data collapse onto the same curve (panel
(c)). The oscillations visible in panels (b) and (c) for long times originate
from sound modes and are a consequence of the finite compressibility of the
MPC fluid combined with the periodic boundary conditions Huang2012 . We
checked that this effect decreases increasing the simulation box size.

The Enskog friction coefficient (19) slightly varies with $n_{s}$; in order to
test the sensibility of the implementation used we fixed $n_{s}=10$

Figure 6: Velocity field induced by an active colloid in the AHP model. Fluid
velocity field around an active colloid for $Pe=20$, in the lab frame (a) and
in the colloid frame (b). The black arrow indicates the direction of the
active force.

and varied the temperature to change the Enskog friction coefficient. Figure
2(a) shows the early time exponential decay of the VACF for the values
$k_{B}T=0.05,0.1$. The measured values of $t_{E}$ are in good agreement with
the theoretical predictions. Also in this case the long-time tail has the
expected $t^{-1}$ slope (panel (b)), and all the curves collapse if time is
rescaled by $t/t_{\nu}$ (panel (c)).

Figure 7: Hexatic order parameter color map in the AHP model (a)-(c) Color
maps of the local hexatic order parameter, $\psi_{6,j}$, as reported in Fig.
3, for $\textrm{Pe}=10$ and $m_{c}=44$, with $\phi={0.78,0.8,0.81}$ from left
to right, for a system of size $L=128\sigma_{c}$.

Using the Green-Kubo relation, and Eq. (21), the diffusion coefficient can be
approximated at long times, assuming that $D_{c}\ll\nu_{s}$ and that the
Enskog and hydrodynamic contributions to the VACF can be separated, as

$\begin{split}D_{c}(t)&=\int_{0}^{t}<u(t)u(0)>dt\\\
&\approx\int_{t_{\nu}}^{t}\frac{k_{B}T}{8\pi\rho_{s}\nu_{s}t}\approx\frac{k_{B}T}{8\pi\rho_{s}\nu_{s}}[\ln
t]_{t_{\nu}}^{t}\ \ .\end{split}$ (23)

Figure 2(d) shows the temporal evolution of the diffusion coefficient computed
from the VACF. On the time scales of the simulation, we observe a behavior
consistent with $D_{c}\simeq ln(t)$, as expected from the $t^{-1}$ tail of the
VACF.

##  3 Hydrodynamic and variable mass effects on hexatic liquid transition

In this Section, we discuss the effects of changing the particles mass for the
2D ABP model and the role of hydrodynamics in the AHP model, using the
numerical framework illustrated in the previous Section. In particular, we
will focus onto characterizing the presence and location of the liquid-hexatic
transition, by varying the system density in a region of the phase diagram at
intermediate active forces where MIPS does not occur for over-damped ABP. The
latter undergo the transition at $\phi_{c}=0.795$ for $\textit{Pe}=10$ and at
$\phi_{c}=0.83$ for $\textit{Pe}=20$ digregorio2017 .

The transition can be characterized by measuring the hexatic order parameter,
$\psi_{6}(r_{i})=\frac{1}{N_{i}}\sum_{j=1}^{N_{i}}e^{i6\theta_{ij}}$, with
$N_{i}$ the number of nearest Voronoi neighbours for particle $i$, and
$\theta_{ij}$ the angle formed between the segment connecting particles $i$
and $j$ and the $x$-axis. From $\psi_{6}(r_{i})$ we can compute the hexatic
correlation function, defined as:

$g_{6}(r)=\frac{\langle\psi_{6}(\bm{r}_{i})\psi_{6}(\bm{r}_{j})\rangle}{\langle\psi_{6}^{2}(\bm{r}_{i})\rangle}\
,$ (24)

where $r=|\bm{r}_{i}-\bm{r}_{j}|$. The transition between hexatic and liquid
phases can be observed by the change in the functional dependence of
$g_{6}(r)$ from exponential decay for short-range order, $g_{6}(r)\sim
e^{-r/l_{c}}$, where $l_{c}$ is the correlation length, to algebraic for
quasi-long-range order, $g_{6}(r)\sim r^{-\beta}$.We use henceforth this
criteria to distinguish between the liquid and the hexatic phase in our
system. In Sec. 3.2.1 we also discuss from a dynamical perspective how
macroscopic flow properties emerge when hydrodynamics is considered.

###  3.1 Effects of different colloids mass in the ABP model

Here we characterize the evolution of ABP following the model described in Sec
2.1 at $\rm Pe=5,10,20$ and compare the results with the ones obtained in Ref.
digregorio2017 . In particular, while in Ref. digregorio2017  only the value
$m_{c}=1$ was considered, here we will study the system with various masses
ranging from $m_{c}=5\text{ to }50$. Thus, the main difference is that here we
are increasing the inertial time $t_{I}=m_{c}/\gamma$, ranging from
$t_{I}=0.5\text{ to }5$ while maintaining the persistence time
$t_{p}=1/D_{\theta}\approx 67$ filyABPs  constant, so that $7\times
10^{-3}<pn<7\times 10^{-2}$. The use of large masses will allow a direct
comparison with the AHP model (where $m_{c}=44$) that will be used in the
following.

We will focus on measuring approximately the value of the critical density
$\phi_{c}$ where the liquid-hexatic transition occurs, computing the hexatic
correlation at a fixed Pe within intervals of $\phi$ ranging from 0.05 to 0.1.

Fig. 3(a)-(c) show typical configurations at $\rm Pe=10$, $m_{c}=44$ and three
different densities. Configurations are colored according to the local hexatic
parameter $\psi_{6,j}$, projected onto its average value. In panel (a)
($\phi=0.71$) we do not observe the appearance of any macroscopic hexatic
domain, while in panel c ($\phi=0.76$) we observe a fully hexatically ordered
system. Panel (b), with $\phi=0.73$, is an intermediate density where
macroscopic and orientationally ordered domains emerge, suggesting that this
density is close to the transition point.

In order to locate the liquid-hexatic transition point at a fixed activity, we
resort to study the hexatic correlation functions, finding the density at
which these functions change from exponential to algebraic decay. Fig.4 shows
these functions for $m_{c}=44$ and $\rm Pe=10,20$. At densities below
$\phi=0.72$ for $\rm Pe=10$ and $\phi=0.74$ for $\rm Pe=20$, we find that the
correlations have an exponential decay, while for larger values the behaviours
that best fits the decay is that of an algebraic function. Thus, we find that
at both activity considered the values where the liquid-hexatic transition
occurs are lowered with respect to the ones at $m_{c}=1$ reported in Ref.
nostroDIFETTI , suggesting that the increase in mass enhances the
orientational ordering at fixed activity. In particular, we estimate
$\phi_{c}=0.730\pm 0.01$ and $0.760\pm 0.01$ for $\rm Pe=10,20$, respectively.

We also checked that this ordering effect occurs while fixing the system
density and activity, and increasing the colloids mass. The correlation
functions in Fig. 5, left side, at $\rm Pe=10$ and $\phi=0.74$, show that by
increasing the mass the system crosses from a liquid state to a hexatic one.
We summarize these measurements in the right panel of Fig. 5, where we show
the location of the critical density for different Pe and different $m_{c}$.
It is evident that the critical density of the liquid-hexatic transition
continuously decreases increasing the value of the mass for the different Pe
considered. Interestingly, the data fit with the function
$\phi_{c}(m_{c})=a+be^{-m_{c}/c}$, with coefficients reported in the caption.

To summarize, the results showed here point out that an enlarged mass, and
therefore an increase in the inertial time $t_{I}$, has an effect of enhancing
the orientational ordering of the system. It is important to note that this is
a non-equilibrium effect not present in the passive system. Indeed, we checked
(not shown) that in the absence of activity the transition density value is
independent of the mass value.We also observed that the asymptotic values for
large $m_{c}$ (coefficient $a$ in the fitting function) are close to the
transition density at $\rm Pe=0$digregorio2017 . When the persistence number
is $pn\gtrsim 10^{-2}$ ($m_{c}\approx 10$), the system behaves closer to the
passive case. On the other hand, when $pn\lesssim 10^{-2}$ the active force
has a disordering effect in the hexatic ordering.

Figure 8: Hexatic order correlation function for the AHP model.  (a)-(b)
Hexatic order correlation function $g_{6}(r)$ at $\textrm{Pe}=10$ (d) and
$\textrm{Pe}=20$ (e) for different global packing fractions given in the keys
and $m_{c}=44$.

###  3.2 Hydrodynamics effects

We now turn our attention to the role of hydrodynamics by studying the AHP
model. To do so, we employ the hybrid mesoscopic approach presented and tested
in section 2, where the MPC solvent is coupled with the active colloids to
account for hydrodynamic interactions. It is important to stress that in our
numerical model no tangential flow velocity is imposed to colloids (they are
not squirmers), thus the resulting velocity field is the result of collisions
between moving colloids and fluid particles. In Fig. 6 we show the velocity
field of our active colloid immersed in a fluid. The flow field strongly
resembles that of a neutral swimmer.

The parameters of AHP, chosen in order to fulfil the constraints discussed in
Sec. 2.2.3, fix the colloid mass to $m_{c}=44$ and $\gamma=10$. In this way,
the AHP simulation results can be directly compared with the ones of ABP with
the same $m_{c}$. We will scan values of $\phi$ between $0.5$ and $0.85$.

####  3.2.1 Liquid-hexatic transition

We start by looking at how the ordering properties are affected by
hydrodynamics. Fig.7(a)-(c) show, for three different densities at
$\textrm{Pe}=10$, the color map of the local hexatic parameter $\psi_{6,j}$
projected onto its average value. In panel (a) ($\phi=0.78$) we do not observe
the appearance of macroscopic hexatic domains, but locally we still observe
small orientationally ordered regions. These regions appear to become larger
upon increasing the density (panel (b), $\phi=0.8$), although global ordering
is not observed. At $\phi=0.81$, panel (c), a single fully hexatically ordered
system is observed. Thus, also AHP present a transition between liquid and
hexatic phases.

Fig. 8 shows the hexatic correlation functions varying the density for
$\textrm{Pe}=10,20$, to be compared with the results presented in Fig. 4 for
the ABP system. For both values of activity, we find that the hexatic order
correlation function shift from an exponential decay to a power-law decay at
substantially higher values of packing fraction $\phi$. More precisely the
transition is located at $\phi_{c}\approx 0.805\pm 0.01$ for $Pe=10$, and
$\phi_{c}\approx 0.840\pm 0.01$ for $Pe=20$.

The increase in value of the transition density $\phi_{c}$ with respect to the
ABP model suggests that the addition of hydrodynamic interactions has a
disordering net effect regarding the global orientational order. This is
opposite to the effect of increasing the particles mass, which instead
promotes hexatic ordering. Indeed, if we measure the average global hexatic
parameter $\psi_{6}=\frac{1}{N}|\sum_{i}^{N}\psi_{6,i}|$ as a function of the
global packing fraction $\phi$ (Fig. 9 ), which increases from 0 to 1 as the
liquid-hexatic transition is crossed we find that the transition is
significantly shifted. Note that both curves converge to almost the same
values for very high densities, suggesting that for densely packed systems
hydrodynamic does not disrupt the ordering properties of colloids.

We also checked (not shown) that in the absence of activity, the transition
density values that limit the coexistence region of the liquid-hexatic
transition of passive colloids nostroDISCHI  are not affected by the presence
of hydrodynamic interactions. However, hydrodynamics produces other relevant
effects which will be now discussed.

Figure 9: Liquid-Hexatic transition in the ABP and AHP model. Global hexatic
parameter as a function of the global packing fraction for $\textrm{Pe}=10$
and $m_{c}=44$,. The orange and blue curves correspond to simulations with and
without hydrodynamics, respectively, for active colloids with the same mass.
Figure 10: Self-sustained active flow. (a)-(c) Color maps of the local hexatic
order parameter for $\textrm{Pe}=10$ and $\phi=0.60,0.80,0.86$, respectively.
Panels (d)-(f) show the corresponding steady state fluid velocity field. The
color code is the same as the one in Fig. 3.

####  3.2.2 Self sustained motion at high density

We now want to better understand the behavior and the role of the fluid
velocity field, which for AHP can be locally organized, while the ABP model
has no such feature, and rely only on hard-core repulsion. Thus, we will have
a deeper look into the velocity field of AHP, and if it can trigger a coherent
motion of small clusters of particles.

Fig. 10 shows the coarse-grained steady state velocity fields of the fluid,
$\mathbf{v(r)}$ (panels d-f) along with associated snapshots of the
configurations colored according to the hexatic parameter (panels a-c), for
AHP with $\textrm{Pe}=10$ and three different values of packing fraction
$\phi$. Coarse-grained velocity fields of the fluid are realized by averaging
the velocity of fluid particles inside blocks of size
$4\sigma_{c}=20\sigma_{s}$ (such large coarse-graining cells are chosen for
the sake of visualization; similar profiles can be obtained with smaller
cells). The first density, $\phi=0.60$ (panels (a), (d)), is characterized by
the absence of any orientational order. At the same time, however, its
corresponding velocity field presents the formation of vortexes along with
regions where flow is both not correlated and lower in magnitude. The
associated velocity field for the active colloids (not shown) has a matching
profile, while the local average direction of the active force is random, and
thus not coherent with the velocity field.

Fig.10b,e show instead a larger density $\phi=0.80$. We observe, here, a case
close to the hexatic transition point, with locally formed fluctuating hexatic
domains with their typical size remaining stationary over time. Along with
these clusters, the flow becomes more coherent than at $\phi=0.60$, with fluid
and colloids having again a similar velocity field. Again, we do not observe a
local average direction of the active force coherent with the flow field. The
same behaviour becomes even more pronounced upon increasing the density
($\phi=0.860$ panels (c) and (f)), where the system is fully orientationally
ordered. In this case, the associated flow field becomes an unidirected self-
sustained flow, with particles moving typically on the same direction, and
with the global direction of the flow slowly changing over time.
Interestingly, this behaviour is similar to travelling bands occurring in
Vicsek-like models vicsek ; levis2018 , where an additional alignment
interaction of active force directions is introduced, which allows particles
to move coherently. Velocity correlations between particles have also been
found in systems of ABP with different persistence times capriniPRL ;
caprini2020 ; Bettolo_Marconi_2021 ; capriniINERTIASM , flowing crystals made
of spontaneously aligning self-propelled hard disks douchot  and self-
sustained spontaneous flows in active gels review2019 ; softmatter ;
giordano2021 ; negro2020 ; NEGRO_2019 .

We do not have at the moment a full theoretical understanding of the emergence
of the coherent motion, which occurs even when there is no orientational
ordering. We can only try to interpret the phenomenology in the following way:
the self-propulsion force of colloids continuously injects energy into the
fluid, setting it into motion. Fluid particles can later self-organize their
motion in a coherent form, and drag colloids along their direction of motion,
which is not necessarily the same direction of the active force of each
particle.

A quantitative measure of this transition to unidirected self-sustained flow,
as a synergetic effect of self-propulsion and hydrodynamic interactions, can
be obtained by measuring the spatial velocity correlation function for the
fluid velocity:

$C_{\textrm{v}}(\mathbf{r})=\frac{\langle\mathbf{v(r)v(0)}\rangle}{\langle\mathbf{v(0)}^{2}\rangle}\
\ .$ (25)

Fig. 11 shows $C_{\textrm{v}}(\mathbf{r})$ for different values of $\phi$, for
$Pe=10$. For low values of $\phi$ (see e.g. $\phi=0.600$) the curve shows an
exponential decay. This corresponds to the case shown in Fig. 10(d),
characterized by the presence of isolated vortexes. When we increase the
density, we observe that the velocity correlation has a slower decay, or a
longer correlation length. Above density $\phi\simeq 0.730$, the correlation
becomes almost constant. We note that the transition in the velocity
correlations between exponential and algebraic decay does not manifest itself
at the liquid-hexatic transition, since the latter appears at higher values of
$\phi$. In the inset of Fig. 11 the velocity correlation function for ABP in
the hexatic phase ($Pe=10$ and $\phi=0.760$) is also shown for comparison. It
shortly decays to zero, while for AHP, even in the liquid case (yellow curve),
the decay is much slower. To gain more insights on the effects of hydrodynamic
interactions we report the radial distribution function $g(r)$ in Fig. 12, for
$\phi=0.750$ both for ABP and AHP at $Pe=10$. We observe that the presence of
fluid in AHP does not considerably change the position of the peaks in the
radial distribution function. However, we notice that the intensity of the
peaks is enhanced in the ABP case, meaning that the fluid interferes with
ordering, thus shifting the hexatic transition to higher densities.

As a last check, we switched off/on hydrodynamics by just removing/adding the
solvent particles and adding/removing the Langevin friction and noise terms in
the colloids equation of motion (3). This enables us to check if a stationary
AHP configuration is naturally able to relax to a stationary conformation of
the ABP when hydrodynamics is switched off. We choose $Pe=10$ and
$\phi=0.795$, a density where the system is hexatically disordered/ordered
with/without hydrodynamics. The results are shown in Fig. 13. We start with
AHP in a fully ordered configuration; after an equilibration time of $10^{4}$
simulation time units, the system forms fluctuating ordered domains which
change over time but do not grow in size (panel (a)). We then turn off
hydrodynamics, and the system gradually sets after $t=10^{5}$ simulation time
units to a almost fully hexatically ordered conformation (panel (b)). The
corresponding colloids velocity field is shown in panels (c)-(d). Note that
the configuration is still not fully ordered only due to the large time
required to relax to the fully ordered state; however we observe that the
global hexatic parameter is steadily growing over time. Switching on
hydrodynamics again the system returns to the configuration shown in Fig.
13(a).

Figure 11: Spatial velocity correlations. Spatial velocity correlation
functions $C_{\textrm{v}}(r)$, for different values of $\phi$, for
$\textrm{Pe}=10$. In the inset the velocity correlation function of ABP at
$Pe=10$ and $\phi=0.760$ is compared with $\phi=0.60$ for AHP at the same
$Pe$.  Figure 12: Radial distribution function. Radial distribution function
for ABP and AHP models with $\phi=0.750$ at $Pe=10$.

##  4 Conclusions

We have studied with extensive simulations the role of particles mass and
hydrodynamics in active colloids, and showed how they affect the liquid-
hexatic transition in an intermediate activity regime in which MIPS does not
occur yet ( $\textit{Pe}=10,20$).

We have first characterized the ABP by changing their mass, while maintaining
the same Pe and $D_{\theta}$, so that we have a non-trivial interplay between
the inertial time and the persistence time $t_{p}=1/D_{\theta}$. We showed
that the critical density of the transition is shifted to a lower density upon
increasing the colloid mass. This critical density is close to the one found
at $\textit{Pe}=0$, suggesting that inertia has an orientational ordering
effect on the system bringing the system closer to equilibrium behaviour and
counteracting the disordering role of self-propulsion.

When hydrodynamic interactions are taken into account, we found instead that
the liquid-hexatic transition moves towards higher values of packing fraction
$\phi$, thus suggesting that hydrodynamics has a net effect of orientationally
disordering the system. We also analyzed the fluid velocity field of AHP, and
found at $Pe=10$ two results: i) the formation below $\phi\approx 0.72$ of
small correlated velocity field regions, characterized by the presence of
vortices, that are not associated to any local orientational ordering; ii) the
arisal above $\phi\approx 0.720$ of a self sustained motion, with the fluid
particles moving in one direction. This change in behavior has been
characterized by measuring the spatial velocity correlation which was found to
change from an exponential to an algebraic decay.

Regarding the role of inertia, it will be interesting in the future to
reconstruct a phase diagram similar to the one of Ref.nostroDISCHI , by
characterizing in more detail the hexatic phase and the location of the solid
phase. Regarding AHP, instead, it will be necessary to better describe the
physical mechanisms producing the vortices at smaller densities and the
transition to a self-sustained motion at larger densities. It remains an open
question whether such a scenario is still encountered in quasi 2D and 3D
geometries as well as in experiments with wet active colloids. It would also
be of interest to investigate the effect of no-slip boundary conditions, which
would completely determine the colloid angular diffusion and could induce
additional cooperative effects, the effects of changing colloidal mass, and
investigating in more details particle-particle flow interactions and local
velocity field effects.We hope that our results can boost further research in
this direction.

Figure 13: Switch between ABP and AHP model. (a)-(b) Snapshots of the system
at different times $t=10^{4},10^{5}$ of the local hexatic order parameter, at
$\textrm{Pe}=10$ and $\phi=0.795$, before (panel (a)) and after (panel (b))
switching off hydrodynamics. The corresponding colloids velocity fields are
shown in panels (c)-(d).

## Acknowledgements

This work was possible thank to the access to Bari ReCaS e-Infrastructure
funded by MIUR through PON Research and Competitiveness 2007-2013 Call 254
Action I and MARCONI at CINECA (Project INF22-fieldturb) under CINECA-INFN
agreement. We acknowledge funding from MIUR Project No. PRIN 2020/PFCXPE.

## Author contribution statement

All authors contributed equally to conceptualizing the research, analysing the
results and writing the paper. GN and CBC carried out the simulations.

## Data Availability Statement

The datasets generated during and/or analysed during the current study are
available from the corresponding author on reasonable request.

## References

  * [1] Y. Fily and M.C. Marchetti.  Athermal phase separation of self-propelled particles with no alignment.  Phys. Rev. Lett., 108:235702, Jun 2012.
  * [2] P. Romanczuk, M. B’́ar, W. Ebeling, B. Lindner, and L. Schimansky-Geier.  Active brownian particles. from individual to collective stochastic dynamics.  Eur. Phys. J. Spec. Top., 202:1–162, 2012.
  * [3] J. Bialké, T. Speck, and H. Löwen.  Crystallization in a dense suspension of self-propelled particles.  Phys. Rev. Lett., 108:168301, Apr 2012.
  * [4] 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, Nov 2016.
  * [5] W. Bialek, A. Cavagna, I. Giardina, T. Mora, E. Silvestri, M. Viale, and A.M. Walczak.  Statistical mechanics for natural flocks of birds.  Proceedings of the National Academy of Sciences, 109(13):4786–4791, 2012.
  * [6] A.G. Be’er.  A statistical physics view of swarming bacteria.  Mov Ecol, 7, 2019.
  * [7] D. Needleman and Z. Dogic.  Active matter at the interface between materials science and cell biology.  Nat. Rev. Mater., 2, 2017.
  * [8] P. Digregorio, D., 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, Aug 2018.
  * [9] 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, Oct 2020.
  * [10] M.E. Cates and J. Tailleur.  Motility-induced phase separation.  Annual Review of Condensed Matter Physics, 6(1):219–244, 2015.
  * [11] G. Gonnella, D. Marenduzzo, A. Suma, and A. Tiribocchi.  Motility-induced phase separation and coarsening in active matter.  Comptes Rendus Physique, 16(3):316–331, 2015.
  * [12] G.S. Redner, M.F. Hagan, and A.Baskaran.  Structure and dynamics of a phase-separating active colloidal fluid.  Phys. Rev. Lett., 110:055701, Jan 2013.
  * [13] I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger, and T. Speck.  Dynamical clustering and phase separation in suspensions of self-propelled colloidal particles.  Phys. Rev. Lett., 110:238301, Jun 2013.
  * [14] F. Ginot, I. Theurkauff, D. Levis, C. Ybert, L. Bocquet, L. Berthier, and C. Cottin-Bizonne.  Nonequilibrium equation of state in suspensions of active colloids.  Phys. Rev. X, 5:011004, Jan 2015.
  * [15] M. Rubenstein, A. Cornejo, and R. Nagpal.  Programmable self-assembly in a thousand-robot swarm.  Science, 345(6198):795–799, 2014.
  * [16] P. Chvykov, T. A. Berrueta, A. Vardhan, W. Savoie, A. Samland, T. D. Murphey, K. Wiesenfeld, D.I. Goldman, and J.L. England.  Low rattling: A predictive principle for self-organization in active collectives.  Science, 371(6524):90–95, 2021.
  * [17] S. Li, R. Batra, D.Brown, H. Chang, N. Ranganathan, C. Hoberman, D. Rus, and H. Lipson.  Particle robotics based on statistical mechanics of loosely coupled components.  Nature, 567(7748):361–365, 2019.
  * [18] R. Kapral.  Perspective: Nanomotors without moving parts that propel themselves in solution.  The Journal of Chemical Physics, 138(2):020901, 2013.
  * [19] Y. Katz, K. Tunstrøm, C.C. Ioannou, C. Huepe, and I.D. Couzin.  Inferring the structure and dynamics of interactions in schooling fish.  Proceedings of the National Academy of Sciences, 108(46):18720–18725, 2011.
  * [20] A. Celani K. Gustavsson, L. Biferale and S. Colabrese.  Finding efficient swimming strategies in a three-dimensional chaotic flow by reinforcement learning.  Eur. Phys. J. E, 40(110), 2017.
  * [21] A. Doostmohammadi and J.M. Yeomans.  Coherent motion of dense active matter.  Eur. Phys. J. Spec. Top., 227:2401–2411, 2019.
  * [22] E. Tjhung and L. Berthier.  Analogies between growing dense active matter and soft driven glasses.  Phys. Rev. Research, 2:043334, Dec 2020.
  * [23] R. Mandal, P.J. Bhuyan, P. Chaudhuri, C. Dasgupta, and M. Rao.  Extreme active matter at high densities.  Nat Commun, 11, 2020.
  * [24] S. Henkes, K. Kostanjevec, J.M. Collinson, R. Sknepnek, and E. Bertin.  Dense active matter model of motion patterns in confluent cell monolayers.  Nat Commun, 11, 2020.
  * [25] D. Bi, X. Yang, M.C. Marchetti, and M.L. Manning.  Motility-driven glass and jamming transitions in biological tissues.  Phys. Rev. X, 6:021011, Apr 2016.
  * [26] J. Yan, M. Han, J. Zhang, C. Xu, E. Luijten, and S. Granick.  Reconfiguring active particles by electrostatic imbalance.  Nature Mater, 15:1095–1099, 2016.
  * [27] J.U. Klamser, S.C.Kapfer, and W. Krauth.  Thermodynamic phases in two-dimensional active matter.  Nat Commun, 9, 2018.
  * [28] P. Digregorio, D. Levis, L.F. Cugliandolo, G. Gonnella, and I. Pagonabarraga.  Unified analysis of topological defects in 2d systems of active and passive disks.  Soft Matter, 18:566–591, 2022.
  * [29] E.P. Bernard and W. Krauth.  Two-step melting in two dimensions: First-order liquid-hexatic transition.  Phys. Rev. Lett., 107:155704, Oct 2011.
  * [30] S.C. Kapfer and W. Krauth.  Two-dimensional melting: From liquid-hexatic coexistence to continuous transitions.  Phys. Rev. Lett., 114:035702, Jan 2015.
  * [31] A.L. Thorneywork, J.L. Abbott, Dirk G. A. L. D.G.A.L. Aarts, and R.P.A. Dullens.  Two-dimensional melting of colloidal hard spheres.  Phys. Rev. Lett., 118:158001, Apr 2017.
  * [32] J.M.Kosterlitz and D.J.Thouless.  Ordering, metastability and phase transitions in two-dimensional systems.  Journal of Physics C: Solid State Physics, 6(7):1181–1203, apr 1973\.
  * [33] B.I. Halperin and D.R. Nelson.  Theory of two-dimensional melting.  Phys. Rev. Lett., 41:121–124, Jul 1978.
  * [34] A.P. Young.  Melting and the vector coulomb gas in two dimensions.  Phys. Rev. B, 19:1855–1866, Feb 1979.
  * [35] P. Nie, J. Chattoraj, A. Piscitelli, P. Doyle, R. Ni, and M.P. Ciamarra.  Stability phase diagram of active brownian particles.  Phys. Rev. Research, 2:023010, Apr 2020.
  * [36] L. Caprini and U.M.B. Marconi.  Inertial self-propelled particles.  The Journal of Chemical Physics, 154(2):024902, 2021.
  * [37] A.K. Omar, K. Klymko, T. GrandPre, P.L. Geissler, and J.F. Brady.  Tuning nonequilibrium phase transitions with inertia.  arXiv:2108.10278, 2021.
  * [38] I. Petrelli, L.F. Cugliandolo, G. Gonnella, and A. Suma.  Effective temperatures in inhomogeneous passive and active bidimensional brownian particle systems.  Phys. Rev. E, 102:012609, Jul 2020.
  * [39] H. L’́owen.  Inertial effects of self-propelled particles: From active brownian to active langevin motion.  The Journal of Chemical Physics, 152(4):040901, 2020.
  * [40] M. Theers, E. Westphal, K. Qi, R.G. Winkler, and G. Gompper.  Clustering of microswimmers: interplay of shape and hydrodynamics.  Soft Matter, 14:8590–8603, 2018.
  * [41] R. Matas-Navarro, R. Golestanian, T.B. Liverpool, and S.M. Fielding.  Hydrodynamic suppression of phase separation in active suspensions.  Phys. Rev. E, 90:032304, Sep 2014.
  * [42] G. Gompper, R.G. Winkler, T. Speck, A. Solon, C. Nardini, F. Peruani, H. Löwen, R. Golestanian, U.B. Kaupp, L. Alvarez, T. Kiørboe, E. Lauga, W.C.K. Poon, A. DeSimone, S. Muiños-Landin, A. Fischer, N.A. Söker, F. Cichos, R. Kapral, P. Gaspard, M. Ripoll, F. Sagues, A. Doostmohammadi, J.M. Yeomans, I.S. Aranson, C. Bechinger, H. Stark, C.K. Hemelrijk, F.J. Nedelec, T. Sarkar, T. Aryaksama, M. Lacroix, G. Duclos, V. Yashunsky, P. Silberzan, M. Arroyo, and S. Kale.  The 2020 motile active matter roadmap.  Journal of Physics: Condensed Matter, 32(19):193001, feb 2020.
  * [43] A. Zöttl and H. Stark.  Hydrodynamics determines collective motion and phase behavior of active colloids in quasi-two-dimensional confinement.  Phys. Rev. Lett., 112:118101, Mar 2014.
  * [44] M. Theers, Westphal E, K. Qi, R.G. Winkler, and G. Gompper".  Clustering of microswimmers: interplay of shape and hydrodynamics.  Soft Matter, 14:8590–8603, 2018.
  * [45] A. Malevanets and R. Kapral.  Mesoscopic model for solvent dynamics.  The Journal of Chemical Physics, 110(17):8605–8613, 1999.
  * [46] G. Gompper, T. Ihle, D.M. Kroll, and R.G. Winkler.  Multi-Particle Collision Dynamics: A Particle-Based Mesoscale Simulation Approach to the Hydrodynamics of Complex Fluids, pages 1–87.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2009.
  * [47] D.S. Bolintineanu, G.S. Grest, J.B. Lechman, F. Pierce, S.J. Plimpton, and P.R. Schunk.  Particle dynamics modeling methods for colloid suspensions.  Computational Particle Mechanics, 1(3):321–356, 2014.
  * [48] 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, Nov 2014.
  * [49] S. Plimpton.  Fast parallel algorithms for short-range molecular dynamics.  Journal of Computational Physics, 117(1):1–19, 1995.
  * [50] A. Malevanets and R. Kapral.  Solute molecular dynamics in a mesoscale solvent.  The Journal of Chemical Physics, 112(16):7260–7269, 2000.
  * [51] N. Kikuchi, C.M. Pooley, J.F. Ryder, and J.M. Yeomans.  Transport coefficients of a mesoscopic fluid dynamics model.  The Journal of Chemical Physics, 119(12):6388–6395, 2003.
  * [52] T. Ihle, E. Tüzel, and D.M. Kroll.  Resummed green-kubo relations for a fluctuating fluid-particle model.  Phys. Rev. E, 70:035701, Sep 2004.
  * [53] E. Tüzel, T. Ihle, and D.M. Kroll.  Dynamic correlations in stochastic rotation dynamics.  Phys. Rev. E, 74:056702, Nov 2006.
  * [54] L. Bocquet and J. Barrat.  Hydrodynamic boundary conditions, correlation functions, and kubo relations for confined fluids.  Phys. Rev. E, 49:3079–3092, Apr 1994.
  * [55] Y. Inoue, Y. Chen, and H. Ohashi.  Development of a simulation model for solid objects suspended in a fluctuating fluid.  Journal of Statistical Physics, 107(1):85–100, 2002.
  * [56] A. Lamura, G. Gompper, T. Ihle, and D.M. Kroll.  Multi-particle collision dynamics: Flow around a circular and a square cylinder.  Europhysics Letters (EPL), 56(3):319–325, nov 2001.
  * [57] J.T. Padding, A. Wysocki, H. Löwen, and A.A. Louis.  Stick boundary conditions and rotational velocity auto-correlation functions for colloidal particles in a coarse-grained representation of the solvent.  Journal of Physics: Condensed Matter, 17(45):S3393–S3399, oct 2005\.
  * [58] J.T. Padding and A.A. Louis.  Hydrodynamic interactions and brownian forces in colloidal suspensions: Coarse-graining over time and length scales.  Phys. Rev. E, 74:031402, Sep 2006.
  * [59] M. Hecht, J. Harting, T. Ihle, and H.J. Herrmann.  Simulation of claylike colloids.  Phys. Rev. E, 72:011408, Jul 2005.
  * [60] H. Chien-Cheng, A. Varghese, G. Gompper, and R.G. Winkler.  Thermostat for nonequilibrium multiparticle-collision-dynamics simulations.  Phys. Rev. E, 91:013310, Jan 2015.
  * [61] M. Ripoll, K. Mussawisade, R.G. Winkler, and G.Gompper.  Low-reynolds-number hydrodynamics of complex fluids by multi-particle-collision dynamics.  Europhysics Letters (EPL), 68(1):106–112, oct 2004.
  * [62] M. Ripoll, K. Mussawisade, R.G. Winkler, and G. Gompper.  Dynamic regimes of fluids simulated by multiparticle-collision dynamics.  Phys. Rev. E, 72:016701, Jul 2005.
  * [63] O.M. Braun and M. Peyrard.  Friction in a solid lubricant film.  Phys. Rev. E, 63:046110, Mar 2001.
  * [64] B. Liu, J. Goree, and O.S. Vaulina.  Test of the stokes-einstein relation in a two-dimensional yukawa liquid.  Phys. Rev. Lett., 96:015005, Jan 2006.
  * [65] A. Lamura and G. Gompper.  Numerical study of the flow around a cylinder using multi-particle collision dynamics.  Eur. Phys. J. E, 9(5):477–485, 2002.
  * [66] E. Guyon, J.P. Hulin, L. Petit, and C.D. Mitescu.  Physical Hydrodynamics.  Oxford University Press, 2015.
  * [67] J.T. Padding J. Sané and A.A. Louis.  Hydrodynamics of confined colloidal fluids in two dimensions.  Phys. Rev. E, 79:051402, May 2009.
  * [68] B.J. Alder and T.E. Wainwright.  Decay of the velocity autocorrelation function.  Phys. Rev. A, 1:18–21, Jan 1970.
  * [69] M.H. Ernst, E.H. Hauge, and J.M.J. van Leeuwen.  Asymptotic time behavior of correlation functions.  Phys. Rev. Lett., 25:1254–1256, Nov 1970.
  * [70] C. Huang, G. Gompper, and R.G. Winkler.  Hydrodynamic correlations in multiparticle collision dynamics fluids.  Phys. Rev. E, 86:056711, Nov 2012.
  * [71] 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, Dec 2017.
  * [72] V. Tamás and A. Zafeiris.  Collective motion.  Physics Reports, 517(3):71–140, 2012.  Collective motion.
  * [73] A. Martín-Gómez, D. Levis, A. Díaz-Guilera, and I. Pagonabarraga.  Collective motion of active brownian particles with polar alignment.  Soft Matter, 14:2610–2618, 2018.
  * [74] L. Caprini, U.M.B. Marconi, and A. Puglisi.  Spontaneous velocity alignment in motility-induced phase separation.  Phys. Rev. Lett., 124:078001, Feb 2020.
  * [75] L. Caprini, U.M.B. Marconi, C. Maggi, M. Paoluzzi, and A. Puglisi.  Hidden velocity ordering in dense suspensions of self-propelled disks.  Phys. Rev. Research, 2:023321, Jun 2020.
  * [76] U.M.B. Marconi, L. Caprini, and A. Puglisi.  Hydrodynamics of simple active liquids: the emergence of velocity correlations.  New Journal of Physics, 23(10):103024, oct 2021.
  * [77] L. Caprini and U.M.B. Marconi.  Spatial velocity correlations in inertial systems of active brownian particles.  Soft Matter, 17:4109–4121, 2021.
  * [78] G. Briand, M. Schindler, and O. Dauchot.  Spontaneously flowing crystal of self-propelled particles.  Phys. Rev. Lett., 120:208001, May 2018.
  * [79] L.N. Carenza, G. Gonnella, A. Lamura, G. Negro, and A. Tiribocchi.  Lattice boltzmann methods and active fluids.  Eur. Phys. J. E, 42(6):81, 2019.
  * [80] G. Negro, L.N. Carenza, A.Lamura, A. Tiribocchi, and G. Gonnella.  Rheology of active polar emulsions: from linear to unidirectional and inviscid flow, and intermittent viscosity.  Soft Matter, 15:8251–8265, 2019.
  * [81] M.G. Giordano, F. Bonelli, L.N. Carenza, G. Gonnella, and G. Negro.  Activity-induced isotropic-polar transition in active liquid crystals.  EPL (Europhysics Letters), 133(5):58004, mar 2021.
  * [82] 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.
  * [83] G. Negro, A. Lamura, G. Gonnella, and D. Marenduzzo.  Hydrodynamics of contraction-based motility in a compressible active fluid.  EPL (Europhysics Letters), 127(5):58001, oct 2019.