Phase Separation and Coexistence of Hydrodynamically Interacting Microswimmers

A striking feature of the collective behavior of spherical microswimmers is that for sufficiently strong self-propulsion they phase-separate into a dense cluster coexisting with a low-density dis- ordered surrounding. Extending our previous work, we use the squirmer as a model swimmer and the particle-based simulation method of multi-particle collision dynamics to explore the influence of hydrodynamics on their phase behavior in a quasi-two-dimensional geometry. The coarsening dynamics towards the phase-separated state is diffusive in an intermediate time regime followed by a final ballistic compactification of the dense cluster. We determine the binodal lines in a phase diagram of P\'eclet number versus density. Interestingly, the gas binodals are shifted to smaller densities for increasing mean density or dense-cluster size, which we explain using a recently introduced pressure balance [S. C. Takatori et al., Phys. Rev. Lett. 113, 028103 (2014)] extended by a hydrodynamic contribution. Furthermore, we find that for pushers and pullers the binodal line is shifted to larger P\'eclet numbers compared to neutral squirmers. Finally, when lowering the P\'eclet number, the dense phase transforms from a hexagonal"solid"to a disordered"fluid"state.


Introduction
Collective dynamics due to the active motion of microorganisms and artificial microswimmers has received a lot of attention among physicists 1-6 as it is relevant both to real world applications 2,6-8 as well as for posing fundamental questions in non-equilibrium statistical physics 4,[9][10][11][12][13] . From the perspective of physics, a unifying feature of microorganisms and artificial microswimmers is that they propel themselves autonomously through a surrounding fluid. As they are constantly driven out of equilibrium, understanding the collective properties of systems comprised of many swimmers is an important field of ongoing research in statistical physics.
One of the most striking features of swimmer systems is motility-induced phase separation 14 , where stable clusters of swimmers form due to their active motion alone and without any attractive forces. The question of motility-induced phase separa-tion has been treated theoretically using so-called active Brownian particles, which perform a persistent random walk and interact sterically only 11,[15][16][17][18][19] . However, real microswimmers such as ciliated microorganisms 20 , catalytic Janus particles 7,8,[21][22][23] , or active emulsion droplets 2,24-27 employ propulsion mechanisms reliant on hydrodynamics. Thus nearby microswimmers can affect each other's velocity and orientation, which can potentially alter their collective behaviour depending on the details of their hydrodynamic interactions 28 .
Biological microswimmers employ non-reciprocal cell-body deformations or the collective dynamics of flagella or cilia in order to propel themselves 29 , whereas active colloids and droplets move forward by creating a slip velocity close to their surface 1 . The essential features of both of these self-propulsion mechanisms are captured by the so-called "squirmer" model [30][31][32][33][34] , which we shall use to generate our hydrodynamic propulsion.
For phase separating systems, it is common to construct the binodal line in order to determine the phase coexistence region. While there have been many studies quantifying the phasecoexistence regime for active Brownian particles, 11,13,15,[35][36][37][38] to list just a few of them, none of these have considered hydrodynamic interactions between the swimmers. Several recent works studied the collective dynamics of hydrodynamic swimmers [39][40][41][42][43] , but did not look at phase separation.
In our previous study, we examined the collective dynamics and clustering of hydrodynamic squirmers in a quasi-two-dimensional geometry using the method of multi-particle collision dynamics (MPCD) to explicitly simulate the fluid environment 28 . Due to computational limitations, we were restricted to 208 swimmers, which was insufficient to quantify the first-order phase-separation transition. Hence, we built on this study and extended the simulations to thousands of squirmers by implementing a parallelized version on 1008 processor cores. This allows us to investigate the dynamics of a sufficiently large number of squirmers in order to quantitatively resolve phase separation into a dilute gas and a dense phase.
In the following we show phase-separated squirmer systems and explain how we analyze them to obtain the area fractions of the dilute and dense phases. After illustrating the coarsening dynamics towards phase separation, we quantify the coexistence region bounded by the binodal line in a diagram Péclet number versus density. Interestingly, the binodal lines differ for varying mean density. We explain this phenomenon as a true hydrodynamic signature using a pressure balance. In particular, we use the swim pressure introduced in Ref. 9 and extend the approach of Refs. 9,10 by including a hydrodynamic pressure generated by squirmers at the rim of the dense cluster. The binodal line is shifted towards larger Péclet numbers, when the force-dipole contribution of the squirmer, characterizing pushers and pullers, is increased. Finally, moving the dense phase from large Péclet numbers to the critical value, one observes a melting transition between a hexagonal and a fluid phase. We start with summarizing computational details.

