# Instabilities in liquid crystal elastomers

L. Angela Mihai† Alain Goriely∗

†School of Mathematics, Cardiff University, Cardiff, UK
∗Mathematical Institute, University of Oxford, Oxford, UK

ORCID: L.A.M. (0000-0003-0863-3729); A.G. (0000-0002-6436-8483)

###### Abstract

Stability is an important and fruitful avenue of research for liquid crystal
elastomers. At constant temperature, upon stretching, the homogeneous state of
a nematic body becomes unstable, and alternating shear stripes develop at very
low stress. Moreover, these materials can experience classical mechanical
effects, such as necking, void nucleation and cavitation, and inflation
instability, which are inherited from their polymeric network. We investigate
the following two problems: First, how do instabilities in nematic bodies
change from those found in purely elastic solids? Second, how are these
phenomena modified if the material constants fluctuate? To answer these
questions, we present a systematic study of instabilities occurring in nematic
liquid crystal elastomers, and examine the contribution of the nematic
component and of fluctuating model parameters that follow probability laws.
This combined analysis may lead to more realistic estimations of subsequent
mechanical damage in nematic solid materials.

Impact statement: Due to their complex material responses in the presence of
external stimuli, liquid crystal elastomers have many potential applications
in science, manufacturing, and medical research. The modeling of these
materials requires a multiphysics approach, linking traditional continuum
mechanics with liquid crystal theory, and has led to the discovery of
intriguing mechanical effects. An important problem for both applications and
our fundamental understanding of nematic elastomers is their instability under
large strains, as this can be harnessed for actuation, sensing, or patterning.
The goal is then to identify parameter values at which a bifurcation emerges,
and how these values change with external stimuli, such as temperature or
loads. However, constitutive parameters of real manufactured materials have an
inherent variation that should also be taken into account. Hence, the need to
quantify uncertainties in physical responses, which can be done by combining
the classical field theories with stochastic methods that enable the
propagation of uncertainties from input data to output quantities of interest.
The present study demonstrates how to characterize instabilities found in
nematic liquid crystal elastomers with probabilistic material parameters at
the macroscopic scale, and paves the way for a systematic theoretical and
experimental study of these fascinating materials.

Key words: polymer, liquid crystals, modeling, uncertainty quantification,
instabilities.

##  1 Introduction

Liquid crystal elastomers (LCEs) are advanced multifunctional materials that
combine the flexibility of polymeric networks with the nematic structure of
liquid crystals [24, 32]. Due to their complex molecular architecture, they
are capable of exceptional responses, such as large spontaneous deformations
and phase transitions, which are reversible and repeatable under certain
external stimuli, namely: heat, light, solvents, electric or magnetic fields.
These properties render them as promising candidates for future ‘animate
materials’, and could be harnessed for a range of technological applications
including soft actuators and soft tissue engineering. Nevertheless, a better
understanding of these materials is needed before they can be exploited on an
industrial scale [26, 58, 59, 64, 67, 89, 107, 108, 111, 112, 118, 120, 121,
122].

Since the early discovery of liquid crystalline solids, probing their
intriguing material properties has been the focus of research laboratories
around the world, and the importance of such essential work is hard to
overstate. However, their accurate description can only be useful if fully
integrated in a multiphysics framework combining elasticity and liquid crystal
theories. Many nematic solids are synthesized as polydomains, where the liquid
crystal mesogens are separated into different domains, and in every domain,
they are aligned along a preferred direction, known as the local director [19,
20, 61, 93, 104, 109]. Depending on the fabrication process, polydomains may
have very different material properties and behaviors. Monodomains, where
mesogen molecules are uniformly aligned throughout the material, can be formed
from polydomains through mechanical stretching or by cooling an isotropic
material under an external stress field to reach the nematic phase.

An ideal continuum model for monodomains is provided by the neoclassical
strain-energy function proposed in [15, 114, 117]. This is a phenomenological
model based on the molecular network theory of rubber elasticity [106]. The
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 [115, 116]. Since elastic stresses dominate over
Frank elasticity induced by the distortion of mesogens alignment [8, 83, 84],
Frank effects are generally ignored [25, 36]. The neoclassical formulation has
been extended to polydomains by assuming that every single domain has the same
strain-energy density as a monodomain [13, 14]. These descriptions have been
generalized to include nematic strain-energy densities based on
phenomenological hyperelastic models (e.g., Mooney-Rivlin, Gent, Ogden) that
better capture the nonlinear elastic behavior at large strains [3, 2, 30]
(molecular interpretations of the Mooney-Rivlin and Gent constitutive models
for rubber are presented in [37, 54]). Further generalizations are proposed in
[4, 125].

Another important characteristic of most materials is that physical properties
are subject to random variations. Typically, average properties are used and
any variations around the mean are neglected when computing the response of a
material. However, as we show here, one can employ information theory [94] and
the maximum entropy principle [55, 56, 57] to include some stochastic
variations of the material parameters, then propagate their uncertainty to
output physical responses. Comprehensive reviews on the information-theoretic
approach in elasticity are presented in [50, 49, 97].

Our main focus is on large strain instabilities of LCE bodies acted upon by
external loads. In addition to the recurring phenomenon of _soft elasticity_ ,
where alternating shear stripes develop at very low stress if a nematic body
is stretched, we explore theoretically a set of classical instabilities
inherited from parent polymeric networks, namely: _necking_ under tensile
load, _cavitation_ of a nematic sphere where a void nucleates at its center
when uniform tensile traction is imposed, and _inflation instability_ of an
internally pressurized shell where the pressure increases, decreases, and then
increases again. The aim is to determine conditions for the onset of
instability and show how nematic materials perform compared to their purely
elastic analogue. Moreover, for these problems, the propagation of stochastic
variation from input material parameters to output mechanical behavior is
mathematically traceable, making our stochastic approach both mathematically
and mechanically transparent. The main effect of random variations is to
replace a well-defined bifurcation point by a probability that a material will
undergo an instability as a function of the bifurcation parameter. Such
fundamental problems are important in their own right and may stimulate
related mechanical testing of nematic materials.

##  2 General set-up

Following the classical work of Flory [35] on polymer elasticity, we use the
stress-free state of a virtual isotropic phase at high temperature as the
reference configuration [8, 18, 27, 28, 29, 30], rather than the nematic phase
in which the cross-linking might have been produced [4, 15, 119, 113, 114,
117, 125]. Within this theoretical framework, the material deformation due to
the interaction between external stimuli and mechanical loads can be expressed
as a composite deformation from a reference configuration to the current
configuration, via an elastic deformation followed by a natural (stress free)
shape change. The multiplicative decomposition of the associated gradient
tensor is similar in some respects to those found in the constitutive theories
of thermoelasticity, elastoplasticity, and growth [47, 63] (see also [46,
92]), but is different on one major aspect, namely: the stress-free
geometrical change is superposed on the elastic deformation, which is applied
directly to the reference state. This difference is important since, although
the elastic configuration obtained by this deformation may not be observed in
practice, it may still be possible for the nematic body to assume such a
configuration under suitable external stimuli. The resulting elastic stresses
can then be used to analyze the final deformation where the particular
geometry also plays a role [74, 75].

Figure 1: Multiplicative decomposition of the deformation gradient for a
nematic elastomer. Unlike traditional theories of anelasticity, the elastic
deformation is applied first to the reference configuration and then followed
by an anelastic strain deformation.

To describe an incompressible nematic material, we combine isotropic
hyperelastic and neoclassical strain-energy density functions as follows [73],

$\overline{W}^{(nc)}(\textbf{F},\textbf{n})=W\left(\textbf{F}\textbf{G}_{0}^{-1}\right)+W^{(nc)}(\textbf{F},\textbf{n}),$
(1)

where, on the right-hand side, the first term is the energy of the ‘parent’
elastic matrix, and the second term is the neoclassical-type function.
Specifically: n is a unit vector for the localized direction of uniaxial
nematic alignment in the present configuration;
$\textbf{F}=\textbf{G}\textbf{A}$ is the deformation gradient tensor with
respect to the reference isotropic state (see Figure 1 and also Figure 1 of
[30]), with
$\textbf{G}=a^{1/3}\textbf{n}\otimes\textbf{n}+a^{-1/6}\left(\textbf{I}-\textbf{n}\otimes\textbf{n}\right)$
the ‘spontaneous’ (or ‘natural’) deformation tensor and A the (local) elastic
deformation tensor;
$\textbf{G}_{0}=a^{1/3}\textbf{n}_{0}\otimes\textbf{n}_{0}+a^{-1/6}\left(\textbf{I}-\textbf{n}_{0}\otimes\textbf{n}_{0}\right)$
is the spontaneous deformation tensor with $\textbf{n}_{0}$ the director
orientation at cross-linking (which may be spatially varying); and $a>0$ is a
temperature-dependent shape parameter, which we assume to be spatially
independent (i.e., no differential swelling). We denote by $\otimes$ the usual
tensor product of two vectors, and by $\textbf{I}=\text{diag}(1,1,1)$ the
identity second order tensor.

Recognizing that some uncertainties may arise in the mechanical responses of
liquid crystal elastomers, and inspired by recent developments in stochastic
finite elasticity [34, 68, 69, 70, 71, 78, 79, 80, 81, 98, 99, 100, 101, 102],
we assume that the model parameters are defined as random variables drawn from
given probability distributions. In practice, material parameters take on
different values, corresponding to possible outcomes of experimental tests.
The maximum entropy principle then allows us to explicitly construct prior
probability laws for the model parameters, given the available information.
Explicit derivations of probability distributions for the elastic parameters
of stochastic homogeneous isotropic hyperelastic models calibrated to
experimental data for rubber-like material and soft tissues were presented in
[34, 78, 100]. Intuitively, such a stochastic body can be regarded as an
ensemble (or population) of bodies that are equal in size and have the same
geometry, and each body in the ensemble is made from a single homogeneous
material with the parameters not known with certainty, but distributed
according to probability density functions that are calibrated to macroscopic
experimental measurements. These models reduce to the usual deterministic ones
when the parameters are single-valued constants.

##  3 Soft elasticity and stress plateaus

Many macroscopic deformations of nematic liquid crystal elastomers induce a
re-orientation of the director with a general tendency for the director to
become parallel to the direction of the largest principal stretch. This re-
orientation is typically uniform across the material. However, non-uniform
behaviors are also possible. In particular, under appropriate uniaxial tension
or biaxial stretch, bifurcation to a pattern of stripe domains is generated,
where adjacent stripes deform by the same shear but in opposite directions.
Early experimental investigations of this phenomenon, known as _soft
elasticity_ [87, 110, 113], were reported in [33, 60, 103, 126]. Its
theoretical explanation is that, for these materials, the energy is minimized
by passing through a state exhibiting a microstructure of many homogeneously
deformed parts [23, 28, 30, 38, 39, 40, 60]. A natural question is then: _How
does soft elasticity depend on the material parameters?_

