Asymptotics and Summation of the Effective Properties of Suspensions, Simple Liquids and Composites

Abstract

We review the problem of summation for a very short truncation of a power series by means of special resummation techniques inspired by the field-theoretical renormalization group. Effective viscosity (EV) of active and passive suspensions is studied by means of a special algebraic renormalization approach applied to the first and second-order expansions in volume fractions of particles. EV of the 2D and 3D passive suspensions is analysed by means of various self-similar approximants such as iterated roots, exponential approximants, super-exponential approximants and root approximants. General formulae for all concentrations are derived. A brief introduction to the rheology of micro-swimmers is given. Microscopic expressions for the intrinsic viscosity of the active system of puller-like microswimmers are obtained. Special attention is given to the problem of the calculation of the critical indices and amplitudes of the EV and to the sedimentation rate in the vicinity of known critical points. Critical indices are calculated from the short truncation by means of minimal difference and minimal derivative conditions on the fixed points imposed directly on the critical properties. Accurate expressions are presented for the non-local diffusion coefficient of a simple liquid in the vicinity of a critical point. Extensions and corrections to the celebrated Kawasaki formula are discussed. We also discuss the effective conductivity for the classical analog of graphene and calculate the effective critical index for superconductivity dependent on the concentration of vacancies. Finally, we discuss the effective conductivity of a random 3D composite and calculate the superconductivity critical index of a random 3D composite.

1. Introduction

Real physical systems can almost never be studied in terms of exact solutions to their governing dynamics. However, it is often possible to find asymptotic expansions of solutions in the vicinity of some characteristics to the problem boundary points. The ultimate goal is to construct a good approximation for the sought function, representing the solution to the problem, valid on the whole domain of its variable(s), having as an input only its asymptotic behavior near the boundaries. A typical complication is associated with availability most often, only of a few terms of the asymptotic expansions.

Sometimes, we have to derive asymptotic power laws that are ubiquitous in nature. Sometimes, following Feynman, one would attempt to estimate in a given order the result for the expansion coefficient(s) without the brute force evaluation. Sometimes the general structure of the expansion can only be guessed. Many instances of such problems were discussed in the books [1,2,3,4]. Most recently, the disordered hyperuniform composites were reviewed in [5].

The transport coefficients in random and regular media could be expressed as expansions in a concentration of particles f. Despite persistent efforts of such outstanding researchers as Einstein, Batchelor, Brady, Jeffrey, Milton, McPhedran, Torquato, Wajnryb and Mityushev, most recently, the problem still exists of finding long enough and correct expansions. When it is only possible to find short truncated series, their true value may still need clarification. Such expansions are relevant for small f.

The transport coefficients such as effective viscosity, permeability, sedimentation rate, elastic moduli, effective diffusion coefficients and conductivity behave critically when approaching critical points in concentration. The relation between the two types of expressions is difficult to establish. However, it has to be looked for if the ultimate goal of a complete understanding of the whole region of concentrations is to be achieved.

The problems of finding critical indices for the effective transport properties are decidedly non-universal, meaning that there is no simple formula and (or) explanation of the value. This is in sharp contrast with the theory of critical phenomena in thermodynamics, where the critical indices are universal [6], meaning that the formulas exist, covering different physical systems and explaining the indices in terms of spatial dimension and number of components of the order parameter. Thus, for transport properties, the challenge consists in developing the solution in each and every case separately, almost from scratch, when existing methods fail.

However, certain universality exists in the methodology. In all cases, we have a short expansion in concentration f of particles (or close in meaning) for small f. On the other side of high concentrations (or equivalent), we are aware of the power-law behavior. In some cases, the critical properties have to be calculated from the expansions on the other end. In some cases, the two asymptotic tails are known, and we have to merge them into a single formula. The location of criticality is also known from geometrical considerations most often, or the power-laws can be attributed to infinity. Another common salient feature is the shortness of the expansions and ensuing paucity of coefficients in the expansions.

To approach the problem, we had to master or even develop some special tools for the two closely related problems, such as

(1)

The derivation of critical properties, such as indices and amplitudes from the short truncated versions of the putative complete expansions for the effective quantities, or critical problems;

(2)

Approximation methods for deriving formulas from the sought quantities from their known asymptotic expressions, or crossover problems.

Thus, the main goal is to attack the group of problems, including critical behavior of the effective transport properties and develop special methods for their solution. All the systems and critical problems considered in the review are being studied extensively and are important on their own merits. However, putting them together and analyzing them by similar methods is intended to stress their commonality. Of course, the commonality is due to the interactions between inclusions taken into account in low orders, as well as criticality arising at some common geometrical threshold.

In Section 2 and Section 3, we discuss the effective viscosity (EV) of random, active and passive suspensions by means of a special algebraic renormalization approach. It is applied to the first and second-order expansions in volume fractions of particles suggested by Einstein, Batchelor et al. The EV of the two-dimensional and three-dimensional passive suspensions is analyzed by means of various self-similar approximants.

A brief introduction to the rheology of micro-swimmers is given. Microscopic expressions for the intrinsic viscosity for the active system of puller-like micro-swimmers are obtained. It turned out that the methods developed for passive suspensions could be applied to the active system of the pullers. For active systems, we have to actually solve the two crossover problems in a row. The first is with respect to their activity, and the second of interpolation between the two concentration regimes.

In Section 4, special attention is given to the problem of calculation of the critical indices and amplitudes of the shear modulus of random elastic composites.

In Section 5, we consider the sedimentation problem and discuss the sedimentation rate in the vicinity of critical points. Critical indices are calculated from the short truncation by means of minimal difference and minimal derivative conditions imposed directly on the critical properties.

In Section 6, crossover expressions are presented for the non-local diffusion coefficient of a simple liquid (another random system) in the vicinity of a critical point, dependent on the wave-number of density fluctuations in the liquids. Extensions and corrections to the celebrated Kawasaki formula are discussed.

In Section 7, we discuss the effective conductivity for the classical analog of graphene with the random addition of vacancies and calculate the effective critical index for superconductivity dependent on the concentration of vacancies.

Finally, in Section 8, we discuss the effective conductivity of a random 3D composite with perfectly conducting spherical inclusions and calculate the superconductivity critical index of such a random 3D composite.

Methods required in the course of the analysis are discussed when necessary and to the extent required. Relevant references are always given for more details to be learned. Since no one can cover the immense subject in its completeness, the presentation will be anchored to the topics most familiar to the author. Only a few problems of transport coefficients where short truncations are known for quite some time or were obtained only recently, were analyzed in regard to their critical properties or crossover behavior. However, they are all of significant interest in the theory of composites, suspensions and simple liquids. New methods developed for critical and crossover phenomena bring new insights to the effective transport coefficients.

The four effective transport coefficients, such as effective viscosity, sedimentation rate, elastic modulus and effective conductivity, are the most popular in the field of effective properties of composite media and suspensions. A similar problem of non-local diffusion coefficients is one of the central points for research in the theory of simple liquids.

However, there are some other important problems, including criticality in various contexts. Permeability in wavy channels also behaves critically or demonstrates a power-law behavior. Several examples of the effective permeability in wavy channels can be found in [7]. The rather long truncated series is given in terms of waviness [4]. Their analysis is more traditional since one can observe a good numerical convergence for the critical indices and amplitudes for the effective permeability. Striking power-laws were also discovered in the gap sensitivity in optimized photonic crystal and disordered networks as the function of the dielectric contrast [8]. The diffusion process in complex media is found to exhibit different regimes, including power-laws, dependent on the length scale [9].

2. Viscosity of Passive Suspensions

Keep in mind that in isotropic Newtonian fluid viscosity, ${\mu }_0$ is a scalar quantifying the rheology of the fluid, establishing the relation between the macroscopic strain-rate $e_{ij}$ and stress ${\sigma }_{ij}$ in the fluid. The effective interactions among the colloidal particles in suspensions are described almost perfectly by the hard-spheres model. The phase behavior of hard spheres includes liquid, crystal, and metastable states, with associated phase transitions and crossover.

Thus, hard-sphere fluid is the simplest and most widely used reference model for describing the behavior of real fluids and concentrated colloidal suspensions. The particles (colloids) can be viewed as large atoms. There is an understanding that the hard spheres equation of state, as well as various effective transport properties, may be extended smoothly, avoiding the first-order freezing transition to describe the metastable fluid branch. In the most interesting three-dimensional case, the metastable branch ends at its first singularity at a random close packing (RCP) volume fraction $f_c$ [4,10].

The fundamental idea of Einstein (1906) was to quantitatively express the notion that the viscous energy dissipation rate of the suspension must be equal to the dissipation rate of the effective homogeneous fluid. Such a condition leads to the definition of the effective viscosity applicable for any volume fractions up to the threshold where lubrication forces dominate [11]. Such an approach naturally includes both passive and active suspensions of cells and bacteria considered in the subsequent Section 2 and Section 3. The effective viscosity could be measured experimentally and may become negative in open bacterial systems as well [12].

While our main interest is in the most physically and biologically relevant case of three-dimensional (3D) suspensions, models of two-dimensional (2D) suspensions can serve as important additional tests of our method. Furthermore, understanding the rheological properties of 2D suspensions may carry considerable importance for the understanding of the function of lung surfactant monolayers, even in developing possible treatments, since the divergent surface viscosity plays a crucial role in breathing [13,14].

Since both two-dimensional and three-dimensional contributions to the viscosity are present for thin films [15], the importance of studying idealized 2D systems is apparent. Thus, we first review the case of the 2D suspensions of impenetrable hard disks immersed in a fluid of viscosity ${\mu }_0$.

2.1. Effective Viscosity of 2D Passive Suspensions

The effective viscosity problem can be approached as a peculiar critical phenomenon. The expansion of effective viscosity in the volume fraction of hard disks $f=\pi r^{2}n$, where r is the disks radius and n is the number density of particles, is available theoretically only up to the first-order term as $f\to 0$ [16,17],

$$\mu \left(f\right)/{\mu }_0\simeq 1+2f,$$

(1)

where ${\mu }_0$ is the viscosity of the solvent. Clearly, such truncation is not sufficient to develop any meaningful and consistent approach to high-volume fractions. We have to recognize that, in addition, there is a critical point characterized by a power-law divergence of the viscosity in the vicinity of the maximum volume fraction value $f_c=\frac{\pi }{\sqrt{12}}\approx 0.9069$,

$$\mu \left(f\right)/{\mu }_0\simeq A{(f_c-f)}^{-\mathcal{S}},f\to f_c.$$

(2)

The value of the critical exponent for viscosity $\mathcal{S}\approx 4/3$ has not been explained theoretically but is known experimentally [18,19]. The value of critical amplitude A has to be estimated in the course of calculations. From simulations with soft disks in two dimensions, using a bi-disperse mixture with equal numbers of disks of two different radii, it is found that $\mathcal{S}=1.7\pm 0.2$ [20].

In our approach to the derivation of the analytical expressions for the effective viscosity of 2D suspensions, we are going to employ the renormalization group (RG)-inspired technique of [21]. The next component of our approach in our manipulations is to dwell on the intuition of the great scientists. Feynman thought to avoid tedious calculations if possible and replace them with rather simple but direct estimates for the coefficients in the perturbative theories.

A couple of analytical expressions for the effective viscosity have been developed over the years. Here, we include several that we find the most useful for comparison with the results derived in this paper. The most often used is the well-known Krieger–Dougherty (KD) equation, applicable in 2D as well as in 3D.

$$\begin{array}{c}\mu \left(f\right)/{\mu }_0={\left(1-\frac{f}{f_c}\right)}^{-a_1{f}_c}.\hfill \end{array}$$

(3)

This is a crossover expression that links polynomial and power–law behaviors at opposite ends of the concentration interval. It can be written explicitly as

$$\mu \left(f\right)/{\mu }_0={(1-1.10266f)}^{-1.8138},$$

so that one can note that the critical exponent from (3) is much larger than $4/3$ but is more compatible with [20].

As the benchmark, we will also use the 2D variant of the simple crossover equation suggested in [11]

$$\mu \left(f\right)/{\mu }_0={\left(1-\frac{f}{f_c}\right)}^{-\mathcal{S}}\left(1-q\frac{f}{f_c}\right),$$

(4)

where $q=\mathcal{S}-a_1{f}_c$. From (4), one can evaluate the second-order coefficient in the expansion $a_2=2.67023$, as well as the critical amplitude $A=1.2996$.

The main shortcoming of the renormalization group method employed in [22] is that the authors only use a perturbative theory for small f and exclude critical point behavior. Below, we specialize the formulae to

$$\Phi =\mu ,x=f,x_c=f_c,$$

with the critical index denoted by $\beta =-\mathcal{S}$. By definition, in the vicinity of the critical point, $\Phi$ has the following asymptotic behavior as $x\to x_c$

$$\Phi \left(x\right)\simeq A{(x_c-x)}^{\beta },$$

(5)

where A is the critical amplitude.

The critical exponent $\beta$ is the limit $x\to x_c$ of the effective critical exponent $\beta \left(x\right)$ defined in a neighborhood of $x_c$:

$$\beta \left(x\right)=(x_c-x)\frac{d}{d(x_c-x)}ln\Phi \left(x\right).$$

(6)

Let us note that with help of the Euler transformation,

$$z\left(x\right)=\frac{x}{x_c-x},x\left(z\right)=\frac{x_{c}z}{1+z},$$

(7)

the critical behavior can always be moved to infinity.

If we define a new quantity ${\mathcal{B}}_a\left(x\right)={\partial }_{x}log\left(\Phi \left(x\right)\right),$ then

$$\beta \left(x\right)=(x-x_c){\mathcal{B}}_a\left(x\right).$$

On the other end of the interval at small x the function $\Phi$ has a polynomial truncation, with the k-th order approximant given by the series

$${\Phi }_k\left(x\right)=1+\sum _{n=1}^k{a}_n{x}^n,$$

(8)

with integer $k\ge 1$. In the k-th order effective exponent can be defined here analogously by

$${\beta }_k\left(x\right)=-(x_c-x)\frac{d}{dx}ln{\Phi }_k\left(x\right),$$

(9)

which, for $k=1$, has an explicit expression valid for small x,

$$\frac{{\beta }_1\left(x\right)}{x-x_c}\simeq a_1.$$

(10)