System and Methods
Our system consists for N spherical squirmers of radius R, where N varies from 3328 to 4624 depending on mean area fraction or densityφ . Squirmers propel themselves by applying an axisymmetric surface-velocity field 33,34 , where v s is the slip velocity, r s points from the squirmer centre to its surface,r s is the unit vector along r s , andê is the squirmer's orientation. The constant B 1 determines the amplitude of the slip velocity and thus its swimming speed v 0 . For a single squirmer in a bulk fluid, v 0 = 2B 1 /3. We use B 1 to control the Péclet number, Pe = Rv 0 /D T , as the appropriate thermal diffusion coefficient D T is kept constant for all simulations. For the MPCD fluid and squirmer parameters presented below, we obtain Pe ≈ 3546.4B 1 . The additional parameter β breaks the symmetry about the squirmer equator. For β < 0, the surface field is more concentrated on the southern hemisphere and the squirmer's far field in bulk is that of a pusher. For β > 0 a puller is realized and a swimmer with β = 0 is called a "neutral squirmer".
Motivated by experiments 2,24,44 and following our previous study 28 , the swimmers are bounded by two parallel flat plates. As in Ref. 28 we choose strong confinement with the distance between the plates being H = 8R/3. We implement the full hydrodynamic interactions between squirmers and confining walls using the MPCD simulation technique 45,46 . The fluid is represented by of the order of 10 7 point-like effective fluid particles with mass m 0 , velocity v i (t), and position r i (t). The velocities v i (t) are ther-mostatted to have a temperature T 0 . The fluid particles perform alternating streaming and collision steps. During the streaming step, fluid particles move ballistically over a time ∆t without interacting with each other and reach the position They do interact with confining boundaries and squirmers, transferring linear and angular momentum. Interactions with the confining boundaries implement no-slip boundary conditions via the so-called bounce-back rule 47 , while interactions with squirmers generate the surface field of Eq. (1). For each collision step, particles are sorted into a cubic grid of cell length a 0 and with a random offset for each new step. The average velocityv = v i cell and centre of massr = r i cell of the particles in each cell is determined. Within each cell they "collide" according to the MPCD-AT+a rule 45 , which is reminiscent of an Andersen thermostat and also conserves linear and angular momentum. Random velocities increments δ v i for each particle are sampled from a Gaussian distribution with variance k B T 0 /m 0 . The velocity increments would give an overall change in momentum m 0 ∆v = m 0 δ v i cell and angular momentum ∆L = m 0 (r i −r) × δ v i cell . Hence, we subtract these changes in the velocity increment for the collision step, where I −1 is the inverse of the moment of inertia tensor for the point-mass distribution of the particles in the cell's centreof-mass frame. Thus, the motion of the fluid particles on the length scale of the grid corresponds to solutions of the Navier-Stokes equations 45,48,49 including thermal noise. Furthermore, this scheme accurately reproduces the hydrodynamic flow field of a squirmer 34,50 and their interactions including near-field lubrication effects 28,50 .

System Parameters
For all simulations discussed here, the fluid and squirmer mass densities are the same, ρ = 10m 0 /a 3 0 , and the squirmer radius is R = 3a 0 .
Together with the streaming-time step ∆t = 0.02a 0 m 0 /k B T 0 , the fluid viscosity becomes η = 16.05 √ m 0 k B T 0 /a 2 0 taken from 28,48 . For the bulk fluid and the parameters stated above, the self-diffusivity of a passive colloid is In the presence of walls, the self-diffusivity becomes a position-dependent tensor 51 . We measure the self-diffusivity at the centre between the walls using the Green-Kubo formula. We find that D T,2D ≈ 0.5D T,B , comparable to the result of Ref. 51 . For a colloid close to the walls this value can decrease to approximately 0.25D T,B due to the enhanced friction. In the following, we use D T,2D for a single particle (particle density φ = 0) to define the Péclet number, The orientation of a squirmer also undergoes rotational diffusion, for which MPCD gives the rotational diffusion coefficient also approximately valid in our quasi-two-dimensional geometry. There, rotational diffusion strongly depends on the density of squirmers due to the flow field they create and increases by as much as a factor of 15 in very dense systems 28 . The rotational diffusion coefficient is used to define the persistence number Pe r = v 0 /RD R , which compares persistence or run length To summarise, in the following we work with Péclet numbers in the range [190,355], which corresponds to persistence numbers Pe r in the range [120, 220] or typical run lengths of single squirmers, l = v 0 /D R , ranging from 120R to 220R. For high-density simulations without phase separation, where rotational diffusion is considerably enhanced, we have Pe r 8 or l 8R.

