Rheology of a concentrated suspension of spherical squirmers: monolayer in simple shear flow

A concentrated, vertical monolayer of identical spherical squirmers, which may be bottom-heavy, and which are subjected to a linear shear flow, is modelled computationally by two different methods: Stokesian dynamics, and a lubrication-theory-based method. Inertia is negligible. The aim is to compute the effective shear viscosity and, where possible, the normal stress differences as functions of the areal fraction of spheres $\phi$, the squirming parameter $\beta$ (proportional to the ratio of a squirmer's active stresslet to its swimming speed), the ratio $Sq$ of swimming speed to a typical speed of the shear flow, the bottom-heaviness parameter $G_{bh}$, the angle $\alpha$ that the shear flow makes with the horizontal, and two parameters that define the repulsive force that is required computationally to prevent the squirmers from overlapping when their distance apart is less than a critical value $\epsilon a$, where $\epsilon$ is very small and $a$ is the sphere radius. The Stokesian dynamics method allows the rheological quantities to be computed for values of $\phi$ up to $0.75$; the lubrication-theory method can be used for $\phi>0.5$. A major finding of this work is that, despite very different assumptions, the two methods of computation give overlapping results for viscosity as a function of $\phi$ in the range $0.5<\phi<0.75$. This suggests that lubrication theory, based on near-field interactions alone, contains most of the relevant physics, and that taking account of interactions with more distant particles than the nearest is not essential to describe the dominant physics.


Introduction
Active matter has become a popular area of research in recent years, among both fluid dynamicists and soft matter physicists. The canonical examples of fluid active matter are suspensions of swimming micro-organisms, which may exhibit fascinating phenomena, ranging from steady, regular patterns, as in bioconvection (Wager 1911;Platt 1961;Childress et al. 1975;Kessler 1986;Pedley & Kessler 1992) among others, to random coherent structures, sometimes referred to as bacterial turbulence (Dombrowski et al. 2004), with many variants in between.
There have been few experimental studies on the effective viscosity of suspensions of microscopic swimmers. Notable among those few are the measurements of Rafaï et al. (2010) on suspensions of motile algae (Chlamydomonas reinhardtii ), which are close to spherical, are bottom-heavy, and pull themselves through the fluid (equivalent in the squirmer model to β > 0; see Eq. (1.3) below). The effective viscosity in a horizontal shear flow was found to increase with volume fraction much more dramatically than for a suspension of dead cells. On the other hand, Sokolov & Aranson (2009) measured the effective viscosity of a suspension of bacteria (pushers: β < 0), and found a significant decrease in shear viscosity with swimming speed. The latter finding was reinforced by Lopez et al. (2015), using a Couette viscometer in which the flow is essentially simple shear. They found that the presence of pusher cells (E. coli bacteria) caused the effective viscosity to fall below that of the ambient fluid at sufficiently low values of the shear rate γ, and even to approach zero as the cell volume fraction was increased. The mechanism for this phenomenon was first set out by Hatwalne et al. (2004) and depends on the cells being elongated (prolate spheroids or rods): in the absence of swimming, such rods describe Jeffery orbits and align preferentially with the directions of extensional strain rate. The contractile nature of the flow generated by a puller opposes this extension and therefore the effective viscosity is increased, but the extensile nature of pusher-generated flow enhances it. The mechanism is operative for dilute suspensions; it is very clearly explained in the review by Saintillan (2018). However, the same mechanism does not apply to a spherical organism, which does not have a preferred orientation in simple shear flow.
Observations of phenomena such as those described above have stimulated an equally large range of theoretical models to see if the observations can be explained by physical processes alone, without requiring an understanding of biological (or chemical) signalling or intracellular processes. Continuum models have been very successful for dilute suspensions, in which the cells interact with their environment but not with each other: bioconvection results from either a gravitational instability when the upswimming of dense cells leads to a gravitationally unstable density profile, or a gyrotactic instability in which the cell's non-uniform density or geometric asymmetry causes them to be reoriented in a shear flow (Childress et al. 1975;Kessler 1986;Pedley & Kessler 1992;Roberts & Deacon 2002). Even when gravity is unimportant the stresses applied by the cell's swimming motions (a swimmer acts as a force dipole or stresslet) lead to instability and random bulk motions (Pedley & Kessler 1990;Simha & Ramaswamy 2002;Saintillan & Shelley 2008).
The rheology of a suspension in an incompressible fluid is represented by the bulk stress tensor Σ, which needs to be calculated or measured in terms of the (changing) configuration of the particles in the suspension and of the velocity field. Rational scientific analysis of the rheology of non-dilute suspensions began with the work of Batchelor (1970), who showed that Σ, for force-free particles in a (quasi-)steady linear flow with strain-rate tensor E, could be expressed as (1.1) where the first term is the isotropic part of the stress, P being the effective pressure and the second term is the Newtonian viscous stress (η 0 being the fluid viscosity). The third term is the particle stress tensor, defined as the average over all spheres of the stresslet for a single particle: where σ is the stress tensor and u is the velocity. A p is the surface of the particle with outward normal n. An effective viscosity can then be defined for each off-diagonal component of the (ensemble) average stress tensor by dividing it by the corresponding component of the bulk rate of strain tensor, which shows that in general the effective viscosity is itself non-isotropic. Moreover it will in general also depend on the mean velocity field, which means the suspension is non-Newtonian. Further non-Newtonian effects that are of practical importance and are often calculated are the differences between the diagonal stress components ('normal stress differences').
There is a considerable literature on the computation of the effective viscosity and normal stress differences in suspensions of passive particles, in particular identical, rigid, noncolloidal spheres at low Reynolds number. The first contribution was from Einstein (1906), who calculated the first correction, of O(φ) where φ is the particle volume fraction, to the fluid shear viscosity; the suspension was dilute and there was no interaction between particles. Semi-dilute suspensions, in which pairwise hydrodynamic interactions were included, were considered by Batchelor & Green (1972), and they could compute the O(φ 2 ) term, at least in linear flows with no closed streamlines. Higher concentrations in general require substantial computations; great progress has been made using the method of Stokesian dynamics, introduced by Brady and his colleagues, e.g. (Bossis & Brady 1984;Brady & Bossis 1985); application to the rheology of concentrated suspensions was made by Sierou & Brady (2002). Singh & Nott (2000) applied Stokesian dynamics to a monolayer of rigid spherical particles. An excellent recent review of the rheology of concentrated suspensions is given by Guazzelli & Pouliquen (2018).
In this paper we wish to consider concentrated suspensions of active particles, in which the flow is dominated by cell-cell interactions. Continuum models fail because there is no agreed way of incorporating cell-cell interactions into the model equations -in particular the particle stress tensor Σ (p) -even when the interactions are restricted to near-field hydrodynamics together with a repulsive force to prevent overlap of model cells. We seek to understand the rheological properties of an idealised suspension from direct numerical simulations. The cells are modelled as identical steady spherical squirmers (Lighthill 1952;Blake 1971;Ishikawa et al. 2006;Pedley 2016): spheres of radius a which swim by means of a prescribed tangential velocity on the surface, where θ is the polar angle from the cell's swimming direction, V s is the cell swimming speed and β represents the stresslet strength; inertia is negligible. The suspension is taken to be a monolayer, in which the sphere centres and their trajectories are confined to a single plane, which is vertical in cases for which gravity, g, is important. The monolayer is here taken to be confined to the narrow gap between stress-free planes, spaced a distance 2.1a apart. In previous work on semi-dilute suspensions, the monolayer was taken to be embedded in an effectively infinite fluid ; we show that the results of the two models do not differ greatly. The spheres may be bottom-heavy, so that when the swimming direction of a sphere, p, is not vertical the sphere experiences a gravitational torque T , where and υ, h are the cell volume and the displacement of the centre of mass from the geometric centre; ρ is the average cell density, assumed in this case to be the same as that of the fluid. The monolayer is taken to be driven by a simple shear flow in the same, x-y, plane: U = (γy, 0, 0), with shear-rate γ; we will also take g = −g(sin α, cos α, 0), (1.5) so the flow is horizontal if α = 0. We have previously studied the rheology of semi-dilute, three-dimensional, suspensions of squirmers -volume fraction φ 0.1 -in which pairwise interactions between the cells were summed simply, neglecting interactions involving more than two cells at a time (Ishikawa & Pedley 2007b). General pairwise interactions were computed exactly using the boundary element method, supplemented by lubrication theory when cells were very close together. Cells were computationally prevented from overlapping by the inclusion of a repulsive interparticle force wherer is the unit vector along the line of centres and is the minimum dimensionless spacing between two cells (Brady & Bossis 1985). Ishikawa & Pedley (2007b) took τ equal to 1000 and F 0 was varied. It was found that the swimming activity made very little difference to the effective shear viscosity when the spheres were non-bottom-heavy, but the suspension showed significant non-Newtonian behaviour, such as anisotropic effective shear viscosity and normal stress differences, when the cells were bottom-heavy. Moreover, horizontal and vertical shear flows (α = 0 or π/2) gave very different results.
Here we use two different methods of simulation for very concentrated monolayer suspensions, with areal fraction φ up to 0.8. One is a full numerical simulation using Stokesian dynamics , while in the other we assume that every cell interacts only with its immediate neighbours; in that case the interactions are described using lubrication theory alone. This model was recently used to investigate the stability of a regular array of bottom-heavy squirmers, swimming upwards in the absence of an imposed shear flow (Brumley & Pedley 2019). The point of this is to see if lubrication theory alone is good enough to account for all the interesting rheological behaviour of a concentrated suspension, or if some effect of more distant particles is necessary, in the case of squirmers. As well as illuminating the physics, this might make related computations significantly cheaper than the full Stokesian dynamics. For inert and force-free spheres, Leshansky & Brady (2005) reported in a footnote that suppressing all far field interactions made less than 5% difference to the quantities they were computing.
These methods are explained in more detail in the first parts of the next two sections, and their respective results are presented and discussed in the second parts. The findings are compared in the final section, with further discussion and an outline of intended future work.

Problem settings and numerical method
In this study, we consider a monolayer suspension of squirmers in a thin fluid film. The film is assumed as infinitely periodic, and the two surfaces are assumed to be flat and stress free. The film thickness L z is set as L z = 2.1a throughout the study. The reason why we consider a monolayer suspension, instead of a 3-dimensional suspension, is because it allows us to handle a larger system size with a limited number of particles. Moreover, a film suspension can in principle be generated experimentally by using a soap film. As shown in Fig. 1, in-plane shear flow is induced in the monolayer suspension of squirmers; the x-axis is taken in the flow direction, the y-axis is taken in the velocity gradient direction, and the z-axis in taken perpendicular to the film plane. The background shear flow can be expressed as (u x , u y , u z ) = (γy, 0, 0). Gravity also acts in the x − y plane, and α is the angle the gravitational acceleration g makes with the −y axis, as given by Eq. (1.5) and shown in the figure. We assume that body centres and orientation vectors of squirmers are placed in the centre-plane of the fluid film. Thus, due to the symmetry of the problem, the squirmers remain in the same centre plane all the time.
The hydrodynamic interactions among squirmers in an infinite periodic monolayer suspension are calculated in the same manner as in our former studies , so we explain the methodology only briefly here. The Stokesian dynamics method ) is employed. The hydrodynamic interactions among an infinite suspension of particles are computed by the Ewald summation technique. By exploiting the Stokesian dynamics method, the force F , torque T and stresslet S of squirmers are given by: where R is the resistance matrix, U and Ω are the translational and rotational velocities of a squirmer, u and ω are the translational and rotational velocity of the bulk suspension, and E is the rate of strain tensor of the bulk suspension. Q sq is the irreducible quadrupole providing additional accuracy, which is approximated by its mean-field value (cf. ). Index far or near indicates far-or nearfield interaction, and 2B or Sq indicates interaction between two inert spheres or two squirmers, respectively. The computational region is a square of side L, and a suspension of infinite extent is modelled with periodic boundary conditions in the x and y directions. In order to express the stress free surfaces of the fluid film, the monolayer is periodically replicated also in the z-direction with 2.1a intervals. For the simple shear flow in the x, y-plane, the periodic conditions in the x and z-directions are straightforward. In the y-direction, the periodicity requires a translation in the x-direction of the spheres at y = L, relative to those at y = 0, by an amount Lγt in order to preserve the bulk linear shear flow, where t is time. The number of squirmers N in a unit domain is set as N = 128. The areal fraction of squirmers is φ (not the volume fraction c), and it is varied in the range 0.1 φ 0.75. Parameters in Eq. (1.6) for the repulsive force are set as F 0 = 10 and τ = 100. The initial configuration of the suspended squirmers is generated by first arranging them in a regular array and then applying small random displacements and random orientations.
The dynamic motions are calculated afterwards by the fourth-order Adams-Bashforth time marching scheme. The computational time step is set as ∆t = 5 × 10 −4 /γ.
The apparent shear viscosity η of the film suspension can be calculated from the suspension averaged stresslet S . The excess apparent viscosity is given by where c is the volume fraction. The areal fraction φ and c satisfy the relation c = 4aφ/(3L z ), where L z (= 2.1a) is the film thickness. When squirmers are torque-free, the stresslet is symmetric and η xy = η yx = η. When squirmers are bottom-heavy, on the other hand, the stresslet becomes asymmetric and η xy = η yx . The dimensionless first and second normal stress differences are defined as In order to obtain suspension-averaged properties, stresslet values are averaged over all particles during the time interval 15 tγ 30, given that the ensemble values reach the steady state. We introduce a dimensionless number Sq, which is the swimming speed, scaled with a typical velocity in the shear flow, i.e.
In simulation cases with Sq > 1, parameters are modified as F 0 = 10Sq, ∆tγ = 5 × 10 −4 /Sq, to prevent overlapping of particles. Stresslet values are averaged during the longer time interval 75/Sq tγ 150/Sq, to obtain the steady state values.

