Liquid crystals (LCs) display many of the flow characteristics of liquids but exhibit long range orientational order. In the nematic phase, the coupling of structure and flow leads to complex hydrodynamic effects that remain to be fully elucidated. Here, we consider the hydrodynamics of a nematic LC in a hybrid cell, where opposite walls have conflicting anchoring boundary conditions, and we employ a 3D lattice Boltzmann method to simulate the time-dependent flow patterns that can arise. Due to the symmetry breaking of the director field within the hybrid cell, we observe that at low to moderate shear rates, the volumetric flow rate under Couette and Poiseuille flows is different for opposite flow directions. At high shear rates, the director field may undergo a topological transition which leads to symmetric flows. By applying an oscillatory pressure gradient to the channel, a net volumetric flow rate is found to depend on the magnitude and frequency of the oscillation, as well as the anchoring strength. Taken together, our findings suggest several intriguing new applications for LCs in microfluidic devices.

## I. INTRODUCTION

In recent years, liquid crystals (LCs) have found an increasing range of applications that go beyond traditional display technologies to colloidal self-assembly,^{1,2} biosensing,^{3} active biomaterials,^{4} and state-imprinted memory effects.^{5} Their static properties in the bulk, such as topological defects and optical characteristics, have been studied extensively. Less is known, however, about their behavior under confinement and their hydrodynamic response to applied fields.

Advances in microfabrication have led to sophisticated technologies that include ink-jet printing and lab-on-a-chip devices. In these technologies, microfluidic devices can collect, separate, and transport a variety of heterogeneous materials, such as bacteria, colloids, or proteins.^{4,6–8} The liquid medium is typically provided by a simple Newtonian fluid. In this work, we examine the usage of a nematic LC as the medium, where the flow and transport properties may be considerably different from those of the traditional, simple liquids that have been considered to date.

Experimental studies of the flow of nematic LCs have primarily focused on issues related to alignment and tumbling.^{9–16} Jewell *et al.* found two steady state director fields for Poiseuille flow in a homeotropic channel.^{17} Earlier experiments had been conducted by Fishers and Fredrickson,^{18} and the two-steady state phenomenon had been proposed by Ericksen^{19} and simulated by Marenduzzo *et al.*^{20} and Batista *et al.*^{21} Horn and Winter and Holmes *et al.* discovered that a small surface pretilt can induce distinct director profiles for shear and Poiseuille flow.^{22–24} Sengupta *et al.* found that a temperature gradient across the channel can tune the flow shape of Poiseuille flow.^{25} A review of recent experiments on LC microfluidics can be found in Ref. 26.

Theories of nematic LC hydrodynamics date back to the 1960s, when Ericksen, Leslie, and Parodi introduced a nematodynamics model in the director representation, the so-called ELP theory.^{27–30} The average orientation of the local LC molecules is represented by a unit vector **n**, which is the central quantity of any director model. The ELP theory has been widely used to interpret the flow behavior of nematic LCs. It does, however, suffer from limitations, most notably the fact that it cannot describe topological defects. Subsequently, Beris and Edwards proposed a tensorial model based on the Poisson bracket formalism^{31–34} in which a symmetric traceless tensor, **Q**, was introduced to include information about the local average molecular orientation **n** and the degree of order. The advantage of the so-called **Q** approach is that it can describe order parameter variations in the LC. For completeness, we note that other tensorial treatments of nematodynamics are available,^{35,36} some of which were also derived from the Poisson bracket formalism.^{37,38}

Within the structure of a tensor-based description, one can resort to various numerical approaches to describe LC hydrodynamics. A radial basis function method and the Stark-Lubensky formalism have been used to investigate cavity flow.^{39} A finite element method has been applied to solve the Qian-Sheng model.^{40} The lattice Boltzmann method has also been used to solve the Beris-Edwards, Qian-Sheng, or even ELP equations in the context of various LC fluid mechanics problems.^{20,21,41–46} A review of the lattice Boltzmann method for liquid crystal hydrodynamics can be found in Ref. 47.