Figure 2: (a)-(b) The strain-energy function $w(\lambda,\varepsilon,\theta)$,
for $(\varepsilon,\theta)=(0,0)$,
$(\varepsilon,\theta)=(\varepsilon_{0},\theta_{0})$, with $\varepsilon_{0}$
and $\theta_{0}$ given by (7), and $(\varepsilon,\theta)=(0,\pi/2)$, when
$a=2$ and (a) $\mu^{(1)}=4.05$ or (b) $\mu^{(1)}=0$, while $\mu^{(2)}=4.05$.
The two vertical lines correspond to the lower and upper bounds on $\lambda$,
given by equations (5) and (6), respectively. Between these bounds, the second
solution, with $(\varepsilon,\theta)=(\varepsilon_{0},\theta_{0})$, minimizes
the energy. The lower and upper bounds on $\lambda$, given by (c) (5)-(6), and
(d) (8)-(9), respectively, as functions of the parameter ratio
$\eta=\mu^{(1)}/\mu^{(2)}$ when $a=2$. For $\eta\to\infty$, the model
approaches a purely elastic form where the homogeneous deformation is always
stable. For $\eta=0$, the bounds are the same, and correspond to the
neoclassical form.

Figure 3: The stochastic shear parameter $\varepsilon_{0}$ and director angle
$\theta_{0}$, given by (7), when: (a)-(b) $a=2$, while the shear modulus $\mu$
follows a Gamma distribution with shape and scale parameters $\rho_{1}=405$
and $\rho_{2}=0.01$, respectively, and $R^{(1)}=\mu^{(1)}/\mu$ is drawn from a
Beta distribution with hyperparameters $\xi_{1}=100$ and $\xi_{2}=100$;
(c)-(d) $\mu=4.05$ and $R^{(1)}=\mu^{(1)}/\mu=0.5$, while the shape parameter
$a$ follows a Gamma distribution with $\rho_{1}=200$, $\rho_{2}=0.01$. In
(a)-(b), the dashed black lines correspond to the deterministic solutions
based only on the mean values $\underline{\mu}=\rho_{1}\rho_{2}=4.05$ and
$\underline{R}_{1}=\xi_{1}/(\xi_{1}+\xi_{2})=0.5$, whereas the red versions
show the arithmetic mean value solutions. In (c)-(d), the dashed black lines
correspond to the deterministic solutions based only on the mean value
$\underline{a}=\rho_{1}\rho_{2}=2$, whereas the red versions show the
arithmetic mean value solutions. In this case, there is a significant
difference between the mean value solutions and the deterministic solutions.

For simplicity, we select an incompressible neo-Hookean-type strain-energy
function [105],

$W(\textbf{A})=\frac{\mu^{(1)}}{2}\left[\text{tr}\left(\textbf{A}\textbf{A}^{\text{T}}\right)-3\right],$
(2)

where the superscript “T” represents the transpose operator, “tr” denotes the
trace operator, and $\mu^{(1)}\geq 0$ is constant, together with the
neoclassical strain-energy function [18, 28, 30],

$W^{(nc)}(\textbf{F},\textbf{n})=\frac{\mu^{(2)}}{2}\left\\{a^{1/3}\left[\text{tr}\left(\textbf{F}\textbf{F}^{\text{T}}\right)-\left(1-a^{-1}\right)\textbf{n}\cdot\textbf{F}\textbf{F}^{\text{T}}\textbf{n}\right]-3\right\\},$
(3)

with $\mu^{(2)}\geq 0$ constant. For the composite model function defined by
(1), the shear modulus at infinitesimal strain is $\mu=\mu^{(1)}+\mu^{(2)}>0$
[72]. However, our results can be easily extended to other choice of strain-
energy density functions [73].

We analyze shear striping under biaxial stretch, and assume that the nematic
director can only rotate in the biaxial plane. To achieve this, we set
$\textbf{n}_{0}=[0,0,1]^{\text{T}}$ and
$\textbf{n}=\left[0,\sin\theta,\cos\theta\right]^{\text{T}}$, where
$\theta\in[0,\pi/2]$, in a Cartesian system of reference, and examine small
shear perturbations of biaxial extensions, with gradient tensor [73]

$\textbf{F}=\left[\begin{array}[]{ccc}a^{-1/6}&0&0\\\ 0&\lambda&\varepsilon\\\
0&0&a^{1/6}\lambda^{-1}\end{array}\right],$ (4)

where $a>1$ is the nematic shape parameter, $\lambda>0$ is the stretch ratio,
and $\varepsilon>0$ is the small perturbation. Denoting
$w(\lambda,\varepsilon,\theta)=\overline{W}^{(nc)}(\textbf{F},\textbf{n})$,
for sufficiently small values of $\varepsilon$ and $\theta$, we find that,
when $\mu^{(1)}>0$ and $\mu^{(2)}=0$ (purely elastic case), the equilibrium
state with $\varepsilon=0$ and $\theta=0$ is stable. If $\mu^{(2)}>0$, then
this equilibrium state is unstable for

$a^{-1/6}\leq
a^{1/12}\left(\frac{\eta+1}{\eta+a}\right)^{1/4}<\lambda<a^{1/12},$ (5)

where the _elasto-nematic ratio_ $\eta=\mu^{(1)}/\mu^{(2)}$ gives the relative
magnitude of the elastic and neoclassical contributions. There is also an
equilibrium state with $\varepsilon=0$ and $\theta=\pi/2$, where the nematic
director is fully rotated so that it aligns uniformly with the direction of
macroscopic extension. By symmetry arguments, this state is unstable for

$a^{1/12}<\lambda<a^{1/12}\left(\frac{\eta+a}{\eta+1}\right)^{1/4}\leq
a^{1/3}.$ (6)

When $\lambda$ satisfies (5) or (6), we obtain

$\varepsilon_{0}=\pm\frac{\lambda(a-1)\sin(2\theta_{0})}{2\sqrt{\left(\eta+1\right)\left(\eta+a\right)}},\qquad\theta_{0}=\pm\arccos\sqrt{\frac{a^{1/6}\sqrt{\left(\eta+1\right)\left(\eta+a\right)}}{\lambda^{2}(a-1)}-\frac{\eta+1}{a-1}}.$
(7)

For the resulting strip pattern, the gradient tensors of alternating shear
deformations in two adjacent stripe domains are $\textbf{F}_{\pm}$ with
$\varepsilon=\pm\varepsilon_{0}$, respectively. The two deformations are
geometrically compatible in the sense that there exist two non-zero vectors q
and p, such that the Hadamard jump condition
$\textbf{F}_{+}-\textbf{F}_{-}=\textbf{q}\otimes\textbf{p}$ is satisfied,
where p is the normal vector to the interface between the two phases
corresponding to the deformation gradients $\textbf{F}_{+}$ and
$\textbf{F}_{-}$ [11, 12, 76, 77]. In other words, $\textbf{F}_{+}$ and
$\textbf{F}_{-}$ are rank-one connected, i.e.,
$\mathrm{rank}\left(\textbf{F}_{+}-\textbf{F}_{-}\right)=1$. If $\eta=0$, then
the above equilibrium states are unstable for
$\lambda\in\left(a^{-1/6},a^{1/12}\right)$ and
$\lambda\in\left(a^{1/12},a^{1/3}\right)$, respectively. Thus, soft elasticity
is always presented by the purely neoclassical model [23, 30]. When
$\eta\to\infty$, there is no shear striping since the material is practically
elastic.

The strain-energy function $w(\lambda,\varepsilon,\theta)$ is illustrated in
Figure 2(a)-(b). Note that, for $\lambda$ with values between the lower and
upper bounds given by (5) and (6), respectively, the minimum energy is
attained for $(\varepsilon,\theta)=(\varepsilon_{0},\theta_{0})$. Assuming
that loading is applied in the second direction, the first Piola-Kirchhoff
axial stress in this direction is equal to
$P^{(nc)}_{2}=\mathrm{d}w(\lambda,\varepsilon,\theta)/\mathrm{d}\lambda$.
Figure 2(b) then suggests that, when $\mu^{(1)}=0$, i.e., for the purely
neoclassical form, the director rotates and alternating shear stripes develop
for $\lambda\in\left(a^{-1/6},a^{1/3}\right)$, at zero load, since the slope
of the curve is equal to zero within this interval. In contrast, if
$\mu^{(1)}>0$, then from Figure 2(a), we infer that the applied load increases
with deformation and is almost constant but non-zero while the director
rotates.

Due to the geometric compatibility, and since the intervals for stretch ratios
$\lambda$ where shear striping occurs are at a maximum length when
$\textbf{n}_{0}=[0,0,1]^{\text{T}}$, the bounds (5)-(6) give also the maximum
interval for shear striping when $\textbf{n}_{0}$ is not uniformly aligned.
The minimum length of those intervals is attained for monodomains with
$\textbf{n}_{0}=[0,1,0]^{\text{T}}$. Experimental results for monodomains
where the tensile load forms different angles with the initial nematic
director are reported in [82, 86]. If $\textbf{G}_{0}=\textbf{I}$ (see [13]),
then the solution with $\varepsilon=0$ and $\theta=0$ is unstable for

$a^{-1/6}\leq
a^{1/12}\left(\frac{\eta+a^{-2/3}}{\eta+a^{1/3}}\right)^{1/4}<\lambda<a^{1/12},$
(8)

and that with $\varepsilon=0$ and $\theta=\pi/2$ is unstable for

$a^{1/12}<\lambda<a^{1/12}\left(\frac{\eta+a^{1/3}}{\eta+a^{-2/3}}\right)^{1/4}\leq
a^{1/3}.$ (9)

For example, when $a=2$, the bounds given by (5)-(6), and by (8)-(9),
respectively, are plotted as functions of the parameter ratio $\eta$ in Figure
2(c)-(d).

In addition, when the model parameters are random variables, characterized by
probability distributions, we are interested in the probability that shear
striping occurs under a given stretch. Here, we show the effects of
fluctuations in the shear modulus $\mu$ and the shape parameter $a$
separately, i.e., we first set $\mu$ as a random variable while $a$ is a
single-valued constant, then keep $\mu$ constant and let $a$ fluctuate. In
Figure 3(a)-(b), the stochastic shear parameter $\varepsilon_{0}$ and director
angle $\theta_{0}$, given by (7), are represented when the stochastic material
parameters were chosen such that their mean values are the same as the
deterministic values in Figure 2(a). For both stochastic versions, to increase
the probability of homogeneous deformation, one must increase the value of
$\eta$, whereas shear striping is certain only if the model reduces to the
neoclassical one (i.e., when $\eta=0$). However, while $\eta>0$, the inherent
variability in the probabilistic system means that there will always be
competition between the homogeneous and inhomogeneous deformations.

##  4 Cavitation