Phase Separation and Coarsening Dynamics
For a range of total densities and Péclet numbers, the system decomposes into a dilute gas and a dense phase. Figure 1: top, left shows a typical snapshot of a phase-separated system with mean area fractionφ = 0.64 and at Péclet number Pe = 355 (B 1 = 0.1). The two phases are clearly visible: many squirmers have come together to form a densely packed and system-spanning cluster, leaving a few squirmers in a much less dense "gas" phase. For lowerφ , the dense phase consists of fewer squirmers and thus might not necessarily form system-spanning clusters. At even lowerφ phase separation does not occur and only the pure gas phase is visible. The supplemental material † contains videos M1 and M2b showing a non phase-separating and phase-separating squirmer system, respectively. As well as the onset of phase separation M2a.
To quantify the local density, we construct around each squirmer a Voronoi cell and determine its area A i . The local density within this Voronoi cell is therefore φ i = 1/A i . A standard analysis, such as the one used in 35 , would be to determine the probability distribution p(φ i ) of local densities φ i and to identify the densities of the two coexisting phases from the expected bimodal distribution. The probability distribution corresponding to the squirmer snapshot in Fig. 1 is shown in Fig. 1: bottom by the blue line. The dense phase produces a very pronounced peak with a very broad tail, which obscures the signature of the gas phase. This broad tail is due to encounters between a small number of swimmers in the gas phase forming small temporary clusters. Thus, some form of temporal averaging is needed, in order to keep only those dense regions that remain intact and relatively fixed in place.
Hence we divide the system area into a grid with cell spacing equal to the squirmer radius R. Each cell is assigned a local density using a weighted average of the local densities from those Voronoi cells, which intersect each grid cell. We then average over time, choosing the averaging interval to be longer than the life time of temporary clusters, yet short enough so that permanent clusters seem relatively unchanged. Hence, we arrive at a position-resolved density φ (x, y) such as the one plotted in Fig. 1: top, right. A histogram of the values of φ (x, y) is plotted in Fig. 1: bottom,green line, and recovers the peak for the gas phase.
The color-coded density φ (x, y) in Fig. 1 also shows a pronounced interfacial region between dense and gas phase. In order to obtain their densities and distinguish them from the interface, we determine the interfacial density profile φ (d), where |d| is the distance to the nearest point on the interface. A preliminary position of the interface is determined by requiring that the timeaveraged number of neighbours is three. Variations in this choice simply shift the density profile along the d-axis, but do not alter its shape. Thus the resulting φ (d) is independent of fixing the exact interface position beforehand. For each grid cell the distance to each point of the estimated interface is determined and the minimal value is chosen as |d|. Collecting the statistics for each grid cell then gives φ (d). The resulting profile corresponding to the color-coded density φ (x, y) can be seen in the inset of Fig. 1: bottom. Inspired by Ref. 10 and Cahn-Hilliard theory 13 , we fit to the profile in the inset (cf. solid line). This provides values for the respective densities of gas and dense phase, φ gas and φ den . The good quality of the fit is not surprising, as an effective Cahn-Hilliard equation was formulated for the local density in systems comprised of active Brownian particles 11,17 . Before we introduce the resulting phase diagram, we discuss the coarsening dynamics from the uniform initial state to the final phase-separated system. Figure 2: top illustrates the time evolution for a system withφ = 0.48 and Pe = 355 (B 1 = 0.1) on different time scales. Initially (first row) clusters of squirmers develop in the unstable uniform phase characteristic for spinodal decomposition. They grow and become denser until a few clusters remain (second row). These clusters already have the final density φ den of the dense phase. They show the typical coarsening dynamics of a phase-separating system due to diffusive exchange of squirmers between the clusters. Ultimately, one single cluster remains (third row), the shape of which becomes more compact over the course of time.
Following Ref. 52 , we quantify coarsening by determining the area a(t) of all dense regions and their total perimeter l(t) as a function of time. The mean domain size then becomes ξ (t) := a(t)/l(t). Figure 2: bottom compares ξ (t) to the different regimes of the coarsening dynamics for a system with conserved order parameter including hydrodynamics 52 . After the initial formation of the clusters (the magenta dots refer to the frames in the first row), the total area a(t) remains constant and all further changes in ξ are due to changes in l(t) alone. In the intermediate regime, the mean domain size roughly obeys the scaling expected when diffusive transport dominates, ξ ∼ t 1/3 (frames in the second row). This regime lasts less than a decade, since our system size is relatively small. Former work on the coarsening dynamics of pure active Brownian particles report a deviation from the scaling exponent 1/3, measuring 0.255 35 or 0.28 53 . However, our system is too small to make a clear statement about the exponent. Finally, the compactification of the single cluster gives rise to clearly visible, jumplike increases of ξ (t). For comparison, we give the expected scaling ξ ∼ t due to viscous hydrodynamic transport. This compactification has not been seen for pure active Brownian particles. Also, in our simulations not allφ enter this regime. Small clusters leave the diffusive transport regime with fairly circular clusters already. Whereas forφ ≥ 0.50 clusters spanning the whole simulation box are formed. These are stable, either as "slabs" forφ ≈ 0.50 or the gas phase forms a "bubble" at even larger mean densities. The rapid increase of ξ (t), shown in Fig. 2: bottom towards the end, therefore requiresφ 0.50 and probably is a finite size effect.