Although simple (Couette and Poiseuille) and pulsatile flows of nematic LCs have been studied experimentally, numerically, and analytically,^{16,48–51} several key issues remain to be addressed. First, systematic studies of flow in hybrid (splay-bend) cells have not been carried out. Second, systematic numerical studies of unsteady flows are not available. In this manuscript, we develop a surface evolution equation and solve the Beris-Edwards equation, coupled with the Navier-Stokes equation, within the framework of a hybrid lattice Boltzmann/finite difference method.^{52,53} We find that by applying opposite flows to a hybrid cell, the resulting field is asymmetric at low Ericksen numbers. For high Ericksen numbers, the steady-state flow is symmetric, but topological transitions may occur for certain initial director orientations. The anchoring strength of the surface determines the flow-rate-difference of the opposing flows, as well as the transition Ericksen number. We further simulate an oscillatory pressure-driven flow, namely, pulsatile flow, in a hybrid cell and find that a net volumetric flow rate develops, even for very small Reynolds number. The net flow is a function of the pressure gradient, oscillation frequency, and surface anchoring.

## II. METHODS

As mentioned above, our calculations rely on solution of the Beris-Edwards and Navier Stokes equations. The hydrodynamics of the system are described by the density *ρ*, the velocity field **u**, and the pressure *P*. The LC phase is described in a tensorial form. An order parameter is introduced, $Q=\u3008nn\u221213I\u3009$, where the unit vector **n** represents the molecular orientation, **I** is the identity tensor, and 〈〉 denotes an ensemble average. By introducing the velocity gradient *W _{ij}* = ∂

_{j}

*u*, the Beris-Edwards equation reads

_{i}^{34}

which governs the behavior of the LC phase under flow. The part of advection term for rod-like molecules **R** is defined as

where **D** = (**W** + **W ^{T}**)/

**2**and

**Ω**= (

**W**−

**W**)/

^{T}**2**are the symmetric and antisymmetric parts of

**W**, respectively. Here,

*ξ*is a material constant related to the molecular aspect ratio, and Γ is related to the rotational viscosity

*γ*

_{1}of the LC by Γ = 2

*S*

^{2}/

*γ*

_{1},

^{41}where

*S*is the scalar order parameter of the nematic LC, obtained from the eigenvalues of

**Q**. The molecular field

**H**, i.e., the free energy driven term is given by

where $\u2026st$ refers to the symmetric and traceless projection, and $F$ is the total free energy of the system, given by

where *f* and *f _{surf}* are the bulk and surface energy densities, respectively. Here,

*f*takes the Landau-de Gennes form (Einstein notation assumed)

^{54,55}

where *κ* is the elastic constant of the LC (here, we adopt a single-elastic-constant approximation) and is related to the Frank elastic constant *K* by *K* = 2*S*^{2}*κ*. Parameter *U* controls the magnitude of *S* of a homogeneous static system via

The nematic coherence length is $\xi N=K/A0$ and determines the size of a defect in the nematic LC. It is also the natural length scale of the system, as will be seen later.

The Navier-Stokes (momentum) equation for LCs reads

The stress tensor **Π** is decomposed into a viscous stress **Π**^{v} and an elastic stress **Π**^{e},^{56}

where *f*_{0} is the bulk free-energy density of an undistorted nematic and the temperature *T* is related to the speed of the sound *c _{s}* by $T=cs2$. Isotropic pressure is

*P*

_{0}=

*ρT*−

*f*+

*f*

_{0}.

Two types of anchoring can be introduced, namely, non-degenerate or degenerate planar. Non-degenerate surfaces exhibit a preferred **Q** field, **Q**_{surf} = *S*(**n**_{s}**n**_{s} − **I**/3), with a preferred surface molecular orientation **n**_{s}. In the so-called Rapini-Papoular form,^{57} the energy penalty associated with deviations of **Q** from **Q**_{surf} is given by

Degenerate planar surfaces penalize any off-plane distortion of the **Q** tensor; in the Fournier-Galatola form,^{58} they are quantified through

where $Q\u0303=Q+(S/3)I$ and $Q\u0303\u22a5=ZQ\u0303Z$. Here, **Z** is the projection operator associated with the surface normal *ν* as **Z** = **I** − *νν*. The anchoring strength *W* has values that range from 10^{−7} to 10^{−3} N/m,^{59} depending on surface treatments. Similar to the derivation of Beris-Edwards equation (1), one can show that the governing equation of the evolution of the surface **Q**-field is given by