Cavitation in solids is the formation of a void within a solid under tensile
loads, by analogy to the similar phenomenon observed in fluids. For rubber-
like materials, following early studies of damage under loads [16, 124], this
phenomenon was first reported in [43] where experiments showed that rubber
cylinders ruptured under relatively small tensile dead loads by opening an
internal cavity. The nonlinear elastic analysis in [10] provided the first
systematic theory for the formation of a spherical cavity at the center of a
sphere of isotropic incompressible hyperelastic material under a radially
symmetric tensile load. Mathematically, cavitation was treated as a
bifurcation from the trivial state at a critical value of the surface traction
or displacement, at which the trivial solution became unstable. This work
paved the way for numerous applied and theoretical studies devoted to this
inherently nonlinear mechanical effect, which is not captured by the linear
elasticity theory. In addition to the well-known stable cavitation post-
bifurcation at the critical dead load, such that the cavity radius
monotonically increases with the applied load, it was shown in [71] that
unstable (snap) cavitation is also possible for some homogeneous isotropic
incompressible hyperelastic materials where Backer-Ericksen inequalities [9,
66] hold. In [81], static and dynamic cavitation in radially-inhomogeneous
spheres with stochastic parameters were analyzed.

Figure 4: Probability distribution of the applied dead-load traction causing
cavitation of radius $c$ at the center of an initially unit sphere of
neoclassical LCE material when: (a) $a=2$, while the shear modulus $\mu$
follows a Gamma distribution with hyperparameters $\rho_{1}=405$,
$\rho_{2}=0.01$, and (b) $\mu=4.05$, while the shape parameter $a$ follows a
Gamma distribution with $\rho_{1}=200$, $\rho_{2}=0.01$. The dashed black
lines correspond to the deterministic solutions based only on the mean values
of the parameters, whereas the red versions show the arithmetic mean value
solutions. The probabilistic interval for the critical load at which the
cavity nucleates is found when $c\to 0$. As the bifurcation from trivial
solution is supercritical, cavitation is stable, i.e., the cavity radius
increases as the applied dead load increases.

To study the cavitation of a nematic sphere, we consider an initially unit
sphere described by the neoclassical model (3), where either the shear modulus
$\mu$ or the nematic parameter $a$ is a random variable. Expressing all
tensors in the spherical coordinates $(R,\Theta,\Phi)$ in the reference
configuration, we assume that the sphere is deformed by radially symmetric
inflation with deformation gradient
$\textbf{F}=\mathrm{{diag}}\left(\lambda^{-2},\lambda,\lambda\right)$, while
the natural deformation tensor is
$\textbf{G}=\mathrm{diag}\left(a^{-1/3},a^{1/6},a^{1/6}\right)$, and
$\lambda>a^{1/6}>1$. When this deformation is due to a radial dead-load
traction applied uniformly on the sphere surface in the reference
configuration, we are interested in the critical load that will cause a
spherical cavity to open at its center. The radial traction at the outer
surface necessary for an inner cavity of radius $c$ to form is illustrated in
Figure 4. For the onset of cavitation, the critical load is found by letting
$c\to 0$. Since the bifurcation is supercritical, cavitation is stable, i.e.,
the cavity radius increases as the applied dead load increases. A comparison
with [81] shows that, for the nematic sphere, cavitation nucleates at a larger
critical load, $P^{(nc)}_{0}$, than the corresponding load, $P_{0}$, for the
elastic sphere with the same shear modulus, and $P^{(nc)}_{0}=a^{1/3}P_{0}$.
We recall that, for the neo-Hookean sphere, $P_{0}=5\mu/2$. When the LCE
strain-energy function includes an additional elastic component, as in (1), so
that the elasto-nematic ratio is $\eta>0$, while the shear modulus $\mu$
remains the same, the critical dead load decreases towards $P_{0}$. For
example, if $\textbf{G}_{0}=\textbf{I}$, then the critical applied load is
equal to
$\overline{P}^{(nc)}_{0}=\left[a^{1/3}+\left(1-a^{1/3}\right)\mu^{(1)}/\mu\right]P_{0}$,
and $P_{0}<\overline{P}^{(nc)}_{0}<P^{(nc)}_{0}$.

##  5 Shell inflation

Figure 5: The normalized internal pressure for an inflated spherical shell of:
(a) Mooney-Rivlin hyperelastic material, defined by (10), when
$\textbf{A}=\mathrm{{diag}}\left(\alpha^{-2},\alpha,\alpha\right)$, with
$\alpha>1$; (b) Mooney-Rivlin LCE material, defined by (11), when
$\textbf{F}=\mathrm{{diag}}\left(\lambda^{-2},\lambda,\lambda\right)$ and
$\textbf{G}=\mathrm{diag}\left(a^{-1/3},a^{1/6},a^{1/6}\right)$, with $a=2$
and $\lambda>a^{1/6}\approx 1.1225$; (c) Mooney-Rivlin LCE material with
$a=2$, while the shear modulus $\mu$ follows a Gamma distribution with
$\rho_{1}=405$, $\rho_{2}=0.01$, and $R_{1}=\mu_{1}/\mu$ is drawn from a Beta
distribution with $\xi_{1}=900$, $\xi_{2}=100$; and (d) Mooney-Rivlin LCE
material with $\mu=4.05$ and $R_{1}=\mu_{1}/\mu=0.9$, while the shape
parameter $a$ follows a Gamma distribution with $\rho_{1}=200$,
$\rho_{2}=0.01$. In (a)-(b), when the parameter ratio $\mu_{1}/\mu$ is
sufficiently small, the required stress changes from increasing to decreasing,
i.e., the material displays inflation instability. For the nematic model,
instability is expected at larger deformation than for the underlying
hyperelastic model. In (c)-(d), the dashed black lines correspond to the
deterministic solutions based only on the mean values of the parameters,
whereas the red versions show the arithmetic mean value solutions. When the
material parameters are stochastic, the critical load for instability is found
in a probabilistic interval where both stable and unstable states are observed
with a given probability.

Internally pressurized hollow spheres and tubes are relevant in many
biological and engineering structures [47]. For rubber spherical and tubular
balloons, the first experimental observations of inflation instabilities under
internal pressure were reported in [65]. Cylindrical tubes of homogeneous
isotropic incompressible hyperelastic material subject to finite symmetric
inflation and stretching were theoretically analyzed for the first time in
[91]. Finite radially symmetric inflation of elastic spherical shells was
initially investigated in [48], then in [1, 95]. For both elastic tubular and
spherical shells, it was shown in [17] that, depending on the material model,
the internal pressure may increase monotonically, or increase and then
decrease, or increase, decrease, and then increase again. These classical
results were extended to elastic materials with stochastic parameters in [68,
69, 70]. Recent theoretical investigations of inflated nematic cylindrical
balloons were presented in [45, 62].

To compare _inflation instabilities_ in nematic and in purely elastic spheres,
we consider the hyperelastic Mooney-Rivlin model [85, 90], given by

$W(\textbf{A})=\frac{\mu_{1}}{2}\left[\mathrm{tr}\left(\textbf{A}\textbf{A}^{\text{T}}\right)-3\right]+\frac{\mu_{2}}{2}\left\\{\mathrm{tr}\left[\mathrm{Cof}\left(\textbf{A}\textbf{A}^{\text{T}}\right)\right]-3\right\\},$
(10)

where $\mu=\mu_{1}+\mu_{2}>0$ is the shear modulus at infinitesimal strain. A
Mooney-Rivlin-type neoclassical strain-energy function for the nematic
material then takes the form

$W^{(nc)}(\textbf{F},\textbf{n})=W(\textbf{G}^{-1}\textbf{F}).$ (11)

Taking a spherical coordinates system with coordinates $(R,\Theta,\Phi)$ in
the reference configuration, we assume that the sphere is deformed by radially
symmetric inflation with deformation gradient
$\textbf{F}=\mathrm{{diag}}\left(\lambda^{-2},\lambda,\lambda\right)$, while
the natural deformation tensor is
$\textbf{G}=\mathrm{diag}\left(a^{-1/3},a^{1/6},a^{1/6}\right)$, and
$\lambda>a^{1/6}>1$. We denote
$\mathcal{W}^{(nc)}(\lambda,\textbf{n})=W^{(nc)}(\textbf{F},\textbf{n})$, and
further assume that the shell is thin, i.e., $0<\epsilon=(B-A)/A\ll 1$, where
$A$ and $B$ represent the inner and outer radii of the reference shell,
respectively. When the deformation is due to a radial pressure applied
uniformly on the inner surface in the present configuration, the corresponding
radial Cauchy stress at the inner surface can be approximated as
$T^{(nc)}=\epsilon
a^{1/2}\lambda^{-2}\mathrm{d}\mathcal{W}^{(nc)}/\mathrm{d}\lambda$. The
relation between the Cauchy stress at the inner surface in the nematic and
purely elastic shell with the same shear modulus is $T^{(nc)}=T$. Then Figure
5(a)-(b) shows that, if the parameter ratio $\mu_{1}/\mu$ is sufficiently
small, the required stress changes from increasing to decreasing, i.e., the
material displays inflation instability, and for the nematic model,
instability is expected at larger deformation than for the hyperelastic model.
However, this value will decrease if the model is modified to include an
additional elastic energy as in (1), so that the elasto-nematic ratio is
$\eta>0$, while the shear modulus $\mu$ remains the same. For the nematic
shell, in Figure 5(c)-(d), either the shear modulus $\mu=\mu_{1}+\mu_{2}$ or
the shape parameter $a$ is a random variable. In both cases, the critical load
for instability resides in a probabilistic interval where the stable and
unstable states compete. To decrease the chance of stable deformation, one
must increase the value of $\mu_{1}/\mu$, and only when $\mu_{1}=\mu$ unstable
deformation is certain.

##  6 Necking

Figure 6: The effect of changing the value of parameter ratio $\mu_{1}/\mu$ on
the normalized tensile first Piola-Kirchhoff stress for: (a) the Gent-Thomas
model, defined by (12), when
$\textbf{A}=\mathrm{{diag}}\left(\alpha^{-1/2},\alpha,\alpha^{-1/2}\right)$
with $\alpha>1$; (b) the Gent-Thomas-type LCE model, defined by (13), when
$\textbf{F}=\mathrm{{diag}}\left(\lambda^{-1/2},\lambda,\lambda^{-1/2}\right)$
and $\textbf{G}=\mathrm{diag}\left(a^{-1/6},a^{1/3},a^{-1/6}\right)$, with
$a=2$ and $\lambda>a^{1/3}\approx 1.2599$; (c) the LCE model when $a=2$, while
the shear modulus $\mu$ follows a Gamma distribution with $\rho_{1}=405$,
$\rho_{2}=0.01$, and $R_{1}=\mu_{1}/\mu$ is drawn from a Beta distribution
with $\xi_{1}=12$, $\xi_{2}=100$; and (b) the LCE model when $\mu=4.05$ and
$R_{1}=\mu_{1}/\mu=0.11$, while the shape parameter $a$ follows a Gamma
distribution with $\rho_{1}=200$, $\rho_{2}=0.01$. In (c)-(d), the dashed
black lines correspond to the deterministic solutions based only on the mean
values of the parameters, whereas the red versions show the arithmetic mean
value solutions. In all cases, when the parameter ratio is sufficiently small,
the required dead load changes from increasing to decreasing, and the material
displays necking instability. For the nematic model, necking is expected at
larger deformation and lower maximum load than for the underlying hyperelastic
model. When the material parameters are stochastic, both stable and unstable
states are presented with a given probability.

