Large deformation analysis of spontaneous twist and contraction in nematic elastomer fibers with helical director

John Harrison invented the bimetallic thermostat in 1759. Ever since, scientists and engineers have been deploying spatial patterns of spontaneous deformation to induce complex and dramatic actuation in solid materials. The thermal strains in metals and shape memory alloys are limited to a few percent, but, in recent years, the soft-matter community has demonstrated several systems in which patterns of geometrically large strains can be programmed into soft solids. Prominent examples include patterns of swelling in gels,^1–4 patterns of contraction in liquid crystal elastomers (LCEs),^5–8 and patterns of inflation in “baromorphs.”^9,10 These large and exquisitely programmable shape changes can appear to bring the matter to life,^4,11,12 and this is no coincidence, as they strongly resemble the patterns of muscular contraction that drive biological locomotion and the patterns of growth that underpin biological development.¹³ 

This special issue focuses on the programming of such spontaneous shape changes in liquid crystal elastomers (LCEs). These are rubbery networks of rod-like mesogens, which spontaneously align along a director, $n$⁠, to form a nematic phase.¹⁴ On heating or illumination, the nematic order can be disrupted [reflecting the nematic to isotropic transition in conventional liquid crystals, Fig. 3(a)], and, in LCEs, this transition is accompanied by a dramatic and reversible contraction by a factor of $λs∼0.5$ parallel to the alignment director $n$⁠.^14,15 LCEs are thus promising artificial muscles and soft actuators.^16,17

Shape programming in LCEs is typically achieved by fabricating an elastomer in which the director $n$ is spatially varying, generating a corresponding pattern of contraction on heating. Director programming can be implemented by using surface-anchoring to pattern the director field in a nematic liquid sheet and then cross-linking to form an elastomer.^5,6 Alternatively, one may use an aligning stress field to orient the director during cross-linking. This latter strategy was used to create the original globally aligned monodomain LCEs,¹⁵ and, much more recently, has been deployed to generate patterns of alignment during extrusion based 3D printing^18–20 and during the direct shape programming of dual-network LCE sheets.⁸ 

The significant majority of work on LCE shape programming has focused on 2D sheets that bend and morph into curved surfaces on heating.^6,7,21,22 However, recently, Nocentini et al. demonstrated an LCE fiber that was twisted and stretched during cross-linking to imprint a helical director field.²³ On heating, the resulting helical contraction caused the fiber to spontaneously twist and contract. The design^24,25 and mechanics²⁶ of other torsional artificial muscles has been the focus of much attention recently, as they offer a miniaturizable version of a conventional rotary engine. However, these previous torsional muscles are fabricated by twisting multiple component fibers together^27–31 (for example, by twisting passive in-extensible fibers around an inflating core^28,29) and their action relies on slip within the resultant fiber bundle. In contrast, the LCE torsional muscle is a monolithic cross-linked solid and must be understood within the framework of misfit elasticity.

In this paper, we seek to construct such an elastic theory to predict and explain the spontaneous twist/stretch actuation of Nocentini et al.’s LCE fibers. At first sight, the natural starting point is Timishenko’s paradigmatic calculation of the curvature of a bimetallic strip.³² Indeed, analogs of this small-strain analysis are frequently deployed to model LCE bilayers, and other spontaneously bending sheets and strips. Moreover, very recently, several authors have developed corresponding theories of spontaneous bend and twist in elastic rods,^33–37 and these certainly offer considerable insight into Nocentini et al.’s LCE fibers. However, all such theories are only valid for small spontaneous strains (or, more precisely, small incompatibilities of spontaneous strains), leading to the complete decoupling of stretch, twist, and bend in the resultant elastic energies and making such theories formally inapplicable to the large strains generated in LCEs.