where Γ_{s} = Γ/*ξ _{N}*. The above evolution equation is similar to Eq. (9) in Ref. 43. For steady flow, it reduces to the mixed boundary condition given in Ref. 21.

We use a finite-difference method to solve evolution equations (1) and (7). The lattice Boltzmann method is employed to solve Navier-Stokes equation (3) over a D3Q15 grid^{60} (15 discrete velocities in a 3D regular lattice). Stress implementation follows the approach proposed by Guo *et al.*^{61} A no-slip boundary condition is assumed at all surfaces. Model validation was performed by comparing our simulation results to predictions using the ELP theory.^{27–30}

The characteristic time scale is set by the relaxation time of the liquid crystal, *τ* = *γ*_{1}/*A*_{0}. Typically, *τ* ≈ 1 *μ*s for 4-Cyano-4’-pentylbiphenyl (5CB). The characteristic length scale is set by the nematic coherence length *ξ _{N}*, which corresponds to the size of a LC defect. For 5CB,

*ξ*= 6.63 nm.

_{N}^{59}The two equations are non-dimensionalized according to these time and length scales.

## III. RESULTS

As illustrated in Fig. 1, we consider a channel bounded by two walls in the *z* direction. To form a hybrid cell, we impose homeotropic and planar anchoring at the top and bottom walls, respectively. In the absence of flow with single-elastic-constant approximation, minimization of the elastic energy yields^{30}

Thus, *θ* is linear in terms of the channel depth *z*. There are two degenerate states associated with the above solution for any *ϕ*_{0}, namely, a counterclockwise state (para-field), with *θ* decreasing from 90° at the bottom wall to 0° at top wall, and a clockwise state (dia-field), with *θ* increasing from 90° at the bottom wall to 180° at the top wall. The flow is introduced by either moving the top wall in the *x*-direction at velocity *v*_{0} (Couette flow) or by imposing a pressure gradient ∂_{x}*P* in the *x*-direction (Poiseuille flow). Applying the same flow type to the two degenerate states is equivalent to applying two opposite flow directions to a single state. Therefore, we focus on the counterclockwise director field (para-field) here and apply flow in the two opposite directions. An Ericksen number *Er* is introduced to describe the ratio of the viscous force and the elastic force

where *η* is the viscosity and *v* is the characteristic velocity. When *Er* ≪ 1, the director field is close to that of a static system. When *Er* ≫ 1, viscous effects dominate and the directors align with the flow. The Reynolds number *Re* = *ηvL _{z}*/

*ρ*in our simulations ranges from 10

^{−6}to 0.2. We are particularly interested in measuring the orientation of the directors described by

*θ*and

*ϕ*, as this information can be directly translated into optical images under cross polarized light that could then be compared to experiments. An additional, experimentally relevant parameter is the volumetric flow rate Φ, which is defined as

Here, *A* is the cross section of the channel (a plane whose normal is in the *x* direction), and 〈〉_{t} represents a time average (for unsteady flow). In this work, we assume that the LC is nearly incompressible (Mach number *Ma* ≪ 1; thus, density *ρ* is a constant). Therefore, the volumetric flow rate is equivalent to the mass flow rate. In the following, we occasionally refer to them simply by “flow rate.” We say that a flow is symmetric if the opposite flow directions yield identical flow fields (flow rates); otherwise, we say that the flow is asymmetric.

Simulation parameters are chosen so that the system is in the flow-aligning regime. The channel height is *L _{z}* = 51

*ξ*, the molecular constants are

_{N}*ξ*= 0.8 and

*U*= 3.5, implying that the equilibrium order parameter is

*S*≃ 0.62. The LBGK (lattice Boltzmann Bhatnagar–Gross–Krook) single relaxation time is chosen to be

_{bulk}*τ*= 0.84. Surface anchoring sets an additional length scale, namely, the Kleman-de Gennes length

_{f}*ξ*≡

_{KD}*K*/

*W*.

^{30}To avoid uniform director fields in our hybrid cell,

*ξ*<

_{KD}*L*requires that we work with anchoring strengths of magnitude

_{z}*W*≥ 10

^{−4}N/m.

### A. Couette flow in a hybrid cell

