Bending of Thin Liquid Crystal Elastomer under Irradiation of Visible Light: Finsler Geometry Modeling

Abstract

In this paper, we show that the 3D Finsler geometry (FG) modeling technique successfully explains a reported experimental result: a thin liquid crystal elastomer (LCE) disk floating on the water surface deforms under light irradiation. In the reported experiment, the upper surface is illuminated by a light spot, and the nematic ordering of directors is influenced, but the nematic ordering remains unchanged on the lower surface contacting the water. This inhomogeneity of the director orientation on/inside the LCE is considered as the origin of the shape change that drives the disk on the water in the direction opposite the movement of the light spot. However, the mechanism of the shape change is still insufficiently understood because to date, the positional variable for the polymer has not been directly included in the interaction energy of the models for this system. We find that this shape change of the disk can be reproduced using the FG model. In this FG model, the interaction between $\sigma$, which represents the director field corresponding to the directional degrees of LC, and the polymer position is introduced via the Finsler metric. This interaction, which is a direct consequence of the geometry deformation, provides a good description of the shape deformation of the LCE disk under light irradiation.

1. Introduction

Liquid crystal elastomer (LCE) is a material in which LC molecules are chemically bonded to polymers, and its most remarkable property is that the anisotropy in the direction of nematic ordering of the LC molecules is reflected in the macroscopic shape [1,2,3,4,5]. Photoinduced bending of liquid-crystalline gel, which includes azobenzen monomers, shares the same mechanism [6,7]. In other words, the macroscopic shape of the LCE is sensitive to the microscopic orientation of the LCs and vice versa. Indeed, LCE elongates into a spontaneously chosen direction at low temperature if no constraint is imposed on the LC molecules, while at high temperatures, the elongation is suppressed due to the phase change from nematic to isotropic. External fields such as electric field and mechanical stresses also elongate the LCEs, and these field-driven elongations have been studied intensively [8,9,10,11,12].

Because of these interesting properties, a great deal of attention has been paid to LCEs, and studies of their technological applications, such as artificial muscles, have been conducted [13,14,15]. A polymer-dispersed LCE has recently been proposed as a functional and shape-programmable material for additive manufacturing [16]. A cellulose liquid crystal motor has also been proposed [17]. In the statistical mechanical studies of LCEs, the variable $\sigma$ for directors and strains are introduced, and the obtained results successfully describe the experimental data [2,3,4,5]. However, the strains themselves are not always identical to the LCE shape, so the position vector is better or more straightforward for the analyses of material shape such as polymers [18,19].

In Refs. [20,21], the interaction between the LC molecules and polymers is coarse-grained and implemented in the 3D Finsler geometry (FG) model, and the elongation phenomenon is successfully explained. The soft elasticity, in which the stress-strain diagram includes a plateau, is also reproduced by the FG model. This soft elasticity arises because the director orientation and the shape of material interact with each other. The so-called J-shaped stress-strain diagrams of soft biological tissues such as animal skin and muscles are also reproduced by the 2D FG model [22]. In the FG models for the soft-elasticity and J-shaped diagram, we discard the detailed information of the interaction between the LCs and polymers, and instead, we modify the underlying geometry such that the internal metric function of the materials is changed from Euclidean to Finsler. In this sense, the interaction is determined in a more basic and constructive manner from the geometries inside the material and the space ${\Re }^3$. As a result of this procedure, the interaction between the direction $\sigma$ of the directors and the position $\mathbf{r}(\in {\Re }^3)$ of the polymers is automatically introduced in the model. This is in sharp contrast to the ordinary modeling techniques, in which an interaction energy is explicitly introduced in the Hamiltonian. Therefore, this FG modeling technique is completely new and should hence be checked thoroughly by application it to a variety of experimentally observed phenomena.

In this paper, we apply the 3D FG model to the interesting experimental observation that a thin and small LCE disk moves on the water due to light irradiation [1]. This phenomenon has attracted a great deal of attention in engineering applications because it is considered a model for a micro robot actuated by a non-mechanical stimulus [6,7,17]. The driving force that moves the disk on the water comes from the shape change. By moving a spot irradiated by light from the center to the radial direction on the disk surface, the position of the bending on the edge of the disk can be controlled. This shape change at the edge moves in the same direction as the spot; at the same time, the water is pushed in the same direction that the bending edge moves, and then, by the action-reaction principle, the disk moves back or swims in the opposite direction on the water. We expect that such a shape change can be studied by using the FG model because the position variable $\mathbf{r}$ is directly included as a variable. Therefore, it is worthwhile to study whether the experimentally reported shape change of the disk is consistent with the results of the 3D FG model.

We should comment that the FG modeling technique can also be applied to a deformation of LCE piece in more general situation even for composite materials such as polymer-dispersed LCE [16]. Such an applicability of FG model is naturally expected because the deformation of the variable $\mathbf{r}$ is automatically determined from the variable $\sigma$, and the deformation of $\mathbf{r}$ is independent of the reason for distorting $\sigma$ because of the implemented interaction between $\mathbf{r}$ and $\sigma$.

2. Reported Experiment

We describe the experiment reported in Ref. [1] and summarize the current understanding of the shape deformation by light irradiation. A dye-doped nematic LCE is used in the experiment; the disk size is 5 mm in diameter and 0.32 mm thick, so the ratio $D/H$ of thickness H and diameter D is approximately $D/H=16$. This small piece of LCE floats nearly motionless on the surface of a water reservoir, which has a depth of 2 cm. The direction of the nematic director, which is the constituent LC molecule, is parallel to the surface.