Results for non-bottom-heavy squirmers
We first calculated the apparent viscosity η of the suspension, for both inert and squirming spheres. The results are plotted as a function of areal fraction φ in Fig. 2(a). We see that η of a suspension of squirmers with Sq = 1 and β = +1 (i.e. pullers) increases rapidly with φ. By comparing with the results for inert spheres in the present study, it is found that the motility of a squirmer considerably increases the viscosity. In the case of inert spheres, layers of particles are formed along the flow direction, and nearfield interactions between particles are reduced. In the case of squirmers, on the other hand, such layers are destroyed by the motility of the squirmers, which move around irregularly, leading to frequent near-field interactions. Near-field interactions generate large lubrication forces, which result in large particle stresslets as well as the large apparent viscosity.
In Fig. 2, Einstein's equation in the dilute limit as well as the former numerical results of Singh & Nott (2000), for inert spheres, are also plotted for comparison. We see that all results agree well with Einstein's equation in the dilute regime. The main differences between the present study and Singh & Nott (2000) are the boundary condition and the repulsive interparticle force. In Singh & Nott (2000) a monolayer suspension is sheared between two walls, while in the present simulation the shear is generated without a wall boundary. The coefficients of repulsive force given by Eq. (1.6) are η 0 a 2 γF 0 = 10 −4 [N ], τ = 100[−] in Singh & Nott (2000), while F 0 = 10[−], τ = 100[−] in the present simulation. We see that the apparent viscosity of a suspension of inert particles reported by Singh & Nott (2000) is larger than that in the present study. The discrepancy may come from the differences in the boundary condition and the repulsive interparticle force. By introducing the wall boundary, motions of spheres in the y-direction are restricted. Figure 1. Problem setting of the simulation. An infinitely periodic monolayer suspension of squirmers in a thin fluid film is sheared in the x − y plane. The unit domain is a square with side length L, and contains 128 squirmers. The film is assumed to be flat with thickness Lz = 2.1a, and has two stress free surfaces. α is the angle of gravitational acceleration g taken from the −y axis.
The restricted motions may enhance the near-field interactions between particles, which results in the larger apparent viscosity. Moreover, the repulsive force may be larger in the present study than in Singh & Nott (2000) (it is hard to find the value of the dimensionless coefficient analogous to F 0 in that paper). This is because squirmers easily overlap if the repulsive forces are not strong enough, given that lubrication flow between two squirming surfaces cannot mathematically prevent the overlapping. The strong repulsive force may reduce the apparent viscosity, which may be another reason for the discrepancy.
The stresslet can be decomposed into the hydrodynamic contribution and the repulsive contribution, as shown in Fig. 2(b). The repulsive contribution to the particle bulk stress can be calculated as (Batchelor 1977 where V is a unit volume, r ij is the centre-centre separation of squirmers i and j, and F ij is their pairwise interparticle force given by Eq. (1.6). We see that the hydrodynamic contribution is dominant in the low φ regime. When φ 0.6, on the other hand, the repulsive contribution becomes larger than the hydrodynamic contribution, which is the reason why the apparent viscosity rapidly increases in the high φ regime.
First and second normal stresses also appear in the suspension of inert spheres and squirmers, as shown in Fig. 3. N 1 is negative in sign, so particles are basically compressed in the flow direction. |N 1 | increases as φ is increased, similar to the apparent viscosity.  (2000) and Einstein's equation of (η − η0)/η0 = 2.5c = (10/6.2)φ are also plotted for comparison. (b) Present results of squirmers with Sq = 1 and β = 1 decomposed into the hydrodynamic contribution and the repulsive contribution.
The value of |N 1 | in a suspension of inert particles, as reported by Singh & Nott (2000), is larger than in the present study. The discrepancy may again come from the differences in the boundary condition and the repulsive interparticle force.
N 2 in the suspension of inert spheres is also negative in sign, so the particle stresses satisfy Σ zz . In the case of squirmers, the sign of N 2 changes at around φ = 0.55. The hydrodynamic contribution to N 2 is positive and steadily increases with φ, while the repulsive contribution to N 2 is negative and rapidly increases with φ. Since the repulsive contribution overwhelms the hydrodynamic contribution when φ 0.6, N 2 becomes negative in the large φ regime.
The difference between pushers and pullers can be investigated by changing the swimming mode parameter β, which is positive for pullers and negative for pushers. Figure 4(a) shows the effective viscosity η for Sq = 1 and various values of β. We see that η increases as the absolute value of β is increased. This is probably because strong squirming velocity with large |β| enhances near-field interactions between squirmers, which induces a strong stresslet. The effect of ±β is not symmetric, but pullers generate larger viscosity than pushers. In order to confirm these tendencies, we calculated the normalised probability density function of squirmers, defined as where P (r 0 |r 0 + r)dr is the conditional probability that, given that there is a squirmer centred at r 0 , there is an additional squirmer centred between r 0 + r and r 0 + r + dr. The results are plotted in Fig. 4(b). We see that I(r) with β = 3 has a considerably larger value than that with β = 0 in the small r regime, while I(r) with β = −3 does not. On the other hand, I(r) with β = −3 has a larger peak than that with β = 0. Thus, near-field interactions between squirmers are increased as |β| is increased. The difference between pushers and pullers is obvious in Fig. 4(b): pullers tend to come closer to each other than pushers. A former study reported that the face-to-face configuration is stable for puller  squirmers while unstable for pusher squirmers, though the face-to-face configuration was expressed as a stress-free surface in the paper (Ishikawa 2019). Thus, puller squirmers facing each other can come closer more easily than pushers. We also investigated the effect of Sq, i.e. the ratio of swimming velocity to the shear velocity. The results are shown in Fig. 5. We see that the apparent viscosity η increases as Sq is increased. This is because large swimming velocity, relative to the shear velocity, destroys the formation of layers by the particles, and enhances their near-field interactions. This tendency can be confirmed by Fig. 5(b), in which I(r) with Sq = 10 has larger values than with Sq = 0.1 and 0.5. Thus, increase in the apparent viscosity can be understood as coming from the increase of near-field interactions.

Results for bottom-heavy squirmers
When squirmers are bottom-heavy, cells tend to align in the direction opposite to gravity. To discuss the effect of bottom-heaviness, we introduce a dimensionless number G bh , defined as G bh is proportional to the ratio of the time to swim a body length to the time for the axis orientation p of a non-swimmer to rotate from horizontal to vertical. When G bh is sufficiently large, the gravitational torque due to the bottom-heaviness balances the hydrodynamic torque due to the background shear flow, and squirmers tend to align in a particular direction. The effect of G bh and the gravitational angle α on the apparent viscosity is shown in Fig. 6. Since an external torque due to the bottom-heaviness is exerted on a squirmer, the partcle stress tensor becomes asymmetric, and the xy and yx components become different. We see that G bh considerably affects the value of η. For β = 0 and −3, in horizontal flow (α = 0), η is decreased by the bottom-heaviness, and η yx with β = −3 even becomes negative when G bh 50. The decrease in the apparent viscosity may be caused by two mechanisms: (a) squirmers with large G bh tend to swim in the same direction, and cell-cell collisions are suppressed; and (b) aligned squirmers induce a net squirming stresslet when β = 0, which directly contributes to η xy and η yx . The effect of α is also significant: η with β = 3 takes its maximum values around α = π/8, whereas that with β = −3 takes its minimum values around α = π/8. So the tendencies are almost opposite between the pusher and the puller.
In order to clarify the mechanism of bottom-heavy effects, we here discuss the orientation of bottom-heavy squirmers in shear flow. Figure 7 shows normalised probability density distribution as a function of angle ζ that is defined from the x-axis in the counter clockwise direction as shown in Fig. 7(b). When I(ζ) = 1 for all ζ, the orientation distribution is isotropic. We see that the peak of I(ζ) is higher, and at a larger value of ζ, in the G bh = 100 case than in the G bh = 30 case due to the strong bottom-heaviness. Squirmers with G bh = 100 are oriented with approximately ζ = 0.38π regardless of β. When β = −3, i.e. pushers, the stresslet with α = 0 is directed as in Fig. 7(b) bottom-left, and the apparent viscosity decreases. If the direction of gravity is rotated to α = π/2, as in Fig. 7(b) bottom-right, the direction of the stresslet becomes opposite, and the apparent viscosity increases. These schematics can qualitatively explain the results in Fig. 6. When the squirmers are pullers, the sign of the stresslet is opposite to that shown in Fig. 7(b). Thus, the apparent viscosity increases with α = 0 whereas it decreases with α = π/2, which is again consistent with Fig. 6. We note that the internal configuration, at α = 0, is more regular for β = −3 than for β = +3; this is indicated by the higher peak for β = −3 in Fig. 7(a). For β = +3, large numbers of the squirmers are aggregated and almost jam the system. The first and second normal stress differences are affected by the angle of gravity α, as shown in Fig. 8. The first normal stress difference with β = 3 increases as α is increased, whereas that with β = −3 decreases. The second normal stress difference with β = 3 takes its minimum values around α = 3π/8, whereas that with β = −3 takes its maximum values around α = 3π/8. So the tendencies are again almost opposite between the pusher and the puller. The opposite tendency can be explained by the sign and rotation of the stresslet, as schematically shown in Fig. 7(b).
Last, we investigate the rheology under constant Sq · G bh conditions. This product represents the ratio of the gravitational torque to the shear torque, independently of the swimming speed V s . Thus, a solitary squirmer under constant Sq · G bh conditions is expected to have the same orientation angle relative to gravity. The results for excess apparent viscosity and normal stress differences under the condition of Sq · G bh = 30 are shown in Fig. 9. We see that η is considerably increased in the small G bh regime. The effect of swimming is larger than that of the background shear in the small G bh regime. The large swimming velocity enhances near-field interactions between squirmers, which results in the large apparent viscosity. The normal stress difference N 2 for β = −3 is  also enhanced in the small G bh regime. This is because the active stresslet, caused by the squirming velocity, plays a dominant role compared to the passive stresslet, caused by inert spheres, in the small G bh regime. Hence, the rheology under constant Sq · G bh conditions is strongly affected by Sq especially when it is large. Figure 9. Rheology under the condition of Sq · G bh = 30, in which the orientation angle of a solitary squirmer relative to the gravity is expected to be the same (φ = 0.7, β = ±3, α = 0). (a) Excess apparent viscosity and (b) Normal stress differences.

Problem setting and numerical method
In this Section, we present a complementary method for calculating the rheological properties of a suspension of steady, spherical squirmers. The methodology for measuring the bulk viscosity is similar to that of a rheometer; a suspension of squirmers situated between flat parallel plates is sheared using Couette flow, and the drag force on the plates is measured. The overall dynamics of the active suspension are solved by summing pairwise lubrication interactions between closely-separated squirmers, along with hydrodynamic forces and torques associated with solitary squirmers in shear flow. This approach utilises elements of previous work (Brumley & Pedley 2019), but with several key advances. Firstly, in the present formulation, an arbitrary areal fraction of squirmers is permitted, so that cells are no longer confined to a hexagonal crystal array. Secondly, the entire suspension will be subject to a background shear flow. The dependence of the rheology on a number of key suspension and external parameters are presented. Despite the key differences between this approach and that of Section 2, the results are in good quantitative agreement at sufficiently large areal fractions.
Two infinite no-slip plates situated at y 1 = +H/2 and y 2 = −H/2 move with velocities γy 1 e x and γy 2 e x , respectively, so that the fluid between the plates is subject to a uniform shear with rate γ. That is, the fluid velocity is given by u = (u x , u y , u z ) = (γy, 0, 0). A suspension of N identical squirmers is introduced into the fluid (see Fig. 10).
The position and orientation vectors of all squirmers are restricted to lie in the xy plane, so that each squirmer has 2 + 1 degrees of freedom. Moreover, the domain is subject to periodic boundary conditions in the x-direction, with period L such that the squirmer areal fraction φ = N πa 2 /LH, where N = n x n y . This can also be expressed in the following way: where 0 a is the spacing between adjacent squirmers in the case where they are arranged in a regular hexagonal array (as depicted in Fig. 10). The total dimensional force F i and torque T i on the i th squirmer are composed of several contributions, as outlined below. As in Brumley & Pedley (2019), these will be scaled by η 0 πa and η 0 πa 2 respectively, so that they have units of velocity. As in Section 2, the parameter η 0 is the viscosity of the solvent around the squirmers. The resistance formulation for the spheres in Stokes flow is given by where the resistance matrix is given by R = R sq-sq + R sq-wall + R drag . The vectorsF = F /(η 0 πa) andT = T /(η 0 πa 2 ) are of length 2N and N respectively, and are composed of the forces [F Explicit hydrodynamic coupling between squirmers is limited to lubrication interactions. For two squirmers i and j, with minimum clearance ij = (|x i − x j | − 2a)/a (subject to periodic boundary conditions), lubrication interactions occur only if ij < 1. This value is chosen because log ij , the functional dependence of these forces and torques, is equal to zero at that point. These terms therefore increase continuously from zero as squirmers approach one another from afar. The matrix R sq-sq captures the hydrodynamic forces and torques acting on every pair of spheres sufficiently close to one another, arising from their linear and angular velocities. These expressions, evaluated for no-slip spheres, can be found in Kim & Karrila (2005) and Brumley & Pedley (2019). In a similar fashion, R sq-wall captures the forces and torques due to motion of the spheres in close proximity to the bounding walls. The matrices R sq-sq and R sq-wall , composed of blocks for each pair of squirmers, are non-zero only for pairs that are sufficiently close to permit lubrication interactions. The sparsity pattern of these resistance matrices therefore depends on the physical configuration of the system, and must be computed at every time step. Conversely, the final component of the resistance matrix is always diagonal, and is given by where I n is the n × n identity matrix. This captures the drag on a solitary translating and rotating sphere in Stokes flow. Inclusion of this matrix ensures that even widely separated spheres that do not experience lubrication interactions, are subject to solitary Stokes drag.
There are a number of contributions to the forces and torques that are independent of the squirmer velocities. First, there are the contributions which depend on the arrangement of cells with respect to one another. The squirming motions generate contributions for all pairs that are sufficiently close to one another. We refer the reader to Brumley & Pedley (2019) for detailed expressions of these forces and torques,F sq andT sq , respectively, noting that a change of reference frame must be made for each pair, in order for the expressions to apply. Closely separated pairs of squirmers and squirmers near the no-slip plate also experience a repulsive forceF rep parallel to the vector joining their centres (or perpendicular to the wall). This follows the same functional form as in Eq. (1.6), but in principle we allow squirmer-squirmer and squirmer-wall interactions to have different strengths F 0 and interaction ranges τ −1 .
There are several forces and torques which do not require pairwise geometries. Every squirmer is subject to a propulsive force parallel to its orientation vector, p i , according toF where V s = 2B 1 /3 is the swimming speed of a solitary squirmer. Each squirmer also experiences a gravitational torque in the z-direction due to bottom-heaviness: where g = −g(sin α, cos α, 0) and g is the strength of gravity. Finally, the effect of the background shear flow is to exert a hydrodynamic force and torque on each squirmer as follows: where y i is the y coordinate of the i th squirmer. Before proceeding, it is instructive to consider the dynamics of the system if interactions between squirmers and the bounding plates are completely neglected. Under these conditions, the matrix system Eq. (3.3) reduces to the following: Since R drag is diagonal, it is easily inverted, and the following results are obtained: (3.10b) Equation (3.10a) reveals that a solitary squirmer will swim at speed V s in the direction of its orientation, and be advected by the background shear flow in the x-direction. Similarly, the orientation of the squirmer will evolve according to the gravitational torque and background vorticity (Jeffery 1922). In addition to these solitary dynamics (Eq. (3.10)), the complete resistance formulation in Eq. (3.3) includes hydrodynamic and steric effects through pairwise interactions, but only for sufficiently close squirmers. The ability for "solitary" squirmers -i.e. squirmers with no neighbours within lubrication range -to propel themselves is a critical advancement of the present model. This maintains realistic behaviour at low areal fractions of the suspension, when it is quite feasible that a squirmer i will have ij > 1 ∀ j, and ensures that the matrix system in Eq. (3.3) is well conditioned for all configurations.
The dynamics of the sheared squirmer suspension are calculated numerically by solving Eq. (3.3) with the Matlab solver ode15s. The squirmers, initially distributed on a hexagonal close-packed array with mean spacing 0 , are each subject to a random translational perturbation within a disk of radius 0 /2, ensuring that cells are nonoverlapping. Initially, the squirmers' orientations are taken to be random.

Calculation of shear viscosity
For the configuration presented in Fig. 10, determining the effective suspension viscosity requires calculating the wall shear stress on each of the no-slip plates at y = ±H/2. We emphasise that only squirmers which are sufficiently close to the walls to facilitate lubrication interactions will contribute to the wall shear stress. Of the full set of squirmers S = {1, 2, . . . , N }, the following subsets can be identified: The sets S T and S B identify squirmers that have a clearance of less than a with the top and bottom plates respectively (i.e. < 1), and therefore whose behaviour contributes to the lubrication forces and torques. The x-component of the force on the bottom plate is given byF where the three terms represent the x-component of the force on the plate, due to the squirming motion of the i th sphere, the difference between the squirmer velocity, V i , and the plate velocity, −(Hγ/2)e x , and the angular velocity of the squirmer, respectively. In Eq. (3.14a) we have made use of the fact that B 1 = 3 2 aγSq and B 2 = βB 1 . Similar expressions exist to calculate the tangential force,F x T , on the top plate. Here, the value B i a represents the minimum clearance between the i th squirmer and the bottom plate. The repulsive force acts normal to the wall, and therefore does not contribute to the shear force. Although the background shear and gravitational torques influence the motion of the squirmers, their combined effects are encapsulated in Eqs. (3.14b) and (3.14c).
In dimensional form, the tangential shear stress on the top and bottom plates are given by σ T = η 0 πaF x T /(Lδ) and σ B = η 0 πaF x B /(Lδ), respectively, where L = (2 + 0 )an x and δ = 2.1a is the thickness of the monolayer (see Section 2). The effective bulk viscosity is therefore equal to The above expression utilises the mean value of the shear stress across both plates. The excess apparent viscosity is therefore equal to Upon first glance, it appears that the viscosity depends on the system size through the denominator in Eq. (3.16). Although the number of squirmers interacting with the wall would depend on suspension micro-structure, to leading order the number of terms in Eq. (3.13) is proportional to n x . Moreover, the terms in Eq. (3.13) have the dimensions of velocity, matching the dimensions of aγ in Eq. (3.16). In what follows, we will compare the results of Section 2 using Eq. (2.2), with the present lubrication formulation Eq. (3.16).
We should note here that the lubrication-theory-based (LT) method does not give an obvious way to compute normal stresses, so the results that follow will concentrate on predicting the shear viscosity.

Results for non-bottom-heavy squirmers
We first calculate the excess apparent viscosity in the absence of gravity. Suspensions of N = 132 spheres with an areal fraction of φ = 0.80 were simulated over t ∈ [0, 150], for both active squirmers (Sq = 1, β = 1) and passive spheres (Sq = 0). The suspension viscosity can be calculated at every time-step using Eq. (3.16), the results of which are plotted in Fig. 11(a). It is evident that the squirmer suspension (blue) has a higher mean viscosity than the suspension of passive spheres (black), and also exhibits greater excursions from the mean value.
The contribution to the mean squirmer suspension viscosity (blue curve in Fig. 11(b)) arising purely from linear and angular velocities is η/η 0 − 1 = 20.6, approximately 90% of the total mean viscosity (η/η 0 − 1 = 22.8). This is still considerably higher than the mean value for the passive sphere case η/η 0 − 1 = 12.0. Despite the propensity for active swimmers to redistribute themselves, the interior interactions still give rise to a higher suspension viscosity.
The areal fraction of cells was systematically varied in the range 0.5 < φ < 0.8 for both squirmers and passive spheres. A sufficiently long averaging window (10 < t < 150) was taken when calculating the time-averaged suspension viscosity. Figure 11( b) shows the results of the lubrication simulations (circles) together with the findings using Stokesian dynamics (dashed, cf. Fig. 2). There is good quantitative agreement between the complementary simulation methods across all values of φ studied.
The effect of the squirmers' swimming properties was also studied. The two dimensionless parameters are Sq = V s /aγ, the swimming speed relative to background shear rate, and β = B 2 /B 1 , which effectively controls the stresslet sign and strength. Figure 12(a) for β = 1 shows that the suspension viscosity increases for small Sq before plateauing around Sq = 1. Increasing Sq further does not result in a noticeable shift in the viscosity, on average.
The dependence of the suspension viscosity on the swimming mode β is shown in Fig. 12(b) for two different areal fractions. The viscosity is higher for φ = 0.7 than for  φ = 0.6 for all values of β studied. The viscosity for pushers (β < 0) does not vary significantly with β, whereas for pullers (β > 0) the viscosity increases dramatically with β. This increase is much more marked than the Stokesian dynamics findings (see Fig. 4(a)). Cross channel probability distributions for the squirmer positions reveal that pullers spend more time near the boundaries than pushers do, and this likely has a corresponding effect on the suspension viscosity.

Results for bottom-heavy squirmers
The effect of bottom-heaviness is to provide an external torque on each squirmer, which tends to reorient the cell in the gravitational field (see Fig. 10). In the absence of hydrodynamic interactions or external shear, the orientation of each squirmer p i will become anti-aligned with the gravity vector g over a timescale that is inversely proportional to the normalised strength of gravity, G bh . A series of simulations were performed in which the strength of gravity and the angle with respect to the shear flow, were independently varied. Figure 13(a) illustrates the excess apparent viscosity as a function of G bh for three different squirming modes, β (all with Sq = 1, φ = 0.7, α = 0). For small G bh , the ordering matches the results of Fig. 12(b), since hydrodynamic effects alone result in the strongest accumulation of pullers near the wall. As G bh is increased, the viscosity differences become more pronounced, with both neutral squirmers and pushers exhibiting η/η 0 −1 ≈ 0, i.e. no enhancement over the solvent viscosity. The internal squirmer motions are similar to those computed by the SD method: pullers with high G bh generate nearjamming aggregates. As the gravitational angle, α, is varied (for fixed G bh = 100), the curves exhibit cross-over points (see Fig. 13(b)). For −π/4 α π/4, puller suspensions exhibit the greatest viscosity. However, for α π/4 and α −π/4, pusher suspensions have a greater viscosity. In Fig. 13(b), the large peak in the β = 3 curve around α = −π/4, corresponds to a situation where the (clockwise) viscous torques due to background shear balance the (counterclockwise) gravitational torques for upward-swimming squirmers. This results in a system which is, to leading order, equivalent to the G bh = 0 case (1% difference between respective viscosities).
The dependence on α can be understood in terms of the fluid mechanics of propulsion in each case. For α ≈ 0, squirmers will tend to align vertically in the shear flow (i.e. perpendicular to the channel), so that their poles are closest to the walls. Moreover, since squirmers swim upwards, the lubrication interactions will be dominated by anterior poles interacting with the upper plate. Under these conditions, pullers will be drawn closer to the wall than pushers, and so will impart a greater lubrication drag on the top plate. For α ≈ π/2, gravity will tend to align squirmers parallel to the walls, so that equatorial regions of the squirmers are in the lubrication zones. The movement of fluid away from the poles in the case of pushers, means that under these circumstances, they will be drawn closer to the translating plates than their puller counterparts, and therefore exhibit a higher viscosity.

Discussion
The main focus of this paper has been to calculate the particle stress tensor in a concentrated suspension of spherical squirmers (modelling swimming cells) in a planar monolayer exposed to a uniform shear flow, over a wide range of parameter values, in the hope that the union of such results could act in place of an analytical constitutive relation, which appears unlikely to be achievable. For non-bottom-heavy squirmers, the particle stress tensor is symmetric and the rheology can be represented in terms of a single effective shear viscosity η, together with relatively small normal stress differences (Fig. 3). The suspension is non-Newtonian because of these and the fact that η depends on the value of the shear rate γ. However, for bottom-heavy squirmers the particle stress tensor is asymmetric (η xy = η yx ) and exhibits significant normal stress differences: the suspension is clearly non-Newtonian. We have principally investigated how the time average of η (= (1/2)(η xy + η yx ) when the two off-diagonal components are different) varies with the areal fraction φ, the squirming parameter β, the ratio Sq of swimming speed to a typical speed of the shear flow, the bottom-heaviness parameter G bh (Eq. (2.7)), the angle α that the shear flow makes with the horizontal, and the two parameters F 0 and τ that define the repulsive force that is required computationally to prevent the squirmers from overlapping when their distance apart is too small.
We have employed two different numerical methods, Stokesian dynamics (SD), which takes account of all cell-cell interactions, including those between distant cells although high-order multipoles are neglected, and a lubrication-theory-based method (LT) that ignores all cell-cell interactions except those between a cell and its closest neighbours, which are calculated using the approximations of lubrication theory. Both computations are started from given initial conditions and results are taken when the effective viscosity has reached a statistically steady state. Both also use periodic boundary conditions in x. In the SD method, the underlying shear flow is imposed by applying a y-dependent translation to all the spheres, so that those at y = L, say, are displaced in the x-direction by an amount Lγt relative to those at y = 0; periodic boundary conditions in y can then be applied. The particle stress tensor is given by the average over all spheres of the stresslet of an individual sphere (Eq. (1.2)). In the LT method, the shear is applied by translating two infinite planes at y = ±H/2 parallel to each other with velocity ±Hγ/2, and the suspension is set into motion by the viscous shear stresses σ B (bottom wall) and σ T (top wall) exerted on them by the spheres nearest to those planes. The effective viscosity is η = (σ B − σ T )/2γ. Thus the method of computing η is very different between the two methods. They would give the same values only if most of the relevant physics took place in the interior of the suspension and not at the y-boundaries. In other words, although the only forces to be calculated in the LT method are the shear stresses at the walls, the viscosity depends on the internal interactions which drive changes to the configuration of the whole suspension by which the behaviour of the cells near the boundary are determined.
As shown in Fig. 11(b), there is good agreement between the values of η derived by the two methods for non-bottom-heavy squirmers, over the areal fraction range 0.5 < φ < 0.75. Since the SD method takes into account multi-body interactions on top of the pairwise lubrication interactions, the agreement indicates that the relevant interactions are basically pairwise in the regime 0.5 < φ < 0.75 and confirms the basic validity of the lubrication theory approach. The variation of η with various parameters also shows qualitative agreement between the two methods, though not good quantitative agreement: compare Figs. 5(a) and 12(a) for η(Sq), Figs. 6(a) and 13(a) for η(G bh ), and Figs. 6(b) and 13(b) for η(α). In particular, the variation of η with β for pullers (positive β) is much more marked according to the LT method ( Fig. 12(b)), even for non-bottom-heavy squirmers. To understand this, consider the respective interactions of pushers (β < 0) and pullers (β > 0) with a translating wall. The cells oriented towards the wall will have a tendency to reside there longer than cells pointing away from the wall, though this interpretation can obviously break down in various situations. It follows that Table 1. Effect of the repulsive force on the excess apparent viscosity (η − η0)/η0. Particles are non-bottom-heavy and the areal fraction is φ = 0.5. Results are shown for inert spheres and squirmers (with β = 1, Sq = 1), for both SD and LT methods. the anterior hemispheres of the squirmers are most likely to constitute the lubrication interaction with the wall. As squirmers are "tilted" clockwise by the translating boundary, pullers will more strongly oppose the motion of the plate, thereby leading to a higher viscosity. This mechanism partly overrides the similarity of the internal configurations; it is a new physical mechanism, which, unlike those considered in previous studies, is not based on cell elongation. It was seen in Fig. 2(b) that the contribution to the effective viscosity of the repulsive forces between squirmers becomes larger at high volume fraction (φ 0.6) than the contribution from hydrodynamic forces. The results in that figure were obtained for β = 1, Sq = 1, using our 'standard' values of repulsive force parameters, F 0 = 10, τ = 100 (F 0 controls the strength of the repulsive force, whereas 1/τ represents the distance over which the force plays a role). The dependence of the excess effective viscosity on the repulsive forces, calculated using the SD and LT methods, is shown in Table 1. The viscosity is only slightly affected by the parameters F 0 and τ . On the other hand, the difference between the squirmers and inert spheres is pronounced, and of approximately the same magnitude, across all force combinations studied, indicating that the results are not critically dependent on the values of these parameters. The difference between squirmers and inert spheres is presumably induced by two mechanisms: a) the surface squirming velocity generates a direct contribution to the stresslet, and b) the squirming motion determines the suspension microstructure in the first place.
The magnitude of the repulsive force has a similar moderate influence on the apparent viscosity in the LT simulations. Both the magnitude F 0 , and characteristic length scale 1/τ , of the repulsive force influence the spacing, , between the bounding walls and the squirmers adjacent to them. Since the effective viscosity is determined solely by lubrication interactions which scale as log (see Eq. (3.13)), the repulsive force does influence the measured viscosity slightly. However, as in the SD method, we emphasise that the differences between the passive spheres and the squirmers prevail regardless of the specific choices of F 0 and τ .
Although the present paper is mainly concerned with the mean rheological properties, we also analysed time-dependent viscosity of the suspensions (see for example Fig. 11(a)). The squirmer suspensions always exhibit a higher mean viscosity than the passive sphere case. Furthermore, the fluctuations about the mean, shown explicitly in Fig. 11(a) and displayed as confidence intervals in subsequent figures, increase with areal fraction φ. It is appropriate to note here that the maximum packing fraction for spheres in a monolayer is φ ≈ 0.91.
The model active suspension analysed in this paper is extremely idealised: a pla-nar array, of identical, non-colloidal, spherical cells, which swim through the fluid by squirming with an unchanging distribution of tangential velocity on their surfaces. Their concentration (areal fraction φ) is high (up to 0.8). It follows that there is not much previous literature against which the results can be compared, either computational or, especially, experimental. As discussed in Section 2.2, the only monolayer predictions of suspension rheology that we could find are those of Singh & Nott (2000), who considered passive spheres and predicted a greater increase with φ of effective viscosity than found in this paper, as shown in Fig. 2(a), and our own earlier work on squirmers at lower values of φ (Ishikawa & Pedley 2007b). There have been predictions of suspension rheology in three-dimensional flows, using concentrated suspensions of rigid spheres (Sierou & Brady 2002), and experiments on dilute suspensions of swimming bacteria (Lopez et al. 2015), but even in three dimensions there appear to be no findings on concentrated suspensions of micro-swimmers, apart from Ishikawa & Pedley (2007b). The computational work referred to here was conducted using Stokesian dynamics; analysing the same system using lubrication theory alone has been tried only by Leshansky & Brady (2005), for inert spheres (in three dimensions) -and this was referred to only in a footnote -and Brumley & Pedley (2019), for squirmers, but not in order to predict the rheology. The absence of previous relevant work, however, means that there is plenty of scope for future research. We would like to extend the LT method to three dimensions, where it might be quicker and easier to use than SD, enabling predictions of 3D suspension behaviour for a range of parameter values, building up a graphical representation of a suspension's constitutive relation, as we have tried to do here in 2D. A more complete description of the rheology, however, either in 2D or 3D, will require a quantitative representation of pressure gradients, which are absent in our current simulations since the particle stress tensor, derived from Eq. (1.2), is traceless and does not contribute to the pressure P in Eq. (1.1), though it does determine the normal stress differences. The presence of walls, for a suspension in a channel or pipe, introduces a particle pressure distinct from the average bulk pressure, because of the forces exerted by the particles that interact hydrodynamically with the walls. A clear explanation of this effect is given by Guazzelli & Pouliquen (2018). Singh & Nott (2000) computed normal stresses and pressure in their SD simulations of a suspension of inert spheres in a channel, and we would expect to follow their lead for squirmers.
Other possible extensions to this work include consideration of particles of different shapes, such as prolate spheroids, to represent bacteria (Saintillan 2018), or squirmers that rotate about their axes, to represent Volvox . Either of these extensions would affect the dynamics of solitary squirmers, and hence their trajectories in shear flows and suspensions. At higher areal fractions reminiscent of granular media, shape would influence the propensity of cells to pack together, and control their nematic order. Moreover, lubrication interactions with the bounding walls would also be influenced by organismal shape (Manabe et al. 2020) or rotation.
The current model assumes that all particles are identical and move deterministically. In any real suspension, however, stochastic factors play a part, because of (possibly small) differences between individuals' shapes and locomotory apparatuses. As long as the differences are small, perturbation theory may be feasible, but since the motions of the particles in a suspension become effectively random after a small number of collisions (see Ishikawa & Pedley (2007a)) this is unlikely to be profitable.
Finally, although G bh is designed to mimic gravitactic micro-organisms such as Volvox, the analysis could also be applied to other systems in which an external field tends to align the cells (e.g. phototaxis, magnetotaxis).