Phase diagram of neutral squirmers
In Fig. 3: top we plot the densities φ gas and φ den of the coexisting gas and dense phases in a Pe − φ diagram, which we determined for different mean densitiesφ . Interestingly, the resulting binodal lines do not agree. In particular, the gas binodal is shifted to smaller values φ gas for increasingφ . In contrast, for pure active Brownian particles, the densities of the two coexisting phases are independent ofφ 11,35 . For our geometry we present them in appendix 6.2. Thus, we suspect hydrodynamic interactions between the squirmers to cause this behaviour. For constant Péclet number Fig. 4 demonstrates how the gas density φ gas decreases for larger φ . Thus, the cluster of dense phase, which grows in size with increasingφ , extracts additional squirmers from the gas phase. Below we will attribute this behaviour to an additional hydrodynamic pressure caused by the squirmers at the interface between the dense and the dilute phase. The gas binodals in Fig. 3: top roughly run parallel to each other. Thus, by introducing an effective Péclet number with fit parameter a ≈ 0.65, we are able to collapse them on a single master curve for the binodal line as Fig. 3: bottom demonstrates. It is tempting to explain the coefficient a by defining an effective Péclet number using a density-dependent diffusion coefficient D T,2D (φ ) = D T,2D (φ = 0)/(1 + 2.027φ ) in the definition (4) of the Péclet number. We determined it from the velocityauto-correlation function and the Green-Kubo formula (see also paragraph after Eq. (8) and Appendix 6.3). However, the prediction of 2.027 for the coefficient a cannot explain the reported data collapse. In order to estimate the position of the critical point (φ * , Pe * ), we fit 11 φ to the gas-binodal in Fig. 3: bottom. The supplemental material † contains videos M3 and M4 of the dynamics of the positionresolved density φ (x, y) inside the phase-coexistence regime and outside of it close to the critical point.
For a system in mechanical equilibrium the total pressure in the coexisting gas and dense phases has to be the same. We now apply this strategy to our system of phase-separating squirmers in order to understand why the binodals depend on the mean density. For active Brownian particles without any hydrodynamic interactions, the pressure balance was first proposed in Ref. 9 by introducing a total pressure consisting of a steric and an active contribution: p = p (s) + p (a) . This pressure balance was then used in the context of phase separation in Refs. 9,10,16,[54][55][56] . The steric pressure results from direct interactions between the squirmers, p (s) = ∑ i j r i j ·F i j /(6V ) 57 . Here, r i j = r i −r j is the distance vector from swimmer j to swimmer i, F i j the force with which squirmer j acts on i, V is the volume of a uniform system, and ... means temporal average. Swimmers moving with constant velocity v 0 along the unit vectorê i also generate a swim pressure 9 : where γ 0 is a friction coefficient connected to energy dissipated by the swimmer in its fluid environment. We note here that the swim pressure is only defined over times, which are long compared to 1/D R , therefore the average · · · needs to be taken over sufficiently long times 58 . As swimmers move, while this averaging is taking place, p (a) becomes "smeared out" and therefore describes the active pressure in the volume V and not at a point. * Due to hydrodynamic interactions with bounding walls and other squirmers, the friction coefficient in Eq. (8) depends not only on the squirmer's position but also the position of every other squirmer. Hence, it is not possible to know γ 0 precisely. To arrive at a sort of mean-field approximation for γ 0 (φ ), we use the Green-Kubo formula and determine a self-diffusion coefficient D T,2D (φ ) from the velocity-auto-correlation function of a passive particle with radius R using the velocity components in the plane of the quasi-two-dimensional geometry. The result is plotted in Fig. 11 of Appendix 6.3 and can be well approximated by D T,2D (φ ) = D T,2D (φ = 0)/(1+2.027φ ). We then use the Einstein relation, γ 0 (φ ) = k B T /D T,2D (φ ), to arrive at a density-dependent friction coefficient. However, the following results only change little compared to working with a constant γ 0 = k B T /D T,2D (φ = 0) for all densities.
In a dilute gas of non-interacting Brownian particles, orientation decorrelates exponentially on the decorrelation time 2τ r and one obtains from Eq. (8) 9 p (a) id = ρ gas τ r γ 0 v 2 0 /6 ∼ ρ gas Pe 2 , where ρ gas is number density. This active gas pressure is controlled by the Péclet number. It holds the dense phase of active Brownian particles together by preventing them from swimming apart.
We are able to calculate steric and the swim pressure of Eq. (8) for the local neighbourhood of each squirmer (see Appendix in section 6 for details). Figure 5: top shows these local pressure values as a function of the distance from the gas-cluster interface. In light of the discussion about the swim pressure following Eq. (8), its pressure values close to the interface have to be treated with caution. But they are not important for the following. Outside the dense phase, d < 0, the steric pressure is low due to very few contacts between squirmers or a squirmer and a bounding wall. And when they do occur, they persist for only a short time. On the other hand, the local swim pressure in the gas phase is large as swimmers can move unhindered. Inside the * One can argue that the linear exention of the volume V should be larger than the persistence length of the swimmer. However, if one averages over a sufficiently long time so that a sufficient number of swimmers enter and leave V , this requirment might be relaxed. This needs to be checked more carefully in a future investigation.  cluster, d > 0, squirmers are in constant contact with their neighbours and the confining walls. This leads to a significantly higher steric pressure. On the other hand, the squirmer's swimming motion is hindered by the neighbours. Swimming motion and actual displacement are much less correlated and the swim pressure is greatly reduced. However, contrary to what is expected for mechanical equilibrium, the total pressure p = p (s) + p (a) plotted in Fig. 5: top is not the same for the gas and the dense phase. Furthermore, in Fig. 5: bottom we see that this total pressure difference ∆p = p den − p gas across the phase interface increases with increasingφ . Thus, there must be an additional negative pressure stabilizing this interface. When introducing the swim pressure, the authors of 9 already remarked that for swimmers swimming in a fluid environment, the hydrodynamic stresslet induces an additional hydrodynamic pressure. In the following we elaborate on this idea for our special geometry and demonstrate that indeed we are able to justify a negative pressure difference between the dense and the gas phase.
It is not possible to give a complete quantitative theory for the negative hydrodynamic pressure produced by the dense cluster of squirmers. Instead, we present here an approximate treatment, which contains the essential ideas. A squirmer at position r 0 and confined between two plates generates the flow field of a source dipole with the approximate dipole moment σ σ σ = 2πR 2 v 0ê 40,59 . It is accompanied by the pressure field where G ∼ H 2 /η measures the system's confinement with H the plate distance and η viscosity. To be concrete, we use the value for a Poiseuille flow in MPCD units, G = H 2 /12η ≈ 1/3. Now, we consider a circular cluster of squirmers and replace the dipole moment σ σ σ by a density of source dipoles. Furthermore, we find that on average the squirmers have a tendency to point radially inwards towards the center, especially at the rim of a cluster. Therefore, we introduce a source dipole densityσ σ σ = σ σ σ ·r c r c , wherer c is a unit vector pointing from the squirmer to the cluster center. The radial dipole density σ σ σ ·r c across the interface from the gas to the dense phase is plotted in Fig. 6 for differentφ . This confirms the preferred radial alignment at the rim of the dense cluster, whereas in the interior the squirmers point towards the bounding walls, on average, since σ σ σ ·r c is zero. We calculate the pressure generated by the dipole density in two steps. First, we integrate along a circle of radius r 0 to obtain pressure per unit length: With r 0 = (r 0 cos ϕ, r 0 sin ϕ), cos ϕ = r 0 · r, and r 0 = −r 0rc one obtains which can be integrated to Thus inside a ring of source dipoles pointing radially inwards, pressure is constant and negative, while outside of the ring it is zero. Second, we integrate over the radial coordinate and obtain the constant hydrodynamic pressure difference between cluster interior and exterior, where R c is the cluster radius. This negative pressure should cancel the difference of the total pressure between dense and gas phases, which is illustrated in Fig. 5: top.
In summary, our argument is as follows: any net polar order of swimmers is accompanied by a jump in solvent pressure, a fact that has also been reported by Ref. 60 . In our case, this is due to the squirmers at the edge of the cluster: the hydrodynamic pressure difference between cluster interior and exterior is necessary to stop the flow initiated by the squirmers, which try to pump fluid out of the cluster interior.
The inset in Fig. 6 shows the hydrodynamic pressure using the radial source-dipole density σ σ σ ·r c from the main plot. Its absolute value increases linearly with cluster size orφ , which results from the stronger radial alignment as the cluster size increases. The linear increase of |∆p (h) | coincides with the linear increase of ∆(p (s) + p (a) ) plotted in Fig. 5: bottom, however the slope is by a factor of 1.1 too low. This is understandable since we approximate the squirmers in the dense cluster by point-like source dipoles with an approximate dipole moment. We only calculated ∆p (h) up toφ = 0.48, where the dense clusters are nearly circular. Beyond this mean density, they become elongated and the radial alignment of the squirmers cannot clearly be defined.
Thus our claim is that by introducing the negative hydrodynamic pressure generated by squirmes swimming into the cluster at its rim, we are able to maintain mechanical equilibrium. Namely, by defining the total pressure p = p (s) + p (a) + p (h) with the hydrodynamic pressure included, this total pressure should be constant across the interface of gas and dense phase. Our observation that |∆p (h) | increases with cluster size and thus withφ , explains the different binodals in Fig. 5 and why we can collapse all of them on a single master curve with the help of the effective Péclet number from Eq. (6). For an increased |∆p (h) |, a lower ∆p (a) ∝ φ gas is required to maintain the interface and thus φ gas is reduced.
For phase-separating systems, interface curvature leads to a Laplace pressure and can result in a shift of the coexistence densities 61 . We point out that this effect cannot be responsible for the shifted binodals as well as the pressure jump in Fig. 5: bottom. First, the Laplace pressure changes signs as one goes from a droplet to a bubble configuration. In particular, the pressure jump in Fig. 5: top was determined for the "slab" configuration with nearly straight interfaces, where the Laplace pressure should vanish. Second, the coexistence densities of the dilute and dense phase in the middle row of Fig. 2: top do not change while coarsening. However, if Laplace pressure were important in our system, the densities would change since the curvature of the interfaces changes.
Until now we have only examined neutral squirmers. In the following we present binodal lines for different squirmer parameter β to demonstrate how pushers and pullers effect the phasecoexistence region of microswimmers.