We first examine Couette flow in a hybrid cell. The bottom wall is stationary and the top wall moves with velocity *v*_{0} or −*v*_{0} in the *x*-direction (*v*_{0} > 0 is assumed). We only consider rigid (infinite) anchoring for Couette flow. We measure the flow rates Φ_{+} for para-field and Φ_{−} for dia-field and plot the difference ΔΦ = |Φ_{−}| − |Φ_{+}| against Ericksen number *Er*. Our results are shown in Fig. 2(a). All flow rates are scaled by Φ_{0}, which is the flow rate of an isotropic liquid with the same viscosity in the same channel. At moderate and low Ericksen numbers, shearing the LC along opposite directions with identical shear rate yields different flow rates, i.e., ΔΦ≠0, and it is linear with respect to *Er*. ΔΦ is largest when the azimuthal angle is *ϕ*_{0} = 0°. Before |Φ_{−}| and |Φ_{+}| become identical, at high *Er* number, |Φ_{+}| surpasses |Φ_{−}| in the transient regime. Fig. 2(b) shows the velocity fields with different azimuthal angles *ϕ*_{0} at *Er* = 5. The velocity profiles of the two states are different: the dia-field has higher velocities in the lower region of the channel. The difference stems from the different responses of the director field to the flow. In Fig. 2(c), we plot the angle profile *θ* in terms of channel depth *z* at *Er* = 5. One observes that *θ* is no longer linear with respect to channel depth. When the flow is weak compared to the elastic response (*Er* < 1), the directors of the para-field tend to align with the flow by increasing *θ*, whereas for the dia-field, the directors rotate against the flow and *θ* decreases. For infinite anchoring conditions, we see super-linear behavior of *θ* for the para-field and sub-linear behavior for the dia-field. When *ϕ*_{0} = 0, the effective shear viscosity can be determined from *θ* via

where *α*_{1} is a Leslie viscosity coefficient and *η ^{b}* and

*η*are the Miesowicz viscosities with

^{c}*η*<

^{b}*η*. The apparent shear viscosity,

^{c}*η*, is plotted in the inset of Fig. 2(b). Due to the viscous anisotropy of the LC, the shear viscosity is inhomogeneous across the channel. In a hybrid cell filled with 5CB at room temperature (25 °C), the shear viscosity can vary by a factor of 5.

_{eff}^{62}A greater

*θ*leads to a lower viscosity, which in turn supports a larger velocity gradient. Thus, the velocity drops from

*v*

_{0}at the top wall to 0 at the bottom wall earlier in the para-field than in the dia-field, leading to a smaller flow rate in para-field than that in dia-field. In Figs. 2(b) and 2(c), we also calculate the velocity and angle profiles from ELP theory (solid and dashed lines). Our model is in agreement with the ELP nematodynamics theory.

When viscous effects overcome elastic effects (*Er* ≥ 10), the dia-field becomes unstable and transits to para-field, which leads to identical flow rates. The time evolution of the directors is shown in Fig. 2(d). It is known that at high *Er*, the directors adopt the Leslie angle^{30}

where *α*_{2} and *α*_{3} are Leslie viscosity coefficients. For our system, *θ _{Leslie}* ≈ ± 76° (compared to

*θ*≈ ± 78° for 5CB at room temperature

_{Leslie}^{62}). At the early stage of evolution, directors near the top wall adopt the Leslie angle. They conflict with the directors near the bottom. Thus, a high elastic distortion appears and gradually approaches the bottom wall. Later, the high-elastic region becomes unstable. The order parameter drops and the directors adopt Leslie angles thereafter. At steady state, the order parameter recovers the bulk value and the channel becomes defect free. We should notice that dia-field and para-field states are topologically different. The defect has to form during the transition from a dia-field to a para-field. The critical Ericksen number

*Er*≈ 50 is consistent with a similar experimental result, which reports that a uniform homeotropic liquid crystal (5CB) subjected to Poiseuille flow shows a topological transition at

_{c}*Er*≈ 25.

_{c}^{17}

### B. Poiseuille flow in a hybrid cell

We next consider Poiseuille flow in a hybrid cell. Here, we limit our analysis to *ϕ*_{0} = 0° with finite anchoring. Again we plot the flow-rate-difference ΔΦ = (|Φ_{−}| − |Φ_{+}|)/Φ_{0} between para-field and dia-field states against Ericksen number *Er* in Fig. 3(a). ΔΦ increases to a maximum at low Er and decreases to 0 at high *Er*. The threshold Ericksen number, *Er _{c}*, at which ΔΦ vanishes is dependent on the anchoring strength

*W*: higher

*W*leads to higher

*Er*. As shown in Fig. 3(b), velocity fields at ∂

_{c}*P*/∂

*x*= 10

^{5}Pa/m (

*Er*≈ 5) are plotted for both para- and dia-fields. Anchoring strength