This LCE sample is illuminated from above by an argon-ion laser with a beam width of 3 mm. The light spot is located on the center of the disk in the beginning and is then moved toward an edge along the direction perpendicular to the LC alignment direction. Then, the sample moves opposite to the direction in which the light spot moves (see Figure 1).

The authors of Ref. [1] discussed why the LCE sample moves on the water surface. One possible reason that they suggested was the shape deformation, as shown in Figure 1. A part of sample edge is lifted and strongly deformed compared to the other part of the edge, and this deformed position moves with the movement of the light spot. This deformed edge pushes water in the same direction as the light movement, and by the action-reaction principle, the LCE sample moves in the opposite direction. This is a rough outline of the mechanism of the swimming mentioned in the Introduction.

The effect of light on the dye-doped nematic LCE samples is also discussed in Ref. [1]. Based on these discussions, there are two possible effects: the first is dye-mediated heating due to the light absorption, with the increase in temperature causing a transition from the nematic ordered phase to the disordered or isotropic phase. The second also reduces the orientation order due to the effect of the photoisomerization of the dissolved dye. The dye-doped LC molecules play the role of an impurity if they change from the trans state to the cis state via irradiation. These two effects are expected to reduce the orientation ordering of directors, and this change in the director orientation leads to a shape change. We should note that these two effects have been throughly investigated and are now well understood [23,24,25,26,27,28].

However, it is unclear why a change in the orientation order causes shape deformation. This is actually very difficult to study from the microscopic perspective because the interaction between the LCs and polymers is unclear. Here, “unclear” means that the corresponding Hamiltonian in the statistical mechanical models remains unknown because the polymer position $\mathbf{r}$ is not directly used. One of the reasons for the lack of a Hamiltonian that includes $\mathbf{r}$ is that the interaction is very difficult to assume due to the hierarchical structure of polymer, which is composed of monomers and forms a crystalline/non-crystalline structure, with these structures also forming domains. All of these hierarchical structures are connected to the shape deformation; for this reason, a simple interaction energy is not expected as the Hamiltonian [29,30,31].

3. 3D FG Model of LCE

To understand why a dye-doped nematic LCE deforms under light illumination, we employ a model that is completely different from the models utilized in previous statistical mechanical studies of materials, as mentioned in the Introduction. The model in this paper is an FG model and is identical to the model used in Ref. [21]; however, we describe the model for the readers’ convenience.

First, the 3D lattices for the simulations are shown in Figure 2a–c. For the LCE sample in the experiment, the ratio $D/H$ of height H and diameter D is $D/H=16$, as mentioned in the previous section. Thus, we use three different lattices of $D/H=12$, $D/H=16$, and $D/H=20$, as shown in Figure 2a–c.

We show the data for the lattices, including the lattice of $N=$ 10,346, $D/H=8$, in

Table 1. The four different lattices used in this paper. The first lattice of size $N=$ 10,346 is used only for the calculation of Gaussian curvature.

N	${\mathit{N}}_{\mathit{B}}$	${\mathit{N}}_{\mathit{T}}$	${\mathit{N}}_{\mathbf{tet}}$	$\mathit{D}/\mathit{H}$	${\mathit{N}}^{\mathit{U},\mathit{L}}$	${\mathit{N}}_{\mathit{B}}^{\mathit{U},\mathit{L}}$	${\mathit{N}}_{\mathit{T}}^{\mathit{U},\mathit{L}}$
10,346	69,964	116,041	56,422	8	1344	3916	2573
10,566	70,503	116,037	56,099	12	1718	5023	3306
10,980	72,049	117,623	56,553	16	2088	6120	4033
10,563	67,752	109,460	52,270	20	2322	6814	4493

. The size of the lattice is given by $(N,N_B,N_T,N_{\mathrm{tet}})$, where N, $N_B$, $N_T$, and $N_{\mathrm{tet}}$ are the total numbers of vertices, edges (or bonds), triangles, and tetrahedrons, respectively. The Euler numbers $N-N_B+N_T-N_{\mathrm{tet}}(=1)$ of these lattices are identical to those of a tetrahedron, which is characterized by $(N,N_B,N_T,N_{\mathrm{tet}})=(4,6,4,1)$, because the cylinders are topologically identical to a tetrahedron. The symbols $N^U,N^L$ denote the total number of vertices on the upper and lower surfaces, respectively, and $N_B^U,N_B^L$ and $N_T^U,N_T^L$ respectively denote the total number of bonds and the total number of triangles on these two-dimensional surfaces. These data also satisfy $N^{U,L}-N_B^{U,L}+N_T^{U,L}=1$ because a disk is topologically identical to a sphere with a hole, e.g., a triangle.

To describe the model, we start with the continuous Gaussian bond potential such that

$$\begin{array}{cc}\hfill S_1(\mathbf{r},\sigma )& =\int \sqrt{g}{d}^{2}x{g}^{ab}\frac{\partial \mathbf{r}}{\partial x_a}·\frac{\partial \mathbf{r}}{\partial x_b},\end{array}$$

(1)

