A Rod Theory for Liquid Crystalline Elastomers

We derive a general constitutive model for nematic liquid crystalline rods. Our approach consists in reducing the three-dimensional strain-energy density of a nematic cylindrical structure to a one-dimensional energy of a nematic rod. The reduced one-dimensional model connects directly the optothermal stimulation to the generation of intrinsic curvature, extension, torsion, and twist, and is applicable to a wide range of liquid crystalline rods subject to external stimuli and mechanical loads. For illustration, we obtain the shape of a clamped rod under uniform illumination, and compute the instability of an illuminated rod under tensile load. This general framework can be used to determine the shape and instabilities of nematic rods with different cross-sections or different alignment of the nematic field.


Introduction
Nematic elastic rods are slender, flexible structures made of liquid crystal elastomers (LCEs) [35,48,65,71,72]. Due to their stimuli-responsive material properties, they can convert light or heat into mechanical work, without the need for batteries, electric wires or gears [58,67].
Rods are filamentary structures capable of two main reversible deformations causing large displacements, namely, bending and twist. Classically, they are modeled as long, thin elastic bodies defined by a central curve and acted upon by external loads [3,24,36,50]. An extension to anelastic rods where non-elastic deformations and changes in microstructure occur was considered in [19]. A generic one-dimensional strain energy for rod models derived from large-strain finite elasticity was recently proposed in [4] where many related models were also reviewed.
For a nematic rod, significant displacements can be achieved by natural deformations caused by external stimuli. A fundamental problem is then to establish how such spontaneous deformation may lead to the generation of curvatures. Experimental observations for A. Goriely 1 Mathematical Institute, University of Oxford, Oxford, UK 2 School of Mathematics, Cardiff University, Cardiff, UK nematic elastomer rods that bend upon exposure to UV light were reported in [35,72]. Rods with an in-printed helical director field under a remotely controlled rotation were presented in [48]. Light-controlled bending and twisting of nematic rods were investigated in [65,71]. Theoretically, three-dimensional nematic cylinders under combined stretch and torsion were analyzed in [17], and a nonlinear beam model for photoresponsive structures was proposed in [30].
For ideal monodomain nematic elastomers, with the liquid crystal mesogens uniaxially aligned throughout the material, a simple continuum model is provided by the neoclassical strain-energy function proposed in [5,62,63]. This is a phenomenological model based on the molecular network theory of rubber [54]. The constitutive parameters appearing in the neo-Hookean-type strain energy can be obtained through statistical averaging at microscopic scale or derived from macroscopic shape changes at small strain [60,61]. Here, we adopt the neoclassical strain-energy function, and assume the isotropic phase at high temperature as the reference configuration [8,[11][12][13][14], instead of the cross-linking nematic phase [2,5,56,[62][63][64]74]. We exploit theoretically the multiplicative decomposition of the deformation gradient from the reference configuration to the current configuration into an elastic distortion followed by a natural shape change [20,[41][42][43][44]. This multiplicative decomposition is similar to those found in the constitutive theories of thermoelasticity, elastoplasticity, and morphoelasticity [19,37] (see [18,53] as well), but it is also different in the sense that the stress-free geometric change is superposed on the elastic deformation, which is directly applied to the reference state.
We consider the macroscopic deformation of a nematic rod as the product of a spontaneous deformation inducing curvatures and a finite elastic deformation preserving cylindrical symmetry. We are interested in answering the following main questions: (i) If the shape of a filamentary nematic elastomer structure is known, what is the internal nematic stretch field created by an external stimulus? (ii) Given the local nematic stretch field of a filamentary structure under a stimulus, what are the intrinsic curvatures and what is the shape of a rod subject to external loads and boundary constraint?
As we show, these questions are particularly difficult to resolve since a change of shape also changes the position of the rod with respect to its stimulus. In particular cases, we can decouple the two problems but in general, the two questions cannot be answered in isolation. Our approach consists in reducing the three-dimensional (3D) strain-energy density function of a nematic cylindrical structure to the one-dimensional (1D) energy of a nematic rod.
To achieve this, we adapt the strategy developed for morphoelastic rods in [46] (see also [34,45,47,49]), and apply the resulting model to specific scenarios leading to the formation of rods with non-trivial intrinsic curvatures. In particular, we consider rods with square or circular cross-section, and a nematic director uniformly aligned at a given angle relative to the longitudinal axis [65,71], or a helical director field [48]. In Sect. 2, the 3D neoclassical strain-energy function for ideal nematic elastomers is briefly reviewed. The 1D nematic rod model is derived in Sect. 3. Examples of light-induced deformations are presented in Sect. 4.

A Continuum Model for Ideal Nematic Elastomers
To describe an ideal nematic LCE, we assume the following strain-energy density function [20,[40][41][42][43][44], where F represents the deformation gradient from the isotropic state, n is a unit vector, known as the director, for the orientation of the nematic field, and W (A) denotes the strainenergy density of the isotropic polymer network, depending only on the (local) elastic deformation tensor A. The tensors F and A are related through the identity where is the 'spontaneous' deformation tensor defining a change of frame of reference from the isotropic phase to a nematic phase. Here, we assume the existence of an isotropic phase that can be reached from a suitably imprinted defect field in the original fabrication process. The exact conditions under which such an intermediary configuration can be mapped to a stressfree reference configuration that is isotropic has been discussed in the general framework of anelasticity in [18]. In (3), a > 0 represents a temperature-dependent stretch parameter, ⊗ denotes the tensor product of two vectors, and I = diag(1, 1, 1) is the second-order identity tensor.

The Nematic Rod Model
Our goal is to devise a 1D model for a slender nematic solid which is sufficiently long and thin so that it can be approximated by a rod equation. For the derivation of the constitutive equations applicable to nematic rods, we rely on the following conditions assumed a priori: (R1) For the undeformed 3D structure, the dimension in one direction (the axial length) is much larger than in the orthogonal directions (the cross-sectional length scales); (R2) Any curvature-inducing deformation is of the same order of magnitude as the aspect ratio between the cross-sectional dimensions and the axial length.
These kinematic assumptions from the classic rod theory [32] are appropriate for filamentary nematic solids, and allow to reduce the dimensionality of the problem, which can be useful in analytical treatments and numerical simulations. Henceforth, we proceed in a similar manner as in [46] where morphoelastic rods were modeled. Apart from the different order of the multiplicative decomposition between elastic growth and LCE theories, we also bear in mind that, unlike in morphoelasticity, the natural deformation tensor for LCEs is always symmetric, i.e., G = G T , where T denotes the transpose operator, and this tensor changes if the nematic director rotates. In addition, as LCEs are generally incompressible, their deformations are isochoric, and in particular, det G = 1.

Model Reduction
We consider a circular cylinder occupying a domain = S × [0, L] ∈ R 3 in the reference configuration, with Cartesian coordinates (X, Y, Z) ∈ , where are cross-sections with centroids at Z ∈ [0, L]. In a Cartesian coordinate system, each material point (X, Y, Z) ∈ is identified with its position given by the vector Xe 1 + Y e 2 + Ze 3 , where (e 1 , e 2 , e 3 ) is the usual right-handed orthonormal basis. We assume that the typical length scale (Z) of every individual cross-section divided by the rod's length L is of order O(ε), where 0 < ε 1, and that its shape is a slowly varying function of Z, so that, locally, the structure is cylindrical. We define the scaled section The potential energy takes the form where W (A) is the elastic strain-energy function on which the nematic strain energy W (nc) (F, n) is based, and p (det A − 1) enforces the incompressibility condition det A = 1. We analyze deformations that map a right circular cylinder, with the longitudinal axis (centerline) in the Z-direction, to a filamentary body with the centerline r(Z), and set a local director basis (d 1 (Z), d 2 (Z), d 3 (Z)), such that r (Z) = (1 + εξ ) d 3 (Z), where denotes differentiation with respect to Z, and ξ is the small axial (longitudinal) strain. We define the Darboux curvature vector to describe the evolution of the director basis along the filament, satisfying where × is the usual vector product. A finite elastic deformation of the cylindrical body from the reference configuration B 0 to the current configuration B is then described in terms of its centerline and director basis by the following one-to-one, orientation-preserving mapping χ : → R 3 , where ρ i = ρ i (X, Y, Z), i = 1, 2, 3, are functions to be determined that describe the local deformation of each cross-section. Each of these functions satisfies the condition ρ i (0, 0, Z) = 0, so that the Z-axis is mapped to the centerline r(Z). Taking into account the small cross-sectional length scales, we denote ρ i = εα i , i = 1, 2, 3.
The corresponding deformation gradient tensor takes the form where the indices x and y denote the derivatives with respect to x and y, respectively. As the rod is slender, we assume a spontaneous deformation tensor of the form where g is an incremental remodeling tensor that may induce curvature, twist, or torsion. We recall that, for a straight rod in torsional deformation, each transverse cross-section is rotated by some angle while the longitudinal axis remains straight, so that the generators on the sides of the rod, which are initially parallel to the axis, become helical. By writing the Taylor expansions: and denoting we have the approximation The potential energy is then written as follows, where

Equilibrium Equations
On each cross-section S(Z), we define a sequence of minimization problems for the functionals The associated Euler-Lagrange equations take the form: with the natural boundary conditions: where n = [ n 1 , n 2 ] T is the outward unit normal vector to the boundary of S(Z). In these equations, at each order ε l , the choice of k ≤ l will depend on the solution at previous order. We note that the above system of equations is linear in α (k) 1 , α (k) 2 , α (k) 3 , and p (k) , for any given strain-energy density W (A) and spontaneous deformation tensor G. For l = 1, the equations are automatically satisfied, hence, α (1) 1 , α (1) 2 , α (1) 3 , and p (1) are only solved at order ε 2 . These variables take a form similar to the ones given in [46] and lead to the potential energy where

Examples of Rod Deformations
Since the strain-energy density of a nematic rod to the approximation considered here only depends on the linear behavior of the 3D strain-energy density [46], without loss of generality, we can use the form where H = A − I and μ is the shear modulus at small strain. We assume a nematic rod with the director taking the following general form in a reference system with cylindrical polar coordinates (R, , Z), n = sin cos e R + sin sin e + cos e Z , (25) where ∈ [−π/2, π/2] and ∈ [0, π] are given angles. The handedness of the fiber is given by the sign of with right-handedness obtained with angles strictly between 0 and π/2. Helical fibers are defined for = π/2 and purely radial fibers are obtained when = ±π/2. Equivalently, in the rectangular coordinate system, the nematic director is equal to n = sin cos( + )e 1 + sin sin( + )e 2 + cos e 3 .

Nematic Rod with Axial Director Field
First, we assume that the nematic field is aligned parallel to the longitudinal axis, i.e., = 0, so that n = e 3 . The corresponding spontaneous deformation tensor is equal to with diag(·, ·, ·) denoting the diagonal second order tensor. We set a 1/3 = 1 + εg(x, y). Hence, When A = G −1 F, we can follow the same steps as in [46] to reduce, after solving the Euler-Lagrange equations, the three-dimensional energy functional to a one-dimensional energy of the form where ζ = 1 + εξ and the stiffnesses are with E = 3μ representing the Young's modulus. Here, is the so-called warping function, and is a solution of where n is the unit vector normal to the boundary. The unstressed extension is with G(X, Y ) = εg(X/ε, Y /ε). Similarly, the unstressed curvatures are given by where The bending curvature and the torque of the deformed rod are, respectively (see [19, p. 103] for details), where denotes the first derivative with respect to Z. Taking the arc length S in a stress-free reference configuration and the twist as the angle between the normal vector to the centerline and d 1 , the quantity represents the rotation of local basis with respect to the Frenet frame as the arc length increases.

Photoresponses
Following [9], we assume that the anelastic photostrain p is due to the conversion of straight (trans) to bent (cis) forms of the dye molecules present in the LCE. In the simplest case, we can take the photostrain to be proportional to the cis fraction of chromophores n c , so that where λ p is the anelastic photostretch and A is the proportionality constant (taken to be 1 in all the particular examples and figures). When λ p = λ p (X, Y ) is given in the cross-section, the rod deformation can be found explicitly.
We can now couple the mechanical theory to the light intensity distribution by using a generalized Beer-Lambert law [27] for light absorbance. We define η to be the arc length along an optic path measured from the point η = 0 where the ray enters the material. We are interested in the light intensity I = I (η) in the material. Defining J = J (η) = I (η)/I (0), the generalized Beer-Lambert law reads where β is a material constant, n t is the proportion of trans molecules with the property n t + n c = 1, and where α = I (0)/I c is a measure of incident intensity, with I c the characteristic intensity [9]. Therefore, (41) reads which integrates to where W 0 denotes the Lambert function. We conclude that, along a single path, we have In general, at a given point (X, Y ) in the cross-section, multiple light rays contribute to the light intensity. Then, the total light intensity at that point is given by the sum of these different contributions from which the photostretch can be found.

Uniform Illumination: Axial Director, Square Cross-Section
As a simple first example, we consider a square cross-section and assume that the illumination enters a single side, as shown in Fig. 1(a). Moreover, on the illuminated surface, we assume that I 0 = I (0) is constant and that light enters the material along its normal direction, as illustrated in Fig. 1(c). In this case, we have η = X + 1 and each light ray only contributes to a single point in the cross-section. Using (35) with G(X, Y ) = An c , the only non-vanishing unstressed curvature reads Then, using (42) and (44), we find after simplification, the bending curvature (see Appendix A for details) where J 2 = J (2). Note that, after a suitable identification of variables, we recover the main result of [9] (viz., Equation 6, with d = 1/β and w = 2). Thus, the resulting deformation is pure bending, without torsion or twist. In particular, the principal bending stiffnesses K 1 and K 2 , and the torsional stiffness K 3 are the same as for the undeformed rod. Experimental observations of bending in LCE rods due to illumination are presented, for example, in [35,65,72]. Figure 2(a) shows the bending curvature κ scaled by β (note the penetration depth is 1/β, so κ/β is dimensionless) as a function of the thickness to penetration depth ratio 2β for various incident reduced light intensities α (see also Fig. 4 of [9]).

Uniform Illumination: Axial Director, Partially Coated Circular Cross-Section
We also take a circular cylindrical rod and assume again that on the illuminated surface (see Fig. 1(b)), I 0 = I (0) is constant, and that light enters the material along its normal direction, i.e., along the radial direction of the cross-section (see Fig. 1(c)). Then, if the cross-section is a unit circle coated in the opening with angle , we can write X = (η − 1) cos and Y = (η − 1) sin , with η ∈ [0, 2] and ∈ [π − /2, 2π − /2]. In this case, using (35) with G(X, Y ) = An c , the non-zero unstressed curvature is equal to (details in Appendix A) We now obtain the bending curvature where J 1 = J (1) and J 2 = J (2). Figure 2(b) shows the curvature κ scaled by 1/β as a function of 2β for various incident reduced light intensities α. Comparison with the mechanical behavior of a rod with a square cross-section suggests that there is very little difference between the two cases (see Figs. 2(a) and (b)).

Non-uniform Illumination: Optical Effects
The light paths in the above scenarios are simple, due to the assumption that the light penetrates the material only in the normal direction. To demonstrate that our modeling approach is also amenable to computing the curvature when the optical paths must be resolved, we consider horizontal light rays entering a material with circular cross-section. In this case, the light rays are refracted into the material at an angle φ 2 from the normal to the material according to Snell's law [6] where n 1 where n = n 1 /n 2 is the ratio of the refractive indices. For n close to 1, the light is refracted by a very small amount and will continue to travel nearly horizontally, with increasing deflection for decreasing n. Generally, we will restrict our attention to the range n < 1; in particular, if the external environment is air, then n 1 ≈ 1.00029, while a typical value for LCE is n 2 ≈ 1.64128 [7,33,51]. For given n, the light paths satisfy Noting that p 2 + q 2 = 1, t is the arc length along the path and may be identified with η in (41). We can now compute the curvature by numerically integrating (41) along the light paths, while simultaneously integrating for u 2 by converting equation (36) into an integral over s and t (details in Appendix B). In Fig. 3(a), we plot the scaled curvature κ/β against the thickness to penetration depth ratio 2β for different values of n 2 > 1 (with n 1 = 1), for fixed incident intensity α = 0.5. Also included in this plot for comparison are the cases of radial rays with = π , and horizontal rays on a square cross-section; these appear respectively as the black dashed and Fig. 3 Curvature generation for horizontal rays refracted into a circular cross-section. In (a), the scaled curvature is plotted against the thickness to penetration depth ratio 2β for varying refractive index of the nematic beam (with n 1 = 1): n 2 = 1.01 (red), 1.1 (orange), 1.6 (green), 2 (blue), and 4 (purple). Also included are the cases of radial rays (black dashed) and horizontal rays for a square cross-section (dotted black). The optical paths are illustrated for each case. In (b), the integrated curvature κ int is plotted against refractive index ratio n 2 /n 1 for the indicated values of incident intensity α. Light paths are shown (at evenly spaced points around the circle) for indicated values of n 2 /n 1 , with caustic curves highlighted in red dotted curves, almost indistinguishable. The first clear difference when comparing the cases of refracted and unrefracted light is that in the refracted cases the quantity κ/β diverges as β tends to zero. This is due to the existence of regions that are not reached by the refracted light (unlike the case of Fig. 1(a) or Fig. 1(b) with = π ). Therefore, even if J ≡ 1 across the entire path of the light rays, n c will remain 0 in the hidden regions so that a curvature still develops with β = 0, and thus κ/β → ∞ as β → 0.
It also appears from Fig. 3 that low refraction generates less curvature: the lowest curvatures occur for n 2 = 1.01 (red) and n 2 = 1.1 (orange), for which the light paths remain nearly horizontal. With increasing n 2 , there are two competing effects for curvature generation: (i) the light becomes more focused, with the light rays tending to the radial direction as n 2 → ∞, and (ii) the size of the hidden domain for which no light paths enter first increases and then decreases with increasing n 2 . The hidden domain is largest for n 2 around 2 (this is more clearly evident in the light paths shown in Fig. 3(b)). When β is small, this second effect is dominant, and the highest curvature is generated when the hidden region is largest. For larger β, the light intensity decays too quickly for the second effect to have much consequence; instead the curvature is higher when the light is more focused: thus the case of n 2 = 1.6 (blue) provides highest curvature for small β, while n 2 = 4 (purple) gives a higher curvature for larger β.
These observations demonstrate a strong dependence of curvature generation on both the optics and the depth of the material, and suggest the potential to balance the competing effects noted above in order to optimize the optical properties for curvature generation. To investigate further, we define an integrated scaled curvature which provides a measure of the ability of the material to generate curvature over a range of thickness to penetration depth ratios. In Fig. 3(b), we plot κ int as a function of 1/n = n 2 /n 1 for a range of values of incident intensity α. We verified numerically that the divergence as β → 0 is slower than 1/x for all values of n 2 /n 1 , so that κ int remains finite. These plots show an optimal value of n 2 ≈ 2.2. We see that the optimal value has almost no dependence on Fig. 4 In a cross-section S, the light intensity at a point p depends on the path-length η the intensity α, and moreover, the maximum value seems to saturate for α 0.5, suggesting that, once the light intensity is high enough, generating curvature is purely a function of the material geometry and optical properties. The dashed vertical line corresponds to n 2 = 1.64, showing that the value of refractive index typically cited for LCE is not far from the optimal value.

Non-uniform Illumination: Bending Towards the Light
Now, we consider a uniform inextensible transversely isotropic rod with a square crosssection that is clamped at one end. We choose the axes so that the sides of the cross-section at the clamped end are oriented along the rectangular X-and Z-directions, respectively, and, in the absence of illumination, the rod's longitudinal axis is parallel to the Y -direction. We are interested in finding its configuration when light rays are produced so that they are parallel to the X-direction, and the rod is illuminated from the negative X-direction. For simplicity, and to avoid any end effects, we also assume that the top end of the rod is opaque, which prevents the light from entering the material. This basic setup is inspired by the experiments in [35] showing how a clamped rod bends into the light. In this scenario, the light rays are not oriented with respect to the normal at each point on the rod surface, rather the light distribution will depend on the rod's inclination that itself depends on the rod's curvature, dictated by light. The problem can be solved by computing the contribution of different light rays to a given cross-section S, as depicted in Fig. 4, following Snell's law.
The light intensity at a point p ∈ S depends on the light decay along a path of length η. Hence, we can use equation (44) to obtain the light intensity J at p as follows, where η is now the length along the section S from the light entry point. Following Kirchhoff's assumptions on the geometry of a rod, the curvature and the incidence angle do not vary rapidly with arc length. Therefore the light paths impacting on a given cross-section will be approximately parallel. With this simplification, basic geometry gives η = η/ cos θ 2 , and therefore, (54) is equivalent to We observe that this relationship for the light distribution inside the material is the same as that given by (44), except that the parameter β is now replaced by β = β/ cos θ 2 (θ 1 ). The result on the curvature provided by (47) can then be readily adapted by replacing β with β in (47). Since a cross-section with θ 1 = π/2 does not receive any light, the curvature vanishes identically on any cross-section that reaches that particular value of the angle. We obtain the following unstressed curvature depending on the incidence angle, where From the curvature, we find the unstressed shape by first recalling that the tangent vector is equal to r (s) = − sin θ 1 e x + cos θ 1 e y and, second, by connecting the incidence angle to the curvature: θ 1 (s) = κ. We then derive the shape by integrating numerically the following system of ordinary differential equations, with initial conditions x(0) = y(0) = θ 1 (0) = 0. An illustration of this integration is shown in Fig. 5, where we see that there is a critical length L crit for the rod when θ 1 (L crit ) = π/2. For L < L crit , a rod of length L only bends partially towards the light, and for L > L crit , the rod has an horizontal segment with zero curvature of length L − L crit .

Nematic Rod with Helical Director Field
We further consider LCE rods with a helical nematic field (see [48,71] for relevant experimental tests, and also [17] for a theoretical discussion). In cylindrical polar coordinates, the nematic director takes the form n = sin e + cos e Z .
Equivalently, in Cartesian coordinates, the director is equal to n = − sin sin e 1 + sin cos e 2 + cos e 3 .
Then G takes the form Performing similar calculations to those for a cylindrical rod with axial nematic field, we obtain the unstressed extension and curvatures where, as before, A relates the concentration of cis molecules to their stretch (introduced in (40)).

Uniform Illumination: Helical Director, Partially Coated Circular Cross-Section
As an application of the curvatures generated by helical nematic directors, we consider a uniform rod with circular cross-section and a constant helical director field, with > 0, that only exists in a sector R ∈ [R 1 , R 0 ], ∈ [0, 2π]. The rod is partially coated so that light only enters in a section of the boundary circle as shown in Fig. 1(b) following the setup of Sect. 4.1.3. We also assume that the illumination is uniform. In this case, we can adapt the computation of Sect. 4.1.3 for the case of helical fibers to obtain the curvatures Fig. 6 Changes in intrinsic curvature u 3 and intrinsic twist u 3 , as a function of the inverse of light penetration depth β, for a nematic beam with circular cross-section and R 0 = 1, R 1 = 1/2 (ω = A = 1, = π/2, = π/4), under various incident reduced light intensities α. We observe the non-monotonic behavior of the twist as a function of α where we have taken R 0 = 1 and J 2 = J (2), J 12 = J (1 + R 1 ), J 11 = J (1 − R 1 ). An example of the behavior of this twist under different light penetration depths and reduced light intensities is shown in Fig. 6. We note that the intrinsic torsion is negative, creating a left-handed intrinsic helix even though the directors are right-handed. This is due to the fact that the molecules contract, hence applying an anticlockwise rotation of the cylinder that is compensated by a torsional response.
We can now use this exact form to compute the critical values at which a straight LCE rod would become unstable when illuminated. For example, we assume that such a straight rod is pulled with an end tension T and there is no twist before illumination. If the ends are prevented from turning, then, upon illumination, its intrinsic curvatures are given by (66)-(68), and will build elastic energy. Depending on the parameters of the system, this elastic energy can be sufficient to create an instability [23,39].
To compute the critical value of the parameters where such an instability occurs, we follow the general perturbation method given in [19, pp 136-137]. The Kirchhoff equations for an inextensible and unshearable rod, written in local director basis, are [46] n 3 u 2 − n 2 u 3 + n 1 = 0, (69a) Fig. 7 Critical tension at which a straight pulled (infinite) rod becomes unstable as a function of the illumination intensity α, for β = 1, . . . , 10. We assume a nematic rod with circular cross-section and R 0 = 1, R 1 = 1/2 (ω = A = 1, = π/2, = π/4, E = 1). For a given set of parameter and illumination, the straight rod is stable when the tension is above the critical value, given by equation (70) where {n 1 , n 2 , n 3 } and {u 1 , u 2 , u 3 } are the local components of the force and the curvature, respectively (see equation (7)).
One can readily verify that the straight configuration, with u 1 = u 2 = u 3 = n 1 = n 2 = 0, n 3 = T , satisfies equations (69a)-(69f) for all tension values T . We start in a regime of large tension, then decrease the tensile force until one reaches the critical value T crit . To find this limit, we linearize equations (69a)-(69f) around the straight solution and determine the largest tension that leads to the existence of a non-trivial solution. After simplification, we obtain where = K 1 /K 3 ∈ [2/3, 1]. We recognize a correction of the well-known twisting instability due to curvature and recover the twisting instability [21,22] as u 2 → 0. Since we know the values of the intrinsic curvatures from equation (66)- (68), we can compute the critical tension necessary to destabilize the system, as a function of the illumination for different penetration depths (see Fig. 7). We note that the main difference between an infinite and a finite length rod is a delay of the bifurcation, so that the critical tension for a finite length rod is slightly lower than the one given here for an infinite length. This difference vanishes as the rod increases in length.

Conclusions
We have presented in this paper a general method to construct an LCE rod model from a 3D material model. As long as the internal changes due to the excitation of the nematic directors remain moderate, so that the intrinsic curvatures are of the same order as the rod's slenderness ratio, the approximate strain-energy density of the LCE rods can be obtained through dimensional reduction.
We have demonstrated the utility of this approach through a series of examples. In the case of uniform illumination, we have obtained explicit formulas for the light-induced curvature as a function of the incident light intensity and penetration depth of the material. We have shown that our approach is also compatible with cases of non-uniform illumination, in which the light paths are non-trivial to resolve. The approach proposed here thus represents a major step in our ability to model these materials and use them in actual structures that respond to external stimuli such as heliotracking devices [25]. Indeed, the change of shape of a structure can be systematically recomputed by updating continuously the effect of the stimuli on the nematic directors. Rather than using a computationally expensive 3D finiteelement approach, the large deformations of the rod can be easily computed as an updated boundary-value problem at each time step as done, for instance, in the case of plant tropism [47].
At the same time, this work hints at a wealth of complexity that exists in the fully general case. Indeed, the curvature depends on the light paths, which depend on both the geometry of the material and the orientation in the light stimulus, which depends on the curvature. Even when solved in an iterative fashion, resolving the paths of refracted light across different cross-sections in a 3D geometry can pose a formidable challenge, particularly when shadowing effects are important. Nevertheless, in principle, our approach is amenable to any scenario, and in most cases analytic progress may be made under reasonable approximations such as uniform rays impacting on a given cross-section.
More generally, LCE materials open a new area of research where mechanics can be coupled with other physical effects affecting in real time the internal microstructure of the material. While great progress has been made in understanding the intrinsic geometry of such materials, much work remains to be done to fully couple geometry to mechanics and other external fields. We are now in a position to use the full power of nonlinear anelasticity to explore the behavior and shape these materials can acquire.

Appendix A: Curvature Calculations for Square or Circular Cross-Section
We provide below some detailed calculations for the curvature in the case of a square or circular cross-section. When the cross-section is a square, we first evaluate the integral in (46), Then, using this expression, we obtain (47) as follows, When the cross-section is circular, the integrals in (48) are evaluated as Then, (49) is calculated as follows,

Appendix B: Curvature Calculation with Refracted Light in Cross-Section
In this appendix, we outline the calculation of curvature for refracted light paths within a cross-section. According to Snell's law, the light will enter the material at an angle φ 1 from the normal to the boundary satisfying n 1 sin φ 1 = n 2 sin φ 2 , where φ 2 is the incident angle between the incoming light and the normal vector, and n 1 , n 2 are the refractive indices of the external environment and LCE material, respectively. In the case of a circular cross-section (X(s), Y (s)) = (cos s, sin s), and light rays traveling in the negative horizontal X-direction, the light will intersect the boundary for s ∈ [−π/2, π/2], and simple geometry gives the connection cos s = cos φ 1 .
Note that the plots in the main text show the case of rays traveling in the positive X-direction, but the curvature developed is equivalent up to a sign, and the calculation is somewhat cleaner with the light impacting on the right side of the circle, so that is what we present here. The light rays follow the path X(s, t) Y (s, t) = cos s + pt sin s + qt , with the unit vector (p, q) = − cos φ 2 e r + sin φ 2 e θ , according to Snell's law, where e r = (cos s, sin s) and e θ = (− sin s, cos s) are the unit radial and circumferential vectors, respectively. From Snell's law we compute also the relation cos φ 2 = n 2 2 − n 2 1 sin 2 s s .
Combining the above, and denoting n = n 1 /n 2 , we obtain p q = − 1 − n 2 sin 2 s cos s − n sin 2 s − 1 − n 2 sin 2 s sin s + n sin s cos s .
To compute the intrinsic curvature u 2 , we take equation (35) and the connection G(X, Y ) = An c = AαJ / (1 + αJ ), and convert the integral over the cross-section into an integral following the light paths. This gives αJ (s, t) 1 + αJ (s, t) ∂(X, Y ) ∂(s, t) dtds.
In the above expression, the ending value of t is defined by the point when each light ray reaches the other boundary of the cross-section. This point can be found by solving x 2 + y 2 = 1 for s, with X(s, t), Y (s, t) as defined above, which leads to T (s) = √ 2 2 − n 2 + n 2 cos(2s).