Nematic liquid crystalline elastomers are aeolotropic materials

Continuum models describing ideal nematic solids are widely used in theoretical studies of liquid crystal elastomers. However, experiments on nematic elastomers show a type of anisotropic response that is not predicted by the ideal models. Therefore, their description requires an additional term coupling elastic and nematic responses, to account for aeolotropic effects. In order to better understand the observed elastic response of liquid crystal elastomers, we analyse theoretically and computationally different stretch and shear deformations. We then compare the elastic moduli in the infinitesimal elastic strain limit obtained from the molecular dynamics simulations with the ones derived theoretically, and show that they are better explained by including nematic order effects within the continuum framework.


Introduction
Liquid crystalline solids are responsive multifunctional materials that combine the flexibility of polymeric networks with the nematic order of liquid crystals [1,2]. Owing to their molecular architecture, they can exhibit dramatic spontaneous deformations and phase transitions, which are reversible and repeatable under heat, light, solvents and electric or magnetic fields [3][4][5][6][7][8][9][10][11][12]. However, their physical behaviour under combined mechanical loading and external stimuli still needs to be fully elucidated.
For ideal monodomain nematic elastomers, with the mesogens uniaxially aligned throughout the material, a continuum model is given by the so-called neoclassical strain-energy density introduced in [13][14][15][16]. This is a phenomenological strain-energy function that extends the neo-Hookean model for rubber [17], where the parameters are derived from macroscopic shape changes at small strain or through statistical averaging at microscopic scale [18,19]. Since elastic stresses dominate over Frank elasticity induced by the distortion of mesogens alignment [20][21][22], Frank energy [23,24] is generally neglected. Extensions to polydomains where every domain has the same strain-energy density as a monodomain are provided in [25,26]. These descriptions have been generalized by employing other hyperelastic models, such as Mooney-Rivlin, Gent and Ogden, which better capture the nonlinear elastic behaviour at large strains [27][28][29] (see [30,31] for molecular interpretations of the Mooney-Rivlin and Gent strain energies in rubber elasticity). Further generalizations can be found in [32,33]. Numerical studies of liquid crystal elastomers (LCEs) are presented, for example, in [33][34][35], where the finite-element method is used, and in [36,37], where molecular Monte Carlo simulations are employed.
Usually, when the elastic properties of a material are investigated, uniaxial deformations, which are easier to reproduce experimentally, are examined first [38][39][40][41]. For a purely elastic isotropic material, the shear modulus is then inferred from a universal relation between elastic moduli from the classical theory. To study the elastic responses of nematic monodomains, uniaxial deformations were assumed in [38] (see also [42]), where it was found that, if only the elastic energy was considered, then the stretch moduli in the direction parallel to the director and in a perpendicular direction were equal. However, experiments clearly show an aeolotropic effect; namely, the stretch moduli depend on the direction in which they are measured. To capture this experimentally observed response of the material, the nematic energy [43] was then also taken into account. Experimental results for monodomains where the tensile load formed different angles with the initial nematic director were reported in [39,44]. In [45], measurements of five independent elastic constants derived from three uniaxial tests, with the director parallel, perpendicular or at an angle of 45°relative to the loading direction, respectively, were obtained for a nematic monodomain treated as a classical transversely isotropic material. However, for many complex materials, shear deformations can reveal important additional mechanical effects, which may not be observed or inferred from uniaxial tests. In particular, to assess LCE materials, shear deformations with the direction of shear either parallel or perpendicular to the nematic director need to be considered independently of uniaxial stretches [46][47][48][49][50][51][52]. For example, on the one hand, it was found experimentally in [38] that, if E and E ⊥ are the stretch moduli in a direction parallel or perpendicular to the nematic director, respectively, then E /E ⊥ > 1 at low temperature, E /E ⊥ < 1 at high temperature and E /E ⊥ = 1 at the transition point. On the other hand, if μ and μ ⊥ denote the shear moduli in a direction parallel or perpendicular to the nematic director, respectively, then it was reported in [47,48,51] that μ /μ ⊥ ≈ 1 in the isotropic phase and μ /μ ⊥ < 1 at temperatures below that for the nematic-isotropic phase transition. Therefore, from a symmetry point of view, monodomain LCEs are transversely isotropic materials with five independent elastic constants and the distinguished direction given by the nematic director. However, despite the constitutive symmetry about the direction given by the nematic field, the mechanical responses of LCEs differ from the known elastic behaviours in traditional transversely isotropic elastic materials where, typically, E > E ⊥ and μ > μ ⊥ [53][54][55][56].
The aim of this study is to develop an explicit approach for the derivation of elastic moduli that captures the aeolotropy of liquid crystalline elastomers. This approach represents an extension of the general theoretical framework by which similar elastic moduli were obtained for hyperelastic materials [57]. In the case of nematic solids, these moduli include information about both the elasticity of the polymeric network and the mechanical responses of the liquid crystal molecules. In §2, we recall the neoclassical model for ideal nematic elastomers, with the isotropic phase at high temperature as the reference configuration [29,[58][59][60][61], instead of the nematic phase at cross-linking [14][15][16]32,33,62,63]. Phenomenologically, this choice is motivated by the multiplicative decomposition of the effective deformation into an elastic distortion, followed by a natural stress-free shape change [64][65][66][67]. This multiplicative decomposition is similar to those found in the constitutive theories of thermoelasticity, elastoplasticity and growth [68,69] (see also [70,71]), but it is also different in the sense that the elastic deformation is directly applied to the reference state. The elastic stresses can then be used to study the final deformation where the stress-free geometrical change also plays a role [65]. In §3, we calculate the two stretch moduli under small elastic uniaxial tension and finite natural deformation when the nematic director is either parallel or perpendicular to the tensile direction, respectively. In §4, we further consider three shear deformations where the elastic component is a simple shear, while the nematic director is either parallel to the shear direction, perpendicular to the shear direction or perpendicular to the shear plane, respectively. When the elastic shear strain is small and the natural deformation is finite, we obtain effective shear moduli with the relative ratio equal to the natural anisotropy parameter of the nematic material. Note that we use the uniaxial and simple shear deformation, respectively, to describe the elastic contribution to the deformation rather than the overall deformation, which also contains a natural deformation component. To derive the elastic moduli, we then take the limit of small elastic strain, while the natural deformation remains finite. This enables us to rigorously adapt the elasticity theory to nematic elastomers (for a review on elastic moduli, see [57]). To account for the physical aeolotropy of real nematic solids, in §5, we extend the continuum model by incorporating a nematic energy, and show how the stretch and shear moduli corresponding to the ideal case are modified by the additional information. The parameters entering the LCE model characterize either the change of the microstructure due to nematic effects or the behaviour of the material in large natural deformations. Nevertheless, in the limit of small elastic deformations and for fixed nematic parameters, the system behaves indeed like a transverse isotropic material with five independent elastic constants. In physics, aeolotropy refers to materials exhibiting different properties depending on the direction in which they are measured, or simply defined by Lord Kelvin as 'That which is different in different directions' [72, p. 122]. While traditional anisotropic elastic materials also exhibit aeolotropy, we refer to nematic solids as aeolotropic materials, and reserve the characterization of isotropic or anisotropic for the elastic part of the energy. In §6, we present a molecular dynamics simulation of a nematic elastomer, and analyse its response under similar stretch and shear deformations as for the continuum model to illustrate the aeolotropic mechanical responses. In § §3-5, physical quantities are treated symbolically, and units of measure only appear in §6, where datasets are used to illustrate the theory. The final section contains concluding remarks.

