Hydrodynamic effects on the liquid-hexatic transition of active colloids

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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi $$\end{document}ϕ 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.


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 [1][2][3][4] 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 [5], to swarming in bacterial colonies [6] and dynamics in cells' cytoskeleton [7].
Active Brownian Particles (ABP) models are very popular among SPP models [1,8,9]. 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) [1,[10][11][12][13] 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 microswimmers [13,14], opening the perspectives for a sysa e-mail: claudio.caporusso@ba.infn.it (corresponding author) tematic control of active systems and collective motion, with the purpose to exploit some of their unique features for technological uses, for instance in robotics [15][16][17], realisation of biological machines [18], or understanding of flocking intelligence [19,20].
Of particular interest is the characterization of ABP in the dense regime, see e.g. spontaneous flow [21] or glassy behaviour [22,23] in biological tissues, biofilms, cell mono-layers [24,25], and can be considered a target for the development of new materials [26]. In two dimensions (2D), ABP present ordering phase transitions when the density of the system is increased [8,9,27,28], which are connected to those encountered for passive hard colloids [29][30][31]. 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 scenario [32][33][34]. If activity is high enough, instead, MIPS takes place, as a phase separation between a dense phase and a gaseous one [8].
The aforementioned features of the ABP phase diagram have been well established in the context of overdamped motion and without an explicit underlying thermo-hydrodynamic bath. At the same time, there 75 Page 2 of 15 Eur. Phys. J. E (2022) 45 :75 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 [35,36]. It has been pointed out in [37] 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 [38], and to highly inhibit phase segregation [39]. 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 [40][41][42], as hydrodynamics favour reorientation of particles' self propulsion direction, while in quasi-2D systems MIPS has been observed for low-density fluids [43] and not when the fluid was made incompressible [42,44]. For elongated colloids, steric alignment and hydrodynamics show highly non-trivial interplay, such that MIPS is enhanced for pullers and suppressed for pushers [40].
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 [45,46], which seamlessly integrate with the dynamics of active Brownian particles [47]. 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 Sect. 2 we discuss the numerical methods and parameter choice for the ABP model and for the AHP model, with Sect. 2.3 providing several tests for implementation of the latter model. In Sect. 3.1 we discuss how the liquidhexatic scenario changes by varying the active colloids mass, while in Sect. 3.2 we discuss the effects due to hydrodynamics interactions. Finally we draw some conclusions discussing the main findings.

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.

Active Brownian particles (ABP)
We consider a two-dimensional system with N c disks of mass m c and diameter σ c in a square box size of side L. Each disk i has also an associated axis n i = (cos θ i (t), sin θ i (t)), where θ i is the angle between the axis and the x-axis and which evolves over time. n i represents the direction in which the self-propulsion occurs.
The particles interact with each other via a shortranged repulsive potential: where r is the inter-particle distance between the center of masses of each colloid, Θ(r) is the Heaviside function (Θ(r) = 0 for r < 0 and Θ(r) = 1 for r ≥ 0), and σ = 2 −1/32 σ c . The evolution of the centre of mass of disks is described by a Langevin equation, with activity modelled as a force F act of constant magnitude acting along the particle axis n i , while the propulsion axis changes its direction in time through a diffusion equation: where i = 1, .., N c , r ij = |r i − r j | and γ is the damping coefficient. The terms ξ i and η i are Gaussian white noises that mimic the interaction with a thermal bath, with average zero and variance fixed by the fluctuationdissipation theorem: where α, β = 1, 2 are the indices of the spatial coordinates, T the temperature of the system, k B the Boltzmann constant and D θ the rotational diffusion coefficient. We express all the quantities in units of mass, length and energy (m,σ and ε, respectively), with the time unit expressed as τ = (mσ 2 / ) 1/2 . Note that we fix σ c = 1σ, while m c is varied with respect to the mass unitm. From now on we will drop the units for simplicity. The density of the system is expressed in terms of the packing fraction φ = πσ 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 σ c and the typical thermal energy k B T , is the Péclet number Pe = F act σ 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 = mcFa σcγ 2 [48]. The typical time scales for a single ABP are the inertial time t I = m c /γ and the persistence time t p = 1/D θ , 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 γ = 10, as previously done for ABP [8] 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 γ and vary the disk mass m c to consider different inertial contributions, k B T = 0.05 and D θ = 3k B T /(σ 2 c γ) = 0.015. We fix N c = 16384 and vary L in order to obtain the correct packing fraction φ. We use LAMMPS [49] to integrate numerically the equations of motion, using a timestep Δt c = 0.001 and periodic boundary conditions. We fix the Péclet number to Pe = 5, 10 and 20, and vary the packing fraction φ 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.

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 multiparticle collision (MPC) dynamics, first introduced in [45]. 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 Sect. 2.3.

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 position r i and velocity v i , both of which are continuous variables. In this algorithm, the time is discretized in units Δt s , and the evolution of the system is composed by two steps, propagation and collision, which are applied consecutively for each Δt s .
In the propagation step, particles are freely streamed according to their velocities as In order to perform the collision step, the system is partitioned into cells of a square lattice with mesh size σ s . Each cell is the scattering area where a MPC occurs, which updates particles velocities according to the rule [45,50] v where u = ( m i=1 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. Ω is a rotation matrix with angle ±α (0 < α < π). The angle α 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 ν s and the self-diffusion coefficient D s will be useful. In 2D the viscosity is equal to [51,52]: while the coefficient D s is [53]: where n s = N s σ 2 s /L 2 is the average number of particles per cell and λ = Δt s k B T /m s is the mean-free path.

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 [46]; here we adopted the one implemented in the LAMMPS software [47].
In this implementation colloids are evolved for n timesteps Δt c , following the equation of motion (2) without the force terms ξ i and γṙ 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 Δt s = nΔt c . Note that both Δt c and Δ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 σ 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 particles [54], while for slip BC only linear momentum is transferred as in the case with radial interactions [50]. Several implementations of the BC are available, such as the so called thermal BC [55] and the bounce-back collision rule [56]. 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 [57,58].
In the no-slip thermal BC, when a solvent particle of velocity v overlaps with a disk, it is moved back to the disk surface along the shortest vector r d and then streamed for a distance v Δt s ε, where v is the updated velocity and ε is a uniformly distributed random number in the interval [0, 1] [59]. The new velocity v 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 centred around the local velocity v d of the colloid sur- 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 Δt s . Ignoring such multiple collisions would cause an attractive depletion-like force between disks [59]. This effect can be kept under control allowing a maximum number N M of multiple collisions. It was found empirically that N M 10 is the best choice to optimize computational speed and accuracy. 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 v is sampled from a Gaussian distribution according to the distribution of Eq. (11) which is centered around the disk velocity V (the angular velocity is irrelevant since collisions are now treated as central) [47].
The choice between no-slip and slip BC is directly connected to the way the axis of colloids 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 θ 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 v i relative to the centre of mass velocity u for each cell by a proper factor that enforces the correct temperature [60]. We do not expect that this approach may alter flow profiles since, as later shown, we will adopt a very small cell size σ 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 [44].

Parameter choice
In the case of the MPC fluid, an additional set of simulation parameters has to be set -n s , m s , σ s , α, Δt swhich will be expressed in terms of the colloids unitsm,σ, . 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 10 2 − 10 3 , typical of liquids [61]. 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 [62]. Sc is defined (10) can be obtained by requiring small values of λ and large rotation angle [62]. Note that the choice λ < σ s is known to break the Galilean invariance [63], although this problem is cured by implementing the random shift procedure [63] which is here implemented. By using the expressions of ν s and D s in the limit of λ/σ s 1, we find that the Schmidt number depends only on the mean-free path and takes the simple form [61]: where the dependence on n s and α has been omitted since the dominant contribution is with λ.
The second criterion is that we want to have the same value of the friction γ as in ABP simulations, where γ has the same role as in the Langevin equation. For the MPC dynamics, this formula is: where ρ s = n s m s /σ 2 s is the solvent density. The coefficient C 2D depends on dimensionality [64] and the MPC model considered [40]. 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 ± 0.1, using six different values of forces and averaging over ten realizations. It is evident that also the choice of γ depends directly only on λ, 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 [40,58]. 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 v s due to the external forces (in our case due to activity) and the sound velocity v sound = 2k B T /m s inside the fluid: Its value depends directly on flow velocity. In order to reduce compressibility effects of the MPC fluid it should be Ma < 0.2 [65,66]. The Pumping number P u, instead, is the ratio between the active stationary colloid velocity F act /γ and the fluid self-diffusion:  and should be less than 1 [40] in order for the fluidparticle 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 σ s = 0.2σ c . This guarantees that there is a sufficiently large number of cells covering a colloid [59]. We fix α = π/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 /σ 2 s = 4m c /(πσ 2 c ). This choice provides a good compromise between avoiding compressibility effects [44], which for example arises if we choose lower n s , and computational cost, which arises with higher values of n s . We use as Δt c = 10 −4 and Δt s = 410Δ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 γ = 10.04 (ν s = 0.061 and ρ s = 56.24), Sc = 99.48, Ma = 0.1 and P u = 0.9 for the highest Pe = 20 value considered. We note that the Reynolds number of the fluid is given by 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 Pe = 10 and 20, and vary the packing fraction φ between 0.60 and 0.88, where the hexatic- liquid transition was found to be critical for m c = 1 [8]. 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σ c , unless otherwise specified.

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. [67], 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 100. We considered large systems, with L = 900σ c to reduce periodic boundary effects.
At very short times, when hydrodynamics effects can be neglected, the main contribution to the overall dif- where u is a Cartesian components (either x or y) of colloids velocity, t E = m c /ξ is the Enskog time, that is the typical velocity decorrelation time, and ξ the Enskog friction coefficient given in two spatial dimensions by [67] The integral of the VACF is related to the diffusion coefficient D c through the Green-Kubo relation, However, as well known [68,69], 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 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 ν = σ 2 c /ν s , that is the time required by the kinematic viscosity ν s to diffuse over the colloid radius. We validate the slip coupling method introduced in Sect. 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 from equipartition, the long-time tails should all scale onto the same curve if time is rescaled by t ν . Figure 1 shows the VACF for three different values of n s and k B T = 0.1. As shown in panel (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 ν , 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 [70]. 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 and varied the temperature to change the Enskog friction coefficient. Figure 2a 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 ν (panel (c)).
Using the Green-Kubo relation and Eq. (21), the diffusion coefficient can be approximated at long times, assuming that D c ν s and that the Enskog and hydrodynamic contributions to the VACF can be separated, as (23) Figure 2d 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 ln(t), as expected from the t −1 tail of the VACF.

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 φ c = 0.795 for Pe = 10 and at φ c = 0.83 for Pe = 20 [71].
The transition can be characterized by measuring the hexatic order parameter, ψ 6 (r i ) = 1

Ni
Ni j=1 e i6θij , with N i the number of nearest Voronoi neighbours for particle i, and θ ij the angle formed between the segment connecting particles i and j and the x-axis. From ψ 6 (r i ) we can compute the hexatic correlation function, defined as: where r = |r i −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) ∼ e −r/lc , where l c is the correlation length, to algebraic for quasi-long-range order, g 6 (r) ∼ r −β . We use henceforth this criteria to distinguish between the liquid and the hexatic phase in our system. In Sect. 3.2.2 we also discuss from a dynamical perspective how macroscopic flow properties emerge when hydrodynamics is considered.

Effects of different colloids mass in the ABP model
Here we characterize the evolution of ABP following the model described in Sect 2.1 at Pe = 5, 10, 20 and compare the results with the ones obtained in Ref. [71]. In particular, while in Ref. [71] only the value m c = 1 was considered, here we will study the system with various masses ranging from m c = 5 to 50. Thus, the main difference is that here we are increasing the inertial time t I = m c /γ, ranging from t I = 0.5 to 5 while maintaining the persistence time t p = 1/D θ ≈ 67 [1] constant, so that 7 × 10 −3 < pn < 7 × 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 φ c where the liquid-hexatic transition occurs, computing the hexatic correlation at a fixed Pe within intervals of φ ranging from 0.05 to 0.1. 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. Figure 4 shows these functions for m c = 44 and Pe = 10, 20. At densities below φ = 0.72 for Pe = 10 and φ = 0.74 for 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 activities considered the values where the liquid-hexatic transition occurs are lowered with respect to the ones at m c = 1 reported in Ref. [28], suggesting that the increase in mass enhances the orientational ordering at fixed activity. In particular, we estimate φ c = 0.730 ± 0.01 and 0.760 ± 0.01 for 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 Pe = 10 and φ = 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 φ c (m c ) = a + be −mc/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 Pe = 0 [71]. When the persistence number is pn 10 −2 (m c ≈ 10), the system behaves closer to the passive case. On the other hand, when pn 10 −2 the active force has a disordering effect in the hexatic ordering.

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 Sect. 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 Fig. 9 Liquid-Hexatic transition in the ABP and AHP model. Global hexatic parameter as a function of the global packing fraction for Pe = 10 and mc = 44,. The orange and blue curves correspond to simulations with and without hydrodynamics, respectively, for active colloids with the same mass 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 Sect. 2.2.3, fix the colloid mass to m c = 44 and γ = 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 φ between 0.5 and 0.85.

Liquid-hexatic transition
We start by looking at how the ordering properties are affected by hydrodynamics. Figure 7a-c show, for three different densities at Pe = 10, the color map of the local hexatic parameter ψ 6,j projected onto its average value. In panel (a) (φ = 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), φ = 0.8), although global ordering is not observed. At φ = 0.81, panel (c), a single fully hexatically ordered system is observed. Thus, also AHP present a transition between liquid and hexatic phases. Figure 8 shows the hexatic correlation functions varying the density for 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 shifts from an exponential decay to a power-law decay at substantially higher values of packing fraction φ. More precisely the transition is located at φ c ≈ 0.805 ± 0.01 for Pe = 10, and φ c ≈ 0.840 ± 0.01 for Pe = 20.
The increase in value of the transition density φ 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 75 Page 10 of 15 Eur. Phys. J. E (2022) 45 :75 opposite to the effect of increasing the particles mass, which instead promotes hexatic ordering. Indeed, if we measure the average global hexatic parameter ψ 6 = 1 N | N i ψ 6,i | as a function of the global packing fraction φ (Fig. 9 ), which increases from 0 to 1 as the liquidhexatic transition is crossed, we find that the transition is significantly shifted. Note that both curves converge to almost the same value 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 [8] are not affected by the presence of hydrodynamic interactions. However, hydrodynamics produces other relevant effects which will be now discussed.

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. Figure 10 shows the coarse-grained steady state velocity fields of the fluid, v(r) (panels d-f) along with associated snapshots of the configurations colored according to the hexatic parameter (panels a-c), for AHP with Pe = 10 and three different values of packing fraction φ. Coarse-grained velocity fields of the fluid are realized by averaging the velocity of fluid particles inside blocks of size 4σ c = 20σ s (such large coarsegraining cells are chosen for the sake of visualization; similar profiles can be obtained with smaller cells). The first density, φ = 0.60 (panels (a), (d)), is characterized by the absence of 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. Figure 10b, e show instead a larger density φ = 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 φ = 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 (φ = 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 along 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 [72,73], 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 [74][75][76][77], flowing crystals made of spontaneously aligning self-propelled hard disks [78] and self-sustained spontaneous flows in active gels [79][80][81][82][83].
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: (25) Figure 11 shows C v (r) for different values of φ, for Pe = 10. For low values of φ (see e.g. φ = 0.600) the curve shows an exponential decay. This corresponds to the case shown in Fig. 10d, 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 φ 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 φ. In the inset of Fig.  11 the velocity correlation function for ABP in the hexatic phase (Pe = 10 and φ = 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 φ = 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 show the corresponding steady state fluid velocity field. The color code is the same as the one in Fig. 3 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 φ = 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 an 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. 13a.

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 (Pe = 10, 20). We have first characterized the ABP by changing their mass, while maintaining the same Pe and D θ , so that we have a non-trivial interplay between the inertial time and the persistence time t p = 1/D θ . 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 Pe = 0, suggesting that inertia has an orientational ordering effect 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 Pe = 10 and φ = 0.795, before (panel (a)) and after (panel (b)) switching off hydrodynamics. The corresponding colloids velocity fields are shown in panels (c-d) 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 φ, 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 φ ≈ 0.72 of small regions of correlated velocity field, characterized by the presence of vortices, that are not associated to any local orientational ordering; ii) the arisal above φ ≈ 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 changes 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. [8], 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 selfsustained 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 to study in more details particle-particle flow interactions and local velocity field effects. We hope that our results can boost further research in this direction.

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.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.