Sheared Nematic Liquid Crystal Polymer Monolayers

We provide a comprehensive study on the planar (2D) orienta- tional distributions of nematic polymers under an imposed shear flow of ar- bitrary strength. We extend previous analysis for persistence of equilibria in steady shear and for transitions to unsteady limit cycles, from closure models (21) to the Doi-Hess 2D kinetic equation. A variation on the Boltzmann distri- bution analysis of Constantin et al. (3, 4, 5) and others (8, 22, 23) for potential flow is developed to solve for all persistent steady equilibria, and characterize parameter boundaries where steady states cease to exist, which predicts the transition to tumbling limit cycles.


Introduction.
Macromolecular materials (for example, nematic liquid crystal polymers (LCP)) exhibit large flexibility and can be easily processed into fibers with high strength, or thin films [1,6,26]. The bulk properties of these materials are closely related to the processing flows. So it is very important to understand the dynamic behavior of the materials in the presence of flows.
In the 2-D case, it is well-known that in the absence of flow the isotropic-nematic phase transition occurs at U = 2, where U is the normalized polymer concentration.
In the presence of an imposed weak shear there is a threshold (U 0 ≈ 2.41144646) for U : When U < U 0 , steady state solution exists; otherwise there is no steady state and the orientational pdf is temporally periodic ("tumbling") and can never reach a steady state.
The goal of this work is to conduct a comprehensive study of the 2-D nematic liquid crystal polymers under an arbitrary shear. In [21] monolayer films of liquid crystal polymers have been modeled with a mesoscopic two-dimensional (2D) analogue of the Doi-Hess kinetic model. The weak-shear steady and unsteady selection criteria for 2D nematic polymers using various second-moment closures was derived. A simple proof was given based on the Poincare-Bendixon Theorem to show that limit cycles ("tumbling orbits") exist beyond the parameter boundary for the steady-unsteady transition. Finally, it was revealed that the shear-perturbed 2D phase diagram is significantly robust to closure approximations than the 3D system. Our approach here differs from [21] in that we use the kinetic Smoluchowski equation instead of the Doi tensor model. There are two minor differences between our work and [21]. 1) In [21], a bistable region (where there are two stable solutions) was observed for the Doi closure tensor model. In our study using the full Smoluchowski equation there is no bistable region. Note that the bistable region in [21] is extremely small and it is difficult to pinpoint what small differences between the full Smoluchowski model and closure-based tensor models might cause this small bistable region. 2) In [21], Bogdanov-Takens bifurcation was observed whereas there is no Bogdanov-Takens bifurcation in our study. Again note that in [21] the Bogdanov-Takens bifurcation occurs in a very narrow region in the phase diagram. Our results are in general consistent with those in [21] in several aspects: 1) Both models identify critical values of polymer concentration, beyond which there is no steady state; 2) Both models yield the fold structure in the phase diagram. However, a detailed fold structure is resolved in this work. Specifically, the fold is not a kink, rather it is a smooth turn with large curvature. Furthermore, we find that fold structure in the phase diagram (the plot of order parameter versus U ) exits for P e < 0.746. For P e > 0.746, the fold structure disappears. At the critical value of P e ≈ 0.746, the boundary separating the steady state region and the tumbing region in the P e − U plane appears to have a singularity. In [21] it was mentioned that unstable state can be further classified as tumbling or wagging. Tensor models were originally motivated by the fact that in the Doi-Hess kinetic model with the Maier-Saupe potential, an equilibrium state is completely specified by the orientation tensor. As a result, tensor models can provide good approximations at equilibrium or near equilibrium. In 2-D tensor models the system state is represented by the orientation direction and the order parameter. Thus, it is natural to try to classify time-periodic unsteady state as tumbling or wagging. If the director rotates continuously in one direction, we may say it is tumbling; if the director moves back and forth in a certain range, we may say it is wagging. Even in the tensor models, such a classification may be problematic since the director angle is undefined at the isotropic state. If the time-periodic evolution ever goes through the isotropic state, such a classification may be vague. It is worthwhile to note that "tumbling" refers to the tumbling of the orientation probability density. As for individual polymer rods, under a shear they are always tumbling in the direction of shear even if the orientation pdf is in a steady state. With a full Smoluchowski equation (i.e. the Doi-Hess model),we find that under a weak shear the orientation probability density function (pdf) is a traveling wave with the location-dependent velocity, that is, the shape of the orientation pdf is maintained while the whole function may shift in the orientation dimension. So in the case of weak shear, the classification of tumbling and wagging makes perfect sense. But in the case of weak shear, it was found that there is only tumbling and there is no wagging [31]. As P e gets larger (i.e. the shear getting stronger), the unsteady evolution of the orientation pdf becomes more complicated. As the peak of the pdf moves, the shape of the pdf changes dramatically for large P e. We feel that the classification of tumbling and wagging is no longer adequate for providing a good picture of the time-evolution of the orientation pdf. We will pursue a comprehensive classification of the pdf in the future work. Here we shall refer the time-periodic unsteady state simply as tumbling.
In this paper we present and discuss phase diagrams for a very wide range of the Peclet number. For small Peclet number, the result of the current analytical/numerical study confirms the conclusions of the previous asymptotic study [31]: When the concentration parameter U is greater than a critical value U 0 (≈ 2.41144646), there is no stable steady state; when U is smaller than U 0 , there is one stable branch of steady state; the stable branch and an unstable branch form a fold and are connected at U 0 . The current study resolves the details of the transition from stability to instability along the phase diagram: the transition occurs around the folding point of the phase curve. As the Peclet number increases, the stable branch and the unstable branch of the fold are peeled off from each other and are separated. For P e > 0.746, the fold disappears completely, that is, the phase curve no longer has a folding point. Furthermore, for P e > 0.746, there is only one steady state (stable or unstable) for each value of U , and the transition from stability to instability occurs at a Hopf bifurcation point.