Like other rubber-like materials, LCEs may suffer from _necking instability_
under stretch [19, 20]. When homogeneous isotropic incompressible hyperelastic
materials are subject to large tension, necking occurs if there is a maximum
load, or in other words, if there exists a critical extension ratio, such that
the force required to extend the material beyond this critical value changes
from increasing to decreasing [5, 6, 7, 21, 31, 53, 88]. The relation between
the onset of necking and the maximum load was originally analyzed for ductile
materials in [22]. For a class of hyperelastic materials where the load-
displacement curve does not possess a maximum, in [96], it was proved that the
homogeneous deformation is the only absolute minimizer of the elastic energy.
In particular, incompressible neo-Hookean or Mooney-Rivlin hyperelastic models
do not exhibit necking, and this property is inherited by the associated
nematic LCE models. For example, in [51], a necking instability observed
experimentally under uniaxial tension could not be captured by the
neoclassical LCE model based on the neo-Hookean strain-energy function alone.
In that case, since necking was initiated during director rotation, a
composite model consisting of a purely neoclassical form and an additional
elastic form, as presented in Section 3, should be used to predict the non-
zero stress plateau associated with the neck formation.

To further explore necking instability that may occur when the director is
parallel to the applied tensile force, and compare that with the same behavior
in a purely elastic material, we consider the hyperelastic Gent-Thomas model
[44] defined by

$W(\textbf{A})=\frac{\mu_{1}}{2}\left[\mathrm{tr}\left(\textbf{A}\textbf{A}^{\text{T}}\right)-3\right]+\frac{3\mu_{2}}{2}\ln\frac{\mathrm{tr}\left[\mathrm{Cof}\left(\textbf{A}\textbf{A}^{\text{T}}\right)\right]}{3},$
(12)

where $\mu=\mu_{1}+\mu_{2}>0$ is the shear modulus at small strain. A Gent-
Thomas-type neoclassical strain-energy function for the nematic material then
reads

$W^{(nc)}(\textbf{F},\textbf{n})=W(\textbf{G}^{-1}\textbf{F}).$ (13)

We take
$\textbf{F}=\mathrm{{diag}}\left(\lambda^{-1/2},\lambda,\lambda^{-1/2}\right)$,
while $\textbf{G}=\mathrm{diag}\left(a^{-1/6},a^{1/3},a^{-1/6}\right)$ and
$\lambda>a^{1/3}>1$. In this case, denoting
$\mathcal{W}^{(nc)}(\lambda,\textbf{n})=W^{(nc)}(\textbf{F},\textbf{n})$, the
uniaxial tensile load is given by the first Piola-Kirchhoff stress
$P^{(nc)}_{2}=\mathrm{d}\mathcal{W}^{(nc)}/\mathrm{d}\lambda$, and necking
occurs when the ratio $\mu_{1}/\mu$ is sufficiently small. The relation
between the first Piola-Kirchhoff tensile stress in the nematic and purely
elastic case with the same shear modulus is $P^{(nc)}_{2}=a^{-1/3}P_{2}$. As
shown in Figure 6(a)-(b), for the nematic model, necking is expected at larger
deformation and lower maximum dead load than for the hyperelastic model.
However, the maximum load will increase if the model is modified to include an
additional elastic energy as in (1), so that the elasto-nematic ratio is
$\eta>0$, while the shear modulus $\mu$ remains the same. In Figure 6(c)-(d),
the stochastic tensile load in the deformed LCE is illustrated.

##  7 Conclusion

Instabilities in liquid crystalline solids can be of potential interest in a
range of applications. Here, we present various situations to understand some
new possibilities offered by LCEs. Specifically, we combine for the first time
nonlinear elastic instabilities with nematic elastomer theory, within a
stochastic setting where the material parameters are probabilistic, and
compare the results with those from classical nonlinear elasticity. While
shear striping instabilities are specific to LCEs, as they do not occur in
rubber, instabilities such as cavitation, shell inflation and necking have
been widely studied in the context of purely hyperelastic materials. In recent
years, we have reviewed such results and extended them to a stochastic
elasticity framework by building directly on the deterministic nonlinear
theory. For LCEs, similar interesting phenomena can occur, and we offer a
perspective on the general methodology to analyze them. Moreover, each of
these instability problems could be further formulated and treated under
slightly different hypotheses, where the external forces or the type of
nematic material might change to reflect different experimental setups and
observations. Our scope is to present some key theoretical ingredients to
study this class of problems.

The particular choice of numerical values in our calculations are for
illustrative purposes only. To compare the stochastic results with the
deterministic ones, we sampled from distributions where the parameters have
mean values corresponding to the deterministic system. The stochastic results
mostly follow the deterministic ones, but transform a single critical value
for instability into a large probability of an instability taking place close
to that value. However, both stable and unstable states have a quantifiable
chance to be observed with a given probability, and small variations in the
input model parameters can have a significant impact on whether an instability
occurs or not at a certain load. Moreover, we observed that, as deformation
progresses, the solution variance tends to change non-uniformly around the
mean value, suggesting that the average value may be less significant from a
physical point of view if fluctuations become large.

The results presented here are universal in the sense that they hold for
families of LCE models for which their classic hyperelastic counterparts
exhibit a similar instability. To gain significant insight into the
macroscopic mechanical properties of these materials, it is imperative that
they are analyzed and tested under more complex multiaxial deformations and
loads before they can be incorporated into real industrial systems. Therefore,
we hope that the general field of instabilities for LCE materials may serve as
an inspiration for new devices or for systematic testing. For example, in [52]
inflation instabilities were reported for elongated nematic balloons. The
inflation of cylindrical elastic balloons is often accompanied by different
types of instabilities, such as bulging and necking. These phenomena have been
analyzed theoretically for hyperelastic balloons in [41, 42, 123], but remain
virgin territory for LCEs. The possibility of activation of LCEs may finally
fulfill the early promises of the rational mechanics pioneers who first
demonstrated the existence of nonlinear elastic instabilities and understood
that they could be used to design new devices.

##  Appendix A Stochastic model parameters

We refer to [50, 49, 97] for extensive reviews on the information-theoretic
approach in stochastic elasticity, and to the various papers cited in the main
text for numerous examples and applications of continuum models with
stochastic parameters. Here, we adopt the following hypothesis required by the
stochastic models for nematic liquid crystal elastomers (LCEs) presented in
[73]: For any given finite deformation, at any point in the material, the
shear modulus $\mu>0$, the shape parameter $a>0$, and their inverse, $1/\mu$
and $1/a$, respectively, are second order random variables, i.e., they have
finite mean value and finite variance. For the shear modulus $\mu$ (and
similarly for $a$), to construct a prior probability law, we note that this
assumption is guaranteed by setting the following mathematical expectations:

$\displaystyle\begin{cases}E\left[\mu\right]=\underline{\mu}>0,&\\\
E\left[\log\ \mu\right]=\nu,&\mbox{such that $|\nu|<+\infty$}.\end{cases}$
(A.1)

The first constraint in (A.1) specifies the mean value for the random shear
modulus $\mu$, while the second constraint provides a condition from which it
follows that $1/\mu$ is a second order random variable. Then, by the maximum
entropy principle, the shear modulus $\mu$ with mean value $\underline{\mu}$
and standard deviation $\|\mu\|=\sqrt{\text{Var}[\mu]}$ (defined as the square
root of the variance, $\text{Var}[\mu]$) follows a Gamma probability
distribution with shape and scale parameters $\rho_{1}>0$ and $\rho_{2}>0$
respectively, such that

$\underline{\mu}=\rho_{1}\rho_{2},\qquad\|\mu\|=\sqrt{\rho_{1}}\rho_{2}.$
(A.2)

The corresponding probability density function takes the form

$g(\mu;\rho_{1},\rho_{2})=\frac{\mu^{\rho_{1}-1}e^{-\mu/\rho_{2}}}{\rho_{2}^{\rho_{1}}\Gamma(\rho_{1})},\qquad\mbox{for}\
\mu>0\ \mbox{and}\ \rho_{1},\rho_{2}>0,$ (A.3)

where $\Gamma:\mathbb{R}^{*}_{+}\to\mathbb{R}$ is the complete Gamma function

$\Gamma(z)=\int_{0}^{+\infty}t^{z-1}e^{-t}\text{d}t.$ (A.4)

The word ‘hyperparameters’ is often used for $\rho_{1}$ and $\rho_{2}$ to
distinguish them from $\mu$ and other material constants.

When $\mu=\mu_{1}+\mu_{2}$, setting a fixed constant value $b>-\infty$, such
that $\mu_{i}>b$, $i=1,2$ (e.g., $b=0$ if $\mu_{1}>0$ and $\mu_{2}>0$,
although $b$ is not unique in general), we define the auxiliary random
variable

$R_{1}=\frac{\mu_{1}-b}{\mu-2b},$ (A.5)

such that $0<R_{1}<1$. Then, the random model parameters can be expressed
equivalently as follows,

$\mu_{1}=R_{1}(\mu-2b)+b,\qquad\mu_{2}=\mu-\mu_{1}=(1-R_{1})(\mu-2b)+b.$ (A.6)

It is reasonable to assume

$\displaystyle\begin{cases}E\left[\log\ R_{1}\right]=\nu_{1},&\mbox{such that
$|\nu_{1}|<+\infty$},\\\ E\left[\log(1-R_{1})\right]=\nu_{2},&\mbox{such that
$|\nu_{2}|<+\infty$},\end{cases}$ (A.7)

in which case, the random variable $R_{1}$ follows a standard Beta
distribution, with hyperparameters $\xi_{1}>0$ and $\xi_{2}>0$ satisfying

$\underline{R}_{1}=\frac{\xi_{1}}{\xi_{1}+\xi_{2}},\qquad\text{Var}[R_{1}]=\frac{\xi_{1}\xi_{2}}{\left(\xi_{1}+\xi_{2}\right)^{2}\left(\xi_{1}+\xi_{2}+1\right)},$
(A.8)

where $\underline{R}_{1}$ is the mean value and $\text{Var}[R_{1}]$ is the
variance of $R_{1}$. The associated probability density function is

