Heat Driven Flows in Microsized Nematic Volumes: Computational Studies and Analysis

: The nematic ﬂuid pumping mechanism responsible for the heat driven ﬂow in microﬂuidic nematic channels and capillaries is described in a number of applications. This heat driven ﬂow can be generated either by a laser beam focused inside the nematic microvolume and at the nematic channel boundary, or by inhomogeneous heating of the nematic channel or capillary boundaries. As an example, the scenario of the vortex ﬂow excitation in microsized nematic volume, under the inﬂuence of a temperature gradient caused by the heat ﬂux through the bounding surface of the channel, is described. In order to clarify the role of heat ﬂux in the formation of the vortex ﬂow in microsized nematic volume, a number of hydrodynamic regimes based on a nonlinear extension of the Ericksen–Leslie theory, supplemented by thermomechanical correction of the shear stress and Rayleigh dissipation function, as well as taking into account the entropy balance equation, are analyzed. It is shown that the features of the vortex ﬂow are affected not only by the power of the laser radiation, but also by the duration of the energy injection into the microsized nematic channel.


Introduction
Impact of liquid crystal (LC) materials on modern technology is very impressive [1]. The recent technological revolution has been brought by these LC materials in the field of displays. By now, for instance, flat-panel LC displays (LCDs) are widely used in consumer electronics, personal computers, and mobile devices, as well as in many types of medical, transportation, and industrial equipment. With the further development of the LC-display market, the question then arises of what are the next areas of the LCs application? Perhaps there is no more appropriate directions for application of the LC materials than the LC sensors and actuators [2]. 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 due to the fact that liquid crystals are extremely sensitive to external stimuli and therefore can be used for the construction of stimuli-responsive devices, such as sensors or actuators. There is another suitable direction for application of the LC materials, which is related to the precise manipulation of complex liquids on a microscale. The problem of electrically driven manipulation of complex liquids such as LCs has brought an increasing number of integrated small-scaled microdevices for chemical and biological applications. These are known as micrototal analysis systems or biological lab-on-a-chip systems [3][4][5]. Electro-osmosis, dielectrophoresis, and electrowetting have been explored for controlling microflows [3,4,6]. Thus, understanding the flow of LCs through micron-scale channels is of increasing significance. Nematic drops of appropriate size confined in the capillary are microdevices, whose molecular orientations also can be manipulated by a temperature gradient ∇T, generated, for instance, by a laser beam focused both in the microfluidic volume or on the boundary of the LC channel. Fluid pumping principle has been developed utilizing the interaction, on the one hand, between the gradient of director field ∇n and a temperature gradient ∇T, excited, for instance, by a laser beam, focused both on the boundary of the LC channel or inside the LC microvolume, and, on the other hand, between the ∇n and the velocity field v, excited in the microfluidic channel by the laser irradiation [7]. The problem of motion of an ultra-small (a few microliters) isotropic and LC drops, under the influence of the temperature gradient, caused by a laser beam, has drawn considerable attention [8][9][10][11][12][13][14]. Understanding of how the liquid crystal material can be manipulated under the influence of the temperature gradient is also a matter of great fundamental interest, as well as an essential part of knowledge in the field of soft materials science.
Despite the fact that the possibility of the formation of hydrodynamic flows in nematic channels under the influence of temperature gradients has been theoretically described since [15], only detailed numerical simulations performed within the framework of the extended Ericksen-Leslie theory [16,17] allowed us to recreate the complete picture of the formation of flows in nematic microchannels and capillaries. One of the aims of this review 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 review is devoted to the latest results describing the possibilities of computational methods implemented in the framework of the nonlinear extension of the Ericksen-Leslie theory, with accounting the entropy balance equation [18]. Another purpose of our review is to show some routes in describing the laser-excited motion of nematics enclosed in a microsized channel with a free surface. This is the first such review that describes in detail the role of thermomechanical forces in the formation of hydrodynamic flow in microsized nematic channels and capillaries. It is based on the nonlinear extension of the Ericksen-Leslie theory, supplemented by thermomechanical correction of the shear stress and Rayleigh dissipation function, and also takes into account the entropy balance equation. The fact that the main results were obtained by numerical methods indicates that experimenters still have a lot of work to do in order to create a more complete picture of the formation of hydrodynamic flows in microsized nematic channels and capillaries.
The layout of this article is as follows. In the next section, we give the review of theoretical background to the description of the physical mechanism responsible for the heat driven flow in the microfluidic rectangular nematic channel. The heat driven nematic flow in a cylindrical microfluidic channel is described in Section 3. Our conclusions are given in Section 4.