Two-Dimensional formulation.
We start with reviewing the mathematical formulation of the Doi-Hess kinetic theory for two-dimensional homogeneous flows of rigid rodlike molecules immersed in a viscous solvent subject to an imposed flow field [7,17]. The orientation of a polymer rod is described by an angle θ and the orientational direction of each polymer rod is denoted by u = (cos θ, sin θ).
We consider the steady state solution of the Smoluchowski equation in the form [7] ∂ρ where ρ(θ, t) is the orientational probability density function of the ensemble, ψ(θ) is a periodic function with period π, and f is a constant torque. In the case of nematic polymers under shear, ψ(θ) contains both the Maier-Saupe interaction and the effect of the elongational component of the shear whereas f is caused by the rotational component of the shear. The steady state solution satisfies where J is the steady state probability flux. Multiplying (2) by the integrating factor exp[−f θ + ψ(θ)], we have Integrating from θ to θ + π yields Solving for ρ(θ) gives The value of the probability flux J is determined from the condition π 0 ρ(θ)dθ = 1.
Note that once we know f and ψ(θ), the probability density ρ(θ) is conveniently determined from (5). An efficient method of implementing (5) based on the Fast Fourier Transform (FFT) is discussed in Appendix A.
When an external shear flow with arbitrary strength is imposed, the Smoluchowski equation (1) can be written as where P e = γ/Dr is the Peclet number, a nondimensional parameter measuring the relative strength of viscosity (γ) and rotational diffusivity (Dr). In equation (7), V SH (θ) is the periodic part of the potential induced by the imposed shear, and V MS (θ) is the Maier-Saupe short-range interaction potential, where the brackets · denote the ensemble average with respect to the probability density function ρ(θ) and the angle α is selected such that Let r ≡ U cos 2(θ − α) . With the introduction of r, the Maier-Saupe potential can be expressed as V MS (θ) = −r cos 2(θ − α) ( 1 1 ) and the order parameter can be written as Note that equation (1) reduces to (7) if we use the substitutions f = P e 2 and It is clear that once the values of P e, r and α are given, the total potential is completely specified and the steady state probability density (if a steady state exists) is given by (5). The existence of steady state solution of the Smoluchowski equation (7) depends on the existence of solution of the following nonlinear system: where the pdf used in averaging is the one given in (5). In the nonlinear system (15), P e and U are parameters whereas r and α are unknowns. Alternatively, q and α can be used as unknowns where q and r are related by (14). The goal of this paper is to study the existence and stability of steady solutions of the nonlinear system (15) for P e > 0 and U > 0. It is worth noticing that with the introduction of r, the pdf given in (5) is independent of U . So we can avoid the parameter U in our analysis until the very end. These are the big advantage and mathematical convenience of dealing with r instead of U .

