Shear-Driven Mechanism of Temperature Gradient Formation in Microﬂuidic Nematic Devices: Theory and Numerical Studies

: The purpose of this paper is to show some routes in describing the mechanism responsible for the formation of the temperature difference ∆ T at the boundaries of the microﬂuidic hybrid aligned nematic (HAN) channel, initially equal to zero, if one sets up the stationary hydrodynamic ﬂow v st or under the effect of an externally applied shear stress (SS) to the bounding surfaces. Calculations based on the nonlinear extension of the classical Ericksen–Leslie theory, supplemented by thermomechanical correction of the SS σ zx and Rayleigh dissipation function while accounting for the entropy balance equation, show that the ∆ T is proportional to the heat ﬂux q across the HAN channel and grows by up to several degrees under the inﬂuence of the externally applied SS. The role of v st = u st ( z ) ˆ i with a sharp triangular proﬁle u st ( z ) across the HAN in the production of the highest ∆ T is also investigated.


Introduction
Consisting of anisotropic molecules, liquid crystal (LC) materials were called curious soft matter until their impressive impact on modern technology. The primary technological revolution was brought by these LC materials in the field of displays. With the development of the LC display market, the question arises about the following areas of application of LC materials. Perhaps there is no more suitable direction for the application of LC materials than LC sensors (LCSs) and LC actuators (LCAs) [1]. They have various advantages in comparison with other types of microsensors and microactuators; simple structure, high shape adaptability, easy downsizing, and low driving voltages. This is because LC materials are extremely sensitive to external disturbances and can be used for the construction of stimuli-responsive devices, such as LCSs or LCAs [1]. Nematic liquid crystal (NLC) channels or capillaries of appropriate size are microdevices in which the molecular orientations can be manipulated by forces applied macroscopically or can be generated locally within the microfluidic nematic channel [2] or capillary [3]. A challenging problem in all such systems is the precise handling of nematic microvolume, which in turn requires self-contained micropumps with small package sizes exhibiting either a very small displacement volume (displacement pumps) or a continuous volume flow (dynamic pumps). One of the pumping principle in the microsized LC channel confined between two infinitely long boundaries is based on the coupling between the tangential component of the shear stress (SS) σ zx and the director fieldn, which describes the preferred orientation of the molecules, together with accounting the effect of the temperature gradient ∇T [4,5].
It has been shown that, due to the temperature gradient ∇T, the horizontal hybrid aligned nematic (HAN) microfluidic channel, being initially at rest, starts moving in the horizontal direction if heated from below or above [6][7][8]. In the case when the directorn is anchored homeotropically to the cooler (T lw ), lower boundaries and homogeneously to the hotter (T up ) upper boundaries, due to coupling between ∇T ∼ ∆T d and ∇n, the hydrodynamic flow v = v xî = uî [7] in the horizontal direction is excited. Here, ∆T = T up − T lw is the temperature difference on the HAN boundaries, d is the thickness of the HAN channel, andî is the unit vector taken parallel to the horizontal boundaries of the HAN channel (see Figure 1). The magnitude of the flow v x is proportional to ∼ d η σ tm zx [6][7][8], where σ tm zx ∼ ξ ∆T d 2 is the tangential component of the thermomechanical stress tensor σ tm ij , η is the viscosity, and ξ is the thermomechanical constant [6]. The direction of the hydrodynamic flow v is influenced by the character of the preferred anchoring of the average molecular directionn to the boundaries of the HAN channel and the heat flux q across the bounding surface [7,8]. Measurements of the temperature of the induced flow were performed on the HAN cell [9], and the main result of that experimental study is the estimation of the thermomechanical constant ξ ∼ 10 −12 J/K m [9].
Despite the fact that the possibility of formation of hydrodynamic flows in nematic channels under the influence of temperature gradients has been theoretically described since [6][7][8], only detailed numerical simulations performed within the framework of the extended Ericksen-Leslie theory [10,11] allowed us to recreate the complete picture of the formation of flows in nematic microchannels and capillaries. One of the aims of this paper is to describe the various regimes of hydrodynamic flow formation due to the interaction between the temperature and director field gradients obtained by numerical modeling of these processes. This paper is devoted to the possibilities of computational methods implemented in the framework of the nonlinear extension of the Ericksen-Leslie theory while accounting for the entropy balance equation [12]. Another purpose of our paper is to show some routes in describing the mechanism responsible for formation of the temperature difference ∆T on the boundaries of the HAN channel, being initially equal to zero, if one sets up the stationary hydrodynamic flow or under the effect of the externally applied shear stress to the bounding surfaces.
It should be noted that a whole class of confined active fluids under external shear stress [13][14][15] as well as the problem of thermomechanical coupling in cholesterics [16] are not included in the scope of our research interests. And this is despite the fact that research in active matter has been fostered by many unexpected finding, ranging from spontaneous flow to symmetry breaking in exotic active emulsions. This is the first review to describe the role of thermomechanical force in the formation of hydrodynamic flow in microsized nematic channels and capillaries in detail. It is based on the nonlinear extension of the Ericksen-Leslie theory, supplemented by thermomechanical correction of the shear stress and the Rayleigh dissipation function, and takes into account the entropy balance equation. The fact that the main results were obtained by numerical methods indicates that the experimenters still have a lot of work to do in order to create a more complete picture of formation of hydrodynamic flows in microsized nematic channels and capillaries.
The layout of this section is as follows. The following section presents the system of hydrodynamic equations describing both the reorientation of the director field and the flow of liquid crystal in a microfluidic HAN channel containing the temperature gradient under the action of the external shear stress. The numerical results for the number of hydrodynamic regimes, caused by both the shear stress and the heat flux across the bounding surfaces of the HAN channel, describing the orientational relaxation of the director, velocity, and temperature are given in Section 2. The role of flow in temperature gradient formation across a hybrid aligned nematic channel is given in Section 3. The conclusions are summarized in in Section 4.