In contrast, in the field of rubber elasticity, there is a classic result, discovered by Poynting³⁸ in 1913, that a rubber fiber that is twisted substantially will also stretch in response. This large strain effect cannot be captured by the small-strain approaches but is captured by a simple and exact large-deformation solution for a twisted and stretched rubber cylinder.^39,40 Here, we derive an analog of this exact large deformation solution for a rubber cylinder subject to an axisymmetric pattern of spontaneous distortion. We find that these spontaneous distortions introduce simple modifications to the energy so that it is minimized by an overall spontaneous twist and contraction. We then compute these spontaneous twists and contractions for various LCE fibers with different director fields, highlighting how a helical field with both azimuthal and longitudinal components is required to produce a spontaneous twist. We validate our results with full 3D finite element simulations. Finally, we estimate the spontaneous twist and stretch expected in fibers created by applying twist and stretch during cross-linking, as reported by Nocentini et al.²³ Counterintuitively, we find that the maximum spontaneous twist is achieved by an optimal finite twist during cross-linking, with both too little and too much twist yielding lower performance. In conclusion, we discuss how these twisting fibers could be used to create coiling artificial muscles, with greatly amplified stroke compared to the intrinsic actuation of the LCE itself.

We start by recalling the classic large deformation solution for a twisted and stretched passive rubber cylinder.^39,40 More precisely, we consider a long cylindrical rubber fiber with undeformed radius $R0$ and length $L≫R0$⁠, which is subject to an angular twist $Δθ$ between the two ends and an overall extension by a factor of $λ$⁠. Working in cylindrical coordinates, if this deformation maps the material point initially at $R=(R,Θ,Z)$ to the point $r=(r,θ,z)$⁠, then the resultant deformation gradient is simply

In general, the hyper-elastic energy of a deformed solid may be written as

where $V$ is the volume in the reference configuration and $W$ is an energy density that depends on the deformation gradient. For an incompressible rubber, the simplest (neo-Hookean) energy density is given by

where $μ$ is the shear modulus and the second term arises to impose volume conservation, $det(F)=1$⁠, with $p$ being a spatially dependent Lagrange multiplier describing the pressure field in the material.

Our challenge is to minimize this energy for the fiber subject to an overall imposed twist $τ=Δθ/L$ and longitudinal stretch $λ$⁠. Since we are looking for states of uniform stretch and twist along the length of the cylinder, the outer surface must have the form

These forms clarify that $τ=dθdZ$ corresponds to the angular twist per unit reference state length.

Minimizing the elastic energy with respect to variations in $r(R)$ leads to the traditional bulk equations of mechanical equilibrium

and free boundary condition

where the first Piola–Kirchhoff stress tensor is given by

Finally, minimizing with respect to variations in $p$ returns the expected bulk condition of incompressibility

In the case of the twisted stretched fiber, the simplest possible fields are those with a homogeneous twist and (isochoric) stretch

which do indeed solve all the relevant bulk and boundary conditions. The associated deformation gradient is simply

and, upon substituting this back into the elastic energy, we obtain the total energy of a twisted and stretched fiber as

We note that the first two terms in the parenthesis are simply the familiar uni-axial stretching energy for a neo-Hookean rubber, while the final term determines the energy cost of twisting the fiber. The energy is trivially minimized by $τ=0$ for all values of $λ$⁠, but the minimum stretch is given by

with $λm=1$ only when $τ=0$⁠, as illustrated in Fig. 1(b). This asymmetric coupling between the two quantities is a direct manifestation of the Poynting effect⁴¹ in nonlinear elasticity. To clarify its origin, we note that the deformation gradient can be achieved as a pure twist followed by a pure stretch, as illustrated in Fig. 1(a). Mathematically, this corresponds to decomposing the deformation gradient as

where

and $δ$ is the identity matrix. From the form of $T$⁠, we see that a smaller radius implies a smaller deformation for the same twist. Indeed, the twist energy $R0τ2/λ$ can be written as $Rf2τ2$⁠, where $Rf=R0/λ$ is the final radius of the fiber. These considerations suggest that twist energy is partially relieved by stretch as this reduces the radius via Poisson effects. The resultant twist–stretch coupling is a paradigmatic example of the inherent geometric nonlinearity of large strains.