Prerequisites
The strain-energy density describing an ideal monodomain nematic liquid crystalline (NLC) solid takes the general form [64][65][66][67]73] where F represents the deformation gradient from the isotropic state, n is a unit vector, known as the director, for the orientation of the nematic field and W(A) denotes the strain-energy density of the isotropic polymer network, depending only on the (local) elastic deformation tensor A. The tensors F and A satisfy the following relation: . is the spontaneous deformation tensor describing a change of frame of reference from the isotropic to a nematic phase (e.g. [74]). In (2.3), a > 0 represents a temperature-dependent stretch parameter, ν is the optothermal analogue to the Poisson ratio [74] and relates responses in directions parallel or perpendicular to the director n, ⊗ denotes the tensor product of two vectors and I = diag (1, 1, 1) is the identity tensor. It is assumed here that a and ν are spatially independent. The ratio r = a 1/3 /a −ν/3 = a (ν+1)/3 represents the anisotropy parameter, which, in an ideal nematic solid, is the same in all directions. In the nematic phase, both the cases with r > 1 (prolate molecules) and r < 1 (oblate molecules) are possible, while when r = 1 the energy function reduces to that of an isotropic hyperelastic material. Nematic elastomers have ν = 1/2, while for nematic glasses ν ∈ (1/2, 2) [10,75]. Natural strains in NLC glasses are typically of up to 4%, whereas for NLC elastomers these may be up to 400%. The nematic director n is an observable (spatial) quantity. Denoting by n 0 the reference orientation of the local director corresponding to the cross-linking state, n may differ from n 0 both by a rotation and by a change in r. In nematic elastomers, which are weakly cross-linked, the director can rotate freely, and the material exhibits isotropic mechanical effects. In nematic glasses, which are densely cross-linked, the director n cannot rotate relative to the elastic matrix, but changes through convection due to elastic strain and satisfies [21,22,58,74,76] This constraint enables patterning of the director field at cross-linking and guarantees that the 'written-in' pattern remains virtually the same during natural shape changes [21,22,77]. For a hyperelastic material described by the strain-energy density W = W(A), the Cauchy stress tensor is equal to where the so-called hydrostatic pressure p denotes the Lagrange multiplier for the incompressibility constraint det A = 1, β 1 = 2∂W/∂I 1 and β −1 = −2∂W/∂I 2 are material parameters, B = AA T is the left Cauchy-Green elastic deformation tensor and I 1 , I 2 are its first two principal invariants (I 3 = 1 owing to incompressibility). The corresponding first Piola-Kirchhoff stress tensor is equal to P = TCof(A), (2.6) where Cof(A) = (det A)A −T . For an NLC solid characterized by the strain-energy density given by (2.1), the stress tensors when the director is 'free' to rotate relative to the elastic matrix and when the nematic director is 'frozen' and satisfies condition (2.4) are presented in appendix A (see also [65]).
In the next sections, we obtain two stretch moduli under small elastic uniaxial tension and three nonlinear shear moduli under elastic simple shear deformations, respectively, which combine elastic and nematic effects. In our calculations, the nematic director is 'frozen', but the case where the director is 'free' can be treated similarly, provided that the elastic deformation is small. In particular, when the director is free to rotate and a tensile force is applied perpendicular to the director, experimental results show that there is a range of strains, up to 10% (e.g. [78,79]), before the director rotates in response to the applied force. However, the local nematic order might be altered. We further note that, in finite elasticity, the strain-energy density W(A) and the stress relationship (2.5) characterize an isotropic material. Yet, the response of a nematic solid also depends on the director orientation. For instance, we show that the ideal model exhibits an anisotropic response under stretch if the natural deformation is finite. There are, in fact, two possible contributions to aeolotropy. First, in the ideal case, the elastic energy of the system is isotropic but the full energy depends on the orientation of the nematic field as well [59]. Second, there is a component of the energy that depends directly on the nematic order through the socalled Q-tensor [23]. In this case, we need to extend the discussion by including this nematic energy density as developed later in the paper.

