## Abstract

Effective viscosity (EV) of suspensions of puller-like microswimmers (pullers), for example *Chlamydamonas* algae, is difficult to measure or simulate for all swimmer concentrations. Although there are good reasons to expect that the EV of pullers is similar to that of passive suspensions, analytical determination of the passive EV for all concentrations remains unsatisfactory. At the same time, the EV of bacterial suspensions is closely linked to collective motion in these systems and is biologically significant. We develop an approach for determining analytical EV estimates at all concentrations for suspensions of pullers as well as for passive suspensions. The proposed methods are based on the ideas of renormalization group (RG) theory and construct the EV formula based on the known asymptotics for small concentrations and near the critical point (i.e. approaching dense packing). For passive suspensions, the method is verified by comparison against known theoretical results. We find that the method performs much better than an earlier RG-based technique. For pullers, the validation is done by comparing them to experiments conducted on *Chlamydamonas* suspensions.

## 1. Introduction

Suspensions of *puller-like microswimmers* such as *Chlamydamonas* algae possess remarkable rheological properties [1,2], but are relatively difficult to investigate experimentally: at meaningful concentrations, it can be challenging to keep the algae from clumping or forming a kind of precipitate. Pullers typically propel themselves by executing a breaststroke-like motion with a pair of flagella attached at the front of the body (e.g. [3]). This should be compared to the propulsion of *pushers* such as *Escherichia coli*, which use a rotating flagellar bundle pushing the fluid behind the bacterial body (e.g. [4]). The rheological behaviour of pushers is rather complicated and at present appears to be beyond the reach of analytical methods except very small concentrations. Pullers, on the other hand, are an important biophysical system that exhibits a range of unique macroscopic behaviours, while being in many ways simpler than pushers. In particular, their rheology is similar to that of passive suspensions. In this paper, we develop a unified analytical approach to approximating the effective viscosity (EV) of suspensions of passive particles and pullers as a function of particle concentration *φ* for *all* physical values of *φ*. Furthermore, our analytical expressions contain only physical parameters, compare favourably with prior analytical approaches and match experimental data for passive suspensions well. Finally, the obtained formulae are in good agreement with experimental data for suspensions of *Chlamydamonas* algae, confirm the similarity of pullers and passive systems, and can be used to probe the rheology of pullers at concentrations inaccessible to experiment, suggest additional experiments and serve as a guide where experimental data are lacking (figure 1).

### 1.1. Effective viscosity

Early experiments on pushers revealed a number of novel macroscopic properties of their suspensions. In particular, in [5,6], pushers were found to exhibit a superdiffusivity for small times. This enhancement in diffusion was attributed to the advection generated by coherent bacterial motion. The same effect was found for pullers, but there it was shown to be much smaller [7].

The effect of different types of propulsion is even more pronounced in EV. For pushers, EV can decrease dramatically with the increasing particle concentrations [8]. For pullers, on the other hand, there is a substantial enhancement in EV for increasing concentrations [1,2]. In experiments on *Chlamydamonas* [1], this enhancement was found to be of smaller magnitude than the corresponding EV decrease for pushers. An increase in the EV also is expected for the model of bottom-heavy *squirmers* oriented by the gravity field [9]. This explanation was proposed for *Chlamydamonas*. However, as detailed in [2,10], it does not seem to apply there.

While there is a relative scarcity of experimental EV measurements for pullers, there are good reasons to expect [1] that pullers behave like passive suspensions with a strongly enhanced EV. Even for passive suspensions, however, the picture of the dependence of EV on the concentration of passive particles is rather incomplete.

### 1.2. Rheology of suspensions

In an isotropic Newtonian fluid, viscosity *η*_{0} is a scalar quantifying the rheology of the fluid—the relation between the macroscopic strain rate *e _{ij}* and stress

*σ*in the fluid. In a suspension of microscopic particles, both the strain rate and the stress are modified by the disturbance flow around the particles, and the EV

_{ij}*η*quantifies the relationship between these volume-averaged quantities. While in general the presence of particles requires a tensor to specify the rheology, operationally we can define a scalar

*η*as relating fixed components of

*e*and

*σ*(e.g. the

*xy*-components for particles suspended in a planar shearing flow). Given the low Reynolds numbers typical of suspensions, one has to solve the Stokes equation in the exterior of the microscopic particles, which amounts to balancing the forces and torques on the fluid. In the simplest case of dilute suspensions, the inclusions do not interact and are statistically independent. In this situation, the disturbance flow of a single typical particle is sufficient to determine the EV, and the normalized correction is linear in the volume fraction of inclusions

*φ*. For prolate spheroids with the Bretherton constant

*B*=

_{r}*b*

^{2}−

*a*

^{2}

*/b*

^{2}+

*a*

^{2}quantifying their aspect ratio of the particles in terms of its minor and major semi-axes

*a*and