In light of what we learned from the simple twisting case, we now turn our attention to a cylindrical rubber fiber that undergoes a heterogeneous spontaneous distortion, such that, locally, the energy minimizing deformation is given by $F=G(R,Θ,Z)$⁠. Since our ultimate motivation is to understand the spontaneous twist and stretch of nematic LCE fibers, we restrict consideration to $G$ that is axisymmetric, isochoric, and independent of $Z$⁠. However, the pattern of spontaneous distortion is allowed to be incompatible so that the cylinder cannot attain $F=G$ throughout but will instead relax to an internally stressed state that minimizes the total elastic energy.

If the actual local deformation from the original state prior to spontaneous distortion is $F$⁠, then the elastic deformation from the local relaxed state is simply $F⋅G−1$⁠, where the second term reverses the effect of spontaneous deformation, and the first applies the actual deformation. The new elastic energy of the fiber, after spontaneous distortions, is thus

This “multiplicative decomposition” form⁴² was first introduced for elasto-plastic deformations⁴³ and now pervades and unifies the study of solids with spontaneous deformations, including growing tissues,^44,45 swelling gels,^46,47 thermal expansion,^34,48 and, as we shall clarify later, nematic elastomers heated to the isotropic state. Importantly, the spontaneous deformation only affects the energy via the combination $g=G−1⋅G−T$ (corresponding to the Finger tensor of $G$⁠) which, in our fibers, takes the symmetric and axisymmetric forms

The isochoric condition on all rubbers, including LCEs, requires $det(g)=det(G)=1$⁠.

Minimizing (15) with respect to variations in $r(R)$ and $p(R)$ leads to the same bulk and boundary equations as before [(5), (6), and (8)], but with the new Piola–Kirchhoff tensor given by

To solve for the actual deformation field, we first observe that, given axisymmetry, $Z$ independence, and incompressibility, the deformation must take the simple form

where $fθ$⁠, $fz$⁠, and $fp$ are as-yet unknown functions of $R$⁠, which modify the solution from the original one for a passive rubber cylinder. The $fθ$ term introduces rotation as a function of radius and allows radii in the reference configuration to become curves in the final configuration. Similarly, the $fz$ term encodes warping of each cross section into an identical surface of revolution, while $fp$ allows the pressure to vary with radius.

At this stage, one could methodically substitute these fields into the bulk and boundary equations and then solve for $fθ$⁠, $fz$⁠, and $fp$⁠. However, one can substantially simplify this process by noting that the resultant deformation field can now be decomposed as

where $λ$ and $T$ are again pure stretch and twist [Eq. (14)], while

describes the deformation of the cross section in the absence of a twist or stretch and is the only part that depends on $fθ$ and $fz$⁠. This decomposition is illustrated in Fig. 2.

Furthermore, since only derivatives of the two fields $fθ$ and $fz$ appear in $F$⁠, we can substitute this deformation into the energy and minimize directly with respect to variations in $fz′(R)$ and $fθ′(R)$⁠, leading to the conditions,

Evaluating these yields the simple uncoupled differential equations,

which can be directly integrated to obtain the fields,

These solutions then imply the full form of the deformation fields,

Finally, if we now substitute these back into the original bulk and boundary equations [(5) and (6)], we can confirm they are fully solved (e.g., in mathematica) provided the pressure is taken as

We now substitute the fields into the elastic energy [Eq. (15)] to evaluate the energy for a rubber rod with imposed twist $τ$⁠, stretch $λ$⁠, and a pattern of spontaneous deformation $g$⁠. To do this, it is convenient to use the decomposition in (19). We then see immediately that the effects of $g$ in the energy are entirely contained within the symmetric tensor,

which captures the residual part of $g$ after allowing the cross sections to relax via $Φ$⁠, but in the absence of twist or stretch. By inspection, we may then multiply out all the terms in the energy to find

with the coefficients given by

We note that these terms are related to the zeroth, first, and second moments of aspects of the spontaneous distortion, reminiscent of those derived via Gamma convergence for linear elastic rods.³⁴ 