Binodal Lines of Pushers and Pullers
The flow field around squirmers is strongly controlled by the dipole strength β used in our swimmer model. For β = 0, the flow field of a single squirmer in a bulk fluid is reminiscent of pushers (β < 0) and pullers (β > 0) 1,62 . It has a strong effect on how swimmers reorient due to nearby surfaces or other squirmers as demonstrated in Refs. 1,28,62 . The qualitative effect of increasing |β | has been shown in our earlier work 28 to largely suppress phase separation, requiring systems of increasingφ in order for stable clusters to form.
The binodal lines for different β are compared to the neutralswimmer case in Fig. 7 for pushers (top panel) and pullers (bottom panel). The overall trend is that for increasing |β | the coexistence regime becomes smaller. For pushers, one realizes very nicely how the binodals are mainly shifted up, thereby increasing the critical Péclet number. The same applies to pullers, however, the gas binodals are also deformed compared to β = 0. We observe that puller squirmers are more mobile in the dense phase compared to pushers. Therefore, the corresponding binodals are shifted to smaller φ .
For pushers and pullers in full three dimensions spontaneous polar order was observed for small but non-zero β 39 . This results in a net motion of swimmers along a common direction. We do not find such a net motion for small β in our simulations. This might be due to the fact that in our quasi-two-dimensional geometry the far fields of neutral and non-neutral squirmers both decay as 1/r 2 40,59 . Thus, squirmers with small β = 0 do not have a sufficiently different far field compared to neutral squirmers.