*b*, the dynamics of particle orientations

*τ*in a linear ambient flow was obtained by Jeffery [11]. When the particles are subjected to the Brownian motion, a unique long-time distribution of orientations can be obtained, predicting an increase in EV. For self-propelled particles modelled as prolate spheroids with a rigidly attached point force (figure 2), the situation is different. Here, the additional propulsion force must be incorporated into the balance of forces and torques acting on the fluid and the particle, when determining the disturbance flow and the propulsion velocity. Elongated particles, in accordance with Jeffery's equations, will spend more time aligned to the stable principal axis of the ambient planar shear flow. In [4], we developed asymptotic expressions for the EV of active suspensions in various limiting cases, including small asphericity as well as diffusion-dominated particle dynamics.

### 1.3. Beyond the dilute limit

The theory developed in [4] produces the leading order expansion in *φ*. Higher order coefficients reflect interactions between the suspended particles, accounting for the influence of pairs, triples, and higher interactions, on the disturbance flow. At the other end of the spectrum lies the behaviour in the densely packed limit—at critical *φ*_{max} the EV can be expected to diverge. We seek to elucidate the regime in between these two extremes.

Initially, we consider the case of passive suspensions. Even here the knowledge of asymptotic expansions is far from complete. The expansion in the particle volume fraction *φ* is available up to the second-order term inclusively, with the result for the EV of spheres owing to Einstein [12]. Estimates of the second-order term in the low-concentration asymptotic expansion for passive suspensions go back to Batchelor & Green [13], where the weak Brownian limit was considered.

Beyond the second-order low-concentration expansion, various analytical expressions have been suggested to capture the enormous volume of experimental material for passive suspensions in the broad range of concentrations. Rutgers [14] studied hundreds of different expressions for the EV and concluded that those with a critical point compared with experiment are much better than others. Such approaches assume the existence of a power-law divergence of the viscosity [15] in the vicinity of the maximum volume fraction value [16]. The broad consensus of what *φ*_{max} should be emerged only recently and is now believed to reflect the density of *maximally random jammed state* [17] not the dense packing of spheres, which is virtually unattainable in a real suspension.

As far as the intermediate range is concerned, it is accessible only via various empirical or analytical approximation techniques. We propose a systematic approach for constructing analytical estimates of the EV in the intermediate range of concentrations that match the known asymptotic behaviour for small concentrations as well as near the critical point. The method is based on the ideas of renormalization group (RG) theory. To fully use the information contained in the asymptotic expansions of EV, we *renormalize* them or ‘restore’ the expression for arbitrary intermediate concentrations based on the information ‘hidden’ within the limiting expansions.

## 2. Algebraic renormalization method

The general idea behind the RG is that it uncovers a hidden symmetry, which is used to reconstruct an unknown function in the vicinity of its critical points [18]. Essentially, the same idea can be applied to the problem of functional reconstruction of the EV, not only at the critical point, but also in a broader region of concentrations. Previously such an approach was applied to the determination of the EV of a hard-sphere suspension in an important publication [19], which followed the general idea of Shirkov [20]. The realization of the idea has remained incomplete until now: there is no final transparent analytical result and no clear physical interpretation of the RG in the context of EV reconstruction. Additionally, numerical results obtained in Altenberger & Dahler [19] are not accurate enough (figure 3).

The main shortcoming of the method employed in Altenberger & Dahler [19] is that the authors only used a perturbative expansion and excluded critical point behaviour. Here, we develop a variant of the RG method, which allows for the natural inclusion of the behaviour in the vicinity of the critical point, making it different from the RG of Altenberger & Dahler [19].

All of the derivations done in this section are presented in a general form, motivated by the specific case of the EV of suspensions. Consider a function *f* defined on interval [0, *x _{c}*). By definition, in the vicinity of the critical point,

*f*has the following asymptotic behaviour as

*x*→

*x*: 2.1where

_{c}*B*is the

*critical amplitude*. The

*critical exponent β*is the limit

*x*→

*x*of the

_{c}*effective critical exponent β*(

*x*) [21] defined in a neighbourhood of

*x*2.2Thus, near the critical point,

_{c}*f*has a divergent behaviour with

*β*< 0. At the other end of the interval,

*f*has a polynomial perturbative expansion, with the

*k-*th order approximation given by the series 2.3with integer

*k*≥ 1. The

*k*-th order effective exponent can be defined here analogously by 2.4which for

*k*= 1,2 has explicit expressions 2.5and 2.6The essential idea is to extend

*β*(

*x*) to the whole interval [0,

*x*) so as to connect the polynomial behaviour (2.3) of

_{c}*f*near

*x*= 0, to the power-law behaviour (2.1) near

*x*=

*x*. According to (2.2), this is equivalent to reconstructing

_{c}*f*and can be viewed as an appropriate resummation of (2.3) in the whole interval.

