Multi-particle collision dynamics with a non-ideal equation of state. I

The method of multi-particle collision dynamics (MPCD) and its different implementations are commonly used in the field of soft matter physics to simulate fluid flow at the micron scale. Typically, the coarse-grained fluid particles are described by the equation of state of an ideal gas, and the fluid is rather compressible. This is in contrast to conventional fluids, which are incompressible for velocities much below the speed of sound, and can cause inhomogeneities in density. We propose an algorithm for MPCD with a modified collision rule that results in a non-ideal equation of state and a significantly decreased compressibility. It allows simulations at less computational costs compared to conventional MPCD algorithms. We derive analytic expressions for the equation of state and the corresponding compressibility as well as shear viscosity. They show overall very good agreement with simulations, where we determine the pressure by simulating a quiet bulk fluid and the shear viscosity by simulating a linear shear flow and a Poiseuille flow. © 2021 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0037934., s


I. INTRODUCTION
Since their introduction in 1999, 1 algorithms belonging to the method of multi-particle collision dynamics (MPCD) have become a standard tool to simulate fluid flows in the field of soft matter physics. [2][3][4] In particular, MPCD algorithms have commonly been used to model solvent dynamics in the context of microswimmers, [5][6][7][8][9][10][11][12][13][14][15][16][17][18] where we can only cite a few examples. Further studies address colloidal suspensions, [19][20][21][22][23] polymers, [24][25][26] blood cells, 27 the African trypanosome as the causative agent of the sleeping sickness, 28 and even fish schools. 29 Also, extensions to binary and ternary fluid mixtures, [30][31][32] liquid crystals, [33][34][35] and chemically reacting systems 36 exist. MPCD methods are particularly suited to simulate solvent flow on the microscopic scale because they solve the Navier-Stokes equations but also incorporate the omnipresent thermal fluctuations. 1,37 The particle-based strategy of MPCD makes the implementation of no-slip boundary conditions in complex geometries very straightforward. 38 Furthermore, the collision rules for the coarse-grained fluid particles are well suited for the implementation on parallel computer hardware 39,40 so that extensive simulations can also be performed on desktop computers with graphic cards.
Although MPCD methods are often used to simulate the dynamics of incompressible solvents, one has to be aware that the coarse-grained fluid particles follow the equation of state of an ideal gas. 1,32,41 Therefore, the fluid is rather compressible and has a low speed of sound cs. 21 This is tolerable for typical flow velocities well below cs. In contrast, in the presence of large pressure gradients, pronounced inhomogeneities in the fluid density can occur due to the high compressibility. For example, such a situation has recently been observed in strongly clustered microswimmers, 12 where the overlapping flow fields of many microswimmers are responsible for strong pressure gradients. While variations in fluid density are, in principal, necessary to generate pressure gradients, these variations need to be small to stay close to the limit of an incompressible fluid. Thus, the compressibility needs to be sufficiently small. For the MPCD fluid with its ideal-gas equation of state, this can be achieved by increasing the number n 0 of fluid particles per collision cell and thereby density. 12 However, such an approach causes an immense increase in the simulation time proportional to the square of the fluid density n 2 0 if the system size should be kept constant at an equal Péclet number.
In this paper, we follow a different strategy to decrease compressibility. All the MPCD algorithms consist of a sequence of collision and streaming steps. Here, we propose a new collision rule that results in a non-ideal equation of state for the MPCD fluid. Note that such non-ideal equations of state are required and have already been introduced in the context of simulating fluid mixtures within MPCD. 31,32 Thus, compressibility is reduced for the constant particle number n 0 , and the computational efficiency is enhanced compared to conventional MPCD algorithms, which need to employ a larger particle density. Our approach extends ideas of Tüzel, Ihle, and collaborators, who included geometric properties of hard-core particles in two dimensions into the collision rule to control momentum transport in the fluid. 42,43 This approach has also been extended to the simulation of fluid mixtures. 30 In contrast to conventional MPCD algorithms, where collisions take place at a fixed rate 1/Δt, collisions instead occur stochastically with a probability that depends on the local density and velocities. In the present work, we extent the approach of Refs. 42 and 43 to three dimensions and strongly modify the geometric rules of the collision so that they can be implemented in an existing MPCD code more easily. Furthermore, our new collision rule allows us to keep the typical canonical thermostat and also to take care of angular momentum conservation during collisions, which is particularly important for the simulation of colloids and active particles. 44 This article is structured as follows: In Sec. II, we introduce the extended MPCD method with its new collision rule including three possible collision probabilities. Then, we derive approximate analytic expressions for the equation of state and the associated compressibility in Sec. III as well as the shear viscosity in Sec. IV. For the shear viscosity, we consider both contributions that arise from the streaming and collision step of the extended MPCD method. In Sec. V, we compare these analytic expressions with the results from simulations and obtain very good agreement for the equation of state. In particular, we demonstrate the reduction of the compressibility for reasonable particle densities. We measure the shear viscosity by determining the collisional and streaming viscosities in a linear shear flow geometry. The total viscosity agrees very well with values determined from simulating a Poiseuille flow and also with the analytic expression above a density of ∼n 0 = 20. We close with conclusions and an outlook in Sec. VI.