One would strive to extend $\beta \left(x\right)$ to the whole interval $[0,x_c)$. To such an end, one can derive a self-similar approximation [21,23], for the effective critical exponent ${\beta }_k^{*}\left(x\right)$. It should also be taken into account that

$${\beta }_k^{*}\left(x\right)\to \beta ,x\to x_c.$$

(11)

The value of $\beta$ may be assumed to be known in advance, and its value is finite. However, it can be calculated as well from the expansion at small x.

Finally, from the knowledge of ${\beta }_k^{*}\left(x\right)$ one can approximate the sought function, ${\Phi }_k^{*}\left(x\right)$, by integrating (6), obtaining

$${\Phi }_k^{*}\left(x\right)=exp{\int }_0^{x}dX\frac{{\beta }_k^{*}\left(X\right)}{X-x_c}.$$

(12)

The critical amplitude A can be evaluated from (12) as

$$A=\underset{x\to x_c}{lim}{\Phi }_k^{*}\left(x\right){(x_c-x)}^{-\beta }.$$

(13)

The formal method of self-similar approximants may, in general, lead to multiple solutions to ${\beta }^{*}$ and, hence, to ${\Phi }^{*}$. Indeed, the formally constructed fixed point depends on unknown exponents and amplitudes, which parametrizes the most general expression for ${\beta }_k^{*}$. With sufficiently many terms in the perturbative expansion, it may be possible to determine the parameters uniquely by selecting them so that all the $a_k$, as well as the critical threshold and exponent, are reproduced. When only low-order expansions are at our disposal, as is frequently the case, additional choices have to be made to select appropriate parameter values for self-similar approximants.

The most frequently used here are what we refer to below as the root, iterated root, factor, and super-exponential approximants [3,4]. Some precise definitions of the relevant constructions will be given a bit later in Section 2.1.1. For convenience, most important notions will be collected in the Nomenclature section. However, we think that for many simpler cases have been discussed in the first order of perturbation theory, and even in the more advanced second-order, emerging formulas are transparent enough and are self-explanatory. In our opinion, a clever researcher could have guessed them if needed.

It is important to be able to make the extrapolation by several methods and then compare the results. If these results yield close values, this suggests that the extrapolation (or interpolation) is reliable.

In line with this idea, we employ the self-similar approximants within the single framework and strive to show that they are close to each other. The latter observation would indicate that the basic approach is robust and not overly sensitive to the uncertainty arising from working with small-concentration expansions or incomplete knowledge of the critical point.

We refer to our own results as obtained by the Algebraic renormalization method (ARM). ARM dwells on studying and drawing conclusions from various expressions for the effective critical index ($\beta$-function) in the lowest orders of perturbation theory.

Using a change in variables $z=f_c-f$ and the simplest root approximant [24], one can obtain the simplest expression for the effective critical index from the first-order perturbative expansion:

$${\beta }_{1,1}^{*}\left(f\right)=-z_0(f_c-f)\left(1+z_1\frac{1}{f_c-f}\right),$$

(14)

where

$$z_0=a_1-\frac{\mathcal{S}}{f_c},z_1=\frac{\mathcal{S}}{z0}.$$

This is a linear expression, which after integration of (12) with ${\beta }_1^{*}={\beta }_{1,1}^{*}$ yields the following formula

$${\mu }_{1,1}^{*}\left(f\right)/{\mu }_0=exp\left(z_{0}f\right){\left(\frac{f_c}{f_c-f}\right)}^{\mathcal{S}}=\frac{{\pi }^{4/3}}{2\sqrt[3]{2}{3}^{2/3}}{e}^{\left(2-\frac{8}{\sqrt{3}\pi }\right)f}{\left(\frac{1}{\frac{\pi }{2\sqrt{3}}-f}\right)}^{4/3}.$$

(15)

Asymptotically close to the critical point, we recover a power-law, while at small concentrations, the exponential term dominates. The second-order coefficient and the critical amplitude are easily evaluated from (15) and found to be

$$a_{2,1}=2.81,A_{1,1}=1.419.$$

One can also look for the solution in the form of an iterated root. Iterated roots of the type

$${\beta }_{1,2}^{*}\left(f\right)=a_1{\left({(1+v_{1}f)}^2+v_2{f}^2\right)}^{-1}(f-f_c),$$

(16)

were designed to find all the parameters iteratively (see [25] and Section 2.1.1 for a complete definition of iterated roots) in such a way that

$${(1+v_1{f}_c)}^2+v_2{f}_c^2=0,$$

(17)

so that

$$v_1=\frac{a_1}{2\mathcal{S}}-\frac{1}{f_c},v_2=-\frac{a_1^2}{4{\mathcal{S}}^2}.$$

(18)

The result of integration can be expressed in closed form by means of a simple formula

$${\mu }_{1,2}^{*}\left(f\right)/{\mu }_0={\left(1+\frac{bf}{f_c-f}\right)}^{\mathcal{S}},$$

(19)

where $b=\frac{\sqrt{3}\pi }{4}$. The intermediate expressions occurring in the course of integration are not shown here. They are rather cumbersome [21], but the simple formula can be established through the direct comparison of the expansions in arbitrary order. From the expression (19), we easily analytically evaluate

$$a_{2,2}=2.705,A_{1,2}=1.323,$$

again, larger than what follows from (4).

However, one can also try to approximate the effective index by the following simple exponential approximant [26],

$${\beta }_{1,3}^{*}\left(f\right)=-a_1{f}_{c}exp\left(\tau f\right),$$

(20)

where $\tau$ is obtained from the condition ${\beta }_{1,3}^{*}\left(f_c\right)=-\mathcal{S}$.

In this case, the effective viscosity can be written using an exponential integral

$$Ei\left(x\right)={\int }_{-\infty }^x\left(e^t/t\right)dt,$$

so that

$${\mu }_{1,3}^{*}\left(f\right)/{\mu }_0=2.{71828}^{1.33333(-Ei(0.339334(0.9069-f))-0.268144)}.$$

(21)

Directly, one can obtain the second-order coefficient in the expansion and the amplitude. It equals

$$a_{2,3}=2.763,A_{1,3}=1.369.$$

Let us also try to estimate the third-order coefficient, assuming that $a_2=2.75$, the value chosen between the three estimates for the second-order coefficient just derived above. Note that one is not obliged to include the trivial zero-order term into the resummation directly, but simply can leave it outside and just re-use the results of ARM. Then, there are two simple expressions for EV, which respect one exact and one calculated coefficient,

$${\mu }_{2,1}^{*}\left(f\right)/{\mu }_0=1+2fexp\left(b_{1}f\right){\left(1-\frac{f}{f_c}\right)}^{-\mathcal{S}}=1+2e^{-0.0952104f}f{\left(1-\frac{2\sqrt{3}f}{\pi }\right)}^{-4/3}.$$

(22)

It generates the third-order coefficient in the expansion and the critical amplitude

$$a_{3,1}\approx 3.51,A_{2,1}\approx 1.46.$$

Based on the same idea, the second expression of such type can be readily written down,

$$\begin{array}{c}{\mu }_{2,2}^{*}\left(f\right)/{\mu }_0=1+2f{\left(\frac{d_{1}f}{f_c-f}+1\right)}^{\mathcal{S}}=1+2f{\left(\frac{0.93524f}{\frac{\pi }{2\sqrt{3}}-f}+1\right)}^{4/3}.\hfill \end{array}$$

(23)

It generates the corresponding third-order coefficient in the expansion and the critical amplitude

$$a_{3,2}\approx 3.5,A_{2,2}\approx 1.46.$$

In the above expressions, unknown coefficients $b_1$ and $d_1$ are found from asymptotic equivalence with the truncated expansions.

Finally, using the third suggested form of the $\beta$-function (20), we find after the same shift yet another expression in terms of the incomplete gamma function,

$${\mu }_{2,3}^{*}\left(f\right)/{\mu }_0=1+0.107521f{e}^{1.33333\Gamma (0,0.0669518-0.0738249f)}.$$

(24)

Formula (24) produces the third-order coefficient in the expansion and the critical amplitude

$$a_{3,3}\approx 3.51,A_{2.3}\approx 1.$$

All of our above solutions work well compared with the numerical data up to $f=0.25$. For example, at $f=0.25$, (15) predicts the EV of $1.75$ against numerical $1.8$, with similar accuracy for smaller f. However, one should take into account that it is a much smaller value for maximum volume fraction $f_c=\pi /4$ used in computations of [27], which are available up to $f=0.25$. The numerical agreement is expected to become better when the higher threshold is used.

Another convenient expression for the effective viscosity can be suggested by directly applying the factor approximant,

$${\mathcal{F}}^{*}\left(f\right)=\frac{{(0.228655f+1)}^{2.31698}}{{(1-1.10266f)}^{4/3}},$$

(25)

with the third-order expansion coefficient and the critical amplitude

$$a_3\approx 3.44,A\approx 1.36,$$

appearing close to the other estimates. In general form, the factor approximants will be presented a bit later in the (34). However, it is quite straightforward to observe that formula (25) is a multiple of two factors, each of them being the simplest root-form [3,23].

Thus, in the 2D passive case, we suggest, for instance, formula (23) for the whole range of concentrations. The predictions of the formula are in good agreement with the results of direct numerical simulation [27]. We also estimated the second-order expansion coefficient, the third-order expansion coefficient and the critical amplitude,

$$a_2\approx 2.75,a_3\approx 3.5,A\approx 1.46.$$

In comparison, in [28], a slightly larger result for the second-order coefficient is obtained. It was found that $a_2=3.6$ by extrapolation from the confined geometry for athermal disks. Such a difference in $a_2$ with our estimates is to be expected since, in our study, we do not specify the strength of thermal fluctuations but draw additional information from high-concentration limits, where thermal and athermal cases look very similar.

To check the consistency of such estimated values for $a_2,a_3$, let us try to estimate the critical index along the different route, as suggested in [29,30]. Assume also that in the vicinity of a threshold,

$$\mu \sim {(f_c-f)}^{-1-\mathsf{v}},$$

with unity to be expected as the usual contribution from the radial distribution function at the particle’s contact $G(2,f)$, [29,30], and the value of $\mathsf{v}$ coming from the hydrodynamic interactions in the suspension. With the total value of the critical index $\mathcal{S}\approx 4/3$, we expect that $\mathsf{v}\approx 1/3$. One can suggest, for the effective viscosity, a simple root approximant

$$\mathcal{R}\left(z\right)=1+b_{1}z{(1+b_{2}z)}^{-\mathsf{v}},z=\frac{f}{f_c-f},$$

which, by design, takes the unity contribution into account. After imposing the asymptotic equivalence with the polynomial

$${\mu }_3\left(f\right)=1+2f+2.75f^2+3.5f^3,$$

in the third-order, we find

$$\begin{array}{c}{\mu }^{*}\left(f\right)/{\mu }_0=1+\frac{1.8138f{\left(\frac{0.689576f}{0.9069-f}+1\right)}^{0.358173}}{0.9069-f},\hfill \end{array}$$

(26)

with an estimated value for $\mathcal{S}\approx 1.358$, reasonably close to expected. The amplitude could be estimated from (26) as well, bringing $A\approx 1.39$.

2.1.1. Comment on Critical Index and Diff-Log Borel Summation

Consider the case when $\Phi$ has the following asymptotic behavior as $x\to \infty$

$$\Phi \left(x\right)\simeq A{x}^{\beta },$$

(27)

Furthermore, let the critical index $\beta$ be known in advance. One can also calculate the critical amplitudes using the method of self-similar Borel summation introduced in [31]. Mind that the problem of critical behavior with a finite threshold (2) can always be mapped to the problem (27) by applying Euler transformation (7), or similar, as in (57).

The Borel summation method could be devised through the Borel transform,

$$B_k\left(x\right)=\sum _{n=0}^k\frac{a_n}{\Gamma (1+n)}{x}^n.$$

(28)

The intention here is to simplify the problem by removing the superficial factorial growth of the coefficients $a_n$. Working on the transformed series with some other approximants may bring about more accurate results. However, returning to the original problem is not so easy. To perform an inverse transformation, one should compute the integral numerically. However, recently, some progress was achieved in [31,32], allowing us to find analytical expressions for the critical amplitude A and index $\beta$. In the case of the Borel summation, one can easily avoid numerics and use only formulas. Calculations are particularly easy to perform with the iterated root approximants.

Summing up the truncated series $B_k(x,u)$ by means of iterated self-similar root approximants [3,4,32] yields the approximant

$$B_k^{*}\left(x\right)={\left({\left({\left({(1+{\mathcal{P}}_{1}x)}^2+{\mathcal{P}}_2{x}^2\right)}^{3/2}+{\mathcal{P}}_3{x}^3\right)}^{4/3}+\dots +{\mathcal{P}}_k{x}^k\right)}^{\beta /k}.$$

(29)

The inverse Borel transformation

$${\Phi }_k^{*}\left(x\right)={\int }_0^{\infty }{e}^{-t}{B}_k^{*}\left(tx\right)dt,$$

(30)

gives the approximation ${\Phi }_k^{*}\left(x\right)$ for the sought function $\Phi \left(x\right)$. At large values of x, function (30) is expressed as follows,

$${\Phi }_k^{*}\left(x\right)\simeq A_k{x}^{\beta }(x\to \infty ),$$

(31)

where the critical amplitude

$$A_k=C_k\Gamma (1+\beta ).$$

(32)

The marginal amplitude $C_k$

$$C_k={\left({\left({\left({\mathcal{P}}_1^2+{\mathcal{P}}_2\right)}^{3/2}+{\mathcal{P}}_3\right)}^{4/3}+\dots +{\mathcal{P}}_k\right)}^{\beta /k},$$

(33)

follows from the definition of iterated root.

One can also apply for the purpose of the Borel summation the self-similar factor approximants [3,4] in the form

$$B_k^{*}\left(x\right)=\prod _{j=1}^{N_k}{(1+{\mathcal{P}}_{j}x)}^{n_j},$$

(34)

where