There are several reasons for the introduction of the *β*-function in addition to *f*, which leads to a seemingly redundant differentiation, only to undo it by integration later. The physical reason lies in the use of the techniques and the tradition of critical phenomena, where *β* has a more transparent meaning as the critical index. The mathematical reason lies in the fact that *β*, although equivalent to *f*, satisfies simpler identities, which are frequently obscured by integration. In particular, *β* has the constant asymptotic at the critical point *β* → *−t*.

In order to achieve the resummation, *f* has to be *renormalized* to so that each subsequent is obtained from a lower order expansion via an appropriate *similarity* transformation—an element of RG. On approaching the transformation's fixed point, convergence of is achieved. This idea leads to the method of *self-similar approximants*, based on the so-called algebraic self-similar RG [22] and self-similar bootstrap [23]. Applying this method, one can derive a self-similar approximant for the effective critical exponent denoted , which possesses the correct limit
2.7

Within the context of this paper, the value of *β* is *known* and finite. Finally, from the knowledge of one can approximate the sought function, by integrating (2.2) to find
2.8

The critical amplitude *B* can be evaluated from (2.8) as
2.9

It is important to note that the formal method of self-similar approximants may in general lead to multiple solutions to *β**, because the formally constructed fixed point depends on unknown exponents and constant factors, which parametrize the most general expression for . With sufficiently many terms in the expansion, it may be possible to determine all of these parameters uniquely by selecting them so that all the *a _{k}* as well as the critical exponent are reproduced. When only low-order expansions are at our disposal, additional choices have to be made to select appropriate self-similar approximants. Various choices have been proposed—detailed in the literature cited—as we encounter the specific self-similar approximants. The most frequently used here are what we call the root, iterated root, factor and super-exponential approximants. These approximants can be obtained systematically from the general self-similar fixed-point expression discussed above. Following the work [22] that introduced the original version of this approach, we refer to it as the algebraic renormalization method (ARM).

There exists a principal problem in extrapolation/interpolation of functional dependencies from low-order asymptotic expansions. The problem is in the reliability or stability of the obtained expressions with respect to the variations in the coefficients of asymptotic expansions. It is important, therefore, to be able to make the extrapolation by several methods, and then compare the results. If these results are close, this suggests that the extrapolation is reliable. In line with this idea, we employ various self-similar approximants from the ARM framework and show that their results are close to each other, indicating that the basic machinery of ARM is robust.

## 3. Passive suspensions

We first apply the ARM method to the case of passive suspensions of spherical particles. Here, we specialize the formulae of the previous section to *f* = *η*, *x* = *φ*, *x _{c}* =

*φ*

_{max}, with the critical index

*t*=−

*β*. The expansion in the particle volume fraction

*φ*is available up to the second-order term inclusively as

*φ*→ 0: 3.1

Einstein determined that for hard spheres *a*_{1} = 5/2 (see [12]). For the value of the second-order coefficient, several *estimates* are available. In particular, Batchelor [24] gives *a*_{2} = 6.2, and Altenberger and Dahler [19] give *a*_{2} = 5.9148, both results including the hydrodynamic and Brownian contributions. We will use Batchelor's estimate unless stated otherwise.

Furthermore, there is a critical point characterized by a power-law divergence of the viscosity in the vicinity of the maximal volume fraction value 3.2

This corresponds to the *random closest packing* of hard spheres at [16]. The value of the critical exponent has been explained both theoretically and experimentally [15,25].

To test our results, we compare them to the analytical expressions for the EV assembled in appendix A. The most often used is the well-known Krieger–Dougherty (KD) equation (A 1), an example of a crossover expression that links polynomial and power-law behaviours at opposite ends of the concentration interval [26]. The empirical Thomas equation (A 2) is probably the most accurate fit *outside* of the region close to *φ*_{max} (see [27]). This result will be considered the main benchmark against which to measure the accuracy of all analytical formulae away from the asymptotic regimes. We present also a very transparent semi-empirical analytical expression for the EV (A 3) suggested by Brady [15]. Finally, the earlier numerical RG-based result (A 4) from Altenberger & Dahler [19] is discussed.

### 3.1. The first-order interpolation

We have constructed EV expressions for 3*D* passive suspensions of spheres from the first-order coefficient *a*_{1} as well as the critical density and exponent values. These analytical expressions were then used to obtain the second-order coefficient *a*_{2} which will be compared against a benchmark, establishing an accuracy test for *φ* close to 0 and *φ*_{max}, respectively.

#### 3.1.1. Transformation of independent variable

Using a change of variables *z* = *φ*_{max} − *φ* and the simplest root approximant [28], one can easily obtain the simplest expression for the effective critical index from the first-order expansion
3.3

This is a linear expression, which after integration of (2.8) with yields the following formula: 3.4which is probably our most important result for passive suspensions. Asymptotically close to the critical point we recover a power-law, while at small concentrations the exponential term dominates.