II. ALGORITHM OF THE EXTENDED MPCD METHOD
Our method shares the common features typical for the group of MPCD algorithms. 1,3,45 Like all MPCD algorithms, it considers point-like particles that represent the fluid at a mesoscopic level of description. They perform a sequence of streaming and collision steps. Since the latter conserves linear momentum, the resulting hydrodynamic flow fields fulfill the Navier-Stokes equations. 1 While we perform the steaming step as in other MPCD algorithms, we alter the collision step as already mentioned in the Introduction. We now explain the extended MPCD method in more detail.
During the streaming step (i), the point particles with masses m 0 , positions xi(t), and velocities vi(t) move ballistically during time Δt, They collide with confining walls or moving objects such as model microswimmers called squirmers 4,8,9,[16][17][18] and thereby transfer both linear and angular momentum to these moving objects. By applying the so-called bounce-back rule, 1,38 the collisions either enforce the no-slip boundary condition at confining walls and passive colloids or the slip-velocity field, which are present at squirmer surfaces. For the collision step (ii), we suggest an alternative algorithm compared to the original SRD method 1 or the collision operator based on the Andersen thermostat. 4 As in all MPCD algorithms, the simulation volume is divided by a cubic lattice and the fluid particles are grouped into the cubic unit cells of linear size a 0 centered around ξ and with volume V ξ . Each cell then contains n ξ particles with the mean velocity v ξ and center-of-mass position x ξ . Additionally, each cell is divided into two halves A and B by a plane P x ξ ,n through the center-of-mass position x ξ 46 and with an orientation defined by the unit normal vectorn [see Fig. 1(a)]. For each cell, n is randomly drawn from a discrete set of 13 possible orientations at each collision step. By definition,n always points to region A. The number of particles on each side of the plane is denoted by nA and nB, respectively. Their mean velocitiesvA andvB, respectively, along the normal vectorn are given relative to v ξ .
The main idea of the new collision step is that the particles in region A and B only collide when they move toward each other. Then, they stochastically exchange a momentum m 0 δvi alongn both with particles in the same half A, B and also on the other side of the plane. The latter mechanism generates momentum flux across the randomly oriented plane and thereby contributes to pressure, which belongs to the isotropic part of the stress tensor.
The collision step can be summarized by where xi,c denotes the position vector of particle i relative to the center-of-mass position x ξ . As we explain below, the collision between the particles in region A and B occurs with a certain probability. To initiate a collision, the stochastic variable χ(Δu) is set to one; otherwise, it is zero. The term following χ(Δu) in the square brackets sets the normal velocity components of all particles i to the normal component of the center-of-mass velocity,n ⋅ v ξ . Then, new values for the relative velocity component δvi are assigned as explained below. They all add up to zero in order to preserve the total momentum. The second term in the curly brackets is added to conserve the angular momentum. Thus, the value L ξ = m 0 ∑ x i ∈V ξ xi × vi  before the collision is preserved. Here, I ξ is the moment-of-inertia tensor of all particles in the cell relative to the center of mass. We have already introduced the mean values of the normal velocity componentsvA andvB on either side of the collision plane. Then, collisions between the particles of region A and B occur, when, on average, they move toward each other, meaning that the relative velocity is positive. Furthermore, collisions between the particle clouds in A and B occur with the rate cΔunAnB, where c quantifies the scattering cross section. A similar term has been used for the collision rate of two clouds of hard-core particles or in chemical reactions of the second order 47 and can be motivated by the collision term in the Boltzmann equation. 48 Thus, the probability that a collision occurs or that the stochastic variable χ(Δu) is set to one becomes Here, Θ(Δu) is the Heaviside step function so that collisions only occur for Δu > 0. In the second line, assuming a sufficiently small c, we have introduced the exponential that guarantees pχ(Δu) ≤ 1.
Another possibility to fulfill this constraint using Eq. (4) is 30 We will explore also this form in Sec. V B when we calculate the pressure in the MPCD simulations. Finally, we introduce the changes δvi in the velocity component along the normaln. It consists of two contributions: δvi = δv t i + δv s i . The first term transfers momentum from region B of the cell to particles i in the region A and vice versa, Here, the first indices apply to particles i in region A that take over the momentum m 0vB from region B and the second indices apply to particles i in region B. The ratios nB/nA and nA/nB guarantee the overall momentum conservation, meaning the total momenta from regions A and B are just swapped. The second contribution, assigns each particle a random velocity δv MB i drawn from a Maxwell-Boltzmann distribution at temperature T, which serves as a thermostat for the fluid. We subtract the mean random velocity to preserve total momentum in both regions A and B, separately.
In particular, the introduction of the momentum transfer in Eq. (7) and the transfer rate Eq. (5) defines the equation of state. As shown in Sec. III, it contains a term proportional to n 2 ξ resembling a virial expansion and thus extends the ideal gas term.
As in other MPCD algorithms, immersed boundaries are represented by the so-called "ghost" particles during the collision step. 38 These are added to the collision cells to interact with the other fluid particles. In simulations with squirmers, the ghost particles are assigned the local velocity of the translating and rotating squirmers plus a random thermal velocity drawn from a Boltzmann distribution. Then, the changes in linear and angular momentum of the ghost particles following from step (ii) are assigned to the relevant squirmer, which ensures that linear and angular momentum are conserved. Finally, before performing each collision step, the lattice is randomly shifted to ensure Galilean invariance. 49