As in the simple case explored in Sec. II, twist and stretch are coupled via the $λ−1$ in the twist energy. However, with the addition of the linear term in $τ$⁠, when $b≠0$⁠, the energy is minimized by a non-zero twist $τm$ and a finite stretch $λm$ given by

We can use these and $λ¯=λ/λm$ to rewrite energy with the same structure as in Eq. (11), but now including the source terms,

This is our main result and captures the emergent twist/stretch nature of fibers with a spontaneous deformation field.

Thus, our calculation reveals that an axisymmetric spontaneous deformation can induce 4 key effects in a fiber:

1. spontaneous length change,

2. spontaneous twist,

3. cross-section warping, via the function $fz$⁠, so initial disks become surfaces of revolution, and

4. cross-section twisting via $fθ$⁠, so initial radii become curves.

These effects stem from different aspects of the spontaneous deformation and may or may not co-occur in a given case.

In Sec. IV, we concentrate on spontaneous deformation fields arising from a LCE undergoing nematic–isotropic transition.

We want to apply our theoretical machinery to nematic elastomer cylinders, encoded with a spatially variable director pattern $n$⁠. On heating through the nematic–isotropic transition schematically shown in Fig. 3(a), the elastomer will spontaneously contract by a factor $λs$ along the director and, to preserve volume, elongate by $1/λs$ in the two perpendicular directions, corresponding to a spontaneous deformation $G=diag(λs−1/2,λs−1/2,λs)$ in a frame aligned with the director.¹⁴ We can thus write $g$ for such spontaneous actuation as

As discussed in the  Appendix, the resulting pre-strained neo-Hookean exactly reproduces the “trace-formula” energy commonly encountered in the LCE literature¹⁴ for an LCE cross-linked in the nematic state then heated to the isotropic state. This form neglects stress induced changes in the degree of nematic alignment, which is appropriate in the isotropic phase, where an LCE is simply a classical rubber provided one is not too close to the transition temperature. Nocentini et al.’s experimental LCE undergoes its transition in a 10 K window around 373 K and is heated to 393 K for actuation, allowing this approximation. We also note that the nematic actuation strains during the transition, $δL/L∼0.3$⁠, vastly exceed other thermal effects such as (linear) thermal expansion which has characteristic size $δL/L∼10−5K−1$⁠.

We again limit attention to axisymmetric director patterns. In cylindrical coordinates, the (unit) director field is then described by two angles $α(R)$ and $β(R)$ [Fig. 3(b)] so that

In what follows, we first consider the case of $α$ and $β$ independent of $R$ and show that spontaneous twist requires an oblique director field between azimuthal and longitudinal directions while warping of the cross sections requires a director that is oblique between the longitudinal and radial directions. Finally, we consider some more realistic $R$ dependent director fields for the fibers created by twisting and stretching during cross-linking²³ and show that spontaneous twisting is maximized by a finite optimum degree of twisting at genesis.

We start by considering director patterns in which the angles $α$ and $β$ are constants, independent of $R$⁠. In this case, the integrals in our solutions can be conducted analytically and yield the displacement fields,

From the linear dependence on $R$ in $z$⁠, we note that flat disk cross sections turn into cones, while $θ(R)$ shows that radii turn into spirals, which tend to conical-spirals at the center of the fiber. Though such singular structures may seem surprising at first glance, a more careful look reveals that integral curves of the director field $n$ are also conical spirals making a constant angle $β$ with the radial direction. Furthermore, in planar LCEs encoded with a constant angle +1 defects, the integral curves also form planar log spirals (the planar projection of a conical spiral) and it is well established that the resultant actuation transforms the sheet into twisted conical surfaces in which the radii transform into conical spirals.^22,49,50

In the simple case in which $β=π/2$⁠, the integral curves become simple helices, as shown in Fig. 3(b), and are helpful to understand how twist is developed during activation. In the most simple sketch of the mechanics, heat/illumination drives a contraction by $λs$ along the integral curve as well as an increase in its radius by a factor of $λs$⁠. The twisting thus occurs to reduce the length of the integral curve, with the integral curve playing the same role of a sub-fiber in twisted fiber bundle torsional muscles.^27–31 

