# Effective viscosity of a two dimensional passive suspension in a liquid crystal solvent

S. Dang, C. Blanch-Mercader, L. Berlyand

► To cite this version:

S. Dang, C. Blanch-Mercader, L. Berlyand. Effective viscosity of a two dimensional passive suspension in a liquid crystal solvent. European Physical Journal E: Soft matter and biological physics, 2024, 48 (4-5), pp.23. 10.1140/epje/s10189-025-00479-2. hal-04779444

HAL Id: hal-04779444

https://hal.science/hal-04779444v1

Submitted on 13 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Effective viscosity of a two dimensional passive suspension in a liquid crystal solvent

S. Dang,¹ C. Blanch-Mercader,²,\* and L. Berlyand³,†

¹Department of Mathematics, The Pennsylvania State University, USA

²Laboratoire Physico-Chimie Curie, Institut Curie, Université PSL, Sorbonne Université, CNRS UMR168, F-75248 Paris, France

³Department of Mathematics and Huck Institute for Life Sciences, The Pennsylvania State University, USA

Suspension of particles in a fluid solvent are ubiquitous in nature, for example, water mixed with sugar or bacteria self-propelling through mucus. Particles create local flow perturbations that can modify drastically the effective (homogenized) bulk properties of the fluid. Understanding the link between the properties of particles and the fluid solvent, and the effective properties of the medium is a classical problem in fluid mechanics. Here we study a special case of a two dimensional model of a suspension of undeformable particles in a liquid crystal solvent. In the dilute regime, we calculate asymptotic solutions of the perturbations of the velocity and director fields and derive an explicit formula for an effective shear viscosity of the liquid crystal medium. Such effective shear viscosity increases linearly with the area fraction of particles, similar to Einstein formula but with a different prefactor. We provide explicit asymptotic formulas for the dependence of this prefactor on the material parameters of the solvent. Finally, we identify a case of decreasing the effective viscosity by increasing the magnitude of the shear-flow alignment coefficient of the liquid crystal solvent.

## INTRODUCTION:

Liquid crystals are anisotropic fluids where particles align on average along a common direction, called the director field [1]. In nematic phases, the director field is equivalent to its opposite and therefore characterises an axis in space. Due to their unique properties, liquid crystals are of considerable technological importance and have attracted interest in a broad range of fields. Some recent examples are from ferroelectric liquid crystals to thin liquid-crystal films or from active nematics to biological physics [2–5].

The motion of colloidal particles can be influenced by properties of the liquid crystal solvent [6]. For example, the director field in the surrounding region of particles can be distorted due to the anchoring of the director field on the interface of particles. These distortions can, for instance, induce additional long-ranged interactions between particles [7–11]. Moreover, the motion of colloidal particles can be used to measure rheological properties of a solvent. A method that has been widely used in biological systems like actin or microtubule suspension, [12–17], bacteria suspensions [18–21] or tissues [22–25]. Here we are interested in understanding how the bulk properties of liquid crystals are modified by a suspension of particles.

Several mathematical approaches have been developed for coarse-graining microscopic models. Among them, homogenization theory offers mathematical tools to deal with well-separated scales both numerically and analytically. In the context of fluid mechanics, this line of research dates back to the seminal work by Einstein 1905 [26] who calculated the homogenized viscosity of a dilute suspension of particles in Stokesian fluids. Subsequently, this work was extended to non-dilute regimes, requiring the inclusion of pairwise particle interactions [27].

Since then, numerous extensions of the problem have been explored, including suspensions of particles with different shapes [28–30] or suspensions with particles that self-propel [31–34]. Other approaches such as the BBGKY hierarchy has been widely used to coarse-grain assemblies of self-propelled particles [35–41].

Various experimental and theoretical techniques have been developed to measure the effective viscosities of liquid crystals. The pioneering work of Miesowicz [42] first reported on the anisotropy of viscosity coefficients in pure liquid crystals under shear flow. In [43, 44], the authors computed the frictional drag on a sphere moving through a liquid crystal solvent with a uniform director alignment parallel to the flow direction at infinity. More recently, [45] introduced the shear horizontal wave method, which measures an effective viscosity averaged over the depth of the liquid crystal solvent, offering a novel experimental approach to obtain the bulk viscosity.

In this work, we first introduce a 2D mathematical model for describing a suspension of particles immersed in a liquid crystalline medium based on the Ericksen-Leslie model. Next, we perform a two-scale expansion to obtain the local problem for a single particle, which is known in homogenization theory as the cell problem. Then, we determine the effective viscosity in the dilute regime up to first order in area fraction. Finally, we discuss a case of decreasing the effective viscosity by changing parameters of the liquid crystal solvent.

## ERICKSEN-LESLIE MODEL FOR SUSPENSIONS

In this section, we present the theoretical framework for describing a suspension of solid particles in a liquid crystal. This suspension is modeled by an $\varepsilon$-periodic square array of undeformable identical discs $P^\varepsilon = \bigcup_n P_n^\varepsilon$
immersed in a liquid crystaline medium occupying a domain $\Omega^\epsilon = \Omega \setminus \mathcal{P}^\epsilon$ with $\Omega \subset \mathbb{R}^2$, see Fig. 1a. The parameter $\epsilon$ is the spacing between centers of particles (i.e. period of the square array) and the parameter $a^\epsilon$ is the radius of the discs. In the subsequent sections, we consider the dilute regime by taking the limit $a^\epsilon/\epsilon \to 0$. Note that because the domain $\Omega$ is divided into an $\epsilon$-periodic array, it can be partitioned into a periodic array of cells where each cell $Y_f$ contains a single particle $\mathcal{P}_n^\epsilon$, see Fig. 1b.

<figure><img src="image_2.png" /><figcaption>FIG. 1. Schematic of a periodic suspension of particles in a liquid crystal. a) $\epsilon$-periodic square array of undeformable identical discs (black) $\mathcal{P}^\epsilon = \bigcup_n \mathcal{P}_n^\epsilon$. The radii of the particles is $a^\epsilon$. The particles are immersed in a liquid crystalline medium occupying a domain $\Omega^\epsilon = \Omega \setminus \mathcal{P}^\epsilon$ with $\Omega \subset \mathbb{R}^2$ being the domain of the container. The vector $\hat{s} = \{\cos(\theta_s), \sin(\theta_s)\}$ is the anchoring on the boundary of the container $\partial\Omega$ with phase $\theta_s$. b) Schematic of the unit cell problem on domain $Y_f$. The dimensions of the periodic cell are rescaled such that the size of the domain $Y_f$ is 1 and the size of the particle is $a$. The vector $\hat{t} = \{\cos(\theta_t), \sin(\theta_t)\}$ is the anchoring of the director on the particle boundary $\partial\mathcal{P}^\epsilon$ with phase $\theta_t$. Dotted black lines represent the director field $\mathbf{n}^\epsilon$.</figcaption></figure>

While, in the literature, there are several descriptions for the liquid crystal solvent, in this work, we choose the Ericksen-Leslie model (1)-(4) because it is amenable to analytical studies. The liquid crystal is described by a velocity field $\mathbf{u}^\epsilon$ and a director field $\mathbf{n}^\epsilon$. Because the liquid crystal has nematic symmetry, the system of continuum equations governing the hydrodynamics of the liquid crystal are invariant under the transformation $\mathbf{n}^\epsilon \to -\mathbf{n}^\epsilon$ and it takes the form

$$ \nabla \cdot \boldsymbol{\sigma}^{\epsilon} = 0 \qquad x \in \Omega^{\epsilon} \quad (1) $$

$$ \nabla \cdot \mathbf{u}^{\epsilon} = 0 \qquad x \in \Omega^{\epsilon} \quad (2) $$

$$ \partial_t \mathbf{n}^{\epsilon} + (\mathbf{u}^{\epsilon} \cdot \nabla) \mathbf{n}^{\epsilon} + A(\mathbf{u}^{\epsilon}) \cdot \mathbf{n}^{\epsilon} = \Gamma \mathbf{h}^{\epsilon} - \nu D(\mathbf{u}^{\epsilon}) \cdot \mathbf{n}^{\epsilon} \quad x \in \Omega^{\epsilon} \quad (3) $$

$$ |\mathbf{n}^{\epsilon}| = 1 \qquad x \in \Omega^{\epsilon} \quad (4) $$

$$ |\mathbf{n}^{\epsilon}| = 1 \qquad x \in \Omega^{\epsilon} \quad (4) $$

where $\boldsymbol{\sigma}^\epsilon$ is the total stress tensor [1] of the form

$$ \boldsymbol{\sigma}^{\epsilon}(x) = 2\eta D(\mathbf{u}^{\epsilon}) - p^{\epsilon} + \frac{\nu}{2}(\mathbf{n}^{\epsilon}\mathbf{h}^{\epsilon} + \mathbf{h}^{\epsilon}\mathbf{n}^{\epsilon}) + \frac{1}{2}(\mathbf{n}^{\epsilon}\mathbf{h}^{\epsilon} - \mathbf{h}^{\epsilon}\mathbf{n}^{\epsilon}) - K(\nabla\mathbf{n}^{\epsilon})(\nabla\mathbf{n}^{\epsilon}) \quad (5) $$

In the limiting case where inertial effects are negligible and in the absence of external forces, the conservation of momentum reduces to the force balance equation (1). In Eq. (5), the first term on the right hand side (RHS) is the viscous stresses with the shear viscosity $\eta$, where $D(\mathbf{u}^\epsilon)$ is the symmetric part of the velocity gradient tensor $D(\mathbf{u}^\epsilon) = (\partial_\alpha \mathbf{u}_\beta^\epsilon + \partial_\beta \mathbf{u}_\alpha^\epsilon)/2$. The second term $p^\epsilon$ is the pressure field, which acts as a Lagrange multiplier to enforce the incompresibility constraint (2). The next three terms are stresses due to distortions of the director field which depend on the field $\mathbf{h}^\epsilon = -\delta\mathcal{F}_d/\delta\mathbf{n}^\epsilon$ often called the molecular field in the context of liquid crystals [1]. Here $\delta/\delta\mathbf{n}^\epsilon$ refers to the functional derivative of $\mathcal{F}_d$ which is the part of the total free-energy that depends on distortions of the director field. Here, $\mathcal{F}_d$ is given by the Frank free-energy in the one-constant approximation [1]

$$ \mathcal{F}_d = \int_{\Omega^{\epsilon}} \left( \frac{K}{2} (\nabla \mathbf{n}^{\epsilon})^2 + \lambda^{\epsilon} |\mathbf{n}^{\epsilon}|^2 \right) da \quad (6) $$

The elastic coefficient $K$ is the reduced Frank constant in two dimensions and $da$ is an infinitesimal element of area. The field $\lambda^\epsilon$ is a Lagrange multiplier that enforces the constraint (4), which means that the liquid crystal is deep into a nematic phase. Using Eq. (6), the expression of the molecular field becomes

$$ \mathbf{h}^{\epsilon} := -\frac{\delta\mathcal{F}}{\delta\mathbf{n}^{\epsilon}} = K\Delta\mathbf{n}^{\epsilon} - \lambda^{\epsilon}\mathbf{n}^{\epsilon}, \quad (7) $$