III. EQUATION OF STATE
To calculate the equation of state, we use the definition of pressure as the normal component of the momentum flux through an arbitrarily oriented plane. 50 In the extended MPCD method, both streaming (i) and collision step (ii) contribute to the pressure, During the streaming step (i), particles do not interact and simply transport momentum across a plane. This results in the ideal gas contribution Pstr = nk B T/a 3 0 , which we already know from the conventional MPCD methods. 32 To evaluate the contribution P coll from the collision step (ii), we consider the momentum flux across a plane with area a 2 0 that lies in a single collision cell. Without loss of generality, we choose the plane Pŷ ,y 0 perpendicular to theŷ axis at position y 0 and then average over all y 0 (see Fig. 2). During the collision step, momentum is transported from the region y < y 0 across the plane Pŷ ,y 0 into y > y 0 during time Δt. Thus, for the pressure as momentum transfer per area and time, we obtain Here, i is restricted to all particles above Pŷ ,y 0 and m 0ŷ ⋅ (v new i − vi) is the change in the normal momentum component of particle i during

FIG. 2.
To derive the equation of state, we consider the momentum transferred to the region above the plane Pŷ ,y 0 with unit normal vectorŷ and at position y = y 0 . Note that to evaluate sign(x ⋅n ) in Eq. (15), one has to distinguish particles that are located in the green region as part of region A relative to the collision plane with normaln and particles in the blue region as part of B.

ARTICLE
scitation.org/journal/jcp collision and given in Eq. (2). The average goes over all possible collisions, particle configurations, orientationsn of the collision planes, and positions y 0 . The term added to Eq. (2) in the second line to preserve angular momentum does not contribute to P coll since it vanishes when averaging over all possible collisions. Furthermore, since particle i is either in region A and B and we average over all particle velocities relative to v ξ with identical velocity distributions, we can ultimately replacen ⋅ (v ξ − vi) in Eq. (2) by the mean velocities −v A/B in region A or B of the collision cell. Note thatv A/B are given relative ton ⋅ v ξ . The choice of index A or B depends on the location of particle i. Thus, we can simplify Eq. (11) to The stochastic contribution δv s i of δvi given in Eq. (8) obeys a Gaussian distribution with zero mean and therefore vanishes on average. The remaining part δv t i given in Eq. (7) becomesv B/A using nA /B ≈ n ξ /2. Thus, with the definition of the collision velocity Δu in Eq. (3), we can finally The factor sign(xi ⋅n) comes in since the first index in v B/A −v A/B applies if particle i is in region A, while the second index refers to a particle i in B (see Fig. 2). Noting also that Δu and xi are independent stochastic variables, we can factorize the average in Eq. (12) and rewrite it as where Here, p(Δu) is the probability distribution for Δu and pχ(Δu) is the probability for a collision to take place, as introduced in Sec. II. We now calculate the two averages of Eq. (13). In the second average, we replace the conditional sum by a volume integral introducing the factor n ξ /a 3 0 Θ(y − y 0 ) and average over y 0 . We write the second average as n ξ αP, where we identify the purely geometrical factor It is the difference between the green and blue volume in Fig. 2 averaged over alln and y 0 and weighted by the projection ofn onŷ. The integrals can be calculated for each of the 13 normal vectorsn so that we obtain in total For the second average, we need the probability distribution for Δu =vB −vA. Since the components of the single-particle velocities are Gaussian distributed with variance k B T/m 0 , also Δu is Gaussian distributed with variance 4k B T/(m 0 n ξ ), as shown in Appendix A using nA = nB ≈ n ξ /2. Taking the collision probability from Eq. (4), we then have Thus, in total, we obtain from Eq. (13) for the pressure contribution of the collision step, which is quadratic in the particle density n ξ . Hence, up to second order in density, the full equation of state reads This gives a compressibility where the ideal gas contribution from the streaming step, β id = a 3 0 /(n ξ k B T), is diminished by the second-order contribution from the collision step. This means that the MPCD fluid is less compressible.
We add two comments. First, if we take for the collision probability pχ(Δu) the expression from Eq. (5), which we use as one option in the simulations, one can still evaluate ⟨χ(Δu)Δu⟩ and then obtain the pressure contribution from the collision step, We will use this form when comparing the pressure in the simulations to the analytic result. Second, in deriving Pc, we have always set nA = nB ≈ n ξ /2, thus neglecting fluctuations in the particle numbers in regions A and B. For sufficiently large particle numbers, these fluctuations are small. When we compare our analytic results to simulations, we obtain good agreement and the approximation seems to be reasonable.

IV. SHEAR VISCOSITY
To derive an expression for the dynamic shear viscosity η, we consider the linear shear flow v(y) =γyx (22) with constant shear rateγ (see Fig. 3). We also note that the nonvanishing component of the viscous stress tensor, σxy = η∂yvx = ηγ, describes the negative flux of the x component of momentum along the y direction. Similar to the derivation of the equation of state in Sec. III, the viscosity consists of two contributions from the collision and streaming step, respectively, While the derivation of the collisional viscosity η coll requires similar steps used in calculating the collisional contribution of pressure, of Chemical Physics the streaming viscosity ηstr needs special attention. We start with the derivation of η coll .

