Elementary Flow Field Profiles of Micro-Swimmers in Weakly Anisotropic Nematic Fluids : Stokeslet , Stresslet , Rotlet and Source Flows

Analytic formulations of elementary flow field profiles in weakly anisotropic nematic fluid are determined, which can be attributed to biological or artificial micro-swimmers, including Stokeslet, stresslet, rotlet and source flows. Stokes equation for a nematic stress tensor is written with the Green function and solved in the k-space for anisotropic Leslie viscosity coefficients under the limit of leading isotropic viscosity coefficient. Analytical expressions for the Green function are obtained that are used to compute the flow of monopole or dipole swimmers at various alignments of the swimmers with respect to the homogeneous director field. Flow profile is also solved for the flow sources/sinks and source dipoles showing clear emergence of anisotropy in the magnitude of flow profile as the result of fluid anisotropic viscosity. The range of validity of the presented analytical solutions is explored, as compared to exact numerical solutions of the Stokes equation. This work is a contribution towards understanding elementary flow motifs and profiles in fluid environments that are distinctly affected by anisotropic viscosity, offering analytic insight, which could be of relevance to a range of systems from microswimmers, active matter to microfluidics.


Introduction
Microswimmers are biological or artificial entities that are characterized by locomotion at microscopic scales, such as bacteria, protozoa, spermatozoa, algae, or Janus colloids, self-propelling droplets, shape-shifting swimmers, and rotating colloids [1][2][3].The locomotion of the swimmers is affected by their relatively small size (order of magnitude 10 µm) and relatively slow propulsion (up to few 100 µm/s) [1,3], which determines that their dynamics is described by the Stokes equation (rather than the full Navier-Stokes), as having very small Reynolds numbers (Re 1).As recognized by Purcell, the swimming mechanism of such microswimmers requires swimming strokes that are irreversible in time [4].While swimming of biological microswimmers is an essential part of their survival and reproduction, artificial microswimmers not only help to understand the biomechanics behind swimming organisms [5], but could possibly be also utilized to drive gears (produce work) [6,7], enhance the diffusion of a solution [8], increase stirring by the process of particle entrainment [9] or produce frictionless fluids [10,11].
Microswimmers are often faced with non-Newtonian environment.For example, biological fluids, such as blood, mucus, and saliva show viscoelastic behavior [12], and cytoskeleton of animal cells shows orientational order [13].Furthermore, many biomaterials can form a nematic phase [14].To investigate the influence of anisotropic background on swimming bacteria, studies of bacterial solutions in lyotropic liquid crystals were conducted, investigating the behavior of single bacteria [15][16][17], their pairwise interaction [18] or their collective motion [19,20].Bacterial suspensions in liquid crystals can be used for transport of microcargo [21].While investigating swimming bacteria in liquid crystals is particularly challenging due to high toxicity of many nematics, this problem is avoided by studying artificial active particles in liquid crystals [22], or, for example, properties of a Taylor's swimming sheet [23,24].Even passive systems, such as sedimentation of spheres shows an intertwined behavior of fluid velocity and internal structure of nematics.Flow profiles in a nematic liquid crystal were discussed typically in the context of a moving sphere inside a nematic for a case of homogeneous [25,26] or inhomogeneous [27][28][29] director profile.In [30], the solution for the flow profile and the drag force on a moving sphere were determined in the weakly anisotropic limit.In [31], the problem of nematic flow around a sphere was investigated by the use of Green functions for the Stokes equation, which is a method that is applied also in this article.Biological and artificial swimmers in non-Newtonian fluids, possibly with a anisotropic structure, represent a further complexity in their underlying physics.Finally, understanding the effects of the internal structure of the complex fluids on the microswimmers is an open challenge in both living systems as well as designed devices.
The flow fields caused by the microswimmers can be characterized by the moments of force that microswimmers apply on a fluid [1][2][3].If a microswimmer exerts an effective point force on a fluid, we speak of a Stokeslet flow.Such flows have been observed in the case of Volvox carteri [32] and externally driven colloids [5].A microswimmer that exerts an effective force dipole on a fluid produces a stresslet or rotlet flow.Stresslet is a typical flow profile for swimming microorganisms, such as Escherichia coli [33] or self-driven colloids [5].Although multipolar approach ignores the details of the flow field close to the swimmers, for example around the flagella, its particular advantage is that higher order multipoles decay faster with the distance from the swimmer.Therefore, it is especially efficient in semi-dilute suspensions, where multipolar expansion of the flow field was used to investigate problems such as fluid mixing by the microswimmers [34], interaction between microswimmers [35,36], and collective behaviour in microswimmer suspensions [37,38].The swimmer flow field can be additionally affected by the confinement, as, for example, in thin films [39] or droplets [40], where the orientation and position of the swimmer with respect to the boundary determines the flow field and the behavior of microswimmers.In dense suspensions, the dominating interactions are steric, electrostatic, van der Waals and close-field hydrodynamic [3].In such high-density suspensions, swimmers tend to move in swarms, or can even produce an active nematic phase [41].
In this paper, we determine the analytical formula for flow fields of micro-swimmers-represented as point forces, point torques, or stresslets-in anisotropic nematic fluids, under the assumption of uniform nematic alignment and within the first order expansion of the anisotropy of the viscosity tensor.Specifically, we first calculate analytical approximations of the Green function to the Stokes equation with anisotropic viscosity, included via the nematic-type stress formulation, and then use these solutions to construct the flow fields of point force, point torque, and stresslet in uniform nematic fluid.Using a similar approach, we also determine analytically the flow field of point sinks and sources in a uniform nematic.Distinctly, we show how the direction of the applied forces and torques with respect to the nematic director affects the calculated flow profiles.The two key assumptions of the demonstrated analytical solutions-uniformity of nematic field and anisotropy expansion in viscosity-are discussed.Anisotropy expansion in viscosity is tested against precise numerical solutions.More generally, the results of this paper give analytical insight into elementary flow field profiles applicable to micro-swimmers in anisotropic fluids and could be used to further explore more complicated flow profiles and swimming dynamics in various complex environments.