Stretch moduli
The stretch modulus of a homogeneous isotropic elastic material is obtained under uniaxial tension with the gradient tensor in a Cartesian system of coordinates taking the form [57] A = where λ > 1 is the stretch ratio in the direction of the applied tensile force. Assuming that the only non-zero component of the associated first Piola-Kirchhoff stress, given by (2.6), is in the tensile direction, i.e. P 22 > 0, it follows that The Young modulus at small strain is then defined as [57] E = lim where P 22 is the tensile first Piola-Kirchhoff stress given by (3.2), λ − 1 is the corresponding tensile strain (for different definitions of an elastic strain, see [57]) and μ = lim λ→1 (β 1 − β −1 ) is the shear modulus at small strain.
To derive stretch moduli for the nematic material, we apply the tensile force in a direction that is either parallel or perpendicular to the reference nematic director [45]. In each case, we assume an overall deformation where the elastic component is a uniaxial tension with the deformation tensor given by (3.1). Then, we calculate the stretch moduli by taking the ratio between the first Piola-Kirchhoff tensile stress and the associated strain in the limit of small elastic tensile strain, while the natural deformation remains finite.

(a) Nematic director parallel or perpendicular to the tensile direction
When the tensile force is acting in the second Cartesian direction and the reference director is parallel to the tensile force, the nematic director in the current configuration and the associated spontaneous deformation tensor take the following form, respectively (figure 1a): When the reference director is perpendicular to the tensile force, the nematic director in the current configuration and the associated spontaneous deformation tensor are, respectively (figure 2a), Assuming that the elastic deformation tensor A is of the form given by (3.1), for the deformation gradient F and the associated first Piola-Kirchhoff stress, the principal components in the direction of the tensile load are, respectively, and F where P 22 is the first Piola-Kirchhoff stress given by (3.2). These are illustrated in figures 1b and 2b, respectively, for μ = 1 and ν = 1/2. We define the following stretch moduli for the nematic material under the above two deformations, respectively: and where E = 3μ is Young's modulus defined by (3.3). Therefore, the stretch moduli given by (3.8) and (3.9) satisfy the following relation: We infer that: if r > 1, then