*W*has the opposite effect on the two fields. For para-field, the weakest anchoring leads to the lowest flow rate, whereas for dia-field, lower

*W*significantly increases the velocities. The underlying reason can be understood by relying on Fig. 3(c). For the para-field,

*θ*decreases throughout the channel due to flow alignment, and the effective shear viscosity also increases. The resulting flow rates have decreased. Dia-field director fields respond to the flow in an opposite way: the flow increases

*θ*and the effective shear viscosity has therefore decreased, leading to increased flow rates. When

*W*= 10

^{−4}N/m, the relative weak anchoring exacerbates the above flow alignment effect and, consequently, the flow-rate-difference becomes the greatest. Note that the shape of the velocity fields tilts to the bottom wall in most cases, due to the fact that the shear viscosity near the bottom wall is lower than that near the top wall. The position of the velocity peak varies with anchoring strength. This implies that one can shape the velocity field by tuning surface chemistry.

^{25}

We next address the origin of the peak in the flow-rate-difference plot of Fig. 3(a). When *Er* ≤ 1, the flow rate increases as the pressure gradient increases and, therefore, the flow-rate-difference between opposite flowing directions depends monotonically on the Ericksen number. When *Er* approaches the transition value *Er _{c}*, the relation between

*Er*and Φ

_{±}is no longer linear. As the directors approach the Leslie angle, they have less space to rotate, and the difference between the two fields decreases. Recall that, for Couette flow, the dia-field is stable only at low

*Er*and, for para-field, it is stable for the entire range of

*Er*’s. In Fig. 3(d), we show the stability diagram for Poiseuille flow. The dia-field is unstable in a narrow window of

*Er*numbers. The transition between both fields being stable and para-field being stable depends on the anchoring strength

*W*. This transition is what we observe from Fig. 3(a) when the flow-rate-difference vanishes. However, the transition between the para-field being stable and the dia-field being stable at a higher

*Er*is independent of

*W*. In Fig. 4(a), we show how the directors of a para-field transit to a dia-field at

*Er*≈ 500, where transient defects emerge and disappear in the channel center. The upper and lower regions of the channel adopt Leslie angles of opposite sign. The conflict happens at the lowest velocity gradient region (near the channel center), which becomes unstable and eventually rotates the directors from vertical to horizontal. In Fig. 4(b), we show the time evolution of the directors of a dia-field at

*Er*≈ 250. For this smaller velocity gradient, the directors in the lower region of the channel also adopt the Leslie angle. Frustration happens near the channel center, where the order parameter drops from its bulk value of 0.62 to 0.5 and stabilizes at this specific value. For the dia-field state at even lower

*Er*, the lower region of the channel does not adopt a Leslie angle and remains at an opposite direction with respect to the velocity gradient. The anchoring of the walls controls the angle at steady state and therefore controls at what shear rate the directors conflict with the surface preferred angle. Therefore, the first transition, between two stable fields and one stable para-field, is related to the anchoring strength. However, for the second transition, the stability of the frustration region is only dependent on the local velocity shear and thus independent of the anchoring strengths of the walls.

We have examined the temporal velocity field and the pressure profile before the flow becomes steady for both Couette and Poiseuille flows. When the system is free of defects, the 2D velocity field gradually approaches steady-state configuration. Whereas when the system undergoes a topological transition, the velocity near the transient defect shows a sudden jump, followed by a decaying oscillation about the steady solution. Provided that our calculations are focused on a creeping flow condition, we do not see any evidence of cross flow or other complex flow patterns in our results. The associated pressure profile behaves correspondingly. The low hydrostatic pressure region arises near the domain where the LC exhibits a large distortion, which leads to high local elastic stress.

### C. Pulsatile flow in a hybrid cell