$$\begin{array}{c}\hfill N_k=\left\{\begin{array}{cc}k/2,\hfill & k=2,4,6,\dots \hfill \\ (k+1)/2,\hfill & k=1,3,5,\dots \hfill \end{array}\right..\end{array}$$

(35)

The condition on the exponent

$$\beta =\sum _{j=1}^{N_k}{n}_j.$$

should be respected as well. For the factor approximants, marginal amplitude for the transformed series is given by the expression

$$C_k=\prod _{j=1}^{N_k}{\mathcal{P}}_j^{n_j},$$

(36)

Calculations could be performed with relative ease using formulas, strictly following the structure of the example presented in (Section 7, [31]). To the Borel-transformed series, we apply factor or iterated root approximants and simply find the leading term at infinity, behaving as a power law. The power calculated in such a manner is nothing else but the sought critical index [31].

As for the amplitude, it consists of the critical amplitude of the factor approximants applied to the transformed series, scaled by a Gamma-function appearing in the course of inverse transformation and dependent only on the found critical index [31]. Thus we find $\mathcal{S}\approx 1.245$, $A\approx 1.518$.

With the value of the critical index fixed to $1.3$ and factor approximants given by (34), we can recalculate $A\approx 1.417$ by using the same method of self-similar Borel summation [31]. If, instead of a factor approximants, we employ iterated roots, the result remains quite close, i.e. $A\approx 1.433$. In both situations, we observed reasonably good numerical convergence to the quoted numbers.

Critical indices can also be found by using the techniques developed for critical amplitudes. Let us deal now with a function with power-law asymptotic behavior

$$\Phi \left(x\right)\simeq A{x}^{\beta }(x\to \infty ).$$

The critical exponent can be expressed as a diff-log function of the limit

$$\beta =\underset{x\to \infty }{lim}x\frac{d}{dx}ln\Phi \left(x\right)\equiv \underset{x\to \infty }{lim}x\psi \left(x\right),$$

(37)

where $\psi \left(x\right)\equiv \frac{d}{dx}ln\Phi \left(x\right),$ as shown, e.g., in the books [3,4].

When the small-variable expansion for the original function is given by the sum ${\Phi }_k\left(x\right)$, we can find the small-variable expression for the diff-log function $\psi \left(x\right)$,

$${\psi }_k\left(x\right)=\frac{d}{dx}ln{\Phi }_k\left(x\right),(x\to 0),$$

which can be expanded in powers of x, leading to the truncated series

$${\psi }_k\left(x\right)=\sum _{n=0}^k{d}_n{x}^n.$$

(38)

However, for $x\to \infty$, we find that simply

$$\psi \left(x\right)\simeq \beta x^{\delta },$$

where the “ critical amplitude” is the sought critical index $\beta$, and the “critical index” is fixed to $\delta \equiv -1$ to make the limit (37). We can also apply the technique of iterated roots described above for the critical amplitudes and calculate the critical index $\beta$. Then, it remains to simply apply the methods to be developed to calculate also the critical index $\beta$. However, the case of $\delta =-1$, appears to be divergent, or indeterminate, when treated by the self-similar Borel summation technique of [31]. In such cases, some other types of summation could be used, similarly to in the paper [32]. Such summations, called Borel–Leroy and Mittag–Leffler summations, include some control parameters [32]. The parameters could be found from the minimal-difference, or minimal derivative conditions imposed on critical amplitudes [4].

There is also some simple way to avoid divergence. Let us consider an inverse,

$$\psi {\left(x\right)}^{-1}\simeq \gamma x^{-\delta }.$$

With such a critical index, there is no divergence in the formulas of [31]. Now the value of critical amplitude $\gamma$ can be found following [31]. After resummation, say with iterated roots [32], we will arrive at the resummed amplitude ${\gamma }^{*},$ and

$$\frac{1}{{\gamma }^{*}}\equiv {\beta }^{*},$$

where ${\beta }^{*}$ gives the sought value of the critical index $\beta$.

The technique just discussed can be readily applied to the case of critical index $\mathcal{S}$. The first-order iterated root properly conditioned at infinity gives ${{\mathcal{S}}_1}^{*}\approx 1.374$, while in the second-order ${{\mathcal{S}}_2}^{*}\approx 1.39$. Application of the second-order factor approximant leads to the estimate very close to the latter value.

Both methods of Borel summation [31,32], also explained in Section 2.1.1, especially with iterated roots, are very easy to work with since the computations are performed on explicitly given formulas.

2.1.2. Comment on Critical Index and Elasticity-Viscosity Analogy

Consider special 2D composite material with an incompressible matrix and perfectly rigid, incompressible inclusions. It represents an elastic analog of the viscous suspension [4].

The snapshot of a random media studied below is displayed in the book [4]. Lengthy but straightforward Monte-Carlo computations result in the polynomial formula for the effective viscosity of 2D random suspension valid for small f,

$$\frac{{\mu }_e}{{\mu }_0}\simeq 1+2f+2f^2+3.81137f^3+8.27305f^4+5.33447f^5+20.5096f^6.$$

(39)

The elasticity-viscosity analogy, in principle, allows us to take into account all elastic & hydrodynamic interactions between inclusions and map them into effective interactions of the pairs, triplets, etc. of groups of the inclusions. Of course, there is no trace of a Brownian motion in such an approach. On the other hand, less sophisticated estimates by ARM presented in the main body of the section allows us to potentially take them into account indirectly.

Application of the self-similar Borel summation with factor approximants from Section 2.1.1, generates single real, non-trivial estimates in the third-order of the perturbation theory, ${{\mathcal{S}}_3}^{*}\approx 1.073$. The higher-order estimates produce only non-physical, complex solutions. Application of the diff-log Borel summation with iterated roots and inverse transformation discussed in Section 2.1.1 brings again only a single non-trivial estimate in, de-facto, the same order, ${{\mathcal{S}}_1}^{*}\approx 1.301$. Again, only complex results are generated in higher orders.

Let us employ the Euler transformation (7),

$$z\left(f\right)=\frac{f}{f_c-f},f\left(z\right)=\frac{f_{c}z}{1+z}.$$

In terms of the transformed variable, let us construct the iterated root approximants introduced in Section 2.1.1,

$${\mathcal{R}}_k^{*}\left(z\right)={\left({\left({\left({(1+{\mathcal{P}}_{1}z)}^2+{\mathcal{P}}_2{z}^2\right)}^{3/2}+{\mathcal{P}}_3{z}^3\right)}^{4/3}+\dots +{\mathcal{P}}_k{z}^k\right)}^{{\mathcal{S}}_k/k},$$

(40)

and define the parameters ${\mathcal{P}}_j$ from the asymptotic equivalence with the original truncated series for the effective viscosity. This gives the large-z asymptotic form

$${\mathcal{R}}_k^{*}\left(z\right)\simeq A_k{z}^{{\mathcal{S}}_k}(z\to \infty ),$$

(41)

where the amplitudes $A_k=A_k\left({\mathcal{S}}_k\right)$ are

$$A_k={\left({\left({({\mathcal{P}}_1^2+{\mathcal{P}}_2)}^{3/2}+{\mathcal{P}}_3\right)}^{4/3}+\dots +{\mathcal{P}}_k\right)}^{{\mathcal{S}}_k/k}.$$

(42)

In order to calculate the critical index $\mathcal{S}$, we analyze the differences for the critical amplitudes in k-th order, as prescribed in [4,33],

$${\Delta }_{kn}\left({\mathcal{S}}_k\right)=A_k\left({\mathcal{S}}_k\right)-A_n\left({\mathcal{S}}_k\right).$$

(43)

The control parameters are defined as the solution to the minimization problem

$$\mid {\Delta }_{kn}\left({\mathcal{S}}_k\right)\mid ={\min }_{\mathcal{S}}\mid {\Delta }_{kn}\left(\mathcal{S}\right)\mid .$$

The results of calculations are presented in

Table 1. Critical indices for the viscosity ${\mathcal{S}}_k$ obtained from optimization conditions on ${\Delta }_{kn}\left({\mathcal{S}}_k\right)$.

${\mathcal{S}}_{\mathit{k}}$	${\mathit{\Delta }}_{\mathit{k}+1}\left({\mathcal{S}}_{\mathit{k}}\right)=0$	${\mathit{\Delta }}_{\mathit{k}6}\left({\mathcal{S}}_{\mathit{k}}\right)=0$
${\mathcal{S}}_1$	0.907	1.37
${\mathcal{S}}_2$	n.a	n.a
${\mathcal{S}}_3$	n.a	1.19
${\mathcal{S}}_4$	1.23	1.317
${\mathcal{S}}_5$	1.16	1.16

. As shown in [4], the most reasonable estimate for the index equals $\mathcal{S}=1.26\pm 0.09$,

The latter value for the index appears to be compatible with the experimentally found number $1.24$, quoted in the book [1], justified by the elasticity-viscosity analogy [4]. It is smaller than the otherwise measured value of $1.3$ since only an elastic (hydrodynamic) contribution to the index is considered. Currently, we are not aware of the consistent estimates for the contributions from the 2D Brownian motion, unlike in the 3D case to be discussed in what follows.

2.2. 3D Viscous Suspensions

We now apply the ARM method to the case of passive suspensions of spherical particles. The effective viscosity of random suspensions of the hard spheres with the stick boundary condition is always of great interest, as the researchers strive for a long time to obtain an accurate and compact formula for arbitrary concentrations. The expansion in the volume fraction of hard spheres $f=\frac{4}{3}\pi r^{3}n$, where r is the spheres radius and n is the number density, is available up to the second-order term as $f\to 0$:

$$\mu \left(f\right)/{\mu }_0\simeq 1+a_{1}f+a_2{f}^2,$$

(44)

and for non-interacting hard spheres, the result $a_1=5/2$ is due to Einstein (1906).

Finding the second-order coefficient turned out to be a great challenge. Batchelor [34] found that $a_2=6.2$, and Wajnryb and Dahler [35] obtained $a_2=5.9147$. The result for $a_2$ takes into account the two different mechanisms. It consists of the hydrodynamic contribution $\approx 5$, with the addition of the smaller contribution from the Brownian motion of the order of unity. We tend to use Batchelor’s estimate. However, we will also recall the Wajnryb and Dahler estimate when required. The effective shear viscosity behaves critically,

$$\mu \left(f\right)/{\mu }_0\simeq A{(f_c-f)}^{-\mathcal{S}},f\to f_c.$$

(45)

The threshold here corresponds to the random closest packing of hard spheres at $f_c\simeq 0.64$. The value of the critical exponent for the viscosity in 3D, $\mathcal{S}\simeq 2$ can be satisfactorily deduced both theoretically and experimentally [19,29].

For practical purposes, for three-dimensional situations, even with varying threshold values and model assumptions, it was suggested to consider such divergence as a universal power-law with the value of $\mathcal{S}=2$ [36].

2.2.1. First-Order Interpolation

Thus, near the critical point, the EV in the 3D case has a divergent behavior (45), with $\beta =-\mathcal{S}$. At the other end of the interval for small f, it is expressed as the second-order polynomial. Let us proceed similarly to the 2D case studied above and in accordance with the paper [21], where more technical details could be looked up. Our primary target is a subsequent application to the case of active suspensions in Section 3 since in that case, only first-order information is available.

In addition to ${\beta }_1\left(x\right)$, we are now able to construct in the next Section 2.2.2${\beta }_2\left(x\right)$. For $k=1$ and small x, we have an explicit expression (10),

$$\frac{{\beta }_1\left(x\right)}{x-x_c}\simeq a_1.$$

In terms of the variables $z=f_c-f$, one can easily obtain the simplest expression for the effective viscosity [24],

$${\mu }_{1,1}^{*}\left(f\right)/{\mu }_0=exp\left(z_{0}f\right){\left(\frac{f_c}{f_c-f}\right)}^{\mathcal{S}}=\frac{0.4096e^{-0.625f}}{{(0.64-f)}^2}.$$

(46)

One can also look for the $\beta$-function in the form of an iterated root

$${\beta }_{1,2}^{*}\left(f\right)=a_1{\left({(1+v_{1}f)}^2+v_2{f}^2\right)}^{-1}(f-f_c),$$

(47)

where

$$v_1=\frac{a_1}{2\mathcal{S}}-\frac{1}{f_c},v_2=-\frac{a_1^2}{4{\mathcal{S}}^2}.$$

After integration in (12) with ${\beta }_1^{*}={\beta }_{1,2}^{*}$, we arrive, after some transformations, at the surprisingly simple explicit solution:

$${\mu }_{1,2}^{*}\left(f\right)/{\mu }_0={\left(1+\frac{0.8f}{0.64-f}\right)}^2.$$

(48)

Another way to approach the problem is simply to approximate the effective index by the exponential approximant [26],

$${\beta }_{1,3}^{*}\left(f\right)=-a_1{f}_{c}exp\left(\tau f\right),$$

(49)

where the control parameter $\tau$ has to be obtained from the condition

$${\beta }_{1,3}^{*}\left(f_c\right)=-\mathcal{S},$$

In this case, integral (12) can be written in the general form

$${\mu }_{1,3}^{*}\left(f\right)/{\mu }_0=exp\left(a_1{f}_{c}exp\left(\tau f_c\right)\right)\left(Ei(\tau (f_c-f)-Ei\left(\tau f_c\right))\right)$$

(50)

Directly from the expressions (46), (48), (50), one can plausibly estimate the second-order coefficient in the expansion. It equals, respectively,

$$a_2=5.566,a_2=5.469,a_2=5.514.$$

2.2.2. Second-Order Interpolation

Let us now consider the second-order approximation to the effective critical index,

$$\frac{{\beta }_2\left(x\right)}{x-x_c}\simeq a_1+(2a_2-a_1^2)x(k=2).$$

(51)

Transformation of the independent variable. In the paper [21], the following expression for the effective critical index in the form of a second-order iterated root [25],

$${\beta }_2^{*}\left(f\right)=-g_0(f_c-f){\left({(1+g_1\frac{1}{f_c-f})}^2+g_2{\left(\frac{1}{f_c-f}\right)}^2\right)}^{1/2},$$

(52)

was suggested. The parameters should follow from the asymptotic equivalence in the region of small f as well as from the critical regime. In general form, iterated roots are given by formula (40). Numerical analysis of the integral (12) for effective viscosity allows us to find the critical amplitude, $A=0.342.$

Iterated Root. However, another expression was suggested in [21]

$${\beta }_2^{*}\left(f\right)=-a_1(f_c-f){\left({({(1+k_{1}f)}^2+k_2{f}^2)}^{3/2}+k_3{f}^3\right)}^{-1},$$

(53)

utilizing the higher-order iterated root [25]. Imposing the condition on critical points together with usual asymptotic conditions, one can obtain all parameters in (53) explicitly. Again, the EV is given numerically by the integral (12). One can also find the corresponding critical amplitude, $A=0.398.$

Super-exponential and root approximants. One can also construct the so-called super-exponential approximant [26],

$${\beta }_2^{*}\left(f\right)=-a_1{f}_{c}exp({\tau }_{1}fexp\left({\tau }_{2}f\right)),$$

(54)

as shown in [21]. The parameters $\tau$ can be found in [21]. After numerical integration (12), one can find $A=0.34.$

Similarly, but using a root approximant [3,24]

$${\beta }_2^{*}\left(f\right)=-a_1{f}_c{(1+{\mathcal{P}}_{1}f)}^m,$$

(55)

with

$$m=0.134,{\mathcal{P}}_1=\frac{(2a_2-a_1^2)f_c-a_1}{f_c{a}_{1}m},$$

after integration (12), one can then find the critical amplitude $A=0.309.$

Averaged approximants. A different approach to the second-order problem could be based on pairwise weighted averages of the expressions obtained in the first-order, such as (46), (48), (50). In total, three such pairs are available. For example,

$${\mu }_{1,2}^{**}\left(f\right)=p_1{\mu }_{1,1}^{*}\left(f\right)+(1-p_1){\mu }_{1,2}^{*}\left(f\right).$$

(56)

The weight $p_1$ is determined through the asymptotic condition on the second-order coefficient $a_2$. By analogy, one obtains ${\mu }_{1,3}^{**}$ and ${\mu }_{2,3}^{**}$. Averaged approximants will be employed in Section 3.

2.2.3. Factor Approximants and Estimates for the Critical Index

As explained in the paper [21], one can take into account some additional proposition on the expansion at small f for finite-size systems [37] by means of the following variables transformation

$$z=\frac{{(f/f_c)}^{1/2}}{1-{(f/f_c)}^{1/2}}$$

(57)

which admits an inverse

$$f=f_c\frac{z^2}{{(1+z)}^2}.$$

The factor approximant is then obtained from the trial expansion in the parameter $f^{1/2}$,

$${\mathcal{F}}^{*}\left(f\right)={\left(1+{\mathcal{P}}_1\frac{{(f/f_c)}^{1/2}}{1-{(f/f_c)}^{1/2}}\right)}^{\mathcal{S}-m}{\left(1+{\mathcal{P}}_2\frac{{(f/f_c)}^{1/2}}{1-{(f/f_c)}^{1/2}}\right)}^m,$$

(58)

as suggested and explained in more detail in [4,21]. Here,

$${\mathcal{P}}_2=\frac{3\mathcal{S}+{\left(\mathcal{S}(9\mathcal{S}-8a_1{f}_c)\right)}^{1/2}}{2\mathcal{S}},m=-\frac{2a_1\mathcal{S}}{{\mathcal{P}}_2^2\mathcal{S}-2a_1},{\mathcal{P}}_1=\frac{{\mathcal{P}}_{2}m{f}_c}{m-\mathcal{S}}.$$

(59)

From (58), one can evaluate the second-order coefficient in the expansion at small concentrations,

$$\begin{array}{ccc}{a}_2^f\left(f_c\right)=\frac{6f_c{a}_1^2{\mathcal{S}}^2-4f_c{a}_1^2\mathcal{S}+9a_1{\mathcal{S}}^2}{12{\mathcal{S}}^2{f}_c}& +& \frac{\frac{27a_1\mathcal{S}}{4}-9a_1^2{f}_c\frac{a_1{(9{\mathcal{S}}^2-8a_1\mathcal{S}{f}_c)}^{3/2}}{4{\mathcal{S}}^2}+\frac{8a_1^3{f}_c^2}{3\mathcal{S}}}{3f_c\sqrt{9{\mathcal{S}}^2-8a_1\mathcal{S}{f}_c}-4a_1{f}_c^2+9\mathcal{S}{f}_c}.\end{array}$$

(60)

This coefficient $a_2^f$ is a function of the critical concentration $f_c$. Evaluating at the random close packing fraction $f_c=0.64$ yields $a_2^f=5.859$, remarkably close to the result of [22,35], obtained after exceptionally tedious calculations. With a slight decrease in the value of the threshold to $f_c=0.636$, we estimate $a_2\approx 5.9.$ In reverse, with known $a_2=5.9147$ [22,35], we estimate $\mathcal{S}\approx 2.036$ in close agreement with the accepted value.

With the threshold $f_c=0.637$ and two known coefficients in the expansion, one can develop the factor approximant

$$\mu \left(f\right)/{\mu }_0={(1-1.56986f)}^{-2.00267}{(1+f)}^{-0.643908},$$

(61)

and almost perfectly deduce the critical index, $\mathcal{S}=2.00267$.

Assimilating the value of the index, threshold and the two coefficients for the viscosity, we arrive at the following factor approximant

$$\mu \left(f\right)/{\mu }_0={(1-1.56986f)}^{-2}{(1.01684f+1)}^{-0.629126}.$$

(62)

The critical amplitude found from such a factor approximant is equal to $0.3$. Such value can be considered most reasonable based on our various approximations.

The method of index function developed in the book [4] can be applied to calculate $\mathcal{S}$ as well. Recall first the well-known Krieger–Dougherty (KD) Equation (3),

$$\mu \left(f\right)/{\mu }_0={\left(1-\frac{f}{f_c}\right)}^{-a_1{f}_c}.$$

which has the form of the simple root-approximant [3]. One can see that the critical index can be expressed as

$$\mathcal{S}=-\underset{f\to f_c}{lim}\frac{log(\mu \left(f\right)/{\mu }_0)}{log(1-f/f_c)},$$

(63)

see also [38], where such a definition is introduced in the general case of phase transition. Let us look for the solution in a more general form,

$$\mu \left(f\right)/{\mu }_0={\left(1-\frac{f}{f_c}\right)}^{-\mathcal{S}\left(f\right)},$$

with

$$\mathcal{S}\left(f\right)=\frac{a_1{f}_c}{1+\frac{f\left({a_1}^2{f}_c+a_1-2a_2{f}_c\right)}{2a_1{f}_c}},$$

and we were able to calculate $\mathcal{S}\approx 2.02$. Overall, the index function monotonically increases from the KD value away from $f_c$ to the calculated critical index in the vicinity of $f_c$. More details on the index function can be found in the book [4].

2.2.4. Reference Formulae for 3D Passive Suspensions

We would also bring up a few useful formulas for the effective viscosity. The empirical Thomas equation [39] gives an accurate fit

$$\begin{array}{c}\mu \left(f\right)/{\mu }_0=1+2.5f+10.05f^2+0.00273exp\left(16.6f\right),\hfill \end{array}$$

(64)

valid outside of the neighborhood of $f_c$.

In [29], Brady suggested an explicit semi-empirical expression for the EV (see also [21] for more details):

$$\mu \left(f\right)/{\mu }_0={\mu }_H/{\mu }_0+\frac{12}{5}{f}^2\frac{G(2,f){\mathcal{D}}_0}{{\mathcal{D}}_s^0\left(f\right)}.$$

(65)

Several physical parameters, which are functions of f, are required: the short-time self-diffusion coefficient of suspended particles

$${\mathcal{D}}_0^s\left(f\right)={\mathcal{D}}_0(1-1.56f)(1-0.27f)$$

(66)

and the high-frequency limiting viscosity ${\mu }_H$,

$${\mu }_H=\frac{1+\frac{3}{2}f(1+f-0.189f^2)}{1-f(1+f-0.189f^2)}.$$

(67)

Both expressions above represent the experimental data of Pearson and Shikata in the form convenient for calculations. They conform to the dilute limits and vanish or diverge at $f_c=0.64$. The theory behind (67) allows us to track experimental data nicely for $0.20<f<0.5$.

For $f<0.5$, the radial distribution function $G(2,f)$ at particle contact is given approximately as

$$G(2,f)=\frac{1-f/2}{{(1-f)}^3}.$$

For higher concentrations, $f\ge 0.5$$G(2,f)$ is found approximately through simulations,

$$G(2,f)=\frac{0.78}{0.64-f}.$$

(68)

Remarkably, such considerations lead to the critical index of $\mathcal{S}=2$! More details and discussions can be found in [21,29]. Perhaps, the value $\mathcal{S}=2$ is even universal [36]. Nevertheless, such value could be employed for practical purposes, notwithstanding the threshold value, as conjectured in [36].

Brady’s approach can be conveniently extended to the whole region of concentrations [40],

$$\mu /{\mu }_0=1+\frac{2.5f}{1-\frac{f}{f_c}}+\frac{3.4f^2}{{\left(1-\frac{f}{f_c}\right)}^2}.$$

(69)

The earlier numerical RG results from [22] can be summarized as follows:

$$\mu \left(f\right)=\mu \left({\psi }_0\right)exp{\int }_1^{f/{\psi }_0}dx{\rm Y}\left(\overline{g}(x,{\psi }_0)\right),$$

(70)

with ${\rm Y}\left(g\right)$ given explicitly [22]:

$$\overline{g}(t,g)={\left(ln(exp(1/g)t^{a_1-2\frac{a_2}{a_1}})\right)}^{-1},{\rm Y}\left(g\right)=a_{1}g+(2a_2-a_1^2)g^2,$$

with ${\psi }_0=0.085$.

From the analysis of the most typical curves presented in Figure 1, one can conclude that for a typical suspension of spherical particles [39], our ARM formula matches it with good accuracy, substantially exceeding the accuracy of prior RG results [22]. More details and discussion can be found in [21]. At present, it seems that the following formula derived from the factor approximant given by the expression (58),

$$\mu \left(f\right)/{\mu }_0=\frac{{\left(\frac{0.863446\sqrt{f}}{1-1.25294\sqrt{f}}+1\right)}^{2.84988}}{{\left(\frac{2.89537\sqrt{f}}{1-1.25294\sqrt{f}}+1\right)}^{0.84988}},$$

(71)

is in better agreement with the Thomas interpolation (64) than all other formulas derived from all kinds of self-similar approximants.

3. Effective Viscosity of Suspensions of Pullers

Pullers are an important biophysical system that exhibit a range of interesting macroscopic behaviors, while being in many ways simpler than pushers [21]. Their rheology could be similar to that of passive suspensions. It is believed that suspensions of puller-like microswimmers, such as Chlamydomonas algae cells, possess interesting rheological properties [41,42].

However, the pullers are difficult to investigate experimentally. At meaningful concentrations, the algae tend to form 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 (see, e.g., [43]). This should be contrasted to the propulsion of pushers, see Figure 2, which illustrates the two types of swimmers. The pushers, such as Escherichia Coli bacteria, use a rotating flagella bundle to push the fluid behind the bacterial body (see, e.g., [44]).

The rheological behavior of pushers at present appears to be beyond the reach of analytical methods. Early experiments on pushers revealed a number of novel macroscopic properties of their suspensions. In particular, in [45,46], pushers were found to exhibit a super-diffusion for small times, with a crossover to classical diffusion on longer time intervals. This initial dramatic modification in diffusion was attributed to the advection generated by coherent bacterial motion. The second effect corresponds to the long-time diffusion enhancement compared to the passive case. The same effect was found for pullers, but in that case, it was shown to be much smaller [47,48].

The effect of different types of propulsion is even more pronounced in effective viscosity. For pushers, EV can decrease dramatically with the increasing particle concentrations but still remain in the dilute regime $f<0.3\%$ [12,49]. Even leading to a viscous-less, “super-fluid” regime with $\mu \propto 0$ or even negative [49].

Experimental and theoretical models have been developed to examine fundamental aspects of a collective motion exhibited by various biological systems. Following the papers [12,50,51,52] and works cited therein, it was suggested that hydrodynamic interactions between the swimmers lead to the collective motion when every bacterium (pushers) interacts with others through the viscous environment. Implicit evidence of the collective motion can be found in a reduction in the effective viscosity [12]. The theoretical investigations of the collective motion were based on the framework of mechanical dynamical systems [53,54,55,56,57].

Experimental EV measurements for pullers are scarce, yet there are good reasons to expect [41] that pullers behave similarly to passive suspensions with a strongly enhanced effective viscosity [41,42]. In experiments on Chlamydomonas [41], this enhancement was found to be of a smaller magnitude than the corresponding EV decrease for pushers.

In a suspension of microscopic particles, both the strain-rate and the stress are modified by the disturbance flow around the particles, and the effective viscosity $\mu$ quantifies the relationship between these volume-averaged quantities. Operationally, we can define a scalar $\mu$ as relating fixed components of e and $\sigma$ (e.g., the $xy$-components for particles suspended in a planar shearing flow) [44]. Given the low Reynolds numbers typical of suspensions, one has to solve the Stokes equation in the exterior of the microscopic particles. The solution amounts to balancing the forces and torques on the fluid. In the simplest case of dilute suspensions, the inclusions are statistically independent and do not interact. In this situation, the disturbance flow of a single typical particle is sufficient to determine the EV, and the normalized correction appears to be linear in the volume fraction of inclusions f.

The aspect ratio for prolate spheroids is quantified by the Bretherton constant $B_r=\frac{b^2-a^2}{b^2+a^2}$, in terms of its minor and major semi-axes a and b. The dynamics of particle orientations $\tau$ in a linear ambient flow were elucidated by Jeffery (1922). When the particles are subject to the Brownian motion, a unique long-time distribution of orientations can be obtained, predicting only an increase in EV. For self-propelled particles modeled as prolate spheroids with a rigidly attached point force (see Figure 3), the situation is different. Elongated particles, according to Jeffery (1922), will spend more time aligned to the stable principal axis of the ambient planar shear flow.

Having established the ARM for passive suspensions, we would like to apply the same technique to suspensions of swimming cells. Specifically, it is important to suggest the concentration dependence of EV for suspensions of a puller cell organism. The interesting characteristic of pullers is that their presence tends to monotonically increase the effective viscosity of the suspension. The latter property is in contrast with pushers. The effective viscosity of pushers may depend on concentration non-monotonically and even exhibit a significant lowering of the effective viscosity for low to moderate concentrations.

However, one would expect all swimming cells or bacteria to affect the effective viscosity of the suspension similarly. All swimmers inject energy into the fluid, thus decreasing the effective rate of viscous dissipation and leading to an effective decrease in viscosity. The usual qualitative explanation of the different behavior of pullers 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 effective viscosity is increased. When a puller, approximated by a contractile dipole, is aligned with the extensional axis, the rate of strain is decreased. Since the viscosity is inversely proportional to the rate of strain, one can think that under a constant external force, the effective viscosity should increase.

Let us keep in mind that reconstruction of EV for pullers should be accomplished only on the basis of the (effective) intrinsic viscosity,

$$a_1=\underset{f\to 0}{lim}\frac{\mu -{\mu }_0}{{\mu }_{0}f}.$$

(72)

Even then, the problem is far more complicated than for passive systems since: (1) For comparison with experiments, $a_1$ has to be calculated first as a function of the puller’s internal parameters and of the flow. Second, (2) it should be supplemented with the information or insight on critical behavior at the dense packing limit, which we argue to be similar to that of passive suspensions.

Only then, would we be in a position to apply the techniques developed for passive suspensions and justify all the extensive groundwork on passive systems with only $a_1$ available from the “microscopic” approach of Stokes equations. To reconstruct the dependence of the intrinsic viscosity on the parameters of the background flow and the shape of the swimmers, we are going to rely on the results of [44]. In the paper [44], it was shown how the case of pullers can be obtained from the same formulae as for pushers just by changing the rigidly attached propulsion force $F_p$. The analytical expressions for the intrinsic viscosity as a function of the arbitrary strain rate $\gamma$ of the background flow, as well as the Bretherton constant, were obtained in [44], in the form of another two asymptotic pieces.

For the small parameter $\xi =\frac{\gamma B_r}{{\mathcal{D}}_r}$, applicable when the background flow is much smaller than the rotational diffusion of pullers (weak flow regime), the following asymptotic is valid,

$$\frac{\mu -{\mu }_0}{{\mu }_{0}f}\simeq {\alpha }_0+{\alpha }_1\xi +{\alpha }_2{\xi }^2+O\left({\xi }^3\right),\xi \to 0,$$

(73)

with

$$\begin{array}{c}{\alpha }_0=M_0;{\alpha }_1=\frac{F_{p}K}{160\pi a^2{\mu }_0\gamma };{\alpha }_2=M_2.\hfill \end{array}$$

(74)

Here, ${\mathcal{D}}_r$ is a tumbling constant specifying the strength of rotational diffusion.

For the large $\xi$, the background flow dominates the rotational diffusion (strong flow regime), and a different asymptotic is valid,

$$\frac{\mu -{\mu }_0}{{\mu }_{0}f}\simeq L=S+\frac{F_{p}K}{16\pi a^2{\mu }_0\gamma }+O(1/\xi ),\xi \to \infty ,$$

(75)

Coefficients $M_0,M_2,K$, S and N are the shape factors, dependent only on the geometry of the puller [44]. In the aligned limit of a strong flow, the particles tend to align along the extensional axis of the flow.

It just so happens that the experimental data from [41] are outside the asymptotic regimes of (73) and (75). The only way to obtain estimates for the viscosity at real parameters is to find a good or just reasonable interpolation formula between these asymptotic regimes by means of yet another interpolation step with respect to the parameter $\xi$. It seems convenient to separate the passive and active contributions to the effective intrinsic viscosity.

Therefore, we start with the construction of the suitable interpolation in the form

$$a_1^{\left(1\right)}\left(\xi \right)={\alpha }_0+{\alpha }_1\xi {\left(1+\frac{{\alpha }_1}{L-{\alpha }_0}\xi \right)}^{-1},$$

(76)

which is the simplest two-point Padé approximant (see, e.g., book [3]). The second approximant

$$a_1^{\left(2\right)}\left(\xi \right)=L+\left({\alpha }_0-L\right)exp\left(\frac{{\alpha }_1}{{\alpha }_0-L}\xi \right),$$

(77)

is the simplest exponential interpolation.

We next form a mixture of the two approximants using the method of weighted averaging, according to the (56),

$$a_1^{\left(p\right)}\left(\xi \right)=p{a}_1^{\left(1\right)}\left(\xi \right)+(1-p)a_1^{\left(2\right)}\left(\xi \right)$$

(78)

using the remaining condition—matching ${\alpha }_2$—to determine p:

$$\begin{array}{c}p=\frac{{\alpha }_2-u_2}{u_1-u_2},u_1=\frac{{\alpha }_1^2}{{\alpha }_0-L};u_2=u_1/2,\hfill \end{array}$$

(79)

where $u_1$ and $u_2$ are second-order terms in the expansion of $a_1^{\left(1\right)}$ and $a_1^{\left(2\right)}$ in the parameter $\xi$, respectively.

A curve drawn for pullers with positive $F_p$ (see Figure 4) [21] demonstrates that the effective intrinsic viscosity always remains positive for pullers. It approaches the aligned limit from above and another limit of small $\gamma$ from below.

The asymptotic expressions in [44] were derived in the absence of interactions between particles and yield only the intrinsic viscosity $a_1$. Instead, just as is the case of passive suspensions, we apply the ARM to the first-order expansion in concentration with $a_1^{\left(p\right)}\left(\xi \right)$.

To complete the application of the ARM, the information on critical point behavior is required. Intuitively, the densely packed regime with power-law divergence can be invoked here. Indeed, for large concentrations, the algae cells become immobilized and nearly spherical [41], apparently falling dormant in the untenable conditions of very high concentrations. Thus we can apply the techniques developed for passive spherical particles, such as factor approximant ${\mathcal{F}}^{*}\left(f\right)$ from (58) using the intrinsic $a_1^{\left(p\right)}\left(\xi \right)$ as the first-order coefficient.

The approximant is a function of both f and $\xi$, describing the two crossovers, (1) with increasing $\xi$, from zero to infinity, inside of the expression for the intrinsic viscosity, and (2) with increasing f from the active system at small f to the passive system at $f\to f_c$.

The first crossover was already formalized above by means of weighted averages of the two simple approximations. In the course of formalizing the second crossover, we have to rely on the approximants already tested above in the passive case.

For instance, the factor approximant ${\mathcal{F}}^{*}\left(f\right)$ provides a lower bound on our estimates of the EV of pullers. The very same expression as derived for the passive case is to be applied in the active case of pullers, but with $a_1$ replaced by $a_1^{\left(p\right)}\left(\gamma \right)$ for arbitrary $\gamma$. After the two crossovers are taken into account, the final expression ${\mu }^{f_a}(f,a_1^{\left(p\right)})$ should work, in principle, for arbitrary f and $\gamma$.

To assess the quality of the resulting prediction, we calculated the predicted enhancement in the effective viscosity relative to the passive spherical case. In practice, one can use the most transparent ARM approximant given by (46), with coefficients from (14). The only difference is that $a_1$ is replaced with $a_1^{\left(p\right)}$ from (78), so that it becomes ${\mu }_{1,1}^{*}(f,a_1^{\left(p\right)})$. The enhancement ratio for this approximant becomes particularly transparent

$$\frac{{\mu }_{1,1}^{*}\left(f,a_1^{\left(p\right)}\right)}{{\mu }_{1,1}^{*}\left(f,a_1\right)}=exp\left(\left(a_1^{\left(p\right)}-a_1\right)f\right).$$

(80)

To calculate the enhancement, we used the following parameters of the puller model (see Figure 3): the major semi-axis of the approximating spheroid $b=5\mathsf{\mu }\mathrm{m}$, force location characterized by $\lambda =0.5$, bacteria’s swimming speed $v=100\mathsf{\mu }\mathrm{m}/\mathrm{s}$, rotational diffusion ${\mathcal{D}}_r=0.4{\mathrm{s}}^{-1}$ taken from [41,43]. Finally, the Stokes drag law holds, i.e., $F_p=6\pi b{\mu }_{0}vN$, where an expression for a factor N can be found in [44].

Let us finally make the desired estimates of the enhancement for the typical values of parameters of the problem, following [21]. For $\gamma =5{\mathrm{s}}^{-1}$, $f\approx 0.15$ and $b/a=2.5$ (motivated by an estimate of flagellar’ length, 10–12 $\mathsf{\mu }\mathrm{m}$, see [58]), the enhancement is around $1.757$. The first-order renormalized coefficient $a_1^{\left(p\right)}$ is estimated as $7.122$, and the second-order coefficient $a_2^{\left(f\right)}=17.23$. Both numbers appeared as rather different from the corresponding numbers for passive systems.

When the enhancement ratio is averaged over various models, including the factor approximant, three ARM models and the KD model, it becomes larger and equals $1.956$. The ratio also depends weakly on $\gamma$ in agreement with the results on shear-thinning from [59]. Mind that all models involved in the averaging are based on the first-order expansion in concentration. Hence, the structure of the starting sections is dedicated in large part to such models.

It is interesting to realize, at least in principle, that an increase in the effective viscosity for pullers could have matched the corresponding decrease in the effective viscosity for pushers [44]. If ${\mathcal{D}}_r$ is decreased to $0.1$ in accordance with [42], the average enhancement ratio increases to $2.5$ in good agreement with [41] and also could become considerably larger for smaller $\gamma$. For instance, for $\gamma =1.3{\mathrm{s}}^{-1}$, we are able to estimate an increase of seven times the ratio, matching the corresponding decrease in the effective viscosity for pushers [21].

4. Critical Index for 3D Elasticity, or High-Frequency Viscosity

In the problem of finding the effective shear modulus of perfectly rigid spherical inclusions embedded at random into an incompressible matrix, we have two shear moduli, $G_1$, G, for the particles and matrix, respectively, and the $\frac{G_1}{G}$ is going to be set to infinity. The jamming regime sets in when f tends to a random close packing fraction $f_c\approx 0.637$, and is very hard to simulate numerically. Therefore, the value of theoretical approaches to this limit is increasing.

The elasticity problem is treated by analogy to the problem of high-frequency effective viscosity of a suspension [60]. In the analogous case of suspension, consisting of hard spheres immersed in a Newtonian fluid of viscosity ${\mu }_0$, only hydrodynamic interactions between pairs of suspended particles are considered.

For such suspension (or elastic composite), which is described by the ratio of the effective shear modulus $G_e\equiv {\mu }_e$ to that of the matrix, and for small f, the second-order truncation is available [35],

$$\begin{array}{c}{G}_e\left(f\right)=1+a_{1}f+a_2{f}^2+O\left(f^3\right),\\ a_1=\frac{5}{2},a_2=5.0022.\end{array}$$

(81)

The coefficient $a_2$ was also estimated earlier by Batchelor and Green, who gave $a_2=5.2$ [60]. In the vicinity of the 3D threshold, the effective modulus characteristically behaves as a power-law

$$G_e\left(f\right)\simeq A{(f_c-f)}^{-\mathcal{S}},$$

where $f_c\approx 0.637$, $\mathcal{S}\approx 1.7$[30,33]. Here $\mathcal{S}$ stands for the critical index for 3D (super)elasticity, analogous to the critical index for the viscosity when only hydrodynamic interactions are considered.

The idea of application of optimal conditions in the space of approximations was put forward by V.I. Yukalov [61,62], and the minimal difference condition was formulated in [63]. The minimal derivative (sensitivity) condition was discussed in [64]. Kleinert also employed minimal derivative conditions in the course of his variational-perturbation method [65]. Technically, it is very difficult to extract expressions for the physical quantities in a closed form by Kleinert’s methods.

4.1. Original Suggestion

Let us recapitulate the original suggestion for the calculation of the critical index [66,67]. In the case of a critical quantity divergent at the critical point, one has to apply a resummation technique to the inverse series

$${\left[G_e\left(f\right)\right]}^{-1}\simeq 1-a_{1}f+(a_1^2-a_2)f^2\approx 1+\widetilde{a_1}f+\widetilde{a_2}{f}^2,$$

(82)

where $\widetilde{a_1}=-5/2$, $\widetilde{a_2}=1.2478$. Here, we construct the simplest pair of approximations given below, with m being a control parameter,

$$R_1^{*}\left(f\right)={\left(1-\frac{{\tilde{a}}_1}{m}f\right)}^{-m},R_2^{*}\left(f\right)=1+\widetilde{a_1}f{\left[1-\frac{\widetilde{a_2}f}{\widetilde{a_1}(1+m)}\right]}^{-(1+m)}.$$

(83)

The, we set $m_2=m_1\equiv m,$ to have the same behavior at infinity.

Let us have the two solutions (83) differ minimally in the vicinity of a critical point. The minimal difference condition between the two approximations is reduced to the condition on m,

$$\left|1+\widetilde{a_1}{f}_1^c\left(m\right){\left[1-\frac{\widetilde{a_2}{f}_1^c\left(m\right)}{\widetilde{a_1}(1+m)}\right]}^{-(1+m)}\right|,$$

(84)

where the threshold as a function of m,

$$f_1^c\left(m\right)=\frac{m}{\widetilde{a_1}},$$

(85)

was found from the first-order approximant. With the concrete parameters pertaining to the problem, we calculate $m=-1.57646$. The corresponding critical index $\mathcal{S}=-m$ is found automatically.

4.2. Conditions on Thresholds

One can also select the iterated root as the second approximant, so that

$$\begin{array}{c}{\mathcal{R}}_1^{*}\left(f\right)={\left(1+\frac{a_1}{m}f\right)}^m,{\mathcal{R}}_2^{*}\left(f\right)={({({\mathcal{P}}_1\left(m\right)f+1)}^2+{\mathcal{P}}_2\left(m\right)f^2)}^m,\end{array}$$

(86)

where

$${\mathcal{P}}_1\left(m\right)=\frac{a_1}{2m},{\mathcal{P}}_2\left(m\right)=\frac{-2{a_1}^{2}m+{a_1}^2+4a_{2}m}{4m^2}.$$

The second-order approximant gives an explicit expression for the threshold,

$$f_2^c\left(m\right)=\frac{\sqrt{2{a_1}^2{m}^3-{a_1}^2{m}^2-4a_2{m}^3}+a_{1}m}{{a_1}^{2}m-{a_1}^2-2a_{2}n}.$$

The same threshold is expected for each approximation, i.e.

$$f_1^c\left(m\right)=f_2^c\left(m\right).$$

(87)

From (87), find the stabilizer m,

$$m=\frac{{a_1}^2}{{a_1}^2-2a_2}.$$

The critical index is determined by the control parameter, $\mathcal{S}=1.66471$.

4.3. Conditions on the Critical Index

Let us consider the value of a threshold as a parameter $f_0$. It could be introduced into consideration through the Euler transformation [4,68],

$$z=\frac{f}{f_0-f}\iff f=\frac{z{f}_0}{z+1},$$

applied to the original truncation (81).

Let us begin by deriving the two low-order approximations,

$${\mathcal{R}}_1^{*}\left(z\right)={\left(1+p_1\left(f_0\right)z\right)}^{m_1\left(f_0\right)},{\mathcal{R}}_2^{*}\left(z\right)={({({\mathcal{P}}_1\left(f_0\right)z+1)}^2+{\mathcal{P}}_2\left(f_0\right)z^2)}^{m_2\left(f_0\right)/2}.$$

(88)

After such a task is accomplished, it remains to minimize the difference

$$\left|m_2\left(f_0\right)-m_1\left(f_0\right)\right|,$$

with respect to the control parameter $f_0,$ where $m_1\left(f_0\right)$ and $m_2\left(f_0\right)$ are explicit formulas for the critical index obtained from the two approximations. Their difference can be minimized at a calculated value of threshold $f_0=0.656746$. The first-order approximant

$${\mathcal{R}}_1^{*}\left(f\right)={\left(\frac{1+1.01373f}{0.656746-f}\right)}^{1.61963},$$

gives the critical index $\mathcal{S}=1.61913$.

4.4. Conditions on the Amplitude

Assume now that we know the correct value of threshold $f_c=0.637$ in advance, and try to use the knowledge in the course of estimating the index. Let us apply the Euler transformation (7), with known threshold $f_c$. The original truncation is transformed to the following form,

$$\mathcal{R}\left(z\right)=1+b_{1}z+b_2{z}^2+O\left(z^3\right)b_1=1.5925,b_2=0.437238.$$

(89)

The lower-order approximation ${\mathcal{R}}_1^{*}$ can be written down as follows,

$${\mathcal{R}}_1^{*}={\left(1-\frac{b_{1}z}{{\mathbf{s}}_1}\right)}^{-{\mathbf{s}}_1}.$$

(90)

In order to construct a second approximation, we leave the constant outside of the renormalization procedure and obtain

$${\mathcal{R}}_2^{*}=1+b_{1}z{\left[1-\frac{b_{2}z}{b_1(1+{\mathbf{s}}_2)}\right]}^{-(1+{\mathbf{s}}_2)}.$$

(91)

Let us require that both approximations (90) and (91) have the same power-law behavior as $z\to \infty ,$ so that

$${\mathbf{s}}_2={\mathbf{s}}_1=\mathbf{s}.$$

The parameter $\mathbf{s}$ has to be deduced from the two approximants. Require now the fastest convergence of the two approximations, and impose the minimal difference condition on critical amplitudes. The solution to the minimization problem [4,68] is found so that $\mathcal{S}=1.60483$. As usual, more details can be found in [4,68], where the case of elastic modulus is explained.

4.5. Minimal Derivative

We also may start with the iterated root employed for the second-order approximant,

$$\begin{array}{c}{\mathcal{R}}_2^{*}\left(z\right)={({({\mathcal{P}}_1\left(\mathbf{s}\right)z+1)}^2+{\mathcal{P}}_2\left(\mathbf{s}\right)z^2)}^{\mathbf{s}},\end{array}$$

(92)

where

$${\mathcal{P}}_1\left(\mathbf{s}\right)=\frac{b_1}{2\mathbf{s}},{\mathcal{P}}_2\left(\mathbf{s}\right)=\frac{-2{b_1}^2\mathbf{s}+{b_1}^2+4b_2\mathbf{s}}{4{\mathbf{s}}^2}.$$

The control parameter $\mathbf{s}$ may be interpreted as the “critical” index in the vicinity of the (quasi)threshold $f_e^c\left(\mathbf{s}\right)$. The second-order approximant (92) allows for an explicit expression for the quasi-threshold

$$f_e^c\left(\mathbf{s}\right)=\frac{\sqrt{-{\mathcal{P}}_2\left(\mathbf{s}\right)}-{\mathcal{P}}_1\left(\mathbf{s}\right)}{{\mathcal{P}}_1{\left(\mathbf{s}\right)}^2+{\mathcal{P}}_2\left(\mathbf{s}\right)}.$$

(93)

Let us ask for an independence of the critical index at infinity on the position of a quasi-threshold. One can find the value of $\mathbf{s}$ from the minimal derivative (sensitivity) condition imposed on the quasi-threshold,

$$\frac{\partial f_e^c\left(\mathbf{s}\right)}{\partial \mathbf{s}}=0,$$

and we found the critical index $\mathcal{S}=1.788$ [4,68].

In summary, we found five fairly close estimates for the critical index $\mathcal{S}$. Their average equals $\mathcal{S}=1.65\pm 0.13$, which is in fair agreement with the experiment [30]. Finer details and explanations can be found in [4,68].

4.6. Comment on the Method of Coherent-Anomaly

It is interesting to compare our results for the critical index with other approaches with similar objectives, such as the coherent-anomaly method (CAM) of Suzuki. CAM derives the critical index from a couple of systematically obtained approximations for the sought quantity, whose critical points converge to the known true value of $f_c$. The physical idea behind initial CAM is to first take into account the diluted regime (long-range interactions) and to correct it by taking into account progressively larger clusters with shorter-range interactions. At each step, the problem becomes more difficult, but the result is still a mean-field approximation for the critical index with the expected critical index $\mathcal{S}=1$ for any finite cluster.

Assume that the problem of finding the effective property is solved analytically for a number of finite clusters of varying size, including all the interactions between the growing number n of particles at each step. One can think about the regular or randomly generated clusters or about the snapshots of the real material. In application to elasticity, the problem would consist of the numerical extrapolation of the results for the effective moduli of such clusters to the infinite cluster as n tends to infinity.

Is it possible to move beyond the mean-field approximation within the framework of the cluster approach? The answer is “no” if only finite-size clusters are considered without additional efforts to make an inference on the case of infinite clusters. Suzuki suggested that additional information about the critical index is hidden in the parameters of two successive approximations [69]. In order to extract the information, one has to find the approximate thresholds $f_{c1}$ and $f_{c2}$ and the amplitudes $A_1$ and $A_2$, i.e., the derivatives calculated in the vicinity of the approximate thresholds in each approximation. Then, the critical index value, different from the mean-field value, could be estimated. With increasing complexity of the approximation, the value of the threshold is supposed to approach the true threshold, and the value of the amplitude is supposed to increase, reflecting the so-called coherent anomaly.

The CAM is built in the same spirit as any renormalization group (RG), but the flow of the sought quantity is considered with respect to the approximation parameter, discrete or continuous. The critical amplitude of the mean-field approximation is considered the sought physical quantity. The value of amplitude should be known for at least two successive approximations. In addition, for the amplitude, it is most important that one should introduce a scaling hypothesis, which postulates that the amplitude diverges as a power-law with the index of $\tilde{s}$ in the vicinity of the fixed point with exact threshold $f_c$. The divergence rate should give a non-mean-field contribution to the critical index on top of the mean-field value.

The scaling hypothesis, in turn, should be expressed in terms of the parameter characterizing proximity to the true critical point, or more generally, a fixed point of various RG transformations. To this end, one should know how such a parameter varies from approximation to approximation. It is enough to know the derivative of such transformation in the vicinity of the fixed point. When the ratio of amplitudes for two successive approximations is expressed in two different ways, one arrives at the formula for the critical index to be found as the solution of a simple equation. CAM recommends the simplest and most natural choice of the approximate threshold $f_c$ as the parameter characterizing the different approximations, and the derivative is simply expressed in finite differences between the two thresholds and the exact threshold. From the approximate thresholds and known amplitude ratios, one can express the critical index.

The final expression for the critical index $S=1+\tilde{s}$ becomes explicit and simple, as $\tilde{s}$ is expressed below

$$\tilde{s}=\frac{log\left(\frac{A_1}{A_2}\right)}{log\left(\frac{f_c-f_{c2}}{f_c-f_{c1}}\right)}.$$

(94)

Here, the value of unity comes from the mean-field approximation and $\tilde{s}$ from the unaccounted interactions. This formula written in appropriate variables is standard for most RG calculations of the critical indices [70], except for the approaches of the type of [66], where the critical index is treated as a control parameter and found from some sort of additional condition of a non-perturbative nature.

We applied the two methods suggested by Suzuki to generate systematic mean-field approximations from the series. The first method of [71] is based on an analysis of the series and by studying the zeroes for the first two polynomials of the power-series of the inverse of the relevant physical quantity, giving $f_{c1}=0.4$ and $f_{c2}=0.552187$. Then, the mean-field amplitudes $A_1=2.5$ and $A_2=1.12196$ can be found near each of the points.

Because we are working with inverse quantities, the value of amplitude decreases with an increasing number of terms. The critical exponent is estimated on the basis of the general CAM theory [69,72] and formula (94) with changed sign, $\mathcal{S}=1+\frac{log\left(\frac{A_1}{A_2}\right)}{log\left(\frac{f_c-f_{c1}}{f_c-f_{c2}}\right)},$ with the result $\mathcal{S}\approx 1.7797$, for $f_c=0.637$.

In the second method of [73], one should first construct the continued fractions, approximant asymptotically equivalent to the original series. Then, find the critical point from the position of the simple pole and determine the critical coefficient/amplitude as the corresponding residue so that

$$f_{c1}=0.4,f_{c2}=0.49978,A_1=0.4,A_2=0.62445.$$

The information extracted from the two continued fractions is then combined with the CAM theory to estimate fractional critical exponents [69,72] from formula (94), with the result $\mathcal{S}\approx 1.8151$. The value of the critical coefficient now is increasing with an increasing number of terms, as is anticipated in the framework of CAM, and formula (94) is fully applicable.

The techniques of Padé approximants and their various modifications are rather popular, and their advantages and limitations are understood. Thus, first applying the diagonal Padé approximants or continued fractions as the zero approximation for the critical point should be quite natural. However, in an unexpected turn, the value of a critical index can also be found from such a pedestrian approach. However, it is given only implicitly, as it is hidden or scattered within the different Padé approximants.

The CAM estimates agree with other estimates from this section and are close to the upper bound. Such an estimated critical index could be employed as an initial guess/input for some other technique, which would operate with the critical index explicitly. Such a technique would calculate a correction to the CAM value and reconstruct corresponding improved approximants, giving the expression for the sought quantity for all concentrations.

5. Sedimenatation and Particle Mobility

The sedimentation problem is concerned with a suspension movement under gravity. The suspension is built of small rigid spheres with random positions, and they are falling through Newtonian fluid under gravity. A homogeneous mixture of solid particles and a fluid stands in a container. The rate of particles settling under gravity depends on their concentration, which is dependent on the hydrodynamic interaction between particles.

The quantity of interest is called sedimentation velocity U. It is nothing else but the averaged velocity of suspended particles. It is usually normalized with respect to the velocity $U_0$ of a single particle movement in the absence of any other particles. The ratio $U/U_0$ is also called the collective mobility or sedimentation coefficient. In sedimentation experiments, two different scenarios are played out. The first scenario corresponds to the phase transition [74], but different scenarios with metastable phase extension up to RCP [75] can be observed as well.

Below, we mainly follow the discussion presented in [4,68], present some empirical and semi-empirical observations, and develop formulas based on some low-order expansions for the sedimentation rate. We also provide some estimates for the critical properties of the sedimentation rate. Based on Stokes relation $U_0\sim \frac{1}{{\mu }_0}$, where $\mu$ is an ambient fluid viscosity, one might expect the sedimentation rate to be inversely proportional to the effective viscosity of the suspension

$$U\sim \frac{1}{{\mu }_e\left(f\right)},$$

(95)

One can think that the critical index value for the sedimentation rate is close to the critical index for effective viscosity. Various studies [76,77,78,79,80,81] suggest a critical behavior for the sedimentation rate in the form of formula (7),

$$\frac{U}{U_0}\sim {(1-f)}^{\beta },$$

(96)

with a positive critical index in the range from $\beta \approx 3$ to $\beta \approx 5.5$. Various formulas for the sedimentation rate can be found in [4]. For small f, the following single-term expansion formula was found in [77]

$$\frac{U}{U_0}=1-6.55f+O\left(f^2\right).$$

(97)

The value of $a_2$ will be estimated below, following [4]. It turns out that the estimated value is in rather good agreement with some other microscopical, theoretical estimates.

We are interested in accurate expressions satisfying (97), with power-law behavior (96). Let us follow the method of derivation, which leads to formula (19). Let us find for $U/U_0$ the simplest approximant with fixed position of singularity and given the critical index [4],

$$\frac{U}{U_0}={\left(1+\frac{2.18333f}{1-f}\right)}^{-3}.$$

(98)

From formula (98), one can find the approximate value of the second-order coefficient, i.e., $a_2\approx 22.05$. One can also estimate $a_3\approx -53.425$.

Another simple formula can be obtained similarly to the derivation of ([21], eq. 3.6),

$$\frac{U}{U_0}=e^{-3.55f}{(1-f)}^3,$$

(99)

and its expansion at small f gives a close value for $a_2\approx 19.95$. Both estimates for $a_2$, especially the former, agree with the value of $a_2\approx 21.918$ calculated in [82].

In reverse, accepting the theoretical estimate of $a_2=21.918$, we can find

$$\frac{U}{U_0}={\left(1+\frac{2.15061f}{1-f}\right)}^{-3.04379},$$

(100)

and estimate that $\beta \approx 3.044$. The latter value is exactly the result of the application of the diff-log Padé method from [4,68], after Euler transformation (7), from the finite interval to semi-infinity. After the transformation, one can simply apply the technique sketched in the end of Section 2.1. The new formula (98) agrees better in comparison with the simulation results of [80].

We conclude that the value of $a_2\approx 21.918$ from the [82] is reasonable. We also employed it for more detailed estimates of the critical properties. Following minor modifications to the same approach of the five methods of application, as in Section 4, we found $\beta =3.003\pm 0.323$ [68]. More details are presented below in Section 5.1. Keep in mind that in [33], we estimated the critical index $\beta =3.3049\pm 0.26.$

Practitioners also respect the formula suggested by Barnea–Mizrahi [83],

$$\frac{U}{U_0}=\frac{{(1-f)}^2}{(\sqrt[3]{f}+1)exp\left(\frac{5f}{3(1-f)}\right)}.$$

(101)

It works quite well in the region of high concentrations but gives different results as $f<0.2$ compared to the numerical data of [80].

5.1. Direct Estimates of the Critical Indices

In reverse, when the coefficient $a_2$ calculated in [82] is known in advance, so that for small f

$$\frac{U}{U_0}=1-6.55f+21.918f^2+O\left(f^3\right),$$

one can estimate the value of a critical index at $f_c=1$. Such an estimate is independent of any other experimental or theoretical considerations [68]. Let us again employ the special resummation techniques, just as in the previous section, for the shear modulus. In this case, for obtaining estimates of the critical index $\beta$ for the sedimentation rate, remember that the strategy consists of employing various conditions for the fixed point in order to compensate for a shortage of the coefficients. We present below only the main results obtained in [68].

The method of Section 4.1 gives the critical index $\beta =2.966$.

The method of Section 4.2 produces an unreasonably large value and is discarded.

The method of Section 4.3 brings the critical index $\beta =2.863$.

The minimal difference condition of Section 4.4 gives the index $\beta =2.577.$

The minimal derivative condition of Section 4.5 leads to the critical index $\beta =3.566$.

The index can also be obtained with the diff-log Padé method of [4,68], with the result $\beta =3.044$.

The average taken among all estimates gives the value of $\beta =3.003\pm 0.323$.

The formulas for the sedimentation rate U presented above are applicable up to $f\approx 0.5$ before the freezing point for the hard-spheres suspension.

In the case of a metastable branch of the suspension, the following expression was derived in [4,68],

$$\frac{U}{U_0}\left(f\right)={\left(1+\frac{4.17235f{\left(1+\frac{1.4814f}{0.637-f}\right)}^{0.6937}}{0.637-f}\right)}^{-1}.$$

(102)

In the latter case, we looked for a formula to continue up to the jamming threshold of $f_c\approx 0.637$. The critical index $\beta \approx 1.7$ was found for such a threshold in [4,68]. The latter value of $\beta$ appeared as close to the estimates for the critical index of the high-frequency effective viscosity due to hydrodynamic interactions, obtained in the preceding Section.

6. Non-Local Diffusion in Critical Liquids

The definition of a non-local diffusion coefficient in liquids and suspensions dwells on some other quantities, which are briefly discussed below. First, one has to introduce the space- and time-dependent local number density

$$\rho (\mathbf{r},t)=\sum _{i=1}^N\delta \left(\mathbf{r}-{\overrightarrow{r}}_i\left(t\right)\right),$$

where N stands for the number of particles.

Second, consider the intermediate scattering function, which is the auto-correlation function of the Fourier components ${\rho }_q\left(t\right)$ [84],

$$F(q,t)=\frac{1}{N}\langle {\rho }_q\left(t\right){\rho }_{-q}\left(0\right)\rangle ,$$

where the averaging is performed over a statistical ensemble. $F(q,t)$ depends on time t and on the wave-number q.

Even more convenient for experimental and theoretical purposes is the dynamic structure factor $S(q,\omega )$[85]. It is defined as the Fourier-transform of the $F(q,t)$ in the temporal domain

$$S(q,\omega )={\int }_{-\infty }^{\infty }F(q,t)exp\left(ti\omega \right)dt.$$

The celebrated $S(q,\omega )$ is directly measurable in light scattering experiments. The density fluctuations spectrum ${\omega }_c$ corresponds to the simple pole $\omega =-i{\omega }_c,$ in $S(q,\omega )$, dependent on frequency $\omega$ and on the wave-number q. To obtain the non-local, q-dependent diffusion coefficient, one has to simply divide ${\omega }_c$ by $q^2$ [86].

Near the critical point, temperature is measured by the ratio $|\frac{T-T_c}{T_c}|$. We are interested here only in long-time relaxation properties, dependent on the wave-number q. They are characterized by the characteristic relaxation time $\tau$ of the density fluctuations. Such fluctuations are the inverse of the ${\omega }_c$.

The dynamic scaling hypothesis conjectures that the characteristic frequency can be expressed in the following form

$${\omega }_c(\zeta ,q)=q^z\Omega \left(\frac{q}{\zeta }\right).$$

Here, the dynamical critical index z is introduced, and $\zeta$ stands for the inverse correlation length [38]. In fact, we are looking for an approximation for the function $\Omega$ and would like to estimate the index z. The length scale ${\zeta }^{-1}$ is intermediate between the atomic and macroscopic scales.

The behavior of ${\omega }_c$ can be reconstructed approximately from some general asymptotic information [38] using the crossover approach of [24]. The scheme was extended and developed further in the book [4]. Crossovers are typical of the thermodynamics of a supercritical state [87].

In the hydrodynamic regime $\frac{q}{\zeta }\ll 1$, for isotropic fluid, we have an expansion [38],

$${\omega }_c(\zeta ,q)\simeq D_T{q}^2\left(1+B{\left(\frac{q}{\zeta }\right)}^2+\dots \right),B>0(\frac{q}{\zeta }\ll 1).$$

(103)

In the fluctuation regime $\frac{q}{\zeta }\gg 1$, another expansion is considered,

$${\omega }_c(\zeta ,q)\simeq A{q}^z\left(1+A^{\prime }{\left(\frac{\zeta }{q}\right)}^2+\dots \right),(\frac{q}{\zeta }\gg 1),$$

(104)

The expansion is expounded in [38]. Here, $D_T$ is thermal diffusivity. The value of B is estimated as $B=3/5$ [38,88]. Assume also that the value of A is known.

The simple analytical expression for the characteristic frequency for arbitrary $\frac{q}{\zeta }$ was suggested [24],

$${\omega }_c(\zeta ,q)=D_T{q}^2{\left(1+C{\left(\frac{q}{\zeta }\right)}^2\right)}^{\frac{z}{2}-1},$$

(105)

where C is constant, and

$$C={\zeta }^2{\left(\frac{A}{D_T}\right)}^{\frac{2}{z-2}}.$$

(106)

If we now substitute the dependencies of $D_T\sim {ϵ}^{\gamma -a}$ and $\zeta \sim {ϵ}^{\nu }$ into the expression for C, then the well-known relation

$$z=2+\frac{(\gamma -a)}{\nu },$$

(107)

we simply recover the known relation between the dynamic critical index z and three other critical indices $\gamma \approx 1.239$, $\nu \approx 0.63$ and a. Formula (107) represents one of the central results of the dynamical scaling hypothesis [38]. Remember that $\gamma$ stands for the critical index for heat capacity, while a stands for the critical index for thermal conductivity. Index $\nu$ is the correlation length critical index. The combination

$$a^{\prime }=\gamma -a\approx 0.67$$

stands for the critical index for thermal diffusivity.

From (106), it follows then that

$$A=A\left(ϵ\right)=D_T{\zeta }^{2-z}{C}^{\frac{z}{2}-1},$$

in such a way that requirement $A\left(0\right)$=const is satisfied when the relation between indices (107) is valid.

The index $a^{\prime }$ is determined mainly by the concentration (long time) diffusion coefficient (see, e.g., [89]), which in turn, can be estimated from the analog of Stoke’s formula [4]. The value of $a^{\prime }$ is dominated by the correlation length index, provided the effective viscosity $\mu$ remains finite at the critical point as $ϵ=0$ [4]. Therefore, $a^{\prime }\approx \nu$.

Using the scaling relation (107), the dynamic critical index can be crudely estimated as $z\approx 3$. However, one can also take into account a correction to z, originating from weakly divergent viscosity. If such divergence is characterized by the index $x_{\mu }\approx 0.068$ (measured in units of the critical index $\nu$), then

$$z\equiv z\left(x_{\mu }\right)=2+(1+x_{\mu }).$$

The latter formula means that it takes longer for short wavelength density fluctuations to relax if compared with a long wavelength regime. Different estimates for $x_{\mu }$ can be found in [90].

The unknown C can be estimated from the asymptotic equivalence with the expansion (103) in the hydrodynamic regime. For instance, $C\approx 1.2$, with typical B and z. To calculate the thermal diffusivity $D_T$, one can simply use the diffusion coefficient approximated by the Stokes formula, as explained in [4,86,89].

The non-dimensional part of ${\omega }_c$ can be approximated already by a simple self-similar root approximant

$${\omega }_{nd}=\sqrt{1+1.2x^2},$$

(108)

where $x=\frac{q}{\zeta }$. Its asymptotic behavior as $x\to 0$

$${\omega }_{nd}\left(x\right)\simeq 1+0.6x^2-0.18x^4+\dots ,$$

(109)

respects the coefficient B, while as $x\to \infty$

$${\omega }_{nd}\left(x\right)\simeq 1.09545x+0.456435x^{-1}-0.0950907x^{-3}.$$

(110)

Both expressions are in qualitative agreement with the expressions (103), (104). Expansions (109), (110) are designed to obey the dynamic scaling in all orders.

The non-dimensional part of ${\omega }_c$

$${\omega }_{kaw}\left(x\right)=\frac{3\left(\left(x^3-\frac{1}{x}\right){tan}^{-1}\left(x\right)+x^2+1\right)}{4x^2},$$

(111)

was derived by Kawasaki [88]. Its asymptotic behavior as $x\to 0$

$${\omega }_{kaw}\left(x\right)\simeq 1+\frac{3x^2}{5}-\frac{x^4}{7}+\frac{x^6}{15},$$

(112)

is qualitatively correct in all orders. It is in agreement with the dynamic scaling in the hydrodynamic limit (103).

While as $x\to \infty$

$${\omega }_{kaw}\left(x\right)\simeq \frac{3\pi }{8}x+\frac{1}{x^2}-\frac{3\pi }{8x^3}.$$

(113)

However, the second expansion (113) agrees only with the leading term in (104). This means that it violates the dynamic scaling by means of the $x^{-2}$ term in (113). The term $O\left(x^{-1}\right)$ is altogether absent. We conclude that the Kawasaki approximation [88] is correct only as far as the leading term is considered in the limit-case of $x\to \infty$.

In the book [4], we attempted to find some form of approximant extending (108) based asymptotically on a few terms from the correct Kawasaki limit $x\to 0$ and the leading term as $x\to \infty$. Such approximants should be free from the terms violating dynamic scaling. In other words, we thrive on restoring dynamic scaling in all orders. To this end, we have to come up with some additional considerations.

Already, the simplest factor approximants give a reasonable approximation to the Kawasaki function while restoring the dynamic scaling in all orders [4]. Based on this observation, one can suggest a simple factor approximant to the non-dimensional part of ${\omega }_c$, generalizing the simplest forms presented above to an arbitrary z ($x_{\mu }$). In other words, the second-order factor approximant satisfies just one non-trivial coefficient from (112) and is designed to give the correct critical index,

$$\begin{array}{c}{\mathcal{F}}_2^{*}(x,x_{\mu })={\left(1+\frac{6x^2}{5(x_{\mu }+1)}\right)}^{\frac{x_{\mu }+1}{2}}.\hfill \end{array}$$

(114)

The higher-order approximant asymptotically obeys the three non-trivial coefficients from (112) and also gives the correct critical index z. It also has a closed form but with expressions for the parameters $P_1\left(x_{\mu }\right)P_2\left(x_{\mu }\right),m_1\left(x_{\mu }\right)$ given explicitly in [4],

$$\begin{array}{c}{\mathcal{F}}_4^{*}(x,x_{\mu })=\left(x^2{\mathcal{P}}_1\left(x_{\mu }\right)+1\right){}^{m_1\left(x_{\mu }\right)}\left(x^2{\mathcal{P}}_2\left(x_{\mu }\right)+1\right){}^{\frac{1}{2}\left(x_{\mu }+1\right)-m_1\left(x_{\mu }\right)}.\hfill \end{array}$$

(115)

Approximants (114) and (115) respect the dynamic scaling by design. Of course, they can also be used with a non-zero $x_{\mu }$. In such cases, the Kawasaki function is not rigorously applicable. Remember that in developing the factor approximations, we rely on the fact that the Kawasaki approximation works well at small x, at least in leading orders [38]. However, the Kawasaki function is often used as a part of the formulas that account for the correct value of the critical index z. For instance, the following approximation was developed in [91],

$${\omega }_{pp}\left(x\right)={\left(x^2+1\right)}^{x_{\mu }}{\omega }_{kaw}{\left(x\right)}^{1-x_{\mu }},x_{\mu }=0.07.$$

(116)

Another formula was explored in [86,92,93],

$${\omega }_{bsbg}\left(x\right)={\omega }_{kaw}\left(x\right){\left(1+\frac{1}{4}{x}^2\right)}^{x_{\mu }/2}{x}_{\mu }=0.068.$$

(117)

In Figure 5, factor approximants (114) and (115) are compared with the two non-monotonous approximations based on the Kawasaki function and with reference to the monotonous Kawasaki function itself. They deviate from the Kawasaki (111) in less than two other semi-empirical formulas. Such a difference is potentially above the typical experimental error of $10\%$ [89].

7. Effective Conductivity of a Classical Analog of Graphene

For the subsequent study, it is important to remember that graphene consists of atoms tightly packed into a 2D honeycomb lattice. Graphene has very good electric conductivity, at par with silver and copper. On the other hand, silver (or copper) inclusions, embedded into a very weakly conducting matrix (vacuum, air, water, polyfoam), can be considered perfect conductors [94].

Consider a continuum analog of graphene [68,94]. To this end, one can introduce and study a regular honeycomb array of perfectly conducting (superconducting) disks. The disks are embedded into a matrix with a conductivity normalized to unity. With respect to the matrix conductivity, one can set the inclusions conductivity to infinity. The relevance for establishing the analogy case corresponds to the volume fraction $f\to f_c$. In such a limit-case, one can expect that the effective conductivity of the array goes to infinity.

If some vacancies are introduced at random within the continuum model, the threshold $f_c$, or maximal possible volume fraction of disks in a regular case, will be lowered to the value of $f_c\left(p\right)$, where the parameter p is the number of vacancies per cell and $f_c\left(0\right)=f_c$. The effective conductivity ${\sigma }_e$ is expected to diverge with some critical index $\mathsf{s}\left(p\right)$. For the sake of analogy, with graphene, it only makes sense to consider the limit $f\to f_c\left(p\right)$ [68,94].

The expansion for small f is possible to obtain by means of the Monte Carlo simulations combined with analytical results from the method of basic structural sums, as described in [94], leading to the polynomials representing effective conductivity at small and moderate f. The details of the methodology of deriving the effective conductivity of finite samples for various p can be found in [68,94]. A complete set of polynomials for $p=36$ can be found in [68].

For instance, in some particular cases, the following truncations were computed,

$$\begin{array}{c}{\sigma }_e\simeq 1+2f+2f^2+2.085728259543634f^3+4.503578882551432f^4+\hfill \\ 7.022350274464803f^5+9.889873882233639f^6,p=2;\hfill \\ \\ {\sigma }_e\simeq 1+2f+2f^2+2.4267342319989553f^3+5.952924825710871f^4+\hfill \\ 10.215408230511798f^5+17.000009114788686f^6,p=8.\hfill \end{array}$$

(118)

Smoothing of the polynomials was performed following the methodology explained in [94] in great detail. Concrete sets of polynomials for smaller samples can be found in [94]. Polynomials (118) can be seen as evolving smoothly from the expression (121) for the regular case with increasing p.

The critical index will be found by extrapolation from the small-f polynomials to $f\to f_c\left(p\right)$. The threshold is conveniently given by formula (120). In the vicinity of the threshold, the following family of power-laws is considered

$${\sigma }_e\left(f,p\right)\sim {\left(f_c\left(p\right)-f\right)}^{-\mathsf{s}\left(p\right)}.$$

(119)

Here, the superconductivity critical index $\mathsf{s}$ depends on p. We also expect that as $p=0$, the index should be equal to $1/2$.

The threshold $f_c\left(p\right)$ also depends on p. It is given by the compact formula

$$f_c\left(p\right)=\frac{\pi }{3\sqrt{3}}\frac{p_{max}-p}{p_{max}},$$

(120)

where $f_c\left(0\right)=\frac{\pi }{3\sqrt{3}}$ is the maximum volume fraction attainable for the honeycomb lattice, and the maximal value of p is denoted as $p_{max}$. In simulations of [68], $p_{max}=36$. When the threshold is known, we would first apply the following Euler transformation (7),

$$z=\frac{f}{f_c\left(p\right)-f}\iff f=\frac{z{f}_c\left(p\right)}{z+1}$$

to the original series. The transformation maps the segment to a half-line.

It is instructive to first study the regular case with no defects. The small-f polynomial corresponding to the case of $p=0$ has the following form,

$${\sigma }_e\simeq 1+2f+2f^2+2f^3+4.14933f^4+6.29865f^5+8.44798f^6,$$

(121)

calculated exactly from the formulas of [3], applied to the case of a honeycomb array of disks. Let us calculate the critical index $\mathsf{s}$ applying the technique suggested in Section 2.1.1. After the diff-log transformation, taking an inverse and Borel summation of Section 2.1.1, we arrive at the following sequence of approximations dependent on the number of terms k involved in the construction of the iterated roots,

$${\mathsf{s}}_1=0.6046,{\mathsf{s}}_2=0.5617,{\mathsf{s}}_3=0.5372,{\mathsf{s}}_4=0.523,{\mathsf{s}}_5=0.5142.$$

The latter sequence can be continued beyond the fifth order if we also keep taking the inverse of the sixth order term, i.e., ${\mathsf{s}}_6=0.5089$. The results agree rather well with the expected value of $0.5$. The latter estimate seems better than the result $\mathsf{s}\approx 0.519$ of [68], which was obtained by the technique of the paper [33], as already employed in Section 2.1.2. Note that the method of [31], based on the Borel summation and factor approximants [31], failed to produce a convergent sequence.

Similarly, for the case of $p=2$, we obtain the following sequence

$${\mathsf{s}}_1=0.571,{\mathsf{s}}_2=0.5306,{\mathsf{s}}_3=0.507,{\mathsf{s}}_4=0.4928,{\mathsf{s}}_5=0.4837,{\mathsf{s}}_6=0.4779,$$

The latter estimate agrees with the result $\mathsf{s}\approx 0.4814$ of [68] for the same disordered case.

Again, similarly for the case of $p=8$, we obtain the following sequence

$${\mathsf{s}}_1=0.4702,{\mathsf{s}}_2=0.4367,{\mathsf{s}}_3=0.4157,{\mathsf{s}}_4=0.4023,{\mathsf{s}}_5=0.393,{\mathsf{s}}_6=0.3864,$$

The latter estimate agrees with the result $\mathsf{s}\approx 0.3796$ of [68] for the same disordered case. The dependence of $\mathsf{s}$ on the degree of disorder measured in steps p of randomization achieved through random placement of vacancies on a honeycomb lattice is displayed in Figure 6. The results given by the Borel summation method in the sixth order, as explained in the examples above, are shown with squares. The results found by applying the fifth-order Borel summation method are shown with disks.

We see that both approximations give rather close results for the index for all p.

One can also plot curves for the (inverse) effective conductivity [68,94] for varying p. The dependence is characterized by the steepening of the curves with an increasing fraction of vacancies [68,94], while the threshold simultaneously shifts to smaller values. The index $\mathsf{s}$ is non-universal, smaller than $1/2$, and decreases with increasing concentration of vacancies. Remarkably, the honeycomb array of highly conducting disks has the ability to adapt to randomly placed vacancies by adjusting (lowering) the threshold and critical index, effectively retaining its quasi-regularity. One can think that analogous graphene can also withstand the infusion of vacancies while remaining a good conductor.

The dependence of the critical index on the degree of disorder defies the universality hypothesis [95]. Similarly, the critical index in the so-called eight-vertex model and equivalent Ising model with four-spin interaction [95] depends on temperature [96].

It would also be very interesting to understand if there is a way to find a continuous analog of the order to disorder abrupt transition in quasi-periodic twisted bilayer graphene composites [97]. The interesting effect of order emerging from disorder was discovered in [98]. In both situations, it could be interesting to access the change in critical properties in the course of transitions.

8. Superconductivity Critical Index of Random 3D Composite

It was argued in [4] that in the limit-case of a random array of perfectly conducting spherical inclusions embedded into a matrix with a conductivity normalized to unity, the expression for effective conductivity ${\sigma }_e\left(f\right)$, given by the classical Jeffrey formula [99], could possibly contain the wrong quadratic term in the volume fraction f.

Interesting and challenging analyses of various formulas in the theory of composites can be found in [100,101]. However, the effective medium-type approaches can be useful if applied properly for the qualitative analysis of the complex composite systems [102,103], when modeling is not feasible, when the system is designed particularly to satisfy the certain, generalized Hashin–Shtrikman bounds [104], or when the expansions are converging fast [105]. The method of basic structural sums developed by Mityushev et al. [3,4], serves as the firm foundation for obtaining polynomial truncations of the effective properties, as well as for the direct analysis of a micro-structure of nanocomposites [106].

For macroscopically isotropic composites, the following truncated expansion for effective scalar conductivity was found

$${\sigma }_e=1+3f+3f^2+4.80654f^3+O\left(f^{\frac{10}{3}}\right).$$

(122)

Jeffrey [99], in this case gives, only a quadratic term, $4.51f^2$. Mityshev et al. concluded in [4] that such a disparity is caused by a methodological error, “related to an intuitive physical treatment of the conditionally convergent integral discussed in [99]”.

With a known threshold corresponding to the random close packing (RCP), $f_c=0.637$, the conductivity is given by the power-law of the type (7),

$${\sigma }_e\simeq A{(f_c-f)}^{-\mathsf{s}}.$$

From simulations, the superconductivity critical index $\mathsf{s}$ is expected to have the value of $0.73\pm 0.01$[107]. Let us estimate the value of $\mathsf{s}$ based on asymptotic information contained in (122). Let us use the same method for calculating the critical indices as in Section 2.1.2. To this end, compose the sequences ${\Delta }_{kn}=0$, with $k=1,2$. From the series (122), one can find only three independent estimates for the critical index $\mathsf{s}$, see

Table 2. Critical indices for the superconductivity $s_k$ obtained from the optimization conditions ${\Delta }_{kn}\left({\mathsf{s}}_k\right)=0$.

${\mathsf{s}}_{\mathit{k}}$	${\mathbf{\Delta }}_{\mathit{k},\mathit{k}+1}\left({\mathsf{s}}_{\mathit{k}}\right)=0$	${\mathbf{\Delta }}_{\mathit{k}3}\left({\mathsf{s}}_{\mathit{k}}\right)=0$
${\mathsf{s}}_1$	0.725	0.721
${\mathsf{s}}_2$	0.715	0.715

. One can see that all three are close to $0.72$ and to the expected value of $0.73$. More details can be found in the book [4].

To obtain a final formula, one can simply fix the value of $\mathsf{s}=0.73$ and find the following approximant for the effective conductivity:

$$\begin{array}{c}{\mathcal{R}}_3^{*}\left(f\right)={\left({\left(\frac{2.56706f^2+2.06109f+0.405769}{{(0.637-f)}^2}\right)}^{3/2}-\frac{0.372504f^3}{{(0.637-f)}^3}\right)}^{0.243333}.\hfill \end{array}$$

(123)

From the approximant, we obtain $A\approx 1.456$. Furthermore, the fourth-order coefficient, $a_4\approx 7.48$, could be found from ${\mathcal{R}}_3^{*}\left(f\right)$.

Let us calculate the critical index $\mathsf{s}$, applying the technique suggested in Section 2.1.1. After the diff-log transformation, inversion and Borel summation, we arrive at the following sequence of approximations dependent on the number of terms k involved in the construction of the iterated roots,

$${\mathsf{s}}_1=0.725,{\mathsf{s}}_2=0.72.$$

With the index fixed to the value of $\mathsf{s}=0.72$ and applying the technique of [32] directly to the truncation (122), we obtain the following sequence of amplitudes

$$A_1=1.332,A_2=1.398,A_3=1.43,A_4=1.446.$$

The latter sequence was continued beyond the third order since we kept going after the Euler variable transformation (7) to the fourth-order term.

The method of self-similar Borel summation [31] based on Borel transformation and factor approximants, and described in more detail in Section 2.1.1, produces reasonable values

$${\mathsf{s}}_2=0.84,{\mathsf{s}}_3=0.721,$$

for the index, and in the third order gives $A\approx 1.489$. Various techniques bring very close results for the critical index and amplitude.

9. Concluding Remarks

The conclusion of the paper is straightforward. Even short polynomial expressions derived by the great scientists for the effective properties at small values of relevant parameters are rather informative, as long as their treatment is complemented with the consideration of critical properties for a non-universal critical phenomenon. Such a particular realization of the L. I. Manevitch asymptotic complementarity principle (from the 1960s) [108] is our contribution to the “Big Picture”. Moreover, the two limit-cases are so thoroughly intertwined or entangled in the course of obtaining an approximate solution, that each of them retains hidden information about the other even when they are considered separate or disentangled. Uncovering such information turns out to be possible by means of mostly analytical tools. Thus, an analytical approach to perturbative problems advocated by Feynman and developed by Bogolyubov in the refined form of an RG turns out to be fruitful in application to the effective properties.

One can think about the methodological universality in place of the physical universality of critical phenomena elaborated on by Wilson and Kadanoff [6,70].

In Section 2, we applied the methods of estimation of the coefficients in the truncated expansions based on the first-order and second-order coefficients, supplemented by the information on the maximum volume fraction and critical indices for the effective viscosity. In the three-dimensional case, the critical index is even calculated based on the two non-trivial coefficients and threshold by applying the factor approximants. The results for the index are found to be in good agreement with Brady predictions. In the three-dimensional case, an accurate crossover formula for the conductivity is brought up. They are compared with others and are in good agreement with the reference formula.

In Section 3, the formulas developed in Section 2 for passive suspensions are applied to the active suspensions of puller micro-organisms. The two crossover problems for the intrinsic viscosity dependent on activity, and the effective viscosity dependent on volume fractions are solved approximately. The quality of approximations could be evaluated by estimating the enhancement in the effective viscosity, which turned out correct by an order of magnitude for pullers.

In Section 4, we discussed the main five methods for the problem of the calculation of the critical indices and amplitudes from extremely short truncated series. The critical index for the shear modulus of random elastic composites is calculated by all five methods, and the average (and variance) is found in agreement with expectations from the experiments.

In Section 5, the critical indices for sedimentation are calculated from the short truncation by means of minimal difference and minimal derivative conditions imposed directly on the critical properties. The estimates at various thresholds agree with other theoretical estimates and real experimental results.

In Section 6, we found the crossover expressions for the non-local diffusion coefficient of a simple liquid in the vicinity of a critical point, dependent on the wave-number of density fluctuations in the liquids. The extensions and corrections to the celebrated Kawasaki formula are designed in such a way that the dynamic scaling hypothesis is satisfied in all orders.

In Section 7, we calculated the effective critical index for superconductivity dependent on the concentration of vacancies. The index is predicted to be smaller than the value of $1/2$ for the two-dimensional regular composites.

In Section 8, we calculated the superconductivity critical index of the random three-dimensional composite in good agreement with results from simulations on lattice models.

The problems to be solved are technical but difficult, consisting of improvements of the quality of expansions using some smoothing procedures, such as those applied in the derivation of the expansions in [94]. Hopefully, a more accurate truncated expansion could be obtained for the problems of finding critical indices t for conductivity in 3D and elasticity $\mathcal{T}$ in 2D and 3D cases [1,4].

The method of self-similar Borel summation could be improved further by applying it iteratively, with the number of iterations playing the role of a control. Generally speaking, the quality of Borel summation could be improved through various generalizations, and finding the most effective optimization technique, in addition to those already available [32], would help to improve the accuracy of summations in various problems of physics and applied mathematics.

Another challenge to be solved is envisaging a 2D system with a critical index for superconductivity $\mathsf{s}$ larger than $1.3$, as our work already covered the physics behind the region of $0<\mathsf{s}\le 1.3$ [3,68,94].