Shear moduli
To obtain the shear modulus of a homogeneous isotropic elastic material, a standard deformation is the simple shear, with the following gradient tensor in a Cartesian system of coordinates [57]: where k > 0 is the shear parameter. The non-zero components of the associated Cauchy stress tensor, given by (2.5), are and the non-zero components of the corresponding first Piola-Kirchhoff stress tensor, given by (2.6), are The shear modulus at small shear is then obtained as follows [57]: To obtain suitable shear moduli for the nematic solid, we assume an overall deformation where the elastic component is a simple shear with the deformation tensor given by (4.1). In general, one cannot easily separate the individual contributions to the deformation gradient given by (2.2). However, for the three shear configurations adopted here, it is possible to write these components explicitly. Then, in each case, we can calculate the shear modulus by taking the ratio between the first Piola-Kirchhoff shear stress and the corresponding shear strain in the limit of small elastic shearing strains, while the natural deformation remains finite.
where P 12 is the shear component of the elastic first Piola-Kirchhoff stress given by (4.3). These are represented in figure 3b, for μ = 1 and ν = 1/2. We now define the shear modulus for the nematic material at small shear as follows: For sufficiently small values of k, the corresponding shear strain and first Piola-Kirchhoff shear stress of the nematic material take the form with the components P 12 and P 22 of the elastic first Piola-Kirchhoff stress given by (4.3). These are illustrated in figure 4b, for μ = 1 and ν = 1/2. In this case, the associated shear modulus at small shear is equal to By (2.2) and (A4), the corresponding shear strain and first Piola-Kirchhoff shear stress for the nematic solid are, respectively, F where P 12 is the shear component of the elastic first Piola-Kirchhoff stress given by (4.3). These are represented in figure 5b, for μ = 1 and ν = 1/2. We now define the shear modulus for the nematic material at small shear as follows: A comparison of the shear moduli given by (4.7), (4.10) and (4.13) implies We conclude that: if r > 1, then