Existence of steady state solutions.
In this section we will study (both analytically and numerically) the existence of solutions of the nonlinear system (15) with 0 < U < ∞ for a given value of P e. It is important to keep in mind that the pdf in (5) is independent of U with the use of r.
Since sin 2θ is a periodic function with period π, we can restrict the range of α to [− π 4 , 3π 4 ]. Also note that sin 2θ satisfies That is, the polymer orientation distribution with order parameter cos 2(θ − α) and phase angle α is identical to the polymer orientation distribution with order parameter − cos 2(θ − α) and phase angle α + π/2. So by allowing r ≡ U cos 2(θ − α) to be both positive and negative, we can further restrict the range of α to [− π 4 , π 4 ]. We first make a change of variable θ new = θ − α. For simplicity, we still denote the new variable θ new by θ and the corresponding pdf and potential for the new variable θ new by ρ(θ) and ψ(θ), respectively. In our notation, (15) becomes where the new pdf ρ(θ) and the new potential ψ(θ) are For a given value of P e, the function F 1 (α, q, P e) depends on only q and α. Our extensive numerical calculations suggest that for any value of P e > 0 and any value of α ∈ [− π 4 , π 4 ], there exists one and only one value of q such that F 1 (α, q, P e) = 0.
A rigorous mathematical proof of the uniqueness is still open. In Appendix B, we prove the existence by showing the asymptotic behaviors of F 1 (α, q, P e): Since F 1 (α, q, P e) is continuous and it has opposite sign at positive and negative infinity, there exists a point such that it is zero. In other words, equation (16) can always be satisfied and for a fixed value of P e it defines uniquely a function q(α) such that F 1 (α, q(α), P e) = 0. Notice that in (16)-(19), α appears only in the potential ψ(θ) in the form of cos 2α. Since cos 2α is an even function with respect to α, the whole problem (16)-(19) is even with respect to α. Consequently, q(α) is an even function of α. Now, in potential ψ(θ) given in (19), we replace q by function q(α). It follows that for a fixed value of P e, the pdf ρ(θ) given in (18) is now a function of α only for α ∈ [− π 4 , π 4 ]. For mathematical convenience, let us introduce several functions as follows: As we discussed earlier, the function s(α) is also an even function of α. As we will see later, the nonlinear system (16) and (17) governing the steady state solutions can be conveniently described by these introduced functions. Furthermore, the nonlinear system (16) and (17) in terms of these introduced functions can be simplified significantly.
It is clear from Figures 1-2 that function q(α), together with functions s(α) and w 0 (α) are not sensitive to the value of P e. In contrast, the function w(α) is highly dependent on P e since the second term in the definition of w(α) has a coefficient P e 4 . Figure 2 shows that for small values of P e, the function w(α) initially increases to a local maximum; after the local maximum w(α) decreases to a local minimum; after the local mininum w(α) starts to increase again. As the value of P e increases, the locations of the maximum and minimum get closer. In other words, as the value of P e increases, the zigzag part of w(α) shrinks. At a critical value of P e c = 0.746027, both the value and the location of the maximum coincide with those of the minimum, and the zigzag part of w(α) disappears. For P e > P e c , function w(α) is monotonically increasing for α ∈ [− π 4 , π 4 ]. Note that if we replace q in equation (17) by q(α) (and replace r by r(α)), then equation (17) is an equation of α only where P e and U are treated as parameters. In terms of function w(α), equation (17) becomes Here we write it as w(α, P e) to show explicitly its dependence on P e. The phase diagram (steady states as functions of U or P e or both) can be obtained by solving equation (20). As we discussed above, for P e < P e c , equation (20) has one solution for small values of U ; it has three solutions for U in an intermediate range; and it has one solution for large U . For P e > P e c , the zigzag part of function w(α) disappears. As a result, for P e > P e c , equation (20) has one solution for any value of U > 0. Figure 3 shows the phase diagram (steady state as a function of U ) for P e = 0.01, corresponding to the case of weak shear. As shown in the figure, a cusp (actually it is smooth, see the insert in Figure 3) separates the diagram into two branches. As we argued in [31], the branch shown in solid line, which corresponds to the stable branch in the case of no shear (P e = 0), is stable. The branch shown in dashed line is unstable [31]. The lower part of the dashed line branch corresponds to the unstable isotropic branch for U > 2 in the case of no shear (P e = 0). The stability of the upper part of the dashed line branch was analyzed in our multi-scale asymptotic study [31]. In the case of weak shear (i.e. P e << 1), a nematic equilibrium will keep its shape. But under the influence of weak shear (P e << 1), its phase angle α will change with respect to the slow time scale T 1 = P e · t. The evolution of the phase angle α is governed by the differential equation [31] dα where c 1 > 0. For −1 < c 2 < 0, the phase angle will converge to a steady state. The stability of a steady state phase angle α 0 is determined by the derivative of sin 2 α + c 2 at α 0 : d dα (sin 2 α + c 2 ) = 2 sin α cos α = sin 2α.  If sin 2α 0 > 0, then α 0 is unstable. If sin 2α 0 < 0, then α 0 is stable. As shown in Figure 4, for any −1 < c 2 < 0, there are two steady state phase angles in (0, π). The one in (π/2, π) is stable while the one in (0, π/2) is unstable. The upper part of the dashed line branch in Figure 3 corresponds to the unstable steady state phase angle, and thus is unstable. In the asymptotic study of weak shear (P e << 1) both the stable branch and the upper part of the unstable branch cease to exist for U > U c = 2.41144646. For U > U c , there is no steady state solution. Instead, there is a tumbling solution. In the analytical/numerical study here, the stable branch and the upper part of the unstable branch are connected to each other at the cusp (see the insert in Figure 3).  In Figure 5 we plot the phase diagram (steady state as a function of U ) for several values of P e. As shown in Figure 5, the cusp, which separates the phase diagram into two branches, is less and less well defined as the value of P e increases. For P e > P e c , the cusp disappears completely. Here it is important to point out that the system is not driven by a potential field because the effect caused by the shear flow is not a potential field. As a consequnce, a steady state branch may change its stability at any point even when it remains isolated. In contract, when a system is driven by a potential field, a steady state branch may change its stability only if some other branch bifurcates from it (or a fold occurs) [8].