Nematic Green Function for the Stokes Equation
A strong approach to the formulation of elementary flow profiles as generated by local stimuli-like micro-swimmers-is to first consider point-like forces acting on a fluid, which is known as the Method of Green's functions [42].Specifically, in a distinct point in space, a point force F is taken to act on the fluid (in our case at r = 0), and then one tries to solve the corresponding flow field from the dynamic equation of the fluid.For fluid flows relevant to micro-swimmers, we are interested in the solutions of the Stokes equation with a point force acting on a fluid: where ∂ i is a spatial derivative along the i-th direction, p is the fluid pressure, F a point force, σ the stress tensor, and δ( r) the three-dimensional Dirac delta function.We also assume incompressible flow 1) needs to be solved for both the pressure field p( r) and the velocity field v( r).
The solution of Equation ( 1) is called Stokeslet and is well known for isotropic fluids [42], where the stress tensor is proportional to the strain rate σ 1).Often, the interaction between a microswimmer and the surrounding fluid is not best described by a single point force, but by a distribution of forces.In this case, the main advantage of the Green function approach is that, once the solution for the point force is found, it can be easily expanded to other force distribution, such as a pair of opposite forces (stresslet flow), or a point torque (rotlet flow), solutions which are shown in Figure 1 for the isotropic stress tensor.Here, we are interested in the same elementary solutions, however in the systems with the nematic degree of order, i.e., with anisotropic viscosity and some preferred homogeneous orientation (director) field of the system.The anisotropy of the system is reflected in the nematic stress tensor, which we write in the Ericksen-Leslie formulation [43]: where , n is the director, v the flow field velocity vector, and α i are six Leslie viscosity coefficients (five are independent due to constraint α 6 − α 5 = α 2 + α 3 ).In an isotropic liquid, only α 4 viscosity coefficient is non-zero, giving the standard isotropic fluid stress tensor.
Flow fields may cause realignment in the nematic orientational profile.The dynamics of the director profile as coupled to the velocity field is given by the equation [43] where f is a free energy density of the nematic due to elastic deformations of the director field and due to coupling to the external fields, and Λ is a Lagrange multiplier preserving the unit length of n.
In principle, a solution to the proposed problem of microswimmer flow profiles in nematic fluids could be found by solving the Stokes equation together with the incompressibility condition, and molecular field equation as a set of coupled nonlinear differential equations.Indeed, such solutions could be found by using numerical simulations [44,45], but not analytically.However, a particular question we want to address in this paper is whether some analytic insight can be found, which is quite rare in complex fluid systems.We search for a solution by assuming a homogeneous director field, therefore excluding Equation (3) from our calculations.Such analytic solutions typically prove distinct relevance as they can be used in a variety of experimental or computer modelling setups either for analysis or as elementary building solutions in a much more complex system.For example, one idea would be to take these elementary analytic solutions of individual flow elements and incorporate them as simplest-order approximations in systems of multiple flow objects-like colonies of swimmers.
The posed problem of elementary solutions in Stokes fluid with anisotropic viscosity can be solved analytically-within approximation-if considering the regimes of the anisotropic viscosity coefficients α i .The large variety of nematic materials-for example from molecular to virus and colloidal, thermotropic or lyotropic, passive or active-exhibit a variety of regimes of the actual values of the viscosity coefficients α i .For example, in some lyotropic nematics, α 5 − α 2 can be very large, whereas |α 1 |, |α 3 |, α 4 , and |α 6 | are smaller in size and α 1 , α 3 , and α 6 are negative [46][47][48].Differently, in a standard thermotropic nematic like 5 4-cyano-4'-pentylbiphenyl (5CB) α 2 , α 4 , α 5 are of similar order of magnitude which is slightly larger than the magnitude of α 6 , whereas α 1 and α 3 are smaller [49].However, a common trend in nematic systems is that if the nematic degree of order gets smaller (for example, as when approaching the nematic-isotropic phase transition), the (isotropic) viscosity coefficient α 4 typically can become the leading viscosity coefficient [50].Stimulated by this behaviour of α 4 , we assume in this work that α 4 is the largest coefficient and calculate the solutions to the Stokes equation that are of first order in other viscosity coefficients.From a more theoretical perspective, such approach of taking isotropic viscosity coefficient α 4 as the leading term in the stress tensor can be seen also as searching for solutions that are close to the isotropic fluid solutions but adapted (within first order) for anisotropic viscosity.An example of how to give physical meaning to viscosity coefficients is to compare effective viscosity at different orientations of velocity, velocity gradient, and director field-in the context of Miesowicz geometry [43].For the hierarchy of α 2 , α 3 , α 5 , and α 6 of lyotropic nematics discussed above, the lowest effective viscosity is in the case of a director being parallel to flow direction and perpendicular to flow gradient, whereas the highest viscosity is for a director parallel to flow gradient and perpendicular to flow direction.This regime of viscosity coefficients is used also for the graphic representation of the results throughout the paper, where we use values α 5 − α 2 = 0.6α 4 , α 1 = α 3 = α 6 = −0.15α 4 , although the calculations are for general values α i , as long as α 4 is the leading viscosity coefficient.
We solve the Stokes equation (Equation ( 2)) with the assumption of homogeneous director pointing along the z-axis and write the Stokes equation and incompressibility condition in terms of x, y, and z components: where ∂ ij is a second derivative over the i-th and j-th spatial coordinates.Due to linearity of the Stokes equation, we look for the solution in the form where G ij and P j are the Green functions for velocity and pressure, respectively.Following the standard approach of the method of Green's functions, we insert Equations ( 8) and ( 9) into the Stokes and continuity equations and make the Fourier transform of the equations, using and the fact that Green functions are independent on the force F j , we obtain the Green functions in the reciprocal space: where the caret symbol above Pj and Ĝij indicates solutions in the reciprocal space.We have also used the notation ) to indicate which of the viscosity coefficients are present in the Green functions.While α 4 only rescales the velocity, α1 , ν b , and ν c determine the actual shape of the velocity profile.ν d is included only in the pressure field.
In order to be able to perform the inverse Fourier transform of Equations ( 10), ( 11) and ( 13) we first expand them in the series of αi , where i = 4:

Ĝzz
The expressions for Ĝyz , Ĝyy , and Py can be deduced by simply switching the x-coordinate with y-coordinate in the above equations.To perform an inverse Fourier transform of the above Green functions, we employ a procedure, described in Ref. [51], based on the plane wave expansion of e i k• r and integral properties of spherical Bessel functions.This method allows for the calculation of inverse Fourier transforms of expressions in the form k i 1 k i 2 . . .k i N /k N+2 , among others, for any value of N and indices i by performing the angular momentum decomposition.Following the approach in Ref. [51], we calculate the list of inverse Fourier transforms in Table 1.Finally, we obtain for the Green functions the following expressions: Note that, if α 1 is much smaller than other viscosity coefficients, which is true for quite some nematic fluids (e.g., in nematic 5 4-cyano-4'-pentylbiphenyl (CB) at 30 • C α 1 /α 4 is ∼0.03 [49]), the nematic Green functions can be further simplified neglecting terms with α 1 .If only α 4 is non-zero, a solution reduces to the well known Green functions for an isotropic fluid [42].
Table 1.A list of inverse Fourier transforms used in this article, obtained by the procedure, described in [51].Note that the expressions form self-consistent pairs bound by the relation x 2 z 3 r 5 − 3 x 2 z r 3 − 5 Two notable conclusions about the properties of the flow fields of micro-swimmers can be drawn from the calculated Green functions in k-space (Equations ( 10)-( 13)): (i) Green function for the velocity decays as 1/r and for the pressure as 1/r 2 ; and (ii) Green functions are symmetrical.Furthermore, the calculated profiles exactly satisfy the incompressibility condition, despite being only an approximation to the exact solution of the Stokes equation.More generally, the calculated nematic Green functions now allow for the calculation of the flow and pressure fields at any force distribution and could possibly be also used in the presence of walls, where, in isotropic fluids, a set of mirror flow fields is constructed to accommodate the no-slip velocity condition at the wall [39,42].

Flow Fields of Point Force in Nematic Fluid
We use Green functions from Equations ( 24)- (32) to calculate the velocity profile in a weakly anisotropic homogeneous nematic due to a point force.In isotropic fluids, such Stokeslet profiles are associated with the dynamics of externally driven colloidal particles [3] or distinct microswimmers [32].Differently as in the isotropic fluids, in the case of anisotropic nematic fluids, there are two principal solutions to the flow profile, corresponding to force being parallel or perpendicular to the director.For any other angle, the flow profile can be calculated as a linear combination of the two principal solutions due to the linearity of the Stokes equation.The nematic flow field generated by a point force F can be written as: where G ij are Green functions given in Equations ( 24)-( 28).The calculated flow profiles are shown in Figure 2. Figure 2a,e shows the flow profile of a point force that is aligned with the director.Similar as in the isotropic case (Figure 1), the flow profile retains the rotational symmetry around the force direction; however, the effective viscosity is lower along the director and more fluid is pumped in the direction of the force, meaning that less fluid is pumped from the perpendicular directions, thus reducing the curvature of the velocity field.If the force is aligned perpendicular to the director (Figure 2b To generalise, our results demonstrate that for moderate anisotropy of the stress tensor that was chosen for graphical representation, flow profiles show obviously visible dependence on the direction of the applied force.We show three distinct cases, however, calculated results for the Green function allow for easy calculation of flow and pressure profiles at any angle between director and the driving force.The value of d 0 can be chosen arbitrarily as long it is much larger than a swimmer size, since it only rescales the velocity magnitude.The values are compared to the isotropic case (dashed lines, see also Figure 1).Flow field is drawn for arbitrary values of the length d 0 and force F. Results show that for the given nematic anisotropy of the viscosity tensor, spreading of the momentum in the direction perpendicular to the director is suppressed, while the direction along the director offers much less resistance to the fluid flow.Note that graphs (d,e) are symmetrical with respect to θ = 0 case (or to φ = 0 case).In (f), which is no longer the case-the velocity field is tilted with respect to the direction of the applied force.

Flow Field of Force Dipole
Microswimmers with no external force applied have no Stokeslet contribution to the surrounding flow field.Instead, flow field is generated by higher moments of the distribution of forces exerted on the fluid, with force dipole being often the leading term.Indeed, typically the flow field of a force dipole exerted on the fluid falls as 1/r 2 , where r is the distance from the swimmer.In addition, since higher moments of force distribution decay faster, dipolar flow field is the most pronounced in the far field of the velocity distribution, which was observed both in self-driven colloidal particles [5] and biological microswimmers [33].
More formally, the flow multipoles are introduced by using the linearity of the Stokes equation and writing the flow field of an arbitrary force distribution in the integral form as Equation ( 34) is then Taylor expanded for distances much larger than the size of a microswimmers, which reveals separate contributions of force monopoles, dipoles, quadrupoles and other multipoles to the flow field of a swimmer, where the force dipole contribution equals [42] and D jk = − r k f j d 3 r is the dipolar force moment.Typically, D ij is decomposed into a traceless symmetrical tensor S ij (stresslet) and an antisymmetrical tensor T ij (rotlet): The trace of D ij was subtracted since the zero compressibility condition renders it ineffective to the fluid flow.The stresslet can be related to swimmer-imposed strain of the fluid and the rotlet to the net torque of the force distribution, with notably both mechanisms relevant in the complex swimming strokes of various micro-elements of micro-organisms.

Stresslet Flow Field in Nematics
To demonstrate stresslet flow field in a nematic, we plot the corresponding flow profiles of the far flow field due to two equal and opposite uniaxial point forces with size F acting on the fluid at a separation l (see scheme in Figure 1).Eigenvalues of the S tensor are in that case −S/2, −S/2, and S, where S = 2 3 lF.Negative S is associated with the pusher type of microswimmers, and positive values with puller type.In Figure 3, we plot flow fields of a puller microswimmer at different orientations of the axis of the applied forces with respect to the nematic director.Velocity field of a pusher has exactly the opposite direction of the flow.In order to calculate the appropriate flow field, one must first determine the components of the S tensor.For example, Figure 3 shows a flow field due to two contractile forces along z-direction distance l apart.In that case, the components of the S tensor equal: S zz = S = 2 3 lF, S xx = S yy = −S/2 = −lF/3 and other components equal zero.The stresslet flow field is given by equation where appropriate derivatives of the Green functions (Equations ( 24)-( 28)) are easily calculated.Flow magnitude decays distinctly as 1/r 2 ; therefore, all the main information about the flow field is contained in the angular dependence of the flow components, which we show in the bottom row of Figure 3 and compare them with the solution for the isotropic fluid, which has only a radial flow profile.
If the force dipole is aligned along the director (Figure 3a,d), the flow field retains the radial-only flow profile and only the magnitude of the velocity field is changed: here given in spherical coordinates v( r) = v r e r + v θ e θ + v φ e φ , where the z-axis is defined by homogeneous director field.Radial-only flow profile in the aligned case is a necessary consequence of the symmetry imposed by the direction of the director and the force dipole, 1/r 2 dependence of the flow, and flow incompressibility.However, in the case where the force dipole is perpendicular to the director (Figure 3b,c,f,g), flow field attains polar and azimuthal components.The profiles of the velocity magnitudes once again show rough correspondence with the regime of the lowest fluid resistance along the director.Upon tilting the force dipole at angle 45 • to the director (Figure 3d,h), the symmetry of the velocity profile is reduced and the difference to the radial profile is even more pronounced.We have measured the volume V t around a dipolar swimmer, where the flow magnitude exceeds the value of (as some selected value).For a microswimmer aligned along the director, V t = 0.36d 3 0 .If the force dipole is perpendicular to the director V t is reduced by 4.2%.We expect similar behavior if one would consider the detection volume of a microswimmer [52]-a volume within which the gradients of the velocity field exceed a certain threshold.In nematic fluids, at the same force dipole, a microswimmer has lesser influence on the surrounding flow field if it is oriented at large angles with respect to the director.(a,e) flow field of force dipole aligned with the director has only radial component.Note that radial flow field is characteristic for dipolar flow in isotropic fluids; however, in nematic fluids, there is still a difference between the magnitudes of the flow field between the isotropic case and a dipole aligned along the director due to anisotropic viscosity (see, e).In (a,e), the magnitude of the velocity field falls to zero at angle θ = 48.3• , whereas, in the isotropic case, the transition from inward to outward flow occurs at θ = 54.7 • .At (b,c,f,g) 90 • or (d,h) 45 • , angles between force dipole and the director before radial flow configuration gains additional terms in the azimuthal and polar directions.

Rotlet Flow
Motile parts of a microswimmer can in principle impose torques on the surrounding fluid, for example by rotating flagella [53].However, for an isolated microswimmer, total torque applied on the fluid has to vanish [1].To show the effect of torques in a nematic, we plot a flow field due to a point torque in Figure 4.A point torque in the y-direction was constructed by an antisymmetric force dipole tensor T ij , where T xz = −T zx = T and other components equal zero.Flow field follows from Equation ( 35): where derivatives of the Green functions (Equations ( 24)-( 28)) need to be calculated.Compared to the isotropic fluid rotlet flows (Figure 1), Figure 4 shows concentric vortices around the origin deformed due to the anisotropy of the nematic fluid.The deformation occurs both in terms of magnitude of the velocity field, which is stretched along the director, and in terms of the shape of the velocity profile, which gains a radial component (compared to the isotropic case of only polar flow direction).

Flow Fields of Sources and Sinks in Homogeneous Nematics
Flow sources and sinks represent another type of fundamental microfluidic elements that in various arrangements could also be associated with micro-organisms.More formally, point sources and sinks are another set of elementary solutions to the Stokes and compressibility equations [42].They are characterized by zero force exerted on the fluid, but a non-zero velocity divergence at the position of the source or sink (in our case at r = 0): where A represents volumetric flow rate through a source (positive A) or a sink (negative A).
Point-source dipoles or even more complex flow structures can be achieved using a linear combination of flow sources and sinks, as allowed by the linearity of the Stokes equation.We calculate the velocity and pressure profiles of a point source or sink in the homogeneous nematic and find the exact solutions in the reciprocal space and write the expressions in the real space in the limit of weak viscous anisotropy, using the same procedure as in Section 2.1.Notably, we obtain the velocity profile directly as analytic formulas.We write in terms of x-, y-, and z-components the compressibility condition (Equation ( 40)) and the Stokes equation for homogeneous director along the z-axis (Equations ( 4)-( 6) for zero force) and perform the Fourier transform, similarly as was done in Equations ( 10)- (13).We obtain the solutions for velocity and pressure in the k-space: We expand over and perform the inverse Fourier transform.The results for the velocity field v( r) = v x e x + v y e y + v z e z in the spherical coordinates v( r) = v r e r + v θ e θ + v φ e φ are: The velocity field of a point source is drawn in Figure 5.As shown in Equations ( 45) and ( 46), flow field has only a radial component.Similar to the stresslet flow in Figure 3a,d, the purely radial flow field is a consequence of the symmetry imposed by the director, flow incompressibility outside the source, and 1/r 2 dependence of the flow.However, the anisotropy of the nematic medium results in the azimuthal dependence of the flow magnitude.The flow profile is dependent only on the α 1 /α 4 and α 6 /α 4 ratios of Leslie viscosities, whereas differently the pressure field is dependent on four viscosity coefficients (α 1 , α 4 , α 5 , α 6 ).In the isotropic fluids (α i =4 = 0), the pressure is a delta function at the point source.Other viscosity terms add 1/r 3 terms that also depend on the azimuthal angle θ.  40)) shown for an isotropic and for a nematic fluid as a function of the azimuthal angle.In nematic fluid, the velocity field retains the radial direction; however, its amplitude as a function of azimuthal angle θ shows an along the z-axis-the director-and a decrease perpendicular to the director.Note that the flow field of a point sink is obtained directly by taking the opposite sign of the flow field of the source; (d) flow source and sink can be combined in a source dipole, here shown for isotropic medium.