A. Collisional viscosity
Similar to our derivation of the pressure starting from Eq. (11), we consider the momentum transported during the collision step from the region y < y 0 across the plane Pŷ ,y 0 into y > y 0 during time Δt (see Fig. 3). However, for σxy, we need the transfer of the x component of momentum per area and time, and thus, When evaluating the term in the angular bracket using Eq. (2), we concentrate on the first line of Eq.
(2) and then at the end following Ref. 45 to include angular momentum conservation, which is the origin of the second line. Replacing δvi in Eq. (2) byv B/A as before in Sec. III, we obtain On average, the relative velocity Δu =vB −vA is equally distributed on the mean velocities of regions A and B so that we can use Δu/2 ≈vA = −vB. Furthermore, on average, v ξ can be replaced by vi/n ξ , 51 and we arrive at As explained in Sec. III, sign(xi ⋅n) is necessary to distinguish between particle i being in either region A or B. Now, we need to introduce the shear rateγ from the applied linear shear profile of Eq. (22). It comes in by setting vi equal to its deterministic partγyix and through the Gaussian distribution of Δu. In Appendix B, we show that the conditional distribution for Δu, given the collision vectorn and fixed position xi of particle i, is Gaussian with the mean value The first factor originates from the orientation of the collision plane with collision vectorn relative to the shearing directionx and the second factor from keeping particle i at fixed height yi. The quantity yA is the y coordinate of the center of mass of region A defined by the collision plane. For the different collision vectorsn, we give them in Table I. We need this conditional mean value μ i,n = ⟨Δu⟩ i,n when averaging over Δu since in Eq. (25), we also average over the position of particle i. Now, to evaluate the shear viscosity, it is sufficient to only consider the terms of σxy linear in the shear rateγ. As we demonstrate in the following, they result from either thermal fluctuations of Δu or the deterministic part of vi equal toγyix. Thermal fluctuations of vi can be neglected since they produce higher-order terms inγ. To perform the average over Δu in Eq. (25), we first evaluate the required averages using the conditional distribution p(Δu−μ i,n ) (see Appendix C). Since the shear-induced shift μ i,n is small compared to the width of the distribution, m 0 μ 2 i,n /k B T ≪ 1, we can always linearize in μ i,n ∝γ. First, for the mean conditional collision rate using Eq. (4) for pχ(Δu), we obtain Only the contribution of the zeroth order in μ i,n is required, since the last term in Eq. (25) already contributes the required term linear inγ by setting vi =γyix. For the second necessary mean value, we obtain up to the linear order in μ i,n ∝γ. One can show that the zerothorder term Ξ(n ξ , c) does not contribute to σxy, as it should be. The contribution vanishes for collision planes with yA = 0 or in combination with two collision vectors. Using Eqs. (27) and (28) in the expression (25) and only collecting all terms linear inγ, we arrive at Here, the remaining average ⟨⋯⟩ goes over xi, the offset y 0 of the plane Pŷ ,y 0 , and the collision vectorn.
As in the derivation of the equation of state, we replace the average over all particles and the conditional sum by a volume integral over n ξ /a 3 0 Θ(y − y 0 ) and also average over y 0 . With only the average over the collision vectorn remaining, we obtain where In Table I, we give the values α η,n and yA for all collision vectorsn. Averaging over all of them, we obtain for the collisional viscosity without taking into account angular momentum conservation in the collision rule of Eq. (2), For our choice of c = 1/100 and n ξ = 20, we obtain Ω ≈ 5/14 and Γ = Ω/2. If instead of Eq. (4), we use Eq. (5) to have a bounded collision probability pχ(Δu), we still can evaluate the averages of Eqs. (27) and (28) and expand intoγ. The resulting expressions for Ω and Γ are given in Appendix C. For c = 1/100 and n ξ = 20, we then obtain Ω ≈ 1/4 and Γ ≈ 9/65. So far, we did not consider the term due to angular momentum conservation in our collision rule (2) when evaluating σxy. We follow here Ref. 45 to take into account two additional terms. The essential contribution is the rotational motion of the particles in the collision cell induced by the vorticity of the shear flow, which generates the rotational velocity ω = ∇ × v/2 = −γ/2ẑ. The velocity ω × xi,c of particle i due to this rotational flow is removed during the random collision, and we have to add it to v new i − vi considered so far to preserve angular momentum. More precisely, our collision rule (2) only considers the component normal to the collision plane, and when we average over all collision vectors, we realize that only the x component of ω × xi,c is needed. Hence, in total, we need to add to the last term in Eq. (25) the normal velocity component (n⋅x)γyic/2. When averaging over all particle positions j ≠ i, we can set yi,c = yi − y ξ = yi(1 − 1/n ξ ) following a similar reasoning as in footnote 1. Thus, a careful inspection of Eq. (25) and the following steps show that we have to subtract half of the last term in Eq. (25). Ultimately, this replaces Γ in Eq. (32) by Γ/2.
A minor contribution comes from the random velocity changes δvj during collision. We mention it here since it gives a near perfect agreement with the simulation results we will present in Fig. 5(a). The random changes δvj add the angular momentum ∑ j xj,c ×n δvj to the cell, for which we have to subtract a term in the second line of Eq. (2) in order to restore angular momentum conservation. This velocity also has to be considered inx ⋅ (v new i − vi) when starting with Eq. (24) for σxy, and then, the relevant steps carefully have to be repeated. As before, we only need the z component of the angular momentum. Dividing by the moment of inertia for the relevant z direction, I ξ,zz /m 0 = ∑ j x 2 j + y 2 j , and taking the average over Δu, we obtain the mean angular velocity ⟨ω δv i z ⟩ Δu =x ⋅n ⟨χ(Δu) ∑ j y j,c Δu/2 sign(x j ⋅n) I ξ,zz /m 0 ⟩. Then, evaluating this average and introducing the mean moment of inertia ⟨I ξ ⟩zz/m 0 = a 2 0 (n ξ − 1)/6 gives ⟨ω δv i z ⟩ Δu = −12(x ⋅n) 2 y 2 Aγ Ω/a 2 0 , where we neglected correlations in the product yjΔu and used ⟨χ(Δu)Δu⟩ ≃ −2Ωx⋅nyAγ. Finally, from the angular velocity, one calculates the mean x component of the velocity correction, ⟨⟨ω δv i z ⟩ Δu yic⟩. Performing the remaining averages, one ultimately realizes that this changes the prefactor of Ω in Eq. (32) from 361/567 to 0.5034 ≈ 1/2. Together with the correction from the previous paragraph, we then obtain the final formula for the shear viscosity,