Shear-Driven Flow Regimes in Microfluidic Nematic Devices: Tumbling and Laminar
In the case of the stationary shear stress flow v st = u stî , the directorn is oriented in the shear x − z plane, where the x − z plane is defined by the liquid crystal flow (the direction x coincides with the unit vectorî) and the velocity gradient ∇v in the z direction coincides with the unit vectork; y is the vorticity axis coinciding with the unit vectorĵ (see Figure 1). It has long been thought that, in shear flows, the dynamics of nematics always produces an alignment regime, where the directorn aligns at a stationary angle [17][18][19], with respect to the direction of the flow velocity v st =γyî, when the hydrodynamic torque, exerted per unit LC volume in a shear flow vanishes. Here, γ 21 = −γ 2 /γ 1 , γ 1 = α 3 − α 2 and γ 2 = α 3 + α 2 are the rotational viscosity coefficients (RVCs), α 2 and α 3 are the Leslie coefficients, andγ = ∂v x /∂z is the shear rate. However, it has been found that some LC materials exhibit an unusual type of instability, when the directorn continuously rotates in the shear flow. It is clear from Equation (1) that, if |γ 1 | > |γ 2 | or α 3 > 0 (because, in practice, α 2 < 0), then no real solution for θ st exists. Physically, this means that, in this case, the director tumbles under the shear flow of the nematic. An interesting aspect of the tumbling rotation of the director's field in the shearing flow should be noted here [20]. In this case, there is a certain value of the shear rate, above which the director rotates out from the x − z plane, and there is a steady-state solution without tumbling.
Among the many questions that arise in this connection, we are interested in two. First, how does the viscous torque T vis = γ 1 2 (1 + γ 21 cos 2θγ)ĵ affect the character of the director fieldn (or the polar angle θ) evolution to its stationary orientationn st with respect to the nematic flow v st = u stî in the microsized HAN channel when the director is strongly anchored to both boundaries under the influence of the tangential component of the shear stress σ zx while accounting for the temperature gradient ∇T? This is investigated for two types of nematic phases: first, for the "laminar" case of nematic phase, when γ 21 > 1, and, second, for the "tumbling" case of nematic phase, when γ 21 < 1 [17][18][19]. For instance, the liquid crystal composed of 4 − cyano − 4 − pentylbiphenyl (5CB) molecules belongs to the laminar nematic phase, whereas the liquid crystal composed of 4 − cyano − 4 − octylbiphenyl (8CB) molecules [21] belongs to the tumbling nematic phase. Second, is it possible to set up a temperature difference ∆T across the HAN microfluidic channel, initially being equal to zero, under the action of the tangential component of the shear stress σ zx applied to the boundary of the LC channel?
The answers to these questions are given within the framework of the nonlinear extension of the classical Ericksen-Leslie theory [10,11], supplemented by thermomechanical correction of the shear stress and Rayleigh dissipation function while taking into account the entropy balance equation [12]. It has recently been shown, both experimentally [22] and numerically [23], that only stationary flow u st with a triangular sharp profile and position of the maximum in the vicinity of the restricted boundaries may achieve the highest temperature difference ∆T = T up − T lw in the HAN microfluidic channel of few degrees. By means of other hydrodynamic flow with profiles, which cannot demonstrate the sharp growth of u st across the HAN channel, one can achieve the same result only by using of the "high-speed" hydrodynamic flow ∼0.1 µm/s [22]. However, taking into account that some LC systems driven by external SS exhibit such non-equilibrium phenomena as a tumbling behavior [7,8,23], the mechanical description can shed some light on the problem of temperature gradient formation; when under certain conditions, the thermomechanical force can overcome elastic, viscous, and anchor forces and cause a temperature gradient across the HAN channel.