Source Dipole Flow
Another elementary flow profile is a source dipole, which was shown to be an important contribution to the flow field of microswimmers [5].It is formulated as a flow due to a point source and a point sink (each with strength A), which are separated by a dipole vector d.The solution is obtained by using solutions for flow source and flow sink (as calculated In Section 2.4) and Taylor expanding the formulas for distances r much larger than d, giving where B = 1 + 3 4 α1 + α6 , C = − 9 2 α1 − 3α 6 , and D = 15 4 α1 .Note that since the curl of the source flow is not zero-contrary to the isotropic case-we are not allowed to use the formalism of the potential flow [42] and have derived the above equation directly from Equation (45).For α1 = 0 and α6 = 0, Equation ( 48) is equivalent to the solution for isotropic fluid [42] that is shown in Figure 5d.For non-zero values of α1 and α6 , source dipole flow is dependent on the orientation between the dipole and the director.Figure 6 shows flow profiles for (a,d) dipole parallel to the director, (b,c,f,g) dipole perpendicular to the director, and (d,h) dipole at angle 45 • to the director.Similar to other flow profiles discussed in this article, the magnitude of the velocity field is spread along the director field.In particular, at angle 45 • between the director and the dipole, flow magnitude profile is skewed compared to the isotropic case.The bottom row of Figure 6 gives a direct comparison to the isotropic solution in terms of radial and azimuthal (or polar) component of the velocity field.