Contribution of the nematic free energy
The results of § §3 and 4 imply that, for an ideal nematic material, the effective shear and stretch moduli respect the same inequalities (e.g. if r > 1, then we have both E (1) < E (2) and μ (1) < μ (2) < μ (3) , and so on). However, numerous experimental results have demonstrated that there are significant differences between the behaviour of real nematic solids and that of ideal nematic materials analysed in the previous sections [38,42,45]. In particular, it was found that the ideal behaviour of the stretch and shear moduli ratio does not match experimental results. We therefore extend the strain-energy function by taking into account the nematic free-energy density, in addition to the isotropic strain-energy density W (nc) = W (nc) (F, n), given by (2.1), as where W n is equal to [38,43] (see also [19], ch. 2) Btr(QQQ) + 2 9 Ctr(QQQQ) + · · · , (5.2) with Q the order parameter tensor. This macroscopic tensor parameter is used to describe the orientational order in nematic liquid crystals [23] (see also [80]). For incompressible nematic elastomers subjected to uniaxial stretches, the contribution given by (5.2) to the total strain-energy density described by (5.1) was originally analysed in [38]. Following a similar approach, we restrict our attention to the one-term strain-energy density of the form (2.1) (see also [58]), where 'tr' denotes the trace operator and μ > 0 is the elastic shear modulus at small strain. This strain-energy density can be expressed equivalently as follows [14]:  [19], §2.2), we consider the nematic energy described by (5.2) and approximate it by [38] CQ 4 + · · · + 1 6 (A + 2BQ + CQ 2 + · · · )b 2 + · · · . (5.5) We therefore approximate the total strain-energy density defined by (5.1) as where W (nc) = W (nc) (F, n) and W n = W n (Q, b) are given by (5.4) and (5.5), respectively. For the components of G 2 , the first-order approximation of the Taylor expansions about the backbone order parameters in the initial state (Q, b) = (Q 0 , 0) are [38] where l ,Q = ∂l 1 /∂Q and l ⊥,Q = ∂l 2 /∂Q = ∂l 3 /∂Q denote the first derivatives of l 1 and l 2 (and also l 3 ) with respect to Q at b = 0, respectively, and l 2,b = ∂l 2 /∂b and l 3,b = ∂l 3 /∂b are the first derivatives of l 2 and l 3 , respectively, with respect to b at Q = Q 0 . For a nematic solid with the strain-energy density given by (5.6), we derive the stretch and shear moduli in a direction parallel or perpendicular to the nematic director, as follows.
This contributes to the small variation in total strain-energy density defined by (5.6) where W (lce) The above quadratic function in δQ has a minimum of Denoting E = 3μ, the stretch modulus in a direction parallel to the nematic director is then

) Extension perpendicular to the director
When the director is aligned in the first direction and F = diag(λ 1 , λ 2 , λ 3 ), where λ 1 = λ 3 = 1/ √ λ and λ 2 = λ > 1, assuming λ = 1 + ε and λ −1 = 1 − ε, with ε the infinitesimal strain, the elastic strainenergy density given by (5.4) is approximated by This contributes to the total strain-energy variation QQ is defined as before and W (lce) bb = ∂ 2 W (lce) /∂b 2 | b=0 . Minimizing the above function with respect to δQ and b gives Denoting E = 3μ, the stretch modulus in a direction perpendicular to the nematic director is

(c) Shear parallel to the director
We recall that a simple shear deformation of a hyperelastic material is equivalent to a biaxial stretch ('pure shear' [81]) in the principal directions [82]. We assume that the director is aligned uniformly in the first (or second direction), and F = R diag(λ 1 , λ 2 , λ 3 ), where λ 1 = λ > 1, λ 2 = 1/λ and λ 3 = 1, and R is the rotation by an angle of π /4 in the plane formed by the first two directions. Taking λ = 1 + ε and λ −1 = 1 − ε, with ε the infinitesimal strain, the elastic strain-energy density is approximated as follows: This contributes to the total strain-energy variation bb defined as before. We minimize this function with respect to δQ and b, The linear shear modulus in a direction parallel to the nematic director is then