B. Streaming viscosity
To determine the streaming viscosity based on the linear shear flow of Eq. (22), we follow the work of Kikuchi et al. 50 They determined the shear stress component σxy from the momentum along the x direction transported through the plane y = 0 during the streaming time Δt. They showed that this results in the expression where vx and vy are velocity components of the fluid particles. The average is performed at the beginning of the streaming step. In the steady state, we can immediately use ⟨v 2 y ⟩ = k B T/m 0 . However, as we explain now, the velocity correlation ⟨vxvy⟩ changes when we cycle once through the streaming and collision step. However, in the steady state, it should be back to the value at the start of the cycle. Using this self-consistency condition, one can ultimately determine ⟨vxvy⟩ and therefore σxy.
First of all, if p(vx, vy) is the velocity distribution of the particles at the beginning of the streaming step, it will evolve toward The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp the distribution p(vx +γvyΔt, vy) at the end of the streaming step since particles in the shear flow acquire additional speed along the x axis when moving along the y direction. Based on this altered distribution, the velocity correlation at the end of the streaming step becomes 50 ⟨vxvy⟩ new str = ⟨vxvy⟩ −γΔt⟨v 2 y ⟩.
In other words, the value of ⟨vxvy⟩ decreases by a constant value during the streaming step. Both Eqs. (34) and (35) are common to all MPCD algorithms. 45 In a second step, ⟨vxvy⟩ new str is altered during the subsequent collision step. This depends on the detailed collision rule. As we demonstrate below and in Appendix D, the velocity correlations change by a constant factor during collision. Thus, ⟨vxvy⟩ new coll = (1 − b)⟨vxvy⟩ new str . Inserting Eq. (35) and using the self-consistency condition ⟨vxvy⟩ new coll = ⟨vxvy⟩ as explained above, we can solve We insert this result into the expression (34) for σxy and use ⟨v 2 y ⟩ = k B T/m 0 to finally arrive at Thus, after determining the factor 1 − b for our collision rule, we will have an expression for the streaming viscosity. In order to write ⟨vxvy⟩ new coll in a compact way, we abbreviate in the collision rule of Eq. (2) the term added to restore angular momentum conservation by A and use for the other term Bi =n[n ⋅ (v ξ − vi) + δvi]. Furthermore, right before the collision, the velocity correlation is ⟨vxvy⟩ new str so that we have ⟨vxvy⟩ new coll = ⟨vxvy⟩ new str + ⟨χ(Δu)[vi,xBi,y + Bi,xvi,y + vi,xAy + Axvi,y + AxAy + Bi,xBi,y + AxBi,y + AyBi,x]⟩.
Since the value of χ(Δu) is either 0 or 1, we have set χ(Δu) 2 = χ(Δu). Note that in δvi = δv s i + δv t i , we can drop δv s i since it is zero, on average, and also set δvt =v B/A using nA = nB ≈ n ξ /2 in Eq. (7). Hence, we will always use δvi =v B/A in the following. For the first term in Eq. (38), we demonstrate here how it is evaluated and refer to Appendix D for the evaluation of all the other terms. We obtain with δvi =v B/A ⟨χ(Δu)vi,xBi,y⟩ = ⟨χ(Δu)vi,xny[v B/A +n ⋅ (v ξ − vi)]⟩.
In the last line, we used again the molecular chaos assumption and neglected higher correlations for particle i.
from which we read off the factor b as the second term in the brackets. Using it in Eq. (37) together with ⟨χ(Δu)⟩ ≈ Γ(n ξ , c) from Eq. (27) and dividing by the shear rate, we obtain the streaming viscosity The sum of this equation and the collisional viscosity from Eq. (33) gives the complete shear viscosity in this new version of MPCD,

V. COMPARISON WITH SIMULATIONS
In this section, we compare the derived analytic expressions (21) for the pressure and (33) and (41) for the collision and streaming viscosities with values obtained from simulations. To calculate the collisional contribution to the pressure, we use with Eq. (11) the same formula with which we started the analytic calculations. Likewise, for the collisional and streaming viscosities, we set up the linear shear-flow profile v(y) =γyx and then explore Eqs. (24) and (34), respectively, to evaluate the viscosities from η = σxy/γ. Finally, to test our method in a realistic situation, we also simulate a Poiseuille flow profile and measure the total viscosity from the maximum flow velocity. We start with some computational details.