Melting of the Cluster Phase
In the region below the critical Péclet number in Fig. 3: bottom the phase is isotropic. On the other hand, at large effective Péclet numbers we observe hexagonal order in the dense phase, which hence has to undergo a melting transition if the Peclet number is reduced. It is not our intention to present a thorough discussion of this subtle phenomenon but to only give an indication for melting.
The dense branch of the binodal lines of neutral squirmers and of pushers (see Fig. 7: top) reveal a similar trend. At Pe ≈ 400 and below, the density of the dense phase decreases and stabilizes at a lower value around Pe ≈ 350. In particular, this feature is clearly seen for neutral squirmers and coincides with a loss of hexagonal order in the cluster. As in Ref. 28 we introduce the bond parameter |q 6 | 2 to measure local sixfold bond-orientational order around a squirmer. In Fig 8 we plot the probability distribution for |q 6 | 2 in the dense phase for a range of Pe . For Pe ≥ 402 a peak at |q 6 | 2 ≈ 1 dominates the distribution and indicates hexagonal order in the dense phase. The hexagonal domains are stable over time with only short perturbations when defects and vacancies move through the hexagonal lattice due to the random component of the squirmers' motions. These perturbations slightly reduce the time-averaged |q 6 | 2 and shift the peak position below 1.
Below Pe ≥ 402 the sharp peak in the distribution function disappears indicating a melting transition with a loss of hexagonal order. In the simulations we still observe stable clusters, which now form fluid droplets. During the melting transition we also observe an increase of the mobility of the squirmers by a factor of 5-10, when monitoring the mean-square displacement in the dense phase. The supplemental material † contains videos M5 and M6 showing the hexagonal and fluid dense phase, respectively. In video M6 we see that, although hexagonal patches exist, they are not stable over time. Finally, when lowering Pe towards the critical point, the fluid cluster "evaporates" as the system exits the phase-coexistence regime.