Stability study.
For a given value of P e, a steady state solution is completely specified by the value of α. Once the value of α is known, U and the order parameter s are given by U = w(α), s = s(α). The probability density function (pdf) is calculated as follows: where the potential ψ(θ) is ψ(θ) = P e 4 cos 2α sin 2θ − q(α) cos 2θ.
Putting all these in the right hand side of (24) and simplifying, we get Hence, from (24), it follows that the Fourier coefficients of the perturbation satisfy If we collect all of the Fourier coefficients of ∆ρ(θ, t) into a vector u(t) = [a 1 (t), b 1 (t), a 2 (t), b 2 (t), · · · , a M (t), b M (t)] T , then the above equations can be cast in a compact form where we have used the notation A[ρ(θ)] to show explicitly that the matrix A depends on the unperturbed steady state ρ(θ). If all of the eigenvalues of A have negative real parts, then the perturbation decays to zero as time goes to infinity and the steady state ρ(θ) is linearly stable; if at least one of the eigenvalues of A has a positive real part, then the steady state ρ(θ) is unstable. We use the Matlab package to solve this eigenvalue problem to determine the stability.

Numerical results.
Referring to Figure 3, we showed the phase diagram with stability/instability indicated for the case of P e = 0.01. In the figure, the solid line represents the stable branch and the dashed line represents the unstable branch. The numerical stability result is consistent with that of the asymptotic analysis. Numerical stability analysis shows that the stable and the unstable branches are connected at the folding point of the phase curve where the phase curve is vertical (i.e., the point where ∂s/∂U = ∞). The transition point from the stable branch to the unstable branch can only be calculated from the numerical stability analysis. The asymptotic analysis is not capable of resolving this detail. Figure 6 shows the phase diagram with stability/instability indicated for the case of P e = 0.1. The solid line represents the stable branch and the dashed line represents the unstable branch. Numerical stability analysis shows that the transition point from the stable branch to the unstable branch is the folding point of the phase curve where the phase curve is vertical (i.e. the point where ∂s/∂U = ∞).  In Figure 7 we show the phase diagrams with stability/instability indicated for several intermediate values of P e. The solid line represents the stable branch and the dashed line represents the unstable branch. Numerical stability analysis shows that the transition from the stable branch to the unstable branch is the folding point of the phase curve if such a folding point exists. Note that this qualitative behavior based on the closure models has been observed in [21] even though they differ quantatively. For P e < 0.746, such a folding point exists. For P e > 0.746, there is no such folding point. For P e > 0.746, the transition from the stable branch to the unstable branch occurs at a point beyond the maximum.  Numerical stability analysis indicates that the transition from the stable branch to the unstable branch occurs at a point beyond the maximum. Furthermore, as the value of P e increases, the value (U ) of the transition point increases. Figure 9 describes the region of stable steady state (shaded) and the region of tumbling solution in the (U, P e)-plane. This figure is also inherent in [21], where the locus of the Hopf bifurcation was shown. Here, it is captured from the kinetic theory.