Assumption of Weakly Anisotropic Nematic Fluid
The demonstrated results for the Green function in nematics rely on the first order expansion in the ratio α i =4 /α 4 , where α i are six Leslie viscosities (coefficients), from which α 4 corresponds to isotropic stress.For the variety of nematic materials-i.e., the variety of regimes of values of the Leslie coefficients-this assumption is not necessarily justifiable in general.Therefore, we calculate in full-by numerical integration of the inverse Fourier transforms of the Green functions in Equations ( 10)-( 13)-the corresponding flow profiles and compare them to the analytical flow profiles obtained through first order approximation of small α i =4 /α 4 .Note that in order to achieve good numerical convergence of the Fourier integrals, we have followed the approach in Ref. [31] and multiplied Equations ( 10)-( 13) by a regularization function e −x 2 0 k 2 /25π to eliminate large wavenumbers.We compare angular dependence of a Stokeslet flow for analytical and numerical results at the values of α 1 , η b , and η c used throughout this article and at higher/lower anisotropies.The comparison is shown in Figure 7.At smaller anisotropy (α 1 = −0.075,ν b = 0.3, ν c = −0.15), the full numerical calculations show an excellent match with first order approximations.At higher anisotropy (α 1 = −0.3,ν b = 1.2, ν c = −0.6),there is a clear mismatch between numerical and analytical results, indicating a clear overstep of the first-order approximation in α i =4 /α 4 .Analytical and numerical results for the viscosity coefficients used in the figures of this article (α 1 = −0.15,ν b = 0.6, ν c = −0.3)show some discrepancies, but they still capture similar behaviour.To generalise, the analytic to numeric comparison shows that, although slight discrepancies between the exact flow profiles are possible, flow profiles calculated in this article capture the essential response of flow fields to the anisotropy in the environment.24)-( 29) are compared to the results obtained by numerical inverse Fourier transform of Equations ( 10)-( 13) at distance d 0 from the point of force origin.The comparison is performed for the values of Leslie coefficients used throughout this article (middle column), values twice as small (first column) and twice as large (third column).