Formulation of the Relevant Equations for Dynamical Reorientation in Microsized Nematic Fluids
First, we consider the description of the physical mechanism responsible for the shear-driven nematic flow in microfluidic hybrid aligned nematic (HAN) channels under the action of the external shear stress applied to the upper boundary of this channel (see Figure 1): We consider a hybrid aligned channel composed of both the laminar and tumbling types nematics, which is bounded by two horizontal surfaces at a distance of d on a scale of the order of tens of micrometers. According to this geometry, the director is maintained in the xz-plane (or in the yz-plane), defined by the heat flux q = Q 0k normal to the horizontal boundaries of the LC channel. As we deal with the HAN channel under the influence of both the SS σ 0 zx and the heat flux q perpendicular to the HAN channel, taking into account the fact that the length of the channel L is much greater than the thickness d, it can be assumed that the component of the directorn = n xî + n zk = sin θ(z, t)î + cos θ(z, t)k as well as the rest of the physical quantities depend only on the zcoordinate and time t. Here, θ denotes the angle between the director and the unit vectork (see Figure 1). In order to understand how the viscous T vis , elastic T el and thermomechanical T tm torques as well as the tangential component of the shear stress σ zx affect the character of the director fieldn evolution to its stationary orientation with respect to the nematic flow v, we must formulate the boundary conditions for the temperature T(z, t), velocity v, and the director n fields.
We consider a hydrodynamic regime where the HAN channel is subjected to uniform heating from above, for instance by the laser irradiation [24], while directorn is strongly anchored to both solid surfaces, homeotropically to the lower cooler (T 1 ), and homogeneously to the upper bounding surfaces, where whereas the boundary conditions for the temperature field are Here, λ ⊥ is the heat conductivity coefficient perpendicular to the directorn whereas Q 0 is the heat flux across the upper boundary. As a result, we obtaina picture where there is a balance between the heat flux q; SS σ 0 zx ; and the viscous, elastic, and anchoring forces, and in general, the LC fluid settles down to a stationary flow in the horizontal direction [7,8]. Under the assumption of an incompressible fluid, the hydrodynamic equations describing the orientation dynamics induced by both SS σ 0 zx and q can be derived from the torque, linear momentum, and the entropy balance equations for such a LC system.
Taking into account the micro-size of the HAN channel, we can assume that the mass density ρ is constant across the LC film, and thus, we deal with an incompressible liquid. The incompressibility ∇ · v = 0 implies that there is only one nonzero component of the vector v, viz. v(z, t) = u(z, t)î.
When the temperature gradient ∇T is set across the HAN channel, we expect that the temperature field T(z, t) satisfies the entropy balance equation [7,8,12] where Q = −T δR δ∇T is the heat flux and C P is the heat capacity of the liquid crystal phase. To describe the evolution of the director fieldn (or the polar angle θ(z, t)) to its stationary orientationn st (z) and exciting the velocity field v(z, t) caused by both the heat flux q and the external SS σ 0 zx , we consider the dimensionless analog of the torque and linear momentum balance equations as well as the entropy balance equation.
In order to elucidate the role of both the heat flux q = q 0k and the external SS σ 0 zx on the reorientation process in the microsized HAN channel, we consider the hydrodynamic regime when the directorn is strongly anchored to both solid surfaces, homeotropically to the lower, cooler boundary (χ 1 ), whereas on the upper boundary, it is assumed that the heat flux is vanished or restricted. In this case, the boundary conditions must satisfy the following equations is the dimensionless heat flux across the upper boundary of the HAN channel.
The velocity on the lower boundary must satisfy the no-slip boundary condition, whereas on the upper boundary the SS is applied as Now, the reorientation of the director in the microsized HAN channel confined between two solid surfaces, when the relaxation mode is governed by viscous, elastic, thermomechanical forces and the SS σ 0 zx with accounting for the heat flux q = q 0k , can be obtained by solving the system of nonlinear partial differential Equations (9), (10), and (12), with the appropriate boundary conditions for the polar angle θ(z, τ), temperature χ(z, τ) (Equation (13)), and the velocity u(z, τ) (Equations (14) and (15)) as well as with the initial condition