A. Computational details
To calculate the collisional contribution to the pressure equation of state, we perform MPCD simulations in a box with edge length L using periodic boundary conditions and the parameters introduced further below. The setups of the shear flow profile and the Poiseuille flow need more comments. All the simulations are performed with the bounded collision probability of Eq. (5). For the pressure, we will also show results for the alternative form of Eq. (6).

Linear shear flow profile
To generate a steady shear flow profile with constant shear rate ∂yvx =γ in a cubic simulation box with edge length L, we use the so-called Lees-Edwards boundary conditions. We introduce them shortly. 50,52 In the directions along thex andẑ axes, perpendicular to the shear gradient, regular periodic boundary conditions are applied. However, along the direction of the shear gradient, the boundary conditions are modified such that the periodic images of the system move with velocity ±Lγ. This means that a particle receives a The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp shift in position and velocity when crossing the boundaries along the y direction. If the particle crosses the lower boundary at time t and position (x, y = −L/2, z), it re-enters the system at the upper boundary at position ((x + Lγt) mod L , y = +L/2, z) with velocity (vx + Lγ, vy, vz). If it crosses the corresponding upper boundary at time t and position (x, y = + L/2, z), it re-enters the system at the lower boundary at position ((x−Lγt) mod L , y = −L/2, z) with velocity (vx − Lγ, vy, vz).

Poiseuille flow profile
To generate a Poiseuille flow profile, we do not introduce bounding walls but simulate a driven system with periodic boundaries in all three dimensions and two profiles with opposing flow directions along thex axis. 53 For this, we introduce a pressure difference Δp by acting with a constant body force on all the particles. Particles with positions yi < 0 experience the force −m 0 gx, and for particles with positions yi > 0, the force points in the opposite direction, m 0 gx. This setup with the box dimensions L × 2L × L produces two opposing Poiseuille flow profiles and thereby avoids the implementation of any solid boundaries. With the resulting pressure gradient Δp/L = m 0 n 0 g/a 3 0 , the viscosity then follows from the maximum flow velocity vmax = ΔpL 8η . 54

Parameters
For all simulations, we use the edge length L = 64a 0 , the collision parameter c = 1/100, and, in MPCD units, set k B T = 1 and mass m 0 = 1. For the Lees-Edwards simulations, the shear rate is chosen asγ = ∂yvx = 0.006 25 in units of the inverse MPCD time scale t 0 = a 0 √ m 0 /k B T. Each system is initialized by randomly distributing N = n 0 L 3 particles in the volume L 3 and by choosing their velocities from the Maxwell-Boltzmann distribution. For the Lees-Edwards simulations, a local offset for the mean velocity component along the x direction is chosen, ⟨vx⟩ =γy. To equilibrate the system at the beginning, we simulate it for 10 5 time steps Δt. Then, we sample Eqs. (11), (24), and (34) during a simulation time of 5 × 10 5 Δt. When simulating the Poiseuille flow, we average the flow profile over the same amount of time but use an increased equilibration time of 5 × 10 5 Δt to assure that the flow has reached its maximum velocity.
Our goal is to perform the MPCD simulations with defined values of the parameters, which we keep constant throughout the simulations. The collision parameter c introduced in Eq. (4) has to be sufficiently small so that we can explore the dependence on Δu, nA, and nB. It turns out that c = 1/100 and an average number of n 0 = 20 particles per cell is a suitable choice, which yields a collision rate of ⟨χ(Δu)⟩ ≈ 0.14.
Together with a time step of Δt = t 0 /200, our set of parameters is particularly interesting because it yields a total viscosity η = ηstr+η +A coll ≈ 16a 0 √ k B T/m 0 , which is commonly used for simulating microswimmers with MPCD. 8,9,12,17 Hence, in the following, we focus on densities between n 0 = 7 and n 0 = 35 and investigate how pressure and viscosities behave in this range of densities centered around n 0 = 20. In Fig. 4(b), we plot the corresponding compressibility as a function of n 0 relative to its ideal-gas value β id = 1/P id . Using P = P id [1 + f (n 0 )], the compressibility

B. Equation of state
can be directly related to the deviation of pressure from P id . The dashed and dotted lines in Fig. 4(b) represent the analytic results calculated from the formulas for pressure, while the derivative ∂f /∂n 0 The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp for the numerical results was determined with the standard Python toolchain. Relative to the ideal-gas value, compressibility is further reduced and, in particular, β also decreases with the decrease in Δt.
For example, at Δt = t 0 /200 and with n 0 = 20 as a reasonable density, compressibility is reduced to 0.4β id . Now, applying the bounded collision probability of Eq. (6), the compressibility is down to 0.3β id .
To obtain such a reduction with conventional MPCD methods and the ideal-gas pressure, one would need to increase the particle number per cell by a factor of three. Thus, the new collision rule with its non-ideal equation of state reduces the computational efforts.