Deformations in the Director Profile
We have calculated the nematic Green functions by also assuming a homogeneous nematic director field (and also constant nematic degree of order).However, director deformations may occur due to anchoring of the nematic molecules on the surface of a microswimmer, anchoring at the confining interfaces, or due to backflow coupling to high velocity gradients.In a more complete analysis, simultaneously with the Stokes equation, an equation for the nematic order would have to be solved-for example the Ericksen-Leslie equation for the time derivative of the director alignment (Equation ( 3)).Due to complexity of the problem, the solution would most likely be only numerical.The resulting behavior could in principle include oscillatory solutions or even more complex temporal flows.However, experiments [19,20] and simulations [54] that consider interconnected behavior of nematic order and microswimmer-induced flow profiles typically observe stable trajectories of microswimmer motion and orientation.More complex temporal behavior though can occur at larger concentrations of microswimmers [20].Notice that in real systems-such as swimmers-the centres of these flow multipoles are effectively within the body of the swimmer.Consequently, the exact singularity is always overlayed physically with the swimmer and only the flow fields already at some distance from the singularity affect the nematic.Additionally, surface anchoring at the body of the microfluidic object or possible external fields can further stabilize the nematodynamic behavior.Anchoring induced director field deformations due to a single swimmer and elastic interactions between individual swimmers, mediated by the deformation of the director field, can be assessed within the multipole expansion of the director field components [55].In that case, the corrections to the homogeneous director field typically decay with a square of distance from the swimmer (in a case of a dipolar configuration) or faster.The torque on the nematic director due to material flow decays with the gradient of the velocity field, which in the leading (isotropic) approximation scales as r 3 for a Stokeslet and stresslet, respectively.Director field deformations due to all of the above effects enter the divergence of the stress tensor and alter the velocity profile in the close proximity of a microswimmer.However, if at a sufficient distance from a microswimmer a homogeneous director far field is established, a general set of solutions for the flow field coincides with the results obtained in this article.A homogeneous director field could be additionally imposed by strong electric or magnetic fields.Finally, introducing and exploring the full role of nematic elasticity in the response of microswimmers is a highly exciting and interesting challenge.

