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 [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 micro-swimmers [13, 14], 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 [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 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 [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 liquid-hexatic 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.

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:

$$\begin{aligned} 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) \end{aligned}$$
(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 \ge 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_\mathrm{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:

$$\begin{aligned} m_c\ddot{{\varvec{r}}_i}&= - \gamma \dot{{\varvec{r}}_i} + F_{\text {act}}\mathbf {n}_i - {\varvec{\nabla _i}}\sum _{j\ne i}U(r_{ij}) + {\varvec{\xi }}_i, \end{aligned}$$
(2)
$$\begin{aligned} \dot{\theta }_i&= \eta _i \ , \end{aligned}$$
(3)

where \(i=1,..,N_c\), \(r_{ij}=|\mathbf {r}_i-\mathbf {r}_j|\) and \(\gamma \) is the damping coefficient. The terms \({\varvec{\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:

$$\begin{aligned}&\langle \xi _{i\alpha }(t) \rangle = 0,\; \langle \eta _{i}(t) \rangle = 0, \end{aligned}$$
(4)
$$\begin{aligned}&\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), \end{aligned}$$
(5)
$$\begin{aligned}&\langle \eta _{i}(t_1) \eta _{j}(t_2) \rangle = 2D_\theta \delta _{ij}\delta (t_1-t_2), \end{aligned}$$
(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_BT\), is the Péclet number Pe = \(F_\mathrm{act} {\sigma _c}/(k_BT)\). 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}\) [48].

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 [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 \(\gamma \) and vary the disk mass \(m_c\) to consider different inertial contributions, \(k_BT=0.05\) and \(D_\theta = 3 k_BT/(\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 LAMMPS [49] to integrate numerically the equations of motion, using a timestep \(\varDelta t_c=0.001\) and periodic boundary conditions. We fix the Péclet number to \(\text {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 [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.

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 position \(\mathbf{r}_i\) and velocity \(\mathbf{v}_i\), both of which are continuous variables. In this algorithm, the time is discretized in units \(\varDelta t_s\), and the evolution of the system is composed by two steps, propagation and collision, which are applied consecutively for each \(\varDelta t_s\).

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

$$\begin{aligned} \mathbf{r}_i(t+\varDelta t_s)=\mathbf{r}_i(t)+\mathbf{v}_i(t) \varDelta t_s . \end{aligned}$$
(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 [45, 50]

$$\begin{aligned} \mathbf{v}_i(t+\varDelta t_s)=\mathbf{u}(t)+\Omega [\mathbf{v}_i(t)-\mathbf{u}(t)] , \end{aligned}$$
(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 to [51, 52]:

$$\begin{aligned} \nu _s= & {} \frac{\sigma _s^2}{2 \varDelta t_s} \left[ \left( \frac{\lambda }{\sigma _s} \right) ^2 \left( \frac{n_s}{(n_s-1+\exp {(-n_s)})\sin ^2(\alpha )}-1 \right) \right. \nonumber \\&\left. +\frac{(n_s-1+\exp {(-n_s)})(1-\cos (\alpha ))}{6 n_s} \right] , \end{aligned}$$
(9)

while the coefficient \(D_s\) is [53]:

$$\begin{aligned} D_s=\frac{\lambda ^2}{2 \varDelta t_s} \left( \frac{2n_s}{(n_s-1+\exp {(-n_s)})(1-\cos (\alpha ))} \right) , \end{aligned}$$
(10)

where \(n_s=N_s \sigma _s^2 / L^2\) is the average number of particles per cell and \(\lambda =\varDelta 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 [46]; here we adopted the one implemented in the LAMMPS software [47].

In this implementation colloids are evolved for n timesteps \(\varDelta t_c\), following the equation of motion (2) without the force terms \({\varvec{\xi }}_i\) and \(\gamma \dot{ {\varvec{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 \(\varDelta t_s=n\varDelta t_c\). Note that both \(\varDelta t_c\) and \(\varDelta 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 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 \(\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}' \varDelta t_s \varepsilon \), where \(\mathbf{v}'\) is the updated velocity and \(\varepsilon \) is a uniformly distributed random number in the interval [0, 1] [59]. The new velocity \(\mathbf{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

$$\begin{aligned} p_N(v_N)&= (m_s v_N/k_B T) \exp {(-m_s v_N^2/2k_B T)} , v_N>0 \end{aligned}$$
(11)
$$\begin{aligned} p_T(v_T)&= \sqrt{m_s/2\pi k_B T} \exp {(-m_s v_T^2/2k_B T)} , \end{aligned}$$
(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 \(\varDelta \mathbf{P}=\sum _s m_s (\mathbf{v}-\mathbf{v}')\) and \(\varDelta \mathbf{L}=\sum _s m_s (\mathbf{r}_d-\mathbf{R}) \times (\mathbf{v}-\mathbf{v}')\), and the linear and angular velocities of the colloid are updated as: \(\mathbf{V}'=\mathbf{V}+\varDelta \mathbf{P}/m_c\) and \(\mathbf{\omega }'=\mathbf{\omega }+\varDelta \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 \(\varDelta 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 \simeq 10\) is the best choice to optimize computational speed and accuracy.

Fig. 1
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)

Full size image
Fig. 2
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}\)

Full size image

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}'\) 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) [47].

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 [60]. 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 [44].

Fig. 3
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, \({\varvec{\Psi }} = 1/N \sum _j \psi _{6,j}\), at fixed \(\text {Pe} = 10\) and \(m_c = 44\) for \(\phi ={0.710, 0.730, 0.760}\) respectively, for a system of size \(L=256\sigma _c\)

Full size image

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 \), \(\varDelta 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 [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 as \(Sc=\nu _s/D_s\). Values \(Sc \sim O(10)\) can be obtained by requiring small values of \(\lambda \) and large rotation angle [62]. Note that the choice \(\lambda < \sigma _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 \(\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 form [61]:

$$\begin{aligned} Sc \simeq \frac{1}{12 (\lambda /\sigma _s)^2} \ , \end{aligned}$$
(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:

$$\begin{aligned} \gamma = C_{2D}\pi \nu _s\rho _s (\sigma _c/2), \end{aligned}$$
(14)

where \(\rho _s=n_s m_s/\sigma _s^2\) 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\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 [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 \(\text {v}_s\) due to the external forces (in our case due to activity) and the sound velocity \(\text {v}_\text {sound}=\sqrt{2 k_B T/m_s}\) inside the fluid:

$$\begin{aligned} Ma=\frac{\mathrm {v_s}}{\mathrm {v}_{sound}}. \end{aligned}$$
(15)

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 Pu, instead, is the ratio between the active stationary colloid velocity \(F_{act}/\gamma \) and the fluid self-diffusion:

$$\begin{aligned} Pu=\frac{\sigma _c F_{act}}{6 \gamma D_s}, \end{aligned}$$
(16)

and should be less than 1 [40] 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 [59]. 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=4 m_c/(\pi \sigma _c^2)\). 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 \(\varDelta t_c=10^{-4}\) and \(\varDelta t_s=410 \varDelta 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 \(\mathrm{Pe}=20\) value considered. We note that the Reynolds number of the fluid is given by

$$\begin{aligned} Re=\text {v}_s \sigma _c / \nu _s=F_{act} \sigma _c / \gamma \nu _s=6 Pu D_s / \nu _s = 6 Pu / Sc , \end{aligned}$$
(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_BT/m_s\) and \(k_BT/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 \(\text {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\) [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\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. [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 \simeq 100\). We considered large systems, with \(L=900 \sigma _c\) to reduce periodic boundary effects.

Fig. 4
figure 4

Hexatic order correlation functions for the ABP model. Hexatic order correlation functions \(g_6(r)\) for \(m_c=44\) at \(\text {Pe} = 10\) (left) and \(\text {Pe} = 20\) (right) for different global packing fractions given in the keys

Full size image
Fig. 5
figure 5

Effects of colloids mass on the Liquid-Hexatic transition. On the left panel, orientational correlation functions, \(g_6(r)\), at fixed \(\mathrm 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 \(\mathrm 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 \(\mathrm Pe=5\); \(a=0.72,\, b=0.08,\, c=12.92\) for \(\mathrm Pe=10\); and \(a=0.76,\, b=0.08,\, c=14.30\) for \(\mathrm Pe = 20\). The error bars correspond to the gap \(\varDelta \phi \) between the densities scanned in the simulations for each value of the mass \(m_c\)

Full size image

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

$$\begin{aligned} C_u(t)=<u(t)u(0)>=\frac{k_B T}{m_c}exp(-t/t_{E}) \ \ , \end{aligned}$$
(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 [67]

$$\begin{aligned} \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}\ \ . \end{aligned}$$
(19)

The integral of the VACF is related to the diffusion coefficient \(D_c\) through the Green-Kubo relation,

$$\begin{aligned} D_c=\int _0^{\infty }<u(t)u(0)> dt =\frac{k_B T}{\xi } \ \ . \end{aligned}$$
(20)

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

$$\begin{aligned} C_u(t)=\left( \frac{1}{2 \rho _s}\right) \frac{k_B T}{[4\pi (D_c+\nu _s)t]} \ \ , \end{aligned}$$
(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 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

$$\begin{aligned} C_u(0)=\frac{k_B T}{m_c} , \end{aligned}$$
(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 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_{\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 [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\)

Fig. 6
figure 6

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

Full size image

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_{\nu }\) (panel (c)).

Fig. 7
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 \(\text {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\)

Full size image

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{aligned} \begin{aligned} 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'}dt' \approx \frac{k_B T}{8\pi \rho _s\nu _s}[\ln t]_{t_{\nu }}^{t}. \end{aligned} \end{aligned}$$
(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\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 \(\text {Pe}=10\) and at \(\phi _c=0.83\) for \(\text {Pe}=20\) [71].

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:

$$\begin{aligned} g_6(r) =\frac{\langle \psi _6({\varvec{r}}_i)\psi _6({\varvec{r}}_j) \rangle }{ \langle \psi _6^2({\varvec{r}}_i) \rangle } \ , \end{aligned}$$
(24)

where \(r = | {\varvec{r}}_i - {\varvec{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 Sect. 3.2.2 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 Sect 2.1 at \(\mathrm 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 \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\)  [1] 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.

Figure 3a–c show typical configurations at \(\mathrm 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. Figure 4 shows these functions for \(m_c=44\) and \(\mathrm Pe = 10, 20\). At densities below \(\phi =0.72\) for \(\mathrm Pe = 10\) and \(\phi =0.74\) for \(\mathrm 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 \(\phi _c= 0.730\pm 0.01\) and \(0.760\pm 0.01\) for \(\mathrm 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 \(\mathrm 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 \(\mathrm Pe = 0\) [71]. When the persistence number is \(pn > rsim 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.

Fig. 8
figure 8

Hexatic order correlation function for the AHP model. a–b Hexatic order correlation function \(g_6(r)\) at \(\text {Pe} = 10\) (d) and \(\text {Pe} = 20\) (e) for different global packing fractions given in the keys and \(m_c=44\)

Full size image
Fig. 9
figure 9

Liquid-Hexatic transition in the ABP and AHP model. Global hexatic parameter as a function of the global packing fraction for \(\text {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

Full size image

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 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 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 \(\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. Figure 7a–c show, for three different densities at \(\text {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.

Figure 8 shows the hexatic correlation functions varying the density for \(\text {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 \(\phi \). More precisely the transition is located at \(\phi _c\approx 0.805 \pm 0.01\) for \(\mathrm{Pe}=10\), and \(\phi _c\approx 0.840\pm 0.01\) for \(\mathrm{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 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.

Fig. 10
figure 10

Self-sustained active flow. a–c Color maps of the local hexatic order parameter for \(\text {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

Full size image

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.

Figure 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 \(\text {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 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 \(\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 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:

$$\begin{aligned} C_{\text {v}}(\mathbf {r})=\frac{\langle \mathbf {v(r)v(0)}\rangle }{\langle \mathbf {v(0)}^2\rangle } \ \ . \end{aligned}$$
(25)

Figure 11 shows \(C_{\text {v}}(\mathbf {r})\) for different values of \(\phi \), for \(\mathrm{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. 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 \(\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 (\(\mathrm{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 \(\mathrm{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 \(\mathrm{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 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.

Fig. 11
figure 11

Spatial velocity correlations. Spatial velocity correlation functions \(C_{\text {v}}(r)\), for different values of \(\phi \), for \(\mathrm{Pe}=10\). In the inset the velocity correlation function of ABP at \(\mathrm{Pe}=10\) and \(\phi =0.760\) is compared with \(\phi =0.60\) for AHP at the same \(\mathrm{Pe}\)

Full size image
Fig. 12
figure 12

Radial distribution function. Radial distribution function for ABP and AHP models with \(\phi =0.750\) at \(\mathrm{Pe}=10\)

Full size image
Fig. 13
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 \(\text {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)

Full size image

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 (\(\text {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 \(\text {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 \(\mathrm{Pe}=10\) two results: i) the formation below \(\phi \approx 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 \(\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 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 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 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.