Summary and Conclusion
We have studied the phase separation of model hydrodynamic microswimmers called squirmers in quasi-two-dimensional confinement extending our previous investigations 28 . The full threedimensional hydrodynamics was simulated in order to accurately capture the squirmer interactions. This makes our simulations a realistic representation of real-world artificial and biological microswimmers. By employing a parallelized implementation of the MPCD algorithm on over 1008 CPU cores, we were able to simulate thousands of squirmers, allowing us to quantitatively determine the binodal lines for different mean densitiesφ and squirmer parameters β . The coarsening dynamics towards the phase-separated system shows a characteristic intermediate diffusive regime and a final rather ballistic regime due to compactification of the dense squirmer cluster. However, most strikingly, we found that the position of the binodal lines is strongly influenced byφ . This effect has not been observed for active Brownian particles, which only interact sterically, and thus is of real hydrodynamic origin. We were able to collapse the different binodals on a single master curve by introducing an effective Péclet number. Extending the mechanical pressure balance of Ref. 9,10 by a hydrodynamic pressure due to the squirmer flow field, the reason for the different binodals became apparent.
We also explored the dependence of the binodal line on β , a parameter which allows us to tune the swimmer type and their flow fields between pushers and pullers. We found that in all cases going to strong pushers or pullers suppresses phase separation by shifting the phase-coexistence region to larger Péclet numbers. The shape of the coexistence region does, however, depend on whether the squirmers are pushers or pullers.
Finally, we examined the structure of the dense phase as the squirmers' swimming speed, the Péclet number, is reduced. We find that fast squirmers in the dense phase form a hexagonal cluster, which melts with decreasing Péclet number into a fluid cluster. This is expected since close to the critical Péclet number the swimmer fluid is isotropic.

Acknowledgements
We would like to thank I. Aranson, R. Golestanian, G. Gompper, R. Kapral, K. Kroy, T. Speck, and J. Tailleur for stimulating discussions and the referees for their very helpful feedback. This project was funded under the DFG priority program SPP 1726 "Microswimmers -from Single Particle Motion to Collective Behaviour" (grant number STA352/11). Simulations were conducted at the "Norddeutscher Verbund für Hoch-und Höchstleistungsrechnen" (HLRN), project number bep00050.