Possible Application to Experiments
Recently, thorough investigative work has been conducted in the field of bacteria swimming in lyotropic nematics [19,20].Due to anchoring of the nematic molecules on the surface of bacteria in lyotropic nematics, bacteria tend to align along the nematic director [18,19].This is one of the rationale behind the control of the bacteria swimming direction by using patterned surface anchoring [56], similar to patterned anchoring in colloidal suspensions driven by liquid-crystal-enabled electrophoresis [57].However, in principle, bacteria can change its direction and swim even perpendicular to the nematic director [58].This means that, in principle, all orientations of force dipoles relative to the director field (as shown in Figure 3) are relevant to experimental applications.Ref. [18] shows how bacteria swimming along the nematic director can advect particles in the direction of its way, while particles which are positioned on the side of the swimming trajectory of the bacteria move only a little.This could be explained by the fact that the stresslet velocity field due to bacteria is strong mostly along the director, as shown in Figure 3 and is further affected by the strong anisotropy of the viscosity of the lyotropic nematic used in Ref. [18].In Ref. [19], a collision of two bacteria swimming along the nematic director is discussed, leading to formation of chains.Notice though that if bacteria could be oriented at a notable angle relative to the director, such collision and/or binding events could further gain in complexity.As seen from Figure 3, non-radial terms in the flow field could in principle shift the colliding bacteria to non-colliding trajectories.
Emergent velocity profiles as discussed in this article are highly dependent on the microswimmer orientation with respect to the nematic director.This suggests that the swimming speed and the power consumption will also depend on the microswimmer orientation.Indeed, the existence of an effectively lower viscosity direction for swimming provides another or additional mechanism that affects the direction of swimming of the micro-swimmers [18,19].Simulations for spherical squirmers show that pusher types tend to swim along the director, while pullers tend to swim perpendicular to the director [54].It would be interesting to see how microswimmers would behave if the preferred alignment due to elasticity effects would be at large angles compared to the effective lower viscosity direction of swimming.Possibly, this could be achieved by different types of microswimmers, more complex shapes or anchoring conditions at the surface of a microswimmer, or possibly by using nematics with different material parameters.
The mutual interaction of bacteria depend on their concentration and the strength of force dipoles they impose upon the nematic.At high concentrations, interactions between bacteria are mediated through the deformation of the director field [18].Instabilities due to the coupling of the orientation of bacteria, director field of the nematic, and bacterial flow can induce periodic stripe patterns in the director field [20].However, at more dilute suspensions of bacteria, deformation of the director field might get suppressed by the homogeneous anchoring at the confining walls.In such cases, hydrodynamic interactions between bacteria could become of prime importance, including the effect of anisotropic viscosity as investigated in this article, suggesting that nematic flow profiles that mediate hydrodynamic interaction between microswimmers could significantly alter collective behavior of active particle suspensions.

Conclusions
Elementary solutions of the Stokes equation in nematic fluids are important from the perspective of understanding nematic flows as well as in the context of swimming active particles in anisotropic fluids.We have derived Green functions for the Stokes equation for homogeneous nematic fluid (given by Equations ( 24)-( 32)) under the approximations that (isotropic) viscosity α 4 is the leading viscosity coefficient in the Ericksen-Leslie stress tensor.Specifically, to obtain the solution, we have written Green functions for the Stokes equation for velocity and pressure in the Fourier space, expanded them to first order in anisotropic viscosity coefficients, and converted them back to real space.The calculated Green functions or their derivatives were then used to determine the flow field profiles of point forces (Equation ( 33)), stresslet (Equation ( 37)) and point torques (Equation ( 39)).Similarly, we have derived flow fields of point sources (sinks)(Equations ( 45)-( 46)) and source dipoles (Equation ( 48)) in nematics.The solutions for the corresponding flow profiles are given as direct analytical formulas that could be used in further studies.The general qualitative consequence of anisotropic viscosity is that flow profiles change depending on the angle between the applied force (or stresslet or rotlet) and the nematic director.For viscosity coefficients used in the graphical representation of the derived equations, effective viscosity is smallest if flow direction coincides with the director.Therefore, for Stokeslet flow, the magnitude of the velocity field is the highest if the point force is aligned along the director (Figure 2a).If the alignment is at some angle, we observe a tilt of the velocity field towards the director (Figure 2c,f).The main characteristic of the stresslet flow in nematics is that, unless the force dipole is aligned along the director, flow field gains non-radial terms (Figure 3).Similarly, concentric vortices that form a rotlet flow in isotropic fluids get stretched along the nematic director in nematics (Figure 4).Source (sink) flow field remains only radial also in nematics; however, the flow magnitude shows an oval shape, oriented along the director (Figure 5).The obtained results were also critically assessed with regards to the two main assumptions of our approach-the leading isotropic viscosity and homogeneous nematic director-identifying the regimes of validity of the approach.In conclusion, our work provides an analytical insight into the field of elementary flow elements in nematic fluids-relevant for e.g., microswimmers in experiments in nematic liquid crystals or in biological fluids with anisotropic order-where the derived Green functions and flow profiles can be used to further construct flow fields of other-in principle arbitrary-swimmers and model their dynamics.