Building on the idea that flow in a hybrid cell is asymmetric, we impose an oscillatory pressure gradient ∂_{x}*P* = *G*sin(*ωt*) to the channel, where *G* is the amplitude and *ω* = 2*π*/*T* is the frequency and *T* is the period. It is known that a Newtonian fluid at low Reynolds number (Stokes flow) shows no net flow subjected to slow oscillatory pressure drop. For a viscoelastic liquid, this is not true even at low Reynolds number. The Womersley number, *α*, is defined as $\alpha =Lz/2\nu /\omega $.^{63} The kinematic viscosity *ν* = *γ*_{1}/*ρ*. Here, *α* ranges from 0.12–6.0. In Fig. 5, we plot the net flow rate |Φ| against the oscillation period *T*. For small pressure gradient, |Φ| is inversely dependent on *α*. Rapid oscillations lead to a low net flow, and slow oscillations yield a high net flow. This observation can be understood from the velocity and director profiles in Fig. 6. When *T* is comparable to the LC relaxation time *τ*, the director has little time to respond to the pressure gradient and thus remains linear; the associated velocities are small as well. As the frequency decreases, or when *T* > 10*τ*, one starts to see the asymmetric oscillations of the director, as well as the velocity fields, which contributes to a greater net flow. In the extreme case, when *T* ≫ *τ*, the outlines of the velocity and director angle profiles are essentially the steady-state solutions of the dia and para-fields. The magnitude of the extreme velocity reaches a maximum when *T* ≈ 10*τ*. This is because when the oscillation period equals the viscous time (≈5*τ*), the LC directors resonate with the pressure oscillations. At the limit of high frequency (*T* is small or *α* is large), weaker anchoring leads to higher net flow. This is consistent with what we have found in steady Poiseuille flow, where weaker surface anchoring allows more elastic deformation, resulting in higher asymmetry of the flow fields. For a moderate or high pressure gradient, the net flow only exists at high frequencies. As shown in Figs. 5(b) and 5(c), the transition period at which |Φ| vanishes is dependent on the pressure gradient, but the surface anchoring has little influence. Fig. 7 shows the details of the underlying mechanism. For rapid pressure oscillations, the director has no time to perform a topological transition (as in Fig. 4) and |Φ| remains nonzero. However, when the pressure oscillation slows down, a topological transition occurs and the director of the bulk becomes horizontal. The transition of the directors being vertical to horizontal drastically lowers the effective shear viscosity, and the extreme velocity increases accordingly. As a consequence of the topological transition, only the director near the top wall oscillates according to the pressure, the responding flow becomes symmetric, and the net flow vanishes. As found in the steady flow, the higher the pressure gradient, the shorter time the LC needs to undergo a topological transition. Therefore, for smaller oscillation magnitudes, the transition period is longer than that observed at higher oscillation magnitudes.

## IV. CONCLUSION

In this work, we have developed a surface evolution equation and coupled it with a lattice Boltzmann solver to examine the flow of liquid crystals described by the Beris-Edwards equation. Finite anchoring boundary conditions and time dependent dynamics can be handled by our method. We observe that at low and moderate Ericksen numbers (*Er* ≤ 1), the Couette and Poiseuille flows in a hybrid cell are asymmetric and dependent on the direction of the shear or pressure gradient. The difference in flow rates is linear in terms of the Ericksen number. At high Ericksen numbers (*Er* ≥ 10), the director fields may transit from one state to the other through a topological transition during which a transient defect emerges. The resulting flow becomes symmetric at steady state. The transition Ericksen number is dependent on the surface anchoring strength. Our low Ericksen calculations are in agreement with predictions of the ELP theory, which cannot describe topological transitions. Our high Ericksen number observations are consistent with an experimental result in which a critical Ericksen number of the same order of magnitude was reported. We then extend the method to simulate pulsatile flow in a LC channel. A sinusoidal pressure gradient is applied to the hybrid cell. We show that the net flow depends on the oscillation frequency and the anchoring strength, as well as the pressure magnitude. At high pressures and low frequency, the LC undergoes a topological transition, the resulting flow field becomes symmetric, and the net flow rate vanishes.

Our surface evolution equation enables us to systematically study the temporal behavior of LC hydrodynamics. The asymmetry of the hybrid cell due to flow stems from a spontaneous symmetry breaking of the splay-bend structure. For more complicated director configurations, one should expect similar behaviors. Unlike Newtonian fluids, which have no such internal structure, viscoelastic liquids could find useful applications in future microfluidic devices.

## Acknowledgments

The studies of transient, non-equilibrium behavior of flowing liquid crystals presented here were supported by the National Science Foundation, Grant No. DMR-1410674. The authors are grateful for the support from the University of Chicago and Argonne National Laboratory for development of innovative collaborative efforts. The development of robust algorithms and computational codes for simulation of nanostructured complex fluids was supported by the Department of Energy, Basic Energy Sciences, Division of Materials Research, through the Midwest Integrated Center for Computational Materials (MICCoM). The research of I.S.A. was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Science and Engineering.