Instabilities in liquid crystal elastomers

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. Because of 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, thus 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.


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 [2,3,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 [49,50,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.

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,113,114,117,119,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], 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; F = GA is the deformation gradient tensor with respect to the reference isotropic state (see Figure 1 and also Figure 1 of [30]), with G = a 1/3 n ⊗ n + a −1/6 (I − n ⊗ n) the 'spontaneous' (or 'natural') deformation tensor and A the (local) elastic deformation tensor; is the spontaneous deformation tensor with 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 ⊗ the usual tensor product of two vectors, and by I = 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-71, 78-81, 98-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.

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?
For simplicity, we select an incompressible neo-Hookean-type strain-energy function [105], where the superscript "T" represents the transpose operator, "tr" denotes the trace operator, and µ (1) ≥ 0 is constant, together with the neoclassical strain-energy function [18,28,30], with µ (2) ≥ 0 constant.For the composite model function defined by (1), the shear modulus at infinitesimal strain is µ = µ (1) + µ (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 n 0 = [0, 0, 1] T and n = [0, sin θ, cos θ] T , where θ ∈ [0, π/2], in a Cartesian system of reference, and examine small shear perturbations of biaxial extensions, with gradient tensor [73] where a > 1 is the nematic shape parameter, λ > 0 is the stretch ratio, and ε > 0 is the small perturbation.Denoting w(λ, ε, θ) = W (nc) (F, n), for sufficiently small values of ε and θ, we find that, when µ (1) > 0 and µ (2) = 0 (purely elastic case), the equilibrium state with ε = 0 and given by equations ( 5) and ( 6), respectively.Between these bounds, the second solution, with (ε, θ) = (ε 0 , θ 0 ), minimizes the energy.The lower and upper bounds on λ, given by (c) ( 5)-( 6), and (d) ( 8)-( 9), respectively, as functions of the parameter ratio η = µ (1) /µ (2) when a = 2.For η → ∞, the model approaches a purely elastic form where the homogeneous deformation is always stable.For η = 0, the bounds are the same, and correspond to the neoclassical form.θ = 0 is stable.If µ (2) > 0, then this equilibrium state is unstable for where the elasto-nematic ratio η = µ (1) /µ (2) gives the relative magnitude of the elastic and neoclassical contributions.There is also an equilibrium state with ε = 0 and θ = π/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 When λ satisfies ( 5) or ( 6), we obtain For the resulting strip pattern, the gradient tensors of alternating shear deformations in two adjacent stripe domains are F ± with ε = ±ε 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 F + − F − = q ⊗ p is satisfied, where p is the normal vector to the interface between the two phases corresponding to the deformation gradients F + and F − [11,12,76,77].
In other words, F + and F − are rank-one connected, i.e., rank (F + − F − ) = 1.If η = 0, then the above equilibrium states are unstable for λ ∈ a −1/6 , a 1/12 and λ ∈ a 1/12 , a 1/3 , respectively.Thus, soft elasticity is always presented by the purely neoclassical model [23,30].When η → ∞, there is no shear striping since the material is practically elastic.The strain-energy function w(λ, ε, θ) is illustrated in Figure 2(a)-(b).Note that, for λ with values between the lower and upper bounds given by ( 5) and (6), respectively, the minimum energy is attained for (ε, θ) = (ε 0 , θ 0 ).Assuming that loading is applied in the second direction, the first Piola-Kirchhoff axial stress in this direction is equal to 2(b) then suggests that, when µ (1) = 0, i.e., for the purely neoclassical form, the director rotates and alternating shear stripes develop for λ ∈ a −1/6 , a 1/3 , at zero load, since the slope of the curve is equal to zero within this interval.In contrast, if µ (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 λ where shear striping occurs are at a maximum length when n 0 = [0, 0, 1] T , the bounds ( 5)- (6) give also the maximum interval for shear striping when n 0 is not uniformly aligned.The minimum length of those intervals is attained for monodomains with n 0 = [0, 1, 0] T .Experimental results for monodomains where the tensile load forms different angles with the initial nematic director are reported in [82,86].If G 0 = I (see [13]), then the solution with ε = 0 and θ = 0 is unstable for and that with ε = 0 and θ = π/2 is unstable for For example, when a = 2, the bounds given by ( 5)-( 6), and by ( 8)-( 9), respectively, are plotted as functions of the parameter ratio η 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 µ and the shape parameter a separately, i.e., we first set µ as a random variable while a is a single-valued constant, then keep µ constant and let a fluctuate.In Figure 3(a)-(b), the stochastic shear parameter ε 0 and director angle θ 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 η, whereas shear striping is certain only if the model reduces to the neoclassical one (i.e., when η = 0).However, while η > 0, the inherent variability in the probabilistic system means that there will always be competition between the homogeneous and inhomogeneous deformations.

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.
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 µ or the nematic parameter a is a random variable.Expressing all tensors in the spherical coordinates (R, Θ, Φ) in the reference configuration, we assume that the sphere is deformed by radially symmetric inflation with deformation gradient F = diag λ −2 , λ, λ , while the natural deformation tensor is G = diag a −1/3 , a 1/6 , a 1/6 , and λ > 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 → 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µ/2.When the LCE strain-energy function includes an additional elastic component, as in (1), so that the elasto-nematic ratio is η > 0, while the shear modulus µ remains the same, the critical dead load decreases towards P 0 .For example, if G 0 = I, then the critical applied load is equal to P (nc) 0 = a 1/3 + 1 − a 1/3 µ (1) /µ P 0 , and P 0 < P The probabilistic interval for the critical load at which the cavity nucleates is found when c → 0. As the bifurcation from trivial solution is supercritical, cavitation is stable, i.e., the cavity radius increases as the applied dead load increases.

Shell inflation
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 where µ = µ 1 + µ 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 Taking a spherical coordinates system with coordinates (R, Θ, Φ) in the reference configuration, we assume that the sphere is deformed by radially symmetric inflation with deformation gradient F = diag λ −2 , λ, λ , while the natural deformation tensor is G = diag a −1/3 , a 1/6 , a 1/6 , and λ > a 1/6 > 1.We denote W (nc) (λ, n) = W (nc) (F, n), and further assume that the shell is , when the parameter ratio µ 1 /µ 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.
thin, i.e., 0 < = (B − A)/A 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) = a 1/2 λ −2 dW (nc) /dλ.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 µ 1 /µ 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 η > 0, while the shear modulus µ remains the same.For the nematic shell, in Figure 5(c)-(d), either the shear modulus µ = µ 1 + µ 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 µ 1 /µ, and only when µ 1 = µ unstable deformation is certain.

Necking
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-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 where µ = µ 1 + µ 2 > 0 is the shear modulus at small strain.A Gent-Thomas-type neoclassical strain-energy function for the nematic material then reads We take F = diag λ −1/2 , λ, λ −1/2 , while G = diag a −1/6 , a 1/3 , a −1/6 and λ > a 1/3 > 1.In this case, denoting W (nc) (λ, n) = W (nc) (F, n), the uniaxial tensile load is given by the first Piola-Kirchhoff stress P (nc) 2 = dW (nc) /dλ, and necking occurs when the ratio µ 1 /µ 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 η > 0, while the shear modulus µ remains the same.In Figure 6(c)-(d), the stochastic tensile load in the deformed LCE is illustrated.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.

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.

A Stochastic model parameters
We refer to [49,50,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 µ > 0, the shape parameter a > 0, and their inverse, 1/µ and 1/a, respectively, are second order random variables, i.e., they have finite mean value and finite variance.For the shear modulus µ (and similarly for a), to construct a prior probability law, we note that this assumption is guaranteed and the variances and covariance take the form, respectively,

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 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; F = GA is the deformation gradient tensor with respect to the reference isotropic state, with G = a −1/6 I + a 1/3 − a −1/6 n ⊗ n the 'spontaneous' (or 'natural') deformation tensor and A the (local) elastic deformation tensor; a > 0 is a temperaturedependent, spatially-independent shape parameter; ⊗ denotes the tensor product of two vectors; and I = diag(1, 1, 1) is the identity tensor.The strain-energy function given by (B.1) takes the equivalent form where {λ 2 i } i=1,2,3 are the eigenvalues of the tensor FF T .For the hyperelastic material described by the strain-energy function W (A), the Cauchy stress tensor (representing the internal force per unit of deformed area acting within the deformed solid) is equal to where p denotes the Lagrange multiplier for the incompressibility constraint det 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 where Cof(A) = (det A) 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, where T is the Cauchy stress tensor defined by (B.3), J = det F, and the scalar p (nc) represents the Lagrange multiplier for the internal constraint J = 1.
The principal components T , T , T of the Cauchy stress defined by (B.5) are the solutions of the characteristic equation det T it follows that the principal Cauchy stresses for the underlying hyperelastic model satisfy , T .

(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 T (nc) given by (B.5) is not symmetric in general.In addition, the following condition is required, or equivalently, by the principle of material objectivity, where T (nc) − T (nc) T /2 represents the skew-symmetric part of the Cauchy stress tensor.The first Piola-Kirchhoff stress tensor for the nematic material is equal to 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 where T is the Cauchy stress defined by (B.3), J = det F, p (nc) is the Lagrange multiplier for the volume constraint J = 1, and q is the Lagrange multiplier for the constraint As the Cauchy stress tensor given by (B.12) is not symmetric in general, the following additional condition must hold, ∂ W (nc) ∂n = 0, (B.The corresponding first Piola-Kirchhoff stress tensor for the nematic material is equal to 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 where (R, Θ, Φ) and (r, θ, φ) are the spherical polar coordinates in the reference and current configuration, respectively, such that 0 ≤ R ≤ B, and g(R) ≥ 0 is to be determined.The corresponding deformation gradient is equal to A = diag (α 1 , α 2 , α 3 ), with where α 1 and α 2 = α 3 are the radial and hoop principal stretches, respectively, and dg/dR denotes the derivative of g with respect to R. By (C.2), where c ≥ 0 is a constant to be calculated.If c > 0, then g(R) → c > 0 as R → 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) i.e., the cavitation is stable, regardless of the material parameter µ > 0.

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 < = (B − A)/A 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, where the deformation gradient for radially symmetric inflation is equal to A = diag α −2 , α, α , with α = r/R > 1, and W(α) = W (A). The critical value of α where a limit-point instability occurs is obtained by solving for α > 1 the following equation, 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 F = diag λ −2 , λ, λ and G = diag a −1/3 , a 1/6 , a 1/6 , with λ > 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 Then, the critical value of λ where a limit-point instability occurs is found by solving for λ > a 1/3 the equation dT (nc) dλ = 0, (D.4) with T (nc) given by (D.3).

Figure 4 :
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 µ follows a Gamma distribution with hyperparameters ρ 1 = 405, ρ 2 = 0.01, and (b) µ = 4.05, while the shape parameter a follows a Gamma distribution with ρ 1 = 200, ρ 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 → 0. As the bifurcation from trivial solution is supercritical, cavitation is stable, i.e., the cavity radius increases as the applied dead load increases.

Figure 5 :
Figure 5: The normalized internal pressure for an inflated spherical shell of: (a) Mooney-Rivlin hyperelastic material, defined by (10), when A = diag α −2 , α, α , with α > 1; (b) Mooney-Rivlin LCE material, defined by(11), when F = diag λ −2 , λ, λ and G = diag a −1/3 , a 1/6 , a 1/6 , with a = 2 and λ > a 1/6 ≈ 1.1225; (c) Mooney-Rivlin LCE material with a = 2, while the shear modulus µ follows a Gamma distribution with ρ 1 = 405, ρ 2 = 0.01, and R 1 = µ 1 /µ is drawn from a Beta distribution with ξ 1 = 900, ξ 2 = 100; and (d) Mooney-Rivlin LCE material with µ = 4.05 and R 1 = µ 1 /µ = 0.9, while the shape parameter a follows a Gamma distribution with ρ 1 = 200, ρ 2 = 0.01.In (a)-(b), when the parameter ratio µ 1 /µ 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.

Figure 6 :
Figure6: The effect of changing the value of parameter ratio µ 1 /µ on the normalized tensile first Piola-Kirchhoff stress for: (a) the Gent-Thomas model, defined by(12), when A = diag α −1/2 , α, α −1/2 with α > 1; (b) the Gent-Thomas-type LCE model, defined by(13), when F = diag λ −1/2 , λ, λ −1/2 and G = diag a −1/6 , a 1/3 , a −1/6 , with a = 2 and λ > a 1/3 ≈ 1.2599; (c) the LCE model when a = 2, while the shear modulus µ follows a Gamma distribution with ρ 1 = 405, ρ 2 = 0.01, and R 1 = µ 1 /µ is drawn from a Beta distribution with ξ 1 = 12, ξ 2 = 100; and (b) the LCE model when µ = 4.05 and R 1 = µ 1 /µ = 0.11, while the shape parameter a follows a Gamma distribution with ρ 1 = 200, ρ 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.

7 )
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 formP = x 3 + 1 2/3 T rr (b) = 2µ x 3 + 1 1/3 + 1 4 (x 3 + 1) 2/3 , (C.6)and increases as x increases.The critical dead load for the onset of cavitation is thenPTo 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,