where $\mathbf{r}(\in {\Re }^3)$ is a material position, $x_a(a=1,2,3)$ is a local coordinate inside the material, g is the determinant of the Finsler metric $g_{ab}$, which is a $3\times 3$ matrix, and $g^{ab}$ is the inverse of $g_{ab}$. This continuous $S_1$ is discretized on 3D lattices such as those in Figure 2a–c, which are composed of tetrahedra. Please note that this $S_1$ depends not only on $\mathbf{r}$ but also on $\sigma (\in S^2:\mathrm{unit}\mathrm{sphere})$, which is a director field corresponding to the directional degrees of freedom of the LC molecule. The dependence of $S_1$ on $\sigma$ comes from the fact that $g_{ab}$ depends on $\sigma$. Note also that this $S_1$ is a 3D extension of 2D $S_1$, which is a model for the membranes [32,33,34,35,36,37,38], and 2D $S_1$ is an extension of the 1D model for polymers [39]. Therefore, this type of Gaussian potential $S_1$ in Equation (1) is of fundamental importance in the studies of polymers.

To evaluate the dependence of $S_1$ on $\sigma$, we write the discrete 3D Finsler metric $g_{ab}=\left(\begin{array}{ccc}\hfill 1/v_{12}^2\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1/v_{13}^2\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1/v_{14}^2\hfill \end{array}\right),$ which is obtained by replacing the element 1 of Euclidean metric ${\delta }_{ab}$ with $1/v_{ij}^2$. This Finsler metric is defined on the tetrahedron in Figure 3a, where the local coordinate origin is at the vertex 1. In this $g_{ab}$, $v_{ij}$ is given by

$$\begin{array}{c}\hfill v_{ij}=|{\mathbf{t}}_{ij}·{\sigma }_i|+v_0,{\mathbf{t}}_{ij}={\overrightarrow{\ell }}_{ij}/{\ell }_{ij},{\overrightarrow{\ell }}_{ij}={\mathbf{r}}_j-{\mathbf{r}}_i,\end{array}$$

(2)

where ${\mathbf{t}}_{ij}$ is the unit tangential vector from the vertex i to j (Figure 3b). This $v_{ij}$ corresponds to the unit of the Finsler length [20,40,41,42], and it depends on $\sigma$; hence, $g_{ab}$ depends on $\sigma$. We should note that $v_0$ in $v_{ij}$ plays the role of a cut-off because $1/v_{ij}^2$ is divergent if ${\mathbf{t}}_{ij}$ is vertical to ${\sigma }_i$. We fix $v_0$ to $v_0=0.001$ in this paper.

A discrete version of $S_1$ is obtained by replacing the differentials ${\partial }_a\mathbf{r}$ with differences such as ${\mathbf{r}}_j-{\mathbf{r}}_i$ and the integral $\int d^{2}x$ with the sum ${\sum }_i$. Including several terms, we have the Hamiltonian $S(\mathbf{r},\sigma )$ such that