Given constant $α$ and $β$⁠, the coefficients in our energy [Eq. (31)] are

where

These can be used together with Eqs. (32) and (30), to obtain the spontaneous twist,

The behavior of the twist as a function of $α$ and $β$ is shown in Figs. 4(a) and 4(b) for typical values during a nematic–isotropic (heating) transition (⁠$1>λs>0.5$⁠). Importantly, we see that the twist vanishes when $β=0$ or $α=0,π/2$⁠, corresponding to a purely longitudinal or azimuthal–radial director field where no shear $Θ$–$Z$ is present. On the other hand, $τm$ is maximized when $β=π/2$ (independent of $α$ and $λs$⁠), indicating a director with no radial component. The largest possible twist is then given by choosing $α=12cos−1⁡(λs3−1λs3+1)∼π/4+34(λs−1)+O((λs−1)2)$⁠, showing that maximum twist is achieved at an oblique angle biased towards the azimuthal direction for larger $λs$⁠, as shown in Fig. 4(b).

Similarly, the overall spontaneous stretch, $λm$⁠, is given by Eq. (30). Although the full expression is too complicated to reproduce here, in the simple and twist-maximizing case of $β=π/2$⁠, it reduces to

In Fig. 4(c), we can see how this stretch behaves between the two extremes of a longitudinally aligned (⁠$α=0$⁠) director field, yielding a simple contraction of $λm=λs$⁠, and of an azimuthal director field yielding $λm=(12(1+λs)/λs3)1/3$⁠.

To test our results, we use the open-source finite element software FEBio^51–53 to compute the spontaneous deformations of an LCE fiber encoded with an $R$-independent director field. We used FEBio’s standard prestrained neo-Hookean material on cylindrical fibers with an almost incompressible Poisson ratio of $ν=0.45$⁠. The isochoric prestrain parameter $1/λs$ was set to vary between values of $0.5$ and $2$⁠. After applying the pre-strain, the energy minimizing deformation was found using a static analysis.

In Fig. 5, we compare the predicted spontaneous twist and stretch of a fiber encoded with $β=π/2$ and $α=π/4$ with finite element simulations. The theory shows excellent agreement with the numerical simulations, accurately capturing the nonlinearities in both twist and stretch.

Fibers with $β=π/2$ have cross sections that remain flat during deformation, as coning is driven by the $gRZ$ component of the spontaneous deformation, which is only present if the director has an $RZ$ component. Therefore, to confirm our predictions about cross-section warping, we also computed the deformation of a fiber with $β=π/4$ and $α=π/4$⁠, as shown in Fig. 6(a). Again, comparing our theoretical results with finite element (FE) simulations, we obtain excellent agreement between the numerical and theoretical coning [Fig. 6(b)] and winding [Fig. 6(c)] of the cross section. The logarithmic nature of the spirals implies a theoretically infinite number of rotations at the center of the cross section, although the stress, strain, and energy are all finite. The divergent rotation stems from the line of director discontinuity (disinclination) along the central axis of the fiber. Accordingly, in a real fiber, the rotation would be cut off near the axis by a regularization of the director discontinuity within a defected core,⁵⁴ and in our finite elements, it is cut off by the element size near the axis. However, the director discontinuity line and the associated infinite rotation are an artifact of patterns of constant $α$ and $β$⁠, and, as discussed in Sec. IV B, are not expected in the experimentally generated fibers.

Finally, we consider the twisting LCE fibers reported by Nocentini et al.²³ These fibers were produced by pulling a filament out of a viscous LC monomer mixture while rotating the drawing end and simultaneously cross-linking with a UV light. The director alignment is imprinted through the strains induced during this drawing process, shown in Fig. 7(a). Given there is both twisting and stretching during cross-linking, and twisting strains are larger at larger radii, we expect this fabrication to produce an $R$-dependent director field with azimuthal and longitudinal components.