$\beta(r;\xi_{1},\xi_{2})=\frac{r^{\xi_{1}-1}(1-r)^{\xi_{2}-1}}{B(\xi_{1},\xi_{2})},\qquad\qquad\mbox{for}\
r\in(0,1)\ \mbox{and}\ \xi_{1},\xi_{2}>0,$ (A.9)

where $B:\mathbb{R}^{*}_{+}\times\mathbb{R}^{*}_{+}\to\mathbb{R}$ is the Beta
function

$B(x,y)=\int_{0}^{1}t^{x-1}(1-t)^{y-1}dt.$ (A.10)

Then, for the random coefficients given by (A.6), the corresponding mean
values are

$\underline{\mu}_{1}=\underline{R}_{1}(\underline{\mu}-2b)+b,\qquad\underline{\mu}_{2}=\underline{\mu}-\underline{\mu}_{1}=(1-\underline{R}_{1})(\underline{\mu}-2b)+b,$
(A.11)

and the variances and covariance take the form, respectively,

$\displaystyle\text{Var}\left[\mu_{1}\right]=(\underline{\mu}-2b)^{2}\text{Var}[R_{1}]+(\underline{R}_{1})^{2}\text{Var}[\mu]+\text{Var}[\mu]\text{Var}[R_{1}],$
(A.12)
$\displaystyle\text{Var}\left[\mu_{2}\right]=(\underline{\mu}-2b)^{2}\text{Var}[R_{1}]+(1-\underline{R}_{1})^{2}\text{Var}[\mu]+\text{Var}[\mu]\text{Var}[R_{1}],$
(A.13)
$\displaystyle\text{Cov}[\mu_{1},\mu_{2}]=\frac{1}{2}\left(\text{Var}[\mu]-\text{Var}[\mu_{1}]-\text{Var}[\mu_{2}]\right).$
(A.14)

##  Appendix B Stress tensors for ideal nematic elastomers

We briefly recall the relations between the stress tensors of an ideal nematic
elastomer and those of the underlying hyperelastic model. These relations were
originally obtained in [74]. For an ideal incompressible nematic elastomer,
the neoclassical strain-energy density function takes the general form

$W^{(nc)}(\textbf{F},\textbf{n})=W(\textbf{A}),$ (B.1)

where the right-hand side represents the strain-energy function of a
homogeneous isotropic incompressible hyperelastic material, depending only on
the elastic deformation gradient A. On the left-hand side, n is a unit vector
for the localized direction of uniaxial nematic alignment in the present
configuration; $\textbf{F}=\textbf{G}\textbf{A}$ is the deformation gradient
tensor with respect to the reference isotropic state, with
$\textbf{G}=a^{-1/6}\textbf{I}+\left(a^{1/3}-a^{-1/6}\right)\textbf{n}\otimes\textbf{n}$
the ‘spontaneous’ (or ‘natural’) deformation tensor and A the (local) elastic
deformation tensor; $a>0$ is a temperature-dependent, spatially-independent
shape parameter; $\otimes$ denotes the tensor product of two vectors; and
$\textbf{I}=\text{diag}(1,1,1)$ is the identity tensor.

The strain-energy function given by (B.1) takes the equivalent form

$\mathcal{W}^{(nc)}(\lambda_{1},\lambda_{2},\lambda_{3},\textbf{n})=W^{(nc)}(\textbf{F},\textbf{n}),$
(B.2)

where $\\{\lambda_{i}^{2}\\}_{i=1,2,3}$ are the eigenvalues of the tensor
$\textbf{F}\textbf{F}^{T}$.

For the hyperelastic material described by the strain-energy function
$W(\textbf{A})$, the Cauchy stress tensor (representing the internal force per
unit of deformed area acting within the deformed solid) is equal to

$\textbf{T}=\left(\det\textbf{A}\right)^{-1}\frac{\partial
W}{\partial\textbf{A}}\textbf{A}^{T}-p\textbf{I},$ (B.3)

where $p$ denotes the Lagrange multiplier for the incompressibility constraint
$\det\textbf{A}=1$.

The associated first Piola-Kirchhoff stress tensor (representing the internal
force per unit of undeformed area acting within the deformed solid) is then

$\textbf{P}=\textbf{T}\text{Cof}(\textbf{A}),$ (B.4)

where $\text{Cof}(\textbf{A})=\left(\det\textbf{A}\right)\textbf{A}^{-T}$ is
the cofactor of A.

###  B.1 Free director

When the nematic director is ‘free’ to rotate relative to the elastic matrix,
F and n are independent variables, and the Cauchy stress tensor for the
nematic material with the strain-energy function described by (B.1) is
calculated as follows,

$\begin{split}\textbf{T}^{(nc)}&=J^{-1}\frac{\partial
W^{(nc)}}{\partial\textbf{F}}\textbf{F}^{T}-p^{(nc)}\textbf{I}\\\
&=J^{-1}\textbf{G}^{-1}\frac{\partial
W}{\partial\textbf{A}}\textbf{A}^{T}\textbf{G}-p^{(nc)}\textbf{I}\\\
&=J^{-1}\textbf{G}^{-1}\textbf{T}\textbf{G},\end{split}$ (B.5)

where T is the Cauchy stress tensor defined by (B.3), $J=\det\textbf{F}$, and
the scalar $p^{(nc)}$ represents the Lagrange multiplier for the internal
constraint $J=1$.

The principal components $\left(T^{(nc)}_{1},T^{(nc)}_{2},T^{(nc)}_{3}\right)$
of the Cauchy stress defined by (B.5) are the solutions of the characteristic
equation

$\det\left(\textbf{T}^{(nc)}-\Lambda\textbf{I}\right)=0.$ (B.6)

Since

$\det\left(\textbf{T}^{(nc)}-\Lambda\textbf{I}\right)=\det\left[\textbf{G}^{-1}\left(J^{-1}\textbf{T}-\Lambda\textbf{I}\right)\textbf{G}\right]=J^{-1}\det\left(\textbf{T}-J\Lambda\textbf{I}\right),$
(B.7)

it follows that the principal Cauchy stresses for the underlying hyperelastic
model satisfy

$\left(T_{1},T_{2},T_{3}\right)=J\left(T^{(nc)}_{1},T^{(nc)}_{2},T^{(nc)}_{3}\right).$
(B.8)

Therefore, if the Baker-Ericksen inequalities hold for the hyperelastic model,
then the greater principal Cauchy stress occurs in the direction of the
greater principal elastic stretch for the nematic model. We recall that, for a
hyperelastic material, the Baker-Ericksen inequalities state that the greater
principal stress occurs in the direction of the greater principal stretch [9,
66].

In the presence of a nematic field, the total Cauchy stress tensor
$\textbf{T}^{(nc)}$ given by (B.5) is not symmetric in general. In addition,
the following condition is required,

$\frac{\partial W^{(nc)}}{\partial\textbf{n}}=\textbf{0},$ (B.9)

or equivalently, by the principle of material objectivity,

$\frac{1}{2}\left(\textbf{T}^{(nc)}-\textbf{T}^{(nc)T}\right)\textbf{n}=\textbf{0},$
(B.10)

where $\left(\textbf{T}^{(nc)}-\textbf{T}^{(nc)^{T}}\right)/2$ represents the
skew-symmetric part of the Cauchy stress tensor.

The first Piola-Kirchhoff stress tensor for the nematic material is equal to

$\textbf{P}^{(nc)}=\textbf{T}^{(nc)}\text{Cof}(\textbf{F})=\textbf{G}^{-1}\textbf{T}\textbf{A}^{-T}=\textbf{G}^{-1}\textbf{P},$
(B.11)

where P is the first Piola-Kirchhoff stress given by (B.4).

###  B.2 Frozen director

If the nematic director is ‘frozen’, the Cauchy stress tensor for the nematic
material takes the form

$\widehat{\textbf{T}}^{(nc)}=J^{-1}\textbf{G}^{-1}\textbf{T}\textbf{G}-J^{-1}q\left(\textbf{I}-\frac{\textbf{F}\textbf{n}_{0}\otimes\textbf{F}\textbf{n}_{0}}{|\textbf{F}\textbf{n}_{0}|^{2}}\right)\textbf{n}\otimes\frac{\textbf{F}\textbf{n}_{0}}{|\textbf{F}\textbf{n}_{0}|},$
(B.12)

where T is the Cauchy stress defined by (B.3), $J=\det\textbf{F}$, $p^{(nc)}$
is the Lagrange multiplier for the volume constraint $J=1$, and $q$ is the
Lagrange multiplier for the constraint

$\textbf{n}=\frac{\textbf{F}\textbf{n}_{0}}{|\textbf{F}\textbf{n}_{0}|}.$
(B.13)

As the Cauchy stress tensor given by (B.12) is not symmetric in general, the
following additional condition must hold,

$\frac{\partial\widehat{W}^{(nc)}}{\partial\textbf{n}}=\textbf{0},$ (B.14)

or equivalently,

$\frac{1}{2}\left(\widehat{\textbf{T}}^{(nc)}-\widehat{\textbf{T}}^{(nc)T}\right)\textbf{n}=\textbf{0}.$
(B.15)

The corresponding first Piola-Kirchhoff stress tensor for the nematic material
is equal to

$\widehat{\textbf{P}}^{(nc)}=\widehat{\textbf{T}}^{(nc)}\text{Cof}(\textbf{F}).$
(B.16)

##  Appendix C Cavitation of a nematic sphere

The static and dynamic cavitation of homogeneous and radially inhomogeneous
isotropic incompressible hyperelastic spheres with stochastic material
parameters was analyzed in [71, 81] where up-to-date reviews of the relevant
literature are presented. For an inflated elastic sphere, the radially
symmetric deformation takes the form

$r=g(R),\qquad\theta=\Theta,\qquad\phi=\Phi,$ (C.1)

where $(R,\Theta,\Phi)$ and $(r,\theta,\phi)$ are the spherical polar
coordinates in the reference and current configuration, respectively, such
that $0\leq R\leq B$, and $g(R)\geq 0$ is to be determined. The corresponding
deformation gradient is equal to
$\textbf{A}=\mathrm{diag}\left(\alpha_{1},\alpha_{2},\alpha_{3}\right)$, with

$\alpha_{1}=\frac{\mathrm{d}g}{\mathrm{d}R}=\alpha^{-2},\qquad\alpha_{2}=\alpha_{3}=\frac{g(R)}{R}=\frac{r}{R}=\alpha,$
(C.2)

where $\alpha_{1}$ and $\alpha_{2}=\alpha_{3}$ are the radial and hoop
principal stretches, respectively, and $\mathrm{d}g/\mathrm{d}R$ denotes the
derivative of $g$ with respect to $R$. By (C.2),

$g^{2}\frac{\mathrm{d}g}{\mathrm{d}R}=R^{2},$ (C.3)

hence,

$g(R)=\left(R^{3}+c^{3}\right)^{1/3},$ (C.4)

