Finsler geometry modeling and Monte Carlo study of liquid crystal elastomers under electric fields

The shape transformation of liquid crystal elastomers (LCEs) under external electric fields is studied through Monte Carlo simulations of models constructed on the basis of Finsler geometry (FG). For polydomain side-chain-type LCEs, it is well known that the external-field-driven alignment of the director is accompanied by an anisotropic shape deformation. However, the mechanism of this deformation is quantitatively still unclear in some part and should be studied further from the microscopic perspective. In this paper, we evaluate whether this shape deformation is successfully simulated, or in other words, quantitatively understood, by the FG model. The FG assumed inside the material is closely connected to the directional degrees of freedom of LC molecules and plays an essential role in the anisotropic transformation. We find that the existing experimental data on the deformations of polydomain LCEs are in good agreement with the Monte Carlo results. It is also found that experimental diagrams of strain versus external voltage of a monodomain LCE in the nematic phase are well described by the FG model.

external-field-driven shape transformations. Rubber elastic theory for polymer network with an interaction of LC director and electric field are developed [23][24][25]. In the context of statistical mechanics of LCs combined with the theory of rubber elasticity, the elongation phenomena are well understood [1][2][3][4]. These microscopic theories for LCEs are further extended by incorporating spatial symmetries associated with the director degrees of freedom [26][27][28][29][30][31][32]. From these theoretical studies, the essential part of the elongation phenomena is successfully explained.
However, the interaction between the direction of LC molecules σ(∈ S 2 : unit sphere) and the polymer is still insufficient, because the polymer position r(∈ R 3 ) is not directly included in the interaction energies in those models. Indeed, the interaction is assumed between σ and the strain variable, which is defined by a displacement field u. The displacement field is defined by the displacement from its original position r 0 such that r = r 0 + u, however, this r 0 is unknown in general. Even if r 0 is successfully assumed, the assumed condition is not always applicable to any problem.
In contrast, the position r of polymer is directly treated in the Finsler geometry (FG) [33][34][35][36] modeling, and as a consequence the mechanical response to stimuli can be calculated independent of how the orientation of σ changes [37]. If σ is isotropic, then the bulk shape does not respond to σ. However, once σ aligns into a direction spontaneously by an internal phase transition, or artificially by an external stimulus such as electric fields, then the bulk shape responds to σ. Moreover, as will be described in the next section in more detail, the treatment of the variable r is considered as an extension of that of the Doi-Edwards model for polymers [38]. In this sense, the FG model for LCE is considered as a combined model of 'statistical mechanical model of LC' and '3-dim extension of the linear polymer model', where the dynamical variable of the former part is σ and that of the latter is r. The interaction between these σ and r is introduced by Finsler metric. Note that σ is identified as −σ because of the non-polar interaction of LC molecules.
In this paper, we use an FG model to analyze two types of elongations driven by external electric fields E. The first one is observed in a polydomain cylindrical LCE of side chain type, which undergoes a transition from the random to the oriented state when a constant electric field is applied [8]. The second one is a monodomain sheet-like LCE of side chain type, the size of which changes when an alternating voltage is applied [9,10]. Monte Carlo (MC) results for the strain versus E will be compared to the reported experimental data. In the experimental data for the monodomain LCE, the variation of the strain versus E is relatively smooth. For this reason, we modify the Finsler metric such that the interaction between σ and r becomes weak, and consequently, the MC results can be compared with those experimental data. The alternating voltage in the experiment is effectively replaced by a constant voltage in the MC simulations.
Additionally, note that the real LCE is simplified in the FG model such that the shape deformation is caused only by the electrical orientation of σ, and the effects of electrostriction and the other parameters except for the electrical orientation of σ are all neglected. Moreover, the experimental condition in [8], where the LCE is surrounded by a solvent, is also simplified in the model. The reason for why the electric effect is not as simple is that the applied electric field influences not only LC molecules but also polymers and surrounding solvents [39].

