Q-tensor model for undulatory swimming in lyotropic liquid crystal polymers

Abstract Microorganisms may exhibit rich swimming behaviours in anisotropic fluids, such as liquid crystals, which have direction-dependent physical and rheological properties. Here we construct a two-dimensional computation model to study the undulatory swimming mechanisms of microswimmers in a solution of rigid, rodlike liquid crystal polymers. We describe the fluid phase using Doi's $Q$-tensor model, and treat the swimmer as a finite-length flexible fibre with imposed propagating travelling waves on the body curvature. The fluid–structure interactions are resolved via an immersed boundary method. Compared with the swimming dynamics in Newtonian fluids, we observe non-Newtonian behaviours that feature both enhanced and retarded swimming motions in lyotropic liquid crystal polymers. We reveal the propulsion mechanism by analysing the near-body flow fields and polymeric force distributions, together with asymptotic analysis for an idealized model of Taylor's swimming sheet.

beating to generate directional motions by breaking symmetry (Lauga & Powers 2009). Additionally, extensive research has been conducted to reveal rich dynamics and new propulsion mechanisms that use the complex fluids' intrinsic anisotropy arising from the non-Newtonian physical and rheological properties. For example, while the extra particle stresses from shear-thinning viscoelastic (e.g. Oldroyd-B) fluids generally impose additional hindrance effects, it has been found that undulatory swimmers, such as C. elegans, may leverage both the finite body-length and the polymeric stress relaxation to achieve a higher swimming speed than that in Newtonian fluids (Teran, Fauci & Shelley 2010; Thomases & Guy 2014;Salazar, Roma & Ceniceros 2016). Furthermore, there has been a growing interest in exploring non-equilibrium physics of biosynthetic active materials 'powered' by many-body interactions in complex fluids that are capable of exploiting collective dynamics for useful mechanical work (Ramaswamy 2010;Shelley 2016). Of particular interest is the study of active microparticles (e.g. motile bacteria) in lyotropic liquid crystals (LCs) wherein the extra stress generation is determined by the LC's orientational order, which leads to intriguing swimming behaviours of microswimmers, as well as far-from-equilibrium physical and topological properties of the LCs (Zhou et al. 2014;Lavrentovich 2016 So far, there have been several models proposed to study the dynamics of microswimmers in lyotropic LCs. Zhou et al. (2017) employed a Q-tensor model (here Q-tensor denotes the second-rank orientational order-parameter tensor) derived from the generalized Ericksen-Leslie equation (Sonnet, Maffettone & Virga 2004) to solve for the orientation field of LCs when a rigid, rodlike particle (B. subtilis) moves in a homeotropic nematic cell where the director is perpendicular to the cell wall. They illustrated how the induced shear determines the ordering patterns (or birefringent cloud) around the moving particle. Using an Ericksen-Leslie model described by the nematic director field (Larson 1999) and enforcing the local director orientation to be tangential to the undulatory body (i.e. tangential anchoring), Krieger, Dias & Powers (2015a); Krieger, Spagnolie & Powers (2015b) studied the motion of Taylor's swimming sheet (Taylor 1951), an infinite-length, zero-thickness sheet with imposed small-amplitude travelling waves, in unconfined LCs. The asymptotic solutions they derived suggest that both forward and backward motions can be achieved in different parameter regimes, and the swimmer's mean speed can be either faster or slower than that in a Newtonian fluid. Similar bi-directional swimming motions were reported by the same group when studying the motions of Taylor's swimming sheet in LCs confined by two parallel plates under the tangential anchoring condition (Krieger, Spagnolie & Powers 2019). In addition to asymptotic analysis, they performed nonlinear simulations using the immersed boundary (IB) method to resolve the fluid-structure interactions of an infinite-length sheet with prescribed finite-amplitude travelling waves. Nevertheless, as pointed out by Krieger et al. (2019), the Ericksen-Leslie model cannot be robustly extended to more general cases where the anchoring direction is misaligned with the director (e.g. homeotropic anchoring), and may break down owing to generations of singularities in the director field.
Motivated by the prior studies and to further seek quantitative understanding of swimming mechanisms in anisotropic fluids, we present a computational framework for studying the undulatory motion of a finite-length microswimmer in a solution of liquid crystal polymers (LCPs), a class of rigid, rodlike aromatic polymers that have much larger sizes and higher aspect ratios than small LC molecules (e.g. para-azoxyanisole). Compared with flexible polymers or molecule aggregates (e.g. lyotropic chromonic LCs Zhou et al. 2014;Zhou 2018), LCPs can much more easily re-orient themselves when interacting with external fields (e.g. fluid flows, electric fields) to show large birefringence (Doi & Edwards 1988;Larson 1999). As discussed below, the evolution of their orientation distributions can be described by a Q-tensor model that is coarse-grained from Doi's kinetic model (Doi 1981;Doi & Edwards 1988), a prevalent theoretical model for simulating the flow and rheology of LCPs. We treat the microswimmer as a finite-length elastic fibre whose undulatory motion is activated by imposing travelling waves on the body curvature. An IB method is employed to resolve the fluid-structure interactions for the swimmer. The nonlinear simulations of finite-amplitude swimming motions reveal both enhanced and retarded swimming motions in the nematic regime, which correspond to the scenarios when the swimming direction is parallel and perpendicular to the nematic alignment directions, respectively. Our numerical results are qualitatively confirmed by the asymptotic solutions of Taylor's swimming sheet model, and can be further understood by revealing the characteristic near-body flow fields and polymer force distributions.
The paper is organized as follows. Section 2 is dedicated to the mathematical formulation of the problem. We introduce the governing equations for a minimal Q-tensor model, and construct an IB formulation. In § 3, we systematically explore the parameter space to study and analyse the finite-amplitude undulatory swimming motions via nonlinear simulations and analytical asymptotic analysis. Finally, we summarize and draw conclusions in § 4. The details of the IB numerical scheme and the derivation steps of the asymptotic solutions of Taylor's swimming sheet are provided in Appendices A and B, respectively.

Mathematical model
Doi's kinetic model (Doi & Edwards 1988) adopts a probability distribution function (p.d.f.) Ψ (x, p, t) to describe the ensemble dynamic of concentrated microrods in terms of their centre-of-mass position x and orientation p (|p| = 1). When simulating the macro scale hydrodynamic interactions between rods and ambient flows, it is often preferred to use coarse-grained partial differential equations that are derived via proper moment closure owing to their low computation costs by removing p dependency. Consider rodlike polymers with fore-aft symmetry and uniform length . In the fluid domain Ω f , a two-dimensional (2-D) model can be derived for the second-moment tensor D(x, t) = p ppΨ dp (Doi 1981;Doi & Edwards 1988;Feng & Leal 1997) as is the so-called upper-convected time derivative, E = (∇u + ∇u T )/2 is the symmetric strain-rate tensor and S(x, t) = p ppppΨ dp is the fourth-moment tensor. The first term on the right-hand side represents the steric alignment effects arising from the Maier-Saupe potential written as with ζ as the dimensionless strength coefficient (Maier & Saupe 1958). The diffusivity coefficients d t and d r originate from the translational and rotational Brownian diffusion motion, respectively. The above equation is a 'Q-tensor' model because the trace-free (normalized) order-parameter tensor, the so-called Q-tensor, is defined by Q(x, t) = c −1 D − I/2, with c(x, t) = 1 as the concentration field (i.e. zeroth-moment of Ψ ) that remains as a homogeneous constant. In two dimensions, the tensor Q has a maximal non-negative eigenvalue ρ satisfying 0 ≤ ρ ≤ 1/2. Assuming that ρ is isolated, then we call its associated unit eigenvector the director m, and 0 ≤ S = 2ρ ≤ 1 the orientational order parameter. The incompressible fluid velocity u is incorporated in a mean-field fashion, and is governed by the forced Stokes equation: where μ f is the fluid viscosity and f e (x, t) represents the force exerted by the elastic swimmer. The extra stress of the polymer solution is defined as (Doi 1981;Feng & Leal 1997) where ν ∼ −3 is an effective volume fraction, k B T is the thermal energy scale, β 0 (ν 3 ) −2 represents a crowdedness factor with β 0 being an empirical coefficient. In this study, we will choose a small empirical crowdedness factor β 0 (ν 3 ) −2 ∼ 10 −3 − 10 −4 (Feng & Leal 1997;Feng, Chaubal & Leal 1998).
Here we have several comments on the above Q-tensor model. First, we consider the model to be minimal because (2.2) only describes the nematic elasticity arising from the rods' orientation distribution but ignores spatial heterogeneity. According to Greco & Marrucci (1992), the Maier-Saupe potential can be further improved by taking the distortion elasticity into account by adding a term proportional to pp : D. However, here we consider the scenarios in the nematic regime where the microswimmer undulations only weakly disturb the LCPs, and the concentration field remains uniform across the domain. For simplicity, we adopt Doi's original model with the Maier-Saupe potential in (2.2). This is consistent with the nonlinear simulation results shown below where the polymer distribution remains spatially homogeneous, with small near-body fluctuations in the orientational order.
Second, as pointed out by Feng, Sgalari & Leal (2000), adding distortion elasticity into the above Q-tensor model can mathematically recover the director formulation of the Ericksen-Leslie model in the limit of weak flow and mild spatial distortion, which can be described by the so-called Frank elasticity with equal (Frank) constants (DeGennes & Prost 1993). Nevertheless, the director formulation is more suitable for modelling small-molecule LCs that have short orientation relaxation time and hence admit a uniaxial form for the quasi-equilibrium orientational distributions. In comparison, the LCPs' orientations are easily distorted into a non-uniaxial configuration, especially when imposing various flow conditions (e.g. shear and extensional flows). Therefore, the generalized LC theories, either of Ericksen-Leslie or Doi type, which adopt a tensorial order parameter are preferred in simulating LCPs (Beris & Edwards 1994;Qian & Sheng 1998;Feng et al. 2000;Sonnet et al. 2004;Klein et al. 2007).
Third, our Q-tensor model is essentially non-polar, which has a two-fold meaning. At the micro level, the Doi-Onsager kinetic model adopts Jeffrey's equation to describe the polymers' rotational motion driven by the mean-field fluid velocity gradient and the local alignment torque arising from the Maier-Saupe potential. The resultant rotational flux term hence takes the following form: where ∇ p = (I − pp) · ∂/∂p denotes the gradient operator on the unit sphere of orientations. One can prove there is no net-torque exerting on the unit volume of the fluid phase by taking the configurational average of the torque produced by a single test polymer, i.e. p ∇ p U MS Ψ dp = 0 (Feng et al. 2000). At the macro level, the symmetric extra stress in (2.5) guarantees the conservation of angular momentum in the fluid phase. Moreover, when interior or exterior boundaries are introduced, only the boundary conditions (e.g. no-slip condition) for the velocity field need to be implemented. Solving the coarse-grained evolution equation (2.1) does not require imposing any boundary conditions on D in the limit d t → 0. Instead, the near-boundary LCPs' orientational distribution is determined by the interplay among the resultant fluid shear, mean-field alignment torque and rotational diffusion as reflected by the right-hand side of (2.6) microscopically.
In the Lagrangian frame (Ω L ) for the undulatory swimmer of length L s , its kinematics can be described by the parametric form X (s, t), with s the local arc length s ∈ [0, L s ]. The actual shape X (s, t) is determined by penalizing deviations from actuation imposed on the body curvature κ 0 (s, t) as starting from an initial configuration given by X (s, t = 0) = (s, A 0 (sin(ks))). Equation (2.7) describes the rightward-propagating travelling waves with amplitude A 0 , wavenumber k and angular frequency ω. Imposing actuation in (2.7) produces elastic forces F e (X ) along the body, and effectively drive periodic shape changes (or swimming gaits). Following Peskin (2002), the Lagrangian body force can be derived by performing the variational derivative upon the elastic energy E, i.e. (2.8) The discretized forms for the two components of F e are given in (A2). Here the total elastic energy E[X ] include the contributions from both stretching (denoted by subscript s) and bending (denoted by subscript b) deformation (Fauci & Peskin 1988): where n denotes the local normal direction. Using the IB method, the Eulerian forcing term f e can be calculated by where δ denotes the Dirac delta function that maps between the Eulerian (Ω f ) and Lagrangian (Ω L ) domain. Likewise, the Lagrangian velocity can be interpolated from the nearby Eulerian grids via Here we perform mapping between the Eulerian (x(x, y)) and the Lagrangian (X (X, Y)) domain via a four-point Dirac delta function that is numerically constructed as where h is the (uniform) Eulerian mesh width, and the function φ(r) is constructed as (2.13) This choice of the numerical Delta function guarantees the momentum conservation in both the fluid and the solid phases (Peskin 2002). For non-dimensionalization, we choose the wavelength λ = 2πk −1 as the length scale, wave speed ω/k as the velocity scale and 2νk B T as the LCP's stress scale, which lead to the dimensionless equations in the fluid phase after proper re-scaling: (2.17) Here we define two Péclet numbers, Pe = ω/8πd r and Pe t = 2πω/d t k 2 , which represent the ratio of the time scales for the rod's rotation and transport over that of undulation (i.e. ω −1 ), respectively. Especially Pe characterizes the time evolution of the orientation field. In this study, we vary Pe ∼ O(10 −1 ) − O(10) over three orders of magnitude to fully resolve the non-Newtonian swimming behaviours that are prominent at a finite Pe (100) is chosen to be at least one order of magnitude higher than Pe so that the translational diffusion effect is small or negligible. The coupling between the LCPs and the viscous fluid is characterized by the Ericksen number (Larson 1999;Krieger et al. 2015b). This parameter choice reveals the strong (weak) flow modification on the polymers' configuration at higher (lower) Er values. When Er becomes even higher (e.g. Er = 20), we notice that the numerical solutions may be difficult to converge. Another dimensionless number in the fluid phase is an effective crowdedness coefficient . Equations (2.7)-(2.9) and the two mapping functions (2.10) and (2.11) retain the same mathematical forms, with the two dimensionless stiffness coefficients We fix the wavenumber k = 2π, the swimmer length L = L s k/2π = 1 (same as the wavelength) and the amplitude A = A 0 k/2π = 0.05. We choose a large stretching stiffness σ s = O(10 3 ) to enforce inextensibility and vary the bending stiffness σ b = O(10 −3 ) − O(10 −1 ) to achieve various different swimming gaits. We then use a pseudo-spectral method to solve the above coupled equations in a periodic computation domain. The reader is referred to Appendix A for more details of the numerical algorithm and parameter choice.

Results and discussion
Without the swimmer, the polymer distribution can be described by the equilibrium solution of the original Doi's kinetic model, where the coefficient χ satisfies .
( 3.2) As shown in figure 1(a), χ(ζ ) admits a bifurcation at ζ c = 4, and has two solution branches for the isotropic and nematic states. In figure 1(b), the symmetric form of Ψ 0 (φ) suggests that the rod's primary (or principal) orientation direction is aligned with the x-axis when φ = 0 and π are in the nematic regime (ζ ≥ ζ c ), compared with the constant solution in the isotropic regime (0 ≤ ζ < ζ c ). At each time step in the simulation, we reconstruct Ψ 0 (x, p, t) of form (3.1) locally in the principal coordinates of the D(x, t) field, and then approximate the fourth moment as where we follow a quasi-equilibrium ansatz that D and S co-align in principal directions in the eigenspace (see more details in . This method is also referred to as the Bingham closure (Bingham 1974;Feng et al. 1998), which has proven to be accurate in p.d.f. reconstruction and has been widely used in simulating LCPs (Feng et al. 1998;Gao & Li 2017;).
In the following, we numerically investigate the undulatory motions of a finite-length swimmer in LCPs whose equilibrium configurations are described by Ψ 0 , and focus on the scenarios when the swimmer moves either parallel with or perpendicular to the alignment direction (i.e. x-direction) for both the 'stiff' (σ b = O(10 −1 )) and 'soft' (σ b = O(10 −3 )) cases. In figure 2(a), we show the time sequence of the swimmer's shape-change during Hence, we denote the mean centre-of-mass velocity as −U LCêx withê x a unit basis vector pointing in the positive direction of the x-axis. We measure the velocity magnitude U LC by performing time averaging over five undulation periods after the swimming speed reaches quasi-steady states, and then compare it with U N measured in Newtonian fluids. As shown in figure 2(b), the speed ratio U LC /U N as a function of ζ clearly shows non-Newtonian behaviours, especially in the nematic regime where similar bifurcations occur, which correspond the motions that are parallel (denoted by symbol '//') and perpendicular (denoted by symbol '⊥') to the x-axis. More specifically, we observe enhanced swimming speeds when the swimmer moves parallel in the nematic regime; while retarded motions are seen in the isotropic regime, and in the most nematic regime during perpendicular motions. Additionally, while qualitatively similar behaviours are consistently observed for the stiff and soft cases in the parameter space, it is seen that the stiff swimmer can generally achieve faster mean speeds, which exhibit more significant speed gain and loss. Here it is worthwhile to mention that when the swimming direction is tilted, i.e. with an arbitrary angle between the nematic alignment direction, we observe that (not reported here) the swimmer can perform an entire-body rotation before reaching a steady free swimming. This observation is consistent with the asymptotic solutions of Taylor's swimming sheet in a nematic LC derived by Shi & Powers (2017). It can be explained by the fact that the anisotropic elastic stress will produce a net torque on the plate when there is any misalignment between the nematic director and the swimming direction. Figure 3 shows the variations of U LC /U N with respect to Pe at various values of Er in both isotropic and nematic regimes. Apparently, for all the cases the speed ratio exhibits stronger non-Newtonian behaviours (U LC /U N / = 1) at larger values of Er when fixing Pe. For the stiff swimmer shown in panels (a-c), at a given Er, U LC /U N is seen to be close to 1 in the small Pe regime, which corresponds to a slow undulation (or equivalently, quick rotational diffusion of LCPs). This is because when Pe → 0, the  rotational diffusion dominates over convection, which not only drives the nematic structure towards equilibrium states (i.e. the right-hand side of (2.17) approximately to be 0) but also produces smaller and smaller τ p . At a high Pe, where convection dominates during fast actuation, the polymer configuration is modulated by flow via constraining stress, i.e. E : S, owing to the rod's rigidity, which leads to more complicated behaviours. For a given Er, in the isotropic regime, U LC /U N monotonically decreases with respect to Pe. Nevertheless, in the nematic regime, U LC /U N varies non-monotonically, with the maximum (minimum) occurring at small Pe (∼ O(1)) during the parallel (perpendicular) swimming motions. Furthermore, in panels (b,c), we mark the extrema locations on the U LC /U N − Pe curves using black dashed lines. It is seen that for both parallel and perpendicular cases, the critical values of Pe that correspond to the maximal and minimal speed ratio become smaller and smaller when increasing Er. In panels (df ), we show the corresponding results for the soft swimmer cases, and similar trends as for the stiff cases are observed. However, the extrema locations in the nematic regime would be at much larger values of Pe where the soft swimmer motions become very slow, and hence are not studied here.
To understand the numerical observations, we build an idealized linear model for Taylor's swimming sheet that has an infinite length (Taylor 1951). With the swimmer moving at the speed −U LCêx , its vertical displacement can be described as Then we simplify the above governing equations by employing a stream function ϕ via withê z as the unit basis vector pointing in the out-of-plane direction, which ensures the incompressibility condition ∇ · u = 0. By following the solution procedure of Lauga (2007) and Krieger et al. (2015bKrieger et al. ( , 2019, we impose the no-slip condition on the wavy sheet, and perform asymptotic analyses by expanding all the variables in terms of ε in both the isotropic and nematic regimes. Moreover, we neglect the crowdedness effect (i.e. β = 0) in τ p and the translational diffusion (i.e. Pe −1 t → 0), and adopt different closure methods for S in the isotropic and nematic regimes to facilitate derivation. After some manipulations (see more derivation details in the supplemental material), we find directional motions only occur in the second order, i.e. U LC = U (2) LC ε 2 , and the speed ratio can be calculated as where D 11 ), and the plus (minus) sign corresponds to parallel (perpendicular) swimming motion with respect to the alignment direction. Unlike the possible bi-directional motions predicted by Krieger et al. (2015b), it is straightforward to show that U LC /U N > 0 for the cases of isotropic (3.6) and parallel swimming in the nematic regime (3.7); while for the perpendicular swimming cases, (3.7) suggests that for a given Pe, U LC /U N < 0 only occurs at very large Er values, which are far above the practical parameter range. Hence the asymptotic model predicts universal unidirectional motions that are in the opposite direction of the prescribed travelling wave, which are also consistent with our numerical results of finite-length swimmers. Moreover, as shown in figure 3(g-i), the asymptotic solutions predict the qualitatively similar behaviours of U LC /U N as the numerical results in figure 3(af ). However, (3.7) suggests the maximal and minimal speed ratios all occur at Pe = ζ − 2 for both the parallel and perpendicular cases and are independent of Er. This analytical prediction is different from the numerical results that depend on both Pe and Er, which is likely owing to the finite-length effect of the swimmer (see more discussion below).
In the isotropic regime, the monotonic decay as a function of Pe at a given Er is reminiscent of the behaviours in viscoelastic (e.g. Oldroyd-B) fluids (Lauga 2007;Shen & Arratia 2011). This can be witnessed by introducing an effective polymer viscosity μ p (Doi & Edwards 1988;Feng & Leal 1997) as μ p /μ f = α(S)ErPe, where ErPe defines the effective polymer contribution to the zero-shear-rate viscosity, and the estimated coefficient α characterizes the dependence on the order parameter S. When ζ → 0,  we estimate S = 0 and α(S) = 1, and hence can rewrite (3.6) as which resembles the asymptotic solution for Taylor's swimming sheet in viscoelastic fluids (Lauga 2007). In the nematic regime, from (3.7), apparently the speed gain and loss are seen to be directly determined by the horizontal (D (0) 11 > 1/2) and vertical (D (0) 11 < 1/2) alignment directions. Moreover, U LC /U N → 1 is observed in the limit of Pe → ∞, which suggests diminishing non-Newtonian behaviours for very fast undulations in sharply aligned LCPs. The reader is referred to Appendix B for more details of the derivation.
Next, we take a closer look at the flow patterns to uncover the nonlinear effect arising from the finite length of the swimmer. As shown in figure 4(a-d), we compare the simulation results of the instantaneous flow fields at the end of each period, after the swimming motions reach quasi-steady states. While the isotropic case in panel (b) resembles the Newtonian one in panel (a), the parallel (perpendicular) swimming case in the nematic regime in panel (c) (panel d) has slightly wider (narrower) circulation zones near the two ends. Indeed, when averaging over three to five periods, we consistently find the flow patterns shown in figure 4(e-h) have distinctive features. Compared with the Newtonian case in panel (e), the isotropic flow in panel ( f ) appears to be slightly weakened. More interestingly, the mean flow field in the nematic regime appears to be extensile (panel g) and contractile (panel h), which correspond to the parallel and perpendicular motions, respectively. Figure 5(a-c) reveal the characteristics of the nematic field. Without imposing particular anchoring conditions for D along the wavy body, we find the time-averaged (denoted by ' ') director distribution m approximately follows the isotropic (panel a) and aligned (panels b,c) configurations. However, owing to the finite length of the swimmer, the resultant structures lack left-right symmetry. The disturbances on the nematic field are concentrated in the areas around the body's rear side without propagating further. High-orientational orders suggests enforcement of local alignment owing to strong steric interactions characterized by large values of ζ . When comparing the two scenarios in panels (b) and (c), it appears that parallel motions can impact larger areas. This is because the vertical oscillations during undulation are a lot faster (∼5-10 times) than the resultant horizontal migration, and can lead to stronger momentum exchange along the y-direction. Hence, it is easier for the shear-induced force to re-orient the LCPs perpendicularly, which leaves large low-order zones during parallel motions. Further examination of the time-averaged polymer force, f p = ∇ · τ p in figure 5, reveals more precise propulsion mechanisms. Because the periodic elastic body force has a zero mean, i.e. f e = 0 from (2.15), the near-body fluid flows are driven by f p at a finite Er. As typically shown in panel (d), isotropic LCPs only permit small f p without showing clear correlations with the swimming direction. In comparison, nematic LCPs lead to much stronger, asymmetrical f p in panel (e, f ), with the largest values near the 'tail.' Additionally, they are clearly seen to serve as the driving forces of the extensile and contractile flow patterns in figures 4(g) and 4(h), respectively. More intriguingly, unlike the isotropic cases, the near-body distributions of f p in the nematic regime highly correlate with the swimming direction. In figure 5(e), the projected force vectors at y = 0 all approximately point leftwards to 'push' the swimmer forward; while in panel ( f ), the projected force along the x-direction, while weaker, appears to be overall pointing rightward to 'pull' the swimmer backward.

Conclusion
This work has built a Q-tensor model to study the nonlinear dynamics of undulatory swimming motions in lyotropic LCPs in the limit of vanishing Reynolds numbers. Combining direct numerical simulations and asymptotic analysis, we have revealed both enhanced and retarded swimming behaviours for a flexible swimmer moving in rigid, rodlike polymers. Compared with those top-down macro models of Ericksen-Leslie or Landau-de Gennes (DeGennes & Prost 1993) for small-molecule LCs, ours is derived bottom-up, and has more accurate microscopic origins and fewer computation parameters. Hence it provides a simple and convenient computation framework to model fluid-structure interactions in lyotropic LCs. Our results suggest that in addition to adopting wavy body undulations to break symmetry, it is possible for a finite-length flexible swimmer to make use of the unique near-body fluid dynamics and polymer force generation to adjust the migration speed. We expect to apply the same methodology to generally study a microorganism's motion, both individually and collectively, in anisotropic complex-fluid environments.
Acknowledgements. Z.L. acknowledges the Fundamental Research Funds for the Central Universities of China grant no. 2020QNA4046. T.G. acknowledges the National Science Foundation grant no. 1943759. The anonymous reviewers' are thanked for their constructive comments that helped improve and clarify this manuscript.
Declaration of interests. The authors report no conflict of interest.

Appendix A. Immersed boundary solver
In the flexible fibre, the elastic body force F e (F x , F y ) is calculated explicitly at each time step. Following the IB approaches (Fauci & Peskin 1988;Lee & Wolgemuth 2016), we discretize the variational derivative of the bending energy functional, and derive the following component form for the ith node at the position X i (X i , Y i ): where the curvature is defined as κ i =ê z · (X i+1 − X i ) × (X i − X i−1 )/ s 3 and s = L/N s . In the fluid phase, we first solve the constitutive equations for the velocity u and the second-moment tensor D. Here we use the method of Vaithianathan & Collins (2003) to enforce the symmetric positive definiteness of D whose diagonal components are bounded, i.e. tr(D) = D 11 + D 22 = 1 and 0 ≤ D 11 , D 22 ≤ 1. We apply the Cholesky decomposition on D to obtain D = L · L T = l 11 0 l 21 l 22 where L is a lower triangular matrix. Then we perform another transformation for the components of L, so that the parameters (η, γ, χ) are now unconstrained. Substituting ((A3) and (A4a-c)) into the evolution equation of D leads to K 11 + l 21 l 22 l 11 In the above a shorthand notation K = 2E : S − (ζ /Pe)(D · D − D : S) is introduced and computed explicitly. When D is obtained, we evaluate the polymer stress τ p , and then directly compute the 2-D velocity field u = (u x , u y ) in the Fourier space as where f total = Er∇ · τ p + f e ,ξ is the wavevector with |ξ | = ξ (ξ = ξ /ξ ). The equations in the fluid phase are solved using a pseudo-spectral method, with the time integration being evaluated using the fourth-order Runge-Kutta method. For all simulations performed in this study, we choose the computation parameters as Ω f = 2 × 2, h = 1/128, t = 10 −5 -10 −4 , Pe = 0.1-10, Pe t = 10-100, Er = 1-10, ζ = 0-10, β = 0.0005-0.005, N s = 32, L = 1, A = 0.05, k = 2π, ω = 2π, σ b = 0.005-0.5, σ s = 3000σ b . As shown in figure 6(a-c) for the time-dependent velocity components, we verified the numerical results by performing convergence tests for a parallel moving swimmer.
Moreover, as shown in figure 7, we performed another benchmark study for the undulatory swimming motions in an Oldroyd-B (OB) fluid where the dimensionless Deborah (De) number, which plays a similar role as Pe in the LCP cases, is defined as the wave frequency by the OB fluid relaxation time. We measured the mean centre-of-mass swimming speed U OB of the swimmer that has a long body length L = 4, and compared  the speed ratio with the numerical data of Salazar et al. (2016) and the asymptotic results for Taylor's swimming sheet of Lauga (2007), where η s and η p respectively represent the solvent and polymer contribution to the viscosity. The Newtonian speed U N can be derived as where k = 2π is the wavenumber. We find our results agree well with the previous studies.
In the third test, we examined the impact of varying bending stiffness when choosing high values of σ b . As shown in figure 8 for the time-dependent velocity, only small  differences are seen when σ b goes beyond 0.5. In these scenarios, the swimmer can quickly respond to the imposed target curvature and well follow the travelling-wave actuation. Hence, in the main text, we choose σ b = 0.5 for the cases of a stiff swimmer. the polymer stress, when neglecting β, is written as By ignoring the translational diffusion, the D evolution equation is written as When applying the curl on both sides of (B4), we have Next, we expand all the variables in terms of ε up to the second order as After some manipulations, we can derive the following governing equations. zeroth-order: First-order: where a 0 = (ErPe/(4 + ErPe))16Pe 2 /((4 − ζ ) 2 + 16Pe 2 ). Eventually we can evaluate the mean swimming speed at the second-order as where U N = U (2) N ε 2 = ε 2 /2 is the mean swimming speed in a Newtonian fluid. Note that in Doi's theory, the polymer contribution to the zero-shear-rate viscosity can be effectively defined as (Feng & Leal 1997) μ where S = √ (D : D − 1/2) is the order parameter, α(S) is a concentration coefficient and μ p represents the polymer contribution to the viscosity. Hence, when ζ = 0, we estimate S = 0 and α(S) = 1, which leads to B.2. Nearly-aligned cases In the deep nematic regime, we adopt a classical quadratic closure to approximate S (Doi & Edwards 1988): which corresponds to the strong alignment configuration with large ζ values when ζ > ζ c . At the zeroth-order, we solve for the equilibrium solution D (0) = diag(D The plus (minus) sign in the above equation represents the scenario when the nematic alignment direction is along the x-(y-) axis. Hence, when fixing the swimmer's moving direction to be horizontal (vertical), the above solution corresponds to the parallel (perpendicular) swimming cases. At the first-order, it is straightforward to show that To facilitate further analytical manipulations, we consider the strong alignment cases when ζ is large (ζ 1), which allows us to neglect the last term in (B44), and obtain the same first-order solution as (B27).