Conclusions.
We have semi-analytically investigated all steady state solutions and their stability for planar 2D nematic polymers under an imposed shear flow with arbitrary strength. We have confirmed the mesoscopic closure models of Lee et al. [21], extending their analysis and numerical studies to the Smoluchowski equation for the orientational probability distribution. A detailed phase diagram is given, consisting of all stable steady states versus the normalized concentration U and normalized shear rate or Peclet number P e. Special attention has been given to the locus of bifurcations in the phase diagram where stable steady states no longer exist. These bifurcations are shown to consist of a fold among stable and unstable branches up to some critical P e, after which the fold disappears and is replaced by Hopf bifurcation. Extension of the current work to 3-D model will be much more challenging mathematically. Appendix A An efficient method for calculating π 0 exp −fθ + ψ(θ +θ) dθ Using the Fourier series expansion of exp(ψ(θ)) = ∞ k=−∞ c k exp(ik2θ), we obtain The numerical procedure for calculating π • Apply FFT on ψ(θ) to calculate c k .
Appendix B Existence of solution of F 1 (α, q, P e) ≡ sin 2θ = 0 We consider several cases. Case 1: q → +∞ Since the pdf ρ(θ) given in (18) is periodic with period π, we can consider any interval of length π. Without loss of generality, let us select the interval [− π 2 , π 2 ]. We have where Z = For positive and large q, the dominant contribution in the above two integrals comes from a small region near θ = 0 andθ = π 2 . Away from this small region, the contribution is exponentially small.
In the inner integral with respect toθ, we make a change of variablesθ new = θ+θ. For simplicity, we still useθ to denote the new variableθ new . Since the dominant contribution comes from a small region near θ = 0 andθ = π 2 , we can extend the integration limits away from the region of the dominant contribution. Specifically, we rewrite (30) as Here T.S.T. stands for transcendentally small term. Let We have We apply the Laplace method to the integral in (33). The dominant contribution comes from a small region near θ = 0. To proceed, we need to find the leading term expansion of the part not containing q around θ = 0. We discuss two cases. Case 1A: cos 2α < 1 In this case, we have sin 2θ sinh P e 2 θ − P e 4 cos 2α sin 2θ = P e(1 − cos 2α)θ 2 + · · · , and F 1 (α, q, P e) = 2C 1 Z Similarly, we get Z = 2C 1 π 2 0 cosh( P e 2 θ − P e 4 cos 2α sin 2θ) exp(q cos 2θ)dθ + T.S.T.