The dynamics of the director field are governed by Eq. (3). The terms on the left hand side (LHS) of Eq. (3) represent the co-rotational convective time derivative of $\mathbf{n}^\epsilon$, where $A(\mathbf{u}^\epsilon)$ is the anti-symmetric part of the velocity gradient tensor $A(\mathbf{u}^\epsilon) = (\partial_\alpha \mathbf{u}_\beta^\epsilon - \partial_\beta \mathbf{u}_\alpha^\epsilon)/2$. The first term on the RHS of Eq. (3) is the elastic torque associated with perturbation of the director field, where $\Gamma$ is the inverse of the rotational viscosity. The second term on the RHS in Eq. (3) couples the dynamics of $\mathbf{n}^\epsilon$ to the local shear flows with the parameter $\nu$ called the shear-flow alignment coefficient. This coefficient controls two regimes of the dynamics of the director field: The tumbling regime $|\nu| < 1$, where the director field rotates when subjected to shear flows, and the aligning regime $|\nu| > 1$, where the director field aligns when subjected to shear flows. It follows from the Onsager relations that the coefficient $\nu$ in Eq. (5) and in Eq. (3) are equal [1].

As shown in Appendix A, Eq. (3) can be further simplified by representing the director field as

$$ \mathbf{n}^{\epsilon} = \{\cos(\psi^{\epsilon}), \sin(\psi^{\epsilon})\} \quad (8) $$

where $\psi^\epsilon$ is the angle of the director field with respect to an arbitrary axis. Note that this expression automatically satisfies the constraint (4). Then, the continuum
equations (1)-(4) become:

$$ \nabla \cdot \boldsymbol{\sigma}^{\varepsilon} = 0 \qquad x \in \Omega^{\varepsilon} \quad (9) $$

$$ \nabla \cdot \mathbf{u}^{\varepsilon} = 0 \qquad x \in \Omega^{\varepsilon} \quad (10) $$

$$ \begin{aligned} \partial_t \psi^{\varepsilon} + (\mathbf{u}^{\varepsilon} \cdot \nabla)\psi^{\varepsilon} + \mathbf{n}_{\perp}^{\varepsilon} A(\mathbf{u}^{\varepsilon})\mathbf{n}^{\varepsilon} & = \Gamma h_{\perp}^{\varepsilon} - \nu \mathbf{n}_{\perp}^{\varepsilon} D(\mathbf{u}^{\varepsilon})\mathbf{n}^{\varepsilon} \quad x \in \Omega^{\varepsilon} \quad (11) \\ \Gamma h_{\parallel}^{\varepsilon} & = \nu \mathbf{n}^{\varepsilon} D(\mathbf{u}^{\varepsilon})\mathbf{n}^{\varepsilon} \quad x \in \Omega^{\varepsilon} \quad (12) \end{aligned} $$

where

$$ h_{\perp}^{\varepsilon} = \mathbf{n}_{\perp}^{\varepsilon} \cdot \mathbf{h}^{\varepsilon} = K \Delta \psi^{\varepsilon} \quad (13) $$

$$ h_{\parallel}^{\varepsilon} = \mathbf{n}^{\varepsilon} \cdot \mathbf{h}^{\varepsilon} = -K |\nabla \psi^{\varepsilon}|^2 - \lambda^{\varepsilon} \quad (14) $$

are the perpendicular and parallel components of $\mathbf{h}^{\varepsilon}$ with respect to $\mathbf{n}_{\perp}^{\varepsilon} = \{-\sin(\psi^{\varepsilon}), \cos(\psi^{\varepsilon})\}$ and $\mathbf{n}^{\varepsilon}$ respectively. Since $\lambda^{\varepsilon}$ is a Lagrange multiplier, then $h_{\parallel}^{\varepsilon}$ in Eq. (14) becomes a Lagrange multiplier that is fixed through Eq. (12).

Next, we introduce boundary conditions for the system (9)-(12). Here, there are two boundaries: the boundary of the particles $\partial \mathcal{P}^{\varepsilon} = \bigcup_n \partial \mathcal{P}_n^{\varepsilon}$ and the exterior boundary of the container $\partial \Omega$. The set of boundary conditions reads

$$ D(\mathbf{u}^{\varepsilon}) = 0 \qquad x \in \mathcal{P}^{\varepsilon} \quad (15) $$

$$ \mathbf{u}^{\varepsilon} = \mathbf{E} \cdot \mathbf{x} \qquad x \in \partial \Omega \quad (16) $$

$$ K(\hat{\mathbf{N}} \cdot \nabla)\psi^{\varepsilon} = W \sin(2(\theta_t(\mathbf{x}) - \psi^{\varepsilon})) \qquad x \in \partial \mathcal{P}^{\varepsilon} \quad (17) $$

$$ \psi^{\varepsilon} = \theta_s(\mathbf{x}) \qquad x \in \partial \Omega \quad (18) $$

The rigidity condition (15) extends the velocity $\mathbf{u}^{\varepsilon}$ to a constant inside $\mathcal{P}^{\varepsilon}$ and so that at the particle boundary $\partial \mathcal{P}^{\varepsilon}$ the velocity field is continuous. At the container boundary, we impose in Eq. (16) that the velocity is a linear shear flow given by $\mathbf{E} \cdot \mathbf{x}$, where $\mathbf{E}$ is a symmetric traceless matrix. In addition, the director field on the particle boundary satisfies the torque balance equation (17), where the term on the RHS of (17) corresponds to the torques generate by misalignment of the director field with respect to a preferred anchoring $\hat{\mathbf{t}} = \{\cos(\theta_t(\mathbf{x})), \sin(\theta_t(\mathbf{x}))\}$ with phase $\theta_t(\mathbf{x})$, Fig. 1b. The strength of the anchoring alignment is controlled by the parameter $W$ in (17). Vector $\hat{\mathbf{N}}$ is the outward normal to the particle boundary. Likewise, on the container boundary in (18), we impose that the director field aligns with a preferred anchoring $\hat{\mathbf{s}} = \{\cos(\theta_s(\mathbf{x})), \sin(\theta_s(\mathbf{x}))\}$ with a fixed given phase $\theta_s(\mathbf{x})$, Fig. 1a.

Note that the model (9)-(18) depends on the following set of parameters: the size of a cell $\varepsilon$, the radii of discs $a^{\varepsilon}$, and the size of the container $|\Omega|$ which does not depend on $\varepsilon$, the shear viscosity $\eta$ and the rotational viscosity $1/\Gamma$, the elastic coefficients $K$ and $W$, two given phases $\theta_s(\mathbf{x})$ on the boundary of the container and $\theta_t(\mathbf{x})$ on the boundary of the particles, the shear flow rate $\mathbf{E}$, which has two independent components, and the dimensionless shear-flow alignment coefficient $\nu$. In the following, we consider the special case where the phases are

equal and constant $\theta_s = \theta_t = \text{constant}$ and without loss of generality this constant phase is set to zero. Moreover, we choose an extensional shear flow at infinity for shear strength $\gamma > 0$

$$ \mathbf{E} = \gamma \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \quad (19) $$

Our focus is on steady state solutions and therefore we consider time-independent physical fields and disregard the time derivative in (11).

Note that in the absence of particles, the liquid crystal in the configuration $\psi^{\varepsilon} = 0$ is linearly stable to long-wavelength perturbations of the phase when $\nu\gamma < 0$, [1]. In this work, we choose to consider the range of $\gamma > 0$ and therefore the stable range is $\nu < 0$.

## TWO-SCALE EXPANSION

In this section, we perform a two-scale expansion of the solutions $(\mathbf{u}^{\varepsilon}, \psi^{\varepsilon})$ to the Eqs. (9)-(18) (see also (8)). We consider the limit $\varepsilon \ll 1$ and derive the so-called cell problem which describes the local flows around a particle. We solve the cell problems under two conditions on the range of the parameters $\nu$ and $a^{\varepsilon}$.

I. Two limiting regimes of the shear-flow alignment coefficient, where cell problem can be solved asymptotically: $|\nu| \ll 1$ and $|\nu| \gg 1$.

II. The dilute regime where $a^{\varepsilon}/\varepsilon \to 0$ and interparticle interactions can be neglected.

Two-scale expansion is an asymptotic technique that is useful in constructing approximate solutions for problems with two well-separated scales. This technique has been used in a wide variety of contexts e.g., [46–48]. For example, it has been used to study suspensions in the Stokes fluids [49] as well as nonhomogeneous materials, including elastic materials and conductive materials [50, 51]. In short, we consider the physical fields, $\mathbf{u}^{\varepsilon}$ and $\mathbf{n}^{\varepsilon}$ (equiv. $\psi^{\varepsilon}$), as vectorial functions of slow varying spatial coordinate $\mathbf{x}$ and “fast” variable $\mathbf{y} = \mathbf{x}/\varepsilon$. This expansion describes two well-separated length scales in the system and the key mathematical idea is to treat $\mathbf{x}$ and $\mathbf{y}$ as independent variables. The slow variable $\mathbf{x}$ represents variations of the physical fields on the scale of the domain $\Omega$ which we normalize to be of order $O(1)$ as $\varepsilon \to 0$. On the other hand, the fast variable $\mathbf{y}$ represents variations of the physical fields on the scale of the interparticle distance of order $O(\varepsilon)$. Then, the gradient operator writes as:

$$ \nabla \mathbf{f}(\mathbf{x}, \mathbf{y}) = \left( \nabla^x + \frac{1}{\varepsilon} \nabla^y \right) \mathbf{f}(\mathbf{x}, \mathbf{y}) \quad (20) $$

where the superscripts in $\nabla$ denote derivatives with respect to either the slow variable $\mathbf{x}$ or the fast variable $\mathbf{y}$.
We expand the unknown physical fields in a power series in $\epsilon$ for $\mathbf{u}^{\epsilon}$, $\psi^{\epsilon}$, $p^{\epsilon}$, $h_{\perp}^{\epsilon}$ and $h_{\parallel}^{\epsilon}$:

$$ \left\{ \begin{aligned} \mathbf{u}^{\epsilon}(\mathbf{x}, \mathbf{y}) &= \sum_{n=0}^{\infty} \mathbf{u}^{(n)}(\mathbf{x}, \mathbf{y}) \epsilon^n \\ \psi^{\epsilon}(\mathbf{x}, \mathbf{y}) &= \sum_{n=0}^{\infty} \psi^{(n)}(\mathbf{x}, \mathbf{y}) \epsilon^n \\ p^{\epsilon}(\mathbf{x}, \mathbf{y}) &= \sum_{n=-1}^{\infty} p^{(n)}(\mathbf{x}, \mathbf{y}) \epsilon^n \\ h_{\perp}^{\epsilon}(\mathbf{x}, \mathbf{y}) &= \sum_{n=-2}^{\infty} h_{\perp}^{(n)}(\mathbf{x}, \mathbf{y}) \epsilon^n \\ h_{\parallel}^{\epsilon}(\mathbf{x}, \mathbf{y}) &= \sum_{n=-2}^{\infty} h_{\parallel}^{(n)}(\mathbf{x}, \mathbf{y}) \epsilon^n \end{aligned} \right. \qquad (21) $$

The superscripts in the coefficient fields denote the order of the term in the two-scale expansion (21). Because the domain $\Omega$ is partitioned into an $\epsilon$-periodic array of cells $Y_f^{\epsilon}$, we consider all terms in the RHS of Eq. (21) to have the same $\epsilon$-periodicity in $\mathbf{y}$.

Substituting Eqs. (21) into the system (9)-(12) and (15)-(18) (see Appendix B to D for calculations), we obtain at the lowest order $\epsilon^{-2}$

$$ \psi^{(0)}(\mathbf{x}) = \theta_s \qquad (22) $$

$$ \mathbf{u}^{(0)}(\mathbf{x}) = \mathbf{E} \cdot \mathbf{x} \qquad (23) $$

$$ p^{(-1)} = h_{\perp}^{(-2)} = h_{\parallel}^{(-2)} = 0 \qquad (24) $$

The lowest order terms given by Eqs. (22)-(24) do not depend on the fast variable $\mathbf{y}$ and therefore represent the macroscopic (or coarse-grained) physical fields.

Due to the presence of particles, local perturbations can be generated on the physical fields, which are accounted for by the next order terms in the two-scale expansions Eqs. (21). To order $\epsilon^{-1}$, the perturbations of the director phase and of the two components of the molecular field vanish, $\psi^{(1)} = 0$ and $h_{\perp}^{(-1)} = h_{\parallel}^{(-1)} = 0$ (see Appendix B). Moreover, the perturbations of the velocity and pressure field solve the following "cell problem" (see Appendix C and D)

$$ \eta_1 \partial_{11}^y u_1^{(1)} + \eta_2 \partial_{22}^y u_1^{(1)} - \partial_1^y \bar{p}^{(0)} = 0 \qquad y \in Y_f \qquad (25) $$

$$ \eta_1 \partial_{11}^y u_2^{(1)} + \eta_2 \partial_{22}^y u_2^{(1)} - \partial_2^y \bar{p}^{(0)} = 0 \qquad y \in Y_f \qquad (26) $$

$$ \partial_1^y u_1^{(1)} + \partial_2^y u_2^{(1)} = 0 \qquad y \in Y_f \qquad (27) $$

$$ \mathbf{u}^{(1)} \text{ periodic} \qquad y \in \partial Y_f \qquad (28) $$

$$ \mathbf{u}^{(1)} = -\mathbf{E} \cdot \mathbf{y} \qquad y \in \partial \mathcal{P} \qquad (29) $$

where $\eta_1 = \eta + (\nu - 1)^2/4\Gamma$ and $\eta_2 = \eta + (\nu + 1)^2/4\Gamma$ are two effective viscosities. The subscripts in the derivatives denote the Cartesian coordinates in two dimensions. The pressure field $p^{(0)}$ is expressed in terms of $\bar{p}^{(0)}$, which satisfies the Laplace equation (see Appendix D). Because the ansatz (21) is $\epsilon$-periodic in the fast variable $\mathbf{y}$, it is sufficient to find the solution in a single cell $Y_f$ with periodic boundary conditions on $\partial Y_f$ (28) and the boundary conditions (29) on $\partial \mathcal{P}$. This is the so-called cell problem which depends on the radius of the particles and the material parameters of the liquid crystal solvent, such as

the shear-flow alignment coefficient $\nu$. We note that since the cell problem only depends on $\mathbf{y}$ the domain $Y_f$ is of size 1 and the particle radius becomes $a$, see Fig. 2a

<figure><img src="image_3.png" /><figcaption>FIG. 2. Schematic of the unit cell domain $Y_f$ and the rescaled domain in the dilute limit. a) In the unit cell domain $Y_f$, the particle is of size $a$ in a container of size 1. In this work, we choose the scaling of $a$ with $\epsilon$ in Eq. (32). b) Rescaled domain for the choice of scaling of $a$ with $\epsilon$ in Eq. (32). In this domain, the particle is of size 1 in a container of size $\epsilon^{-1}$. The dotted curves represent the director field and the black circles represent the particles.</figcaption></figure>

The cell problem (25)-(27) describe an anisotropic fluid with two viscosities $\eta_1$ and $\eta_2$. Because the fluid is incompressible, the velocity field can be expressed in terms of a streamfunction $\Psi^{(1)}$, such that $u_1^{(1)} = -\partial_2^y \Psi^{(1)}$ and $u_2^{(1)} = \partial_1^y \Psi^{(1)}$. This streamfunction is a solution of a generalised (anisotropic) biharmonic equation (see Appendix D) with the general form

$$ \Psi^{(1)}(\mathbf{y}) = \mathrm{Re} \left[ f(y_1 + iy_2) + g \left( \frac{y_1}{\sqrt{\eta_1}} + i \frac{y_2}{\sqrt{\eta_2}} \right) \right] \qquad (30) $$

where $f$ and $g$ are analytic functions and $y_1$ and $y_2$ are the Cartesian coordinates of the fast variable $\mathbf{y}$. By incompressibility, the pressure field $\bar{p}^{(0)}$ satisfies the $\Delta \bar{p}^{(0)} = 0$. In order to solve the cell problem (25)-(27) analytically, we focus on the limiting regime I introduced at the beginning of this section controlled by the parameter $\nu$ (see Appendix E for calculations).

In the limit of small shear-flow alignment $|\nu| \ll 1$, we solve the cell problem by expanding the perturbations of the streamfunction and pressure field as a power series in $\nu$.