where $c\geq 0$ is a constant to be calculated. If $c>0$, then $g(R)\to c>0$
as $R\to 0_{+}$, and a spherical cavity of radius $c$ forms at the center of
the sphere, from zero initial radius, otherwise the sphere remains undeformed.

In particular, for a static sphere of neo-Hookean material, if the surface of
the cavity is traction-free, the radial component of the Cauchy stress is
equal to (see [81] for a detailed derivation using the same notation)

$T_{rr}(b)=\frac{2}{3}\int_{x^{3}+1}^{\infty}\mu\frac{1+u}{u^{7/3}}\mathrm{d}u,$
(C.5)

where $x=c/B$. After evaluating the integral in (C.5), the required uniform
dead-load traction at the outer surface, $R=B$, in the reference
configuration, takes the form

$P=\left(x^{3}+1\right)^{2/3}T_{rr}(b)=2\mu\left[\left(x^{3}+1\right)^{1/3}+\frac{1}{4\left(x^{3}+1\right)^{2/3}}\right],$
(C.6)

and increases as $x$ increases. The critical dead load for the onset of
cavitation is then

$P_{0}=\lim_{x\to 0_{+}}P=\frac{5\mu}{2}.$ (C.7)

To analyze the stability of this cavitation, we study the behavior of the
cavity opening, with radius $c$ as a function of $P$, in a neighborhood of
$P_{0}$. After differentiating the function given by (C.6), with respect to
the dimensionless cavity radius $x=c/B$, we have

$\frac{\mathrm{d}P}{\mathrm{d}x}=2\mu
x^{2}\left[\frac{1}{\left(x^{3}+1\right)^{2/3}}-\frac{1}{2\left(x^{3}+1\right)^{5/3}}\right]>0,$
(C.8)

i.e., the cavitation is stable, regardless of the material parameter $\mu>0$.

For a nematic sphere of neoclassical material with the strain-energy function
given by (B.1) derived from the neo-Hookean hyperelastic model, if
$\textbf{F}=\mathrm{{diag}}\left(\lambda^{-2},\lambda,\lambda\right)$ and
$\textbf{G}=\mathrm{diag}\left(a^{-1/3},a^{1/6},a^{1/6}\right)$, with
$\lambda>a^{1/6}>1$, then the Cauchy stress is equal to that of the neo-
Hookean sphere. Hence, $T^{(nc)}_{rr}=T_{rr}$, and the first Piola-Kirchhoff
stress representing the critical dead load for cavitation in the nematic
sphere is equal to $P^{(nc)}_{0}=a^{1/3}P_{0}=5a^{1/3}\mu/2$.

##  Appendix D Inflation of a nematic spherical shell

Static and dynamic inflation instabilities of homogeneous and radially
inhomogeneous isotropic incompressible hyperelastic spheres with stochastic
material parameters were analyzed in [68, 69, 70]. For a thin hyperelastic
spherical shell, such that $0<\epsilon=(B-A)/A\ll 1$, where $A$ and $B$
represent the inner and outer radii of the reference shell, respectively, if
the external pressure is equal to zero, then the internal pressure can be
approximated as follows,

$T=\frac{\epsilon}{\alpha^{2}}\frac{\text{d}\mathcal{W}}{\text{d}\alpha},$
(D.1)

where the deformation gradient for radially symmetric inflation is equal to
$\textbf{A}=\mathrm{{diag}}\left(\alpha^{-2},\alpha,\alpha\right)$, with
$\alpha=r/R>1$, and $\mathcal{W}(\alpha)=W(\textbf{A})$. The critical value of
$\alpha$ where a limit-point instability occurs is obtained by solving for
$\alpha>1$ the following equation,

$\frac{\text{d}T}{\text{d}\alpha}=0,$ (D.2)

where $T$ is the radial component of the Cauchy stress given by (D.1).

For a nematic sphere of neoclassical material with the strain-energy function
given by (B.1) derived from the Mooney-Rivlin hyperelastic model, when
$\textbf{F}=\mathrm{{diag}}\left(\lambda^{-2},\lambda,\lambda\right)$ and
$\textbf{G}=\mathrm{diag}\left(a^{-1/3},a^{1/6},a^{1/6}\right)$, with
$\lambda>a^{1/6}>1$, the Cauchy stress is equal to that of the hyperelastic
sphere. If the shell is thin, assuming zero external pressure, the internal
pressure can be approximated as

$T^{(nc)}=T=\frac{\epsilon}{\alpha^{2}}\frac{\text{d}\mathcal{W}}{\text{d}\alpha}=\frac{\epsilon
a^{1/2}}{\lambda^{2}}\frac{\mathrm{d}\mathcal{W}^{(nc)}}{\mathrm{d}\lambda}.$
(D.3)

Then, the critical value of $\lambda$ where a limit-point instability occurs
is found by solving for $\lambda>a^{1/3}$ the equation

$\frac{\text{d}T^{(nc)}}{\text{d}\lambda}=0,$ (D.4)

with $T^{(nc)}$ given by (D.3).

#### Acknowledgement.

We gratefully acknowledge the support by the Engineering and Physical Sciences
Research Council of Great Britain under research grants EP/R020205/1 to Alain
Goriely and EP/S028870/1 to L. Angela Mihai.