C. Shear viscosity
We first discuss the collisional viscosity. Figure 5(a) shows the simulated collisional viscosity η +A coll in MPCD units m 0 /(a 0 t 0 ) as a function of the density n 0 for three values of Δt. The circle symbols show data points from simulations using the bounded collision rule Eq. (5). Over a wide range of densities, the values are in very good agreement with the analytical result of Eq. (33) shown as dashed lines. Hence, there is a quantitative agreement between simulations and theory. Similar to the pressure, the simulations confirm the scaling η +A coll ∝ 1/Δt. For the two larger time steps Δt = t 0 /100Δt = t 0 /200, we see a deviation at densities n 0 ≤ 10. We attribute this to the following reasons: first, our collision rule is not constructed for small numbers of particles, and second, to derive Eq. (33), we neglected fluctuations of the center of mass position, which also requires a higher number of particles.
We now continue with the streaming viscosity ηstr that we extract from the same simulations. Figure 5(b) shows ηstr in MPCD units m 0 /(a 0 t 0 ) as a function of the density n 0 and for different values of Δt. Again, the circle symbols show data points for the bounded collision rule from Eq. (5), while dotted lines refer to the analytic values given by Eq. (41). Although we observe an approximate quantitative agreement of the simulated values for ηstr with Eq. (41) for larger densities n 0 , there are clear differences. First, Eq. (41) predicts an increase in the streaming viscosity toward smaller n 0 , which then falls sharply to zero at n 0 = 0 (not shown). The simulated streaming viscosities only show a slight increase for Δt = t 0 /100; otherwise, they are roughly independent of n 0 . Second, while we do not reproduce the predicted scaling ηstr ∝Δt, we observe a clear increase in the streaming viscosity with Δt, and thus, the expected trend is reproduced qualitatively. As a main reason for the disagreement of the simulated viscosities with Eq. (41), we consider the approximation ⟨χ(Δu)vi,xvi,y⟩ ≈ ⟨χ(Δu)⟩⟨vi,xvi,y⟩ made during the derivation of Eq. (41). Nevertheless, factoring out the collision rate ⟨χ(Δu)⟩ in the previous expression provides a rough quantitative estimate of the streaming viscosity ηstr as demonstrated.

D. Poiseuille flow
To simulate the Poiseuille flow profiles, we used the time step Δt = t 0 /200. After averaging the velocity field over the time 5 × 10 4 Δt, the final flow profiles are generated by also averaging along the x and z directions. The resulting profiles vx as a function of the y position are shown in Fig. 6(a) for different densities n 0 . The two opposing profiles are clearly visible, and in both regions y < 0 and y > 0, we observe excellent agreement with the expected parabolic shape. The decrease in the flow velocity toward higher densities n 0 already indicates an increase in the total viscosity η with n 0 . This increase is even more pronounced since also the pressure difference Δp increases with n 0 because we always use the same body force per fluid particle.
In Fig. 6(b), we plot the total viscosity η as a function of the density n 0 . The blue circle symbols show data points obtained by extracting vmax from the Poiseuille flow profiles shown in Fig. 6(a) and using η = ΔpL/(8vmax). The red triangle symbols show the total viscosity η +A coll +ηstr consisting of the collisional and streaming viscosity, which we determined in Sec. V C from the simulated linear shear profile. The numerical data are compared to the analytical result of Eq. (42), which is shown as the dashed line.
First of all, the values for the viscosities, determined in the simulations by analyzing either momentum transfer in linear shear flow or the maximum flow speed of the Poiseuille profile, are in excellent agreement over the whole range of densities. In addition, we also observe a good agreement with the analytical expression of Eq. (42) for densities n 0 ≥ 10, where the collisional viscosity η +A coll dominates. For low densities n 0 < 10, the deviations occur due to ηstr as discussed before.

VI. CONCLUSIONS AND OUTLOOK
The new collision rule in our extended MPCD method provides the fluid with a non-ideal equation of state by introducing stochastic collisions between two particle clouds in the collision cell. In short, the collision frequency is quadratic in density and collisions only occur if the particle clouds move toward each other. In contrast to prior approaches, the extended MPCD method is designed for three dimensions, conserves angular momentum, and features a thermostat. The main goal of our method is to guarantee a low fluid compressibility for simulations in which significant pressure gradients occur. Since with the reduced compressibility, we can keep the particle number per collision cell at reasonably small values, our method requires significantly less simulation time compared to raising the fluid density in classical MPCD algorithms. 12 We provide an example in footnote. 55 At the same time, our method saves computer memory necessary to store MPCD particles so that we do not need to reduce the system size. We will explore this more in a planned second publication.
Based on the new collision rule, we have derived the equation of state and also demonstrated the impact of different collision probabilities. Indeed, in the regime where the collision probability is quadratic in density, we observe the nonlinear quadratic variation of pressure with density. For larger densities, where the collision frequency is bounded by the maximum value 1/Δt, the pressure again becomes linear in density, albeit at a higher value, which increases with 1/Δt. For typical values of Δt/t 0 = 1/200 and n 0 = 20 together with the most effective collision probability, compressibility is reduced by a factor three compared to the ideal-gas value at n 0 = 20. Overall, we find very good agreement with values obtained from simulations in the regime where our analytic expressions apply.
Moreover, for the shear viscosity, we have derived analytic expressions for the contributions of the collision and streaming step. For the collisional viscosity, we find very good agreement with the values obtained from simulating a linear shear flow and determining momentum transport, while for the streaming viscosity, the analytic expression only provides a rough estimate. However, for density values n 0 above 10, the collisional viscosity starts to dominate and we obtain a very good agreement with the simulated values. This is also demonstrated by simulating a Poiseuille flow and extracting viscosity from the maximum flow velocity.
In a planned second publication, we will use our extended MPCD method for selected flow problems to demonstrate its applicability. Furthermore, we intend to apply it to dense systems of microswimmers, where large pressure fields arise naturally. Preliminary simulations of such systems show that the extended MPCD method keeps the inhomogeneities in fluid density small. This will help us to obtain reliable insight into how hydrodynamic flow fields influence the collective dynamics of clustering and swarming microswimmers.