$$ \left\{ \begin{aligned} \Psi^{(1)}(\mathbf{y}) &= \sum_{k=0} \nu^k \Psi_s^{(1,k)}(\mathbf{y}) \\ \bar{p}^{(0)}(\mathbf{y}) &= \sum_{k=0} \nu^k \bar{p}_s^{(0,k)}(\mathbf{y}) \end{aligned} \right. \qquad (31) $$

This allows us to transform (25)-(27) into a set of PDEs for the functions $\Psi^{(1,k)}(\mathbf{y})$ and $\bar{p}^{(0,k)}(\mathbf{y})$ which describe an isotropic fluid solvent with shear viscosity $\eta + 1/4\Gamma$ and forcing terms determined by the previous order in $\nu$ (see Appendix E).

Next, we consider the dilute limit of the cell problem where the particle radius $a \to 0$. Since the PDE for $\mathbf{u}^{(1)}$ reduces to an isotropic Stokes equation in the limit of $|\nu| \ll 1$, we choose the scaling [49]

$$ a^{\epsilon} \propto \epsilon^2, \quad a \propto \epsilon \qquad (32) $$
Later we will show that this choice of scaling for the radius $a$ leads to corrections to the rate of total dissipation for the suspension.

The cell problem on domain $Y_f$ can be rescaled so the radius of the particle is of size 1 and the domain is of size $\varepsilon^{-1}$, Fig. 2b. In the dilute limit ($a \to 0$) this rescaled domain approaches $\mathbb{R}^2 \setminus P$ and we replace the periodic boundary conditions on $Y_f$ (28) with decay condition at infinity $\mathbf{u}^{(1)} \to 0$ at $\infty$. Then, we solve for $\Psi^{(1)}$ and $\bar{p}^{(0)}$ on the infinite domain. The solution in this infinite domain is an approximation to leading order of the cell problem on a domain of size 1 with a particle of size $a$, Fig. 2a. Up to order $k=2$ in $\nu$, the particular solution in the unit cell takes the form (see Appendix E)

$$ 
\left\{ 
\begin{aligned}
\Psi_s^{(1,0)} &= \gamma a^2 \left(1 - \frac{a^2}{2r^2}\right) \sin(2\theta) \\
\Psi_s^{(1,1)} &= \frac{\gamma a^2}{2(1+4\Gamma\eta)} \left(1 - \frac{a^2}{r^2}\right)^2 \sin(4\theta) \\
\Psi_s^{(1,2)} &= \frac{\gamma a^2}{6(1+4\Gamma\eta)^2} \left(2 - \frac{5a^2}{r^2}\right) \left(1 - \frac{a^2}{r^2}\right)^2 \sin(6\theta) \\
\bar{p}_s^{(0,0)} &= -\frac{\gamma a^2 (1+4\Gamma\eta)}{\Gamma r^2} \cos(2\theta) \\
\bar{p}_s^{(0,1)} &= 0 \\
\bar{p}_s^{(0,2)} &= -\frac{4\gamma a^2 \eta}{(1+4\Gamma\eta)r^2} \cos(2\theta)
\end{aligned}
\right. 
\quad (33) 
$$

where the Cartesian coordinates $y_1 = r \cos \theta$ and $y_2 = r \sin \theta$ are expressed in the polar coordinates $r$ and $\theta$.

In the limit of large shear-flow alignment $|\nu| \gg 1$, we solve the cell problem in a similar way and expand the perturbations of the streamfunction and pressure field as a power series in $1/\nu$.

$$ 
\left\{ 
\begin{aligned}
\Psi^{(1)}(\mathbf{y}) &= \sum_{k=0} \nu^{-k} \Psi_l^{(1,k)}(\mathbf{y}) \\
\bar{p}^{(0)}(\mathbf{y}) &= \nu^2 \sum_{k=0} \nu^{-k} \bar{p}_l^{(0,k)}(\mathbf{y})
\end{aligned}
\right. 
\quad (34) 
$$

By a similar procedure, up to order $k=2$ in $1/\nu$, the particular solution in the unit cell takes the form (see Appendix E)

$$ 
\left\{ 
\begin{aligned}
\Psi_l^{(1,0)} &= \gamma a^2 \left(1 - \frac{a^2}{2r^2}\right) \sin(2\theta) \\
\Psi_l^{(1,1)} &= \frac{\gamma a^2}{2} \left(1 - \frac{a^2}{r^2}\right)^2 \sin(4\theta) \\
\Psi_l^{(1,2)} &= \frac{\gamma a^2}{6} \left(2 - \frac{5a^2}{r^2}\right) \left(1 - \frac{a^2}{r^2}\right)^2 \sin(6\theta) \\
\bar{p}_l^{(0,0)} &= -\frac{\gamma a^2}{\Gamma r^2} \cos(2\theta) \\
\bar{p}_l^{(0,1)} &= 0 \\
\bar{p}_l^{(0,2)} &= -\frac{4\gamma a^2 \eta}{r^2} \cos(2\theta)
\end{aligned}
\right. 
\quad (35) 
$$

where the Cartesian coordinates $y_1 = r \cos \theta$ and $y_2 = r \sin \theta$ are expressed in the polar coordinates $r$ and $\theta$.

In both limiting regimes of $\nu$, the net force and net torque on the particle boundary associated with either perturbation (33) or (35) vanish, and therefore the particle's linear and angular velocity vanish (see Appendix F for calculations).

## DISSIPATION RATE AND EFFECTIVE VISCOSITY

In this section, we will determine an effective viscosity of a suspension of colloidal particles in a liquid crystal solvent. Using the results from the two-scale expansion in the previous section, we will determine the rate of total dissipation for two systems: a dilute suspension of particles in a liquid crystal solvent and an homogenized liquid crystal. The comparison between the two rates of total dissipation allows us to identify an effective viscosity for the suspension and its dependence on the particle density and material parameters of the liquid crystal solvent.

Finding an equivalent homogenized medium without particles that produces the same rate of total dissipation in a macroscopic volume as a suspension in a solvent is a classical problem in fluid dynamics, [26, 27]. In the case of a suspension of rigid particles in a liquid crystal solvent, the rate of total dissipation is equal to the time variation of the total free-energy which takes the form

$$ \frac{dF}{dt} = - \int_{\Omega^{\varepsilon}} 2\eta D(\mathbf{u}^{\varepsilon}) : D(\mathbf{u}^{\varepsilon}) + \Gamma \mathbf{h}^{\varepsilon} \cdot \mathbf{h}^{\varepsilon} da \quad (36) $$

where $\mathbf{u}^{\varepsilon}$ is defined in (21), $\mathbf{h}^{\varepsilon}$ is defined in (7), and the operator : denotes a contraction of the two indices of the symmetric velocity gradient tensor $D(\mathbf{u}^{\varepsilon})$. In this expression, the first term is the rate of dissipation due to viscous stresses and the second term is the rate of dissipation due to dynamics of the director field, which has no counterpart in Newtonian fluids.

We first focus on the dissipation due to viscous stresses in Eq. (36). We express the area of the container $|\Omega|$ minus the area of the particles of radius $a^{\varepsilon}$

$$ |\Omega^{\varepsilon}| = |\Omega| - N_p \pi (a^{\varepsilon})^2 \quad (37) $$

where $N_p = C\varepsilon^{-2}$ is the number of particles, see Fig. 1a and $a^{\varepsilon}$ is given by Eq. (32). Substituting the two-scale expansion (21) into the first term in the RHS of Eq. (36) to obtain

$$ \int_{\Omega^{\varepsilon}} 2\eta D(\mathbf{u}^{\varepsilon}) : D(\mathbf{u}^{\varepsilon}) d\mathbf{x} = 4\eta\gamma^2 (|\Omega| - N_p\pi(a^{\varepsilon})^2) + \int_{\Omega^{\varepsilon}} 2\eta D^y(\mathbf{u}^{(1)}) : D^y(\mathbf{u}^{(1)}) d\mathbf{x} + \mathcal{O}(\varepsilon^3) \quad (38) $$

where we used that $\int_{\Omega^{\varepsilon}} D^y(\mathbf{u}^{(1)}) da = 0$ and $\gamma$ is the shear strength. The first contribution in the RHS in Eq. (38) is the rate of dissipation due to the macroscopic shear flow in the solvent. The second term in the RHS in Eq. (38) is the rate of dissipation due to the flow perturbations in the surrounding region of particles. As written, this term has contributions of order $\varepsilon^n$ for $n \ge 2$ which comes from the dependence of $\Psi_s^{(1,k)}$ or $\Psi_l^{(1,k)}$ on the radius of the particle $a$ in Eqs. (33) and (35) and the dependence of the domain $\Omega^{\varepsilon}$ on $a^{\varepsilon}$. In the following, we ignore the additional contributions at order $\varepsilon^n$ for $n > 2$.
Because Eqs. (25)-(29) are linear in the perturbation flows $\mathbf{u}^{(1)}$ and $\mathbf{u}^{(1)}$ is $\varepsilon$-periodic in $\mathbf{y}$, the second term in the RHS in Eq. (38) is expressed as the sum of the rate of dissipation by the flow perturbations of a single particle $\mathcal{P}$ in the unit cell, Fig. 2a.

$$ \int_{\Omega^{\varepsilon}} 2\eta D^{y}(\mathbf{u}^{(1)}): D^{y}(\mathbf{u}^{(1)})d\mathbf{x} = N_{p}\varepsilon^{2} \int_{Y_{f}} 2\eta D^{y}(\mathbf{u}^{(1)}): D^{y}(\mathbf{u}^{(1)})d\mathbf{y} \quad (39) $$

We focus on the dilute regime $\varepsilon \to 0$, (32), where inter-particle interactions can be neglected, and approximate the integral on the RHS by using the solutions of the perturbations in the two limiting cases of $\nu$, (33) and (35), where the integration domain is replaced by $\mathbb{R}^2/\mathcal{P}$. The results obtained take the form

$$ -\int_{\Omega^{\varepsilon}} 2\eta D(\mathbf{u}^{\varepsilon}): D(\mathbf{u}^{\varepsilon})d\mathbf{a} = -4\eta\gamma^2\pi R^2 * \left\{ \begin{array}{ll} (1 + \phi\frac{2\nu^2}{(1+4\alpha)^2}) + \mathcal{O}(\phi^2) + \mathcal{O}(|\nu|^3) & \text{if } |\nu| \ll 1 \\ 1 + \mathcal{O}(\phi^2) + \mathcal{O}(|\nu|^{-1}) & \text{if } |\nu| \gg 1 \end{array} \right. \quad (40) $$

where $\phi = N_p(a^\varepsilon)^2/R^2 = \varepsilon^2/R^2$ is the area fraction of particles in the container $\Omega$, and $\alpha = \Gamma\eta$ is a dimensionless parameter that arises from the ratio between the two viscosities $\eta$ and $1/\Gamma$.

Recall that the limit in $\nu$ is performed before the limit in $a$. In a similar way, the rate of dissipation due to dynamics of the director field in Eq. (36) takes the form in the limit $|\nu| \ll 1$,

$$ -\int_{\Omega^{\varepsilon}} \Gamma \mathbf{h}^{\varepsilon} \cdot \mathbf{h}^{\varepsilon} d\mathbf{a} = -\frac{\gamma^2 \pi R^2}{\Gamma} * \left( \nu^2 + 2\phi \left( 1 - \frac{(1 + 8\alpha)\nu^2}{(1 + 4\alpha)^2} \right) \right) + \mathcal{O}(\phi^2) + \mathcal{O}(|\nu|^3) \quad (41) $$

and in the limit $|\nu| \gg 1$,

$$ -\int_{\Omega^{\varepsilon}} \Gamma \mathbf{h}^{\varepsilon} \cdot \mathbf{h}^{\varepsilon} d\mathbf{a} = -\frac{\gamma^2 \pi R^2}{\Gamma} (\nu^2 + 2\phi) + \mathcal{O}(\phi^2) + \mathcal{O}(|\nu|^{-1}) \quad (42) $$

Adding these two expressions yields the rate of total dissipation in the suspension. In the regime of small $|\nu| \ll 1$, dissipation due to viscous stresses (40) is modified by parameters that control the dynamics of the director field such as $\nu$ and $\Gamma$. Conversely, dissipation due to the dynamics of the director field (42) is modified by parameters that control the viscous stresses, such as the shear viscosity $\eta$. This shows that the presence of particles modifies the rate of total dissipation of the suspension by coupling the two sources of dissipation in the liquid crystal solvent.

The presence of a particle in a solvent influences the rate of total dissipation in suspensions in several ways.

On the one hand, particles create perturbations to the macroscopic fields in a region surrounding the particles, increasing the density of dissipation in the solvent. On the other hand, particles occupy space in the container which reduces the area of the solvent which reduces the total dissipation per unit time (see the first term on the RHS of Eq. (38)). The balance between these opposing contributions determines the first order correction in $\phi$ of the various sources of dissipation (40) and (42).

Next, we determine the rate of total dissipation of an homogenized liquid crystal with effective material parameters. For clarity we denote with the superscript $H$, the physical fields and material parameters of the homogenized medium. A priori, all material parameters of the homogenized medium can differ from the material parameters of the solvent. In this work we focus on the computation of the effective viscosity $\eta^H$ via the solution of the cell problem (25)-(29), but we note that the effective rotational viscosity $1/\Gamma^H$ and shear-flow alignment coefficient $\nu^H$ could be computed in a similar manner. For simplicity, we set $\Gamma^H = \Gamma$ and $\nu^H = \nu$. The homogenized medium is subjected to the same extensional shear flow and preferred anchoring at the boundaries of the container as the suspension (16) and (18). The expression of the physical fields of the homogenized medium are determined in Appendix G.

The time variation of the total free-energy in the homogenized liquid crystal is

$$ -\frac{d\mathcal{F}^H}{dt} = \int_{\Omega} 2\eta^H D(\mathbf{u}^H) : D(\mathbf{u}^H) + \Gamma \mathbf{h}^H \cdot \mathbf{h}^H d\mathbf{a} \quad (43) $$

Substituting the expression of the physical fields of the homogenized medium (see Appendix G), the rate of total dissipation up to $\mathcal{O}(\varepsilon^3)$ in (43) takes the form in the two limits of $\nu$

$$ \frac{d\mathcal{F}^H}{dt} = -C_R * \left\{ \begin{array}{ll} (1 - 2\phi) + \mathcal{O}(\phi^2) + \mathcal{O}(|\nu|^3) & \text{if } |\nu| \ll 1; \\ (1 - 2\phi) + \mathcal{O}(\phi^2) + \mathcal{O}(|\nu|^{-1}) & \text{if } |\nu| \gg 1. \end{array} \right. \quad (44) $$

where the constant $C_R = (4\eta^H + \frac{\nu^2}{\Gamma})\gamma^2\pi R^2$.

Balancing the sum of Eqs. (40) and (42) to Eq. (44) allows us to define an effective viscosity up to $\mathcal{O}(\varepsilon^3)$ for the homogenized medium in the limit $|\nu| \ll 1$:

$$ \frac{\eta^H}{\eta} = 1 + \phi \left( \frac{1 + 4\alpha}{2\alpha} - \frac{\nu^2(1 + 8\alpha)}{8\alpha^2(1 + 4\alpha)} \right) + \mathcal{O}(\phi^2) + \mathcal{O}(|\nu|^3) \quad (45) $$

and in the limit $|\nu| \gg 1$:

$$ \frac{\eta^H}{\eta} = 1 + 2\phi \left( 1 + \frac{1}{\nu^2} \right) + \mathcal{O}(\phi^2) + \mathcal{O}(|\nu|^{-1}) \quad (46) $$
Note that in the limit $\Gamma \to \infty$, we recover the well-known viscosity of a two-dimensional suspension in a Stokesian fluid $\eta^H/\eta = 1+2\phi+O(\phi^2)$. The effective viscosity (45) and (46) depends linearly on the area fraction of particle $\phi$ and the dimensionless parameters $\alpha$ and $\nu$.

For sufficiently small values of $|\nu| \ll 1$, the effective viscosity (45) depends on $\alpha$ even in the case of vanishing $\nu$. At a fixed area fraction, the effective viscosity can be reduced by modifying properties of the liquid crystal solvent, such as increasing the parameter $|\nu|$ or the dimensionless parameters $1/\alpha$. For sufficiently large values $|\nu| \gg 1$, we recover the Einstein formula for a Stokesian fluid in Eq. (46), $\eta^H/\eta = 1+2\phi+O(\phi^2)$.

Previous studies of suspensions of active swimmers in a Stokesian fluid have shown both experimentally and theoretically that the homogenized viscosity (i.e., effective viscosity) of the suspension can be decreased by increasing the area fraction $\phi$ [31–34, 52–55]. This effect is captured by a negative term due to activity in the classical Einstein formula. Observe that there is a negative term in (45) (second order in $\nu$ and first order in $\phi$). This term can cause a decrease of the effective viscosity by either increasing the parameter $|\nu|$ or decreasing the parameter $\alpha$ at a fixed area fraction $\phi$. However, it cannot lead to a decrease of the effective viscosity ($\eta^H < \eta$) by increasing the area fraction $\phi$ in the parameter regimes considered above ($|\nu| \ll 1$ and $|\nu| \gg 1$). To illustrate that in our case $\eta^H > 0$, we computed the density of energy dissipation per unit time in the solvent (i.e., the integrand in Eq. (36)), see Fig 3, and found that it is overall negative in the two regimes of $|\nu|$. As the area fraction increases the suspension dissipates more energy per unit time, leading to an increase in the effective viscosity ($\eta^H > \eta$), as shown in Eqs. (45) and (46). An interesting open question is whether there exists an intermediate parameter regime in $0 < |\nu| < \infty$ where this effective viscosity is negative.

DISCUSSION

In this work, we obtained an effective viscosity for a dilute suspension of solid particles in a liquid crystal solvent. A key aspect of our analysis is the consideration of two limits: $|\nu| \ll 1$ and $|\nu| \gg 1$ (shear-flow alignment) and $a^\epsilon/\epsilon \to 0$ (diluteness parameter). The order in which these limits are taken can be critical, as it could affect the resulting behavior of the macroscopic flow and the bulk properties of the homogenized liquid crystal. While $|\nu| \ll 1$ and $|\nu| \gg 1$ are regular perturbations of a material parameter of the liquid crystal solvent, the diluteness parameter $a^\epsilon/\epsilon$ is a geometric parameter which affects the solution at both microscopic and macroscopic scales and thus may not commute with the former limit in $\nu$. By performing the limit in $\nu$ first, we obtain an isotropic Stokes Eqs. (111)-(113) and (117)-(119) for which the

<figure><img src="image_4.png" /><figcaption>FIG. 3. Flow perturbation $\mathbf{u}^{(1)}$ in the limit $|\nu| \ll 1$ (tumbling phase) (a) and in the limit $|\nu| \gg 1$ (aligning phase) (b). For a background shear flow that extends in the horizontal axis and contracts in the vertical axis. Colormap of the density of energy dissipation per unit time (integrand in (36)) in the limit $|\nu| \ll 1$ (c) and in the limit $|\nu| \gg 1$ (d). Note that locally, the density of energy dissipation remains negative (shown in blue). In panels a-d, the parameter set is: $\gamma = 2$, $1/\Gamma = 10$, $\eta = 1$ and $a = 1$. In panels a and c, we used $\nu = 0.01$ and in panels b and d, we used $\nu = 10$. In panels c and d, the dashed circle represents the boundary of a particle with radius $a$.</figcaption></figure>

scaling of $a^\epsilon/\epsilon$ for the dilute regime are known [49]. Understanding the relationship between $\nu$ and $a^\epsilon/\epsilon$ deserves further attention in the future.

We would like to acknowledge that we learned of formula (30) from J. King (University of Nottingham) and we are grateful to him. The work of L. B. was partially supported by NSF grants DMS-2404546 and PHY-2140010. S. D. was partially supported by NSF grant PHY-2140010.

<footnote>* carles.blanch-mercader@curie.fr</footnote>

<footnote>† lvb2@psu.edu</footnote>

<footnote>[1] P. De Gennes and J. Prost, The physics of liquid crystals, 1993, Oxford University Press, New York, Olbrich E., Marinov O., Davidov D., Phys. Rev. E 2713, 48 (1993).</footnote>

<footnote>[2] T. Lopez-Leon and A. Fernandez-Nieves, Drops and shells of liquid crystal, Colloid and Polymer Science 289, 345 (2011).</footnote>

<footnote>[3] M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Hydrodynamics of soft active matter, Reviews of Modern Physics 85, 1143 (2013).</footnote>

<footnote>[4] T. B. Saw, W. Xi, B. Ladoux, and C. T. Lim, Biological tissues as active nematic liquid crystals, Advanced materials 30, 1802579 (2018).</footnote>
[5] C. Blanc, G. Durey, R. D. Kamien, T. Lopez-Leon, M. O. Lavrentovich, and L. Tran, Helfrich-hurault elastic instabilities driven by geometrical frustration, Reviews of Modern Physics **95**, 015004 (2023).

[6] H. Stark, Physics of colloidal dispersions in nematic liquid crystals, Physics Reports **351**, 387 (2001).

[7] F. Brochard and P. De Gennes, Theory of magnetic suspensions in liquid crystals, Journal de Physique **31**, 691 (1970).

[8] S. Ramaswamy, R. Nityananda, V. Raghunathan, and J. Prost, Power-law forces between particles in a nematic, Molecular Crystals and Liquid Crystals Science and Technology. Section A. Molecular Crystals and Liquid Crystals **288**, 175 (1996).

[9] P. Poulin, H. Stark, T. Lubensky, and D. Weitz, Novel colloidal interactions in anisotropic fluids, Science **275**, 1770 (1997).

[10] R. Ruhwandl and E. Terentjev, Long-range forces and aggregation of colloid particles in a nematic liquid crystal, Physical Review E **55**, 2958 (1997).

[11] T. Lubensky, D. Pettey, N. Currier, and H. Stark, Topological defects and interactions in nematic emulsions, Physical Review E **57**, 610 (1998).

[12] L. Hou, F. Lanni, and K. Luby-Phelps, Tracer diffusion in f-actin and ficoll mixtures. toward a model for cytoplasm, Biophysical journal **58**, 31 (1990).

[13] M. Dichtl and E. Sackmann, Colloidal probe study of short time local and long time reptational motion of semiflexible macromolecules in entangled networks, New Journal of Physics **1**, 18 (1999).

[14] T. G. Mason, T. Gisler, K. Kroy, E. Frey, and D. Weitz, Rheology of f-actin solutions determined from thermally driven tracer motion, Journal of Rheology **44**, 917 (2000).

[15] Y.-C. Lin, G. H. Koenderink, F. C. MacKintosh, and D. A. Weitz, Viscoelastic properties of microtubule networks, Macromolecules **40**, 7714 (2007).

[16] R. Zhang, N. Kumar, J. L. Ross, M. L. Gardel, and J. J. De Pablo, Interplay of structure, elasticity, and dynamics in actin-based nematic materials, Proceedings of the National Academy of Sciences **115**, E124 (2018).

[17] I. Vélez-Cerón, P. Guillamat, F. Sagués, and J. Ignés-Mullol, Probing active nematics with in situ microfabricated elastic inclusions, Proceedings of the National Academy of Sciences **121**, e2312494121 (2024).

[18] C. Valeriani, M. Li, J. Novosel, J. Arlt, and D. Maren-duzzo, Colloids in a bacterial bath: simulations and experiments, Soft Matter **7**, 5228 (2011).

[19] Y. Peng, L. Lai, Y.-S. Tai, K. Zhang, X. Xu, and X. Cheng, Diffusion of ellipsoids in bacterial suspensions, Physical review letters **116**, 068303 (2016).

[20] A. Lagarde, N. Dagès, T. Nemoto, V. Démery, D. Bar-tolo, and T. Gibaud, Colloidal transport in bacteria suspensions: from bacteria collision to anomalous and enhanced diffusion, Soft Matter **16**, 7503 (2020).

[21] N. Figueroa-Morales, M. M. Genkin, A. Sokolov, and I. S. Aranson, Non-symmetric pinning of topological defects in living liquid crystals, Communications Physics **5**, 301 (2022).

[22] O. Campàs, T. Mammoto, S. Hasso, R. A. Sperling, D. O'connell, A. G. Bischof, R. Maas, D. A. Weitz, L. Mahadevan, and D. E. Ingber, Quantifying cell-generated mechanical forces within living embryonic tissues, Nature methods **11**, 183 (2014).

[23] P. Guillamat, C. Blanch-Mercader, G. Pernollet,

K. Kruse, and A. Roux, Integer topological defects organize stresses driving tissue morphogenesis, Nature materials **21**, 588 (2022).

[24] J. Barbazan, C. Pérez-González, M. Gómez-González, M. Dedenon, S. Richon, E. Latorre, M. Serra, P. Mariani, S. Descroix, P. Sens, et al., Cancer-associated fibroblasts actively compress cancer cells and modulate mechanotransduction, Nature Communications **14**, 6966 (2023).

[25] A. Vian, M. Pochitaloff, S.-T. Yen, S. Kim, J. Pollock, Y. Liu, E. M. Sletten, and O. Campàs, In situ quantification of osmotic pressure within living embryonic tissues, Nature Communications **14**, 7023 (2023).

[26] A. Einstein, A new determination of molecular dimensions, Annln., Phys. **19**, 289 (1906).

[27] G. Batchelor and J. Green, The determination of the bulk stress in a suspension of spherical particles to order c2, Journal of Fluid Mechanics **56**, 401 (1972).

[28] M. Duerinckx and A. Gloria, Corrector equations in fluid mechanics: Effective viscosity of colloidal suspensions, Archive for Rational Mechanics and Analysis **239**, 1025 (2021).

[29] M. Duerinckx and A. Gloria, Quantitative homogenization theory for random suspensions in steady stokes flow, Journal de l'École polytechnique—Mathématiques **9**, 1183 (2022).

[30] M. Duerinckx, Semi-dilute rheology of particle suspensions: derivation of doi-type models, Archive for Rational Mechanics and Analysis **248**, 94 (2024).

[31] B. M. Haines, I. S. Aranson, L. Berlyand, and D. A. Karpeev, Effective viscosity of dilute bacterial suspensions: a two-dimensional model, Physical biology **5**, 046003 (2008).

[32] B. M. Haines, A. Sokolov, I. S. Aranson, L. Berlyand, and D. A. Karpeev, Three-dimensional model for the effective viscosity of bacterial suspensions, Physical Review E—Statistical, Nonlinear, and Soft Matter Physics **80**, 041922 (2009).

[33] B. M. Haines, I. S. Aranson, L. Berlyand, and D. A. Karpeev, Effective viscosity of bacterial suspensions: A three-dimensional pde model with stochastic torque, Commun. Pure Appl. Anal **11**, 19 (2012).

[34] A. Girodroux-Lavigne, Derivation of an effective rheology for dilute suspensions of microswimmers, SIAM Journal on Mathematical Analysis **55**, 6305 (2023).

[35] E. Bertin, M. Droz, and G. Grégoire, Boltzmann and hydrodynamic description for self-propelled particles, Physical Review E—Statistical, Nonlinear, and Soft Matter Physics **74**, 022101 (2006).

[36] K. Kruse and F. Jülicher, Dynamics and mechanics of motor-filament systems, The European Physical Journal E **20**, 459 (2006).

[37] A. Baskaran and M. C. Marchetti, Enhanced diffusion and ordering of self-propelled rods, Physical review letters **101**, 268101 (2008).

[38] A. Baskaran and M. C. Marchetti, Nonequilibrium statistical mechanics of self-propelled hard rods, Journal of Statistical Mechanics: Theory and Experiment **2010**, P04019 (2010).

[39] K. Doubrovinski and K. Kruse, Self-organization in systems of treadmilling filaments, The European Physical Journal E **31**, 95 (2010).

[40] Y.-L. Chou and T. Ihle, Active matter beyond mean-field: Ring-kinetic theory for self-propelled particles, Physical
Review E 91, 022103 (2015).

[41] E. Bertin, A. Baskaran, H. Chaté, and M. C. Marchetti, Comparison between smoluchowski and boltzmann approaches for self-propelled rods, Physical Review E 92, 042141 (2015).

[42] M. Miesowicz, The three coefficients of viscosity of anisotropic liquids, Nature 158, 27 (1946).

[43] H. Kneppe, F. Schneider, and B. Schwesinger, Axisymmetrical flow of a nematic liquid crystal around a sphere, Molecular Crystals and Liquid Crystals 205, 9 (1991).

[44] H. Heuer, H. Kneppe, and F. Schneider, Flow of a nematic liquid crystal around a sphere, Molecular Crystals and Liquid Crystals Science and Technology. Section A. Molecular Crystals and Liquid Crystals 214, 43 (1992).

[45] R. Ozaki, M. Aoki, K. Yoshino, K. Toda, and H. Moritake, Effective viscosity for nematic-liquid-crystal viscosity measurement using a shear horizontal wave, Physical Review E—Statistical, Nonlinear, and Soft Matter Physics 81, 061703 (2010).

[46] L. Berlyand, V. Rybalko, et al., Getting acquainted with homogenization and multiscale (Springer, 2018).

[47] A. Bensoussan, J.-L. Lions, and G. Papanicolaou, Asymptotic analysis for periodic structures, Vol. 374 (American Mathematical Soc., 2011).

[48] V. A. Marchenko and E. Y. Khruslov, Homogenization of partial differential equations, Vol. 46 (Springer Science & Business Media, 2008).

[49] G. Allaire, Homogenization of the navier-stokes equations in open sets perforated with tiny holes i. abstract framework, a volume distribution of holes, Archive for Rational Mechanics and Analysis 113, 209 (1991).

[50] O. A. Oleinik, A. Shamaev, and G. Yosifian, Mathematic

ical problems in elasticity and homogenization (Elsevier, 2009).

[51] E. Sánchez-Palencia, Non-homogeneous media and vibration theory, Lecture Note in Physics, Springer-Verlag 320, 57 (1980).

[52] J. Gachelin, G. Mino, H. Berthet, A. Lindner, A. Rouselet, and É. Clément, Non-newtonian viscosity of escherichia coli suspensions, Physical review letters 110, 268103 (2013).

[53] H. M. López, J. Gachelin, C. Douarche, H. Auradou, and E. Clément, Turning bacteria suspensions into superfluids, Physical review letters 115, 028301 (2015).

[54] E. Clement, A. Lindner, C. Douarche, and H. Auradou, Bacterial suspensions under flow, The European Physical Journal Special Topics 225, 2389 (2016).

[55] V. A. Martinez, E. Clément, J. Arlt, C. Douarche, A. Dawson, J. Schwarz-Linek, A. K. Creppy, V. Škultéty, A. N. Morozov, H. Auradou, et al., A combined rheometry and imaging study of viscosity reduction in bacterial suspensions, Proceedings of the National Academy of Sciences 117, 2326 (2020).

[56] E. Guazzelli and J. F. Morris, A physical introduction to suspension dynamics, Vol. 45 (Cambridge University Press, 2011).

[57] S. Fürthauer, M. Neef, S. W. Grill, K. Kruse, and F. Jülicher, The taylor-couette motor: spontaneous flows of active polar fluids between two coaxial cylinders, New Journal of Physics 14, 023001 (2012).

[58] F. Jülicher, S. W. Grill, and G. Salbreux, Hydrodynamic theory of active matter, Reports on Progress in Physics 81, 076601 (2018).

## Appendix A: The derivation of the molecular field

In this Appendix, we derive Eqs. (11)-(13), which respectively express the dynamics of the director field $n^\epsilon$ and the components of the molecular field $h_\perp^\epsilon$ in terms of the phase $\psi^\epsilon$ of the director field.

In two dimensions, the director field can be expressed as $n^\epsilon = \{\cos(\psi^\epsilon), \sin(\psi^\epsilon)\}$. Inserting this expression into Eq. (7), one obtains

$$h^{\epsilon} = K \Delta n^{\epsilon} - \lambda^{\epsilon} n^{\epsilon} = n^{\epsilon} (-K |\nabla \psi^{\epsilon}|^2 - \lambda^{\epsilon}) + n_{\perp}^{\epsilon} (K \Delta \psi^{\epsilon}) \quad (47)$$

where $n_\perp^\epsilon = \{-\sin(\psi^\epsilon), \cos(\psi^\epsilon)\}$ is a vector perpendicular to $n^\epsilon$. Therefore the components of the molecular field are

$$h_{\perp}^{\epsilon} = n_{\perp}^{\epsilon} \cdot h^{\epsilon} = K \Delta \psi^{\epsilon} \quad (48)$$

$$h_{\parallel}^{\epsilon} = n^{\epsilon} \cdot h^{\epsilon} = -K |\nabla \psi^{\epsilon}|^2 - \lambda^{\epsilon} \quad (49)$$

Note that the Lagrange multiplier $\lambda^\epsilon$ associated with the constraint $|n^\epsilon|^2 = 1$ can be re-expressed such that $h_\parallel^\epsilon$ becomes a Lagrange multiplier.

Inserting the above expression of the director field into the Ericksen-Leslie equation (3), one obtains

$$n_{\perp}^{\epsilon} (\partial_t \psi^{\epsilon} + \mathbf{u}^{\epsilon} \cdot \nabla \psi^{\epsilon} + n_{\perp}^{\epsilon} A(\mathbf{u}^{\epsilon}) n^{\epsilon}) = \Gamma h^{\epsilon} - \nu D(\mathbf{u}^{\epsilon}) n^{\epsilon} \quad (50)$$

and its projections in the direction parallel and perpendicular to the director field takes the form

$$\partial_t \psi^\epsilon + \mathbf{u}^\epsilon \cdot \nabla \psi^\epsilon + n_\perp^\epsilon A(\mathbf{u}^\epsilon) n^\epsilon = \Gamma h_\perp^\epsilon - \nu n_\perp^\epsilon D(\mathbf{u}^\epsilon) n^\epsilon \quad (51)$$

$$\Gamma h_{\parallel}^{\epsilon} = \nu n^{\epsilon} D(\mathbf{u}^{\epsilon}) n^{\epsilon} \quad (52)$$

which correspond to Eqs. (11) and (12) in the main text. The dynamics of the phase $\psi^\epsilon$ are governed by Eq. (51), where $h_\perp^\epsilon$ is given by Eq. (48). The Lagrange multiplier $h_\parallel^\epsilon$ is fixed by Eq. (52).
## Appendix B: Two-scale expansion of the equations for director

In this Appendix, we perform the two-scale expansion for Eqs. (11)-(13) as well as for the boundary conditions (17)-(18). For more details on the two-scale expansion, we refer the reader to the main text and to references [46, 47].

Inserting the relations (21) into Eqs. (11)-(13) and Eqs. (17)-(18) and treating the derivative as in Eq. (20), one obtains to the lowest order in $\epsilon$ that

$$ \epsilon^{-2}: \Gamma h_{\perp}^{(-2)} = 0 \qquad x \in \Omega^{\epsilon} \qquad (53) $$

$$ \epsilon^{-2}: h_{\perp}^{(-2)} = K \Delta^{y} \psi^{(0)} \qquad x \in \Omega^{\epsilon} \qquad (54) $$

$$ \epsilon^{-2}: \Gamma h_{\parallel}^{(-2)} = 0 \qquad x \in \Omega^{\epsilon} \qquad (55) $$

$$ \epsilon^{-1}: K(\hat{N} \cdot \nabla^y)\psi^{(0)} = 0 \qquad x \in \partial\mathcal{P}^{\epsilon} \qquad (56) $$

$$ \epsilon^{0}: \psi^{(0)} = 0 \qquad x \in \partial\Omega \qquad (57) $$

The solution of this set of PDEs is

$$ \psi^{(0)} = 0 \qquad (58) $$

The corrections to the next order in $\epsilon$ of Eqs. (11)-(13) and Eqs. (17)-(18) are

$$ \epsilon^{-1}: \Gamma h_{\perp}^{(-1)} = n_{\perp}^{(0)} A^{y}(u^{(0)}) n^{(0)} + \nu n_{\perp}^{(0)} D^{y}(u^{(0)}) n^{(0)} \qquad x \in \Omega^{\epsilon} \qquad (59) $$

$$ \epsilon^{-1}: h_{\perp}^{(-1)} = K \Delta^{y} \psi^{(1)} \qquad x \in \Omega^{\epsilon} \qquad (60) $$

$$ \epsilon^{-1}: \Gamma h_{\parallel}^{(-1)} = \nu n^{(0)} D^{y}(u^{(0)}) n^{(0)} \qquad x \in \Omega^{\epsilon} \qquad (61) $$

$$ \epsilon^{0}: K(\hat{N} \cdot \nabla^y)\psi^{(1)} = 0 \qquad x \in \partial\mathcal{P}^{\epsilon} \qquad (62) $$

$$ \epsilon^{1}: \psi^{(1)} = 0 \qquad x \in \partial\Omega \qquad (63) $$

The advection in Eq. (11) does not appear at order $\epsilon^{-1}$ because the lowest order term of the director's phase $\psi^{(0)}$ is independent on the fast variable.

In Appendix D, we will show that the lowest order term of the velocity field $u^{(0)}$ is only a function of the slow variable $x$. Note that in this case Eqs. (59) and (61) simplify to $h_{\perp}^{(-1)} = h_{\parallel}^{(-1)} = 0$, and therefore the solution of the previous set of PDEs is

$$ \psi^{(1)} = 0 \qquad (64) $$

Finally, the corrections to order $\epsilon^0$ of Eqs. (11) and (12) are

$$ \Gamma h_{\perp}^{(0)} = n_{\perp}^{(0)} (A^{x}(u^{(0)}) + A^{y}(u^{(1)})) n^{(0)} + \nu n_{\perp}^{(0)} (D^{x}(u^{(0)}) + D^{y}(u^{(1)})) n^{(0)} \qquad (65) $$

$$ \Gamma h_{\parallel}^{(0)} = \nu n^{(0)} (D^{x}(u^{(0)}) + D^{y}(u^{(1)})) n^{(0)} \qquad (66) $$

In our special case, both the lowest order term and the first correction to the director phase vanish $\psi^{(0)} = \psi^{(1)} = 0$ and thereby the advection in Eq. (11) does not appear at order $\epsilon^0$.

## Appendix C: Two-scale expansion of the stress

In this Appendix, we perform the two-scale expansion for the total stress tensor (5). The results from Appendix B are used to express the molecular field $h_{\epsilon}$ in terms of velocity field gradients. The two-scale expansion of the total stress is performed up to first order corrections in $\epsilon$.

We start from the total stress in index notation and perform the two-scale expansion in component-wise. From (5)

$$ \sigma_{ij}^{\epsilon}(x) = \eta(\partial_i u_j^{\epsilon} + \partial_j u_i^{\epsilon}) + \frac{\nu}{2}(n_i^{\epsilon} h_j^{\epsilon} + h_i^{\epsilon} n_j^{\epsilon}) + \frac{1}{2}(n_i^{\epsilon} h_j^{\epsilon} - h_i^{\epsilon} n_j^{\epsilon}) - K(\partial_i n_j^{\epsilon})(\partial_i n_j^{\epsilon}) - p^{\epsilon} \delta_{ij} \qquad (67) $$

We use the notation $\partial_i^{\alpha}$ where $i = 1, 2$ is the vector index and $\alpha = x, y$ denotes either the slow or fast variable respectively. Correspondingly for a function $u_k^{(l)}, k = 1, 2$ denotes the vector index and $l$ denotes the order of $\epsilon$ in the
two-scale expansion. In addition, plugging (47) into the total stress, the second and third terms on the RHS of (67) are

$$ \frac{\nu}{2}(n_i^{\epsilon} h_j^{\epsilon} + h_i^{\epsilon} n_j^{\epsilon}) = \frac{\nu}{2} (2h_{\parallel}^{\epsilon} n_i^{\epsilon} n_j^{\epsilon} + h_{\perp}^{\epsilon} (n_i^{\epsilon} (n_{\perp})_j^{\epsilon} + (n_{\perp})_i^{\epsilon} n_j^{\epsilon})) \quad (68) $$

$$ \frac{1}{2}(n_i^{\epsilon} h_j^{\epsilon} - h_i^{\epsilon} n_j^{\epsilon}) = \frac{1}{2} h_{\perp}^{\epsilon} (n_i^{\epsilon} (n_{\perp})_j^{\epsilon} - (n_{\perp})_i^{\epsilon} n_j^{\epsilon}) \quad (69) $$

which together with the solution $\psi^{(0)} = 0$ from Appendix B, simplifies (67) at order $\epsilon^{-1}$ to

$$ \sigma_{11}^{(-1)} = 2\eta \partial_1^y u_1^{(0)} + \nu h_{\parallel}^{(-1)} - p^{(-1)} \quad (70) $$

$$ \sigma_{22}^{(-1)} = 2\eta \partial_2^y u_2^{(0)} - p^{(-1)} \quad (71) $$

$$ \sigma_{12}^{(-1)} = \eta (\partial_1^y u_2^{(0)} + \partial_2^y u_1^{(0)}) + \frac{\nu+1}{2} h_{\perp}^{(-1)} \quad (72) $$

$$ \sigma_{21}^{(-1)} = \eta (\partial_2^y u_1^{(0)} + \partial_1^y u_2^{(0)}) + \frac{\nu-1}{2} h_{\perp}^{(-1)} \quad (73) $$

Next, we express Eqs. (59) and (61) from Appendix B in components and obtain the equations for $h_{\parallel}^{(-1)}$ and $h_{\perp}^{(-1)}$ in terms of velocity gradients:

$$ \Gamma h_{\parallel}^{(-1)} = \nu \partial_1^y u_1^{(-1)} \quad (74) $$

$$ \Gamma h_{\perp}^{(-1)} = \frac{(\nu - 1)}{2} \partial_1^y u_2^{(0)} + \frac{(\nu + 1)}{2} \partial_2^y u_1^{(0)} \quad (75) $$

Finally, substituting (74) and (75) into Eqs. (70)-(73), one obtains at order $\epsilon^{-1}$:

$$ \sigma_{11}^{(-1)} = (2\eta + \frac{\nu^2}{\Gamma}) \partial_1^y u_1^{(0)} - p^{(-1)} \quad (76) $$

$$ \sigma_{22}^{(-1)} = 2\eta \partial_2^y u_2^{(0)} - p^{(-1)} \quad (77) $$

$$ \sigma_{12}^{(-1)} = \left(\eta + \frac{\nu^2 - 1}{4\Gamma}\right) \partial_1^y u_2^{(0)} + \left(\eta + \frac{(\nu + 1)^2}{4\Gamma}\right) \partial_2^y u_1^{(0)} \quad (78) $$

$$ \sigma_{21}^{(-1)} = \left(\eta + \frac{(\nu - 1)^2}{4\Gamma}\right) \partial_1^y u_2^{(0)} + \left(\eta + \frac{\nu^2 - 1}{4\Gamma}\right) \partial_2^y u_1^{(0)} \quad (79) $$

Through a similar procedure, one obtains the corrections of the stress components at order $\epsilon^0$, which take the reduced form:

$$ \sigma_{11}^{(0)} = (2\eta + \frac{\nu^2}{\Gamma}) (\partial_1^x u_1^{(0)} + \partial_1^y u_1^{(1)}) - p^{(0)} \quad (80) $$

$$ \sigma_{22}^{(0)} = 2\eta (\partial_2^x u_2^{(0)} + \partial_2^y u_2^{(1)}) - p^{(0)} \quad (81) $$

$$ \sigma_{12}^{(0)} = \left(\eta + \frac{\nu^2 - 1}{4\Gamma}\right) \left(\partial_1^x u_2^{(0)} + \partial_1^x u_2^{(1)}\right) + \left(\eta + \frac{(\nu + 1)^2}{4\Gamma}\right) \left(\partial_2^x u_1^{(0)} + \partial_2^x u_1^{(1)}\right) \quad (82) $$

$$ \sigma_{21}^{(0)} = \left(\eta + \frac{(\nu - 1)^2}{4\Gamma}\right) \left(\partial_1^x u_2^{(0)} + \partial_1^x u_2^{(1)}\right) + \left(\eta + \frac{\nu^2 - 1}{4\Gamma}\right) \left(\partial_2^x u_1^{(0)} + \partial_2^x u_1^{(1)}\right) \quad (83) $$

where we used Eqs. (65) and (66) from Appendix B in components:

$$ \Gamma h_{\parallel}^{(0)} = \nu (\partial_1^x u_1^{(0)} + \partial_1^y u_1^{(1)}) \quad (84) $$

$$ \Gamma h_{\perp}^{(0)} = \frac{(\nu - 1)}{2} (\partial_1^x u_2^{(0)} + \partial_1^y u_2^{(1)}) + \frac{(\nu + 1)}{2} (\partial_2^x u_1^{(0)} + \partial_2^y u_1^{(1)}) \quad (85) $$

Finally, because the pressures $p^{(-1)}$ and $p^{(0)}$ are Lagrange multipliers, we redefine them as follows in order to simplify the stress components in Eqs. (80)-(83):

$$ \bar{p}^{(-1)} = p^{(-1)} - \frac{\nu(\nu+1)}{2\Gamma} \partial_1^x u_1^{(0)} = p^{(-1)} + \frac{\nu(\nu+1)}{2\Gamma} \partial_2^y u_2^{(0)} \quad (86) $$

$$ \bar{p}^{(0)} = p^{(0)} - \frac{\nu(\nu+1)}{2\Gamma} (\partial_1^x u_1^{(0)} + \partial_1^y u_1^{(1)}) = p^{(0)} + \frac{\nu(\nu+1)}{2\Gamma} (\partial_2^x u_2^{(0)} + \partial_2^y u_2^{(1)}) \quad (87) $$
where the equality follows from the incompressibility conditions Eqs. (98) and (104). Then, the stress components in Eqs. (76)-(79)

$$ \sigma_{11}^{(-1)} = \left(2\eta + \frac{\nu^2 - \nu}{2\Gamma}\right) \partial_1^y u_1^{(0)} - \bar{p}^{(-1)} \quad (88) $$

$$ \sigma_{22}^{(-1)} = \left(2\eta + \frac{\nu^2 + \nu}{2\Gamma}\right) \partial_2^y u_2^{(0)} - \bar{p}^{(-1)} \quad (89) $$

$$ \sigma_{12}^{(-1)} = \left(\eta + \frac{\nu^2 - 1}{4\Gamma}\right) \partial_1^y u_2^{(0)} + \left(\eta + \frac{(\nu + 1)^2}{4\Gamma}\right) \partial_2^y u_1^{(0)} \quad (90) $$

$$ \sigma_{21}^{(-1)} = \left(\eta + \frac{(\nu - 1)^2}{4\Gamma}\right) \partial_1^y u_2^{(0)} + \left(\eta + \frac{\nu^2 - 1}{4\Gamma}\right) \partial_2^y u_1^{(0)} \quad (91) $$

and Eqs. (80)-(83) become

$$ \sigma_{11}^{(0)} = \left(2\eta + \frac{\nu^2 - \nu}{2\Gamma}\right) \left(\partial_1^x u_1^{(0)} + \partial_1^y u_1^{(1)}\right) - \bar{p}^{(0)} \quad (92) $$

$$ \sigma_{22}^{(0)} = \left(2\eta + \frac{\nu^2 + \nu}{2\Gamma}\right) \left(\partial_2^x u_2^{(0)} + \partial_2^y u_2^{(1)}\right) - \bar{p}^{(0)} \quad (93) $$

$$ \sigma_{12}^{(0)} = \left(\eta + \frac{\nu^2 - 1}{4\Gamma}\right) \left(\partial_1^x u_2^{(0)} + \partial_1^y u_2^{(1)}\right) + \left(\eta + \frac{(\nu + 1)^2}{4\Gamma}\right) \left(\partial_2^x u_1^{(0)} + \partial_2^y u_1^{(0)}\right) \quad (94) $$

$$ \sigma_{21}^{(0)} = \left(\eta + \frac{(\nu - 1)^2}{4\Gamma}\right) \left(\partial_1^x u_2^{(0)} + \partial_1^y u_2^{(1)}\right) + \left(\eta + \frac{\nu^2 - 1}{4\Gamma}\right) \left(\partial_2^x u_1^{(0)} + \partial_2^y u_1^{(0)}\right) \quad (95) $$

## Appendix D: Computing the force balance PDE

In this Appendix, we use the results in Appendix C to perform the two-scale expansion for Eqs. (9)-(10) as well as for the boundary conditions (15)-(16) up to the first order corrections in $\epsilon$.

Inserting the relations (21) into Eqs. (9)-(10) and Eqs. (15)-(16) and treating the derivative as in Eq. (20), one obtains to the lowest order in $\epsilon$ that

$$ \varepsilon^{-2}: \quad \eta_1 \partial_{11}^y u_1^{(0)} + \eta_2 \partial_{22}^y u_1^{(0)} - \partial_1^y \bar{p}^{(-1)} = 0 \qquad x \in \Omega^{\epsilon} \quad (96) $$

$$ \varepsilon^{-2}: \quad \eta_1 \partial_{11}^y u_2^{(0)} + \eta_2 \partial_{22}^y u_2^{(0)} - \partial_2^y \bar{p}^{(-1)} = 0 \qquad x \in \Omega^{\epsilon} \quad (97) $$

$$ \varepsilon^{-1}: \quad \partial_1^y u_1^{(0)} + \partial_2^y u_2^{(0)} = 0 \qquad x \in \Omega^{\epsilon} \quad (98) $$

$$ \varepsilon^{-1}: \quad \mathbf{u}^{(0)} = \mathbf{E} \cdot \mathbf{x} \qquad x \in \partial \mathcal{P}^{\epsilon} \quad (99) $$

$$ \varepsilon^0: \quad \mathbf{u}^{(0)} = \mathbf{E} \cdot \mathbf{x} \qquad x \in \partial \Omega \quad (100) $$

where we used the stress components (70)-(73) in Appendix C. The parameters $\eta_1 = \eta + (\nu - 1)^2/4\Gamma$ and $\eta_2 = \eta + (\nu + 1)^2/4\Gamma$ are two effective viscosities. The solution of this set of PDEs is

$$ \mathbf{u}^{(0)} = \mathbf{E} \cdot \mathbf{x} \quad (101) $$

Note that because the strain rate $\mathbf{E}$ represents a simple shear flow, it is traceless and symmetric. Therefore $\partial_1^x u_1^{(0)} + \partial_2^x u_2^{(0)} = 0$.

Since the velocity and pressure fields are $\epsilon$-periodic in the fast variable $\mathbf{y}$, the corrections to the next order in $\epsilon$ of Eqs. (9)-(10) and Eqs. (15)-(16) are

$$ \varepsilon^{-1}: \quad \eta_1 \partial_{11}^y u_1^{(1)} + \eta_2 \partial_{22}^y u_1^{(1)} - \partial_1^y \bar{p}^{(-1)} = 0 \qquad y \in Y_f \quad (102) $$

$$ \varepsilon^{-1}: \quad \eta_1 \partial_{11}^y u_2^{(1)} + \eta_2 \partial_{22}^y u_2^{(1)} - \partial_2^y \bar{p}^{(-1)} = 0 \qquad y \in Y_f \quad (103) $$

$$ \varepsilon^0: \quad \partial_1^y u_1^{(1)} + \partial_2^y u_2^{(1)} = 0 \qquad y \in Y_f \quad (104) $$

$$ \varepsilon^0: \quad \mathbf{u}^{(1)} = -\mathbf{E} \cdot \mathbf{y} \qquad y \in \partial \mathcal{P} \quad (105) $$

$$ \varepsilon^1: \quad \mathbf{u}^{(1)} \text{ periodic} \qquad y \in \partial \partial Y_f \quad (106) $$
Because the fluid is incompressible, the velocity field can be expressed in terms of a streamfunction such that $u_1^{(1)} = -\partial_2^y \Psi^{(1)}$ and $u_2^{(1)} = \partial_1^y \Psi^{(1)}$. Taking the curl of (102) and (103), one obtains that $\Psi^{(1)}$ solves an anisotropic biharmonic equation

$$ (\eta_1 \partial_{11}^y + \eta_2 \partial_{22}^y)(\partial_{11}^y + \partial_{22}^y)\Psi^{(1)} = 0 \qquad (107) $$

Taking the divergence of (102) and (103), one obtains that $\bar{p}^{(0)}$ solves the Laplace equation

$$ (\partial_{11}^y + \partial_{22}^y)\bar{p}^{(0)} = 0 \qquad (108) $$

The set of PDEs (102)-(106) is solved using asymptotic expansions in Appendix E.

## Appendix E: Computing the Weak Anisotropic Limit:

In this Appendix, we perform the asymptotic expansion of the physical fields in the limits $\nu \ll 1$ and $\nu \gg 1$ in the anisotropic Stokes Eqs. (102)-(106) obtained in Appendix D. This reduces the PDEs to a set of isotropic Stokes equations with forcing terms.

We start from expanding the velocity and pressure field in the anisotropic Stokes Eq. (25)-(27) as a power series in small shear-flow alignment $|\nu| \ll 1$ (31)

$$ \mathbf{u}_s^{(1)} = \sum_{k=0}^{\infty} \nu^k \mathbf{u}_s^{(1,k)} \qquad (109) $$

$$ \bar{p}_s^{(0)} = \sum_{k=0}^{\infty} \nu^k \bar{p}_s^{(0,k)} \qquad (110) $$

where we denote the $\varepsilon$-order $\alpha$ and $\nu$-order $\beta$ function by $\mathbf{u}^{(\alpha,\beta)}$. We collect terms at each order $\nu^k$

$$ \nu^0: \begin{cases} \bar{\eta} \Delta \mathbf{u}_s^{(1,0)} - \nabla \bar{p}_s^{(0,0)} = 0 \\ \nabla \cdot \mathbf{u}_s^{(1,0)} = 0 \end{cases} \qquad y \in Y_f \qquad (111) $$

$$ \nu^1: \begin{cases} \bar{\eta} \Delta \mathbf{u}_s^{(1,1)} - \nabla \bar{p}_s^{(0,1)} + \frac{1}{2\Gamma} (\partial_{22}^y \mathbf{u}_s^{(1,0)} - \partial_{11}^y \mathbf{u}_s^{(1,0)}) = 0 \\ \nabla \cdot \mathbf{u}_s^{(1,1)} = 0 \end{cases} \qquad y \in Y_f \qquad (112) $$

$$ \nu^2: \begin{cases} \bar{\eta} \Delta \mathbf{u}_s^{(1,2)} - \nabla \bar{p}_s^{(0,2)} + \frac{1}{2\Gamma} (\partial_{22}^y \mathbf{u}_s^{(1,1)} - \partial_{11}^y \mathbf{u}_s^{(1,1)}) \\ + \frac{1}{4\Gamma} (\partial_{11}^y \mathbf{u}_s^{(1,0)} + \partial_{22}^y \mathbf{u}_s^{(1,0)}) = 0 \\ \nabla \cdot \mathbf{u}_s^{(1,2)} = 0 \end{cases} \qquad y \in Y_f \qquad (113) $$

where $\bar{\eta} = \eta + 1/4\Gamma$. We solve this system by a streamfunction $\Psi_s^{(1,k)}$

$$ \Psi_s^{(1,k)} := A^{(k)}(r) \sin(2(k+1)\theta) \qquad (114) $$

Furthermore, we expand Eq. (25)-(27) as a power series in large shear-flow alignment $|\nu| \gg 1$. From Eq. (34)

$$ \mathbf{u}_l^{(1)} = \sum_{k=1}^{\infty} \nu^{-k} \mathbf{u}_l^{(1,k)} \qquad (115) $$

$$ \bar{p}_l^{(0)} = \sum_{k=0}^{\infty} \nu^{-k} \bar{p}_l^{(0,k)} \qquad (116) $$
$$ \nu^2: \begin{cases} \bar{\eta} \Delta \mathbf{u}_l^{(1,2)} - \nabla \bar{p}_l^{(0,2)} = 0 \\ \nabla \cdot \mathbf{u}_l^{(1,2)} = 0 \end{cases} \qquad y \in Y_f \qquad (117) $$

$$ \nu^1: \begin{cases} \bar{\eta} \Delta \mathbf{u}_l^{(1,1)} - \nabla \bar{p}_l^{(0,1)} + \frac{1}{2\Gamma} (\partial_{22}^y \mathbf{u}_l^{(1,1)} - \partial_{11}^y \mathbf{u}_l^{(1,1)}) = 0 \\ \nabla \cdot \mathbf{u}_l^{(1,1)} = 0 \end{cases} \qquad y \in Y_f \qquad (118) $$

$$ \nu^0: \begin{cases} \bar{\eta} \Delta \mathbf{u}_l^{(1,0)} - \nabla \bar{p}_l^{(0,0)} + \frac{1}{2\Gamma} (\partial_{22}^y \mathbf{u}_l^{(1,1)} - \partial_{11}^y \mathbf{u}_l^{(1,1)}) = 0 \\ + \frac{1}{4\Gamma} (\partial_{11}^y \mathbf{u}_l^{(1,2)} + \partial_{22}^y \mathbf{u}_l^{(1,2)}) = 0 \\ \nabla \cdot \mathbf{u}_l^{(1,0)} = 0 \end{cases} \qquad y \in Y_f \qquad (119) $$

where we use the streamfunction $\Psi_l^{(1,k)}$

$$ \Psi_l^{(1,k)} := B^{(k)}(r) \sin(2(k+1)\theta), \quad k = 0, 1, 2 \qquad (120) $$

Using the scaling (32), we pass to the dilute limit in both (111)-(113) and (117)-(119) by taking the particle radius $a \to 0$ as $\epsilon \to 0$, see Fig. 2. In this limit, we rescale the domain $Y_f$ such that it becomes an infinite domain $\mathbb{R}^2 \setminus \mathcal{P}$ with a particle of size 1 and we replace the periodic boundary condition with a decay condition $\mathbf{u}^{(1)} \to 0$ at infinity. Solving for the coefficients $A^{(k)}, B^{(k)}$ in the streamfunction ansatze (114) and (120) in the infinite domain $\mathbb{R}^2 \setminus \mathcal{P}$, we obtain (33) and (35), respectively, after rescaling back to the unit cell domain $Y_f$.

We solve explicitly the velocity field $\mathbf{u}^{(1)}$ by finding for either streamfunction $\Psi_s^{(1,k)}$ or $\Psi_l^{(1,k)}$

$$ u_1^{(1,k)} = -\partial_2^y \Psi^{(1,k)}, \quad u_2^{(1,k)} = \partial_1^y \Psi^{(1,k)}, \quad k = 0, 1, 2 \qquad (121) $$

Due to the incompressibility constraint the pressure satisfies a Laplace equation at each order

$$ \Delta \bar{p}^{(0,k)} = 0, \quad k = 0, 1, 2 \qquad (122) $$

Note that for the 2D isotropic Stokes Eqs. (111)-(113) and Eqs. (117)-(119), the solution to the problem with rigidity condition (15) coincides (see for example [56], Section 2) with the solution with a no-slip condition boundary condition

$$ \mathbf{u} = 0 \qquad (123) $$

where the velocity of the particle is assumed to be 0, see Appendix F.

## Appendix F: Net force and net torque on the particle boundary:

In this Appendix, we determine the net force and net torque on the particle boundary as a result of the perturbations of the physical fields. For these calculations, we will use the result from Appendix B to E. For more details on the derivation of the net force and net torque on the interface of a liquid crystal, we refer to [57, 58].

Because the particles are undeformable, their motion is limited to translations and rotations. The linear and angular velocities of a particle are determined by enforcing that the net force and the net torque on the particle boundary vanishes.

On the interface of a liquid crystal and in the absence of external sources, the expressions of the net force and net torque are

$$ \int_{\partial \mathcal{P}^{\epsilon}} \boldsymbol{\sigma}^{\epsilon} \cdot \hat{\mathbf{N}} dl = 0 \qquad (124) $$

$$ \int_{\partial \mathcal{P}^{\epsilon}} \mathbf{r} \times (\boldsymbol{\sigma}_{\tilde{d}}^{\epsilon} \cdot \hat{\mathbf{N}}) dl = 0 \qquad (125) $$
where $\sigma^\epsilon$ is the total stress tensor, $\hat{N}$ is the outward facing normal vector, and

$$ \sigma_d^\epsilon = 2\eta D(\mathbf{u}^\epsilon) + \frac{\nu}{2}(\mathbf{n}^\epsilon \mathbf{h}^\epsilon + \mathbf{h}^\epsilon \mathbf{n}^\epsilon) + \frac{1}{2}(\mathbf{n}^\epsilon \mathbf{h}^\epsilon - \mathbf{h}^\epsilon \mathbf{n}^\epsilon) \quad (126) $$

and $dl$ is an infinitesimal element of length on the boundary, which in our case is the boundary of a particle.

The two-scale expansion of the stress components of $\sigma^\epsilon$ and (126) can be determined using the results in Appendix C. To order $\epsilon^0$, the components of $\sigma^\epsilon$ are given by Eqs. (92)-(95) and the components of takes a similar form

$$ \sigma_{d,11}^{(0)} = (2\eta + \frac{\nu^2}{\Gamma}) (\partial_1^x u_1^{(0)} + \partial_1^y u_1^{(1)}) \quad (127) $$

$$ \sigma_{d,22}^{(0)} = 2\eta (\partial_2^x u_2^{(0)} + \partial_2^y u_2^{(1)}) \quad (128) $$

$$ \sigma_{d,12}^{(0)} = (\eta + \frac{\nu^2 - 1}{4\Gamma}) (\partial_1^x u_2^{(0)} + \partial_1^x u_2^{(1)}) + (\eta + \frac{(\nu + 1)^2}{4\Gamma}) (\partial_2^x u_1^{(0)} + \partial_2^x u_1^{(1)}) \quad (129) $$

$$ \sigma_{d,21}^{(0)} = (\eta + \frac{(\nu - 1)^2}{4\Gamma}) (\partial_1^x u_2^{(0)} + \partial_1^x u_2^{(1)}) + (\eta + \frac{\nu^2 - 1}{4\Gamma}) (\partial_2^x u_1^{(0)} + \partial_2^x u_1^{(1)}) \quad (130) $$

Using the expression for $\mathbf{u}^{(0)}$ in Eq. (101) and the expression for $\mathbf{u}^{(1)}$ in Eqs. (33) or (35), one can show that the net force and net torque vanish in the two limiting cases controlled by $\nu$. As a consequence the particle velocity vanishes.

## Appendix G: Macroscopic field of an homogenized medium:

In this Appendix, we determine the form of the physical fields of an homogenized liquid crystal medium that is subjected to the same extensional shear flows and preferred anchoring as the suspension (17) and (18).

When the homogenized medium is subjected to an external shear flow, each particle can create a disturbance on the macroscopic scales that can influence the and the integrate effect can generate corrections to the macroscopic flows. Following [26], we consider that the flow field of the homogenized medium is

$$ \mathbf{u}^H = \mathbf{E} \cdot \mathbf{x} + \sum_{n=1}^{N_p} \mathbf{u}^{(1)}(\mathbf{x}, \mathbf{y} = \mathbf{x} - \mathbf{x}_i) \quad (131) $$

where the summation runs over all particles in the medium and $\mathbf{x}_i$ is the center of the i-th particle. Because in our case, particles do not perturb the director field (see Appendix B), the corrections of an ensemble of particles to the macroscopic director field vanishes.

The corrected macroscopic shear flow is defined as

$$ E_{\alpha\beta}^{H} = E_{\alpha\beta} + \left( \sum_{n=1}^{N_p} \partial_{x_{\alpha}} (u_{\beta}^{(1)}(\mathbf{x}, \mathbf{y} = \mathbf{x} - \mathbf{x}_i)) \right) \Big|_{\mathbf{x}=0} \quad (132) $$

where the first term on the RHS represents the macroscopic flows and the second term on the RHS represents the sum of the disturbance flows generated by each particle. Considering that the particles are evenly distributed and taking the continuum limit, we can express (132) as

$$ E_{\alpha\beta}^{H} = E_{\alpha\beta} + \frac{N_p}{|\Omega|} \int_{\Omega} \partial_{x_{\alpha}} (u_{\beta}^{(1)}(\mathbf{x}, \mathbf{y} = \mathbf{x} - \mathbf{x}'))|_{\mathbf{x}=0} d\mathbf{x}' \quad (133) $$

$$ = E_{\alpha\beta} - \frac{N_p}{|\Omega|} \int_{\Omega} u_{\beta}^{(1)}(\mathbf{x}, \mathbf{y} = \mathbf{x} - \mathbf{x}')|_{\mathbf{x}=0} x'_{\alpha} d\mathbf{x}' \quad (134) $$

Substituting the expression of the perturbation flows $u_{\alpha}^{(1)}$ in the two limits of $\nu$, Eqs. (33) or (35), leads to the same expression for the homogenized macroscopic shear flow

$$ E^H = \begin{pmatrix} \gamma(1-\phi) & 0 \\ 0 & -\gamma(1-\phi) \end{pmatrix} \quad (135) $$