#### 3.1.2. Iterated root

One can approach the problem of reconstruction of the effective critical index differently: we look for the solution in the form of an iterated root. Such approximants were designed to find all the parameters iteratively [29] 3.5in such a way that . After integration in (2.8) with , we have the explicit solution 3.6

#### 3.1.3. Exponential approximant

Another approach to interpolation is to look for the solution for the effective critical index in the form of a regular function represented by the simplest (super) exponential approximant [23]
3.7where *τ* is obtained from the condition . In this case, integral (2.8) can be written using exponential integral :
3.8Directly from the expressions (3.4), (3.6) and (3.8), one can obtain the second-order coefficient in the expansion and it equals, respectively, *a*_{2} = 5.566, *a*_{2} = 5.469 and *a*_{2} = 5.514.

We stress that the first-order interpolation schemes are particularly relevant for our study of active systems with only the first-order expressions available.

### 3.2. The second-order interpolation

In the expressions above, *a*_{2} was not used to construct a more accurate interpolation, which we do here. In order to check these second-order interpolants, we have to compare our predictions of *a*_{3} against those obtained from the benchmark formulae, because there is a relative paucity of numerical results in the broader range of concentrations. For future comparison, we also note that the critical amplitude *B* may be estimated as *B* ∼ 0.41 following [15].

#### 3.2.1. Transformation of independent variable

Let us again employ the variable *z* and the second-order iterated root approximant [29], so that one can obtain the following expression for the effective critical index using (3.1):
3.9where all parameters can be easily obtained explicitly from the asymptotic conditions. The EV is then represented by integral (2.8), enabling estimation of the third-order coefficient, *a*_{3} = 12.563 and the critical amplitude *B* = 0.342.

#### 3.2.2. Iterated root

We look for the solution in the form of a higher order iterated root [29] 3.10

Imposing the condition on ‘critical’ point together with usual asymptotic conditions, one can obtain all parameters in (3.5) explicitly. Again, the EV is given by integral (2.8), and one can find an estimate for the third-order coefficient *a*_{3} = 13.382 and for the critical amplitude *B* = 0.398.

#### 3.2.3. Super-exponential and root approximants

Extending the first-order derivation, one can construct the following super-exponential approximant [23]:
3.11where *τ*_{1} and *τ*_{2} can be obtained explicitly from the asymptotic conditions. Correspondingly, follows from (2.8), from which one can find an estimate for the third-order coefficient *a*_{3} = 12.726 and for the critical amplitude *B* = 0.34.

Instead of the super-exponential approximant, one can use a root approximant [28] leading to another expression for the effective index
3.12where *c* can be obtained numerically, *c* = 0.134 and . As usual, the EV is given by integral (2.8). One can then find an estimate for the third-order coefficient *a*_{3} = 11.326 and for the critical amplitude *B* = 0.309.

#### 3.2.4. Averaged predictions

A different approach to the second-order problem could be based on pairwise weighted averages (see [30]), of (3.4), (3.6) and (3.8), where a total of three such pairs are available. For example,
3.13where the weight *p*_{1} is determined uniquely through the asymptotic condition on the second-order coefficient *a*_{2}. and are defined by analogy. Such averages have some advantage over the other second-order interpolations, because they express the EV analytically.

#### 3.2.5. Factor approximant

