Mesoscopic simulations of active nematics

Coarse-grained, mesoscale simulations are invaluable for studying soft condensed matter because of their ability to model systems in which a background solvent plays a substantial role but is not the primary interest. Such methods generally model passive solvents; however, far-from-equilibrium systems may also be composed of complex solutes suspended in an active fluid. Yet, few coarse-grained simulation methods exist to model an active medium. We introduce an algorithm to simulate active nematics, which builds on multiparticle collision dynamics (MPCD) for passive fluctuating nematohydrodynamics by introducing dipolar activity in the local collision operator. Active nematic MPCD (AN-MPCD) simulations not only exhibit the key characteristics of active nematic turbulence but, as a particle-based algorithm, also reproduce crucial attributes of active particle models. Thus, mesoscopic AN-MPCD is an approach that bridges microscopic and continuum descriptions, allowing simulations of composite active-passive systems.


I. INTRODUCTION
While fish shoals, bird flocks and insect swarms [1,2] are magnificent macroscopic examples of active systems that captivate onlookers with their collective behaviours, the majority of biophysical research focuses on active systems composed of microscopic agents.Microscopic, motile particles locally convert free energy from their surroundings into mechanical work [3][4][5] and collective dynamics [6].These active particles are subject to relatively large stochastic fluctuations that play a significant role in their dynamics.Prominent examples of active stochastic systems include suspensions of swimming microbes [7], bacterial colonies [8], tissue monolayers [9][10][11], and mixtures of cytoskeletal filaments and motor proteins [12].
In addition to possessing microscopic activity and stochastic dynamics, each of these examples displays orientational ordering [13][14][15].While shape-anisotropy is not a strict prerequisite of active self-propulsion, an innate direction of self-propulsion typically is.Furthermore, directionality regularly materializes as apolar orientation, even when the microscopic agents possess polar motility [16].Despite the fact that swimming bacteria [14], kinesin motor proteins marching along microtubules [12] and most other constituent self-propelled particles [5] possess distinctly polar behaviour, the hydrodynamic limit of suspensions of many such particles is nematic in nature.This is because the interactionsincluding dipolar forces enacted on the surrounding fluid medium [17] -are principally apolar [17,18].Thus, active nematics have proven a fruitful model for studying intrinsically out-of-equilibrium materials.
Alongside biophysical experiments, numerical simulations of nematic systems have been essential in developing physical understanding of activity's consequences for living materials.Different studies have attempted to model the microscopic details to a greater or lesser extent.For example, simulations of simple self-propelled rods in the dry limit [16] have found many of the same properties as more detailed simulations of substrate-crawling bacilliforms [19].Similarly, self-propelled rods in the wet limit [20,21] agree with simulations of swimming bacteria [22,23].While the microscopic details of active agents can vary significantly often even simple models can reproduce the essential emergent behaviours [5].
Simple microscopic models of stochastic, self-propelled particles (such as nematic variations on the Vicsek model [6,18,24,25]) have been particularly well-studied and can be coarse-grained into kinetic theories, which in turn lead to specific hydrodynamic equations of motion [26].On the other hand, symmetry considerations allow one to write generalised hydrodynamic equations for active fluids in the continuum limit [2,3,27,28].Unlike the microscopic models, these often omit stochastic effects within the fluid.Substantial work has explored the physics of bulk active fluids, such as steady-state creation and annihilation of topological defects [29][30][31] and the emergence of low-Reynolds number active turbulence [32][33][34][35][36].The emergence of complex spontaneous flows and topological singularities from field equations suggests that confining active nematics could produce complex and dynamic self-actualized dissipative structures.However, studies of confined active nematics have mostly been limited to simple, fixed geometries [37][38][39][40][41][42][43][44] and more complicated microfluidic geometries have rarely been considered.
Similarly, passive solutes suspended in active fluids have not received extensive consideration [28,45].Driven particles, such as colloids and disks within active nematics, have exhibited exciting characteristics, such as effective negative viscosity [46] and higher order defects [47].Experimentally, particles embedded in dense suspensions of swimming bacteria have exhibited anomalous diffusion [28,[48][49][50][51].While polymers suspended in intrinsically out-of-equilibrium athermal baths have been studied [52][53][54][55][56], hydrodynamic interactions mediated through the active medium are typically neglected [57,58].This highlights the need for coarsegrained simulation techniques that are capable of simu- lating both active fluids in complex geometries and suspended solutes possessing complex shape or internal degrees of freedom.
Here, we present a novel mesoscopic particle-based algorithm for stochastically simulating active, wet, compressible, nematic fluids.

II. ACTIVE NEMATIC MPCD ALGORITHM
Continuum active nematohydrodynamics are described by three transport equations for mass, momentum and orientational order [3].The method introduced here to simulate these equations of motion extends the multi-particle collision dynamics (MPCD) algorithm to active nematohydrodynamics.The active collision operator locally injects energy, while conserving momentum.This section first introduces the conceptual basis of the passive MPCD framework, summarises the collision operators employed and finally extends these to an activenematic collision operator.

A. Multi-particle collision dynamics
Multi-particle collision dynamics algorithms discretise continuous hydrodynamic field into N point-particles (labelled i).Each MPCD particle of mass m i , position r i and velocity v i is said to stream ballistically for a time δt to a new position before undergoing a multi-particle stochastic collision event [59].Though real molecules that constitute a fluid interact with one another via microscopically-specific pair potentials, the details of these molecular interactions are commonly inconsequential to the resulting continuum equations of motion in the isotrop hydrodynamic limit.The MPCD method leverages this reality to recovered hydrodynamic equations on long-time and length-scales using an artificial and mesoscopic multiparticle collision operations, rather than specific intermolecular potentials between pairs of molecules [60,61].MPCD collision events occur within lattice-based cells (labelled c) defined by a size a.In each cell, the instantaneous population N c of particles stochastically exchange properties through collision operators that respect the relevant conservation laws.The collision operator governs the fluid transport coefficients [60].
To conserve mass and reproduced the continuity equation ∂ t ρ = −∇ • (ρv cm c ), collision operators must simply leave the number of MPCD particles unchanged between collisions events.Like other particle-based hydrodynamic solvers, MPCD fluids are not strictly incompressible [62].
To conserve momentum, the cell's net momentum must be unchanged by collisions, which amounts to an unchanged centre of mass velocity v cm c = v c if all the particles have the same mass m i = m [59].The average • c is over the particles within cell c.The average velocity v cm c (r c ; t) is interpreted as the hydrodynamic velocity field at the position of that cell r c .Similarly, constraints to the stochastic exchange of particle velocities lead to conservation of energy and angular momentum, and their associated hydrodynamic fields.

B. Angular-Momentum Conserving Andersen MPCD
The continuous velocity field of an active nematic is given by a generalised Navier-Stokes equation in which the total hydrodynamic stress is the sum of a viscous stress, the orientational elastic stresses of nematic liquid crystals and an active stress [3].To simulate the stresses MPCD velocities are updated according to where the collision operator Ξ i,c is a non-physical exchange of momenta within each cell c.The collision operator is designed to be stochastic, while constrained to conserve net momentum (and kinetic energy, in passive fluids).The Andersen collision operator [63,64] is a common passive, isotropic collision operator with the form Here, ξ i is a random velocity drawn from the Maxwell-Boltzmann distribution for thermal energy k B T and ξ j c is the cell average.The moment of inertia is Nc j m j r 2 j 1 − r j r j for point-particles in cell c relative to their centre of mass r cm c , where r i = r i − r cm c .The third term in the collision operator corrects any spurious angular momentum δL c = δL vel + δL ori introduced by the collision operator δL vel = Nc j r j × v j − ξ j or liquid crystalline backflow δL ori from the dynamics of the director field.
The conceptual basis of MPCD can be extended from momentum collision operators that reproduce the hydrodynamic velocity field to collision operators for other hydrodynamic fields.In particular, nematic liquid crystals can be simulated via a collision operator to exchange particle orientations [65].

C. Nematic MPCD
In nematic liquid crystals, the tensor order parameter Q c measures the extent and direction of orientational order through its largest eigenvalue S c and associated eigenvector n c .The evolution of Q c (r c ; t) = du j u j − 1 c / (d − 1) can be simulated using passive nematic MPCD by assigning an orientation u i to each point-particle.The identity matrix is denoted 1.In this work we focus on an approach where the orientation is updated through a stochastic nematic multi-particle orientation collision operator [65] based on the local equilibrium distribution for the orientation field The noise η i is drawn from the Maier-Saupe distribu- , where the width of the distribution about n c (t) is controlled by a mean-field interaction potential, with an interaction constant U and inverse thermal energy β = 1/k B T .The interaction constant βU determines whether the fluid is in the isotropic or nematic phase, as well as the Frank elastic coefficient K [65,66].Nematic MPCD reproduces the one-constant approximation [65].
The orientations u i are coupled to the velocity field via a bare tumbling parameter λ and heuristic hydrodynamic susceptibility χ.Backflow coupling is introduced through a viscous rotation coefficient γ R [65].The backflow coupling is mediated through angular momentum transfer from changes to orientation back into linear momentum through δL ori = γ R Nc j u j × uj .Alternative hybrid finite difference/MPCD approaches have been proposed to apply MPCD to nematic liquid crystals [67,68].Nematic MPCD has been employed to study nematohydrodynamic fluctuations and correlations [69,70], electroconvection [71], defects around nanocolloids [66,72,73] and living liquid crystals [74].We now extend the nematic Andersen-thermostatted collision operator to simulate wet active nematics.

D. Active Nematic MPCD
To make MPCD active, an intrinsically out-ofequilibrium term must be included in the collision operation.We seek an algorithm for momentum-conserving (wet) active fluids [3]; therefore, we propose a collision operator that injects energy but not momentum.The form of the introduced active stress should correspond to a force dipole density that can be extensile or contractile [3].As in continuum models of active nematics, this algorithm employs a form for which the local active stress is proportional to the nematic order tensor ∼ Q c [75] or, equivalently, a form for which the force dipole co-aligns with the local director n c .
To account for these considerations, the collision operator is a linear combination of passive and active contributions The passive contribution Ξ 0 i,c is given by Eq. 3. The active nematic MPCD (AN-MPCD) collision operator is composed of two terms: (i) individual impulses (per unit mass) (α c δt/m i )κ i n c on each particle i and (ii) a term to ensure local conservation of momentum −(α c δt/ m c ) κ j c n c .Thus, as in the Andersen collision operator (Eq.3), any residual impulse is removed.
The activity term (α c δt/m i ) κ i n c is composed of four factors (i) δt, (ii) m −1 i , (iii) α c , and (iv) κ i n c = ±n c .(i) The factor of δt ensures activity α c is given per unit MPCD time; rather than per unit iteration.
(ii) The factor of m −1 i gives α c units of force.
(iii) The factor α c (t) represents the local active dipole strength.The cellular activity α c is found by summing over the individual activity value α j of each MPCD particle within a cell.
This represents a local activity that is directly proportional to the local density of active agents within a given MPCD cell.Positive particle activity α j corresponds to extensile active nematics; whereas negative values correspond to contractile activity.
Here, all particles have the same activity α j = α.
(iv) The factor κ i (r i , r cm c , n c ) n c gives the direction of the active force acting on particle i.In an active nematic, the activity is dipolar and is parallel to the local nematic director n c (t).Whether the impulse on each individual particle is parallel or antiparallel to the cell director is set through κ i (r i , r cm c , n c ) = ±1.This parallel/anti-parallel coefficient evenly splits the particles within cell c into those that are driven "forward" (κ i = +1) and those kicked "backward" (κ i = −1).To do this, the collision operator considers the plane passing through the center of mass r cm c with surface normal n c .Generally, half of the MPCD particles within cell c are on either side of this plane.The parallel/anti-parallel coefficient κ i = +1 for particles above the plane and κ i = −1 for particles below the plane.This algorithm is visually illustrated in Fig. 1.
Including the active collision operator (Eq.7) in the nematic MPCD algorithm results in a wet, particle-based, mesoscale, active nematic simulation method.

E. Simulation Units and Parameters
Distances are reported in units of MPCD cell size a ≡ 1 and energy in units of thermal energy k B T ≡ 1.This work considers only a single species, such that all particles have the same mass m i = m ≡ 1 and activity The MPCD time scale is thus τ = a m/k B T ≡ 1. Simulation parameters are expressed in these MPCD simulation units.Simulations are run in 2D in systems of size sys × sys with periodic boundary conditions and sys = 200, except when explicitly stated otherwise.Cartesian axes are denoted x and ŷ.The system-averaged fluid density is represented by the average particle density per cell, N c = 20.The time step is set to δt = 0.01 and, unless otherwise stated, simulations utilize a warmup phase of 10 5 time steps and a data production phase of 5 × 10 5 time steps.MPCD particles are initialized randomly within the control volume with random velocities, and orientations are initialised uniformly along the y-axis.The simulations are performed in the nematic phase by setting U = 10, with a rotational friction γ R = 0.01, and hydrodynamic susceptibility χ = 0.5.The tumbling parameters λ = 2 and λ = 1/2 are simulated as representative of shearaligning and flow-tumbling regimes respectively.The shear-aligning regime (λ = 2) is assumed, except where explicitly stated.The activity |α| < 1 is varied across four orders of magnitude to identify behavioural regimes of the algorithm.The activity is taken to be extensile (α > 0), unless otherwise stated.There exists a minimum activity α eq below which the Andersen thermostat can absorb the active energy injection, yielding effective equilibrium behaviour.Active turbulence, with its associated steady-state population of defects, arises at some greater activity α turb ≥ α eq and exists for a finite activ- ity regime α turb ≤ α ≤ α † , where α † denotes the end of the hydrodynamic active turbulence scaling regime.

A. Characteristic Activity Scales
Active nematic MPCD simulations are expected to possess a number of length scales.The algorithm itself has two; the system size sys and cell size a ≡ 1. Addi-tionally, an incompressible active nematic fluid has two: (i) the passive nematic persistence length which is proportional to the defect core size in the nematic phase; and (ii) the active length scale in fully developed mesoscale turbulence, which arises when the nematic elastic stress balances the active stress From this idealized dimensional analysis, the active length scales as α ∼ α µ with the expected power law µ = −1/2 in the regime of fully developed active turbulence.However, it has been demonstrated that µ can saturate if activity is not sufficiently large and the nematic structure spans a substantial fraction of the system size [76].Dimensional analysis also reveals an expected characteristic active velocity scale which arises when active stress is balanced by the viscous stress, characterized by the viscosity η.This is comparable to the speed of self-propelled +1/2 defects [29].Since Eq. 9 includes both a direct factor of activity and an indirect factor through the active length scale (Eq.8), the velocity scales as v α ∼ α γ , where ideally γ = 1/2.A mesoscale algorithm for simulating the hydrodynamic limit of active nematics should be consistent with µ ≈ −1/2 and γ ≈ 1/2 in the regime of simulating fully developed turbulence (i.e.α turb ≤ α ≤ α † ).
One method to quantify coherent flow structure within the AN-MPCD fluid is to analyze the spatial autocorrelation functions where R is the separation distance between two points in the fluid and r is the radial unit vector.The average is over all points in space and time.The unspecified field x (r 0 ; t) in Eq. 10 might be the director or velocity fields, x ∈ {n, v}.Since this work focuses on the hydrodynamics that result from activity, the microscopic scale R ≤ R kBT within which thermal effects dominate correlation functions are removed from the analysis, and the remaining correlation function is renormalised such that C xx (R = 0) = 1.Hydrodynamic correlations tend to initially decay and possess an anticorrelation well [33,77].To measure of active length scales from autocorrelation functions, the decorrelation length is determined by fitting an exponential decay to the small-R hydrodynamic region.The fit is preformed for the range Enstrophy spectra are employed when studying the spatial structure of traditional turbulence, and comparisons have been made to mesoscale active turbulence [78,79].For example, it has been found that Komogorov's universal −5/3 scaling exponent for inertial turbulence does not hold [79].The enstrophy Ω = |ω • ω| represents the magnitude of the vorticity ω = ∇ × v.In this work, the enstrophy spectra are computed via the radial Fourier transform of the vorticity correlation function where J 0 is the Bessel function of the first kind with order zero.

B. Density Analysis
A common feature of active particle-based models is their tendency to exhibit giant variations in local density [5,18,80,81].Nematic adaptations of the Vicsek model exhibit such giant number fluctuations [16,18,26].Furthermore, these have been observed in experimental and theoretical work on active-nematics [75,82].As a particle-based algorithm for active nematics, AN-MPCD is expected to exhibit significant number fluctuations.
In the typical case of equilibrium systems, density distributions are Gaussian due to the Central Limit Theorem [82].However, continuous energy injection can lead to non-Gaussian distributions.Hence, a measure of a distribution's Gaussianity can be revealing.A non-Gaussianity measure (NGM) can be defined to be where d = 2 is the dimension and ∆r k is the k'th moment of the distribution [19].When a distribution is purely Gaussian χ NGM = 0.However, χ NGM > 0 occurs when the the tails of the distribution stretch, and χ NGM < 0 occurs when the tails contract relative to normal.Physically, χ NGM 0 holds for density distributions of systems in equilibrium, while active-particle models exhibit non-Gaussinity.
To quantify the fluctuations in density, consider how the standard deviation of density, σ Nc , scales with the mean number density N c [82].In equilibrium, one expects with the scaling ν = 1/2 given by the Central Limit Theorem.However, in systems of active particles, ν > 1/2 is typical.As the scaling approaches ν 1, an active system is said to exhibit giant number fluctuations.These fluctuations in density are predicted by theory [75], experiments [82] and observed in simulations [83].The simulation domain is segmented into sub-domains in order to quantify these.By averaging over these sub-domains and computing the mean and standard deviations, the scaling ν is obtained.

IV. RESULTS
We now present the results of AN-MPCD simulations by first demonstrating that a hydrodynamic instability causes a proliferation of topological defects.The dependence of the emergent flows to the activity α are systematically explored by considering the structure of the flow fields.The results show that AN-MPCD possesses an operational active-turbulence-scaling regime, in which fully developed active turbulence exists with a characteristic emergent length scale α that scales with activity.The effects of activity on the local density are characterized and it is found that AN-MPCD exhibits giant-number fluctuations at sufficiently high activity.

A. Hydrodynamic Instability and Defect Unbinding
Experiments [84][85][86] and continuum simulations [87,88] of active nematic systems exhibit activity-induced instabilities, which generate inhomogeneities in the nematic order.This leads to a steady-state population of topological defects and active turbulence [30].The AN-MPCD algorithm reproduces this two-stage development of active turbulence through hydrodynamic instability and defect pair-creation.This processes is depicted by simulation snapshots in Fig. 2 from Movie 1.The simulation is initialised in a fully ordered state with all nematogen orientations parallel to the ŷ-axis, i.e. u i = ±ŷ ∀ i (Fig. 2a) and the warmup period is eschewed.At early times the hydrodynamic instability manifests as kink walls -narrow, initially parallel, well-spaced, alternating regions of high bend deformation (Fig. 2b).For extensile activity, fluctuations in bend produce active forces proportional to ∼ α∇ • Q, which drive spontaneous flows along the kink walls.However, these self-driven streams exacerbate the inhomogeneities through the back flow, which in turn in-creases the active forces causing the instability [89][90][91].
For sufficiently high activity, this hydrodynamic-bend instability is responsible for the formation of topological defect pairs (Fig. 2c).Beyond a certain point, the total distortion free energy cost along a kink wall is higher than the cost of two oppositely charged defects and a pair is spontaneously formed [30].Although there is an elastic restoring force, the +1/2 defect is a self-propelled quasi-particle and it moves along the kink wall, becoming unbound from the non-motile −1/2 defect [30,92] (Fig. 2d).
A pair-creation event is illustrated in detail in the top row of Fig. 3. Along a kink wall, a pair of ±1/2 defects are formed, conserving the net topological charge.The motility of +1/2 defects causes the pair to unbind.As the motile +1/2 defect travels along the kink wall, it alleviates the deformation free energy, leaving the region of uniform order in its trail (Fig. 3; top row) but generating local vorticity, which in turn perturbs nearby kink walls, disordering the orientation of future defect pair unbinding events giving rise to active turbulence (see Fig. 2c-d).Though activity generates defects, they come with an inherent free energy cost due to the deformation they induce in the director field, and as such they annihilate to minimise the deformation cost (Fig. 3; bottom row).As a +1/2 defect approaches a −1/2 defect, they annihilate and ordering the local region.Continual creation and annihilation leads to a steady-state defect population.

B. Steady-State Defect Population
To assess the degree of nematic ordering within the steady-state AN-MPCD fluid, the spatial autocorrelation of the director field C nn (R) is considered.For the lowest activities considered (α < ∼ α eq = 10 −3 ), the correlation remains high even as R → sys (Fig. 4a).In these cases, activities less than α eq are negligible and the injected energy can be absorbed by the thermostat so that the system remains globally ordered.This is seen in Movie 2 where the nematic director field remain oriented primarily in the ŷ-axis, with slight thermal fluctuations.For the low activity regime (α eq < ∼ α < ∼ α turb where α turb 10 −2 ), the director field decorrelates somewhat but does maintain long-range correlations/anti-correlations, which represent the persistent existence of kink walls with a characteristic separation length (Movie 3).Only once α > ∼ α turb do the correlation functions fully decorrelate to lim R→∞ C nn (R) → 2/3 and dip below the longrange value, as observed in continuum simulations of fully formed active turbulence (Movie 4) [33].The activity α turb = 10 −2 represents the lower threshold between fully developed active turbulence and kink walls for this system size.In the active turbulence regime, α turb < ∼ α, we find that the correlations belong to the same class of functions by rescaling the separation by the nematic decorrelation length (Fig. 4b).Rescaling collapses the curves for α > ∼ α turb .Consistent with the measured correlation functions, the number of defects is found to be non-zero only for sufficiently strong activity α > ∼ α turb , where the minimum activity to observe turbulence is α turb = 2.5 × 10 −2 for sys = 200 (Fig. 5).Above α turb the mean defect density ρ d rises rapidly and the root mean squared separation between ±1/2 defects, d = ρ −1/2 d , decreases accordingly (Fig. 5, inset).The scaling d ∼ α µ with µ = −1/2 is expected for a range above α turb (Eq.8) and the simulation results are found to match the expected scaling within the fully formed turbulence region α turb ≤ α ≤ α † (Fig. 5 inset).Thus, continual creation and annihilation leads to a steady-state defect population in AN-MPCD characterized by the competition between activity and nematic elasticity (Eq.8).
To assess system size effects, larger and smaller systems sys ∈ {400, 200, 100, 50, 25} for both tumbling (λ = 1/2 in Fig. 6a) and shear-aligning (λ = 2 in Fig. 6b) behaviours are simulated.In both cases, the system size shifts the activity at which defect pairs unbind and the range of activities for which the active-turbulence scaling law holds is extended for larger system sizes.

C. Spontaneous Flow
The previous section demonstrates that the AN-MPCD algorithm reproduces the hydrodynamic instability, topological defect propagation, the self-motility of +1/2 defects and the expected scaling of the characteristic length scale as measure by the mean separation between defects.This section explores the spontaneous flows and confirms that AN-MPCD reproduces the properties of fully developed active turbulence across a range of activities above α turb .
The magnitude of the active flows is quantified by measuring the steady-state, spatiotemporally averaged, root mean squared velocity v av (Fig. 7a).Since AN-MPCD is a thermalized method, lim α→0 v av = 0 and a smallbut-non-zero value of v av is measured below α eq = 10 −3 (Movie 5 and Fig. 7a).Above α eq the average speed rises rapidly -more rapidly than predicted by Eq. 9 (Fig. 7b).This is consistent with Fig. 4: Between α eq and α turb , the activity is sufficient to generate kink walls (Movie 3) and spontaneous flow (Movie 6) but insufficient to produce defects and active turbulence.Thus, the scale of the characteristic speed rises faster than when active turbulence is fully formed.Active turbulence is fully formed at α turb 10 −2 (Movie 7) and the scaling v av ∼ α γ is seen to be in agreement with the expected value of γ = 1/2 (Eq.9), shown as a dashed line in Fig. 7b, for the region α turb < ∼ α.As activity increases, v av does too, representing the broadening of the distribution of velocities.However, this is not the only effect; the probability density function also changes shape (Fig. 8a).For the lowest activities (α < ∼ α eq ), the distributions of the x-and ŷ-components are Gaussian.For larger activities (α > ∼ α eq ), the distribution exhibits longer tails and the kurtosis is larger than for normal distributions, which is indicative of the out of equilibrium behaviour and correlated flow fields.To quantify the degree of non-equilibrium behaviour, a non-Gaussianity measure (χ NGM ) for the velocity distributions is considered in Fig. 8b.For the lowest activi-ties, χ NGM = 0, corresponding to a true Gaussian distribution for velocity; however, this quickly peaks around α α eq .It is worth noting that our non-Gaussinity measure is directly related to the Binder cumulant [93], and the peak is indicative of discontinuous onset of activityinduced flows.The positive peak in χ NGM corresponds to the single narrower distribution in Fig. 8a.As activity increases past α eq the non-Gaussinity measure drops to negative values, and remains negative throughout the turbulence scaling region α turb < ∼ α < ∼ α † , corresponding to the broadening velocity distributions in Fig. 8a.Finally, for α † < ∼ α, the non-Gaussinity measure increases towards χ NGM → 0, indicating the breakdown of the algorithm.
To quantify the flow structures responsible for the non-Gaussianity, the velocity autocorrelation functions C vv (R) are measured (Fig. 9a).By fitting an exponentially decaying function, the characteristic velocity lengths are obtained (Fig. 9b).For α < ∼ α eq , activity is insufficient to induce hydrodynamic effects.Just as in Fig. 7, for α eq < ∼ α < ∼ α turb hydrodynamic effects are apparent, in that there is a characteristic hydrodynamic length scale.However, there is still insufficient activity to drive active turbulence at these low activities and so the decorrelation length does not yet scale with activity.Only once α turb > ∼ α does the velocity decorrelation length scale decrease with activity.In this regime, the decrease in length scale is seen to correspond to the ideal expectations v ∼ α µ where µ = −1/2 (Eq.8), shown as a dashed line.Comparing to Fig. 5, the director and velocity decorrelation length scale both scale as expected within the fullly formed active turbulence regime, α > ∼ α turb .Correspondingly, the measured length scale collapses the velocity-velocity correlation lengths in the fully formed active turbulence regime of α > ∼ α turb (Fig. 9c).
Active nematics do not possess the scale invariance that inertial turbulence does [32,34,79,94].This is because energy is dissipated at the scale it is injected in, impeding energy cascades.In particular, the length scale from Eq. 8 represents the characteristic vortex size and so considering the amount of enstrophy as a function of wave number (Eq.11), reveals the spatial structure of active turbulence.At small wave numbers, the enstrophy rises with wave number rise ∼ k (Fig. 10); while at intermediate wave numbers, the enstrophy decreases as ∼ k −2 (Fig. 10).This change of sign in the scaling of the enstrophy spectra is expected to occur at k > k α , the wave number representing the characteristic vortex size [35,78].The vortex structure results in a non-monotonic enstrophy spectra with a maximum corresponding to the characteristic vortex diameter.At the largest wave numbers, the signs of the MPCD cell discretization size a are apparent as a secondary peak for k ≥ k a ≈ 2π/3a (Fig. 10).This is because nearestneighbour MPCD cells are used to calculate the local vorticity.
These results demonstrate the AN-MPCD algorithm generates spontaneous flows for sufficiently large activities.The scale of these flows increases rapidly at first, when the ratio of system size and activity are insufficient for fully developed turbulence, but for sufficient activity the characteristic velocity grows according to the theoretically expected scaling.The autocorrelation functions and enstrophy spectra demonstrate that AN-MPCD reproduces the essential properties of fully developed active turbulence with flow structures that statistically match expectations.However, AN-MPCD is a point particlebased mesoscale method and, as such, also possesses properties expected of active particle models.

D. Density Fluctuations
As a mesoscale algorithm between the continuum and microscopic limits, AN-MPCD exhibits both the active turbulence predicted by continuum approaches and the density fluctuations of particle models.Since MPCD is composed of point-particles and obeys an ideal gas equation of state [60][61][62]95], it is a compressible fluid, as is illustrated by the probability distribution of the density (Fig. 11a).In the low activity limit α α eq , the fluid possesses a relatively broad distribution of particles with a mean ρ = N c /a 2 = 20 and a standard deviation σ Nc 4.5 (Fig. 11a, b).At the lowest activities, the mean coincides with the mode at N c = 20, and the distribution is relatively Gaussian, with a slight positive skew.At these low activities, the distribution does not change as a function of activity; however, once α ≥ α turb , the standard deviation begins to rise (Fig. 11b).
The distribution broadens but also changes shape with higher moments of the distribution growing more rapidly (Fig. 11a).Once the standard deviation begins to grow, the non-Gaussianity parameter χ NGM also increases with activity (Fig. 11b; inset).This represents the widening tails of the density distributions in Fig. 11a.At high activities, the likelihood of finding regions with either much higher or much lower density than the mean is increased.This can be seen at high activities α > ∼ α † with not only empty cells, but also high-density bands within a sparse nematic gas (Fig. 12 and Movie 8).As the activity is raised further, corresponding with broader density distributions, these bands appear narrower as a consequence (Movie 9).The bands exhibit spatiotemporal chaos through elongation, splitting and merging, reminiscent of particle-based dry active nematic models [18,80,96].However, the AN-MPCD algorithm is not observed to exhibit as substantial a decrease in nematic ordering in the gas-like phase (Fig. 13).
To better understand the variation of MPCD particles, consider the number fluctuations within sub-domains.From Eq. 13, the fluctuations scale with the mean as σ Nc ∼ N c ν with ν = 1/2 in thermal equilibrium.At low activities (α < ∼ α turb ), the fluctuations obey the Central Limit Theorem with ν = 1/2 (Fig. 11c,d).However, for activities in the fully developed turbulence regime, density fluctuations increase and crossover to anomalous is the instantaneous maximum cell population.AN-MPCD exhibits high density bands within a dilute nematic gas.behaviour.As ν → 1, giant number fluctuations dominate the statistics, indicating highly out-of-equilibrium behaviour within the active turbulence scaling regime α turb < ∼ α < ∼ α † (Fig. 11d).The eminence of nonequilibrium particle number fluctuations in AN-MPCD highlights the mesoscale nature of the algorithm.AN-MPCD possesses properties of both active particle ensembles and active fluids.
The exponent ν has a maximum at α = α † 3 × 10 −1 .The decrease in ν past this point can be explained by considering how the average nematic order parameter S varies with activity (Fig. 13).For the near-equilibrium regime of α < ∼ α eq , the nematic order in the system is at its highest and is effectively constant.The system then transitions to fully-formed active turbulence regime α eq < ∼ α < ∼ α † and the nematic order transitions to another plateau at α α turb .At the highest activities, the nematic order drops off sharply for α † < ∼ α, corresponding to the decrease in giant-number fluctuations, ν (Fig. 11d; .The high activity disorders the orientation partially by causing local regions of low density, which are effectively below the isotropic-nematic transition point, leading to a drop in the nematic order.Conversely, the drop in nematic order causes the active nematic dipole to be more broadly distributed which disperses MPCD particles more randomly, leading to less clustering and slightly more homogeneous density structure.V. CONCLUSIONS We have proposed and quantified a mesoscopic, particle-based algorithm for simulating wet active nematics.It is an extension of the nematic version of the multi-particle collision dynamics (N-MPCD) method to account for active force dipoles.This active-dipole contribution to the MPCD collision operator injects kinetic energy but conserves translational momentum.The activedipole is aligned with the local nematic director.The strength of the force dipole is computed from the particles within each cell, resulting in a density dependent local dipole strength.
This active-nematic momentum collision operator generates spontaneous flows, and hydrodynamic instabilities leading to defect unbinding and active turbulence.Activity is found to be negligible below α eq , because the Andersen thermostat absorbs the injected energy.In this negligible-activity limit, the fluid is indistinguishable from a passive nematic.For weak activity (α eq ≤ α ≤ α turb ), spontaneous flows arise but the characteristic length scale of the activity is comparable to the system size and so active turbulence is not fully developed.In this regime, defect-pair unbinding is rare and flows correlations span the entire system.Only as the activity approaches α turb , does the fluid begin to exhibit the aspects of fully developed active turbulence.In the active turbulence regime, the characteristic length scale, as measured from the defect separation, scales with activity as predicted by theory.Likewise, the magnitude of the fluid velocity scales with activity with an exponent comparable to the ideal prediction.
In addition to exhibiting the traits of active turbulence as expected for a continuum model of active fluids, AN-MPCD also possesses characteristics of active particle models.Most prominently, the local density of MPCD particles has begun to exhibit giant number fluctuations.The concurrence of active nematic fluid properties and active particle properties highlights the mesoscale nature of AN-MPCD -multi-particle collision dynamics is a coarse-grained algorithm that spans the intermediate scales between the microscopic and hydrodynamic limits.
Much work has been accomplished studying the fundamental properties of active fluids [3][4][5].However, soft condensed matter physics is often not principally interested in the dynamics of a background solvent itself.Rather, relaxational dynamics, self-assembly, transport of complex solutes embedded within a fluidic medium, or other multiscale phenomenta are the subjects of physical interest.Active-Nematic MPCD opens the door to simulating complex particles embedded within a spontaneously flowing active medium.Future work could simulate of Janus particles [45] or other anisotropic passive colloids, including flexible filaments or even polymeric materials [58].Similarly studies of driven-particles within active nematics [46,47] could be extended.In fact, AN-MPCD could be used to study suspensions of active particles, such as swimming bacteria, suspended in an active solvent with a different activity or symmetry -e.g.pushers in contractile nematics.By bridging microscopic models, such as active Brownian particles, and macroscopic models, such as the Toner-Tu equation, AN-MPCD offers a pathway for numerical studies of novel active hybrid materials.

FIG. 1 :
FIG. 1: Schematic representation of the active portion of the collision operator (Eq.6).(a) Particles are binned into cells.Each particle has an internal orientation u i shown schematically as an ellipsoidal shape and velocity v i denoted by a blue arrow.(b) The local nematic director n c (red arrow) and centre of mass r cm c of the cell define a plane with approximately half the particles above and half below (yellow plane).(c) Each particle is subject to an impulse away from the plane (orange arrows).(d) Particles stream with their new velocities (green arrows).

FIG. 2 :
FIG. 2: Instabilities lead to defect pair creation and active turbulence.Snapshots of the local director field n c (r; t) coloured by scalar order parameter S c (r; t) for a system size of sys = 30 and extensile activity α = 0.03 from Movie 1.(a) Timestep 0: The director is initialised in a fully aligned configuration and thermal noise perturbs the director field.(b) Timestep 110δt: Extensile bend instability causes high-bend kink walls to form perpendicular to the global director with net force density parallel to the bend.(c) Timestep 160δt: Defect pairs unbind within walls.+1/2 and −1/2 defects (red and blue respectively) are highlighted to guide the eye.(d) Timestep 245δt: Steady-state defect population arises, leading to active turbulence.

FIG. 3 :
FIG.3: AN-MPCD exhibits pair creation and annihilation events.Snapshots of the local director field n c (r; t) depicting creation (top row) and annihilation (bottom row) events in extensile AN-MPCD.Snapshots are taken from the same simulation as shown in Fig.2and Movie 1.The frame rate between snapshots is 5δt.Positive +1/2 (red) and negative −1/2 (blue) defects are highlighted to guide the eye.(Top row): (i) An initial kink wall.(ii-iii)The bend instability results in net force parallel to the bend, exacerbating the deformation.(iv) Spontaneous creation of a ±1/2 defect pair.(v) The +1/2 defect is self-motile and departs from the creation event site.(Bottom row): (i) Proximate +1/2 and −1/2 defects.(ii) The self motile +1/2 defect approaches the −1/2 defect.(iii) As the defects approach, they begin to lose their distinct shapes.(iv) Defect annihilation reduces nematic distortion.(v) Locally ordered nematic remains after an annihilation event.

FIG. 8 :
FIG. 8: Activity broadens velocity distributions in AN-MPCD (a) Probability distributions of the local average of velocity components, rescaled by the standard deviation, for λ = 2 and sys = 200.The red dashed line is a standard Gaussian distribution.(b) Non-Gaussianity measure χ NGM of the distributions in Fig. 8a.

FIG. 10 :
FIG. 10: Enstrophy spectral structure of active turbulence.Enstrophy spectra for sys = 400.Scalings of ∼ k 1 at low wave numbers (blue dashed line) and ∼ k −2 at intermediate wave numbers (red dashed line) guide the eye.At the highest wave numbers there exists a secondary structure due to the MPCD cell size a, which appears because vorticity is computed on the MPCD lattice.The wave number at which vorticity computed is marked as the vertical black dashed line.

FIG. 11 :
FIG. 11: AN-MPCD exhibits giant-number fluctuations.(a) Probability density functions (PDF) of the cell density.The average cell density N c = 20 is shown as a dashed red line.(b) Standard deviation of the cell density as a function of activity.(inset) Non-Gaussinity measure χ NGM of cell density PDF Fig. 11a.(c) Standard deviation of population shown as a function of the average cell population.Scalings of σ Nc ∼ N c ν with ν = 1/2 (black dashed line) and ν = 1 (red line).(d) The scaling exponent ν of Fig. 11c measured at the low cell population limit.

FIG. 12 :
FIG. 12: Density Fluctuations in Active-Nematic MPCD.Snapshot of the cell density field of an AN-MPCD system ( sys = 100), α = 0.08, and average cell density N c = 20 from Movie 8. N max C