In reality, the imprinting of the director field is a complex visco-elastic process involving sticky polymers being cross-linked into a rubber. However, here, we take a simple approximation and assume the deformation is mainly elastic, and the director aligns with the direction of maximum strain. The elastic approximation is clearly appropriate once there is sufficient cross-linking but is probably also applicable to the initial visco-elastic drawing as the strain rate is rather high.

To find the imprinted director pattern, we model the fiber during drawing as an elastic cylinder that is stretched by a factor of $λ0$ and twisted by $Δθ0$ as cross-linking proceeds, resulting in a fiber of length $L$ and radius $R0$ and (final state) twist density $τ0=Δθ0/L$⁠. We use $(R~,Θ~,Z~)$ as the reference state coordinates for this problem, so that we may use $(R,Θ,Z)$ for the final state coordinates, which then become the reference state coordinates in our spontaneous deformation analysis. The elastic deformation follows the simple treatment in Sec. II, leading, in our coordinate system, to the deformation fields

To obtain the direction of maximum strain in the final (post cross-linking) configuration, we use the left-Cauchy deformation tensor

The largest eigenvalue of $b$ identifies the largest component of stretch while its corresponding eigenvector (which is a target state object) is its direction. Since this is the direction along which the director will orient, it can be used to express the values of the angles $α$ and $β$ in the fiber. We trivially obtain that $β=π/2$⁠, since the twisting during the manufacturing of the fiber induces no coupling of the $R$-$Z$ components. For $α$⁠, one obtains

where $τ¯0=Rτ0$⁠. We note that, when the fiber is drawn from a drop, $λ0≫1$⁠, which simplifies the angle to

This implies that the director points along the $Z$ direction in the center of the fiber and tilts in the $Θ$–$Z$ plane as one moves outward. This is a reflection of the fact that, during formation, $Θ$–$Z$ shears grow like $R$ as the filament is drawn and twisted, thus inducing a greater azimuthal component further from the center as shown in Fig. 7(b).

Finally, we can use the form of $α$ and $β$ to obtain $g$⁠. We then use Eqs. (29) and (30) to find the twist and stretch capability of a fiber given the twist imposed at its genesis. The results for values of $λs<1$ are shown in Fig. 8. Remarkably, the output twist does not monotonically grow as a function of $τ0$ but reaches a maximum and decays to zero thereafter. Recall that, in the $R$-independent field, we discussed how the twist is maximized when $α∼π/4$⁠. In this case, a small $τ0R0$ implies the director is on average mainly longitudinally aligned (⁠$<α>∼0$⁠), inducing mostly a contraction by a factor $λs$⁠. On the other hand, a large $τ0R0$ leads to a dominantly azimuthal director orientation on the cross section (⁠$⟨α⟩∼π/2$⁠), inducing mainly stretch by a factor of $(12(λs3+1)/λs3)1/3$⁠. The optimum twist output $τ$ is maximized in between the two, when the coupling between azimuthal and longitudinal components is the greatest.

The fibers produced by Nocentini et al.²³ were made from LCEs capable of a maximum spontaneous contraction of $λs∼0.71$ during heating. The twist imposed at genesis on the fibers was of about ten turns with their diameter and length varying between $50$ and $300μm$ and $1$ and $5cm$⁠, respectively. This suggests that their fibers were fabricated with a (dimensionless) genesis twist of around $0.16<R0τ0<1.8$ and the corresponding experimentally explored region is shaded in Fig. 8. The authors only reported the output twist and contraction for one fiber (of unknown length and radius), when activated in a $∼0.5cm$ long region via light. This fiber generated an overall contraction $∼0.84$⁠, suggesting that in the activated region $0.71<λm<0.84$⁠. On the other hand, activation induced a rotation of about $Δθ∼460o∼8.2rad$ corresponding to an output twist of $0.04<R0τ<0.5$⁠. Both these twist and stretch values fall in the shaded region of Fig. 8, consistent with our theoretical results. We highlight how, in general, these experimental fibers appear to have been generated with too little twist at genesis, yielding a sub-optimal output twist. This could perhaps be improved by increasing the number of turns during fabrication.