To achieve further improvement in the estimates of *a*_{2}, based solely on the given *a*_{1}, *t* and *φ*_{max}, one would clearly need some additional information about the system. A non-trivial change of variables can carry such additional physical information about the EV. It is demonstrated in [31] that the next-order correction to the Einstein formula for a suspension is at least of order *φ*^{3/2}. In order to incorporate this information into our study, we apply the self-similar approximants directly to the expressions for the EV in the dilute limit and close to *φ*_{max}. This will satisfy the limits and also allow the terms of order *φ*^{3/2} to appear in the expansions at least implicitly. As we do not know the amplitude for such a term, we only require it to be very small. The so-called factor approximants [32] are uniquely able to address this issue and supply the required transformation, *z* = (*φ*/*φ*_{max})^{1/2}*/*(1 − (*φ*/*φ*_{max})^{1/2}).

We construct the factor approximant *η ^{f}*(

*φ*) and both factors within it to balance each other and suppress all terms appearing after the re-expansion of the order of

*φ*

^{1/2}and

*φ*

^{3/2}, 3.14where all parameters can be uniquely determined from the asymptotic conditions. From (3.14), is a function of the critical concentration

*φ*

_{max}. Evaluating at the random close packing fraction

*φ*

_{max}= 0.64 yields .

The utility of these approximation can be evaluated by comparing them to the (semi-)empirical benchmarks (A 2) and (A 3). In particular, (A 2) is accurate outside of the asymptotic regimes *φ* → 0 and *φ* → *φ*_{max}, and our predictions are very close to this formula away from the asymptotics. Formula (A 3) has a much wider applicability and our predictions, albeit lower than those of (A 3) are in agreement with this benchmark.

Furthermore, comparing the widely used analytical predicition (A 1) to the same (semi-)empirical benchmarks we see that the ARM results fare quite a bit better. Similarly, our predictions outperform the earlier results of the RG method (A 4). In both cases, the advantage seems to be traceable to a more accurate and systematic incorporation of the critical behaviour into the ARM formulae.

## 4. Active suspensions: pullers

Having established the utility of the ARM for passive suspensions, we now apply it to suspensions of swimming bacteria. Specifically, we seek to predict the concentration dependence of EV for suspensions of *Chlamydamonas* algae, a puller organism discussed in the Introduction. The important characteristic of pullers is that their presence tends to increase monotonically the EV of the suspension, similar to passive suspensions. For pushers, the EV concentration dependence is not monotonic and exhibits a lowering of the EV for low to moderate concentrations.

Naively, one would expect all swimming cells or bacteria to affect the EV of the suspension similarly: after all, all swimmers inject energy into the fluid, decreasing the effective rate of viscous dissipation and leading to an effective decrease in viscosity. Therefore, it may still come as a surprise that the effect on the EV is opposite for pushers and pullers [1,8]. The usual qualitative explanation is as follows: in an imposed shear flow, passive elongated particles approximated by infinitesimal rods—force dipoles—tend to align with the extensional axis of the flow and the EV is increased. When a puller, approximated by a contractile dipole, is aligned with the extensional axis, the rate of strain is decreased. Recalling that viscosity is inversely proportional to the rate of strain, it can be concluded that under a constant external force the EV should increase.

Thus, it is plausible to assume that pullers behave similar to passive particles in their overall effect on the EV. This is the conclusion of experiments in [1], based on which we set out to reconstruct the concentration dependence of pullers, compare this prediction to the experimental results in [1] and offer a plausible explanation for the experimentally observed discrepancy in the magnitude of the effect of activity on the EV between pullers and pushers.

We apply the ARM to reconstructing EV's concentration dependence from the (effective) *intrinsic viscosity*, adopting this from [2] as another name for the first-order perturbative expansion coefficient
4.1

For comparison with experiment, *a*_{1} has to be calculated as a function of the puller's internal parameters and of the flow. It is then supplemented with the critical behaviour at the dense packing limit, which we assume to be similar to that of passive suspensions. With these two pieces of asymptotic information at hand, we apply the ARM formulae to establish the concentration dependence of EV for pullers.

### 4.1. Estimates of intrinsic viscosity

To reconstruct the dependence of the intrinsic viscosity on the parameters of the background flow and the shape of the swimmers, we heavily rely on [4]. The paper is dedicated to pushers, modelled as prolate spheroids with semi-axes *a*, *b*, but in [4], it was shown how the case of pullers can be obtained from the same formulae by changing the rigidly attached propulsion force *f*_{p}. All relevant expressions for the case of pushers for planar extensional flow can be found in [4], specifically see (45) and (46) of that paper.

Using those results, we intend to obtain analytical expressions for the intrinsic viscosity as a function of the arbitrary strain rate *γ* of the background flow as well as the Bretherton constant (§1.2 and figure 2). For the small parameter *μ* = *γB _{r}/D*, applicable when the background flow is much smaller than the rotational diffusion of pullers (weak flow regime), the following asymptotic is valid:
4.2with
where

*D*is a constant specifying the strength of rotational diffusion also called tumbling. For the large

*μ*, when the background flow dominates the rotational diffusion (strong flow regime), a different asymptotic is valid: 4.3coefficients

*M*

_{0},

*M*

_{2},

*K*,

*S*and

*N*are the so-called

*shape factors*, which depend only on the geometry of the puller [4]. As in the strong flow regime, the particles strongly tend to align along the extensional axis of the flow (whenever the background flow is purely straining), we refer to this regime as the

*aligned limit*.

#### 4.1.1. Intrinsic viscosity of pullers: weighted approximant

The experimental data from [1] lie outside the asymptotic regimes of (4.2) and (4.3). In order to estimate the intrinsic viscosity for the suitable value of *μ*, we interpolate between these asymptotics by the following physically motivated method. It is natural to separate the passive and active contributions to the effective intrinsic viscosity. Such a separation is justified intuitively by the presence of passive and active properties of the bacteria. Therefore, we start with the construction of the simplest suitable approximants:
4.4which is the simplest two-point Padé approximant, and the simplest exponential approximant:
4.5Both are designed to match the behaviour of the second-order truncation of expansion (4.2) (*α*_{2} = 0) and the aligned limit (4.3).

We next form a weighted average or a mixture of the two approximants
4.6using, remaining condition—matching *α*_{2}—to determine *p*, following [30]:
4.7where *u*_{1} and *u*_{2} are the second-order terms in the expansion of and in the parameter *μ*, respectively. Finally, expressing the dependence of *α*_{1}, *L* and *μ* on *γ* obtains the intrinsic viscosity as a function of the strain rate.

At this point, the weighted approximant applies to pushers and pullers depending on the sign of *f*_{p}. As discussed above, approximant separates the two contributions to the effective intrinsic viscosity. It turns out that the contribution from the exponential term outweighs in magnitude the contribution from the Padé term, both for pushers and pullers.

A similar curve for *f*_{p} of opposite sign is plotted for pullers (figure 4). The effective intrinsic viscosity remains, in this case, always positive and approaches the aligned limit from above and another limit of small *γ* from below.

### 4.2. Concentration dependence of puller effective viscosity

The asymptotic expressions in [4] were derived in the absence of interactions between particles and yield only the intrinsic viscosity *a*_{1}. Obtaining higher order coefficients *a _{k} k* > 1 does not appear feasible, because incorporating pairwise particle interactions into the method of [4] seems to be intractable. Instead, we apply the ARM to the first-order expansion in concentration with .

To complete the setting for the application of the ARM, the critical point behaviour is needed. At least formally, the densely packed regime with power-law divergence can be invoked here. This uses the fact that for large concentrations *Chlamydomonas* algae cells become immobilized and nearly spherical [1], apparently falling dormant in the untenable conditions of very high concentrations. With this justification, we can apply the factor approximant *η ^{f}* from (3.14) using as the first-order coefficient.

The approximant is a function of both *φ* and *μ*, describing a crossover from the *active* system at small *φ* to the *passive* system of dead, immobilized, near-spherical cells near *φ*_{max} at all values of flow strength *μ*. Various approximants can be constructed for pullers with a positive EV but we will rely more on the approximants already tested above in the passive case. The factor self-similar approximant *η ^{f}*(

*φ*) turned out to be the most consistent for the passive systems. It will also provide a

*lower bound*on our estimates for the EV of pullers. The same expression, but with

*a*

_{1}‘renormalized’ (i.e. replaced by for arbitrary

*γ*), will be applied for pullers. The final expression applies for arbitrary

*φ*and

*γ*.

To assess the quality of the resulting prediction, we calculate the predicted enhancement in the EV relative to the passive spherical case where *a*_{1} is independent of the background flow strength *γ*. Motivated by [1], this enhancement is quantified by the ratio and can be compared to experimental results of [1,2].

#### 4.2.1. Comparison with experiment

To calculate the enhancement ratio, we use the following parameters of the puller model (figure 2): the major semi-axis of the approximating spheroid *b* = 5 μm, force location characterized by *λ* = 0.5, bacteria's swimming speed *v* = 100 μm s^{−1} and rotational diffusion *D* = 0.4 s^{−}^{1} taken from [1,3]. Finally, the Stokes drag law holds: *f*_{p} = 6*π*b*η*_{0}*vN*, where an expression for *N* can be found in [4].

For *γ* = 5 s^{−}^{1}, *φ* ≈ 0.15 and *b*/*a* = 2.5 (motivated by an estimate of flagella length, 10−12 μm, see [33]), the enhancement ratio is around 1.757. The first-order renormalized coefficient is estimated as 7.122, and the second-order coefficient , rather different from their counterparts for passive systems. To make our estimates more robust, i.e. less sensitive to the approximate nature of interpolation formulae, we calculated a similar enhancement ratio for other than factor approximants, constructed and tested for passive suspensions and discussed at length above, with *a*_{1} replaced by . One can relatively easily construct the KD approximant (A 1) which provides an *upper bound* for the EV. Finally, we construct EV with all three ARM approximants obtained in closed form by transformation of variables (3.4), iterated root (3.6) and exponential approximation (3.8) applied to their corresponding *β*-functions.

The enhancement ratio calculated and then averaged over all five models based on the first-order expansion in concentration (factor, three ARMs and KD), becomes larger and equals 1.956. The ratio also depends weakly on *γ* in agreement with the results on shear-thinning from [10]. On the other hand, if *D* is decreased to 0.1 in accordance with [2], then the average enhancement ratio increases to 2.5 which is in good agreement with [1], and also could become considerably larger for smaller *γ* (for *γ* = 1.3 s^{−}^{1,} we estimate an increase of seven times the ratio). Such an increase in the EV for pullers matches the corresponding decrease in the EV for pushers [4].

The parameter values required for an explanation of the experiment cannot be explained on a purely geometrical basis as shown in [2], but rather point to their effective nature. These quantities are probably generated by the flagellar shape and beating movements, in accordance with [2]. A refinement of Haines *et al*. [4] would include such a description of the flagellar dynamics. In practice, it is sufficient and most convenient to use the most transparent ARM approximant given by (3.4), with coefficients from (3.3). The only difference is that *a*_{1} is replaced by from (4.6), so that it becomes . The enhancement ratio for this approximant becomes particularly transparent
4.8

## 5. Discussion

The EV formula is reconstructed by ARM from two asymptotic expressions by treating EV as a peculiar critical phenomenon. Why does the ARM work?

(1) The first-order coefficient

*a*_{1}, contains the information about the interaction of the inclusion with the background flow alone, while dilute binary interactions are encoded in the second-order coefficient*a*_{2}.(2) The asymptotic behaviour near the random dense packing limit encodes the many-particle interactions.

(3) The geometry of the problem enters through the value of random closest packing of hard spheres, which characterizes the typical particle arrangements in the dense limit.

Therefore, any higher terms in the concentration expansion carry little information not already contained in (1) or (2). As long as no other critical phenomena arise in the intermediate range of concentrations, the two asymptotic expansions can be used to effectively reconstruct EV in the full range of concentrations.

Comparison with experiment and other analytical results is rather encouraging.

For pushers, however, collective modes arising at intermediate concentrations result in a phase transition well before the packing transition, and additional information about the asymptotic behaviour near the coherence transition is necessary to recover the EV in the full range of concentrations. In particular, active suspensions of pushers exhibit at low concentrations a non-monotonic behaviour of the EV making them qualitatively different from the passive case.

Pullers, on the other hand, rheologically behave as renormalized passive particles of a certain larger effective size. Perhaps, this can be attributed to the fact that pullers activity, being the opposite of that of pushers, enhances the EV and is similar to the passive effect. Thus, there is no additional transition at intermediate concentrations in the EV as a function of concentration.

In case of passive suspensions, we suggested a new formula (3.4) for the EV of suspensions for all concentrations, using only Einstein's coefficient as input from the low-concentration regime and the values of the critical exponent and critical density from the densely packed regime. This formula is more accurate than the prior renormalization-based result of Altenberger & Dahler [19]. This is primarily owing to the careful incorporation of the asymptotic information from the high-concentration regime as well as using a carefully chosen variable transformation.

Concerning suspensions of puller-like microswimmers, our primary goal was to derive working analytical expressions, most important see (4.8), for their EV. This task is particularly important in light of the considerable difficulty of experimental investigation of pullers. We combined the microscopic derivation of the first-order asymptotic expansions of EV for pullers [4] with the ARM developed for and tested on passive suspensions. The result was a single formula that relates the EV of active systems of pullers with the background rates of strain and particle concentrations. The main prediction of [4], the increase in the EV for pullers, now becomes quantitative and in good agreement with the experiment of [1].

Finally, using our quantitative EV results for pullers, we propose an explanation for the relatively modest experimental enhancement of puller EV [1] when compared with the more pronounced EV reduction for pushers [8]. We show that this discrepancy can be remedied by adjusting experimental parameters. We can obtain, at least in principle, an enhancement similar to the reduction for pushers. For the typical values of *D*, quantifying rotational diffusion owing to tumbling or other random effects, EV has a very weak dependence on strain rate *γ*. By decreasing *D*, the EV dependence on *γ* becomes stronger, predicting a more pronounced shear-thinning effect. Thus, by lowering *γ* at the simultaneously reduced *D* a significant enhancement of the EV is obtained. Thus, the relatively modest EV enhancement observed in [1] can be explained by the relatively high rotational diffusion and background strain rates used in the experiment. This suggests additional, albeit non-trivial, experiments that can test this hypothesis.

## Acknowledgements

The work of all three authors was supported by NIH/NIGMS R01 grant no. 1R01GM104978-01. The authors express their gratitude to I. Aronson for his explanations of puller experiments, to P. Peyla for useful discussions and to S. Ryan for multiple useful suggestions.

## Appendix A. Reference formulae for passive suspensions

### A.1. Krieger–Dougherty equation

The most often used is the well-known KD equation [26] A1

It links the polynomial and power-law behaviours at the opposite ends of the concentration interval. It can be used to evaluate the second-order coefficient and the result is 5.078.

In table 1, the numerical values for the second-order expansion coefficient obtained by different approaches are compared.

### A.2. Phenomenological formula

The empirical Thomas equation [27] is probably the most accurate fit outside the neighbourhood of *φ*_{max}
A2

Thomas excluded many non-universal experimental features and presented by means of (A 2) an average picture of the EV based on rich experimental data. This result can be considered a benchmark against which to measure the accuracy of all analytic formulae. Another important such test involves evaluating the third-order coefficient *a*_{3} resulting from an analytical expression obtained from the second-order expansion. In this situation, the empirical expression (A 2) is not very useful.

### A.3. Brady's formula

In [25], Brady suggested an explicit semi-empirical expression for the EV (see also [34]): A3

Several physical parameters, which are functions of *φ*, are required: the short-time self-diffusion coefficient of suspended particles and the high-frequency limiting viscosity *η*_{H}. The radial distribution function *G*(2, *φ*) at particle contact, was obtained for *φ* < 0.5 from the Carnahan–Starling equation of state for rigid spheres. For higher concentrations, *φ* ≥ 0.5 the formally obtained number of particles ‘at contact’ has to be obtained from numerical simulations, hence *G*(2, *φ*) = 0.78/(0.64 − *φ*), and it diverges as random close packing is approached. Using the expressions above, we are able to make estimates *a*_{2} = 7.4 and *B* = 0.5, which will be used in comparisons against our results below.

Brady's formula implies that the critical index *t* equals 2. This is non-trivial and encompasses two critical phenomena—the divergence of *G*(2, *φ*) with an exponent of 1 and a similar divergence of (inverse) . The critical exponents add to give the value of *t* = 2.

### A.4. RG formula

The earlier numerical RG results from Altenberger & Dahler [19] can be summarized as follows:
A4with given explicitly [19]
and *ψ*_{0} = 0.085 is chosen to minimize the difference between the RG solution (A 4) and expansion (3.1).

### A.5. Comparison with references

From the analysis of experimental data in [35], one can only conclude that the *a*_{3} is in the range of 10–16. When only three-particle hydrodynamic interactions are taken into account the second- and third-order coefficients are known theoretically, namely *a*_{2} = 5 and *a*_{3} = 9.09 [36]. With such a value of *a*_{2}, using the ARM approximant corresponding to the *β*-function (3.10), we obtain a rather reasonable estimate, *a*_{3} = 9.53. Coefficient *a*_{3} can be also estimated from Brady's expression (A 3). Using the input parameters we are able to estimate the second-order coefficient in the expansion for small concentrations as *a*_{2} = 7.4 and the third-order coefficient *a*_{3} = 17.42, while the critical amplitude estimate is *B* = 0.5. In order to compare this against our results, we use the theoretical value of *a*_{2} = 7.6 from [13] obtained in the so-called weak Brownian limit. Averaging then over all second-order approximants (except the root approximant) we obtain *a*_{3} = 15.81 and *B* = 0.46, in reasonable agreement with Brady.

In table 2, we compare the numerical values for the third-order expansion coefficient obtained by different approaches.

We next turn to pairwise averages of first-order interpolants: compares a little better than other ‘pairs’ with Brady's results for EV, while the corresponding results from the RG of Altenberger & Dahler [19] are inferior. For the third-order coefficient in the expansion is evaluated as *a*_{3} = 17.079 and the critical amplitude is *B* = 0.558. For the typical large value of *φ* = 0.4, gives values lower than (A 3), deviating from it by 8%, while (A 4) is lower by 19.8%. Our simple formula is also in good agreement with [34], which also systematically gives results lower than Brady's formula (A 3).

In the critical region, as *φ* approaches *φ*_{max}, is only lower than Brady's formula by 2.8% (*φ* = 0.6), while (A 4) is lower by 48%.

To establish a comparison of our results with other analytical predictions in the broader range of concentrations, we compare against the main experimental benchmark (A 2). One of the most interesting comparisons is with the RG formula (A 4), an earlier representative of the renormalization approach.

Specifically, we compute the deviation *ɛ* from (A 2) for the typical large value *φ* = 0.4, as predicted by the following ARM results: (i) formula (3.4), (ii) the results of numerical integration for the *β*^{*}-function (3.12) and (iii) factor approximant (3.14).

The deviation or percentage error *ɛ*, is a conventional measure defined as the difference of our ‘prediction’ with the ‘exact’ result, divided by the exact result and measured in percentages. We had chosen the typical non-small value for concentration where Thomas equation is valid. It is also well outside of the dilute regime and of the critical region, to see how well ARM predicts/interpolates.

We find that the ARM results are all quite accurate, with the corresponding errors for (i) *ɛ* = −2.78%, for (ii) *ɛ* = −2.92% and for (iii) *ɛ* = 5.4%.

In comparison, for the widely used KD formula (A 1), the error is larger, *ɛ* = −15.68%. The largest error is found for the RG formula (A 4) *ɛ* = −22.97%, but is still much better than that from the second-order expansion (3.1), which gives *ɛ* = −47.48%.

In figure 3, we present the two typical curves, (3.4) and the results of numerical integration for the *β*^{*}-function (3.12), both given by our ARM. They are compared with Thomas empirical formula (A 2), Altenberger–Dahler RG (A 4), and the expansion (3.1). Clearly, the ARM formulae bring a substantial improvement on the important original RG method.

Only the best first- and second-order ARM approximants are used in these comparisons. We also calculated the error for all other ARM approximants and found it of the same order.

- Received August 5, 2013.
- Accepted September 3, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.