Heat Driven Nematic Flow in Rectangular Microfluidic Channel
We are primarily concerned with the description of the physical mechanism responsible for the heat driven flow in the microfluidic rectangular nematic channel. This heat driven flow v, generated by the laser beam focused both on the restricted surface or inside the nematic microvolume can set up the temperature gradient [9,12,13]. A challenging problem in all such LC systems is the precise handling of nematic microvolume which requires self-contained micropumps of small volume flow (dynamic pumps). This the nematic pumping principle is based on the coupling between the velocity and director fields, together with accounting the effect of the temperature gradient. This problem will be treated in the framework of the extended Ericksen-Leslie theory [16,17], supplemented by the thermomechanical correction of shear stress [7,15], and the entropy balance [18] equation.
If a liquid crystal drop placed between two horizontal plates is at rest and if it is heated both from below or above, then due to the coupling of the director and temperature gradients, a thermomechanical force develops that can overcome the viscous, elastic, and anchoring forces, and the LC drop begins to move in the horizontal plane. Such horizontal motion of the LC drop due to uniform heating from below was first studied by Lehmann [19].
In the past decade, some research activity has been focused on the problem of formation the hydrodynamic flows in the microsized LC volumes excited by the thermomechanical force σ tm zx [7,[20][21][22][23]. In the case when the temperature gradient ∇T in the hybrid aligned nematic (HAN) volume is proportional to ∆T/d, the hydrodynamic flow v excited by the temperature gradient is proportional to v ∼ d η σ tm zx , where ∆T = T 2 − T 1 > 0 (the range [T 2 , T 1 ] falls within the stability region of the nematic phase) is the temperature difference on the LC volume boundaries, η is the viscosity, and σ tm zx ∼ ξ ∆T d 2 is the tangential component of the thermomechanical stress tensor σ tm ij , and ξ is the thermomechanical constant. Any physical effect that induces flow in the LC volume leads to the coupling of that flow with the director fieldn. While having a subordinate role, the hydrodynamic flow v has been found to qualitatively change the orientational behavior of the directorn in temperature driven flow.
Therefore, we are primary concerned here with the description of the physical mechanism responsible for the flow in the LC channel confined between solid surfaces, which is excited by the temperature gradient ∇T, and the magnitude of that flow is proportional to ∆T, whereas the direction of v influences both the direction of the heat flux q and the character of the preferred anchoring of the average molecular directionn to the restricted surfaces. We focused on the analysis of the problem of formation of hydrodynamic flows in microsized nematic channels under the influence of directed heat fluxes from the position of numerical methods. The analysis of the dynamics of Lehmann-type effects in chiral LCs is far from our research interests. For more information see, for instance, ref. [22].
In an attempt to better understand the dissipation processes in microfluidic nematic systems confined in the various microsized geometries, under the influences the temperature gradient, we will review a number of numerical studies of these nematic systems. With this aim, a number of theories which include the hydrodynamic equations describing both the director reorientation and flow of nematic fluid, as well as the redistribution of the temperature field, will be described. These equations also will be supplemented by appropriate anchoring conditions for the director field on the bounding surfaces, as well as the no-slip condition for the velocity field. Since the geometry of the microsized nematic channel or capillary affects the nature of the hydrodynamic flow that is formed, two cases will be considered. The first is when dealing with a rectangular channel, and second is the cylindrical capillary.

Formulation of the Balance of the Momentum and Torque Equations and Conductivity Equation for Nematic Fluids Confined in Rectangular Channel
The aim of this section is to show some useful routes for further examining of the validity of theoretical treatment of the orientational dynamics in the microsized rectangular nematic channel under the effect of externally directed heat flux q across the bounding surface. It will be done in the framework of the extended Ericksen-Leslie theory [16,17], supplemented by the thermomechanical correction of shear stress [7,15], and the entropy balance equation [18]. The first part of the section is devoted to the hydrodynamic model, which is the basis for describing the formation of flows in microsized hybrid aligned nematic (HAN) channels and capillaries under the action of the thermomechanical force. This force is the result of interaction of the directorn and temperature T gradients, and it is responsible for the formation of the hydrodynamic flow v in these channels. In our case, the temperature gradient ∇T in the HAN channels is formed by focused laser radiation. In the second part of the section, a numerical analysis of two modes of heating the LC channels is given. These modes are characterized by different laser radiation power and duration: slow and fast heating modes caused by the laser irradiation focused on the lower boundary of the HAN channel.
With this aim we consider the microsized HAN channel delimited by two lower and upper horizontal solid surfaces, located at z = 0 and z = d, respectively, and two lateral solid surfaces at distance L on scale on the order of micrometers. The coordinate system defined by this task assumes that the directorn = (n x , 0, n z ) = sin θî + cos θk, where θ ≡ θ(t, x, z) denotes the polar angle, i.e., the angle between the direction of the directorn and the normalk to the horizontal boundary surfaces, is in the XZ plane. Here,î is the unit vector directed parallel to the horizontal boundaries of the HAN channel, andĵ =k ×î (see Figure 1). Figure 1. The coordinate system used for theoretical analysis. The x-axis is taken as being parallel to the director directions on the upper surface, θ(z, t) is the angle between the directorn and the unit vectork, respectively. Both the heat flux q and the unit vectork are directed normal to the horizontal boundaries of the liquid crystal (LC) channel.
We shall be considering the HAN with the heat flux q = qk across the upper boundary whereas on the remaining boundaries the temperature is kept constant Here λ ⊥ is the heat conductivity coefficient perpendicular to the directorn. Therefore, the hybrid aligned nematic state, where at the upper boundary of the channel, the nematic phase is anchored homogeneously and the bulk of the nematic sample contains a gradient in θ(z) = π 2d z from homeotropic orientation at the lower surface to planar orientation at the upper surface, i.e., or θ 0<x<L,z=0 = θ x=0,0<z<d = θ x=L,0<z<d = 0, Taking into account that the width of the nematic channel d ∼ 1 − 5 µm, one can assume that the mass density ρ m = const across the HAN channel, and one deals with an incompressible fluid. In turn, the incompressibility condition ∇ · v = 0 assumes that where u ≡ v x (t, x, z) and w ≡ v z (t, x, z) are the components of the vector v = uî + wk, and u ,α = ∂u ∂α (α = x, z). The hydrodynamic equations describing the reorientation of the nematic phase in 2D case, when the heat flux across the upper restricted surface exists, whereas the temperature on the rest surfaces is kept constant, can be derived from the balance of elastic, viscous, and thermomechanical torques T el + T vis + T tm = 0, the Navier-Stokes equation for the velocity field v, excited by the heat flux q, and the equation for the heat conduction.
To be able to observe the formation of the hydrodynamic flow v in the HAN channel, under the effect of the heat flux q, when the heating occurs during some times τ in , let us consider the dimensionless analog of the torque balance equation [7,16,17,21] n z n x,τ − n x n z, where τ = t/t T is the dimensionless time, t T = ρC p d 2 /λ ⊥ , C p is the heat capacity, γ = γ 2 /γ 1 , γ 1 and γ 2 are the rotational viscosity coefficients, A 0 = n x,x + n z,z , χ ≡ χ(τ, x, z) = T(τ, x, z)/T N I is the dimensionless temperature, T N I is the nematic-isotropic transition temperature, χ ,α = ∂χ ∂α (α = x, z), z = z d is the dimensionless distance away from the lower boundary of the HAN channel, x = x L is the dimensionless space variable corresponding to x-axis, respectively,ψ = t T d 2 ψ is the scaled analog of the stream function ψ for the velocity field v = uî + wk = −∇ ×ĵψ (see the ref. [21]), n z,τ = ∂n z ∂τ , N x = n x n z,x − n z n x,x , N z = n z n x,z − n x n z,z , L x = n x n z,x − 3 2 n z n x,x + 1 2 n x n x,z , and L z = −n z n x,z + 3 2 n x n z,z − 1 2 n z n z,x , respectively. Here, δ 1 = t T K 1 γ 1 d 2 and δ 2 = are two parameters of the nematic system, and ξ is the thermomechanical constant [7,15].
The dimensionless Navier-Stokes equation for the velocity field v = uî + wk = −∇ ×ĵψ takes the form [21] δ 3 ψ ,xzτ = a 1 ψ ,zzzz + a 2 ψ ,xzzz + a 3 ψ ,xxzz + a 4 ψ ,xxxz + a 5 ψ ,xxxx + a 6 ψ ,zzz + a 7 ψ ,xzz + a 8 ψ ,xxz + a 9 ψ ,xxx + a 10 ψ ,zz + a 11 ψ ,xz + a 12 ψ ,xx + F , whereas the dimensionless entropy balance can be written as [21] χ ,τ = χ ,x λn 2 where λ = λ /λ ⊥ is the dimensionless parameter, λ is the heat conductivity coefficient parallel to the director. It should be noted that the function , the coefficients a i (i = 1, . . . , 12), the functions σ tm ij (i, j = x, z) and σ el ij (i, j = x, z) are given in the ref. [21], whereas σ ij = P δ ij + σ el ij + σ vis ij + σ tm ij (i, j = x, z) is the stress tensor (ST) of the nematic system, P is the hydrostatic pressure in the HAN channel, and σ el ij , σ vis ij , and σ tm ij are the dimensionless ST components corresponding to the elastic, viscous, and thermomechanical forces, respectively. The set of the rest parameters of the nematic system, corresponding to the case of 4 − cyano − 4 − pentylbiphenyl (5CB), are: , respectively. In the following we are focused primary on the heat conduction regime in the HAN channel which assumes that across the upper surface the heat flux is set up (see Equation (1)), whereas on the rest surfaces the temperature is kept constant (see Equation (2)). Physically, this means that across the LC channel the temperature gradient ∇T, directed from the cooler to the warmer surfaces, excited by the heat flux q across to the upper boundary, may be built up.
Notice that the dimensionless ST σ ij can be obtained directly from the elastic contribution to the energy and Rayleigh dissipation function as σ el = − ∂F el ∂∇n · (∇n) T , σ vis = δR vis δ∇v , and σ tm = δR tm δ∇v , for the elastic, viscous, and thermomechanical contributions, respectively. Straightforward calculations give the following expressions for ST components σ el ij , σ vis ij , and σ tm ij which are listed in the ref. [21]. Now the evolution of the director in the HAN channel confined between two horizontal and two vertical solid surfaces, under the influence of viscous, elastic, and thermomechanical forces and taking into account the flow, can be obtained by solving the system of the nonlinear partial differential Equations (6)-(8) with the appropriate dimensionless boundary conditions for the directorn or the polar angle θ temperature field χ(τ, x, z) and the initial condition taken in the form respectively. Here, χ 0 = T 0 T N I , and Q = − d q is the dimensionless heat flux across the upper restricted surface. The laser-induced heating was used to inject energy Q across the bounding surface at the microscopic scale [8,13]. In that case the nematic channel was heated by a laser beam, focused, for instance, on the upper bounding surface (z = 1, x = x 0 ), where P 0 is the laser power, and ω 0 is the Gaussian spot size. Taking into account that the total absorbed laser power is P in = αI (x), the heat flux across the upper restricted surface can be written as where α is the absorption coefficient, is the Heaviside step function, and τ in is the duration of the energy injection into the nematic sample. Here, and everywhere else in this section, Q 0 is equal to |Q 0 |. For the case of 4 − cyano − 4 − pentylbiphenyl (5CB), at temperature corresponding to nematic phase, the first four parameters δ 1 = that are involved in Equations (6)-(8) have the following values [21]: Using the fact that δ 3 1, the Equation (7) can be considerably simplified and takes the form [21] a 1 ψ ,zzzz + a 2 ψ ,xzzz + a 3 ψ ,xxzz + a 4 ψ ,xxxz + a 5 ψ ,xxxx + a 6 ψ ,zzz + a 7 ψ ,xzz + a 8 ψ ,xxz + a 9 ψ ,xxx + a 10 ψ ,zz + a 11 ψ ,xz + a 12 ψ ,zz + F = 0, where a i (i = 1, . . . , 12) and F are functions which have been defined in ref. [21].

Orientational Relaxation of the Director, Velocity, and Temperature Fields in Rectangular Microsized Nematic Channels
With the evolution of the directorn to its equilibrium orientationn eq , which is described by the polar angle θ(τ, x, z) from the initial θ(τ = 0, x, z) = θ elast (x, z) to the equilibrium θ eq (x, z) conditions, both the horizontal and vertical components of the velocity field v = uî + wk = −∇ ×ĵψ, and the temperature field χ(τ, x, z) redistribution, can be obtained by solving the system of the nonlinear partial differential Equations (6), (8) and (14), together with the boundary conditions (9)- (11), and the initial condition taken in the form of (12). The calculations have been carried out by using both the relaxation [24] and the sweep [25] methods. The initial distribution of the director fieldn el (x, z) has been obtained from Equation (6) by means of the relaxation method with ψ ,x = ψ ,z = (χ) ,x = (χ) ,z = 0, and with the boundary (n x ) x=0,0≤z≤1 = (n x ) x=1,0≤z≤1 = (n x ) z=1,0≤x≤1 = 0, (n x ) z=0,0≤x≤1 = 1, and initial n x = 1−z 2 , for 0 < x < 1, and n x = 0, for x = 0, 1, conditions. Having obtained the initial distribution of the director fieldn el (x, z), the initial distribution of the temperature field χ(∆τ, x, z), corresponding to the first time step ∆τ, has been obtained from Equation (8), by means of the relaxation method with ψ ,x = ψ ,z = 0 and with the boundary and initial conditions in the form of Equations (9)- (12). Then the initial distribution of the stream function ψ(∆τ, x, z), corresponding to the first time step ∆τ can be calculated by using the initial distributions of the directorn el (x, z) and temperature χ(∆τ, x, z) fields, as well as the function F , which is involved in Equation (14). The next time step ∆τ for the velocity and temperature fields, as well as for the director s distribution across the nematic sample is initiated by the sweep method. The stability of the numerical procedure for Equations (6), (8) and (14) was defined by the following conditions [24]: where ∆x and ∆z are the space steps in the x and z directions, and a 1 and a 5 are the coefficients defined in Equation (14). In the calculations, the relaxation criterion = |(χ (m+1) (τ, x, z) − χ (m) (τ, x, z))/χ (m) (τ, x, z)| was chosen to be equal to 10 −4 , and the numerical procedure was then carried out until a prescribed accuracy has been achieved.
Here, m is the iteration number and τ R is the relaxation time of the system.

The maximum of both the absolute magnitudes of the dimensional velocities
in the rectangular HAN channel, at the initial stage of the slow heating mode are equal to 1.2 mm/s and 0.88 mm/s, for the horizontal and vertical components (see Figure 5a and Figure 7a), at ∆χ max = 0.03 (∼9.3 K), respectively.
The distance z dependence of the temperature field χ(τ, x = 0.5, z) in the middle part of the dimensionless HAN channel (x = 0.5), during its evolution to the equilibrium distribution χ eq (x = 0.5, z), under the effect of the dimensionless heat flux Q = 0.44 (∼3.7 × 10 −3 mW/µm 2 ), caused by the laser beam, at different times, from τ 1 = 0.00042 [curve (1)] to τ 6 = τ in = 0.0134 [curve (6)], respectively, is shown in Figure 8a.  Figure 8. (a) The distance z dependence of the temperature field χ(τ, x = 0.5, z) during its evolution to the equilibrium distribution χ eq (x = 0.5, z) = χ(τ = τ 6 , x = 0.5, z) in the middle part of the HAN channel [21], during the slow heating mode at different times, from The same as in (a), but for the cooling mode at different times, from τ 7 = 0.014 [curve (7)] to τ 11 = τ R = 0.22 [curve (11)], respectively. Figure 8b shows the same as in Figure 8a, but for the cooling mode at different times, from τ 7 = 0.014 [curve (7)] to τ 11 = τ R = 0.22 [curve (11)], respectively. The distance z dependence of the temperature field χ(τ, x = 0.5, z) is characterized by the temperature growth on the upper boundary of the HAN channel (z = 1), from χ z=1 = 0.97 (∼298 K) to χ z=1 = 1.0 (∼307.3 K), during the slow heating mode (curves, from (1) to (6)) (see Figure 8a), with the following temperature decrease, from χ z=1 = 1.0 (∼307.3 K) to χ z=1 = 0.97 (∼298 K), after switching off the laser power (see Figure 8b). Here, the value of τ in = τ 6 is equal to ∼0.0134 (∼1.9 ms). The second part of the relaxation process of χ(τ, x = 0.5, z) (τ ≥ τ 6 ) is characterized by much faster cooling down of the upper restricted surface than the rest bulk of the nematic sample (see Figure 8b). The growth of the temperature field χ(τ, x, z = 0.97) to its equilibrium distribution χ eq (x, z = 0.97) along the width of the HAN channel (0 ≤ x ≤ 1) in the vicinity of the upper warmer boundary (z = 0.97), during the slow heating mode at different times, from τ 1 = 0.00042 [curve (1)] to τ 6 = τ in = 0.0134 [curve (6)], respectively, is shown in Figure 9a. Figure 9b shows the same as in Figure 9a, but for the cooling mode at different times, from τ 7 = 0.014 [curve (7)] to τ 11 = τ R = 0.22 [curve (11)], respectively. This slow heating mode is characterized by symmetric growth of the profile χ(τ, x, z = 0.97) with respect to the middle part (x = 0.5) of the HAN channel. Such symmetric dependence of the temperature field along the width of the channel can be explained by much faster evolution of χ(τ, x, z) to χ eq (x, z), than the evolution both of the director and velocity fields to their equilibrium distributions across the HAN channel, and, as a result, the weak effect of these fields on χ [3].  Figure 9. The same as in Figure 8a,b, but the distance x dependence of the temperature field χ(τ, x, z = 0.97) during its evolution to the equilibrium distribution χ eq (x, z = 0.97) = χ(τ = τ 6 , x, z = 0.97) along the width of the HAN channel (0 ≤ x ≤ 1), in the vicinity of the upper warmer restricted surface (z = 0.97), both during the slow heating (a) and cooling (b) modes [21].

B. Fast heating mode
The equilibrium distributions of the polar angle θ eq (x, z) = θ(τ = τ 5 , x, z), both across the middle part θ eq (x = 0.5, z) = θ(τ = τ 5 , x = 0.5, z) and in the vicinity of the upper boundary of the HAN channel θ eq (x, z = 0.97) = θ(τ = τ 5 , x, z = 0.97), under the effect of the dimensionless heat flux Q = 3.54 (∼2.95 × 10 −2 mW/µm 2 ) (which is 8 times greater than Q in the first case A), caused by the laser beam focused on the upper boundary, are shown in Figure 10a,b, respectively. Here, the duration τ in of the energy injection into the HAN channel across the upper boundary is equal to 0.0006 (∼90 µs), and during that time the temperature on the upper boundary grows from Figure 10a), with the following cooling down to χ z=1 (τ = τ R = 0.064) = 0.97 (298 K). Note that the duration of the laser pulse, at fixed power P 0 = 115 mW, is limited by condition that the higher temperature on the upper boundary of the HAN channel χ z=1 must fall within the nematic stability range. The growth of the temperature field χ(τ, x = 0.5, z) in the middle part of the dimensionless HAN channel (x = 0.5) to its equilibrium distribution χ eq (x = 0.5, z), under the effect of the dimensionless heat flux Q = 3.54 (∼2.95 × 10 −2 mW/µm 2 ), caused by the laser beam at different times, from (1)] to τ 5 = 0.0006 [curve (5)], respectively, corresponding to the fast heating mode, is shown in Figure 11a.  Figure 11. (a) The evolution of the temperature field χ(τ, in the middle part of the HAN channel, under the effect of the dimensionless heat flux Q = 3.54 (∼2.95 × 10 −2 mW/µm 2 ) [21]. The fast heating mode is characterized by the sequence of times, from τ 1 = 0.00004 [curve (1)] to τ 5 = 0.0006 [curve (5)], respectively, whereas the cooling mode (b) is characterized by the sequence of times, from Here, τ in = τ 5 ∼ 0.0006 (∼90 µs) is the duration of the energy injection into the HAN channel across the upper boundary by the infrared laser with the power P 0 = 115 mW. The sequence of times, from τ 6 = 0.001 [curve (6)] to τ 12 = τ R = 0.064 [curve (12)], corresponding to the cooling mode, is shown in Figure 11b. The distance x dependence of the temperature field χ(τ, x, z = 0.97) during its evolution to the equilibrium distribution θ eq (x, z = 0.97) = θ(τ = τ 6 , x, z = 0.97) [21], along the width of the HAN channel (0 ≤ x ≤ 1) in the vicinity of the upper warmer boundary (z = 0.97), during the fast heating mode at different times, from τ 1 = 0.00004 [curve (1)] to τ 5 = 0.0006 [curve (5)], respectively, is shown in Figure 12a, whereas the cooling mode with the sequence of times, from τ 6 = 0.001 [curve (6)] to τ 12 = τ R = 0.064 [curve (12)], respectively, is shown in Figure 12b.  Figure 12. (a) The same as in Figure 11a, but the evolution of the temperature field χ(τ, x, z = 0.97) to its equilibrium distribution χ eq (x, z = 0.97) = χ(τ = τ 5 , x, z = 0.97) along the width of the HAN channel (0 ≤ x ≤ 1) in the vicinity of the upper warmer boundary (z = 0.97) [21] at different times, from τ 1 = 0.00004 [curve (1)] to τ 5 = 0.0006 [curve (5)], respectively, is shown for the fast heating mode. (b) The same as in (a), but the sequence of times, from τ 6 = 0.001 [curve (6)] to τ 12 = τ R = 0.064 [curve (12)], corresponds to the cooling mode.
The distance x dependence of the horizontal component u(τ, x, z = 0.97) of the velocity v during its evolution to the equilibrium distribution u eq (x, z = 0.97) = u(τ = τ 5 , x, z = 0.97), along the width of the HAN channel (0 ≤ x ≤ 1) in the vicinity of the upper warmer restricted surface (z = 0.97) [21], for the fast heating mode at different times, from τ 1 = 0.00004 [curve (1)] to τ 5 = 0.0006 [curve (5)], respectively, is shown in Figure 14a. Figure 14b shows the same as in Figure 14a, but the sequence of times, from τ 6 = 0.001 [curve (6)] to τ 12 = τ R = 0.064 [curve (12)], respectively, corresponds to the cooling mode. Both the fast heating ( Figure 14a) and cooling (Figure 14b) modes are characterized by near symmetric profile of u(τ, x, z = 0.97) with respect to the middle part (x = 0.5) of the HAN channel.
In turn, the distance x dependence of the vertical component w(τ, x, z = 0.97) of the velocity v during its evolution to the equilibrium distribution w eq (x, z = 0.97) = w(τ = τ 5 , x, z = 0.97), along the width of the HAN channel (0 ≤ x ≤ 1) in the vicinity of the upper warmer boundary (z = 0.97), both during the fast heating and cooling modes, are shown in Figure 16a,b, respectively.
The fast heating mode is characterized by exciting the vertical component w(τ, x, z = 0.97) of the velocity v only in the vicinity of the vertical boundaries of the HAN channel (see Figure 16a), with the following velocity decreasing to 0, after switching off the laser power (see Figure 16b), respectively. Note that under the effect of the heat flux across the upper boundary, only approximately 15% of the HAN channel (see Figure 15a) close to the upper warmer boundary is involved in the moving process, whereas the rest amount of the nematic sample is kept unmoved during the full relaxation term. These numerical studies show that in the hybrid aligned nematic material under the effect of q, in the vicinity of the warmer boundary, the vortical flow may be excited. The maximum for both the absolute magnitudes of the dimensional velocities v x (τ, x, z) = ( and v z (τ, x, z) = ( γ 1 d K 1 )w(τ, x, z) in the HAN channel, at the initial stage of the fast heating mode, are equal to 0.2 mm/s and 0.05 mm/s, for horizontal and vertical components (see Figure 15a and Figure 16a), at ∆χ max = 0.03 (∼9.3 K), respectively.  (6)] to τ 12 = τ R = 0.064 [curve (12)], corresponds to the cooling mode.
In summary, we present the review of numerical results describing the evolution of not only the directorn(t, x, z) and velocity v(t, x, z) fields, but also the redistribution of the temperature field χ(t, x, z) in the 2D hybrid aligned nematic channel to their equilibrium values, under the effect of the heat flux directed normal to the upper bounding surfaces, when the rest bounding surfaces of the HAN channel are kept at constant temperature. In this case, the upper nematic boundary is aligned homogeneously and heated by an infrared laser beam, and the dynamics of heating occurs with two distinct time scales: (i) a fast time scale, with the duration of the dimensional heat flux q ∼ 2.95 × 10 −2 mW/µm 2 over time t in ∼ 90 µs, and the laser power in P 0 = 115 mW, and (ii) a slow time scale, with the duration of the dimensional heat flux q ∼ 3.7 × 10 −3 mW/µm 2 over time t in ∼ 1.9 ms, and the laser power in P 0 = 14.3 mW, respectively [21]. In both these cases, the duration of the laser pulse, at fixed power in P 0 = 14.3 mW and P 0 = 115 mW, is limited by condition that the higher temperature on the upper bounding surface must fall within the nematic stability range. These numerical studies show that in such the hybrid aligned nematic volume under the effect of the heat flux, across the upper boundary, the vortical flow may be excited only in the vicinity of the warmer boundary. After some time more than the relaxation time, for instance, of the director field in the HAN channel, the vortical flow settles down to the rest state regime, where the horizontal u and vertical w components of velocity vector are equal to zero, and the temperature field across the nematic channel is finally falling to the value of temperature on the lower and both lateral bounding surfaces. Note that in the case of the fast heating mode (i) one deals with the heating which is characterized by "shallow" heating of the nematic microvolume in the vicinity of the upper bounding surface up to 40% of the full nematic sample, whereas in the case of the slow heating mode (ii), one deals with the "deeper" heating of the nematic microvolume up to 60% of the nematic sample far from the upper warmer restricted surface.
Thus, it can be concluded that for the formation of a hydrodynamic flow in a microsized HAN channel using laser radiation focused on the channel boundary, such a radiation power is necessary to pump energy into the liquid crystal phase for as long as possible. Another condition for the formation of a more powerful flow in the HAN channel is the planar anchoring condition at the boundary where the laser radiation will focus. Therefore, a new scenario of the vortical flow formation in the HAN channel under the effect of the focused laser irradiation on the lower homeotropically aligned boundary will be analyzed in the next paragraph.

Laser Excited Vortical Flow in Microsized Nematic Channels
In the previous paragraph, we analyzed the process of formation and relaxation of the velocity field v in the microsized LC volume under the effect of the heat flux q directed orthogonally to the upper boundary of the HAN channel. The purpose of this paragraph is to show how the direction of the external heat flux q affects the process of formation of the vortex flow in the microsized HAN channel. In this case, the heat flux q directed across the lower bounding surface at an angle α with respect to the unit vectorî. Here,î is the unit vector parallel to the horizontal boundaries of the HAN channel, whereas the unit vectork is directed perpendicular to both the horizontal boundaries. This problem will be treated in the framework of the extended Ericksen-Leslie theory [16,17], supplemented by the thermomechanical correction for both the shear stress and the Rayleigh dissipation function [7,15], as well as the entropy balance equation [18]. The hydrodynamic model in which the above problem will be considered is the same as in Section 2.1, and is based on the interaction effect of the directorn and temperature T gradients with the velocity v = uî + wk field. In this case, the magnitude of the hydrodynamic flow v is proportional to v ∼ d η σ tm zx , where σ tm zx is the tangential component of the thermomechanical stress tensor, η is the viscosity of the nematic phase, and d is the thickness of the HAN channel. It will be shown below that the tangential component of the thermomechanical stress, the tensor plays a crucial role in the formation of the thermally excited vortical flow in the HAN channel. In our case, the temperature gradient ∇T in the microsized nematic volume is formed by focused laser radiation. At the beginning of the section, we give a detailed description of the boundary conditions for both the directorn and temperature T(x, z, t) fields, as well as for the velocity field v. The second part of the section presents a numerical analysis of the formation of the thermally excited vortical flow in the HAN channel under the influence of the heat flux q directed across the lower bounding surface at an angle α with respect to the unit vectorî.
With this aim, we consider the 2D nematic fluid composed of polar molecules, such as cyanobiphenyls, at density ρ, and delimited by two horizontal and two lateral surfaces at mutual distances 2D and 2d, respectively. We also consider the effect of the orientational defect on the lower bounding surface on the vortex flow excited in the nematic microvolume under the effect of the focused laser beam. This defect is characterized by the continuous change in the orientation of the director s along the length of the lower restricted surface, from the homeotropic to the planar, and again to the homeotropic orientation.
Therefore, our 2D aligned nematic state with the orientational defect contains a gradient of ∇θ on the lower boundary, i.e., and the planar orientation on the upper and lateral surfaces, i.e., , L = l d , and 2l is the length of the orientational defect on the lower boundary with the director s orientation changing continuously from the homeotropic to planar, and again to homeotropic orientation, whereas on the rest length of that surface there is the homeotropic director s orientation (k n z=−1 ). Notice that in these calculations the ratio D/d is equal to 10. Moreover, we will consider dimensionless spatial variables x = x/d and z = z/d, and in the following equations the overbars will be eliminated.
In order to elucidate the role of both the orientational defect and the heat flux q directed at an angle α across the lower boundary of the HAN channel in formation of the vortical flow we consider the hydrodynamic regime with the heat flux q = q zk = Q sin αk, where α is the angle between vectors q andî. In the dimensionless form it can be written as [13,26] (χ ,z (τ, x, z)) z=−1 = q z − (λ − 1)n x n z χ ,x whereas on the rest boundaries the temperature is kept constant Here, Q = Q(x, z = −1) is the injected energy across the lower boundary, q z = Q sin α, and τ = t t R is the dimensionless characteristic time needed for reorientation of the director n, respectively.
The velocity v = uî + wk = −∇ ×ĵψ on these surfaces has to satisfy the no-slip boundary condition where u ≡ v x (τ, x, z) and w ≡ v z (τ, x, z) are the components of the vector v = uî + wk = −∇ ×ĵψ. The reorientation of the director in the HAN microvolume confined between two horizontal and two vertical solid surfaces under the effect of the viscous, elastic, and thermomechanical forces and taking into account the flow can be obtained by solving the system of the nonlinear partial differential Equations (6), (8) and (14) with the appropriate dimensionless boundary conditions for the polar angle θ (Equations (15) and (16)), or for the director fieldn velocity field v = uî + wk = −∇ ×ĵψ (see Equation (19)), temperature field and the initial condition taken in the form respectively. Here,n eq elast (x, z) is the equilibrium distribution of the director field over the nematic volume obtained from Equation (6), with u ,x = u ,z = w ,x = w ,z = χ ,x = χ ,z = 0, and with the boundary conditions in the form of Equation (20), whereas the initial condition is chosen in the form θ elast (τ = 0, x, z) = 0, at x = ±10, and −1 ≤ z ≤ 1; π 4 (z + 1), at −10 < x < −L, and L < x < 10; π 2 , at −L ≤ x ≤ L, and z = −1, respectively. Notice that the initial distribution of the director fieldn eq elast (x, z) over the microsized volume has been obtained from Equation (6) by means of relaxation method [24]. In calculations, the relaxation criterion = |(θ (m+1) (τ, x, z) − θ m (τ, x, z))/θ (m) (τ, x, z)| was chosen to be equal to 10 −4 , and the numerical procedure was then carried out until a prescribed accuracy was achieved. Here, m is the iteration number. Figure 17 shows a fragment of initial distribution of the director fieldn eq elast (x, z) near the orientational defect located at −L < x < L, z = −1 [13,26]. Here, L = ±0.5 are the right and left ends of the orientational defect on the lower boundary. According to these calculations the highest value of |∇n(x, z)| is reached in the vicinity of the orientational defect and its effect extends up to 3/4 of the thickness of the nematic volume.  Figure 17. The fragment of initial distribution of the director fieldn eq elast (x, z) near the orientational defect located at −L < x < L, z = −1. Here, L = ±0.5 [13,26].
The evolution of both the velocity v = uî + wk = −∇ ×ĵψ and the temperature χ(τ = τ in , x, z) fields over the nematic microvolume with the orientational defect located on the lower boundary of the HAN channel as the function of the angle α are shown in Figures 18-23, respectively. The distribution of the velocity field v in the microscopic nematic volume with the orientational defect located at −L ≤ x ≤ L, z = −1, when the laser beam is directed at the values of the angle α equal to 20 • and 40 • (Figure 18a,b), 60 • and 80 • (Figure 19a,b), 90 • (Figure 20), 100 • and 120 • (Figure 21a,b), 140 • and 160 • (Figure 22a,b) are shown in Figures 18-22 [13,26,27], respectively. The heating occurs during the dimensionless time τ in = 1.6 × 10 −4 (∼0.29 ms), whereas the value of dimensionless heat flux coefficient Q 0 = d T NI λ ⊥ q 0 is equal to 0.05.      The distribution of the velocity v in the microscopic HAN channel with the orientational defect located at −L ≤ x ≤ L, z = −1, when the heat flux q is directed at two values of the angle α, 20 • and 40 • with respect to the lower bounding surface, is characterized by maintaining of vortices as shown in Figure 18a,b [13,26,27], respectively. Here, 1 mm of the arrow length is equal to 1.8 µm/s. The direction and magnitude of the hydrodynamic flow v(τ, x, z) is influenced by both the direction of the heat flux q across the lower bounding surface and the character of the orientational defect. According to these calculations of the director fieldn(τ, x, z) across the nematic volume, the highest value of |∇n(x, z)| is reached in the vicinity of the orientational defect, and, as a result, the biggest thermally excited velocities occur in the vicinity of the lower hotter surface. In that case, the self-sustaining vortical flow is excited in the negative sense (anti-clockwise) around its center, as the value of the angle α reaches 60 • . Further increase in the angle value of more than 60 degrees leads to the formation of two vortices, one, larger, in the left-hand side of the nematic volume, another, smaller, in the right-side of the HAN channel, as shown in Figure 19a,b, Figure 20 and Figure 21a,b, respectively. With increasing the value of the angle α, the larger vortex begins to decrease in size, while the smaller vortex begins to increase. The size of the smaller vortex reaches its maximum size at the angle α = 120 • (see Figure 21b). Notice that the bigger self-sustaining vortical flow in the left-hand side of the HAN channel is thermally excited in the negative sense (anti-clockwise), whereas the smaller vortical flow in the right-hand side of the nematic volume is excited in the positive sense (clockwise) around their centers. With increasing the value of the angle α to more than 120 degrees, the larger vortex disappears and we have only one self-sustaining vortical flow in the positive sense (clockwise) around its center (see Figure 22a,b).
These results show that the direction of the heat flux across the hotter bounding surface with the orientational defect plays a crucial role in maintaining the thermally excited vortical flow in 2D HAN channel.
The value of dimensionless heat flux coefficient q 0 = T N I λ ⊥ d Q 0 is equal to 0.05 (∼4.2 × 10 −4 mW/µm 2 ). In general, under the above conditions, the picture of warming is such that only a small part of the nematic volume is involved in the heating process, while a large part of the volume of the fluid were not heated. It should be noted that in both cases the area of the greatest heating was shifted in the direction in which the heat flux was directed due to laser radiation.
We have reviewed the orientational dynamics in a microsized hybrid aligned nematic channel, where the nematic volume is confined by two horizontal and two lateral surfaces, and under the influence of the temperature gradient ∇T, when the nematic material is heated by a laser beam focused on the lower boundary with the orientational defect, whereas the rest bounding surfaces of the LC volume are kept at constant temperature. It has been shown that in the microsized HAN channel, due to interaction between ∇T and the gradient of the director field ∇n, in the LC volume with the orientational defect on the lower hotter boundary the vortical flow can be excited. The direction and magnitude of the hydrodynamic flow is influenced both by the heat flux q across the lower hotter boundary of the HAN channel, directed at the angle α with respect to the unit vectorî, and by the character of the orientational defect. Calculations show that the biggest thermally exited vortical flow occurs in the vicinity of the orientational defect. The analysis showed that at the same power of the laser radiation, with a change in the value of the angle α, the picture of the vortex flow in the HAN channel changes. When the heat flux q is directed at the angle α less than 40 degrees, in the nematic volume one the self-sustaining vortical flow in the negative sense (anti-clockwise) is excited around its center.
The increase in the angle value to more than 60 degrees leads to the formation of two vortices, one, larger, in the left-hand side of the LC volume, another, smaller, in the right-hand side of the LC volume, respectively. With further increasing the value of the angle α up to right angle and more, the larger vortex begins to decrease in size, while the smaller vortex reaches its maximum size at the angle 120 degrees. When the value of the angle α is more than 120 degrees, a single self-sustaining vortex flow is formed in a positive sense (clockwise) around its center [2].
Based on the analysis of the nature of the thermally excited vortical flow in the microsized HAN volume under the influence of the heat flux q directed across the boundary, we can conclude that the formation of the hydrodynamic flow requires such the radiation power that would pump energy into the liquid crystal phase for as long as possible. Another condition for the formation of the more powerful flow in the HAN channel is a planar anchoring condition on the boundary on which the laser radiation would be focused.
It should be noted that the vortical flow in the bulk of the microsized nematic volume is a unique phenomena only exhibited by liquid crystal systems, and expected to be applied for novel opto-thermal tweezers. Now, in the next paragraph, we turn to the description of the formation of vortex flows in microsized HAN channels with a free upper surface under the effect of focused laser irradiation in the bulk of the nematic phase.

Laser Excited Motion of Nematics Confined in Microsized Channel with a Free Surface
Despite the fact that certain qualitative and quantitative advances in a hydrodynamic description of the relaxation processes in the microsized nematic volume under the effect of the temperature gradient have been achieved, it is still too early to talk about the development of a theory which would make it possible to describe the dissipation processes in confined nematic channel with a free upper LC/air interface under the influence of the temperature gradient ∇T [11,12]. Thus, we are primarily concerned here on describing the set of numerical results which show the way how the temperature gradient caused by the laser radiation focused in the interior of the microsized hybrid aligned nematic (HAN) volume with a free upper LC/air interface can produce the hydrodynamic flow and, as a result, how it can deform the free LC/air interface [11,12].
This problem will be treated in the framework of the appropriate nonlinear extension of the Ericksen-Leslie theory [16,17], supplemented by the thermomechanical correction of shear stress [7,15], and the entropy balance equation [18]. The hydrodynamic model in which the above problem will be considered is the same as in Section 2.1, and is based on the interaction effect of the directorn and temperature T gradients with the velocity v = uî + wk field. The magnitude of the hydrodynamic flow v is proportional to v ∼ d η σ tm zx , where σ tm zx is the tangential component of the thermomechanical stress tensor, η is the viscosity of the nematic phase, and d is the thickness of the HAN channel. It will be shown below that the tangential component of the thermomechanical stress tensor plays a crucial role in the formation of the thermally excited vortical flow in the HAN channel with the free upper surface. In this case, the temperature gradient ∇T in the microsized HAN volume is formed by a laser beam focused in the interior of the HAN channel. At the beginning of the section, we give a detailed description of the boundary conditions for both the directorn and temperature T(x, z, t) fields, as well as for the velocity field v. The bounding condition for the velocity on the upper free LC/air interface Γ can be obtained from the linear momentum balance equation transmitted to the surface Γ. The temperature regime without the heat flux q across the free LC/air interface Γ will be analyzed. The second part of the section presents a numerical analysis of two modes of the thermally excited vortical flow in the HAN channel with the free upper surface. These modes are characterized by different laser radiation power and duration: slow and fast heating modes caused by the laser irradiation focused in the interior of the HAN channel with the free upper surface.
With this aim, we consider the HAN channel delimited by one lower horizontal solid surface, located at z = −d, one upper free flat LC/air interface, initially located at z = d, and two lateral solid surfaces at distance 2L on scale on the order of micrometers. The coordinate system defined by our task is the same as in the previous paragraph.
Therefore, the hybrid aligned nematic phase contains a gradient of ∇n from planar orientation on the lower and both lateral surfaces to homeotropic orientation on the upper free LC/air interface Γ, i.e., Here, , 1 is the normal to the free LC/air interface Γ at any time and is directed from the nematic phase into air, H(x, t) is the height of the LC film on the top of the smooth surface, and H ,x = ∂H ∂x . We consider the temperature regime without the heat flux q across the free LC/air interface Γ whereas on the rest boundaries the temperature is kept constant We will assume the no-slip boundary conditions for the nematogenic molecules on these solid bounding surfaces, i.e., v −L<x<L,z=−d = v x=±L,−d<z<d = 0, (26) where v = uî + wk = −∇ ×ĵψ is the velocity vector with the horizontal u ≡ v x (t, x, z) and vertical w ≡ v z (t, x, z) components. The bounding condition for the velocity on the upper free LC/air interface Γ can be obtained from the linear momentum balance equation transmitted to the surface Γ. In our case, that balance leads to the tangential and normal is an additional unit tangent vector, γ is the LC/air is the curvature of free LC/air interface Γ at any time, and σ is the full stress tensor (ST). Taking into account the microsized HAN volume, one can assume the mass density ρ to be constant across the sample, and thus one deal with an incompressible fluid. The incompressibility condition ∇ · v = 0 assumes that The hydrodynamic equations describing the reorientation of the nematic phase in 2D case, when the system is subjected to a temperature gradient ∇T, due to uniform heat flux q, can be derived from the torque balance equation Equation (6), the linear momentum equation for the velocity field v, which can be written in the form of Equation (14), and the heat conduction equation [11,12] where q = −T δR δ∇T denotes the dimensional heat flux in the HAN channel, ∆ 2 is the dimensional heat source coefficient, α is the coefficient of absorption, P 0 is the laser beam power, ∆ is the Gaussian spot size, and t in is the duration of the energy injection into the HAN channel. Now the dynamics of the height H(t, x) of the LC/air interface under the influence of the temperature gradient can be obtained by solving the system of nonlinear partial differential Equations (6), (14) and (30) with the appropriate boundary and initial conditions. Equations (27) and (28), together with the torque balance Equation (6), transmitted to the LC/air interface, can be combined to yield equation for the height H(t, x) in the form [11,12] where u Γ and w Γ are the horizontal and vertical components of the velocity v on the LC/air interface Γ, respectively. Below, the set of dimensionless analog balance Equations (6), (14) and (30) will be considered. The dimensionless entropy balance equation now can be written as [11,12] where the extra one parameter of the LC system is δ 5 = 2α The evolution of the free LC/air interface under the effect of the temperature gradient ∇χ, caused by the laser irradiation focused in the interior of the HAN channel, is governed by Equations (6), (14) and (32), together with the boundary conditions at the solid boundaries [11,12] (n x ) x=±10,−1≤z≤1 = 0, (n x ) −10≤x≤10,z=−1 = 1, and at the flexible free LC/air interface Γ [1] ( n · ∇χ) Γ = 0, ( n · ν) Γ = −1, respectively, and the initial condition in the form of Equation (22). Here, Ψ = (ψ ,xx , ψ ,xz , ψ ,zz ), and both the matrixB and vector C are given in the Appendix, whereas the vectorn el (x, z) is obtained from Equation (6), with ψ ,x = ψ ,z = χ ,x = χ ,z = 0. Now the dimensionless height H(τ, x) = H(τ, x)/d of the nematic-air interface at any time τ can be calculated as [11,12] where w(τ, x) x∈Γ = (ψ ,z ) x∈Γ and u(τ, x) x∈Γ = (ψ ,x ) x∈Γ are the vertical and horizontal components of the velocity vector v = uî + wk = −∇ ×ĵψ on the interface Γ. Notice that the overbar in the function H has been (and will be) eliminated in the last, as well as in the following equations. Thus, when the directorn is strongly homogeneously anchored to the lower boundary and planar to the lateral boundaries of the HAN channel, the value ofn has to satisfy the boundary conditions (33) and (34) and its initial orientation (22), and then, under the action of the viscous, elastic, and thermomechanical forces, is allowed to relax to its equilibrium valuen el (x, z).

A. Slow heating mode when a laser beam is focused in the interior of the HAN channel
Let us initially consider the case when the hybrid aligned nematic microvolume is heated by the laser beam focused in the interior of the HAN channel with the free upper boundary. In this case the value of the dimensional heat flux coefficient is equal to O 0 ∼ 0.5 W, whereas the heating occurs during time t in ∼ 2.0 µs. Figure 24 shows the distance x dependence of both the dimensionless height h(x, τ) = H(x, τ) − 1 (Figure 24a) of the LC/air interface and the dimensionless temperature χ(τ, x) (Figure 24b) on the free LC/air interface Γ [11,12], during the slow heating mode when the laser beam is focused in the interior (x = 0.0 and z = 0.93) of the HAN channel, at different times τ i = 2 i × 10 −5 (i = 1, . . . , 10), respectively, whereas Figure 25 shows the distance x dependence of both the horizontal u(τ, x) (Figure 25a) and vertical w(τ, x) (Figure 25b) components of the vector v = uî + wk = −∇ ×ĵψ on the LC/air interface. According to these calculations [10,11], the evolution of the height h(τ, x) of the LC/air interface is characterized by the wavelike profile along the x-axis (−0.5 < x < 0.5). At the final stage of the evolution process, for τ = τ in , the highest value of |h| ∼ 2 × 10 −2 is reached in the vicinity of points x ∼ ±0.25, whereas the evolution of the temperature χ is characterized by symmetric profile of χ(τ, x) x∈Γ with respect to the middle point (x = 0.0) of the LC/air interface Γ (see Figure 24b). In that case, during the heat step (τ ∼ τ in ∼ 0.01 (∼2 µs)) the evolution of the temperature profile χ(τ, x) x∈Γ is characterized by its strong growth in the vicinity of the middle point x = 0.0, up to the highest value of 0.987 (∼307 K), whereas the evolution of the dimensionless height h of the LC/air interface is characterized by two combs with the highest value of |h| ∼ 0.02 (∼0.01 µm), which are directed in the opposite sense with respect to their center x = 0.0 (see Figure 24a). The thermally excited flow in the case of the slow heating mode, with O 0 = 0.5 W and t in ∼ 2.0 µs, is characterized by maintaining of two vortices as shown in Figures 25a,b and 26 [11,12].  Figure 25. The same as in Figure 24, but the distance x dependence of both the horizontal u(τ, x) (a) and vertical w(τ, x) (b) components of the velocity vector v = uî + wk = −∇ ×ĵψ on the LC/air interface during the slow heating mode [11,12].  Figure 27 shows how three vortices can be maintained in the HAN channel with a free LC/vacuum interface Γ, during the fast heating mode, one biggest vortical flow in the vicinity of the heat source and directed in the negative sense (anticlockwise) around their center x = 0.0, z ∼ 0.93, and two smallest vortices, which are settled down close to the points x = ±0.13 and z ∼ 0.93 [11], respectively. These calculations also show that the range of distance z, counted from the lower boundary of the HAN channel, over which the laser beam cannot disturb the nematic phase, is 0.8 ≤ z ≤ 1.0, i.e., which is approximately 80% of the LC sample (see Figure 27). Notice that the duration of the energy injection τ in into the LC sample is restricted only by the nematic phase stability. Further calculations (cooling mode), based on the nonlinear extension of the Ericksen-Leslie theory, show that the nematic material settles down to the rest during the time term τ 8 ∼ 2.56 (∼0.5 s), after switching off the laser power, where both the horizontal u and vertical w components of the velocity v are equal to zero, and the temperature field χ across the HAN channel finally downfalls to the value on the lower and two lateral boundaries.
The results presented in [11][12][13] indicate that with a change in the nature of laser radiation focusing, from the surface of the nematic channel to its volume, leads to a significant change in hydrodynamic flows. The question of how the depth of laser radiation focus affects the nature of hydrodynamic flows in the hybrid aligned nematic channel will be discussed in the next paragraph.  (Fast heating mode). Distribution of the velocity field v = uî + wk = −∇ ×ĵψ in the HAN channel after heating during τ = τ in [12]. Here, 1 mm of the arrow length is equal to 0.4 µm/s.

How the Depth of Laser Radiation Focus Affects the Nature of Hydrodynamic Flows in Nematic Channel
In this paragraph, we will consider the effect of the depth of focus of laser radiation in the microscopic volume of the HAN channel with the free LC/air interface on the nature of formation of the hydrodynamic flow and temperature distribution [11,12]. Figures 28 and 29 show the distribution of the dimensionless temperature χ(τ, x = 0.0, z) along the z-axis (−1.0 ≤ z ≤ 1.0), when the laser beam is focused in the center (x = 0.0) of the HAN channel, at different depths [11,12]: z 0 = 0.80 (see Figure 28a), z 0 = 0.90 (see Figure 28b), z 0 = 0.94 (see Figure 29a), and z 0 = 0.98 (see Figure 29b), respectively.
The heating mode when the laser beam is focused in the interior of the HAN channel is given at different times τ i = 2 i × 10 −5 (i = 6, . . . , 10), respectively. It has been shown that as the focus of the laser beam is shifted in the depth of the nematic microvolume, the temperature profiles across the nematic channel do not undergo the crucial change. For instance, in the case when the laser beam is focused on the maximum depth (z 0 = 0.8), the heating does not reach the LC/air interface (see Figure 28a [11,12]).  In turn, the velocity profiles across the HAN channel undergo the crucial change. Figures 30-33 show the distribution of the horizontal u(τ, x = 0.0, z) and vertical w(τ, x = 0.0, z) components of the velocity vector v = uî + wk = −∇ ×ĵψ along the z-axis (−1.0 ≤ z ≤ 1.0), when the laser beam is focused in the center (x = 0.0) of the HAN channel, at different depths. Figure 30a,b [11,12] show the distribution of the horizontal u(τ, x = 0.0, z) component of the velocity v along the z-axis (−1.0 ≤ z ≤ 1.0), when the laser beam is focused in the center (x = 0.0) of the HAN channel, but at different depths [11,12]: Figure 30a, at z 0 = 0.80 and Figure 30b, at z 0 = 0.90, respectively. In turn, the Figure 31a,b show the distribution of the horizontal u(τ, x = 0.0, z) component of the velocity v along the z-axis (−1.0 ≤ z ≤ 1.0), at different depths [11,12], Figure 31a, at z 0 = 0.94, Figure 31b, at z 0 = 0.98, respectively. Figure 32a   It has been shown that as the focus of the laser beam is shifted in the depth of the HAN channel in the vicinity of the LC/air interface, the horizontal component of the velocity u(τ, x = 0, z) changes its direction from negative to positive, approximately at the point x 0 = 0.0, z 0 ∼ 0.9, whereas the vertical component of the velocity w(τ, x = 0, z) rapidly drops to zero. It should be noted that the greatest value of u(τ, x = 0, z), directed in the positive sense in the vicinity of the LC/air interface, is achieved in the case when the laser beam is focused on the maximum depth of penetration (z 0 = 0.8) in the LC volume, whereas the greatest value of u(τ, x = 0, z), directed in the negative sense in the vicinity of the LC/air interface, is achieved in the case when the laser beam is focused on the minimum depth of penetration (z 0 = 0.98) in the LC volume. In both of these cases, the vertical component of the velocity vector w(x = 0, z, τ) at the free LC/air interface is almost zero. Therefore, this distribution of components of the velocity field shows that in the area close to the LC/air interface (0.8 < z < 1.0), the vortex flow is excited due to the energy pumping by laser radiation, similar to what is shown in Figure 26. These calculations [11,12] also show that with further penetration of the injecting energy to the bulk of the LC phase, from x 0 = 0.0, z 0 = 0.98 to x 0 = 0.0, z 0 = 0.8, the thermally excited vortical flow changes the direction from anticlockwise, around the point x = 0.0, z = 0.98, to clockwise, around the point x = 0.0, z = 0.8, approximately at the point x 0 = 0.0, z 0 ∼ 0.9.
Based on the numerical results describing the formation of hydrodynamic flows in microsized HAN channels with free surfaces under the effect of laser radiation focused in the volume of the LC phase, the following conclusion can be made. First, the calculations, based on the appropriate nonlinear extension of the classical Ericksen-Leslie theory, show that due to the interaction between ∇T and the gradient of the director field ∇n, in the nematic volume the thermally excited three-vortical fluid flow is maintained. Second, the direction and magnitude of the hydrodynamic flow at a fixed time of pumping energy and the laser output power are affected by the depth of the laser injection. The above calculations also show that the range of distances measured from the lower solid boundary at which the laser beam cannot disturb the nematic phase is approximately 80% of the nematic sample.
It should be noted that the vortical flow in homeotropically oriented LC film doped by chiral molecules, using the circular polarization techniques, recently has been observed [10]. It has been shown that at the beginning of laser irradiation, the thermocapillary radial flow is transformed into a circular flow around the position of the laser spot on the free LC interface. The formation of the circular flow on the top of the LC film has been ascribed to thermocapillary convection in the LC sample.
In turn, in the above-mentioned cases [11,12], the vortical flow occurred in the vicinity of the free LC/air interface and penetrated to the bulk of the LC sample. The mechanism responsible for the occurrence of the vortical flow near the LC/air interface is based on the coupling between the director and the temperature gradients initiated by the laser beam radiation. Thus, this vortical flow is a unique phenomenon that is exhibited only by liquid crystal systems, and is expected to be applied to new optical-thermal tweezers [1].

Heat Driven Nematic Flow in Cylindrical Microfluidic Channel
The objective of this section is to analyze the response of the nematic phase confined in the cylindrical micro cavity between two horizontal coaxial cylinders under the influence of the temperature gradient ∇T directed from the inner (outer) cooler (warmer) to outer (inner) warmer (cooler) cylinders [28]. Therefore, we are primarily interested here in describing how the temperature gradient across the microvolume cavity between two coaxial cylinders can produce the hydrodynamic flow v. It has been treated in the framework of the classical Ericksen-Leslie theory [16,17], supplemented by the thermomechanical correction of shear stress [7,15], and the entropy balance equation [18]. The hydrodynamic model in which the above problem will be considered is the same as in Section 2.1, and is based on the interaction effect of the directorn and temperature T gradients with the velocity v field. The magnitude of the hydrodynamic flow is proportional to v ∼ d η σ tm zr , where σ tm zr is the tangential component of the thermomechanical stress tensor, η is the viscosity of the nematic phase, and d is the thickness of the hybrid aligned nematic (HAN) cavity confined between two infinitely long horizontal coaxial cylinders. It will be shown below that the tangential component of the thermomechanical stress tensor plays a crucial role in the formation of the thermally excited vortical flow in the HAN cavity.
To fix ideas and notations, in ref. [28] was considered the HAN system composed of asymmetric polar molecules, such as cyanobiphenyls, at the density ρ, and confined between two infinitely long horizontal coaxial cylinders with radii R 1 and R 2 . Here, R 1 < R 2 , that imposes a preferred orientation of the average molecular directionn on both bounding surfaces, for instance, homeotropic on the inner cooler (T in = T 1 ), and planar, on the outer warmer (T out = T 2 ) bounding cylinders. Therefore, it was considered the HAN system under the influence of the radially directed temperature gradient ∇T parallel to the unit cylindrical vectorê r along the radius r [28]. The coordinate system defined by this geometry assumes that the directorn lies in the rz plane, whereê z is the unit vector which coincides with the planar director orientation, for instance, on the outer warmer (inner cooler) bounding cylinder (ê z n r=R 2 (=R 1 ) ),ê r denotes the unit vector along the radius r, andê ϕ =ê z ×ê r is the tangential unit vector.
Assuming that the temperature gradient ∇T varies only in the r direction, ∇T = ∂T(t,r) ∂rê r , and the directorn belongs to the rz plane. Further, it was proposed that in the case of the infinitely long cylinders, the components of the directorn = sin θ(t, r)ê r + cos θ(t, r)ê z , as well as the rest of the relevant physical quantities depend only on the radius r and on the time t. Here θ denotes the angle between the direction of the directorn and the unit vectorê z directed parallel to the long cylinder s axis. Moreover, it was assumed the homeotropic and homogeneous strong anchoring conditions on the outer and inner cylinders [28] and no-slip boundary conditions for the nematogenic molecules on both bounding cylindrical surfaces, i.e., v(r) r=R 1 = v(r) r=R 2 = 0, respectively. Upon assuming an incompressible fluid ∇ · v = 0, the hydrodynamic equations describing the reorientation of the HAN system confined in the microvolume between two horizontal coaxial cylinders under the influence of the temperature gradient ∇T directed from the inner cooler to outer warmer cylinders can be derived from the balance of elastic, viscous, and thermomechanical torques T el + T vis + T tm = 0, the Navier-Stokes equation for the velocity field v, excited by ∇T, and the equation for heat conduction. Here, T el = δF el δn ×n, T vis = δR vis δn ×n, T tm = δR tm δn ×n,ṅ = dn dt is the material derivative of the directorn, F el = 1 2 K 1 (∇ ·n) 2 + K 3 (n × ∇ ×n) 2 is the elastic energy, and R = R vis + R tm + R th is the full dimensionless Rayleigh dissipation function composed by the viscous, thermomechanical, and thermal contributions, respectively. The incompressibility condition together with the no-slip condition (37) implies the existence only one nonzero component for the vector v, viz. v(t, r) = v z (t, r)ê z ≡ u(t, r)ê z . In the cylindrical coordinate system the dimensionless torque balance equation takes the form [28] γ 1 (χ)θ ,τ = A(θ)u ,r + (G(θ)θ ,r ) ,r − 1 2 G ,θ (θ)θ 2 ,r + θ ,r r G(θ) where γ 1 = γ 1 (χ)/γ 10 , δ 6 = ξ T N I K 10 is the parameter of the nematic system, θ ,τ = ∂θ(τ,r) ∂τ , θ ,r = ∂θ(τ,r) ∂r ,u ,r = ∂u(τ,r) ∂r , and A(θ) = − 1 2γ 10 [γ 1 (χ) + γ 2 (χ)u ,r (τ, r) cos 2θ] and G(θ) = K 1 (χ) K 10 sin 2 θ + K 3 (χ) K 10 cos 2 θ are the hydrodynamic and elastic functions, and χ(τ, r) = T(τ, r)/T N I is the dimensionless temperature. Here, γ 10 and K 10 are the highest values of the RVC γ 1 (χ) and of the splay constant K 1 (χ) in the temperature interval ∆χ = χ 2 − χ 1 belonging to the nematic phase, χ 1 = T 1 /T N I , χ 2 = T 2 /T N I , τ = K 10 γ 10 d 2 t is the dimensionless time, r = r d is the dimensionless radius, and d = R 2 − R 1 is the capillary gate, respectively.
In the reviewed case [28], the HAN capillary is heated from above with the dimensionless temperature difference ∆χ = 0.0162 (∼5 K). The solution of the system of nonlinear partial differential equations Equations (39)-(41), together with the boundary conditions (42)-(46), has been obtained by means of the numerical relaxation method [24]. The relaxation criterion = |(θ (m+1) (τ, r) − θ (m) (τ, r))/θ (m) (τ, r)| was chosen to be 10 −4 , and the numerical procedure was then carried out until a prescribed accuracy was achieved. Here, m is the iteration number, and τ R is the relaxation time of the HAN system.
All calculations were carried out for a = 1.0.
In that case, the lower cooler surface is kept at constant temperature χ 1 = 0.97 (T 1 ∼ 298 K), and the relaxation of the temperature field χ(τ, r) across the HAN cavity is characterized practically by the linear increasing of the values of χ(τ, r), from χ 1 to χ 2 (see Figure 35).
According to these calculations [28], the evolution of the dimensionless velocity u(τ, r) in the HAN cavity between two infinitely long coaxial cylinders is characterized by the monotonic decrease of |u(τ, r)| upon increasing τ, before getting to the equilibrium distributions u eq (r) = u(τ R , r) = u(τ 10 , r) across the HAN cavity. That distribution is characterized by the minimum near the middle part of the cavity, where the hydrodynamic flow is directed in the negative sense (see Figure 34b, curve (10)).
It should be noted that the equilibrium hydrodynamic flow u eq (r) change direction, from the negative to positive sense, across the full thickness of the HAN cavity, after changing the temperature gradient ∇χ direction, from the cooler (T in = 298 K) inward to the warmer outward (T out = 303 K) cylinders on the the warmer (T in = 303 K) inward to the cooler outward (T out = 298 K) (see Figure 36 and Figure 37, curves (2) and (1), respectively.).  Figure 36. The distance r dependence of the equilibrium velocity u eq (r) across the nematic cavity between two a ≤ r ≤ a + 1 infinitely long coaxial cylinders, with the anchoring hybrid condition in the form of Equation (42), under influence of the ∇χ, directed from the cooler (warmer) inward χ 1 (χ 2 ) to warmer (cooler) outward χ 2 (χ 1 ) cylinders (see curves (2) and (1) (2) Figure 37. The same as in Figure 36, but the value of a is equal to 0.1 [28].
These calculations [28] show that increasing of the nematic cavity size d ≡ R 2 − R 1 between two infinitely long coaxial cylinders, when the temperature gradient ∇χ is directed from the inward to outward cylinders, leads to increase of the maximum of the absolute magnitude of the dimensionless velocity u eq (r), from u max eq ∼ 0.4 (∼11 µm/s), at a = 1 (R 2 = 2R 1 ) (see Figure 36a) to u max eq ∼ 1.5 (∼41.3 µm/s), at a = 0.1 (R 2 = 11R 1 ) (see Figure 37), respectively.
According to these calculations [28], in the case when the hybrid oriented anchoring conditions for director is described both by Equations (42) and (43), with (χ 1 ) r=a > (χ 2 ) r=a+1 , one has arrived to the picture where the nematic fluid settles down to the stationary flow regime in the negative verse (see Figures 36-39). (2) Figure 39. The same as in Figure 38, but for a value of a equal to 0.1 [28].
In the first case (Equation (42)), the maximum of the velocity field u max eq (r) is found in the vicinity of the cooler outward cylinder (see Figure 37), whereas in the second case (Equation (43)), the maximum of u max eq (r) is found in the vicinity of the warmer inward cylinder (see Figure 39), respectively. Physically, this means that the thermomechanical force can overcome the viscous and elastic forces, and the nematic drop, confined between two infinitely long hybrid-oriented cylinders, begins to move in the horizontal direction with the stationary distributed velocity across the nematic cavity. Thus, we arrive at the picture where there is a balance between the applied temperature gradient and the viscous, elastic, and anchoring forces, and, in general, the nematic fluid settles in a stationary flow regime along the long axis of the cylinder s. The magnitude of the hydrodynamic flow u eq is proportional to the tangential component of the thermomechanical stress tensor σ tm zr and the cavity size R 2 − R 1 , and the direction of u eq is influenced by both the direction of the heat flux and the character of the preferred anchoring of the director to the restricted cylinders.
The highest value of the dimensionless velocity u max eq (r) is built up in the hybrid aligned nematic cavity, when the temperature gradient ∇χ directed from the outward cooler (T out = 298 K) to inward warmer (T in = 303 K) cylinders, in the vicinity of the homogeneously aligned inward cylinder (see Equation (43)), and equal to ∼2 (∼55 µm/s) (see Figure 39), and is directed in the negative sense.
It has been found in ref. [28] that the size of the HAN cavity a = R 1 has a pronounced effect on the magnitude of u max eq (a) (see Figure 40).  Figure 40. Dependence of u max eq (a) on the dimensionless size of the nematic cavity a = R 1 R 2 −R 1 , with the anchoring hybrid condition in the form of Equation (43), for two cases: first, (a) when the heating mode is directed from inward to outward bounding cylinders (χ r=a+1 > χ r=a ), second, (b) when the heating mode is directed from outward to inward bounding cylinders (χ r=a > χ r=a+1 ) [28], respectively. Figure 40a shows the effect of a on the magnitude of the highest dimensionless velocity u max eq (a), when the temperature on the inward cylinder χ r=a is greater than the temperature on the outward cylinder χ r=a+1 . In the case of a = 0.03, i.e., when the radius of the outward cylinder R 2 ∼ 34R 1 , the maximum value of the dimensionless equilibrium velocity has the biggest value which is equal to ∼3.82 (∼105 µm/s), and is directed in the negative sense. Note that the further decrease of a leads to decreasing of the highest dimensionless velocity u max eq (a). Figure 40b shows the effect of a on the magnitude of the highest dimensionless velocity u max eq (a), in the case when χ r=a+1 > χ r=a . Here, the biggest value is equal to ∼26 (∼0.7 mm/s), at a = 0.01 (R 2 ∼ 112R 1 ), and ∆χ = 0.0162 (∼5 K), and directed also in the negative sense.
It should be noted here that the highest temperature difference ∆ = T in(out) − T out(in) across the nematic cavity is still such that both T in and T out fall within the stability region of the nematic phase. The temperature difference, for instance, in ∼5 K in the nematic cavity confined between two infinitely long horizontal coaxial cylinders can be built up by using the laser-induced heating [10][11][12]14]. Indeed, since the laser beam can be focused to its diffraction limit, it can be used to inject energy at scales that are difficult to reach with other techniques [4].
In summary, we have reviewed the evolution of the directorn(t, r), velocity v(t, r), and temperature T(t, r) in the hybrid aligned nematic cavity between two infinitely long coaxial cylinders to their equilibrium values, under influence of the radial temperature gradient directed from the cooler to warmer bounding surfaces. These calculations [28], based upon the classical Ericksen-Leslie theory, show that the microvolume of the HAN material under influence of the temperature gradient, settles down to a stationary flow regime in the horizontal direction. It has also been shown that the magnitude of that velocity u eq is proportional both to the tangential component of the stress tensor σ tm zr and the size d of the nematic cavity gap between two infinitely long coaxial cylinders, and the direction of v = uk influences both the direction of heat flux and the character of the preferred anchoring of the director to the bounding surfaces. These calculations also show that the optimal pumping effect in the microvolume HAN cavity between two coaxial cylinders build up under influence of the radial temperature gradient, when a = R 1 R 2 −R 1 ∼ 0.01, or when the radius of the homogeneously anchored warmer outward cylinder is greater than the radius of the homeotropically anchored cooler inward cylinder, approximately, in 100 times. Note that the above-mentioned flows are described by a continuous approach, where the manipulating fluid s length scales less than ten micrometers is reflected in the magnitudes of the parameters δ i (i = 6,7,8,9). It should be pointed out that the pumping effect under the influence of the radially directed thermal gradient is not found in the nematic cavity with the homeotropically or homogeneously aligned molecules at the restricted cylinders.
We believe that the present review can shed some light on the problem of the reorientation process in the microfibers under the influence of the temperature gradient [3,4,29]. We also believe that the precise handling of the liquid crystal microvolumes, which requires self-contained micropumps of small package size exhibiting a continuous volume flow, can be developed utilizing the interactions of both the director and velocity fields with the radially directed temperature gradient. Hence, the possible pumping technique described above appears applicable to various pumping strategies not involving mobile parts.

Conclusions
This review discusses some recent numerical advances in predicting the structural and hydrodynamic behavior of thermally excited vortical flows in microsized nematic channels caused by a laser beam focused both inside the nematic volume and on bounding surfaces. Despite the fact that certain qualitative and quantitative advances have been made in the hydrodynamic description of relaxation processes in microsized nematic volumes under the influence of a temperature gradient, there are still a number of questions concerning dissipative processes in confined liquid crystal phase with complex boundaries of nematic channels or with a free upper boundary of the liquid crystal material/air interface that need to be clarified.
It should be noted that in the case of a thick horizontal layer, with a thickness d of several tens of millimeters, the Rayleigh-Benard instability may occur in a stationary liquid crystal volume heated from below. In that case, the LC volume is driven by maintaining the ∇T, directed from the lower hotter boundary T lw to the upper cooler T up one. This instability occurs when the driving ∆T = T up − T lw is strong enough to overcome the dissipative effects of thermal conduction, thermomechanical effect, and viscosity. The control parameter describing the instability, the Rayleigh number R = αgd 3 ∆T/(α 4 /2ρ)κ ⊥ , where g is the gravitational acceleration, α is the isobaric thermal expansion coefficient, α 4 /2ρ corresponds to the isotropic kinematic viscosity, and κ ⊥ is the perpendicular thermal diffusivity. The instability occurs at value R = R c ∼ 1708, independent of the fluid under consideration [30]. Taking into account that the size of the HAN channel d ∼ 5 − 10 µm, in our case R R c , and the driving force is weak enough to set up of convection via the Rayleigh-Benard mechanism [7]. Therefore, we are primarily concerned here with the description of the physical mechanism responsible for the horizontal flow in the HAN channel confined between two solid surfaces, which is excited by the temperature gradient ∇T, and the magnitude of that flow v is σ tm , whereas the direction of v is influenced both by the direction of the heat flux q and the character of the preferred anchoring of the average molecular directionn to the restricted surfaces. A number of suitable directions for the application of microliter LC materials are associated with thermally controlled manipulations in biology, medicine, and many areas of engineering [3][4][5][6]. This dynamic new field has also drawn attention due to advances in drug delivery [31,32]. There is another promising dynamic field of application of LC materials of microliter volumes. It is the artificial LC lenses with tunable focal length, which are called tunable lenses [33][34][35]. The idea of using a sessile drop of liquid as an optical lens dates back to the late 17th century. Over the last decade, researchers have become interested not only in liquid, but also in LC lenses due to exceptional properties of the LC/air interface [33,34]. In LC lenses, the focal length can be altered by varying the curvature of free LC/air interface in response to external stimuli or actuations. Taking into account the smoothness and flexibility of the LC/air interface, it is possible to solve a number of problems with optimizing these optical systems. For example, using laser-induced heating of microsized liquid crystal droplets, it is possible to manipulate the curvature of the LC/air interface and, as a result, the tunable focal length of lenses with curved surfaces. Another possible application of microfluidics is liquid crystal photonics [36].
Therefore, further study of a wider range of problems related to understanding how elastic soft matter, such as liquid crystals, begins to move under the influence of a temperature gradient, will require additional efforts, which will ultimately lead to an increase in our knowledge in the field of materials science. We also hope that this review will draw attention to the problem of manipulating the microsized volumes of liquid crystals by forming temperature gradients.