$$\begin{array}{c}\hfill \begin{array}{c}S(\mathbf{r},\sigma )=\lambda S_0\left(\sigma \right)+S_1(\mathbf{r},\sigma )+\kappa S_2\left(\mathbf{r}\right)+U_{3D}\left(\mathbf{r}\right)+U_{\mathrm{vol}}\left(\mathbf{r}\right),\hfill \\ S_0\left(\sigma \right)=\frac{1}{2}\sum _{ij}\left(1-3{({\sigma }_i·{\sigma }_j)}^2\right),\hfill \\ S_1(\mathbf{r},\sigma )=\sum _{ij}{\mathsf{\Gamma }}_{ij}{\ell }_{ij}^2,{\mathsf{\Gamma }}_{ij}=\frac{1}{4\overline{N}}\sum _{\mathrm{tet}}{\gamma }_{ij}\left(\mathrm{tet}\right),{\ell }_{ij}^2={({\mathbf{r}}_i-{\mathbf{r}}_j)}^2,\hfill \\ S_2\left(\mathbf{r}\right)=\sum _i\left[1-\cos ({\varphi }_i-\pi /3)\right],\hfill \\ U_{3D}\left(\mathbf{r}\right)=\sum _{\mathrm{tet}}{U}_{3D}\left(\mathrm{tet}\right),U_{3D}\left(\mathrm{tet}\right)=\left\{\begin{array}{cc}0\hfill & (\mathrm{Vol}\left(\mathrm{tet}\right)>0)\hfill \\ \infty \hfill & \left(\mathrm{otherwise}\right)\hfill \end{array}\right.,\hfill \\ U_{\mathrm{vol}}\left(\mathbf{r}\right)=\left\{\begin{array}{cc}0\hfill & (\mid V-V_0\mid \le \Delta V)\hfill \\ \infty \hfill & \left(\mathrm{otherwise}\right)\hfill \end{array}\right..\hfill \end{array}\end{array}$$

(3)

This Hamiltonian $S(\mathbf{r},\sigma )$ does not include the energy term for the light irradiation. A detailed description of how to treat this important term is provided in the next section. Here, we describe the terms included in the basic Hamiltonian $S(\mathbf{r},\sigma )$ in more detail.

The first term $S_0$ is the Lebwohl-Lasher potential, which is always assumed for the nematic transition of LC molecules [43]. Because of this non-polar interaction of $\sigma$ in $S_0$, $\sigma$ is identified as $-\sigma$, and ${\sigma }_i·{\sigma }_j\to 1$ (${\sigma }_i·{\sigma }_j\to 0$) is expected for sufficiently large (small) $\lambda$.

The second term is a discrete version of the continuous $S_1$ in Equation (1). In the discrete $S_1$, $\overline{N}$ is given by

$$\begin{array}{c}\hfill \overline{N}=\left(1/N_B\right)\sum _{ij}{n}_{ij},\end{array}$$

(4)

where $n_{ij}$ denotes the total number of tetrahedra sharing the bond $ij$, and $N_B(={\sum }_{ij}1)$ is the total number of bonds. The symbol ${\sum }_{ij}$ in $S_1$ denotes the sum over all bonds $ij$, “tet” in ${\gamma }_{ij}\left(\mathrm{tet}\right)$ denotes all tetrahedra connected to the bond $ij$, and ${\gamma }_{ij}\left(\mathrm{tet}\right)$ for the tetrahedron in Figure 3a are given by

$$\begin{array}{c}\hfill \begin{array}{c}{\gamma }_{12}=\frac{v_{12}}{v_{13}{v}_{14}}+\frac{v_{21}}{v_{23}{v}_{24}},{\gamma }_{13}=\frac{v_{13}}{v_{12}{v}_{14}}+\frac{v_{31}}{v_{32}{v}_{34}},{\gamma }_{14}=\frac{v_{14}}{v_{12}{v}_{13}}+\frac{v_{41}}{v_{43}{v}_{42}},\hfill \\ {\gamma }_{23}=\frac{v_{23}}{v_{21}{v}_{24}}+\frac{v_{32}}{v_{31}{v}_{34}},{\gamma }_{24}=\frac{v_{24}}{v_{23}{v}_{21}}+\frac{v_{42}}{v_{41}{v}_{43}},{\gamma }_{34}=\frac{v_{34}}{v_{31}{v}_{32}}+\frac{v_{43}}{v_{41}{v}_{42}}.\hfill \end{array}\end{array}$$

(5)

We should note that ${\gamma }_{ij}={\gamma }_{ji}$ and, hence, ${\mathsf{\Gamma }}_{ij}={\mathsf{\Gamma }}_{ji}$ for all $ij$.

$\kappa S_2$ is an energy term that controls the stiffness of the material, and the coefficient $\kappa$ is the stiffness constant. In $S_2$, ${\sum }_i$ denotes the sum of internal angles i of all triangles, and ${\varphi }_i$ is the internal angle as shown in Figure 3a. For sufficiently large $\kappa$, all internal angles of triangles are expected to be ${\varphi }_i\to \pi /3$, and hence the shape of tetrahedra becomes regular. In contrast, for sufficiently small $\kappa$, the tetrahedra can considerably deviate from the regular shape.

The term $U_{3D}\left(\mathbf{r}\right)$ is the constraint potential for maintaining the positivity of each tetrahedron volume. The final term $U_{\mathrm{vol}}\left(V_0\right)$ is also the constraint potential for maintaining a constant total volume, which is given by $V_0$. This $V_0$ is determined from the model without $U_{\mathrm{vol}}\left(V_0\right)$. Details of this point are discussed in the next section.

The discrete partition function is given by

$$\begin{array}{c}\hfill Z=\sum _{\sigma }\int \prod _{i=1}^{N}d{\mathbf{r}}_i\exp [-S(\mathbf{r},\sigma ;U_{\mathrm{vol}})],\end{array}$$

(6)

where $U_{\mathrm{vol}}$ in $S(\mathbf{r},\sigma ;U_{\mathrm{vol}})$ is written to highlight the constraint on the volume. ${\sum }_{\sigma }$ denotes the sum over all possible values of $\sigma (=\{{\sigma }_1,{\sigma }_2,\cdots ,{\sigma }_N\})$, and $\int {\prod }_{i=1}^{N}d{\mathbf{r}}_i$ denotes the $3N$-dimensional multiple integrations.

4. Monte Carlo Technique and Implementation of Light Irradiation Effect

The standard Metropolis technique is used to update the variables $\mathbf{r}$ and $\sigma$ [44,45]. A new variable ${\mathbf{r}}^{\prime }=\mathbf{r}+\delta \mathbf{r}$ is accepted with the probability $\mathrm{Min}[1,\exp (-\delta S\left)\right]$, where $\delta S=S\left(\mathrm{new}\right)-S\left(\mathrm{old}\right)$. The rate of acceptance of ${\mathbf{r}}^{\prime }$ is controlled to be approximately $50\%$ by the radius $R_0$ of the sphere containing the three-dimensional random vector $\delta \mathbf{r}$. On the other hand, a new value ${\sigma }^{\prime }$ is randomly determined on the unit sphere $S^2$ using three different random numbers, and it is independent of the old value $\sigma$. To be more precise, three different uniform random numbers ${\sigma }_{x,y,z}\in (-0.5,0.5]$ are generated with the constraint ${\sigma }_x^2+{\sigma }_y^2+{\sigma }_z^2\le 0.25$. This constraint makes the point $({\sigma }_x,{\sigma }_y,{\sigma }_z)$ uniform in the sphere of radius $0.5$. Then, the vector $({\sigma }_x,{\sigma }_y,{\sigma }_z)$ is normalized such that ${\sigma }_x^2+{\sigma }_y^2+{\sigma }_z^2=1$. Thus, the distribution of this unit vector $({\sigma }_x,{\sigma }_y,{\sigma }_z)$ is expected to be uniform on the unit sphere.

We should note that $S(\mathbf{r},\sigma )$ is invariant under $\sigma \to -\sigma$ transformation, i.e., that $\sigma$ is degenerate, because of the non-polar interaction in $S_0$, as mentioned above. The rate of acceptance of ${\sigma }^{\prime }$ is uncontrollable and depends on the coefficient $\lambda$ of the energy $S_0$ in Equation (3). The total number of iterations called Monte Carlo sweeps (MCS), is typically $5\times {10}^7$ to $1\times {10}^8$ after $5\times {10}^6$ thermalization MCS.

We now describe how to implement “the light irradiation” in the model. As described in the Introduction and in Section 2, the effect of light irradiation is to reduce the orientation ordering of $\sigma$ [23,24,25,26,27,28]. To focus on this effect, we should remind ourselves that the ordering of $\sigma$ can be controlled by the parameter $\lambda$ in the FG model. Indeed, for a sufficiently large $\lambda$, the variables $\sigma$ align spontaneously into a single direction, and the system becomes nematic if no external force is applied; on the contrary, for a sufficiently small $\lambda$, $\sigma$ becomes isotropic under the same conditions for the external forces. Another technique for the implementation of this effect is to introduce temperature into the model as $1/T$ directly in the Boltzmann factor $\exp (-S/T)$ in Equation (6). However, this factor $1/T$ changes all of the coefficients $\lambda$, $\gamma (=1)$ and $\kappa$ to $\lambda /T$, $\gamma /T$ and $\kappa /T$, respectively, where $\gamma$ is the surface tension coefficient and is not included in S of Equation (3). If we assume a sufficiently high temperature such as $T\to \infty$, then the new factors become $\gamma /T\to 0$ and $\kappa /T\to 0$. However, this modification is expected to be too strong for the model because $\gamma /T\to 0$ can make the tetrahedron infinitely oblong under the constraint $U_{\mathrm{vol}}$. Therefore, it is better to use $\lambda$ to change the orientation of $\sigma$.

The problem is how to modify $\lambda$ to implement the effect of light irradiation. Recalling that the interaction ${({\sigma }_i·{\sigma }_j)}^2$ in $S_0$ is defined on the bond $ij$, we find that one possible modification of the model is to change the first term $S_0$ in S of Equation (3) as follows:

$$\begin{array}{c}\hfill S_0^{\prime }=\frac{1}{2}\sum _{ij}{\lambda }_{ij}\left(1-3{({\sigma }_i·{\sigma }_j)}^2\right),{\lambda }_{ij}=\left\{\begin{array}{cc}0\hfill & \left(\mathrm{vertex}i\mathrm{or}j\mathrm{is}\mathrm{irradiated}\right)\hfill \\ 1\hfill & \left(\mathrm{otherwise}\right)\hfill \end{array}\right.,\end{array}$$

(7)

where the irradiated vertex i or j is on the upper surface of the disk. In this definition, the irradiation of bond $ij$ is represented by $\lambda =0$ on this bond. Indeed, it is easy to see that $\lambda S_0^{\prime }$ for ${\lambda }_{ij}=0$ (⇔ bond $ij$ is irradiated) is identical to $\lambda S_0$ for $\lambda =0$, which corresponds to the disordered phase of $\sigma$, and that $\lambda S_0^{\prime }$ is identical to $\lambda S_0$ in the case of ${\lambda }_{ij}=1$ (⇔ bond $ij$ is not irradiated). In Figure 4a, we illustrate the irradiated (⇔ solid square symbol) and not irradiated (⇔ open and solid circles) vertices. The vertices inside the domain with diameter d are irradiated. The vertices that are not irradiated are divided into two groups: one contains the vertices that are connected to the irradiated vertices (⇔ open circle), and the other contains the vertices that are not directly connected to the irradiated vertices (⇔ solid circle). The domains for light irradiation are shown in Figure 4b,c, where the alignment of the nematic director is in the x-direction.

If we consider $1/\lambda$ to be the temperature, the definition of the light irradiation in Equation (7) can be rephrased as follows: the temperature is defined on the bonds through ${\lambda }_{ij}$ in Equation (7) such that ${\lambda }_{ij}=0$ (${\lambda }_{ij}=1$) corresponds to a high (low) temperature. This implies that only bond $ij$, one of the vertices of i and j is irradiated at least, corresponds to the high temperature (see Figure 4a) and that all of the remaining bonds correspond to the low temperature characterized by $1/\lambda$. Please note that $S_0^{\prime }$ is identical to $S_0$ in Equation (3) if $R=0$. Indeed, if $R=0$, then all vertices are not irradiated, and hence, ${\lambda }_{ij}=1$ for all bonds $ij$. To summarize, the introduction of ${\lambda }_{ij}$ in Equation (7) is equivalent to introduce a domain of disordered nematic directors.

We should note that the irradiated vertices inside the region of radius d fluctuate around their original positions and remain irradiated, even when they move outside the region. This point that the irradiated vertices are fixed and independent of their positions is different from the experimental situation, where the irradiated position on the surface is not always identical to the positions of molecules, which also fluctuate thermally. However, our strategy is not so poor because in real LCE samples, the total number of irradiated LC molecules inside the irradiated region is sufficiently large compared to that of the boundary. The light intensity is not always uniform, and the intensity on the boundary is weaker than that of the center of light spot. For these reasons, the modeling for the light irradiation in the FG model in Equation (7) is sufficiently simple and is considered a good approximation for the experiment.

We comment on how to obtain the constant volume $V_0$ for the constraint $U_{\mathrm{vol}}\left(V_0\right)$ in Equation (3). Here, we call the simulation with (without) $U_{\mathrm{vol}}\left(V_0\right)$ a V-fix (V-free) simulation. To obtain $V_0$, we perform V-free simulations under the same parameters $\kappa$ and $\lambda$ as in the V-fix simulations to be performed. In this V-free simulation, the light is not irradiated ($\iff d/D=0$) because the light irradiation does not change the volume of the sample LCE in the experiment [1].

5. Simulation Results

5.1. Simulation Unit and Physical Unit

Before presenting the simulation data, we comment on the simulation/physical unit, especially the length unit, although the simulation results are not compared to the experimental data, except for the shape of the LCE disks. In the simulations, $k_{B}T$, which has the unit Nm, is always fixed to 1, where $k_B$ and T are the Boltzmann constant and the temperature, respectively. Another parameter that is modified for simplicity is the lattice spacing a m, which is also fixed to $a=1$. These parameters are eliminated from the simulation data for simplicity. Therefore, all quantities with units of length should be multiplied by a if these are compared to the corresponding experimental data. The problem is whether the physical value of a is meaningful. If the physical length of a is less than the order of distance between LC molecules such as 10 $D_{\mathrm{vdw}}$, where $D_{\mathrm{vdw}}$ is the van der Waals distance, it is very hard to interpret the simulation results. Therefore, we check this point here.

The problem is that we have no physical quantity that can be compared with the experimental data, except for the shapes of the LCE disk, which are obtained as snapshots. One possible value for a is $a=1\times {10}^{-7}$ m, which was obtained in Appendix B of Ref. [21]. This value can be used for the model in this paper because the models in this paper and in Ref. [21] are the same, except for the following two facts. First, in Ref. [21], a mechanical constraint is imposed on the LCE to calculate the stress strain diagram, while it is not imposed on the model in this paper; second, the constraint for the volume $U_{\mathrm{vol}}\left(\mathbf{r}\right)$ is imposed on the model in this paper, but it is not imposed on the model in Ref. [21]. Another possible value for a can be obtained by comparing the disk size. The experimental disk size is $D_{\exp }=5$ mm, as mentioned in Section 2, while it is approximately given by $D\simeq 10$ in the simulation units, even though this D is slightly dependent on the simulation conditions, such as $\lambda$. To compare this D with $D_{\exp }$, we should multiply D by the constant a to obtain $D_{\exp }=aD$; we have

$$\begin{array}{c}\hfill a=D_{\exp }/D=5\times {10}^{-4}\mathrm{m}.\end{array}$$

(8)

This value is relatively larger than the value above $a=1\times {10}^{-7}$ m; however, it is sufficient because $\sigma$ is a coarse-grained variable. To be more precise, the variable $\sigma$ should be understood as an average direction of many LC molecules, or in other words a group of LC molecules is represented by a single $\sigma$. This coarse-graining is not special to the FG model but it is always assumed in lattice simulations for spin models for example.

Using $D_{\exp }=aD$, $D=10$ and $a\simeq {10}^{-9}\mathrm{m}$, which is approximately 10 times larger than $D_{\mathrm{vdw}}$, we have the minimum size $D_{\exp }={10}^{-8}\mathrm{m}$, to which the model is applied. This size is sufficiently small for technological applications such as actuators. In contrast, we have no limitation on the maximum size of $D_{\exp }$. Therefore, we can actually assume arbitrary size for the disk, because the lattice spacing a is an arbitrary parameter except for the constraint $a>{10}^{-9}\mathrm{m}$.

5.2. Snapshots

Figure 5a–c are snapshots of disks without light irradiation ($\iff d/D=0$), and the sizes of the lattices are (a) $D/H=12$, $N=$ 10,566, (b) $D/H=16$, $N=$ 10,980, and (c) $D/H=20$, $N=$ 10,563. Note that these $D/H$ values of the original disks are different from the actual $D/H$ values of the lattices in Figure 5a–c: these lattices are anisotropically deformed because the system is in the nematic phase, as discussed below. The parameter $\kappa$ is assumed to be $\kappa =0.5$ for all disks. We also performed simulations for $\kappa =0$ and $\kappa =0.8$; however, these are not suitable for our purpose. In fact, the disks are not always flat for $\kappa =0$, even in the absence of light irradiation, because of the large shape fluctuations. By contrast, the bending of the disks is not always sufficient for $\kappa =0.8$, even when the light is irradiated. Therefore, we consider that $\kappa =0$ is too small and that $\kappa =0.8$ is too large. Therefore, we assume that $\kappa =0.5$.

The parameter $\lambda$ should be sufficiently large because the LCE sample is in the nematic phase on the water; therefore, we assume (a) $\lambda =1.2$, (b) $\lambda =2.4$, and (c) $\lambda =1.2$ to realize the nematic configurations. We find in Figure 5a–c that the directors $\sigma$ align in specific directions, as indicated by the arrows on the figures. These directions are determined by the initial configuration for $\sigma$ in MC simulations because the assumed $\lambda$s are sufficiently large for the system to remain in the initial nematic configuration. As mentioned above, the disks deform their shape from the original circle shape to the ellipse shape, in which the diameter along the director direction is slightly larger than the diameter along the perpendicular direction. Moreover, the thickness H of the disks shrinks from their initial values, corresponding to the initial lattices in Figure 2a–c. If $\lambda$ is too small, as in the case of $\lambda =0$, the director $\sigma$ becomes random and does not align, even when the initial value of $\sigma$ is fixed to a nematic configuration.

We should emphasize that the director $\sigma$ becomes random in the isotropic phase if $\lambda$ is fixed to $\lambda =0$, for example, because the tetrahedra of the 3D disks are relatively close to the regular shape in their original configurations, as shown in Figure 2a–c. For these configurations, there is no reason for $\sigma$ to align to a direction even when this direction is given as an initial configuration. On the other hand, if the configuration of $\sigma$ is given by an initially aligned direction at sufficiently large $\lambda$, this makes the shape of the tetrahedra deformed or anisotropic. Therefore, the disks shown in Figure 5a–c obtained at sufficiently large $\lambda$ are understood to be “stressed” by the alignment of $\sigma$, even though these are in the minimum energy states.

We show the snapshots of bending disks irradiated at the central part (Figure 6a–c) and those irradiated at one side between the center and edge (Figure 6d–f), which correspond to Figure 4b,c, respectively. The size of the disk is $D/H=12$, $N=$ 10,566, which corresponds to the disk shown in Figure 2a and Figure 5a. Small spheres on the upper surface of the snapshots in Figure 6 denote the irradiated vertices and the vertices connected to the irradiated vertices (see Figure 4a). The arrows (↔) and the symbol (⊗) denote the direction of the directors $\sigma$, which are shown in the snapshots by the small bars. We should note that the directions of almost all bars are identical to the directions indicated by (↔) and (⊗). The viewing angle of Figure 4b is vertical to that of Figure 4c, and those of Figure 4e,f are vertical to each other. From these snapshots, we find that the shape of the disks is almost the same as that of the experimentally observed ones in Figure 1a,b. Indeed, it is easy to see that the bending part of the edges moves from the center to the direction of the motion of the irradiated region.

The FG model obtains the same shape as experiment because the directors at the irradiated region become randomly oriented or form an isotropic configuration. In the isotropic phase, the 3D disk is restored to its original shape and consists of almost regular tetrahedra, and the direction of the edges of the tetrahedra becomes random. In short, the shape of the tetrahedra changes from oblong to regular as a consequence of the light irradiation. Indeed, the triangles on the upper surface remain oblong along the nematic director direction if it is not irradiated, but under light irradiation, the triangle shape changes from oblong to regular. Therefore, the irradiated region, e.g., in Figure 6a, shrinks along the direction (↔, ⊗) and expands along the direction vertical to the direction (↔, ⊗).

The simulations are performed on disks of size $D/H=16$, $N=$ 10,980, and $D/H=20$, $N=$ 10,563, as shown in Figure 2b,c, respectively. The snapshots are shown in Figure 7a–f and Figure 8a–f. We find that the bending shapes are almost identical to those in Figure 6a–f. This implies that the shape deformation observed in the FG model is stable in that the shape remains unchanged with respect to a small variation in disk size, which is characterized by the ratio $D/H$.

5.3. Gaussian Curvature

In this section, we show the Gaussian curvature, which characterizes the shape of the disks bending similar to the disks shown in the previous section [40]. This quantity cannot always be compared to experimental data because no experimental data have been reported. However, a quantitative characterization of the bending by light irradiation is interesting, and the Gaussian curvature reflects surface bending. For these reasons, we calculate the curvature energy $S_G$ on the upper surface, where the light is illuminated, and on the lower surface in contact with water, and we examine the difference between these $S_G$ values.

The absolute Gaussian curvature $\left|K\right|$ is defined by $\left|K\right|=|{\partial }_1\mathbf{n}\times {\partial }_2\mathbf{n}|/|{\partial }_1\mathbf{r}\times {\partial }_2\mathbf{r}|$ [40], where ${\partial }_i\mathbf{r}$ is a tangential vector along the local coordinate axis $x_i(i=1,2)$, and $\mathbf{n}$ is a unit tangential vector (see Figure 9a). Since $|{\partial }_1\mathbf{n}\times {\partial }_2\mathbf{n}|/|{\partial }_1\mathbf{r}\times {\partial }_2\mathbf{r}|$ and $\int \sqrt{g}{d}^{2}x$ can also be written as

$$\begin{array}{c}\hfill \frac{\left|{\partial }_1\mathbf{n}\times {\partial }_2\mathbf{n}\right|}{\left|{\partial }_1\mathbf{r}\times {\partial }_2\mathbf{r}\right|}=\underset{A_{\Delta }\to 0}{lim}\frac{a_{\Delta }}{A_{\Delta }},\int \sqrt{g}{d}^{2}x=\int dA,\end{array}$$

(9)

we have

$$\begin{array}{c}\hfill \int \sqrt{g}{d}^{2}x\left|K\right|=\int dA\underset{A_{\Delta }\to 0}{lim}\frac{a_{\Delta }}{A_{\Delta }}.\end{array}$$

(10)

On the right-hand side of the first term of Equation (9), $A_{\Delta }$ denotes the area of a smooth triangle, as shown in Figure 9b, and $a_{\Delta }$ is the corresponding area on the Gauss sphere (see Figure 9c). Recalling that $\int dA$ can be replaced by ${\sum }_{\Delta }{A}_{\Delta }$ in the limit of $A_{\Delta }\to 0$, we have a discrete version of the absolute Gaussian curvature, which we write as $S_G$ such that

$$\begin{array}{c}\hfill S_G=\sum _{\Delta }{a}_{\Delta }.\end{array}$$

(11)

This $S_G$ can be calculated on the upper and lower surfaces of the 3D disks because these surfaces are considered as triangulated 2D disks. For each triangle $\Delta$, $a_{\Delta }$ is calculated using the three unit normal vectors at the three vertices (Figure 9c). The unit normal vector ${\mathbf{n}}_i$ at the vertex i is defined by

$$\begin{array}{c}\hfill {\mathbf{n}}_i=\frac{{\sum }_{j\left(i\right)}{A}_{j\left(i\right)}{\mathbf{N}}_{j\left(i\right)}}{\left|{\sum }_{j\left(i\right)}{A}_{j\left(i\right)}{\mathbf{N}}_{j\left(i\right)}\right|},\end{array}$$

(12)

where ${\mathbf{N}}_{j\left(i\right)}$ is the unit normal vector of triangle $j\left(i\right)$ connected to vertex i and $A_{j\left(i\right)}$ is its area.

First, we plot $S_G/N_T^{U,L}$ vs. $\lambda$ of the model without the constraint $U_{\mathrm{vol}}$ in Figure 10a,b to examine how $S_G$ reflects the effect of light irradiation on the upper surface. The size of the disk used for the results in Figure 10a is $D/H=8$ and $N=$ 10,346, which is relatively thick, whereas for Figure 10b, $D/H=12$ and $N=$ 10,566, as shown in Figure 2a. The light is irradiated on all of the vertices of the upper side of the disks, and this condition is expressed by $d/D=1$ in Figure 4b. The symbols $N_T^U$ and $N_T^L$ denote the total number of triangles on the upper and lower surfaces, respectively (Table 1). We find that there is a nontrivial difference between $S_G/N_T^U$ and $S_G/N_T^L$, which are $S_G$ per triangle on the upper and lower surfaces, respectively, for a large $\lambda$ region. For the disk with $D/H=8$, the difference is relatively small, even when $\lambda$ is sufficiently large. By contrast, the difference can be clearly seen for the disk with $D/H=12$ in Figure 10b.

Next, we show the results obtained for the disks whose snapshots are shown in Section 5.2, where the constraint $U_{\mathrm{vol}}$ is imposed. The results of the lattice of size $D/H=12$, $D/H=16$, and $D/H=20$, are plotted in Figure 11a,b, Figure 11c,d, and Figure 11e,f, respectively. On the upper (lower) Figure 11a,c,e, (Figure 11b,d,f), the results $S_G/N_T^U$ ($S_G/N_T^L$) calculated on the upper (lower) surface are plotted. The results at $\lambda =0.5$ (and at $\lambda =0$) are independent of the size of light spot $d/D$ and are almost independent of whether the surface is upper or lower. This is because $d/D=0$ on the light spot has no non-trivial influence on the ordering of $\sigma$, which is actually not ordered under $\lambda \le 0.5$ in the whole other region. By contrast, $S_G/N_T^U$ on the upper surface becomes dependent on the spot size $d/D$ for the large $\lambda$ region. This occurs because the light irradiation defined by $\lambda =0$ on the region of diameter d destroys the nematic order of $\sigma$ for a large $\lambda$. We should note that $S_G/N_T^L$ depends on $d/D$, even on the lower surface, where no light is illuminated. For the large $\lambda$ region, we find that $S_G/N_T^L$ at $d/D=0.5$ is non-trivially larger than $S_G/N_T^L$ at $d/D=0$. This non-trivial difference indicates that the disordered configuration of $\sigma$ on the upper surface for $d/D=0.5$ influences $\sigma$ on the lower surface.

6. Summary and Conclusions

We used MC simulations to demonstrate that the Finsler geometry (FG) model successfully reproduces the experimental results of LCE shape deformation under light irradiation reported in [1]. In the simulations, the irradiated region is given by a circle on the disk, and the center of the irradiated circle is located at the center of the disk and at the midpoint of the center and edge of the disk. For a sufficiently large $\lambda$, which is the coefficient of Lebwohl-Lasher potential $S_0$, the bending shape of the disk is almost identical to the experimentally observed shapes reported in [1].

For the FG model, the reason for the bending is simply understood as the orientation $\sigma$ becoming disordered (⇔ isotropic) on the irradiated region but ordered (⇔ nematic) on the other parts of the LCE disk. This change in the orientation order of $\sigma$ comes from the variation of $\lambda$. The role of $\lambda$ is simply to make the orientation of $\sigma$ ordered (disordered) if it is sufficiently large (small). On the other hand, in the FG model itself, the ordering of $\sigma$ and the macroscopic shape are connected by the interaction between the direction $\sigma$ of director and the position $\mathbf{r}$ of polymer. This interaction is indirectly introduced by a modification of the intrinsic geometry of materials from Euclidian to Finsler. Using the interaction introduced by the FG model technique, we find that (i) the shape of disk is determined only by the change in the director orientation and (ii) the obtained results are consistent with the experimental observations.

We should emphasize that in FG modeling, it is not necessary to delve into the details of the temperature and photoisomerization effects. In this sense, the FG model is a coarse-grained model and, hence, can be applied to many phenomena that are not always elucidated from the microscopic perspective. There is no limitation on the size of sample to which FG modeling is applied, although the sample LCE size targeted in this paper is limited in the range of a few mm.