In conclusion, we have obtained an exact elastic solution for a twisted and stretched cylindrical neo-Hookean fiber subject to an axisymmetric isochoric spontaneous deformation field. The solution yields a full non-linear elastic energy for such a fiber, which is minimized by a spontaneous twist and stretch. The energy also highlights a large-deformation coupling between twist and stretch, as familiar from the classical Poynting effect. Finally, the elastic fields also capture the large deformations of the fiber’s cross section, which is predicted to warp into a surface of revolution and twist such that radii become curves during activation.

When applied to LCEs, our results show that a helical director field oblique in the longitudinal–azimuthal plane is required to induce a twist. The twist output depends on the spontaneous elongation $λs$ and is maximized when no radial director component is present (⁠$β=π/2$⁠) as well as when the azimuthal and longitudinal components are coupled through an angle $α∼π/4+34(λs−1)$⁠.

It is instructive to compare our theory with recent work on the spontaneous bending and twisting of rods via incompatible (aka misfit) spontaneous distortions.^33–37 These treatments go beyond ours in that the spontaneous strain is not assumed to be axisymmetric, and the rod is allowed to bend so that the centerline no longer remains straight. However, these treatments do assume small spontaneous strains and high aspect ratio rods, allowing a linear elastic treatment similar to the original Kirchhoff model. The small-strain thin-rod regime limits the theories to stretch free deformations and leads to simple bend twist energies of the form

where $κ$ and $τ$ are the bend (curvature vector) and the twist of the rod, while $κm$ and $τm$ are their minimizing values. These minimizing values were first estimated by linearizing the spontaneous deformation in a Taylor series about the rod’s central axis.³³ More recently, Gamma convergence^34,36,37 and 3D energy minimization³⁵ have been used to derive rigorous forms, yielding averages and moments of various terms of the spontaneous deformation over the rod’s cross section. Our large-deformation axisymmetric treatment reduces to the twisting portion of these rod-theories in the limit of small spontaneous distortions and little imposed stretch. Indeed, if we Taylor expand our expressions for $λm$ and $τm$ [Eq. (30)] in the limit of small spontaneous deformations $g=δ+ϵδg(R)$⁠, we find that

If we then also assume the applied twist and stretch are small, $λ−λm∼ϵ$ , $τ∼ϵ$⁠, we may expand our energy [Eq. (31)] to $ϵ2$ to obtain

The twisting portion of this energy agrees with that in Eq. (40), and the linearized form of $τm$ above agrees with that given by Kohn and O’Brien.³⁴ Interestingly, although such small-strain and high aspect ratio assumptions appear necessary to resolve bending, our treatment demonstrates they can be avoided entirely when only treating twisting and stretching. The resulting non-linear formulation generates the highly non-linear form $τm(λs)$ and $λm(λs)$⁠, as seen in Fig. 5, and is clearly required for accurate predictions in large strain systems such as LCEs.

Remarkably, even a small strain rod theory can describe large bend and twist displacements in a suitably long rod. This introduces fascinating and rich geometric coupling between twist and bend deformations.^26,55–57 For example, if one straightens a wound headphone wire, it becomes highly twisted. Similarly, if one twists a fiber and then bring the ends together, it will spontaneously untwist into a lower energy spring-like coiled state. In general, torsional bend and twist can be exchanged in a rod or fiber, without rotating the ends, provided the total number of turns is conserved. This twist–torsion coupling is an example of a geometric phase and is key in the winding of DNA molecules,^58–60 and the coiling of plant tendrils.⁶¹ 