Conceptual introduction to Finsler geometry modeling
Herein we briefly summarize the basic part of FG modeling and simulation results reported in [37,40] to express the reason why FG model is used to understand the material properties in condensed matter physics. First of all, FG model is simply obtained by modifying (or extending) the framework of string model. It is well known that the string theory is a model for elementary particles in subatomic scales [41]. Indeed, the string model of Polyakov is understood to be defined by a mapping r from a two-dimensional (2D) manifold M to R d where d = 3 is assumed for simplicity, and R 3 is the threedimensional (3D) Euclidean space. A smooth triangle on M is mapped to a linear triangle in R 3 for example (figures 2(a) and (b)). The assumption that the target space is R 3 is reasonable only in the scale between van der Waals distance and the diameter of earth or more larger scale, i.e. the scale for mat erial science. In this framework of 'string model', the variables to be solved or clarified are the mapping r and Riemmanian metric g ab (a,b = 1,2) of M. This modeling comes from the Hamiltonian where r(∈ R 3 ) is the image of (x 1 , x 2 ) ∈ M in equation (1), the symbol g is the determinant of g ab , the components of which are functions on M, and g ab is its inverse. Thus, in this modeling, two different manifolds r(M) and M are assumed for a description of 2D material in R 3 . The first manifold r(M) is a real curved surface, which is a manifold with a set of local coordinates and the induced metric of R 3 , and the second M is an abstract manifold with a set of local coordinates, corresponding to those for r(M), and the metrices g ab [42]. These g ab are not for the real surface denoted by r(M), but these are metrices of M. The Euclidean metric δ ab , which is constant, is an example of g ab (figure 2(c)). We should note that if g ab is the induced metric then this abstract manifold M is logically identified with r(M). In this case, S 1 of equation (2) becomes the minimal 'energy', which is identical with the area of this real surface [43]. In this sense, g ab originally represents 'geometry fluctuations' of the real surface around the induced metric in the context of string theory. Note that it is possible, at least mathematically, to assume that the metric g ab is not identical to the induced metric for the surface in R 3 . However, r(M) should be regarded as a real material in R 3 , and therefore it is reasonable to assume that the variable g ab is introduced for M.
As we immediately understand, this framework seems not always suitable for describing materials in R 3 because of the extra degrees of freedom g ab for M. Indeed, the information r and its local coordinates with the induced metric are sufficient for describing the shape of materials in R 3 , and there is no room to use the information g ab on abstract M in the modeling of any material shapes in R 3 . This extra degree of freedom simply comes from a mathematical consideration that r is regarded as a 'mapping' from a 2D manifold M to R 3 , or in other words, g ab used in the Hamiltonian such as the one in equation (2) can be naturally extended from a metric of the surface in R 3 to a metric of the abstract surface M.
However, as pointed out in [37], the above mentioned degrees of freedom g ab of M plays an essential role in describing anisotropic phenomena in membranes. If g ab is extended from Riemmanian to Finsler metric, this g ab of M turns from 'unnecessary' to 'necessary' information and is able to give a framework for describing the mechanical properties, especially anisotropic properties, of materials. Indeed, 'differential interactions' (or 'differential') on the surface such as ∂ a r · ∂ a r (or ∂ a r) of S 1 in equation (2) are locally modified by the Finsler metric g ab (σ), which is a function of internal degrees of freedom σ i such as the direction of LC molecule i. This modification of interaction by g ab (σ) changes the mechanical property of material both locally and direction-dependently. If this direction dependence of interaction is simply occurs only locally and randomly, no change is expected in the macroscopic observation of physical quantities, because the mean value of physical quantities over the surface is independent of the local coordinate used to describe the Hamiltonian [37,42]. Nevertheless, if the direction dependence of interaction appears globally on materials, the direction dependence of mechanical properties can be observed, even though the mean value of physical quantities remains unchanged. This point will be discussed further in the following section.
An important point to note is that two different length scales are assumed for a path on the surface r(M) (figure 2(d)). The first is the real length L of the path defined by the sum of the Euclidean length of the tangential vectors such as L = t1 t0 dt (dr/dt) 2 , and this length is necessary for all physical descriptions of material. The second is the Finsler length s of the path and is defined by s = t1 t0 (dr/dt) 2 /|v|dt, where |v| is the unit Finsler length given by a vector v [37]. This is a general definition of Finsler length on a path of r(M). This leads to the discrete Finsler [37]. This metric is considered to be a Finsler metric of the triangle 123 on the manifold r(M) in R 3 , because it is obtained by modifying the induced metric ( This metric should be understood as a Finsler metric of the triangle 123 on M (figure 2(a)), because it is obtained by modifying the Euclidean metric δ ab on M. The role of the Finsler metric is to define a modified 'differential interaction' between r as mentioned above. Moreover, if the directional degrees of freedom σ of LC is regarded as v for example, the interaction of σ and r is introduced by g ab because v ij is defined by using σ and r in this case. The Finsler metric in equation (3) is assumed in the calculations that have been reported [44], because it is relatively simple. The FG modeling technique can be extended for 3D materials and was recently applied to the soft elasticity and elongation phenomena, and the simulation results are consistent with the reported experimental data [40]. The soft elasticity is characterized by a plateau in the stress-strain diagram of LCE, where the direction of LC molecules changes almost π/2 and aligns to the direction along which a tensile force is applied as mentioned in the Introduction. This directional degree of freedom of LC plays an important role in the plateau of the diagram, and the experimentally observed diagram is reproduced by the 3D FG model. The numerical results are qualitatively in good agreement with experimental data [40]. The shape of LCE becomes anisotropic at low temperature and restores to the original isotropic shape at high temper ature. This elongation phenomenon under temperature variation is also reproduced, and the results are qualitatively in good agreement with the reported experimental data, although the temper ature is not explicitly introduced in the simulations of [40].
To summarize, in the FG technique, the geometry inside LCE is deformed by modifying the metric tensor from Euclid to Finsler in the above mentioned way, instead of introducing the interaction between LCs and polymers [40]. As a result of this geometry deformation by Finsler metric, an interaction between the variable σ and the polymer position r is automatically introduced, and hence, the modeling technique is completely new compared to the existing technique, in which the interactions between dynamical variables are directly introduced in the Hamiltonian [37]. For this reason, the FG technique should be evaluated in further detail especially in the field of material science or (soft) condensed matter physics. In fact, it has been reported that such a technique has applications mainly in the field of cosmology and/or quantization theories [45][46][47].

3D cylinders discretized by tetrahedra
In the experiment reported in [8], cylindrical LCEs were studied. Therefore, in our numerical study, a cylinder such as the one shown in figure 3(a), discretized by tetrahedra, is used to define the discrete Hamiltonian, which will be described in the following section. The size of the cylinder is given by (N, N B , N T , N tet ) = (19 784, 140 507, 238 291, 117 567), where N, N B , N T , and N tet are the total numbers of vertices, edges (or bonds), triangles, and tetrahedrons, respectively. The ratio L/D of length L and diameter D is L/D = 1, and the Voronoi tessellation technique is used to construct the 3D lattices [49][50][51]. In this construction, the vertices are randomly distributed not only on the surface but also inside the surface under the constraint that the mean edge length on the surface is almost identical to that inside the surface.
In the experiment in [9,10], thin and rectangular LCE plates were used. However, we use a thin cylinder, as shown in figure 3(b), for simplicity to simulate the results of the experiment. The size of this cylinder is given by (N, N B , N T , N tet ) = (10 346, 69 964, 116 041, 56 422), and the ratio of height/diameter is H/D = 1/8. This thin cylinder is constructed under the same constraint on the mean edge length as that in the oblong cylinder in figure 3(a). We note that the Euler numbers N − N B + N T − N tet (= 1) of these two cylinders are identical to that of a tetrahedron, which is characterized by (N, N B , N T , N tet ) = (4, 6, 4, 1), because the cylinders are topologically identical to a tetrahedron.
We describe the lattice construction procedure in more detail. The oblong and thin cylinders are constructed by the same procedure. As the first step, the total number of vertices on the upper and lower disks and that on the side surface of cylinder are determined such that the area per vertex is constant. The total number of internal vertices is also determined such that the volume per vertex is constant. These constants are approximately evaluated by assuming that the cylinder is composed of regular tetrahedron and that the triangle of this tetrahedron appears on the surface. The next step is to generate the lattice vertices randomly on and inside the surfaces under a constraint that the minimum distance between the new vertex and the existing vertices is larger than a given distance. If this constraint for the new vertex is not satisfied within a finite number of vertex generations, then the lattice construction fails. In this case, the constraints such as the minimum distance between vertices, minimum distance between the surface and the internal vertices, the total numbers of internal vertices, etc are slightly changed. If this vertex generation process is successful, the next step is to construct the tetrahedra by Voronoi tessellation technique [49][50][51]. If this construction process also fails before it completes, we further change or loosen the constraints and repeat the same procedure from the first step. As we will see from the FG model described in the following subsection, the lattice deformation is closely connected to the internal degrees of freedom σ, and therefore the initial lattice seems crucial for the convergence of results. For this reason, we plot the normalized histogram h(L) of the distribution of the normalized bond length L in figure 2(c). We find that the histograms h(L) of oblong cylinders (upper part of figure 2(c)) are relatively uniform, while the distribution of h(L) for the thin cylinder (lower part of figure 2(c)) is relatively wide, although h(L) → 0 is satisfied for very small L .

Finsler geometry model
The details of the 3D FG model were reported in [40]. In this paper, we show the outline of the model in a self-contained manner because the model is slightly modified and the Hamiltonian includes a new term for the external electric field.
The continuous Hamiltonian is given by where r(∈ R 3 ) is the position of 3D material and is understood as a function of the parameters x a (a = 1,2,3); the symbol g is the determinant of 3 × 3 matrix g ab , the components of which are the metric functions on the parameter space; and g ab is its inverse. We note that the Hamiltonian S 1 in equation (4) is defined by the 3D integration √ gd 3 x over the 3D material, and therefore it is different from S 1 in equation (2), which is defined by the 2D integration √ gd 2 x over the 2D surface.
Three-dimensional cylindrical bodies, which are partitioned into tetrahedra (figures 3(a) and (b)), are used to define the discrete Hamiltonian as mentioned above. On the tetrahedron with vertices 1234 (figure 4(a)), the discrete metric g ab is given by [37,40] In this discrete metric, v ij is defined as with the vertex position r i and the variable σ i at the vertex i. Note that this v ij is the Finsler length of the model and defines an interaction between σ and r [37]. The parameter v 0 in v ij is the cutoff in the model of [40,44]; however, in this paper, it plays the role not only of the cutoff but also of the parameter that controls the interaction strength.
If v 0 is sufficiently large compared to the first term |t ij · σ i |, then no anisotropy is introduced in the model. Indeed, if v ij v 0 and |t ij · σ i | is neglected, then g ab is proportional to the Euclidean metric δ ab . By replacing the integral with the finite sum (4), and including several additional terms, we have the Hamiltonian: (8) The first term λS 0 is the Lebwohl-Lasher potential for LC molecules with the interaction coefficient λ [48]. The second term S 1 is the discrete Gaussian bond potential corresponding to equation (4), and the coefficient N is given bȳ where n ij is the total number of tetrahedra sharing the bond ij and N B (= ij 1) is the total number of bonds. The coef- 34 , 34 , where v ij are defined on the tetrahedron in figure 4(a) and tet denotes the sum over tetrahedra, each of which shares the bond ij. The surface tension coefficient γ of S 1 is fixed to γ = 1 and eliminated from S in equation (8).
Note that the interaction between σ and r is actually implemented in S 1 via the coefficients γ ij in equation (10). These γ ij are considered as an effective tension, which depends on the position and direction inside the 3D body, including the surface. If γ ij is constant in the sense that it only depends on the direction, then it can be called anisotropic surface tension. Moreover, the magnitude of this effective tension is controlled by v 0 in equation (6). The reason for why the interaction strength varies with the value of v 0 is because γ ij in equation (10) is a rational function of v ij and strongly influenced by small v ij in the denominator in general. For these reasons, FG modeling is basic in the sense that it explains where the anisotropic interaction in S 1 comes from. This point should be explained in further detail. From the so-called scale invariance of the partition function Z, which will be described below, under the scale change r → βr, the mean value of ij of the FG model remains constant [52]. This constant is identical with S 1 = ij 2 ij of the canonical model defined by S of equation (8) without the terms λS 0 and αS E , where S 1 = ij 2 ij . As a consequence, if Γ ij = 1 for all ij in the FG model, then there appears no difference in the material shape between the two models, and therefore, no anisotropy is expected in the shape of materials in this case. This is also true in the case where Γ ij changes randomly, locally and direction dependently in the FG model [37]. In contrast, if Γ ij changes direction dependently and globally, an anisotropy is expected in the shape of material as mentioned in section 1.2. For example, if Γ ij in X direction inside the material is smaller than Γ ij in Y and Z directions, then the corre sponding 2 ij in X direction becomes larger than 2 ij in Y and Z directions, where (X, Y, Z) denote the canonical axes of R 3 , and 'Γ ij (or 2 ij ) in X direction' implies that the direction of the bond ij is almost or nearly parallel to X direction. We should emphasize that the direction dependence of Γ ij comes from the direction dependent v ij in equation (6), and this direction dependence of v ij is expected if σ i aligns to a direction along which an external electric field is applied.
We should note that the Finsler metric in equation (5) or the Finsler length v ij in equation (6) simply plays a role in modifying S 1 from ij 2 ij to ij Γ ij 2 ij . Note also that this modified interaction term ij Γ ij 2 ij is considered as a 3D extension of the corresponding term in the model of Doi-Edwards for polymers, because ij 2 ij is simply an extension of the model for membranes, which is considered as a 2D extension of the Doi-Edwards model [38].
The third term κS 2 is the energy for the stiffness of the tetrahedron with the stiffness constant κ, and i is the sum of all internal angles φ i of triangles ( figure 4(a)). This term plays a role of maintaining the shape of the 3D body just like a bending energy in the surface models [52][53][54][55][56][57][58][59]. The fourth term −αS E is the interaction energy between σ and the external electric field E. Since the simulation unit is different from the physical unit, as we will discuss below, the coefficient α corresponds but is not identical to 0 ∆ , where 0 is the dielectric constant of free space and ∆ is the dielectric anisotropy: In this expression, and ⊥ are the dielectric constants along and perpendicular to the direction of σ. In this paper, we assume both positive ∆ > 0 and negative ∆ < 0 values as in the experiments in [8].
The fifth term U 3D prevents the tetrahedron volume from being negative. The positivity of the volume is a natural constraint for the model. The final term U vol (V 0 ) is the constraint for the volume V: Since the volumes of LCEs remain unchanged in the experiments in [8][9][10], we simply assume this constraint for V [60][61][62]. This constraint is imposed only on the vertices on the surface of the cylinder because no change is made on the volume by any movement of the vertices inside the surface under the constraint of U 3D (r). In the expression of equation (12), ∆V is given by the mean tetrahedron volume ∆V = tet V(tet)/N tet . A constant volume V 0 is obtained in the simulations without the constraint U vol (V 0 ) under E = 0 such that V 0 = tet V(tet) using the same parameters, which will be assumed in the production simulations. In short, the constant V 0 is the equilibrium volume of the model without the constraint U vol (V 0 ). The discrete partition function is defined as where U vol in S(r, σ; U vol ) is the constraint on the volume as described above.

Monte Carlo technique
Variables r(∈ R 3 ) and σ(∈ S 2 ) are both updated by the standard Metropolis technique [63,64]. The first one is changed as r → r = r + δr with a small random vector δr, and the second one is changed as σ → σ , where the new variable σ is randomly defined to be independent of the old variable σ. The oblong cylinder-shaped LCE studied in [8] is surrounded by a gel inside a rectangular box, and therefore, the cylinder is loosely prohibited from changing its oblong direction in the simulations. Indeed, the center of mass of vertices on the upper (lower) surface (=disk) is fixed inside a small circle, the center of which is at the position (x, z) = (0, 0), on each disk during the MC updates such that the cylinder does not move on one side (see figure 5(a)). The radius D 0 of this small circle is fixed to D 0 = 6 , where is the mean bond length of the cylinder of volume V 0 . The center of mass along the y axis is fixed to the position y = 0.
The LCE in [8] is in the polydomain phase, and the LC molecules are locally anisotropic but globally isotropic for E = 0. Therefore, this type of LCE is considered to be an isotropic LCE, which can be represented by our model in the limit of or close to λ → 0. Indeed, by replacing the locally aligned LCs by a single variable σ (figures 6(a) and (b)), we obtain an LCE model for the polydomain LCE. This is a coarse-grained picture of the polydomain LCE.
In the experiments reported in [9,10], this LCE plate is sandwiched by two parallel electrodes with a small but sufficient space such that the size change of the plate is not influenced by the electrodes. For this reason, a disk of thickness H is placed inside the two parallel plates separated by h in the simulations, and the hard wall potential is assumed for the plates (figures 5(b)). The separation h is fixed to h 2H, where the width H of the disk varies depending on the simulation parameters, including E. This hard wall potential is not included in S of equation (8). The center of mass along the horizontal direction is fixed to its origin.
The LC molecules inside the LCE in [9,10] are in the nematic phase, and the LCs are forced to align parallel to the plate when E = 0. In this LCE, the internal network structure of polymers is anisotropic. This is in sharp contrast to the LCE in the polydomain phase in [8] described above. For this reason, to simulate the LCE in [9,10] we fix the parameter λ to be sufficiently large such that the variable σ aligns to the predicted direction which is given by the initial configuration of σ. In addition to this, a small electric field E 0 is applied to this direction. This will be explained in the presentation section in more detail. Moreover, the parameter v 0 of v ij in equation (6) is fixed to v 0 = 0.1, which is quite larger than the cutoff δ(= 0.001). As mentioned in section 2.2, this increment of v 0 from v 0 = δ weakens the strength of the interaction between σ and r. The reason why the interaction strength should be weakened is that if the strength is not weakened, then the obtained results deviate from the experimental data, as we will discuss in the presentation section in detail.

Simulation unit and physical unit
The lattice spacing a is assumed to be a = 1 in the simulations for simplicity. The unit of energy is also simplified such that k B T = 1. These modifications define the simulation unit. Meanwhile, the lattice spacing a can be restored to a physical value if the energy scale is changed to the physical one. The problem is whether the actual value of a in this case is reasonable. In this subsection, we comment on the formula for estimating a when the unit of energy for the external field is changed from the simulation unit to physical unit. This change of units is necessary for comparing the simulation results and the experimental data.
Note that the unit of energy αS E in equation (8) is . Therefore, the simulation data α(σ i · E) 2 k B T(Nm) should be divided by a parameter cubed that has a unit of (m 3 ). This parameter is given by the lattice spacing a. Thus, using this a, we have From this, we obtain where we use the fact that E = 10 −6 E exp and the data 0 = 8.85 × 10 −12 , k B T = 4 × 10 −21 . Note that the unit of E exp is (10 6 V m −1 ) in the experiments [8]. From this equation in equation (15), we estimate a using the input data α and the material parameter ∆ . The dielectric anisotropy ∆ in equation (11) of LCN-P-100/SP in [8] is given by ∆ = 11.5 [65], and the constant α assumed in the data analyses is of the order from 10 −2 to 10 1 , and hence, α/∆ is larger than 10 −3 at least. Even for this case of smallest α/∆ , the estimated a is still sufficiently larger than the typical van der Waals distance D VdW (∼1 × 10 −10 (m)). This result is sufficient for the models to be called 'coarse-grained models'. In fact, if a was smaller than D VdW , then a (∼mean bond length) becomes a subatomic scale, and in this case, the meaning of bond length and even the position r itself are hardly interpreted.  We should comment here on how to determine α, which is positive for simplicity. The parameter α is actually not used in the simulation, but it is determined in the data analyses, even though it is included in S of equation (8). This α comes from the above-mentioned relation E = 10 −6 E exp . Indeed, by replacing √ αE with E sim , we have αE 2 = E 2 sim , which corresponds to the energy αS E . Using this definition E sim = √ αE, where E sim is the actual value assumed in the simulations, we have E sim / √ α = 10 −6 E exp for the data analyses. In the simulations, we simply use E sim without α for the energy αS E in equation (8). In other words, α is determined such that the strains obtained at E sim / √ α are identical to the experimental ones at E exp . If these two strains, where one is the simulation data at E sim / √ α and the other is the experimental data at E exp , are identical for all regions of E sim / √ α and E exp , then E sim / √ α can be identified as E exp . This is the meaning of the Note that the unit of α is (C 2 m N −1 ), differing from that of 0 ∆ , as mentioned in section 2.2. The unit of E in equation (8) of the model is same as that of E exp in the experiment except the factor 10 6 .
To summarize, we have the relations We should note that E sim is the actual value assumed in the simulations, where α is not used. The parameter α is determined in the data analyses so that the strains calculated at E sim / |α| are comparable with the experimental strains observed at 10 −6 E exp . This will be shown in the following section.
In the case of monodomain LCE reported in [9,10], the factor 10 −6 in the relation between E and E exp is changed to E = 10 −3 E exp , because the voltage in the experiment is relatively lower than of the polydomain LCE in [8]. In the case of E = 10 −3 E exp , the lattice spacing a in equation (15) is determined by a 3 = 4.07 × 10 −19 (α/∆ ), and this a is larger than that determined by equation (15), and hence a is larger than van der Waals distance D VdW .

Polydomain LCE
The strains ε µ (µ = x, y, z) are defined as where D x and D z are the diameters of the cylindrical region of width L 0 in the central part, where the electric field E is applied to the x axis (figures 7(a)-(d)): D x (as well as D z ) is calculated by D x = x max − x min , where x max and x min are the maximum and minimum x of vertices on the thin cylindrical region. The thickness L 0 is defined as L 0 = 6 , where is the mean bond length. The same technique is used for calculating D y , which is the length of the oblong cylindrical region of diameter D 0 at the central part of the cylinder ( figure 7(c)). The diameter D 0 is identical to L 0 and defined as D 0 = 6 , as mentioned in section 2.3. The D 0 µ in equation (17) is given by the D µ of the cylinder of size D y = D, where the strain becomes zero. The reason why the region for calculating D µ is fixed to be sufficiently smaller than D is to reduce the errors that are suspected when the cylinder leans from the z axis even if its angle is small enough (figure 7(d)).
We use two different lattices of size (N, N B , N T , N tet ) = (10 492, 73 589, 124 139, 61 041) and (19 784 , 140 507, 238 291, 117 567) for the simulations in this subsection for the polydomain LCE in [8]. We find that the size dependence of the results is not large. However, the linear size is important for the size dependence, and hence, the size dependence will be commented in the final part of this section. On the other hand, it is also interesting to see the shape dependence. To see the dependence of the results on the ratio L/D, we use the lattice of L/D = 2 for calculating ε µ and find that the results obtained with the same parameters deviate only slightly or remain almost unchanged from the results obtained from the lattices of L/D = 1.
The two different experimental data (Exp) in [8]  where v 0 = 0.1. The reason why we plot such negative results is to emphasize that the parameter v 0 controls the strength of the interaction between the variables σ and r, as described in sections 2.2 and 2.3. This dependence of ε µ on v 0 reflects one  (17), (c) an oblong cylindrical region for calculating D y , and (d) a small deviation of the symmetric axis from the y axis. The field direction is fixed to the x axis such that E = (E, 0, 0). of the interesting properties of the FG model: the interaction strength is reflected in the shape deformation.
The coefficient λ is fixed to λ = 0, and hence, the variable σ is expected to be completely random for E = 0. This isotropic phase in the model corresponds to the nematic networks with zero external field in a real polydomain LCE. The rigidity κ, which is the coefficient of S 2 in equation (8) It is expected that the strain ε µ depends not only on the parameter v 0 but also on the rigidity κ. The dependence of ε µ on κ is simply understood as the shape deformation coming from the deformation of the tetrahedron, of which deformation is controlled by the rigidity κ. In figures 9(a)-(d), v 0 (= 0.001) is increased to v 0 = 0.1 and v 0 = 0.04, and κ is reduced from those assumed in figures 8(a) and (c). The lattice size is N = 10 492. We find that the results are almost consistent with the Exp data for both positive and negative dielectric anisotropy.
The tensor order parameter for σ is defined as The eigenvalues Σ 1,2,3 are plotted in figures 10(a) and (b), which correspond to the configurations of ε µ versus E and E exp in figures 8(a) and 9(a) and those in figures 8(c) and 9(c), respectively. The largest eigenvalue Σ 1 is considered to be identical to the scalar order parameter for µ = x because σ is expected to align along the x axis, which is the field direction. We find that Σ i (open symbols) for figures 8(a) and (c) are almost identical to those (solid symbols) for figures 9(a) and (c). This result implies that the internal configuration of σ is closely connected to the shape of the materials, and therefore it is consistent with the conclusion of [8]. This elongation is called affine deformation in the neo-classical theory of nematic elastomers [1]. In the left (right) column, the side (front) views are shown. Since the parameter λ is fixed to λ = 0 in the simulations, the configuration of σ is expected to be isotropic if E = 0. This is confirmed from figure 11(a): the σ is random, and the front shape is a circle. However, if E is increased, the configuration of σ changes from isotropic to anisotropic, and the direction of σ depends on whether the aniso tropy constant α is positive ( figure 11(b)) or negative (figure 11(c)), where the electric field E is applied to the x axis. In the case of positive α ( figure 11(b)), the σ aligns along the x axis, and the   front shape changes from circle to ellipse ( figure 11(b) right). In contrast, for negative α ( figure 11(c)), the σ appears to align along the vertical direction (figure 11(c) right), and from the side view (figure 11(c) left), we can see that the σ aligns to nonuniform directions vertical to the E direction, as expected. The scales of the snapshots are the same, and therefore, we can also observe the variation of the size in each direction. Indeed, the size of the cylinder for positive (negative) α is increased (decreased) in the x axis direction and decreased (increased) in the other two directions (figures 11(b) and (c)).

Monodomain LCE
In the experiment for a monodomain LCE in [9,10], the authors focused on the shape deformation due to the orientation change of directors by an applied electric field along the z direction ( figure 12). In the material that they studied, the directors align to the direction (1, 0, 0) for the zero field. During the directional variation of the director from (1,0,0) to (0, 0, 1), the y component does not change from zero. For this reason, the size of the LCE in the y direction remains almost constant during this orientation change in the experiment. We will evaluate whether this phenomenon together with the elongation into the other directions are observed in the FG model.
For this purpose, we assume in the model a constant external field (E 0 ,0,0) in addition to the variable field (0, 0, E), and therefore, we have (figure 12) To obtain the strains ε x,y in equation (17) under this E, we calculate the diameters D x,y using the same technique as that for the oblong cylinder in figures 7(a)-(d). Indeed, we calculate D x by D x = x max − x min , where x max and x min are the maximum and minimum x of vertices on the thin cylinder. The D y is also calculated by the same technique. We examine two different conditions for σ in the simulations: the nematic and isotropic phases. In the nematic phase, the variable σ spontaneously aligns even when E is zero. This nematic phase is consistent with the experimental situation described above. The initial σ aligns along (1, 0, 0) in both conditions because of the field (E 0 ,0,0). Note that the reason why we assume a nonzero E 0 in E = (E 0 , 0, E) is because the initial configuration of σ would be random if E 0 = 0 at least in the isotropic condition, realized by λ = 0, in the model. First, we assume the nematic phase, which is obtained by increasing λ to a sufficiently large value λ(> λ c ), where λ c 0.26 is the critical value between the isotropic and anisotropic phases [40]. In this case, a relatively small E 0 in E = (E 0 , 0, 0) makes the configuration of σ close to σ = (1, 0, 0). Indeed, due to the nematic nature, the variable σ aligns to σ = (1, 0, 0) as the initial configuration even for small E 0 , and σ remains unchanged from this configuration when the MC update is processed. If an external E is applied, then the direction of σ is expected to vary from (1, 0, 0) to (0, 0, 1) for sufficiently large E. We also expect that this (0, 0, 1) returns to the original (1, 0, 0) if the field E is switched off. This is also the reason why a small, non-zero, finite E 0 is necessary. It is also possible to assume a sufficiently large E 0 . However, the variation of ε x versus E becomes abrupt if E 0 is too large in the nematic phase for λ > λ c , and this abrupt change of σ is contradictory to the experimental fact in [9,10], where ε x smoothly varies with respect to E. Therefore, only a small value can be assumed for E 0 .  In figures 13(a) and (b), we show the simulation results with two different Exp data from [9,10]. The simulation results are obtained with the parameters λ = 0.6, v 0 = 0.1, and κ = 0.37 in figure 13(a) and λ = 0.6, v 0 = 0.04, and κ = 0.34 in figure 13(b). As shown in the figures, the simulation results are almost in good agreement with the Exp data. Note that many combinations for λ, v 0 , and κ close to these values are possible for producing results that are consistent with the Exp data.
We also performed the simulations in the isotropic phase, as mentioned above, where λ is fixed to λ = 0, which makes σ random for E = (0, 0, 0). The parameter v 0 in equation (6) is varied in the range from the cutoff v 0 = δ to v 0 = 0.6. In this isotropic phase, the configuration of σ becomes close to σ = (1, 0, 0) if E 0 is sufficiently large in E = (E 0 , 0, 0). However, the simulation results ε x , which are not plotted, vary rapidly with respect to E and deviate from the experimental data. We examined many combinations for the parameters v 0 and κ under λ = 0. However, we failed to obtain data that smoothly vary as those in figures 13(a) and (b), at least for the assumed parameter values. This result indicates that the nematic phase is more suitable for obtaining results that are consistent with the experimental data. This implies that the LCE elongations by electric fields are well understood in the context of FG modeling.
The order parameters Σ 1,2,3 and M z , corresponding to the configurations for ε x,y in figures 13(a) and (b), are plotted in figures 13(c) and (d), where M z is defined in equation (20) for µ = z. The Σ 1,2,3 remain constant in this case, because the nematic director is just rotating. This is in sharp contrast to the cases in the previous subsection (see also figures 10(a) and (b)). Note that the fact that Σ 1,2,3 are constant is also true even for E 0 = 0 because σ always aligns along a spontaneously determined axis in the nematic phase. The parameter M z does not remain constant and varies against E (figures 13(c) and (d)). The reason why M z starts with a negative value at E = 0 is because σ is almost close to (1,0,0) if E = 0.
Snapshots of the thin cylinder shown in figures 14(a) and (b) are obtained from the simulations for figure 13(a). The system is in the nematic phase because λ = 0.6 is larger than λ c , and a small E 0 is applied in the x direction; for these reasons, the σ aligns along the x axis when E = 0 ( figure 14(a)). If E is increased in the z axis, then the direction of σ becomes almost parallel to the z axis, as confirmed in figure 14(b).
Finally we should comment on the dependence of the results on the lattice size N. As we see from the snapshots of initial lattice in figure 3(a) or (b), the size of tetrahedron or the edge length of triangle is not always negligible compared to the cylinder size L or D. In addition to this fact, taking into account the thermal fluctuations of tetrahedrons, we consider that the total number N up to N = 2 × 10 4 is very small. On the other hand, the lattice construction technique described in section 2.1 is very time consuming if the total number N is increased, and for this reason the size N is limited up to N = 2 × 10 4 in this paper. Nevertheless, the numerical results are relatively in good agreement with experimental data. To obtain the dependence of the results on the size N correctly, we should use lattices which are at least 10 times larger in the linear scale L and D than the lattice in this paper [66]. Since N is proportional to the volume of lattice, N must be increased by the third order of the linear scale. For this reason, we consider that the dependence of the results on the lattice size is out of the scope of this paper and remains to be studied. It is  also interesting to study the dependence of the results on the uniformity or regularity of tetrahedron.

Summary and conclusion
We have shown that the elongation phenomena of LCE under electric fields are understood in the context of Finsler geometry (FG) modeling. In the FG model, the director (or LC molecule) degrees of freedom are represented by the variable σ, which is identified with −σ, and the σ i interacts with the variable r i , which is considered as the position of the polymer. This interaction is independent of how σ is oriented. Therefore, it is natural to expect that the LCE elongation through the electric orientation of directors can be understood in the context of the FG model. This expectation is confirmed in this paper using reported experimental data of both polydomain and monodomain LCEs.