(d) Shear perpendicular to the director
Next, assuming that the director is aligned in the first direction, we set F = diag(λ 1 , λ 2 , λ 3 ), where λ 1 = 1, λ 2 = λ > 1 and λ 3 = 1/λ. Taking λ = 1 + ε and λ −1 = 1 − ε, with ε the infinitesimal strain, the elastic strain-energy density is approximated by This contributes to the small variation in total strain-energy density where W (lce) bb is defined as before. The above quadratic function in b has a minimum value of The linear shear modulus in a direction perpendicular to the nematic director is equal to

Molecular dynamics simulation
We performed molecular dynamics (MD) simulations to synthesize LCEs, form nematic LCEs and characterize their response under stretch and shear deformations. Given its accuracy and efficiency for modelling mesogen-polymer systems, we used a hybrid force field, including Gay-Berne coarse-graining potentials for mesogen-mesogen interaction and Lennard-Jones (LJ) potentials for united atoms of hydrocarbon groups, CH x , in polymer chains. Our computer simulations [83] can serve as a virtual experiment to observe the macroscopic mechanical behaviour, which can then be compared directly with the continuum theory, and to calculate the physical quantities arising from atomistic movement and gain a mechanistic understanding.
In the MD simulations of LCEs, the potential energy of the whole system includes contributions from bond stretch, angle bending, dihedral rotation, non-bonded LJ interaction between united atoms (a-a), anisotropic non-bonded Gay-Berne interaction between mesogens (m-m) and the extended Gay-Berne interaction between united atoms and mesogens (a-m): and of the non-branched X-CH 2 -X and θ (2) the equilibrium angle of branched X-CH-X; E dihedral is a multiple-term harmonic dihedral style, with φ i the ith torsion angle, C n the non-branched X-CH 2 -CH 2 -X and C (2) n the branched X-CH 2 -CH-X; E (a-a) represents the LJ potential between non-bonded united atoms CH x , with r ij the distance between two united atoms and a i , c i the factorized energy parameters for CH x ; E m−m is the Gay-Berne potential between non-bonded mesogens, with U r the shifted distance-dependent interaction, η, χ the orientation-dependent energy, A i the transformation matrix for mesogen i, r ij the centre-to-centre vector between the ith and jth mesogens and all the rest of the parameters specified as constants in table 1; and E a−m denotes the extended Gay-Berne potential between a non-bonded mesogen and CH x following the standard mixing rule [86]. The cut-off distance is 9.8 Å for the LJ potential and 16.8 Å for the Gay-Berne potential. For the detailed explanation and explicit form of the Gay-Berne potential, readers should refer to [87,88]. First, 64 molecules of side-chain liquid crystal polymers were created, where every molecule has a backbone of 100 hydrocarbon monomers and 50 side chains attaching to the backbone in a syndiotactic way. Among the 50 side chains for each molecule, 20% were randomly selected to be attached with cross-linking sites, and the rest were attached with mesogens, as shown in figure 6a. Thus, every LC molecule has a different configuration. The 64 molecules were mixed by heating at 800 K for 50 ns. Then the system was quenched down to 500 K during 10 ns. Cross-linking of the first step was established while equilibrating the system at 500 K. A weakly cross-linked isotropic LCE was constructed, as shown in figure 6b. The system was found to have its isotropic-nematic phase-transition temperature below 490 K, consistent with other MD studies of similar LCE systems [89]. To form nematic LCEs, the whole system was quenched down to 450 K with an external field U efield i = −1.0 · P 2 (cos θ i ), where P 2 (x) = (3x 2 − 1)/2 is the second Legendre polynomial and θ i is the angle between the long axis of the ith mesogen and the external field. The external field has been experimentally and computationally proved to accelerate the formation of the nematic phase [90,91]. The quenching stage from 500 K to 450 K lasted 10 ns and the equilibrium stage at 450 K lasted 20 ns, for both of which the external field was applied along the z-direction. The external field was removed while the LCE system was equilibrated at 450 K for another 17 ns and the nematic phase was found to be stable. The nematic order S 2 = P 2 (cosθ i ) is shown in figure 6b during the quenching and equilibrium state. Subsequently, the second stage of cross-linking was performed within the nematic LCEs. The nematic system was then quenched down to 300 K and the nematic order was found to be well maintained. Throughout the whole process of constructing        [92]. However, the MD simulations of LCEs using Gay-Berne potentials are known to show thermal expansion effects [93]. To eliminate these non-physical effects, the isotropic LCE at 500 K was quickly quenched down to 300 K and the three dimensions for the isotropic phase were measured as l iso x = 120.3505 Å, l iso y = 118.7118 Å and l iso z = 117.2290 Å. The respective stretch ratios were then calculated as λ x = 0.8280, λ y = 0.7426 and  λ z = 1.5844. Ideally, the biaxial Gay-Berne potential for mesogen-mesogen interaction should yield the same stretch ratio along the X-and Y-directions. Here, however, the small size of the LCE system in the MD simulations gives rise to slightly different λ x and λ y . To eliminate the size effect on the spontaneous deformation, these stretch ratios were averaged to λ x−y = 0.7853. Then the anisotropy ratio can be estimated as r ≈ λ z /λ x−y = 2.0176, Poisson's ratio is ν ≈ 0.5 and a = r 2 ≈ 4.0706. After obtaining the nematic LCE at 300 K, we subjected the system to either stretch or shear deformations as described in § §3 and 4, respectively. For each deformation, we can calculate the evolution of the first Piola-Kirchhoff stresses with respect to the strains. We first calculate the Cauchy stress tensor in LAMMPS based on T ij = Σ N k m k v ki v kj /V + Σ N k r ki f kj /V, where m k represents the mass of the kth atom, v ki denotes the velocity of the kth atom in the ith dimension, N is the total number of atoms, V is the total volume of the system, N is the number of atoms pairs and r ki and f kj are the positions and forces of atom k along the ith and jth dimensions, respectively [94]. Next, given the deformation imposed, we construct the elastic deformation tensor A and, subsequently, the corresponding first Piola-Kirchhoff stress tensor P, following (2.6). Then, imposing the spontaneous deformation tensor G, or simply following (3.6), (3.7), (4.6), (4.9) or (4.12), we obtain the final P μ (2) μ (1) = 1.9450 ≈ r, μ (3) μ (2) = 3.4460 > r. (6.2) These relations are in agreement with (5.24) when E /E ⊥ > 1 and μ /μ ⊥ < 1 (see also fig. 7 of [38], fig. 7 of [45], fig. 6 of [48] and fig. 3 of [51]). Our theoretical results are independent of the manufacturing history. Experimentally, the mechanical properties of main-chain and sidechain LCEs are compared in [51]. In future work, it would be interesting to carry out separate simulations for these different cases to compare their behaviour computationally as well.

Conclusion
We studied theoretically and computationally the mechanical behaviour of nematic LCEs under different stretch and shear deformations. Theoretically, we first examined ideal nematic elastomers characterized by a homogeneous isotropic elastic strain-energy density, then also phenomenological models incorporating an additional nematic energy. We showed that these cases are qualitatively different, and that the generalized model does not necessarily order stretch moduli in the same way as the shear moduli. We also performed molecular dynamics simulations to analyse numerically the responses of simulated systems under a similar set of deformations, and found different mechanical responses in different directions. Therefore, the trifecta of experiments, computations and theory leads us to conclude that the contribution of the nematic free energy cannot be ignored, even in small deformations, and that LCEs are best understood as aeolotropic materials. When Frank effects also play an important role, they need to be taken into account as well. However, the deviation from isotropy is well captured by including the nematic energy, and this constitutes an important step in the constitutive modelling of liquid crystalline solids. is the first Piola-Kirchhoff stress tensor for the nematic solid material, with P the first Piola-Kirchhoff stress given by (2.6). 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 (2.5), J = det F, p (nc) is the Lagrange multiplier for the volume constraint J = 1 and q is the Lagrange multiplier for the constraint (2.4). Then is the first Piola-Kirchhoff stress tensor for the nematic solid.