Numerical Results for the Relaxation Modes in the HAN Channel
First, we focus on the problem of how much the viscous torque T vis = γ 1 (χ)θ ,τ − A(θ)u ,z = γ 1 (χ)θ ,τ − 1 2 (1 + γ 21 cos 2θ)u ,z influences the evolution of the director field n(z, τ) (or the polar angle θ(z, τ)) to its stationaryn st (z) distribution across the microfluidic HAN channel with the temperature gradient. In our case, the ∇χ is set by the heat flux q (see Equation (13)) directed across the microfluidic HAN channel.
Calculations of the temperature dependence γ 21 = −γ 2 /γ 1 as well as a comparison of the RVCs values γ 1 and γ 2 , both for 5CB and 8CB, at temperatures corresponding to the nematic phase, are given in Table 1.
The set of parameter values involved in Equations (9), (10), and (12) is δ 1 ∼ 24, δ 2 ∼ 2 × 10 −6 , δ 3 ∼ 6 × 10 −4 , and δ 4 ∼ 10 −10 . Using the fact that δ 2 1, the Navier-Stokes Equation (10) can be considerably simplified as velocity adiabatically follows the motion of the director. Thus, the whole left-hand side of Equation (10) can be neglected, reducing it to while Equation (12) can also be significantly simplified. Since both parameters δ 3 and δ 4 1, and the entire left part of Equation (12) and the second term can be neglected, so the Equation (12) takes the form: The last equation has a solution From a physical point of view, this means that the temperature field χ(z, τ) across the HAN cell under the above conditions is proportional to the heat flux q 0 across the upper bounded surface when the temperature on the lower surface is kept constant.
In the case when the SS σ 0 zx is equal to 10 (∼5 Pa) and there is the heat flux q 0 = 0.02 (Q 0 ∼ 200 nW/µm 2 ) directed to the bulk of the nematic channel, the evolution of the director fieldn to its stationary orientationn st in the microsized HAN channel, which is described by the polar angle θ(z, τ k ) for different times starting from τ 1 = 0.001 (curve 1) to τ R = τ 7 ∼ 0.4 (∼0.07 s) (curve 7) for both cases 5CB (see Figure 2a) and 8CB (Figure 2b), is shown in Figure 2. All calculations in this paper were carried out by the numerical relaxation method [28], whereas the relaxation criterion = |(θ (m+1) (z, τ) − θ (m) (z, τ))/θ (m) (z, τ)| was chosen to be equal to 10 −4 and, then, the numerical procedure was carried out until a prescribed accuracy was achieved. Here m is the iteration number.
First, the effect of the viscous torque T vis , or γ 21 = −γ 2 /γ 1 , on the evolution of the velocity field u(z, τ) is manifested in the qualitative difference in the velocity profiles for 5CB and 8CB nematics. In the case of 5CB, we have concave profiles (see Figure 3a), while in the case of 8CB, these profiles represent almost linear dependencies at the final stage of evolution, where the velocity u(z, τ k ) increases from zero (u(z = 0, τ k ) = 0) at the lower boundary of the channel to the value u(z = 1, τ k ) ∼ 22 (∼0.7 mm/s) at the upper boundary. In the case of 5CB, the value of velocity u(z = 1, τ k )(5CB) at the upper boundary is equal to ∼ 23 (∼0.73 mm/s). Second, the main effect of the viscous torque T vis , or γ 21 = −γ 2 /γ 1 , is manifested in the character of evolution of the director field n to its stationary orientationn st in the microsized HAN channel, which is described by the polar angle θ(z, τ k ). Indeed, in the case of 5CB, the polar angle θ(z, τ k )(5CB) increases monotonically from 0 to ∼1.57 (π/2), whereas in the case of 8CB, the polar angle θ(z, τ k )(8CB) increases monotonically from 0 to θ(z = 0.64, τ k )(8CB) ∼ 2.48 in the vicinity of the centrum of the HAN channel, with a subsequent decrease to the value of ∼1.57 (π/2) at the upper boundary of the HAN channel. Thus, the main effect of γ 21 = −γ 2 /γ 1 is to influence the nature of the reorientation of the director fieldn to its stationary orientation n st in the microsized HAN channel, which is described by the polar angle θ(z, τ k ). In the case of the tumbling type nematic phase, composed of 8CB molecules, when |γ 1 | > |γ 2 |, the director tumbles under shear flow of the nematic, whereas in the case of the laminar type nematic phase, composed of 5CB molecules, when |γ 1 | < |γ 2 |, the dynamics of nematic liquid crystals produces the alignment regime.
In turn, when the SS σ 0 zx is increased and equal to 20 (∼10 Pa) (see Figure 4a) and 30 (∼15 Pa) (see Figure 4b) and there is a heat flux at 0.02 (Q 0 ∼ 200 nW/µm 2 ) in the case of the tumbling type nematic phase composed of 8CB molecules, when |γ 1 | > |γ 2 |, the evolution of the directors fieldn to its stationary orientationn st in the vicinity of the centrum of the HAN channel undergoes a qualitative change.  Figure 5. (a) The relaxation of the velocity field u(z, τ k ) to its stationary distribution across the HAN microfluidic channel in the case of the tumbling type nematic phase composed of 8CB molecules and under the effect of SS σ 0 zx = 20 (∼10 Pa) for different times starting from τ 1 = 0.001 (curve 1) to τ R = τ 7 ∼ 0.4 (∼0.07 s) (curve 7). Here, τ k = 0.067 × (k − 1) (k = 2, ..., 7). (b) The same as in (a) but with SS σ 0 zx = 30 (∼ 15 Pa) and τ R = τ 7 = 0.6 (∼0.1 s). Here, τ k = 0.1 × (k − 1) (k = 2, ..., 7), and q 0 is equal to 0.02. According to our calculations, the shear stress σ 0 zx produces the velocity field u(z, τ) directed in the positive direction (see Figure 5) and its effect on the director distribution across the HAN microfluidic channel is so strong that, in the middle part of the nematic channel, the biggest value of the polar angle θ(z, τ) is equal to 5.5 (∼315 • ) at σ 0 zx = 30 (∼15 Pa) and the director practically executes a full cycle of rotation (see Figure 4b). That influence decreases with a further decrease in σ 0 zx . However, taking into account that the director field is strongly anchored to both boundaries of the HAN channel, homeotropically to the lower and homogeneously to the upper, the balance of the viscous, elastic, thermomechanical, and anchoring forces and the SS σ zx applied to the upper restricted surface leads to rotation of the director field mainly in the middle part of the HAN microfluidic channel.
The maximum absolute value of the dimensionless velocity u st (z) = γ 10 d K 10 v st x in the microsized HAN channel at the final stage of the relaxation process is equal to ∼75 (2.266 mm/s) at σ 0 zx = 20 (∼10 Pa) (see Figure 5a) and is ∼95 (2.871 mm/s) at σ 0 zx = 30 (∼15 Pa) (see Figure 5b). In the case when the heat flux q 0 = 0.02 (Q 0 ∼ 200 nW/µm 2 ) across the upper boundary is directed to the bulk of the tumbling type nematic phase composed of 8CB molecules whereas the SS σ zx is applied to the upper restricted surface, the relaxation of the temperature field χ(z, τ) to its stationary distribution χ st (z) across the HAN channel is characterized by an almost linear dependence χ(z, τ) from the temperature at the lower boundary χ z=0 = 0.98 (∼307 K) to the temperature at the upper boundary χ z=1 = χ up (see Figure 6). The effects of SS σ 0 zx directed in the negative direction both on the evolution of director fieldn to its stationary orientationn st in the microsized HAN channel composed of 8CB molecules, which is described by the polar angle θ(z, τ k ) (see Figure 7) and the velocity field u(z, τ k ) (see Figure 8) for different times starting from τ 1 = 0.001 (curve 1) to τ R = τ 7 ∼ 0.4 (∼0.07 s) (curve 7), are shown in Figures 7 and 8, respectively.  Fig.8(a)), and ∼ 102 (∼ 3.14 mm/ at σ 0 zx = −30 (∼ 15 P a) (see Fig.8(b)), respectively. In the case when the heat flux acro the upper surface is restricted (q 0 = 0.02 (Q 0 ∼ 200 nW/µm 2 )), we are dealing with t almost linear increase of χ(z, τ ) across the HAN channel from the temperature at the low  Figure 8. (a) The relaxation of the velocity field u(z, τ k ) to its stationary distribution across the HAN microfluidic channel composed of 8CB molecules and under the effect of SS σ 0 zx = −20 (∼10 Pa) for different times starting from τ 1 = 0.001 (curve 1) to τ R = τ 7 ∼ 0.32 (∼ 0.053 s) (curve 7). Here, τ k = 0.053 × (k − 1) (k = 2, ..., 7). (b) The same as in (a) but with SS σ 0 zx = −30 (∼15 Pa) and τ R = τ 7 ∼ 0.6 (∼0.1 s). Here, τ k = 0.1 × (k − 1) (k = 2, ..., 7), and q 0 is equal to 0.02. According to our calculations, SS σ 0 zx produces the velocity field u(z, τ) directed in the negative direction, and its effect on the director distribution across the HAN microfluidic channel is so strong that, in the middle part of the nematic channel, the director fieldn is directed almost orthogonally to both boundaries (the biggest value of the polar angle is 3.14 (∼180 • ) (see Figure 7b). The relaxation process of the velocity field is characterized by the growth of |u(z, τ)| upon increasing τ, before achieving the stationary distribution u st (z) = u(z, τ = τ 7 = τ R ) across the microsized HAN channel. This distribution is characterized by the maximum value of u st (z = 1) on the upper bounding surface (z = 1), and the hydrodynamic flow u st (z = 1) is directed parallel to both bounding surfaces in the negative direction. The maximum value of the dimensionless velocity |u st (z = 1)| = γ 10 d K 10 |v st z (z = 1)| in the HAN channel on the upper bounding surface at the final stage of the relaxation process is equal to ∼60.4 (∼1.9 mm/s) at σ 0 zx = −20 (∼10 Pa) (see Figure 8a) and ∼102 (∼3.14 mm/s) at σ 0 zx = −30 (∼15 Pa) (see Figure 8b). In the case when the heat flux across the upper surface is restricted (q 0 = 0.02 (Q 0 ∼ 200 nW/µm 2 )), we deal with the almost linear increase of χ(z, τ) across the HAN channel from the temperature at the lower (χ z=0 = 0.98 (∼307 K)) to the value at the upper boundary χ(z = 1). The relaxation of the dimensionless temperature at the upper boundary of the HAN microfluidic channel χ z=1 (τ) consisting of 8CB molecules to its stationary value χ st z=1 for three values of SS, σ 0 zx = −10 (curve 1), −20 (curve 2), and −30 (curve 3), is shown in Figure 9.
Thus, the highest temperature difference ∆χ = 0.0145 (∼4.3 K), which initially was equal to zero, is built in the HAN channel under the influence of SS σ 0 zx = −20 (∼−10 Pa). Note that, in all these cases, the dimensionless temperature at the lower boundary is kept constant χ z=0 = 0.98 (∼307 K), and the vertical temperature gradient ∇χ is created over the HAN channel, directed towards the warmer upper boundary.
The effect of SS σ 0 zx , applied both in the positive 10 (∼5 Pa) (see Figure 10a) (case I) and negative −10 (∼−5 Pa) (see Figure 10b) (case II) directions on the evolution of the director fieldn to its stationary orientationn st in the microsized HAN channel, composed of laminar type nematic (5CB), is shown in Figure 10. This evolution is described by the polar angle θ(z, τ k ), and the calculations are given for different times starting from τ 1 = 0.001 (curve 1) to τ R = τ 7 ∼ 0.47 (∼0.08 s) (curve 7).
First, the effect of SS on the evolution of the director fieldn is manifested in the qualitative difference in the polar angle profiles for cases I (see Figure 10a) and II (see Figure 10b). In case I, we have convex profiles, when the polar angle θ(z, τ k )(5CB) increases monotonically from 0 to ∼ 1.57 (π/2), whereas in case II, the polar angle θ(z, τ k )(5CB) decreases monotonically from 0 to θ(z = 0.3, τ k )(5CB) ∼ −0.28, with a subsequent increase in the value of ∼1.57 (π/2) at the upper boundary of the HAN channel.
Second, the effect of SS applied both in the positive (case I) and negative (case II) directions on the evolution of the velocity field u(z, τ) is mainly quantitative (see Figure 11a,b), where the velocity u(z, τ k ) increases from zero (u(z = 0, τ k ) = 0) at the lower boundary of the channel to the value u(z = 1, τ k ) ∼ 22 (∼0.7 mm/s) at the upper boundary in case I and from zero (u(z = 0, τ k ) = 0) at the lower boundary of the channel to the value u(z = 1, τ k ) ∼ −10 (∼−0.32 mm/s) at the upper boundary in case II. to zero, is built in the HAN channel under the influence of SS σ 0 zx = −20 (∼ −10 P a) that in all these cases, the dimensionless temperature at the lower boundary is kept co χ z=0 = 0.98 (∼ 307 K), and the vertical temperature gradient ∇χ is created over the channel, directed towards the warmer upper boundary.
The effect of SS σ 0 zx , applied both in positive 10 (∼ 5 P a) (see Fig.10(a)) (case negative −10 (∼ −5 P a) (see Fig.10(b)) (case II) directions, on the evolution of the d fieldn to its stationary orientationn st in the microsized HAN channel, composed of la type nematic (5CB), is shown in Fig.10. This evolution is described by the polar θ(z, τ k ), and calculations are given for different times started from τ 1 = 0.001 (curv τ R = τ 7 ∼ 0.47 (∼ 0.08 s) (curve 7), respectively. First of all, the effect of SS on the evo of the director fieldn is manifested in the qualitative difference in the polar angle profi the cases I (see Fig.10(a)) and II (see Fig.10(b)), respectively. In the case I, we have profiles, when the polar angle θ(z, τ k )(5CB) increases monotonically from 0 to ∼ 1.57 whereas in the case II, the polar angle θ(z, τ k )(5CB) decreases monotonically from (b) Figure 11. The evolution of the velocity field u(z, τ k ) to its stationary distribution across the HAN microfluidic channel composed of 5CB molecules and under the effect of two values of SS σ 0 zx : (a) the first is equal to 10 (∼5 Pa), whereas (b) the second is equal to −10 (− ∼−5 Pa), respectively. Different times the same as in Figure 10.

A Role of a Flow in a Temperature Gradient Formation across a HAN Channel
The purpose of this paragraph is to show the simple way the temperature gradient can be built up across the HAN channel under the action of the hydrodynamic flow in the framework of the classical Ericksen-Leslie theory [10,11] while accounting for the entropy balance equation [12]. We consider the thermal conductivity regime, which assumes that the temperature at the lower boundary is kept constant whereas, at the upper boundary, where it was assumed that the heat flux is vanished, it must satisfy the boundary conditions As a result, the temperature difference, being initially equal to zero, grows by up to the maximum possible value ∆χ = χ up − χ lw , corresponding to the nematic phase. The answer to the question of which restricted surfaces are cooler or warmer depends on the direction of the hydrodynamic flow v = u(z)î.
We consider a nematic system consisting of asymmetric polar molecules, such as cyanobiphenyl, which are confined between two solid surfaces that impose a preferred orientation of the average molecular directionn on the restricted surfaces, for instance, homeotropic on the lower boundary surfaces and planar on the upper bounding surfaces. Therefore, we describe the HAN channel under the influence of the temperature gradient ∇χ parallel to the unit vectork. Here,k is a unit vector directed from the lower substrate to the upper one (see Figure 1). The coordinate system defined by this task assumes that the directorn(r, t) lies in the xz plane (or in the yz plane) (see Figure 1).
Assuming that the temperature gradient is ∇χ = ∂χ(z,τ) ∂zk due to the growth of the temperature difference at the HAN boundaries under the action of the hydrodynamic flow changes only in the direction z, we can assume that the components of the director n = sin θ(z, τ)î + cos θ(z, τ)k as well as other physical quantities depend only on the z coordinate. Here, θ denotes the polar angle, i.e., the angle between the direction of the directorn and the normalk to the bounding surfaces.
Assuming that we deal with an incompressible fluid, the dimensionless hydrodynamic equations corresponding to the torque balance (see Equation (9)) and the linear moment balance (see Equation (10)) equations as well as the entropy balance equation (see Equation (12)) take the following form [23]: where A(θ) = 1 2 (1 − γ 21 cos 2θ) and G(θ) = sin 2 θ + K 31 cos 2 θ are the hydrodynamic and elastic functions, respectively; σ zx = δR δu ,z is the ST component; the set of the LC parameters δ i (i = 6,7,8,9) is the same as in Section 3; τ = K 1 γ 1 d 2 t is the dimensionless time; and z = z/d is the dimensionless distance away from the lower boundary of the HAN channel. Note that the overbars in the space variable z in the last three Equations (21)-(23) have been eliminated. Now, consider that the HAN system confined between two solid surfaces when the directorn is strongly anchored to these boundaries, homeotropically to the lower and homogeneously to the upper boundaries, whereas the velocity on these boundaries must satisfy the no-slip boundary condition, Now, the temperature field χ(z, τ) in the HAN channel confined between two solid boundaries when the temperature at the lower boundary is kept constant but where, at the upper boundary, it was assumed that the heat flux vanished must satisfy the boundary conditions [22] χ(z) z=0 = χ 1 , (χ ,z (z)) z=1 = 0.
The set of parameters that are involved in Equations (21)-(23) are equal to δ 6 ∼ 24, δ 7 ∼ 2 × 10 −6 , δ 8 ∼ 6 × 10 −4 , and δ 9 ∼ 2 × 10 −9 . Using the fact that δ 7 , δ 8 , and δ 9 1, the linear moment balance (22) and the entropy balance (23) equations can be considerably simplified. Thus, the whole left-hand side of Equations (22) and (23) can be neglected, and these equations take the following form: where γ 1 h(θ) = α 1 sin 2 θ cos 2 θ − A(θ)θ ,τ u ,z + 1 2 α 4 + g(θ), g(θ) = 1 2 α 6 sin 2 θ + α 5 cos 2 θ , C(τ) is the function that does not depend on z and is fixed by the boundary conditions and where To be able to observe the formation of the temperature difference across the HAN channel under the effect of the stationary hydrodynamic flow, the stationary analog of the Equation (21) was considered when θ ,τ = 0. In this case, the dimensionless temperature across the hybrid aligned nematic channel is given by [23] where H(θ, u ,z ) = h(θ)u ,z − A(θ)θ ,τ , I(θ, z) = δ 6 θ ,z sin 2 θ 1 + 1 2 sin 2 θ , and χ 1 = T 1 /T N I . The formation of the temperature difference across the HAN channel under the influence of the stationary flow with a triangular sharp profile was investigated by the standard numerical relaxation method [23], and the results are shown in Figure 12a,b. The relaxation criterion = |(θ(τ R ) − θ eq )/θ eq | for calculating procedure was chosen to be equal to 10 −4 and then, the numerical procedure was carried out until a prescribed accuracy was achieved. When the stationary hydrodynamic flow v = u(z, ζ)î is directed in the positive direction (see Figure 12a), and the temperature at the lower boundary of the HAN channel keeps a constant value χ z=0 = χ 1 = 0.97 (∼298 K), across the nematic sample, a vertical temperature gradient ∇χ directed to warmer upper boundary forms. The highest tempera-ture difference χ max (ζ) ≡ ∆χ = χ up − χ lw = 0.03 (∼9 K), which was initially equal to zero, is built in the HAN channel under the influence of the hydrodynamic flow u(z, ζ), where the magnitude of the factor α is equal to 0.0009 (∼1.2 nm/s) (see Figure 12a, curve (1)). The rest curves (2) and (3) correspond to α = 0.0007 (∼1 nm/s) and 0.0005 (∼0.7 nm/s), respectively. In the case of the reverse direction v = −u(z, ζ)î (see Figure 12b), the temperature on the upper bounded surface remains constant χ z=1 = χ 1 = 0.99 (∼304 K), while the lower surface cools to 0.96 (∼295 K), which is close to the nematic-solid phase transition temperature at α = 0.0009 (curve (1)). In all of these cases, ζ (∼0.98) is close to the upper boundary of the HAN channel. The value of χ max (ζ) has a huge influence on the location of the maximum of u(z, ζ). In the case where ζ is placed in the middle part of the nematic channel, ζ = 0.5 (see Figure 13a), the temperature difference, which was initially equal to zero, increases to ∆χ ∼ 0.0008 (∼0.3 K) at the value of α = 0.0009 and only with the increase in α by up to a two order of magnitude, from 0.0009 up to 0.09 (∼0.1 µm/s); that difference in temperature increases to a few degrees ∆χ ∼ 0.02 (∼6 K). The effect of the ζ position on the value of the maximum temperature difference χ max (ζ) when a constant temperature is kept at χ z=0 = χ 1 = 0.97 (∼298 K) at the lower boundary for a number of values of the hydrodynamic velocities u(z, ζ) is shown in Figure 13b. Notes that the velocity u(z, ζ) and the temperature χ(z) at the boundaries must satisfy the boundary conditions Equations (24) and (25), respectively. We found that the value of χ max (ζ) is sensitive to the position of ζ and shows an increase in the magnitude of the highest temperature difference when ζ is close to both boundaries of the HAN channel. This behavior of χ max (ζ) is dictated by Equation (29). Indeed, in the case of the stationary flow, the value of the shear rate shift ∆u ,z = u + ,z − u − ,z ∼ 1 ζ(1−ζ) increases to infinity in the vicinity of the bounding surfaces, when the position of ζ → 0 or 1. Physically, this means that only the stationary flow with a triangular sharp profile and a maximum position near boundaries can create the highest temperature difference in the hybrid aligned nematic channel of several degrees. With another hydrodynamic flow with profiles that cannot demonstrate the sharp growth of u ,z (z), the same result can be achieved only by using the "high-speed" hydrodynamic flow ∼0.1 µm/s. It was shown by the Brewster angular spectroscopy (BAS) method that the hydrodynamic flow of arachidic (eicosanoic) acid through the narrow channel with the width of ∼0.1 µm and with the triangular velocity profile at a shear rate of more than 0.2 s −1 at different values of surface pressure can be achieved [22]. It has also been shown that, as the flow rate increases, the velocity profile gradually becomes sharper, eventually becoming triangular. In a typical fluid, such a profile would indicate shear thickening. If this is the case, we do not exclude the possibility of extending the BAS method to the case of the abovementioned nematic system.

Conclusions
This section discusses some recent numerical advances in predicting the structural and hydrodynamic behavior of thermally excited flow in microfluidic hybrid aligned nematic (HAN) channels. Despite the fact that certain quantitative and qualitative advances have been made in the hydrodynamic description of relaxation processes in microsized nematic channels under the effect of a temperature field, there are still a number of questions concerning the temperature gradient formation across these channels. It is shown that, under the influence of the temperature gradient, the horizontal nematic layer, initially at rest, starts moving in the horizontal direction when heated both from below and from above. In the case of strong homeotropic and planar anchorings at the boundaries, the equilibrium distribution of the velocity field u st (z) over the HAN channel is characterized by a sharp increase in the absolute value of u st (z) in the vicinity of the boundary with the planar anchoring [7]. This result, in turn, leads to a number of questions. Is it possible for a temperature difference ∆T to occur between two boundaries of the HAN microfluidic channel as a result of the stationary hydrodynamic flow distribution of the flow u st (z) across the channel or by applying the SS σ zx to the boundaries of this channel? Alternatively, in more general terms, how does ∆T depend on u st (z) or σ zx ?
The second set of questions is related to the study of the influence of the microscopic width of the LC channel (confining) on the nature of the orientation dynamics of such an LC system. First, it should be noted that, in our case, the shear flow is formed under conditions of strong anchoring of the director field to the boundaries of the LC channel despite the fact that we deal with a microsized channel. Thus, confining also plays a crucial role in the formation of orientation dynamics in such an LC system. Apparently, with an increase in the width of the LC channel, the nature of the orientation dynamics of the director field changes in the direction of continuous rotation of the director, as we described in the article [19].
The absence of continuous rotation of the director field in the tumbling regime is actually due to the combined effect of both anchoring and the size of the HAN channel. Undoubtedly, further investigation of changes in the orientation dynamics in a microscopic LC channel under the influence of a shear stress applied to one of the boundaries of the LC channel as the channel width increases will be continued.
The above numerical results show that both the stationary flow with a triangular velocity profile u st (z) and SS σ zx applied to the boundaries of the HAN channel can, under certain conditions, overcome the viscous, elastic, thermomechanical, and anchoring forces and cause a temperature gradient across this channel, with the maximum absolute temperature difference being up to several degrees.
This once again shows that the macroscopic description of the nature of the hydrodynamic flow of an anisotropic liquid subtly senses the microscopic structure of an LC material.
We believe that the present investigation can shed some light on the problem of precise handling of microvolume LC drops, which requires self-contained micropumps.

Conflicts of Interest:
The authors declare no conflict of interest.