ACKNOWLEDGMENTS
We thank Josua Grawitter and Christian Schaaf for helpful discussions on the topic of the manuscript. We also acknowledge financial support from the Collaborative Research Center 910 funded by Deutsche Forschungsgemeinschaft.

APPENDIX A: GAUSSIAN DISTRIBUTION FOR Δu
For the derivation of the pressure, the fluid is considered at rest. Assuming molecular chaos, the velocities vj of the individual particles all obey the Maxwell-Boltzmann distribution. For the scalar product with the collision vectorn, this implies the probability density The probability distribution for Δu =vB −vA then follows from where δ(⋯) stands for the delta function. Eliminating it by integrating over one normal velocity component and then carefully performing the rest of the N − 1 integration finally gives the Gaussian distribution It has zero mean, ⟨Δu⟩ = 0, and variance where we used the approximation nA /B ≈ n ξ /2 in the last step. Thus, the sum over Gaussian distributed random numbers follows again a Gaussian distribution.

APPENDIX B: CONDITIONAL DISTRIBUTION OF Δu IN SHEAR FLOW
For deriving the collisional shear viscosity in Sec. IV A, we need the distribution for Δu under the condition that the collision vector n and the position xi of particle i are given. We again start with the single-particle velocity distributions. Relative to the applied shear flow, the velocities still follow the Maxwell-Boltzmann distribution, Following the same reasoning as in Appendix A, this implies that the conditional distribution for Δu is again Gaussian with the same variance as before: ⟨(Δu − ⟨Δu⟩) 2 ⟩ = 4k B T/(m 0 n ξ ). However, the conditional mean of Δu under the applied shear flow and for fixedn and xi is non-zero. Starting from the definition of Δu and averaging over all particle velocities and all positions besides particle i, one obtains where the subscripts i,n indicate the conditions that the particle i resides at xi and the collision vector takes the valuen. When we introduce the center-of-mass in the respective regions ⟨yj⟩ = yA /B for xj ∈ V A/B and use yA = −yB and nA /B ≈ n ξ /2, we can ultimately write the conditional mean as Based on the conditional distribution p(Δu − μ i,n ), we can now calculate the required mean values ⟨χ(Δu)⟩ and ⟨χ(Δu)Δu⟩ in shear flow.

APPENDIX C: MEAN VALUES ⟨χ(Δu )⟩ AND ⟨χ(Δu )Δu ⟩ IN SHEAR FLOW
We start with the unbounded form of the collision rate pχ(Δu) in Eq. (4) and find for the mean collision rate where in the last line, we show the relevant zeroth-order term after expansion in μ i,n . For the second mean value, we obtain and, after expanding to first order in μ i,n , For the bounded form of the collision rate pχ(Δu) in Eq. (5), we can also calculate the mean values. The mean collision rate becomes where the last line shows the relevant zeroth-order term after expansion in μ i,n . The second mean value becomes and after expanding to first order in μ i,n , one has
We begin by applying some transformations on the abbreviations Bi and A that we introduced to write Eq. (38) in a compact way. First, we note that we may drop the stochastic part δv s i of δvi = δv t i + δv s i because it averages to zero. Furthermore, we replace δv t i =v A/B using nA = nB ≈ n ξ /2. With the definitions of v ξ andv A/B , the quantity Bi reads and using that, we are free to add the constant velocity v ξ inside the round brackets.
With these simplifications, we now begin considering the next most simple term ⟨χ(Δu)Bi,xBi,y⟩ of the new velocity correlation ⟨vi,xvi,y⟩ new coll .
In the next term ⟨χ(Δu)vi,xAy⟩, we note that the quantity A also depends on the positions of the particles so that these must be included in the average. We may perform this average over the particle positions separately on A to obtain where xi,c denotes any of the components of xi,c. In the last line, we assumed that the contribution of the single particle i is low so that we can average ⟨I ξ ⟩ = (n ξ − 1)m 0 a 2 0 I/6 separately. Furthermore, we used that ⟨x 2 i,c ⟩ = a 2 0 /12(1 − 1/n ξ ) for any of the components of the position and that different components of the position xi,c are uncorrelated.
The last term to calculate is of which the last is zero because ⟨xj,xx k,y ⟩ = 0 for all j, k. We continue with the averages over the positions inside the sums. In the first summand of Eq. (D6), we recognize that ⟨xj,c,αx k,c,α ⟩ = ⟨x 2 j,c,α ⟩δ kj , and in the second one, we may rewrite ⟨xi,c,αx i,c,β ⟩ = ⟨x 2 i,c,α ⟩δ αβ and ⟨xj,c,αx k,c,y ⟩ = ⟨x 2 j,c,α ⟩δαyδ kj . This follows from the usual molecular chaos assumption that different particles are uncorrelated and the assumption that the components of a position are also uncorrelated. Performing the sums with these replacements, we obtain In the last line, we have furthermore neglected the term 9 20(n ξ −1) 2 , which is small compared to 1 4n ξ .

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.