## References

  * [1] Adkins JE, Rivlin RS. 1952. Large elastic deformations of isotropic materials. IX. The deformation of thin shells, Philosophical Transactions of the Royal Society of London A 244, 505-531.
  * [2] Agostiniani V, Dal Maso G, DeSimone A. 2015. Attainment results for nematic elastomers, Proceedings of the Royal Society of Edinburgh A 145, 669-701 (doi: 10.1017/S0308210515000128).
  * [3] Agostiniani V, DeSimone A. 2012. Ogden-type energies for nematic elastomers, International Journal of Non-Linear Mechanics 47(2), 402-412 (doi: 10.1016/j.ijnonlinmec.2011.10.001).
  * [4] Anderson DR, Carlson DE, Fried E. 1999. A continuum-mechanical theory for nematic elastomers, Journal of Elasticity 56, 33-58 (doi: 10.1023/A:1007647913363).
  * [5] Antman SS. 1973. Nonuniqueness of equilibrium states for bars in tension, Journal of Mathematical Analysis and Applications 44(2), 333-349 (doi: 10.1016/0022-247X(73)90063-2).
  * [6] Antman SS, Carbone ER. 1977. Shear and necking instabilities in nonlinear elasticity, Journal of Elasticity 7(2), 125-151 (doi: 10.1007/BF00041087).
  * [7] Audoly B, Hutchinson JW. 2016. Analysis of necking based on a one-dimensional model, Journal of the Mechanics and Physics of Solids 97, 68-91 (doi: 10.1016/j.jmps.2015.12.018).
  * [8] Bai R, Bhattacharya K. 2020. Photomechanical coupling in photoactive nematic elastomers, Journal of the Mechanics and Physics of Solids 144, 104115 (doi: 10.1016/j.jmps.2020.104115).
  * [9] Baker M, Ericksen JL. 1954. Inequalities restricting the form of stress-deformation relations for isotropic elastic solids and Reiner-Rivlin fluids, Journal of the Washington Academy of Sciences 44(2), 33-35 (www.jstor.org/stable/24533303).
  * [10] Ball JM. 1982. Discontinuous equilibrium solutions and cavitation in nonlinear elasticity, Philosophical Transactions of the Royal Society A 306(1496), 557-611 (doi: 10.1098/rsta.1982.0095).
  * [11] Ball JM, James RD. 1987. Fine phase mixtures as minimizers of energy, Archive for Rational Mechanics and Analysis 100, 13-52.
  * [12] Ball JM, James RD. 1992. Proposed experimental tests of a theory of fine microstructure, and the two-well problem, Philosophical Transactions of the Royal Society of London A 338, 389-450.
  * [13] Biggins JS, Warner M, Bhattacharya K. 2009. Supersoft elasticity in polydomain nematic elastomers, Physical Review Letters 103, 037802 (doi: 10.1103/PhysRevLett.103.037802).
  * [14] Biggins JS, Warner M, Bhattacharya K. 2012. Elasticity of polydomain liquid crystal elastomers, Journal of the Mechanics and Physics of Solids 60, 573-590 (doi: 10.1016/j.jmps.2012.01.008).
  * [15] Bladon P, Terentjev EM, Warner M. 1994. Deformation-induced orientational transitions in liquid crystal elastomers, Journal de Physique II 4, 75-91 (doi: 10.1051/jp2:1994100).
  * [16] Busse WF. 1939. Physics of rubber as related to the automobile, Journal of Applied Physics 9(7), 438-451 (doi: 10.1063/ 1.1710439).
  * [17] Carroll MM. 1987. Pressure maximum behavior in inflation of incompressible elastic hollow spheres and cylinders. Quarterly of Applied Mathematics 45, 141-154.
  * [18] Cirak F, Long Q, Bhattacharya K, Warner M. 2014. Computational analysis of liquid crystalline elastomer membranes: Changing Gaussian curvature without stretch energy, International Journal of Solids and Structures 51(1), 144-153 (doi: 10.1016/j.ijsolstr.2013.09.019).
  * [19] Clarke SM, Terentjev EM. 1998. Slow stress relaxation in randomly disordered nematic elastomers and gels, Physical Review Letters 81(20), 4436-4439 (doi: 10.1103/PhysRevLett.81.4436).
  * [20] Clarke SM, Terentjev EM, Kundler I, Finkelmann H. 1998. Texture evolution during the polydomain-monodomain transition in nematic elastomers, Macromolecules 31(15), 4862-4872 (doi: 10.1021/ma980195j).
  * [21] Coleman BD, Newman DC. 1988. On the rheology of cold drawing. I. Elastic materials, Journal of Polymer Science: Part B: Polymer Physics 26, 1801-1822 (doi: 10.1002/polb.1988.090260901).
  * [22] Considére A. 1885. Mémoire sur l’emploi du fer et de l’acier dans les constructions, Annales des Ponts et Chaussées 6(9), 574-775.
  * [23] Conti S, DeSimone A, Dolzmann G. 2002. Soft elastic response of stretched sheets of nematic elastomers: a numerical study, Journal of the Mechanics and Physics of Solids 50, 1431-1451.
  * [24] de Gennes PG. 1975. Physique moléculaire - réflexions sur un type de polymères nématiques, Comptes rendus de l’Académie des Sciences B 281, 101-103.
  * [25] de Gennes PG, Prost J. 1993. The Physics of Liquid Crystals, Clarendon Press, Oxford.
  * [26] de Haan LT, Schenning AP, Broer DJ. 2014. Programmed morphing of liquid crystal networks, Polymer 55(23), 5885-5896 (doi: 10.1016/j.polymer.2014.08.023).
  * [27] DeSimone A. 1999. Energetics of fine domain structures, Ferroelectrics 222(1), 275-284 (doi: 10.1080/00150199908014827).
  * [28] DeSimone A, Dolzmann G. 2000. Material instabilities in nematic elastomers, Physica D 136(1-2), 175-191 (doi: S0167-2789(99)00153-0).
  * [29] DeSimone A, Dolzmann G. 2002. Macroscopic response of nematic elastomers via relaxation of a class of SO(3)-invariant energies, Archive of Rational Mechanics and Analysis 161, 181-204 (doi: 10.1007/s002050100174).
  * [30] DeSimone A, Teresi L. 2009. Elastic energies for nematic elastomers, The European Physical Journal E 29, 191-204 (doi: 10.1140/epje/i2009-10467-9).
  * [31] Ericksen JL. 1975. Equilibrium of bars, Journal of Elasticity 5(3-4), 191-201 (doi: 10.1007/BF00126984).
  * [32] Finkelmann H, Kock HJ, Rehage G. 1981. Investigations on liquid crystalline polysiloxanes 3, Liquid crystalline elastomers - a new type of liquid crystalline material, Die Makromolekulare Chemie, Rapid Communications 2, 317-322 (doi: 10.1002/marc.1981.030020413).
  * [33] Finkelmann H, Kundler I, Terentjev EM, Warner M. 1997. Critical stripe-domain instability of nematic elastomers, Journal de Physique II 7, 1059-1069 (doi: 10.1051/jp2:1997171).
  * [34] Fitt D, Wyatt H, Woolley TE, Mihai LA. 2019. Uncertainty quantification of elastic material responses: testing, stochastic calibration and Bayesian model selection, Mechanics of Soft Materials 1, 13 (doi: 10.1007/s42558-019-0013-1).
  * [35] Flory PJ. 1961. Thermodynamic relations for high elastic materials, Transactions of the Faraday Society 57, 829-838 (doi: 10. 1039/TF9615700829).
  * [36] Frank FC. 1958. I. Liquid crystals. On the theory of liquid crystals, Discussions of the Faraday Society 25, 19-28.
  * [37] Fried E. 2002. An elementary molecular-statistical basis for the Mooney and Rivlin-Saunders theories of rubber elasticity, Journal of the Mechanics and Physics of Solids 50, 571-582 (doi: 10.1016/S0022-5096(01)00086-2).
  * [38] Fried E, Sellers S. 2004. Free-energy density functions for nematic elastomers, Journal of the Mechanics and Physics of Solids 52(7), 1671-1689 (doi: 10.1016/j.jmps.2003.12.005).
  * [39] Fried E, Sellers S. 2005. Orientational order and finite strain in nematic elastomers, The Journal of Chemical Physics 123(4), 043521 (doi: 10.1063/1.1979479).
  * [40] Fried E, Sellers S. 2006. Soft elasticity is not necessary for striping in nematic elastomers, Journal of Applied Physics 100, 043521 (doi: 10.1063/1.2234824).
  * [41] Fu YB, Liu JL, Francisco GS. 2016. Localized bulging in an inflated cylindrical tube of arbitrary thickness - the effect of bending stiffness, Journal of the Mechanics and Physics of Solids 90, 45-60.
  * [42] Fu Y, Jin L, Goriely A. 2021. Necking, beading, and bulging in soft elastic cylinders, Journal of the Mechanics and Physics of Solids 147, 104250 (doi: 10.1016/j.jmps.2020.104250).
  * [43] Gent AN, Lindley PB. 1959. Internal rupture of bonded rubber cylinders in tension, Proceedings of the Royal Society of London A 249, 195-205.
  * [44] Gent AN, Thomas AG. 1958. Forms for the stored (strain) energy function for vulcanized rubber, Journal of Polymer Science 28, 625-628.
  * [45] Giudici A, Biggins JS. 2020. Giant deformations and soft-inflation in LCE balloons, Europhysics Letters 132(3), 36001 (doi: 10.1209/0295-5075/132/36001).
  * [46] Goodbrake C, Goriely A, Yavari A. 2021. The mathematical foundations of anelasticity: existence of smooth global intermediate configurations, Proceedings of the Royal Society A 477, 20200462 (doi: 10.1098/rspa.2020.0462).
  * [47] Goriely A. 2017. The Mathematics and Mechanics of Biological Growth, Springer-Verlag, New York.
  * [48] Green AE, Shield RT. 1950. Finite elastic deformations in incompressible isotropic bodies, Proceeding of the Royal Society of London A 202, 407-419.
  * [49] Guilleminot J. 2020. Modelling non-Gaussian random fields of material properties in multiscale mechanics of materials, in Uncertainty Quantification in Multiscale Materials Modeling, Y. Wang and F.L. McDowell (Eds.), Elsevier, Cambridge, MA.
  * [50] Guilleminot J, Soize C. 2017. Non-Gaussian random fields in multiscale mechanics of heterogeneous materials, in Encyclopedia of Continuum Mechanics, Springer, Berlin Heidelberg.
  * [51] He X, Zheng Y, He Q, Cai S. 2020. Uniaxial tension of a nematic elastomer with inclined mesogens,Extreme Mechanics Letters 40, 100936 (doi: 10.1016/j.eml.2020.100936).
  * [52] He Q, Zheng Y, Wang Z, He X, Cai S. 2020. Anomalous inflation of a nematic balloon, Journal of the Mechanics and Physics of Solids 142, 104013 (doi: 10.1016/j.jmps.2020.104013).
  * [53] Hill R. Hutchinson JW. 1975. Bifurcation phenomena in the plane tension test, Journal of the Mechanics and Physics of Solids 23(4-5), 239-264 (doi: 10.1016/0022-5096(75)90027-7).
  * [54] Horgan CO, Saccomandi G. 2002. A molecular–statistical basis for the Gent constitutive model of rubber elasticity, Journal of Elasticity 68, 167-176 (doi: 10.1023/A:1026029111723).
  * [55] Jaynes ET. 1957. Information theory and statistical mechanics i, Physical Review 108, 171-190.
  * [56] Jaynes ET. 1957. Information theory and statistical mechanics ii, Physical Review 106, 620-630.
  * [57] Jaynes ET. 2003. Probability Theory: The Logic of Science, Cambridge University Press, Cambridge, New York.
  * [58] Jiang ZC, Xiao YY, Zhao Y. 2019. Shining light on liquid crystal polymer networks: preparing, reconfiguring, and driving soft actuators, Advanced Optical Materials 7, 1900262 (doi: 10.1002/adom.201900262).
  * [59] Kuenstler AS. Hayward RC. 2019. Light-induced shape morphing of thin films, Current Opinion in Colloid & Interface Science 40, 70-86 (doi: 10.1016/j.cocis.2019.01.009).
  * [60] Kundler I, Finkelmann H. 1995. Strain-induced director reorientation in nematic liquid single crystal elastomers, Macromolecular Rapid Communications 16, 679-686 (doi: 10.1002/marc.1995.030160908).
  * [61] Küpfer J, Finkelmann H. 1991. Nematic liquid single crystal elastomers, Die Makromolekulare Chemie, Rapid Communications 12, 717-726 (doi: 10.1002/marc.1991.030121211).
  * [62] Lee V, Bhattacharya K. 2020. Actuation of cylindrical nematic elastomer balloons, arXiv:2012.12873.
  * [63] Lubarda VA. 2004. Constitutive theories based on the multiplicative decomposition of deformation gradient: thermoelasticity, elastoplasticity and biomechanics, Applied Mechanics Reviews 57(2), 95-108 (doi: 10.1115/1.1591000).
  * [64] Mahimwalla Z, Yager KG, Mamiya J-i, Shishido A, Priimagi A, Barrett CJ. 2012. Azobenzene photomechanics: prospects and potential applications, Polymer Bulletin 69, 967-1006 (doi: 10.1007/s00289-012-0792-0).
  * [65] Mallock A. Note on the instability of India-rubber tubes and balloons when distended by fluid pressure, Proceedings of the Royal Society of London 49, 458 (doi: 10.1098/rspl.1890.0116).
  * [66] Marzano M. 1983. An interpretation of Baker-Ericksen inequalities in uniaxial deformation and stress, Meccanica 18, 233-235 (doi: 10.1007/BF02128248).
  * [67] McCracken JM, Donovan BR, White TJ. 2020. Materials as Machines, Advanced Materials 32, 1906564 (doi: 10.1002/adma.201906564).
  * [68] Mihai LA, Alamoudi M. 2021. Likely oscillatory motions of stochastic hyperelastic spherical shells and tubes, International Journal of Non-Linear Mechanics 130, 103671 (doi: 10.1016/j.ijnonlinmec.2021.103671).
  * [69] Mihai LA, Fitt D, Woolley TE, Goriely A. 2019. Likely equilibria of stochastic hyperelastic spherical shells and tubes, Mathematics and Mechanics of Solids 24(7), 2066-2082 (doi: 10.1177/1081286518811881).
  * [70] Mihai LA, Fitt D, Woolley TE, Goriely A. 2019. Likely oscillatory motions of stochastic hyperelastic solids, Transactions of Mathematics and Its Applications 3, 1-42 (doi: 10.1093/imatrm/tnz003).
  * [71] Mihai LA, Fitt D, Woolley TE, Goriely A. 2019. Likely cavitation in stochastic elasticity, Journal of Elasticity 137(1), 27-42 (doi: 10.1007/s10659-018-9706-1).
  * [72] Mihai LA, Goriely A. 2017. How to characterize a nonlinear elastic material? A review on nonlinear constitutive parameters in isotropic finite elasticity, Proceedings of the Royal Society A 473, 20170607 (doi: 10.1098/rspa.2017.0607).
  * [73] Mihai LA, Goriely A. 2020. Likely striping in stochastic nematic elastomers, Mathematics and Mechanics of Solids 25(10), 1851-1872 (doi: 0.1177/1081286520914958).
  * [74] Mihai LA, Goriely A. 2020. A plate theory for nematic liquid crystalline solids, Journal of the Mechanics and Physics of Solids 144, 104101 (doi: 10.1016/j.jmps.2020.104101).
  * [75] Mihai LA, Goriely A. 2020. A pseudo-anelastic model for stress softening in liquid crystal elastomers, Proceedings of the Royal Society A, 20200558 (doi: 10.1098/rspa.2020.0558).
  * [76] Mihai LA, Neff P. 2017. Hyperelastic bodies under homogeneous Cauchy stress induced by non-homogeneous finite deformations, International Journal of Non-Linear Mechanics 89, 93-100.
  * [77] Mihai LA, Neff P. 2017. Hyperelastic bodies under homogeneous Cauchy stress induced by three-dimensional non-homogeneous deformations, Mathematics and Mechanics of Solids (doi: 10.1177/1081286516682556).
  * [78] Mihai LA, Woolley TE, Goriely A. 2018. Stochastic isotropic hyperelastic materials: constitutive calibration and model selection, Proceedings of the Royal Society A 474, 20170858 (doi: 10.1098/rspa.2017.0858).
  * [79] Mihai LA, Woolley TE, Goriely A. 2019. Likely equilibria of the stochastic Rivlin cube, Philosophical Transactions of the Royal Society A 377, 20180068 (doi: 10.1098/rsta.2018.0068).
  * [80] Mihai LA, Woolley TE, Goriely A. 2019. Likely chirality of stochastic anisotropic hyperelastic tubes, International Journal of Non-Linear Mechanics 114, 9-20 (doi: 10.1016/j.ijnonlinmec.2019.04.004).
  * [81] Mihai LA, Woolley TE, Goriely A. 2020. Likely cavitation and radial motion of stochastic elastic spheres, Nonlinearity 33(5), 1987-2034 (doi: 10.1088/1361-6544/ab7104).
  * [82] Mistry D, Gleeson HF. 2019. Mechanical deformations of a liquid crystal elastomer at director angles between 0∘ and 90∘: Deducing an empirical model encompassing anisotropic nonlinearity, Journal of Polymer Science 57, 1367-1377 (doi: 0.1002/polb.24879).
  * [83] Modes CD, Bhattacharya K, Warner M. 2010. Disclination-mediated thermo-optical response in nematic glass sheets, Physical Review E 81, 060701(R) (doi: 10.1103/PhysRevE.81.060701).
  * [84] Modes CD, Bhattacharya K, Warner M. 2011. Gaussian curvature from flat elastica sheets, Proceedings of the Royal Society A 467, 1121-1140 (doi:10.1098/rspa.2010.0352).
  * [85] Mooney M. 1940. A theory of large elastic deformation, Journal of Applied Physics 11, 582-592.
  * [86] Okamoto S, Sakurai S, Urayama K. 2021. Effect of stretching angle on the stress plateau behavior of main-chain liquid crystal elastomers, Soft Matter (doi: 10.1039/d0sm02244f).
  * [87] Olmsted P. 1994. Rotational invariance and Goldstone modes in nematic elastomers and gels, Journal de Physique II, EDP Sciences 4(12), 2215-2230 (doi: 10.1051/jp2:1994257).
  * [88] Owen N. 1987. Existence and stability of necking deformations for nonlinearly elastic rods, Archive for Rational Mechanics and Analysis 98, 357-383 (doi: 10.1007/BF00276914).
  * [89] Pang X, Lv J-a., Zhu C, Qin L, Yu Y. 2019. Photodeformable azobenzene-containing liquid crystal polymers and soft actuators, Advanced Materials, 1904224 (doi: 10.1002/adma.201904224).
  * [90] Rivlin RS. 1948. Large elastic deformations of isotropic materials. IV. Further developments of the general theory, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 241, 379-397.
  * [91] Rivlin RS. 1949. Large elastic deformations of isotropic materials. VI. Further results in the theory of torsion, shear and flexure, Philosophical Transactions of the Royal Society of London A 242(845), 173-195.
  * [92] Sadik S, Yavari A. 2017. On the origins of the idea of the multiplicative decomposition of the deformation gradient, Mathematics and Mechanics of Solids 22(4), 771-772 (doi: 10.1177/1081286515612280).
  * [93] Saed MO, Torbati AH, Starr CA, Visvanathan R, Clark NA, Yakacki CM. 2017. Thiol-acrylate main-chain liquid-crystalline elastomers with tunable thermomechanical properties and actuation strain, Journal of Polymer Science 55(2), 157-168 (doi: 0.1002/polb.24249).
  * [94] Shannon CE. 1948. A mathematical theory of communication, Bell System Technical Journal 27, 379-423, 623-659.
  * [95] Shield RT. 1972. On the stability of finitely deformed elastic membranes. II: Stability of inflated cylindrical and spherical membranes, Zeitschrift für Angewandte Mathematik und Physik (ZAMP) 23, 16-34.
  * [96] Sivaloganathan J, Spector SJ. 2011. On the stability of incompressible elastic cylinders in uniaxial extension, Journal of Elasticity 105(1-2), 313-330 (doi: 10.1007/s10659-011-9330-9).
  * [97] Soize C. 2017. Uncertainty Quantification: An Accelerated Course with Advanced Applications in Computational Engineering, Interdisciplinary Applied Mathematics Book 47, Springer, New York.
  * [98] Staber B, Guilleminot J. 2015. Stochastic modeling of a class of stored energy functions for incompressible hyperelastic materials with uncertainties, Comptes Rendus Mécanique 343, 503-514 (doi: 10.1016/j.crme.2015.07.008).
  * [99] Staber B, Guilleminot J. 2016. Stochastic modeling of the Ogden class of stored energy functions for hyperelastic materials: the compressible case, Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik 97, 273-295 (doi: 10.1002/zamm.201500255).
  * [100] Staber B, Guilleminot J. 2017. Stochastic hyperelastic constitutive laws and identification procedure for soft biological tissues with intrinsic variability, Journal of the Mechanical Behavior of Biomedical Materials 65, 743-752 (doi: 10.1016/j.jmbbm.2016.09.022).
  * [101] Staber B, Guilleminot J. 2018. A random field model for anisotropic strain energy functions and its application for uncertainty quantification in vascular mechanics, Computer Methods in Applied Mechanics and Engineering 333, 94-113 (doi: 10.1016/j.cma.2018.01.001).
  * [102] Staber B, Guilleminot J, Soize C, Michopoulos J, Iliopoulos A. 2019. Stochastic modeling and identification of an hyperelastic constitutive model for laminated composites, Computer Methods in Applied Mechanics and Engineering 347, 425-444 (doi: 10.1016/j.cma.2018.12.036).
  * [103] Talroze RV, Zubarev ER, Kuptsov SA, Merekalov AS, Yuranova TI, Plate NA, Finkelmann H. 1999. Liquid crystal acrylate-based networks: polymer backbone-LC order interaction, Reactive and Functional Polymers 41, 1-11 (doi: 10.1016/S1381-5148(99)00032-2).
  * [104] Traugutt NA, Volpe RH, Bollinger MS, Saed MO, Torbati AH, Yu K, Dadivanyan N, Yakacki CM. 2017. Liquid-crystal order during synthesis affects main-chain liquid-crystal elastomer behavior, Soft Matter 13, 7013 (doi: 10.1039/c7sm01405h).
  * [105] Treloar LRG. 1944. Stress-strain data for vulcanized rubber under various types of deformation, Transactions of the Faraday Society 40, 59-70 (doi: 10.1039/TF9444000059).
  * [106] Treloar LRG. 2005. The Physics of Rubber Elasticity, 3rd ed, Oxford University Press, Oxford, UK.
  * [107] Ube T, Ikeda T. 2014. Photomobile polymer materials with crosslinked liquid-crystalline structures: molecular design, fabrication, and functions, Angewandte Chemie International Edition 53(39), 10290–10299 (doi: 10.1002/anie.201400513).
  * [108] Ula SW, Traugutt NA, Volpe RH, Patel RP, Yu K, Yakacki CM. 2018. Liquid crystal elastomers: an introduction and review of emerging technologies, Liquid Crystals Review 6, 78-107 (doi: 10.1080/21680396.2018.1530155).
  * [109] Urayama K, Kohmon E, Kojima M, Takigawa T. 2009. Polydomain-monodomain transition of randomly disordered nematic elastomers with different cross-linking histories, Macromolecules 42, 4084-4089 (doi: 10.1021/ma9004692).
  * [110] Verwey GC, Warner M. 1995. Soft rubber elasticity, Macromolecules 28, 4303-4306.
  * [111] Wan G, Jin C, Trase I, Zhao S, Chen Z. 2018. Helical structures mimicking chiral seedpod opening and tendril coiling, Sensors 18(9), 2973 (doi: 10.3390/s18092973).
  * [112] Warner M. 2020. Topographic mechanics and applications of liquid crystalline solids, Annual Review of Condensed Matter Physics 11, 125-145 (doi: 10.1146/annurev-conmatphys-031119-050738).
  * [113] Warner M, Bladon P, Terentjev E. 1994. “Soft elasticity” - deformation without resistance in liquid crystal elastomers, Journal de Physique II France 4, 93-102 (doi: 10.1051/jp2:1994116).
  * [114] Warner M, Gelling KP, Vilgis TA. 1988. Theory of nematic networks, The Journal of Chemical Physics 88, 4008-4013 (doi: 10.1063/1.453852).
  * [115] Warner M, Terentjev EM. 1996. Nematic elastomers - a new state of matter?, Progress in Polymer Science 21, 853-891.
  * [116] Warner M, Terentjev EM. 2007. Liquid Crystal Elastomers, paper back, Oxford University Press, Oxford, UK.
  * [117] Warner M, Wang, XJ. 1991. Elasticity and phase behavior of nematic elastomers, Macromolecules 24, 4932-4941 (doi: 10.1021/ma00017a033).
  * [118] Wen Z, Yang K, Raquez JM. 2020. A review on liquid crystal polymers in free-standing reversible shape memory materials, Molecules 25, 1241 (doi:10.3390/molecules25051241).
  * [119] Verwey GC, Warner M, Terentjev EM. 1996. Elastic instability and stripe domains in liquid crystalline elastomers, Journal de Physique 6(9), 1273-1290 (doi: 10.1051/jp2:1996130).
  * [120] White TJ. 2018. Photomechanical effects in liquid crystalline polymer networks and elastomers, Journal of Polymer Science, Part B: Polymer Physics 56, 695-705 (doi: 10.1002/polb.24576).
  * [121] White TJ, Broer DJ. 2015. Programmable and adaptive mechanics with liquid crystal polymer networks and elastomers, Nature Materials 14, 1087-1098 (doi: 10.1038/nmat4433).
  * [122] Xia Y, Honglawan A, Yang S. 2019. Tailoring surface patterns to direct the assembly of liquid crystalline materials, Liquid Crystals Reviews 7(1) 30-59 (doi: 10.1080/21680396.2019.1598295).
  * [123] Ye Y, Liu Y, Fu Y. 2020. Weakly nonlinear analysis of localized bulging of an inflated hyperelastic tube of arbitrary wall thickness, Journal of the Mechanics and Physics of Solids 135, 103804 (doi: 10.1016/j.jmps.2019.103804).
  * [124] Yerzley FL. 1939. Adhesion of neoprene to metal, Industrial & Engineering Chemistry 31(8), 950-956 (doi: 10.1021/ ie50356a007).
  * [125] Zhang Y, Xuan C, Jiang Y, Huo Y. 2019. Continuum mechanical modeling of liquid crystal elastomers as dissipative ordered solids, Journal of the Mechanics and Physics of Solids 126, 285-303 (doi: 10.1016/j.jmps.2019.02.018).
  * [126] Zubarev ER, Kuptsov SA, Yuranova TI, Talroze RV, Finkelmann H. 1999. Monodomain liquid crystalline networks: reorientation mechanism from uniform to stripe domains, Liquid Crystals 26, 1531-1540 (doi: 10.1080/026782999203869).