Figure 1 .
Figure 1.(a) an outline of the problem-a microswimmer exerts forces upon the surrounding nematic fluid with homogeneous director n, driving a flow field v( r).A set of elementary flow profile solutions for isotropic fluids is shown as reference for later comparison: (b) Stokeslet flow due to a point force; (c) stresslet flow due to a pair of opposite forces (a force dipole), and (d) rotlet flow due to a point torque.The magnitude of the flow field as noted in the colorbar is given in basic quantities (see text for more).
,e), the flow profile shows only a mirror symmetry with respect to xy, xz, and yz planes.The spreading of the flow perpendicular to the force is more efficient along the director, showing larger magnitudes in the xz profile (Figure 2b,f opposed to xy profile in Figure 2b,e.At some angle between the force and the director (45 • in Figure 2c,f), the flow field retains only a symmetry with respect to the r → − r transformation.The flow magnitude is generally larger in the direction along the director.

Figure 2 .
Figure 2. Flow field of point force in nematic fluid oriented (a,d) parallel, (b,e) perpendicular, or (c,f) at angle 45• to the director.The flow field decays as 1/r with the distance from the point force.In the bottom row, the velocity field component parallel to the force v (along e ) and perpendicular to the force v ⊥ (along e ⊥ ) at distance d 0 from the centre is shown as a function of azimuthal and polar angle.The value of d 0 can be chosen arbitrarily as long it is much larger than a swimmer size, since it only rescales the velocity magnitude.The values are compared to the isotropic case (dashed lines, see also Figure1).Flow field is drawn for arbitrary values of the length d 0 and force F. Results show that for the given nematic anisotropy of the viscosity tensor, spreading of the momentum in the direction perpendicular to the director is suppressed, while the direction along the director offers much less resistance to the fluid flow.Note that graphs (d,e) are symmetrical with respect to θ = 0 case (or to

Figure 3 .
Figure 3. Flow field of a force dipole at different orientations with respect to the director field.The bottom row shows radial and azimuthal (or polar) component of the velocity at distance d 0 from the centre as a function of azimuthal (or polar) angle compared to the result for isotropic fluids.(a,e)flow field of force dipole aligned with the director has only radial component.Note that radial flow field is characteristic for dipolar flow in isotropic fluids; however, in nematic fluids, there is still a difference between the magnitudes of the flow field between the isotropic case and a dipole aligned along the director due to anisotropic viscosity (see, e).In (a,e), the magnitude of the velocity field falls to zero at angle θ = 48.3• , whereas, in the isotropic case, the transition from inward to outward flow occurs at θ = 54.7 • .At (b,c,f,g) 90 • or (d,h) 45 • , angles between force dipole and the director before radial flow configuration gains additional terms in the azimuthal and polar directions.

Figure 4 .
Figure 4. Flow field of a point torque in nematic.(a,b) flow field due to torque that is perpendicular to the director.Compared to the isotropic case (b, see also Figure 1), the vortex stretches in the direction of the director; (c,d) same flow field in the xy cross-section.In (c), flow lines point in (out) of the plane.

Figure 5 .
Figure 5. Flow field of point source in (a) isotropic and (b) nematic fluid; and (c) radial flow velocity for a source flow (Equation (40)) shown for an isotropic and for a nematic fluid as a function of the azimuthal angle.In nematic fluid, the velocity field retains the radial direction; however, its amplitude as a function of azimuthal angle θ shows an along the z-axis-the director-and a decrease perpendicular to the director.Note that the flow field of a point sink is obtained directly by taking the opposite sign of the flow field of the source; (d) flow source and sink can be combined in a source dipole, here shown for isotropic medium.

Figure 6 .
Figure 6.Flow field of source dipole in nematic (Equation (48)) for the dipole aligned (a) parallel to the director; (b,c) perpendicular to the director, and (d) at angle 45 • to the director.The bottom row shows radial and azimuthal (or polar) components of the velocity field compared to the solution for isotropic fluid (α 1 = α 6 = 0) for each of the cases above.Spreading of the velocity magnitude along the director is observed.

Figure 7 .
Figure 7.Evaluation of small α i /α 4 expansion assumption.Results for the point force along the director (first row) and perpendicular to the director (second row) obtained through the nematic Green function in Equations (24)-(29) are compared to the results obtained by numerical inverse Fourier transform of Equations (10)-(13) at distance d 0 from the point of force origin.The comparison is performed for the values of Leslie coefficients used throughout this article (middle column), values twice as small (first column) and twice as large (third column).