Position-Resolved Time Averaging
In order to associate an average quantity to a position, one commonly divides up the system into a grid. We call each cell of the grid a "pixel". A spatially-resolved average can thus be found by simply averaging over the squirmers in each pixel. Since we have only a few thousand squirmers, each pixel would need to be relatively large in order to capture enough squirmers per pixel. This would prevent the accurate resolution of the phase interface.
We will use the local density to illustrate our method. First we partition the system into its Voronoi cells. Each Voronoi cell contains a single squirmer, and thus its local density is φ i = 1/A i where A i is the area of the i-th Voronoi cell. Our grid of pixels is then assigned a local density where r i, j is the position of the pixel with indices i and j. The sum is taken over those Voronoi cells overlapping with the square boundaries of the pixel. φ k is the local density of Voronoi cell k, and w k = a overlap /a pixel is a weight equal to the Voronoi cell's overlap area a overlap with the pixel normalized by the pixel area a pixel . An example of this is shown in Fig. 9.
All time averages are taken with respect to a specific spatial position. Due to the motion of the squirmers, the time averages of φ (r i, j ) eventually wash out the structure of the voronoi construction, leaving only those density inhomogeneities that are stationary in space. We use the same construction for calculating local pressure and |q 6 | 2 values. Figure 10 compares the binodals from Fig. 3 with those for active Brownian particles in quasi two dimensions. Brownian dynamics simulations where conducted for the same range of parameters and the same geometry as the full hydrodynamic simulations. Similarly to simulations conducted in two dimensions 11 for the range of Péclet numbers simulated, the coexistence densities for active Brownian particles did not vary much. We also note that the binodals for two different mean densitiesφ = 0.41 andφ = 0.70 lie on top of each other. For comparison, the binodals for the hydrodynamically interacting (HI) squirmers from Fig. 3   C v (t) Fig. 11 Self-diffusion coefficient of passive particles plotted versus their density φ confined between two walls with distance H = 8R/3. The full line is a fit to D T,2D (φ ) = D T,2D (φ = 0)/(1 + bφ ) with b = 2.027. Inset: The corresponding in-plane velocity-auto-correlation functions C v (t) = v(t) · v(0) /2k B T for the different densities in the main plot. The differences are rather subtle. Time is given in MPCD units. Fig. 12 Illustration of the swim pressure in the central Voronoi cell (central particle colored in red). The coordinate system is centred at the central particle. Its Voronoi cell has a volume V i . Nearby particles exert a swim pressure γ 0 v 0 ∑ jê j · r j /3V i "on" the Voronoi cell of the central particle, where j runs over all neighboring Voronoi cells.

Position-Resolved Swim and Steric Pressure
Before we explain, how we calculate the pressure values, we present in the inset of Fig. 11 the velocity-auto-correlation function for the in-plane velocity components of a passive particle, determined in MPCD simulations for different mean densities of passive particles. The self-diffusion coefficient D T,2D (φ ) = v(t) · v(0) dt/2k B T is then plotted in the main graph of Fig. 11 and fitted by D T,2D (φ ) = D T,2D (φ = 0)/(1+2.027φ ). Using the Einstein relation, we obtain a density-dependent friction coefficient, Now, we apply Eq. (8) to the individual Voronoi cells and its surround neighbours (cf. Fig. 12). Our system volume now is V i ≈ HA i , i.e., the height of the system H times the 2D-Voronoi cell area A i , and the neighbours are the immediate neighbours in terms of the Voronoi cell. On this level, the force generating the active motion F swim = γ 0 v j of neighbouring squirmers will act on the central Voronoi cell, i.e. inducing a pressure on its "walls". There is an additional subtlety here: the position vectors r j of the neighboring squirmers always have to start from the central swimmer, for which local pressure is calculated, otherwise shear terms would appear in the derivation of Eq. (8) 63 . Furthermore, the mean velocity also needs to be subtracted from the individual swimmer velocities. Thus we are performing calculations in the centre-of-mass frame of the neighbourhood of each squirmer. We reduce the amount of fluctuations by taking the time average for each pixel (cf. previous subsection on spatial averaging). This time average is longer than the decorrelation time D −1 R of the swimmers as required by Eq. (8).
Similarly, we use the Voronoi cells to also introduce a local steric pressure following Ref. 57 . The steric forces acting on one squirmer include forces from its neighbors and also from the bounding walls.