Therefore, although our treatment includes no mention of bend, we can infer from this coupling effect that a spontaneously twisting LCE fiber will bend into a coiled spring-like configuration if it activates under boundary conditions that prevent the ends from twisting relative to each other. Given the large-strain nature of LCE actuation, a formal treatment of this torsional effect appears to be a formidable challenge. However, as shown in Fig. 9, we were easily able to observe it in finite elements simulations of a twisting LCE fiber actuated under a constant longitudinal force but with the constraint that the ends may not rotate. This coiling mechanics is commonly deployed in other twisting artificial muscles to generate linear actuation.^25,26 In contrast, LCEs are intrinsically contractile actuators, and simple linear contraction can be trivially achieved with a monodomain strip.^15–17 However, the coiling mechanism would allow an LCE actuator with much higher stroke and lower stiffness, forming switchable soft-springs that, like plant tendrils,⁶¹ could be used to gently anchor and position an object in 3D space.

A.G. acknowledges J.L. for support and the EPSRC for funding (Project No. 2108804). J.S.B. is supported by a UKRI Future Leaders Fellowship (No. MR/S017186/1).

The data that support the findings of this study are available within the article.

The free energy density of a nematic elastomer¹⁴ has its roots in statistical physics and contains two contributions: an elastomer energy from the polymer network and a nematic energy from the rods,

The rod energy depends on the scalar order parameter $Q$ of the nematic field (i.e., the degree of alignment) but not its direction (unit vector $n$⁠). Appropriate forms for $Wrod(Q)$ are provided by the Landau–De Gennes (phenomenological) or Maier Saupe (microscopic) theories of liquid nematics.⁵⁴ Either way, $Wrod(Q)$ has the characteristic size of $kBT$ per rod, and transitions from having a minimum at $Q=0$ (isotropic) to a finite $Q$ (nematic) below a critical temperature $T∗$⁠.

The polymer free energy is dominated by conformational entropy, like in the statistical theory of conventional rubber. However, in the presence of a nematic field, the polymer random walks are not isotropic but biased along the director $n$ by an amount $r(Q)$ determined by the degree of alignment, as encoded in a step-length tensor $ℓ∝δ+(r−1)nn$⁠. The resultant polymer energy is described by the “trace formula”^14,62

where $ns$ is the density of polymer strands, $F$ is the deformation from the cross-linking state to the final state, $ℓ0$ is the step-length tensor at cross-linking (which depends on the nematic variables at cross-linking, $Q0$$n0$⁠), and $ℓ$ is the step-length tensor in the final state (which depends on the final state nematic variables at $Q$⁠, $n$⁠). The full behavior of the nematic elastomer is now given by minimizing the sum of both energies over elastic deformations (⁠$F$⁠), final state order parameter (⁠$Q$⁠), and final state director (⁠$n$⁠). In general, this minimization gives a two way coupling between the nematic order and the LCE deformation. However, the characteristic size of $Wrod$ is $kBT$ per rod, while the characteristic size of $Wpol$ is $kBT$ per polymer strand. Since in an elastomer there are typically more than ten rods per strand, the nematic energy dominates the elastic one during minimization over $Q$⁠. Therefore, the degree of nematic alignment is essentially that which minimizes the nematic energy alone and is only modestly affected by the polymer network. For example, the shift of the nematic–isotropic transition caused by the presence of the network is typically only a few Kelvin⁶³ (compared to a transition temperature of 350 K), and other mechanically induced changes in $Q$ are similarly small.^14,64 Such effects are more pronounced in some modern LCE compositions, which have more cross-links per rod,^65,66 but, for simplicity, we focus on the traditional case here.

We may thus consider the simpler problem of minimizing the polymer energy, with the magnitude of nematic order, $Q$⁠, effectively fixed as a constraint by the rod’s energy. In general, when we minimize the polymer energy in the nematic state, we must still do it over final state director $n$⁠; indeed, the nematic director can be observed to rotate within an elastomer in response to stretch.^14,67 However, if we further assume (as here) that the final state is isotropic, $Q=0$⁠, then $ℓ∝δ$ must also be isotropic, and there is no final state director to minimize over. In this case, the elastomer energy is simply

which corresponds to the standard neo-Hookean energy with spontaneous distortion $G=ℓ01/2$ encoded by the nematic field in the cross-linking state and shear modulus $μ=nskBT$⁠.