The Anelastic Ericksen Problem: Universal Deformations and Universal Eigenstrains in Incompressible Nonlinear Anelasticity

Ericksen’s problem consists of determining all equilibrium deformations that can be sustained solely by the application of boundary tractions for an arbitrary incompressible isotropic hyperelastic material whose stress-free configuration is geometrically flat. We generalize this by first, using a geometric formulation of this problem to show that all the known universal solutions are symmetric with respect to Lie subgroups of the special Euclidean group. Second, we extend this problem to its anelastic version, where the stress-free configuration of the body is a Riemannian manifold. Physically, this situation corresponds to the case where nontrivial finite eigenstrains are present. We characterize explicitly the universal eigenstrains that share the symmetries present in the classical problem, and show that in the presence of eigenstrains, the six known classical families of universal solutions merge into three distinct anelastic families, distinguished by their particular symmetry group. Some generic solutions of these families correspond to well-known cases of anelastic eigenstrains. Additionally, we show that some of these families possess a branch of anomalous solutions, and demonstrate the unique features of these solutions and the equilibrium stress they generate.


Introduction
Universal deformations in nonlinear elasticity are deformations that exist for all members of a particular class of materials in the absence of body forces. Given any member of a particular class of materials, any universal deformation for that class can be maintained by the application of surface tractions alone. For instance in unconstrained isotropic elastic materials, only homogeneous deformations are universal. However, adding material constraints, i.e., restricting the class under consideration, expands the set of universal solutions. In particular, under the imposition of incompressibility, there are five known families of universal deformations in addition to the universal homogeneous deformations, which is now restricted to isochoric homogeneous deformations in keeping with the material constraint.
The process of obtaining and classifying all universal solutions is a highly nontrivial task. This line of research originates in the seminal work of Jerald Ericksen. In 1954, he made the first systematic attempt to classify all universal deformations in isotropic incompressible elasticity [Ericksen, 1954]. His work revealed four families of universal solutions in addition to homogeneous solutions. In 1955, he completely solved the analogous problem for unconstrained elastic bodies, proving that the only compressible universal solutions are homogeneous deformations [Ericksen, 1955]. In the case of incompressible elasticity, another family of universal solutions was then discovered by Singh and Pipkin [1965], with a special case of this family discovered by Klingbeil and Shield [1966]. Additionally, Fosdick [1966] noted that a different special case of this deformation represented a universal solution with constant invariants, a special case not addressed by Ericksen's initial work. Further contributions and specializations of this problem were made by a number of authors [Fosdick and Schuler, 1969, Kafadar, 1972, Knowles, 1979, Marris, 1982, Martin and Carlson, 1976 and the current conjecture is that no other solution to Ericksen's problem exist but a proof of it remains an outstanding open problem of rational mechanics [Antman, 1995].
Here, we are interested in generalizing Ericksen's problem to nonlinear anelasticity. In anelasticity, we consider geometric deformations combining both elastic deformations and an additional anelastic component to the deformation, i.e., one that does not contribute to the strain energy. Such theories are known to be relevant in many situations that generalize classical nonlinear elasticity [Epstein and Maugin, 1996a] such as thermal effects [Stojanović et al., 1964, Ozakin andYavari, 2010], plastic flows [Kondo, 1949], dislocations and defects [Nye, 1953, Bilby et al., 1955, Yavari and Goriely, 2012, growth and remodeling [Goriely, 2017, Naumov, 1994, Rodriguez et al., 1994, Takamizawa, 1991, Yavari, 2010, and swelling Tsai, 2006, 2005]. Such processes are characterized by the presence of eigenstrains [Mura, 2013] that do not produce corresponding stresses or, equivalently, by changing the intrinsic geometry of the reference configuration from Euclidean to Riemannian. These eigenstrains are generally incompatible, and therefore further elastic strains are typically required to embed bodies with nontrivial eigenstrains into Euclidean space, with the resulting self-equilibrating elastic stresses generated by this strain referred to as residual stresses.
The first step in this research program was to generalize Ericksen's theorem for compressible isotropic materials to anelastic deformations. By using a geometric formulation, it was proved that only covariantly homogeneous deformations are universal [Yavari and Goriely, 2016]. The second step, considered here, is to extend the current classification of isotropic incompressible nonlinear elasticity to isotropic incompressible nonlinear anelasticity. This general problem is more involved than the classical Ericksen problem as we have to determine simultaneously both the anelastic and elastic components that render the solutions universal.
While anelastic deformations can be modeled through a multiplicative decomposition [Sadik and Yavari, 2017], it is equivalent but more appropriate in our context to model them as a stress-free evolution into of a general Riemannian manifold, the material manifold, via some anelastic process [Noll, 1967, Epstein and Maugin, 1996b, Yavari and Goriely, 2012, Goriely, 2017. This non-Euclidean material configuration contains all the information of the anelastic processes. It can then be mapped by an elastic deformation into the current Euclidean configuration. Only the strain induced by the elastic component appears in the strain-energy formulation, which models the notion that the anelastic deformation changes the relaxed geometry of the material, and that only the further elastic deformation stores energy by straining the material. In looking for universal solutions, both the Riemannian metric of the material manifold and the elastic deformation are unknowns to be determined.
In this paper, we extend the currently known families of incompressible universal solutions to the anelastic setting in an appropriately symmetric way, which will be precisely defined shortly. In the process of doing this, we discover that under these symmetry conditions, all families exhibit a branch of generic universal solutions that contain arbitrary functions as parameters, but some families also contain anomalous universal solution branches outside of these, whose form is fixed up to a finite number of constants. Additionally, we find that the six classical families of universal solutions merge into three distinct anelastic families, characterized by their respective symmetry groups. The Cauchy deformation tensors of the generic branches of these families can all be expressed in the forms presented here in Table  1, with the standard forms of the anomalous branches presented in section 8, as they are more involved. In §2, we provide an overview of the geometric features of anelasticity and provide the relevant governing equations of elasticity in this context. The known universal solution families are summarized in §3. Then, in §4 we compute the appropriate symmetry group for each family and impose this symmetry on the metric tensor field a priori. We then formulate the problem of extending the universal families to the anelastic setting and outline the techniques used in our analysis in §5. In §6 we derive the form of the generic solutions for each family and obtain the constant invariant conditions that are necessary for anomalous solutions to exist. In §7 we present the general form of the anomalous solutions, relegating the explicit calculations to the appendix B. In §8, we examine the overlap of these families of solutions. Finally, we present some visualizations of the Riemannian geometry of strains and stresses induced by these anomalous solutions in §9, and summarize our results in §10. We have provided two other appendices, Appendix A containing summaries of the algebraic and group theoretic tools we employ in our analysis, and Appendix C, detailing particular features of the Lie algebra se(3), which plays a key role in our analysis.

Nonlinear Elasticity and Anelasticity
In the geometric formulation of nonlinear elasticity, we define the ambient space to be a Riemannian manifold (S, m), where m is a fixed background metric. Since we are interested in universal deformations that take place in Euclidean space we assume that the ambient space is Euclidean. Then, a body is defined as a Riemannian manifold (B, M). We define a motion as an isotopy ϕ ϕ : B × R → C ⊂ S, parameterized by time t that gives a homeomorphism at each time t between the reference configuration B and the physical configuration at time t defined by ϕ(B, t) (see Figure 1): We use coordinate charts {X A } and {x a } for B and C, respectively. We utilize uppercase Latin letters to denote most quantities  Table 1: Standard forms of the Cauchy deformation tensors for the generic branches of the three symmetric anelastic incompressible isotropic universal deformation families. and indices in the reference configuration, and lowercase Latin letters to denote most quantities and indices in the physical configuration. The homeomorphism at time t, ϕ(B, t) can be interpreted physically as determining the position of a small piece of material at time t, given its position in B, which is interpreted as an initial position, though the important feature of the reference configuration is that it specifies the relaxed geometry of the body; it is not necessarily the initial configuration of the body. While in elasticity, (B, M) is Euclidean, for an anelastic system, (B, M) may not be Euclidean and, in such case, ϕ 0 ≡ ϕ(−, 0) is not the identity map. Since ϕ 0 (B) and ϕ(B, t) correspond to physical configurations, we model them as subsets of Euclidean space, and hence, we identify positions with vectors.
Since we are interested in equilibrium states, we restrict our attention to a finite time t * > 0 and suppress the explicit time-dependence so that we define, with a slight abuse of notation, ϕ(B) = ϕ(B, t * ).

Kinematics
The local properties of deformations are encapsulated in the derivative of the map ϕ that we explore next. The tangent space of B at X ∈ B is denoted by T X B. The tangent space of the corresponding point x = ϕ(X) in the ambient space is denoted by T ϕ(X) C. The deformation gradient F is the derivative map of ϕ: F(X) = T ϕ(X) : T X B → T ϕ(X) C.
With respect to the coordinate charts {X A } and {x a }, F it is defined as follows The adjoint of F is In components, (F T ) A a = M AB m ab F b B . Note that the adjoint explicitly depends on the metrics M and m. The right Cauchy-Green deformation tensor is defined as It is straightforward to see that C = ϕ * m, which has components C AB = F a A F b B m ab [Marsden and Hughes, 1983]. The Jacobian relates the material and spatial Riemannian volume elements as dv(x, m) = JdV (X, M), where J is given by

Equilibrium Equations
The balance of linear momentum in the absence of body and inertial forces in terms of the Cauchy stress tensor reads div σ = 0 , where div is the spatial divergence operator, defined in components as and γ a bc is the Christoffel symbol of the Levi-Civita connection ∇ m of the metric m in the coordinate chart {x a } , defined as ∇ m ∂ b ∂ c = γ a bc ∂ a . The local form of the balance of angular momentum reads σ ab = σ ba , i.e., the Cauchy stress is symmetric.
The first Piola-Kirchhoff stress is defined as P = JσF − , and in components, P aA = Jσ ab (F −1 ) A b . F is the deal of F, and is defined as F = ∂ϕ a ∂X A dX A ⊗ ∂ ∂x a . The balance of linear momentum in terms of P reads Div P = 0. In components where Γ A BC is the Christoffel symbol of the Levi-Civita connection ∇ M of the metric M in the coordinate chart {X A } , defined as ∇ M ∂ B ∂ C = Γ A BC ∂ A . The balance of angular momentum in terms of P reads PF T = FP T .

Constitutive equations
In this paper we restrict ourselves to bodies that are isotropic in the absence of eigenstrains. We also assume that the elastic deformations are incompressible. The left Cauchy-Green stretch tensor is defined as B = ϕ −1 * (m −1 ) and has components B AB = (F −1 ) A a (F −1 ) B b m ab , where m ab are components of m −1 . The spatial analogues of C and B are defined as where b is called the Finger deformation tensor. The two tensors C and b have the same principal invariants, which are denoted by I 1 , I 2 , and I 3 [Ogden, 1984]. In the case of an isotropic solid the energy function W depends only on I 1 , I 2 , and I 3 . If the material is incompressible (I 3 = 1), W = W (I 1 , I 2 ), and the Cauchy stress has the following representation [Simo and Marsden, 1984] where p is a Lagrange multiplier that is associated with the incompressibility constraint J = 1.
1 As p is a scalar field to be determined the stress representation can equivalently be written as

Ericksen's problem
The classical elastic Ericksen problem is stated as follows: Determine all equilibrium deformations ϕ : B → C with B, C ⊂ E 3 that can be sustained by an arbitrary incompressible isotropic hyperelastic solid with suitable boundary tractions. The emphasis of the classical problem is that both configurations are Euclidean. Here, we consider a generalized version of this problem applicable to anelastic systems. The anelastic Ericksen problem relaxes the requirement that B ⊂ E 3 and is stated as follows: Determine all Riemannian manifolds (B, M), and all maps ϕ : B → C with C ⊂ E 3 that can be sustained by an arbitrary incompressible isotropic hyperelastic solid with suitable boundary tractions.
These problems are treated locally in the sense that the equilibrium equations are locally satisfied by these deformations for arbitrary incompressible isotropic hyperelastic materials. We do not address non-local problems such as self-intersection or self-contact beyond requiring that our solutions be locally homeomorphisms, which is guaranteed by the condition det F > 0. This condition ensures that over some domain our solution is an embedding, rather than an immersion.

The Known Universal Deformations
We begin with the known families of incompressible universal solutions in the absence of eigenstrains. We merely present them and direct the reader to the original papers for their derivation [Ericksen, 1955, 1954, Klingbeil and Shield, 1966, Singh and Pipkin, 1965. The corresponding deformation gradients are derived explicitly in Goriely [2017]. The emphasis and novelty here is in the particular type of symmetry associated to each family as they will play a key role in the generalization of the problem to anelastic systems. Expressed in terms of the standard Cartesian coordinates {x, y, z}, cylindrical polar coordinates {r, θ, z}, or spherical polar coordinates {r, θ, φ} (letting capital letters denote reference configuration coordinates, and lower case letters denote current configuration coordinates), we have the following six universal families.
Family 0: Homogeneous Deformations. Using the Cartesian coordinates {x a } = {x, y, z} and {X A } = {X, Y, Z}, this family is most compactly expressed as where F a A is a constant tensor with det F a A = 1, and c a is a constant vector. A deformation of this type is depicted in Figure 2. The form of equation (13) immediately reveals that the deformation gradient is F a A , as evidenced by the induced tangent map dx a = F a A dX A . Since the deformation gradient F is spatially constant, F a A (X A ) = F a A , the transformation leaves the deformation gradient unchanged, and hence, C AB remains unchanged. In terms of symmetry group, we notice that the action X A → X A + C A is precisely the action of T(3) ⊂ SE(3) on E 3 . Family 1: Bending, Stretching, and Shearing of a Rectangular Block. When expressed using cylindrical polar coordinates {r, θ, z} and Cartesian coordinates {X, Y, Z} in the current and reference configurations, respectively, universal deformations in this family take the form though the parameters E and F only correspond to rigid motions, and hence, can be safely disregarded. These parameters generate rotation around and translation along the r = 0 axis, as seen in Figure 3. We can compute the deformation gradient, which when expressed on mixed orthonormal frames, takes the form or in terms of a mixed coordinate bases, we have the components Working in terms of the mixed bases will often be advantageous when computing symmetries and the components of arbitrary metric tensors, as these bases are automatically induced by the coordinate map. Additionally, we demand AB = 0 to ensure the deformation is invertible; the incompressibility condition is automatically satisfied. We compute C AB as and note that this remains unchanged under the transformation This is precisely the action of T(2) ⊂ SE(3) on E 3 . Figure 3: The stretching, shearing, and bending of a rectangular block around a cylinder.
Family 2: Straightening, Stretching, and Shearing of an Annular Wedge. Deformations in this family are most naturally expressed using Cartesian and cylindrical polar coordinates {x, y, z} and {R, Θ, Z} in the current and reference configurations respectively, and take the form An example of one of these deformations is depicted in Figure 4. The corresponding deformation gradient is or, in terms of the induced coordinate bases We demand AB = 0 to ensure the deformation is invertible, and as in the previous case, the incompressibility condition is automatically satisfied. Thus The transformation leaves these components unchanged. This is the action of SO(2) × T(1) ⊂ SE(3) on E 3 . Family 3: Torsion, Extension, and Shearing of an Annular Wedge. When expressed in cylindrical polar coordinates, deformations in this family take the form and an example of a deformation from this family is depicted in Figure 5. The deformation gradient can be naturally expressed on orthonormal cylindrical polar bases as or equivalently on the coordinate bases, it has components We have the incompressibility condition A (CF − DE) = 1, which also ensures invertibility. Thus, C AB is written as and notice that C AB does not depend on Θ or Z. Family 4: Inflation/Eversion of a Sphere. In spherical polar coordinates {r, θ, φ} and {R, Θ, Φ}, deformations in this family take the form An example of one of these deformations is depicted in Figure 6. The deformation gradient on orthonormal bases reads or, in terms of the coordinate bases, we have the components Incompressibility and invertibility are trivially satisfied. Thus We can then represent this tensor on an orthonormal spherical basis (using the standard Euclidean metric) as where I is the identity tensor, and E R = X |X| . This obeys the symmetry transformation which is symmetry under the prolonged action of SO(3) ⊂ SE(3) on TE 3 ⊗ T * E 3 . Family 5: Inflation, Bending, Extension, and Azimuthal Shearing of an Annular Wedge. When expressed in cylindrical polar coordinates {r, θ, z} and {R, Θ, Z}, deformations in this family take the form r = AR, θ = B log R + CΘ + D, z = EZ + F.
An example of one of these deformations is presented in Figure 7. The deformation gradient expressed on orthonormal bases is written as or, on the coordinate bases In order to satisfy incompressibility, we have A 2 CE = 1. This family is peculiar, as the stretch generated by the other inhomogeneous families has an eigenvector along the direction of inhomogeneity, but this one does not. Additionally, the invariants of b for this family are spatially constant. When we generalize these deformations to include an anelastic component, we have to change the incompressibility condition from det F = 1 to to reflect the fact that we are only constraining the elastic component of the deformation to be isochoric. It is easier however, to consider square of this equation in the form This is because when we move to the anelastic setting, the stretch b is a more natural object to work with, since it captures geometric data about the material manifold, but itself lives in Euclidean space. The right Cauchy-Green strain reads which is invariant under changes in Z or Θ.
3.1 Summary of the symmetry groups.
The symmetries we have calculated are all generated by the usual action of Lie subgroups of the special Euclidean group on the reference configuration. For each family, there is some continuous family of rotations and/or translations, which once prolonged, leaves the right Cauchy-Green stretch tensor field unchanged. In a similar manner, we can compute symmetries of the left Cauchy-Green stretch tensor field, which also happen to be Lie subgroups of the special Euclidean group but, acting now on the current configuration. These groups are summarized in Table 2, expressed in terms of the group of n-dimensional rotations SO(n) and the group of n-dimensional translations T(n). Interestingly enough, for each family, the symmetry group of C does not necessarily match the symmetry group of b, but the dimensions of these two groups do match, as their actions are related through the maps relating the coordinates in the two configurations. Additionally, when we impose the symmetries of C on the material manifold, families with 3-dimensional Lie symmetries automatically satisfy the equilibrium conditions, but those only containing a 2-dimensional Lie symmetry require additional constraints to satisfy equilibrium. This seems to suggest that the dimension of the symmetry group plays a role in constraining the material response, and sufficiently high dimensional symmetries are sufficient for guaranteeing equilibrium.
These symmetries can be summarized as a single key property, namely that for a given universal deformation and its associated deformation gradient F and right Cauchy-Green tensor, C = F T F = M −1 ϕ * m, we have that C = ϕ * m is invariant under the prolonged action of a Lie subgroup of the special Euclidean group acting on the reference configuration. We will use this key property when studying the anelastic Ericksen problem, and accordingly refer to the symmetries of C as 'universal symmetries'.

General construction
The previous section demonstrated a remarkable symmetry property of the known universal solutions in the absence of eigenstrain. We can use these symmetries to generalize Ericksen's problem to anelastic systems. The problem is then to find a suitable metric on the material manifold that preserves both the symmetry and the general functional form of the universal deformations. The eigenstrains and metric associated with this new metric are referred to as 'universal'. This can be achieved by the following construction: • First, in the absence of eigenstrains, the body is embedded in the ambient space with an induced metricM. The material manifold is the flat Riemaniann manifold (B,M) and the deformation is a map from this manifold to the ambient space.
• Second, in the presence of eigenstrains the natural configuration of the body is a Riemannian manifold (B, M), where M has non-vanishing curvature [Yavari and Goriely, 2013, 2015, Golgoon and Yavari, 2018. In this case, the deformation is a map from (B, M) to the ambient space (S, m).
• Third, we choose curvilinear coordinates {X A } on B and curvilinear coordinates {x a } on C ⊂ S. These coordinates are not necessarily associated with the metricsM and M. We know of the previously presented classes of universal deformations x a = ϕ(X A ) for (B,M). We fix these functional dependences ϕ on the coordinates and determine the the metrics M compatible with these solutions.
• Fourth, we pull back m under the deformation ϕ, and consider the three manifolds: (B,M), (B, M), and (B, ϕ * m). We have two candidates for determining the symmetry to apply to M:M and ϕ * m.
We use ϕ * m since it captures information about the deformation. We compute its symmetries, and demand M to have these same symmetries. Explicitly, since bothM and m are Euclidean metrics, both are invariant under the full action of SE(3) acting on their respective base spaces. By unknown quantity to be determined is the metric tensor field M of the material manifold, which dictates the stress-free geometry.
considering Euclidean symmetries of ϕ * m, we are identifying Euclidean symmetries in the current configuration that are mapped to Euclidean symmetries in the reference configuration when pulled back under the classical universal deformation in question.
• Fifth, we compute the deformation mapping (B, M) to (C, m), where now M is restricted by the derived symmetries, and compute the specific metrics M that generate universal eigenstrains.

Universal equilibrium equations
We fix the coordinate labels for the anelastic component of the local deformation; if X A are coordinates for the material manifold, we write the anelastic deformation in terms of coordinates as where δ Ā A is the Kronecker delta, and use the metric tensor field to reflect the change in geometry (see Figure 8). In other words, we convect the referential coordinates along with the anelastic motion onto the intermediate configuration. This allows us to trivially pull the material metric tensor back to the reference configuration using and treat MĀB as our unknown quantity. Since MĀB and CĀB live in the same space, we can immediately impose the symmetry of CĀB on MĀB , which naturally imposes a symmetry on M AB via (43). After determining the most general form of a metric tensor field obeying these symmetries, we can compute the general form of the elastic left Cauchy-Green tensor and its inverse in terms of the undetermined components of the metric tensor field. Both of these appear in the Cauchy stress of an incompressible isotropic elastic solid, which has the following representation in terms of W i = ∂W /∂I i , where I 1 and I 2 are the two non-trivial invariants of b, and p is the Lagrange multiplier corresponding to the incompressibility constraint. We seek equilibrium solutions, hence we must satisfy We wish to eliminate the pressure field from the analysis, so we take a second covariant derivative, lower the upper free index, and compute the antisymmetric part. This operation eliminates the pressure identically and yields the condition which must be satisfied for arbitrary choices of the strain energy function W . Because W is arbitrary, we can freely vary its partial derivatives independently, so in order to satisfy (47) for any W , we require each of the terms multiplying a partial derivative of W to vanish independently. This yields nine independent conditions that must be satisfied, corresponding to the nine independent mixed partial derivatives of W that appear: W 1 , W 2 , W 11 , W 12 , W 22 , W 111 , W 112 , W 122 , and W 222 . These universal equilibrium equations are Each of these nine equations is antisymmetric in the two free indices. Therefore, each equation has three independent components providing an overdetermined system of 27 scalar conditions that must be satisfied. Imposing the universal symmetries lead to two different situations depending on the symmetry: Case I: For two universal families (the homogeneous and spherical deformations) these equations are trivially satisfied and do not lead to any new conditions.
Case II: For the other families, these equations will either require particular off-diagonal components of the metric tensor to vanish 2 (Case IIa) or the invariants to be constant (Case IIb). For Case IIa the equilibrium equations are trivially satisfied. This is the so-called generic case for all the remaining families analysed in Chapter 7. Case IIb corresponds to the anomalous solutions. In this case the symmetry condition generates a set of ordinary differential equations constraining these components in terms of a single independent variable. In addition to satisfying these, we also must satisfy three algebraic constraints, namely that the invariants of b are spatially constant. This leaves us with an overdetermined system of four linear differential equations, one linear algebraic equation and two nonlinear algebraic equations for the six unknown components of the metric tensor field. We can integrate the differential equations, and use the linear algebraic equation to express the other two algebraic equations in terms of the following 15 ≤ n ≤ 18 variables: a single unknown component of the metric tensor, the remaining independent variable in space (e.g. radius), the integration constants introduced by our integration of the ordinary differential equations, the deformation parameters, and the two constant invariants. Therefore, the remaining algebraic conditions are polynomial equations in these n variables that are quadratic in the unknown component of the material metric tensor field. We compute the resultant (see appendix A.1) of these polynomials in this component, and demand it vanish for the two equations to have a common solution, since our metric tensor must simultaneously satisfy both. This resultant is itself a polynomial in the dependent variable that must vanish. Because we seek solutions that are universal over an open set, we can send each of the coefficients of the resultant to zero identically. This leaves us with a set of algebraic equations relating the integration constants, the deformation parameters, and the invariants that must be necessarily satisfied in order for these anomalous solutions to exist. These algebraic equations are solved in Chapter 8.

Symmetries of the Material Metric
Family 0: Homogeneous Deformations. Recall that the action X A → X A + C A is the action of T(3) ⊂ SE(3) on E 3 . We seek to impose this symmetry on the material metric tensor field, and hence demand Taking the derivative of this with respect to C D gives and hence, we consider constant metric tensor fields. It is however, more useful to express this condition in terms of the current configuration variables, which we can do via the chain rule Since F is invertible, this implies that M AB is constant when expressed in terms of the current configuration coordinates as well. Additionally, it is more useful to consider the inverse metric tensor field M AB , which also must be constant in order for the identity to hold.
Family 1: Bending, Stretching, and Shearing of a Rectangular Block. Recall, that the symmetry associated with this deformation is the action of T(2) ⊂ SE(3) on E 3 . We therefore require the same invariance for M AB : Taking the derivative with respect to C 1 and C 2 independently gives the conditions Therefore, we assume that the metric tensor is of the form M AB (X), which because of the form of the deformation, can be recast into the form M AB (r), and equivalently M AB (r).
Family 2: Straightening, Stretching, and Shearing of a Sector of a Cylinder. Here, we require that the metric is invariant under the action of SO(2) × T(1) ⊂ SE(3) on E 3 : and hence, by the same reasoning as before one finds Note that this does not imply ∂M ∂Θ = 0, since the basis vectors of C do change with Θ. Rather, the components M AB do not change with Θ when M is represented with respect to a cylindrical polar basis.
We can use the equation Hence, we write M AB (R) or equivalently M AB (x). This symmetry of M AB is precisely the action SO(2) × T(1) ⊂ SE(3) prolonged to TE 3 ⊗ TE 3 , with ⊗ denoting the tensor product bundle, the bundle formed by taking fiberwise tensor products.
Family 3: Inflation, Bending, Torsion, Extension, and Shearing of an Annular Wedge. As with Family 2, we demand that M AB be invariant under the transformation which renders the conditions So we consider metric tensor fields of the form M AB (R), or equivalently, M AB (r).
Family 4: Inflation/Eversion of a Sphere. Here, we demand the invariance of M(X) under the prolonged action of SO(3) ⊂ SE(3) on T * E 3 ⊗ T * E 3 . We then seek the most general positive-definite symmetric tensor field that satisfies Because M(X) is positive-definite and symmetric, we can represent it in spectral form on an orthonormal basis {u a } as We then consider the subgroup of rotations leaving X fixed. Under this one-parameter family, we have the symmetry condition This implies that (suppressing the X dependence) i.e., the rotated vector Q T u i lies in the same eigenspace as the eigenvector u i . For this to hold for all Q in this one-parameter family, the eigenspaces of M at the point X must be unchanged by these rotations, i.e., the swept vector Q T u i remains in the eigenspace. Generally, a rotating vector sweeps out a cone, which not being an affine space, cannot be the eigenspace of a linear operator, as depicted in Figure 9. However, there are two degenerate cases where cones become affine spaces, where the rotation axis and swept vector are coincident, and where they are perpendicular, which means that the eigenspaces of M at X must contain the axis of the rotation, and the plane orthogonal to it, as depicted in Figure 10, forcing M(X) to be of the form because for each X, the axis of rotation is E R (X). Imposing the more general symmetry condition (68), on this spectral form, we get the condition Figure 9: The set of orthonormal eigenspaces at a point X are affine spaces that must be preserved under all rotations fixing the origin and X, but these swept areas are generally not affine spaces, but rather cones.
This ultimately requires that m 1 and m 2 depend on X solely through its norm, R = |X|, since for any two points, X 1 = Rn 1 , and X 2 = Rn 2 , where n 1 and n 2 are unit vectors, we can construct orthogonal transformations such that n 2 = Qn 1 . Hence, because the functions m 1 and m 2 take on the same values whenever their arguments have the same norm, these functions must only depend on their argument through its norm, as shown in Figure 11. We then have which is the general form of the pullback of our material metric tensor. 3 Computing the components of M AB then gives Since m 1 and m 2 are arbitrary, and R is only a function of r, this can be rewritten as and equivalently Figure 11: As we rotate the point X, we see that the eigenspaces are unchanged when expressed on an orthonormal spherical basis, i.e., they rotate with X, hence their associated eigenvalues must be constant on concentric spherical surfaces.
These metrics are precisely the form considered by Ben Amar and Goriely [2005], and Yavari and Goriely [2013], though the first represents this tensor on an orthonormal spherical basis, and the second works with the components of the metric tensor rather than its inverse, as we have done.
Family 5: Inflation, Bending, Extension, and Azimuthal Shearing of an Annular Wedge. As in other cases, we have the following symmetry relations for M: so we ultimately consider inverse metric tensor fields of the form M AB (r).
In conclusion, we note that the application of the symmetry condition leads to a reduction of the independent variables to a single one (either the radial variable r in cylindrical or spherical coordinates, or x for the deformation of rectangular blocks).

Generic Universal Solutions
Having established the symmetry conditions on the material metric, we can express the universal equilibrium equations under these restrictions. For all families (Case I and IIa), these equations will have generic solutions, and in some cases, these solutions are the only ones satisfying the symmetry conditions. Other families also have anomalous solutions (Case IIb)) outside of these generic branches that occur only when the invariants of the tensor b, or equivalently C, are constant; these will be addressed in the next section. The nature of the anomalous solutions differs markedly from the generic solutions found here: generic solutions contain arbitrary functions as free parameters, while the form of the anomalous solutions is determined up to a finite number of undetermined constants. Additionally, for generic solutions, the eigenvectors of b are contained within or perpendicular to the span of the infinitesimal generators of the symmetry group, while for anomalous solutions, this alignment does not occur. While the invariants, and hence, their gradients could be calculated explicitly in terms of the unknown inverse metric components, it is easier to keep these functions unevaluated at the moment, because we will ultimately show that they must be constant for the anomalous solution to exist.
We then fix an orthonormal frame on these constructed intermediate configurations, and express the anelastic deformation required to generate these intermediate configurations in terms of these frames. We note that since we are in reality only determining the elastic component of these motions, the factor A corresponds to the universal motion, while the factor G corresponds to the stretch required to obtain the universal intermediate configurations from the classical reference configurations, which play no dynamic role. Hence after computing A = FG −1 , so long as the geometry of the intermediate configuration is retained, the factor G can be discarded, which corresponds to the fact that we can prepend an arbitrary compatible anelastic deformation onto our universal deformations. As this is supplementary to our main results, we will simply give G on some orthonormal anholonomic frame, leaving further computations to interested readers.
Family 0: Homogeneous Deformations. The deformation mapping written in Cartesian coordinates is given by (13), and the deformation gradient F a A is constant. We compute the left Cauchy-Green tensor as which is also constant, and hence its invariants are constant. The Cauchy stress takes the form and the equilibrium equations read Because m ab is invertible, this implies that p is constant, and the equilibrium equations are satisfied simply because of the assumed form of M AB . The only remaining condition is the incompressibility condition which, in the chosen coordinate systems, reads det b ab = det m ab = 1.
We can express this constraint as a condition on M AB , or as a condition on F a A . In reality, these conditions are equivalent, since we have and one can be freely transformed into the other by changing coordinates. However, for the other families, it will be easier to express this condition as a restriction on the inverse metric tensor, so for consistency we choose the invertible tensor F a A and then enforce which ensures that the volume form in the material manifold agrees with that in the current configuration. Because the material metric is constant, its Levi-Civita connection produces no curvature, and thus the material manifold is Euclidean. This is useful when using a multiplicative decomposition of the deformation gradient into elastic and anelastic factors F = AG, as we can choose a Cartesian frame {e α } in the material manifold and its corresponding coframe {ϑ α }, in which case the anelastic factor must satisfy G α A G β B δ αβ = M AB . Since the matrix of components M AB is positive definite and symmetric, we can take the matrix of components G α A to be its unique positive-definite symmetric square root, in which case we satisfy G α A G β B δ αβ = M AB . Alternatively, we may prescribe the anelastic factor in such a way that the induced material metric is valid. In this case, since M AB is constant, any constant invertible anelastic factor will furnish a valid material metric. The incompressibility constraint becomes det(A a α A b β δ αβ m bc ) = 1, which furnishes a differential equation constraining the volume in the current configuration to agree with that in the material manifold.
Family 1: Bending, Stretching, and Shearing of a Rectangular Block. The deformation for this family is given in (15) with the deformation gradient (16). We compute the quantity m a[l ∇ k] ∇ b σ ab = 0, which for this family, only has two independent nonzero components, and we take the coefficients of the partial derivatives of W to vanish independently. The W 111 coefficient of this equation gives the conditions ABrM 12 (r) The first equation implies that either M 12 (r) = 0, or I 1 is constant, because AB = 0 to ensure the invertibility of the deformation. If I 1 is not constant, we have M 12 (r) = 0. The second component then becomes M 13 (r) = 0. Therefore, if I 1 is not constant, we have that M 12 (r) = M 13 (r) = 0. If I 1 is constant, we can examine the W 122 component, which implies the conditions Therefore, I 2 is constant, or M 12 (r) = M 13 (r) = 0 and we have established that either M 12 (r) = M 13 (r) = 0, or all the invariants of b are constant, which is the condition for the anomalous solution.
Hence, in this section, we take M 12 (r) = M 13 (r) = 0 and we consider the constant invariant case in the next section. With M 12 (r) = M 13 (r) = 0, we have the equilibrium equation satisfied, i.e., all of its terms identically vanish. We only have to satisfy incompressibility. Demanding det b ab = det m ab , we have the condition The current form of M AB automatically captures incompressibility because we imposed a particular form of r(X). If we leave this unspecified, we can simply say that M AB is of the form and use the incompressibility constraint to determine r(X). This is equivalent to introducing a change in coordinates rescaling X. When using a multiplicative decomposition, F = AG, we can choose an orthonormal frame {e α } and its coframe {ϑ α } in the material manifold, and require G α A G β B δ αβ = M AB . Since M AB is block diagonal, positive definite, and symmetric, we can take G α A to be its unique positive definite symmetric square root. Because the components M AB are arbitrary functions of X, we can take G α A to be of the form which will yield a suitable M AB . Additionally, we can multiply G by an arbitrary local rotation Q yielding QG, which may be more useful depending on the particular problem. 4 This is equivalent to choosing a different orthonormal frame in the material manifold, which being non-Euclidean in general, does not possess a preferred orthonormal frame to begin with.
Family 2: Straightening, Stretching, and Shearing of a Sector of a Cylinder. Recall that any deformation in this family is given by (20) with deformation gradient (21). We compute the equilibrium condition m a[l ∇ k] ∇ b σ ab = 0, and to aid computations, we use the incompressibility constraint det b ab = det m ab by evaluating c ab as c ab det (b n m ) = c ab det b np m pm . The W 111 coefficient of the equilibrium equation has two independent components giving the condi- If I 1 is constant, we satisfy these equations, but if If I 1 is constant, we then consider the W 122 component of the equilibrium condition to obtain which again implies that either M 12 (x) = M 13 (x) = 0, or I 2 is constant. Therefore, unless the invariants of b are constant, we have M 12 (x) = M 13 (x) = 0. Setting these components to 0 satisfies equilibrium, and so we compute the incompressibility condition det b ab = det m ab . This becomes which implies that Bringing all of these together we have or in terms of referential coordinates, writing M AB (R) = M AB (x(R)), This is the generic solution, and we have set up the conditions for the anomalous solution, namely that the invariants of b must be constant. As before, we can introduce a coordinate rescaling, treating x as an unknown function of R, which allows the tensor M AB to take the form and turns the incompressibility constraint into a differential equation that can be integrated to determine x(R). If a multiplicative decomposition of F = AG is used, we can express G on an orthonormal frame in the form which guarantees that M AB = G α A G β B δ αβ is of the proper form. As before, an arbitrary local rotation Q can be imposed yielding the factor QG, where G is as above, and this new factorization will yield a material metric of the proper form.
Family 3: Inflation, Bending, Torsion, Extension, and Shearing of an Annular Wedge. This family of deformations can be written using cylindrical polar coordinates in both configurations as given in (25) with deformation gradient (26).
As before, we compute the W 111 coefficient of the equilibrium equation and obtain The matrix on the left-hand side is invertible, since its determinant, CF − DE, being a factor of the determinant of F, is nonzero to ensure invertibility. Therefore, we have either I 1 being constant, or both M 12 (r) and M 13 (r) must be 0. If I 1 is constant, we examine the W 122 coefficient of the equilibrium equation and obtain which as before implies that M 12 (r) = M 13 (r) = 0, or I 2 is constant. Therefore, to satisfy equilibrium, we must have all of the invariants of b being constant, or M 12 (r) = M 13 (r) = 0. The latter of these conditions is also sufficient to guarantee equilibrium. We only have to satisfy incompressibility, which amounts to the equation which we can do by setting This gives the generic solution or in referential variables, writing M AB (R) = M AB (r(R)), As in the other families, we can introduce a coordinate rescaling to express the material metric in the form which means that on some orthonormal frame, the anelastic factor of a multiplicative decomposition F = AG takes the form Doing this turns the incompressibility condition into a differential equation for the unknown function r(R), and as before, any other compatible anelastic factor can be expressed as QG, where Q is an arbitrary local rotation and G is as above.
Family 4: Inflation/Eversion of a Sphere. For this family, the symmetry enforced on the metric tensor automatically satisfies the universal equilibrium equations without additional restrictions. Demonstrating this, under this symmetry, the left Cauchy-Green tensor reads and its inverse is We can compute the invariants of b as r 4 m 2 1 (r) + 2r 2 m 2 2 (r), Notice in particular, that these invariants only depend on r. The Cauchy stress is diagonal with components Taking the divergence of this tensor and setting it equal to zero gives Therefore, the undetermined pressure must only depend on r, and the components of ∇ b σ ab only depend on r. Using ∂p ∂θ = 0 and ∂p ∂φ = 0, and defining v a = ∇ b σ ab , we can compute V a c = ∇ c v a . For simplicity, we will only compute the off-diagonal components of this tensor. Note that Computing the off-diagonal components, we get Therefore, V a c is diagonal. Because m ab is also diagonal, we conclude that V bc is diagonal, and hence, is identically symmetric. Recognizing V bc as m ba V a c = m ba ∇ c ∇ d σ ad , this means that the equilibrium equations are automatically satisfied for an appropriate pressure field, because the antisymmetric part of m ba ∇ c ∇ d σ ad vanishes simply due to the symmetry of the tensor field M AB .
We now only need to satisfy the incompressibility condition det b = 1. Computing this yields or in referential variables, writing m 1 (R) = m 1 (r(R)), R 4 m 2 1 (R)m 4 2 (R) = 1. Therefore, the final form of the inverse material metric tensor for this family is Alternatively, introducing a coordinate rescaling as before, the material metric tensor takes the form and we obtain a differential equation that can be integrated to determine r(R). Then, taking a multiplicative decomposition of F = AG, we can express G using an orthonormal frame in the material manifold yielding which again can be left multiplied by an arbitrary local rotation Q if desired.
Family 5: Inflation, Bending, Extension, and Azimuthal Shearing of an Annular Wedge. For this family, we have the deformation as given in (35) with deformation gradient (36). Again, we compute the equilibrium condition m a[l ∇ k] ∇ b σ ab = 0, and look at the W 111 coefficient. This contains two independent equations: If I 1 is constant, this equation is satisfied, and if I 1 is not constant, we require M 13 (r) = 0, and M 12 (r) = − ABM 11 (r) Cr . If I 1 is constant, we examine the W 122 coefficient in the equilibrium equation and obtain AEM 13 (r) Ar ABM 11 (r) + CrM 12 (r) , which characterizes the generic solution.
The conditions on the components of the metric are sufficient to satisfy equilibrium, so we only have to satisfy incompressibility. The incompressibility condition in this case reads We solve this for M 22 (r) to obtain which gives the generic solution or in referential variables, writing M AB (R) = M AB (r(R)), Unlike the other families, the standard Euclidean inverse metric M AB (R) = diag{1, R −2 , 1} is not a member of the generic solution branch for this family. This is because this Euclidean metric yields a special case of the anomalous solution, having constant invariants.
In principle, we can rescale our coordinates and compute the form of the anelastic factor in a multiplicative decomposition for a member of this family as we have done for the previous families. However we will not concern ourselves with the multiplicative decomposition for this family, because its generic solution branch does not contain the solution without eigenstrain. As such, any continuous process based on this family beginning with zero eigenstrain will not lie in this solution branch, but rather on the anomalous branch, and so the multiplicative decomposition associated with this generic branch is of limited use.

The Anomalous Universal Solutions
As we have covered cases I and IIa in the previous section, we turn our attention to case IIb. The groundwork for this case, namely the spatial constancy of the strain invariants has already been laid in our analysis of case IIa in section 6. The analysis for each family follows the same general pattern, so we will merely outline these steps here in an example appearing in Family 5, then present the results. Details are given in Appendix B.
Step 1: For the anomalous solution, we start with the equations derived from the equilibrium conditions: four second-order linear differential equations for each family, involving the six undetermined components of the inverse metric tensor. By integrating the equilibrium conditions, up to two of these components can be expressed in terms of the other variables.
We take, as an example the deformation for which we compute the components of b ab as The first universal equilibrium equation is m a[l ∇ k] ∇ b b ab = 0, which in these coordinates amounts to the two equations r 5M 12 (r) + M 11 (r) + rM 12 (r) + 3M 12 (r) + 3M 11 (r) = 0.
The general solution of these equations is Step 2: After integrating these equations, we have the three constant invariant conditions for each family to solve. The constant trace condition is linear in the unknown components of the inverse metric, so we can use it to solve for one undetermined inverse metric component in exchange for introducing the trace of b as a parameter.
For the purposes of our example, we will take γ 2 = 0, γ 1 = 0, and M 13 (r) = r. With this, b ab becomes Next, we compute the equilibrium condition m a[l ∇ k] ∇ b c ab = 0, which is simplified by multiplying c ab by the condition det b = 1. This condition puts the remaining two ODEs in the form 3M 12 (r) + 3M 11 (r) + r 5M 12 (r) + M 11 (r) + rM 12 (r) = 0.
Integrating these equations gives the solutions We have the constant trace condition hence, The incompressibility condition det b = 1 then can be written as β 1 r 6 + β 2 r 4 M 11 (r) 2 + r 2 µ 1 + µ 2 2 + r 3 β 1 + rβ 2 2 − I 1 β 1 r 6 + β 2 r 4 M 11 (r) and the constant second invariant condition can be written as Step 3: We are left with two nonlinear algebraic equations. The first is the incompressibility condition det b = 1, and the second is the constancy of the second invariant of b. Both are quadratic equations in the remaining component of the inverse metric tensor, which creates an overdetermined system. We compute the resultant of these two equations in this component, and demand this resultant vanish to ensure that these two equations have a common root. The resultant of these equations is itself a polynomial in the other undetermined integration constants: the invariants of b, and the remaining independent spatial coordinate. Even in our simplified example, the resultant of these equations in M 11 (r) yields the condition a relatively lengthy polynomial equation p(r) = 0. It can be immediately simplified by noticing that one of the coefficients is simply µ 6 2 , so µ 2 = 0 is a necessary condition for there to be a common solution to these two equations. It can be further simplified by noting that after using µ 2 = 0, one of the coefficients becomes β 6 2 , so we demand β 2 = 0. With this, a different coefficient becomes µ 6 1 , and hence µ 1 = 0 as well.
Step 4: Therefore, the resultant is a polynomial equation of the form p (q) = 0, which must hold for all values of the independent variable q (which is either r or x depending on the family). Accordingly, we set each coefficient to zero independently, and obtain an overdetermined system of nonlinear polynomial equations for the undetermined constants. We wish to find all the solutions to these equations, and so we compute a primary decomposition of the radical ideal generated by these equations. These equations are simple enough that this can be done with the assistance of a symbolic algebra package, though even then the computations are rather cumbersome (see appendix B). After we have done this, we are left with a set of conditions on the undetermined constants that are necessary and sufficient for the existence of a common root of the original quadratic equations in the undetermined inverse metric component over an open set. We substitute these constants into these equations and use them to solve for the final component of the inverse metric tensor, which gives us the general form of the anomalous solution. In all of these cases, despite encountering branching conditions in the course of analyzing the conditions on the constants, the separate branches ultimately are redundant, and we are left with a single anomalous solution branch for each family.
For our example, using the conditions µ 2 = 0, β 2 = 0, and µ 1 = 0, the polynomial simply becomes β 3 1 − I 1 β 2 1 + I 2 β 1 − 1 2 r 18 = 0, which demands β 3 1 − I 1 β 2 1 + I 2 β 1 − 1 = 0, because r > 0. We recognize that this is the eigenvalue equation for the tensor b, so we require β 1 to be an eigenvalue of b. We can satisfy this by writing I 1 = β 1 + e 1 and I 2 = β 1 e 1 + 1 β1 , where e 1 = λ 1 + λ 2 is the sum of the other two eigenvalues of b and we have used incompressibility in the form λ 1 λ 2 β 1 = 1. When we substitute these conditions back into the original equations for M 11 , they both become (up to some nonzero constant) This equation can be solved for M 11 and we obtain which gives one example of M AB as This lets us compute the corresponding elastic left Cauchy-Green stretch tensor as We can verify that b ab satisfies the equilibrium conditions and the constant invariant conditions. Completely determining the universal anelastic extensions of these families amounts to doing a similar analysis for each of the remaining families, but in full generality, i.e., not assuming particular values for the parameters appearing in the deformation, nor selecting values for the integration constants a priori.
These computations are included in the Appendix B, only the results are presented next. As we are considering the case where the strain invariants are constant, the set of universal equilibrium equations reduces to (48) and (49). Ordinarily these equations have three independent components, but in our case, one of these components vanishes identically for each equation, hence we have four ordinary differential equations to solve for each anomalous branch, together with the three algebraic equations constraining the strain invariants to be constant.
Family 1: Bending, Stretching, and Shearing of a Rectangular Block. Integrating the equilibrium equations (48) and (49) and solving the constant invariant conditions gives the following anomalous solution branch for this family: Here the constants e 1 and e 2 are the elementary symmetric polynomials in two of the three free eigenvalues of b e 1 = λ 1 + λ 2 , e 2 = λ 1 λ 2 , with the incompressibility condition determining the third eigenvalue as λ 3 = 1 λ1λ2 . The parameters e 1 and e 2 must be positive with e 2 1 > 4e 2 , since b is positive definite. The remaining constants, α 1 , α 2 , γ 1 , and γ 2 are arbitrary, subject to the condition that the choice of e 1 , e 2 , α 1 , α 2 , γ 1 , and γ 2 must yield a positive-definite metric tensor. One can explicitly verify that the invariants of b generated by this metric are I 1 = e 1 + 1 e 2 = λ 1 + λ 2 + λ 3 , Additionally, we can express this in terms of the referential variables by expressing r in terms of X by the relation r = A (2X + D).
Family 2: Straightening, Stretching, and Shearing of a Sector of a Cylinder. With this family, it is prudent to make the substitution ξ = x − D, which allows us to express the anomalous solution branch as Alternatively, we can express this in terms of referential variables using the equation ξ = 1 2 AB 2 R 2 . As in the previous case, e 1 and e 2 are the elementary symmetric polynomials in the free eigenvalues of b, λ 1 and λ 2 .
With this, the third eigenvalue is determined via the incompressibility condition λ 3 = 1 λ1λ2 . This ensures that the invariants of b are The constants e 1 , e 2 , α 1 , α 2 , γ 1 , γ 2 are largely arbitrary, apart from the condition that e 1 > 0, and e 2 1 > 4e 2 > 0, and that the constants are chosen such that the metric tensor is positive definite. Alternatively, we have the case where the anelastic strain is compatible, and we have where {M 11 , M 12 , M 13 , M 22 , M 23 , M 33 } are constants and det M AB (R) = R 2 . At first glance this case appears slightly more general than the previous one under the special case α 1 = γ 1 = 0, because for a fixed overall deformation, there are five independent parameters determining this solution, while setting α 1 = γ 1 = 0 in the other family yields a special case of (165) depending on the free parameters α 2 , γ 2 , e 1 , and e 2 . However, this case causes the stretch tensor b ab to be constant, which requires that the material manifold be Euclidean, and a constant isochoric stretch only depends on two independent stretches, with the remaining degrees of freedom representing a global rotation, which we can freely add or remove. This case appears to not be a special case of the previous branch, because once we eliminate the dependence on ξ, we not longer have a preferred direction, and hence we spontaneously gain additional rotational degrees of freedom that can be removed by the choice of the orientation of our current configuration Cartesian coordinates. Physically, this amounts to the reference configuration deforming anelastically into a parallelepiped, which can be elastically deformed into the desired block, as that elastic deformation is homogeneous. Indeed, the stress required to accomplish this is always constant, and hence equilibrium conditions are trivially satisfied. One can easily verify that the only nonzero Christoffel symbol generated by this metric is Γ 1 11 = 1 R , which generates a vanishing curvature tensor R = 0. In fact, the anelastic strain can be integrated up to an arbitrary rigid rotation and translation to obtain the position vector where ε a is an arbitrary right-handed set of linearly independent vectors spanning a parallelepiped with unit volume. With this, the constants M ab = ε a · ε b , i.e., the arbitrary constants appearing in the metric tensor are given by the Euclidean inner products of the constant basis vectors.
As with the other families, we can use the deformation equation r = √ AR 2 + B, to recast this into the referential variables. Additionally, the parameters e 1 and e 2 are the same as the previous families, and other than demanding that the eigenvalues they determine be positive, we also demand that the choice of variables α 1 , α 2 , γ 1 , γ 2 , e 1 , and e 2 leaves the metric tensor positive definite.
Family 5: Inflation, Bending, Extension, and Azimuthal Shearing of an Annular Wedge.
To facilitate the analysis of this family, it is useful to define the function f (r) = γ 1 + γ2 r 2 . With this, we have Again, the constraints on the constants appearing are as before, and are only necessary to ensure the positive definiteness of b and the metric tensor. We can recast this into referential variables using r = AR, if desired.

Merging of Universal Solution Families
After obtaining the previous results, it is natural to ask if solutions in one family correspond to solutions in another, and if so, to what extent? It is possible that the material manifolds, and the corresponding elastic deformations from two different families differ only by a change of coordinates, or equivalently, by a compatible anelastic deformation connecting the reference configurations of the two total deformations.

Equivalent Universal Solutions
Two universal deformations, ϕ 1 : or, equivalently, the following diagram commutes: For general manifolds, this is a difficult task, as we not only have to determine whether or not two isometries exist, i.e., solve the Riemannian manifold equivalence problem twice, but also whether or not they satisfy equation (180). However, in our case, the current configurations of the universal deformations are both Euclidean, and hence ψ must be an element of SE(3). Additionally, ϕ 1 and ϕ 2 are invertible and in principle known, so if we can find ψ, we can solve for Ψ = ϕ −1 2 • ψ • ϕ 1 . It is then a simple matter of checking whether or not Ψ is an isometry.
We choose coordinates for all of these manifolds, writing current configuration coordinates as {x a }, and material manifold coordinates as {X A }, with each set numbered by the universal deformation pertaining to it. In terms of these coordinates, these maps are where we have used different indices on the different sides of the equations to emphasize that in principle each new coordinate depends on all of the old coordinates. These maps prolong to tangent maps (F 1 ) In terms of these tangent maps, we then have the isometry conditions and the prolongation of equation (180) as Because both current configurations are Euclidean, we can trivially satisfy equation (185) by choosing ψ to be an element of SE (3), and we can then use equation (186) to express equation (184) We can write this expression in terms of the inverse of and hence

Relationships Between the Six Universal Families
We would like to identify which families are likely to contain overlap, and take note of Table 2. Specifically, the left Cauchy-Green tensor of each family is symmetric with respect to the prolonged action of a subgroup of SE(3). Therefore, if two universal deformations are equivalent, their corresponding strain tensors should have isomorphic symmetry groups. Denoting the symmetry group of b 1 as K 1 ⊂ SE (3), and the symmetry group of b 2 as K 2 ⊂ SE(3), we seek ψ ∈ SE(3) such that This immediately identifies a possible correspondence between Families 1, 3, and 5, because their symmetry groups are isomorphic. Additionally, we expect that there might be some universal solutions in Family 2 that are also in Family 0, since the symmetry group of Family 2 is a subgroup of that of Family 0, though we can immediately recognize that there are solutions in Family 2 that are not equivalent to any in Family 0, because not all solutions in Family 2 are invariant under the action of the full symmetry group of Family 0. This observation immediately reveals that, up to an element of these symmetry groups, ψ must be the obvious one implied by our choice of coordinates in each family, because it must send invariant sets of K 1 to invariant sets of K 2 . We recall that if a (sub)group K acts on a manifold M, an invariant set of K is a set S K ⊂ M such that ∀x ∈ S K , and ∀k ∈ K, k • x ∈ S K . Here we consider the smallest nonempty invariant sets: the orbits of a single point under the action of the subgroup K i . The invariant sets of the symmetry groups of Families 1, 3, and 5 are concentric cylinders, hence any potential ψ connecting these two families must map a family of concentric cylinders to another. The coordinates for each family were chosen such that this family of cylinders is centered on the z axis, hence we require ψ to be a Euclidean isometry mapping the z axis to itself. Apart from rotations and translation that leave the left Cauchy-Green tensor fields unchanged, this restricts ψ to either be the identity, or a rotation reversing the orientation of the z axis. We will see that we can freely take ψ to be the identity.
We first show that Family 0 is contained within Family 2. To do this, we must find an equivalent deformation in Family 2 for any choice of deformation in Family 0. Identifying our coordinate systems (i.e., taking ψ to be the identity), we can express the left Cauchy-Green tensor field for any deformation in Family 0 as We choose a universal solution in Family 2 with material inverse metric of the form withM AB being appropriate constants, which is one of the cases where the material manifold is Euclidean. Pushing this forward to the current configuration, we obtain the equations Therefore, for any given element of Family 0, the choices yield an equivalent member of Family 2. Also we note that these compatible material manifolds are contained as special cases of the non-homogeneous branch of Family 2 via the same argument presented in section 7. Denoting U A to be the set of universal deformations corresponding to family A, we conclude that We then seek to establish similar correspondences between the sets U 1 , U 3 , and U 5 . First, we consider an element of U 5 lying in its generic branch. The left Cauchy-Green tensor field of this element is fully determined by specifying three functions of R, hence implicitly of r through R = A r , namely M 11 (r), M 23 (r), and M 33 (r), along with values for the constants A, C, and E. Labeling these choicesM 11 (r), M 23 (r),M 33 (r),Ã,C, andẼ, we seek elements in Families 1 and 3 that generate the same stretch tensor field.
The left Cauchy-Green tensor field for the generic branch of Family 1 depends on three arbitrary functions of X(r) = r 2 2A − D 2 , M 11 (X(r)), M 22 (X(r)), and M 23 (X(r)) as well as the constants A, B, C. If we select the functions and constants such that it is straightforward to verify that the stretch tensor fields generated coincide. Therefore, the generic solution branch of Family 5 is contained in the generic solution branch of Family 1, since we can find universal solutions in Family 1 that are equivalent to any universal solution in Family 5.
Similarly, the generic branch of Family 3 depends on three functions of r through R(r) = r 2 −B A : M 11 (R(r)), M 22 (R(r)), and M 23 (R(r)) as well as the constants A, B, C, D, E, and F . The choice M 11 (R(r)) =Ã 2 r 2M 11 (r) also generates an identical stretch field, hence the generic branch of Family 5 is also contained in the generic branch of Family 3. We have shown that the generic branch of Family 5 is contained in those of both Family 1 and Family 3. To examine the opposite direction, suppose we take an arbitrary member of the generic branch of Family 1, defined by parametersM 11 (X(r)),M 22 (X(r)),M 23 (X(r)),Ã,B, andC, and seek to find a solution in Family 5 that generates the same stretch tensor field. Elements in Family 5 depend on the parameters M 11 (R(r)), M 23 (R(r)), M 33 (R(r)), A, C, and E, and the choice generates the same stretch tensor fields as the member of Family 1. Hence, the generic solution branch of Family 1 is contained in the generic branch of Family 5. Coupled with the previous result, we conclude that the generic solution branches for Families 1 and 5 are equivalent, in that every universal solution in one of these branches has at least one equivalent universal solution in the other. Next, choosing an arbitrary universal solution in the generic branch of Family 3, we seek a universal solution in Family 5 that is equivalent. Choosing parametersM 11 (r),M 22 (r),M 23 (r),Ã,B,C,D,Ẽ, andF determining an arbitrary solution in Family 3, we can choose an element of Family 5 by specifying the parameters A, C, E, M 11 (R(r)), M 23 (R(r)), and M 33 (R(r)), where R(r) = r A . If we choose these such that we obtain a universal solution that is equivalent to the specified solution in Family 3. Hence, the generic solution branch of Family 3 is contained within that of Family 5. Coupled with our previous results, this result means that the generic solution branches of Families 1, 3, and 5 are all equivalent to each other. Next we consider the anomalous solution branches for these families. First, we select an arbitrary member of Family 3 anomalous solution branch by specifying the parametersÃ,B,η =DẼ −CF , 5ẽ 1 , ẽ 2 ,α 1 ,α 2 ,γ 1 , andγ 2 . We seek to find solutions in Family 1 anomalous solution branch, and Family 5 solution branches that generate equivalent solutions. First examining Family 1, we can select values for constants α 1 , α 2 , γ 1 , γ 2 , e 1 , e 2 , A, B, and C. It is straightforward to verify that the choice α 1 =γ 2 , α 2 =γ 1 , γ 1 =α 1 , γ 2 =α 2 + Cγ 1 , B = 1 η Ã , A =Ãη 2 , e 1 =ẽ 1 , e 2 =ẽ 2 , (214) generates an equivalent solution. Likewise for Family 5, we can choose values for the parameters α 1 , α 2 , γ 1 , γ 2 , e 1 , e 2 , A, and E, where the specific choices A =η Ã , E = 1, α 1 =α 1 , α 2 =α 2 , γ 1 =γ 1 , γ 2 =γ 2 , e 1 =ẽ 1 , e 2 =ẽ 2 , generate a solution that equivalent to the arbitrary solution from Family 3. Hence the anomalous branch from Family 3 is contained in both that of Family 1 and Family 5. Conversely, we select an arbitrary member of the anomalous branch of Family 5 by specifying the parametersα 1 ,α 2 ,γ 1 ,γ 2 ,ẽ 1 ,ẽ 2 ,Ã, andẼ. We can verify that the choice of parameters yields a solution from Family 3 that is equivalent to the arbitrary one from Family 5. Finally, we select an arbitrary member of the anomalous branch of Family 1 by specifying the parametersα 1 ,α 2 ,γ 1 ,γ 2 ,ẽ 1 ,ẽ 2 ,Ã,B, andC, and seek an equivalent solution in Family 3. The parameter choices A =Ã 2 , B =B, η = 1, e 1 =ẽ 1 , e 2 =ẽ 2 , generate such a solution. Hence, we deduce that the anomalous branches of Families 1 and 5 are contained in that of Family 3, which combined with our previous results implies that the anomalous branches of all the three families are the same. Therefore, having examined both the generic and anomalous branches of these families, we conclude that U 1 = U 3 = U 5 . Hence, in the anelastic setting, our initial six families of universal solutions have collapsed into three families U 2 , U 3 , and U 4 , one corresponding to each of the three surfaces with constant principal curvatures in 3D Euclidean space: planes, cylinders, and spheres, respectively. These surfaces are the invariant sets of the symmetry groups of the left Cauchy-Green tensor fields, and they played a central role in [Ericksen, 1954], being the level sets of the strain invariants. Here, we see that not only are the invariants of b constant on these surfaces, but b itself is symmetric with respect to these surfaces in the manner induced by the action of the special Euclidean group. This symmetry is present even in the degenerate case when the invariants of b are constant, which is why we can identify the symmetry groups even in the anomalous solution branches. In the classical problem, similar surfaces can be identified in the material manifold, since in the absence of eigenstrains, the material manifold and the reference configuration coincide. These surfaces are invariant sets of the symmetry groups of the right Cauchy-Green tensor fields, and prevent the identification of the classical families with each other, since the only two classical families with matching invariant sets in both configurations are Families 3 and 5. These however cannot be identified with each other because solutions in Family 5 have constant invariants, while those in Family 3 do not. Hence, it is only after the addition of eigenstrains that many of the classical families become redundant.

Standard Forms of the Three Distinct Universal Families
We note that there is some redundancy in the parameterizations we currently have, which is exhibited by observing that the parameter selections we have used to identify the families with each other are not mutual inverses. We can reparametrize to eliminate this redundancy and have a single representation for strain field of each family. Concretely, we can express the left Cauchy-Green stretch field for the anomalous branches of U 3 in the following standard form: where p is the product of the two free eigenvalues of b, and m is the mean of the two free eigenvalues. The inverse of this, c ab , is the push forward of the material metric, and has components The generic branch of this family likewise has a standard expression: with the incompressibility condition being c 11 (r) c 22 (r)c 33 (r) − (c 23 (r)) 2 = r 2 . The positive definiteness condition is equivalent to c 11 (r) > 0, c 22 (r) > 0, and c 33 (r) > 0 in addition to the incompressibility condition (221), or in the anomalous solution, requiring m > 0, and 0 < p < m 2 . An example of one of these generic solutions was investigated by Yavari and Goriely [2015], with the parameter choices c 11 (r) = λ 2 , c 22 (r) = λ 2 r 2 , c 23 (r) = r 2 (ψ (λr) − τ ), and c 33 (r) = 1+λ 2 r 2 (ψ(λr)−τ ) 2 λ 4 . Similarly, the left Cauchy-Green tensor field for the anomalous branch of the family U 2 takes the standard form with p and m defined as previously. Inverting this to obtain c ab , one obtains The left Cauchy-Green tensor for the generic branch of this family also has a standard form with the incompressibility condition becoming The inverse of this then takes the standard form with the incompressibility condition The positive-definiteness condition is equivalent to requiring c 11 (x) > 0, c 22 (x) > 0, and c 33 (x) > 0, in addition to the incompressibility condition (226), or in the anomalous case, requiring m > 0, and 0 < p < m 2 . Finally, the spherically-symmetric family U 4 can be expressed in the standard form through its left Cauchy-Green tensor which has the inverse The incompressibility condition and positive definiteness is automatically satisfied for arbitrary functions g(r) satisfying g(r) > 0. In terms of parameters defined by Goriely [2017] in Chapter 15.1.1, this function is g(r) = α r = α −2 : it is the radial stretch. These standard forms make it clear that universal solutions in anelasticity can be categorized by computing the tensor c , and comparing the result with the standard forms here. As a consequence of this, the symmetry of the elastic strain in the current configuration determines which family any particular universal solution belongs to, as it is this symmetry that is reflected in c .
We have examined particular symmetry groups, namely T(2), T(1)×SO(2), and SO(3). All of these are Lie subgroups of SE(3), and specifically they are generated by two independent generators; choosing two translational generators yields T(2), choosing a rotation and a translation about the axis of that rotation yields T(1)×SO(2), and choosing two rotations fixing a common point yields SO(3). We show in Appendix C that any Lie subgroup of SE(3) generated by two arbitrary independent generators contains at least one of these groups by necessity, hence we have the following theorem: Theorem 8.1 (Classification of Symmetric Universal Solutions). Any universal solution that is equivariant under the action of two independent 1-dimensional Lie subgroups of SE(3) is contained in one of the three universal families U 2 , U 3 , or U 4 . This allows us to precisely state our conjecture regarding the completeness of our classification in terms of symmetry: Conjecture 8.1 (Symmetry Necessity). A deformation must be equivariant with respect to the action of two independent 1-dimensional Lie subgroups of SE(3) in order to be universal, hence our classification is complete.

Graphic Representation
Because the material manifolds are generally non-Euclidean, visualizing them is difficult. A way to overcome this difficulty is to approximate their geometry as "piecewise Euclidean" and examine the deformation of each piece. This approach is similar to the three-dimensional version of approximating a curved surface with a polyhedron, and then representing that polyhedron in the plane by its net. The original surface can then be build up by connecting appropriate edges, but because of the curved nature of the surface, these edges cannot all be connected without distorting the pieces, or lifting them out of the plane. To demonstrate this technique, we will first start with a two-dimensional example, and then move on to a Euclidean three-dimensional example, and then finally apply the techniques to examples of material manifolds obtained from our analysis.

A two-dimensional example
We know that representing spherical geometry in the plane isometrically is an egregiously impossible task [Gauss, 1828]. To get around this, we only do this approximately, and allow for incompatibility by partitioning and separating our domain in multiple pieces. We can then stretch each piece in such a way that the deformed pieces can be approximately recombined in three-dimensional space to form an upper hemisphere. The deformed pieces are individually flat, so they can all be placed in the plane, but not in a way such that they can be pieced together without gaps (see Figure 12). Figure 12: We start with a disk, partition it, and separate the resulting pieces to allow for room for each piece to strain without overlapping its neighbors. We then strain each piece, and recombine the deformed cells to create an approximation of an upper hemisphere. Here each cell [R i−1 , R i ] × [Θ j−1 , Θ j ] is positioned so both the position of the point (r, θ) = ((R i−1 + R i )/2, (Θ j−1 + Θ j )/2), and the orientation of its tangent plane match that of the exact map from the disk to the hemisphere.
Explicitly, we want to take the region r ∈ [0, 1], θ ∈ [0, 2π), where r and θ are polar coordinates in the plane, and map it to the surface z = √ 1 − r 2 in three-dimensional space. The stretch induced by this map is described by the metric tensor with cylindrical components which we can approximate as constant on each piece, while keeping each piece in the plane, by evaluating the metric at r * = (r max − r min )/2. The deformed pieces can then be rigidly translated and rotated in Figure 13: As the partitioning gets finer, the resulting approximation of the anelastic strain becomes more and more accurate.
three-dimensional space to approximate the desired spherical surface, with the approximation becoming better as the partition becomes finer (see Figure 13). This two-dimensional example allows one to see the correspondence between the deformed partitioned approximation and the recombined non-Euclidean configuration, which is important because once we move up to three-dimensional examples, we are no longer able to recombine the strained pieces; we must deduce properties of its geometry from the deformed partitioned approximation alone. Additionally, while we assembled the resulting deformed pieces into a hemisphere by lifting them into a higher dimensional Euclidean space, we could have assembled them into any number of other surfaces that are isometric to the hemisphere. Because we only determine the intrinsic geometry of the material manifold, there is no preferred isometric embedding in some higher dimensional Euclidean space, unless as above, we explicitly specify the embedding.

A three-dimensional Euclidean example
Just as in the two-dimensional example, we can partition a flat three-dimensional body, explode it, and approximate strains on each piece to represent the non-Euclidean geometry of our deformations. The only difference is that, in general, we cannot recombine the distorted pieces into a cohesive whole, because the resulting shape is not globally flat. However, if the strain that we impose is actually induced by a map between Euclidean spaces, we can apply this procedure to the partitioned pieces, and observe the local strain, while separately observing the global deformation. We can then compare the two results to see which features are preserved by this local partitioning approach to better interpret the results of applying this procedure to our derived material metrics.
Consider the following map given in cylindrical polar coordinates This map produces azimuthal shear and angular stretching. Choosing µ = 2 and ν = 1 2 , and mapping the domain R ∈ [2, 3], Θ ∈ 0, 2π 3 , Z ∈ [0, 1], we obtain a transformation shown in Figure 14. If we instead compute the strain tensor, and use it as a material metric for the current configuration, we obtain Applying this stretch to our partitioned domain, we obtain the depiction shown in Figure 15. This side-by-side comparison shows what is happening when we do this piecewise transformation, namely we capture the strain of each piece, but we do not capture any local rotation that is present in the global deformation. This is because local rotations produce no strain, so they do not contribute to the strain tensor, and hence, we cannot expect to be able to capture them through this reconstruction.

Anomalous anelastic strains
For the anomalous families, we will use the partitioning technique to attempt to visualize the deformations. We note that not all choices of parameters are valid over arbitrary domains. In particular, the parameters must be chosen such that the metric is positive definite over the chosen domain, in addition to making sure that the strain tensor b is positive definite, a much simpler task as this requires e 1 and e 2 to be positive with e 2 1 > 4e 2 . We also depict the total overall deformation, coloring the current configuration by the trace of the Cauchy stress required to maintain it for a Mooney-Rivlin solid with strain energy of the form both in the presence and absence of anelastic strain. Because the invariants of the left Cauchy-Green tensor are constant for the anomalous universal eigenstrains, any choice of strain energy is indistinguishable from a Mooney-Rivlin energy, and the only invariant of the Cauchy stress that can potentially vary spatially is the pressure generated due to the constraint stress. Here we choose the two material parameters in the Mooney-Rivlin energy to each be equal to 1, though different choices of parameters would yield qualitatively similar results.
Family 1. For this family we choose the reference domain X ∈ [0, 1] , Y ∈ [0, 6], and Z ∈ [0, 4], and take the deformation parameters A = 3 2 , B = 1, C = 1 4 , D = 2, E = 0, F = 0. To examine the effects of anomalous universal eigenstrain on the equilibrium stress distribution, we consider the same overall deformation, and contrast the stress generated in the presence of eigenstrain with that generated in the absence of eigenstrain. For the anelastic strain parameters, we use , e 1 = 9 4 , e 2 = 9 8 , One can verify that this ensures that M is positive definite over the chosen domain. To visualize the anelastic strain, we subdivide the domain, separate the pieces, and approximate the anelastic strain on each (see Figure 16). This anelastic strain is generally not compatible, i.e., the deformed pieces cannot be reassembled in Euclidean space without further deformation. We map the body into the current configuration, and color it to denote the spherical part of the Cauchy stress generated by a Mooney-Rivlin solid. This requires us to integrate the indeterminate constraint pressure field, both with and without eigenstrain. Without eigenstrain, we have the following differential equations for the constraint pressure which can be easily integrated to obtain Notice in particular that p does not vary with z or θ. Additionally, only the gradient of p affects the motion, which allowed us to ignore the integration constant when integrating the above equations. When there is eigenstrain, we obtain a different set of differential equations determining p(X): where k(R) is an algebraic function of r alone. We can in principle integrate these to obtain p eig (X) = p eig (r, θ, z) = r r0 k (r) dr + 2AB (1 + e 2 ) e 2 α 2 θ + 2 (1 + e 2 ) e 2 γ 1 z.
In contrast with the ordinary case, when we have eigenstrain, we can generate pressure gradients that vary with θ and z. Interestingly enough, even the generic universal eigenstrain cannot generate pressure gradients in these directions; the anomalous universal eigenstrain is the only universal eigenstrain that can create pressure gradients in these directions. This suggests that the measurement of these pressure gradients can be used to partially measure the eigenstrain, and conversely these anomalous solutions can be used to generate pressure gradients in these directions to, for example, counteract the pressure generated by body forces. We then compute the first stress invariant, the trace of the Cauchy stress, or equivalently its spherical part, for the material both in the absence and presence of eigenstrain, and plot the resulting stress invariant in Figure 17.
Family 2. For this family, we choose the reference domain R ∈ [2, 3], Θ ∈ [0, 5], and Z ∈ [0, 4], and take the deformation parameters A = 1, B = 3 4 , C = 1 4 , D = 0, E = 0, F = 0. These parameters define the total deformation, allowing us to examine the effects of eigenstrain on the Cauchy stress. In particular, we take the parameters appearing in the anelastic strain to be , e 1 = 21 10 , e 2 = 9 10 , This ensures that the metric is positive definite over the chosen domain. We subdivide and explode the domain, and apply our anelastic strain to each piece, as shown in Figure 18. We are then left with a set Figure 18: A depiction of the anomalous anelastic strain for Family 2.
of differential equations determining the constraint stress. In the absence of eigenstrain, we have Upon integration, we obtain .
In contrast, when we consider eigenstrain, we have the equations where as before k(x) is an algebraic function of x. We can integrate these equations to obtain As before, the presence of this anomalous universal eigenstrain generates pressure gradients in directions that are not possible in their absence, in this case, the y and z directions. Again, even the generic branch of universal eigenstrains cannot generate pressure gradients in these directions, further highlighting the unique nature of the anomalous solutions. We then compute the trace of the Cauchy stress, and plot the resultant distributions both in the absence and presence of eigenstrain (see Figure 19).  , H = 0. This completely defines the total deformation, allowing us to examine the stress generated with eigenstrain in contrast with that generated without eigenstrain. We take 10 , e 2 = 9 10 , as our anelastic parameters. Over the defined domain, these choices ensure that the anelastic metric tensor is positive definite. We then partition and explode the domain, approximating the eigenstrain on each piece, depicting the result in Figure 20. As before, we obtain differential equations for p, yielding in the elastic case which can be integrated to obtain In contrast, in the presence of eigenstrain, we have the differential equations with k(R) an algebraic expression in r as in other families. This pressure can be integrated to obtain We compute the first invariant of the Cauchy stress and color the deformation according to it in Figure  21. As in other families, the presence of this anomalous branch of universal eigenstrain can generate pressure gradients in directions that do not occur otherwise, specifically pressure that varies with θ and z. This property would allow one to indirectly measure the eigenstrain by measuring the pressure variation required to sustain this deformation. Likewise, if we can specify the eigenstrain, we can create specific pressure gradients that would otherwise be impossible for the generic branch. and Z ∈ [0, 4], and we take the deformation parameters to be A = 4 5 , B = 1, C = 1, D = − π 4 , E = 5 4 , F = 0. Choosing the anomalous eigenstrain parameters as we subdivide and explode our domain, then approximate the eigenstrain on each piece. The result of this is depicted in Figure 22. As with the other families, we are left with a set of differential equations determining the constraint stress. When eigenstrain is absent, we have the equations which can be integrated to obtain Figure 22: A depiction of the anomalous anelastic strain for Family 5.
When we consider the anomalous solution, we have the equations in terms of an algebraic function k(R) which can be integrated to obtain When we compute the pressure gradient in the case of the generic universal solution, we obtain a pressure that only varies with r. We see that the anomalous branch generates pressure gradients that vary with θ and z, unlike the generic solutions. We can compute the trace of the Cauchy stress, both with eigenstrain and without eigenstrain, and use this to color the deformed configurations in Figure 23. Notice that unlike other families, the constraint pressure in the absence of eigenstrains can vary in a direction different than in the generic anelastic situation, specifically, for an eigenstrain in the generic solution branch, we have ∂p/∂θ = 0, but in the absence of eigenstrain, we have ∂p/∂θ = 2B A 4 + 1 C 2 /A 2 . This is due to the fact that the standard Euclidean metric in terms of cylindrical polar coordinates lies on the anomalous solution branch, not within the generic branch. In the elastic case, when the azimuthal shearing term does not vanish, corresponding to B = 0, we have a pressure variation in θ that is necessarily nonzero. In the anomalous case, we can determine the pressure variation in both θ and z by adjusting our choices for γ 1 and α 1 , allowing us to determine the pressure variation in these directions independently of the overall total deformation by choosing the anelastic strain appropriately. Specifically, we can have arbitrarily large azimuthal shearing, while also causing the azimuthal pressure variation to vanish. While the other families also allow us to select the ordinarily absent components of the pressure gradient in a similarly arbitrary way, this family is unique in having one of these pressure gradients present without eigenstrain, so the anomalous universal eigenstrain allows us to both create pressure variations in these directions, but also remove pressure variations that are ordinarily necessary to maintain the overall deformation.
Additionally, with the anomalous anelastic solution branch, the azimuthal pressure variation ∂p/∂θ = 2A (1 + e 2 )γ 1 /e 2 does not depend on the degree of azimuthal shearing, i.e., it is independent of B. While the anomalous eigenstrain itself does depend explicitly on B, if the eigenstrain and the total deformation are simultaneously varied by changing B, the azimuthal pressure gradient should not change. Doing this in practice would be difficult, because fundamentally the parameter B partially determines the overall deformation, hence both the overall deformation and the eigenstrain would have to be simultaneously controlled in precise ways to realize this thought-experiment. Thankfully, this does not have to be done dynamically; a new value of B could be selected, the overall deformation could be controlled, and once it is established, fixed. Then the eigenstrain could be controlled until the universal eigenstrain corresponding to the chosen value of B is obtained. After this is done, the pressure variation could be measured, and this process can be repeated to establish the independence of the azimuthal pressure gradient.

Conclusions
We have generalized the universal solutions of Ericksen's problem to the case of anelasticity. The main idea was to first indentfy the symmetry group associated with each solution of the classical Euclidean With eigenstrain (bottom), an axial pressure gradient can also be sustained.
problem and use this symmetry group in a non-Euclidean setting by finding the possible metrics that guarantee each symmetry group. We used both the structure of existing universal solutions for a given M and their symmetries to find possible material metrics M. This was done by interpreting the classical universal deformations passively as coordinate changes. Then all local changes in geometry can be captured by changing the metrics. In this way, once we moved away fromM to the general M, we recognized that any homeomorphism can be expressed by the particular coordinate maps for each family, since the coordinates themselves no longer hold any specific interpretation. By identifying an appropriate symmetry to impose on M for each family a priori, we accomplished the following: • First, we constrained this problem to a point where it is still nontrivial, but solvable in a systematic algorithmic way.
• Second, because this symmetry depends on the classical family of universal solutions, our construction provided a direct extension and classification of the new anelastic solutions consistent with the elastic ones.
• Third, it is likely that these symmetries play a fundamental role in constraining the classical problem; the known classes of universal deformations all have particular symmetries. Identifying the relevant symmetries to impose on M highlights this explicitly and we conjecture that all possible cases of anelastic universal solutions possess such symmetry. Specifically, all known universal solutions are preserved under the induced action of subgroups of the special Euclidean group. These subgroups are precisely those having two-dimensional invariant sets, which are either parallel planes, concentric cylinders, or concentric spheres.
It should also be noted that the generic solution branches and the anomalous solution branches differ largely in character. The generic solution branches contain arbitrary functions, and as such are infinite dimensional, while the anomalous solution branches are entirely determined up to a handful (∼ 10) of arbitrary constants, 6 of which are not redundant, and have a highly nontrivial structure. In all cases, the different branches of the analysis ultimately yield the same family of anomalous solutions. These anomalous branches also allow for the pressure required to sustain these deformations to vary in directions that would otherwise not be supported, even on the generic anelastic branch. This suggests possible applications of these anomalous branches in manipulating the surface tractions required to sustain these deformations, as well as a way to indirectly measure some of the eigenstrain parameters.
Additionally, symmetry appears to play an important role in these universal solutions. The right Cauchy-Green stretch tensor field for every family is invariant under some subgroup of SE(3), the group of orientation preserving isometries of 3D Euclidean space. The dimension of the Lie symmetry seems to play an important role as well; equilibrium conditions for families with three-dimensional Lie symmetries (Families 0 and 4) are trivially satisfied by imposing that symmetry on the material manifold, while families containing two-dimensional Lie symmetries require further restrictions.
While we have framed this problem in the context of an anelastic deformation from Euclidean space with the usual metric in the chosen coordinates to some Riemannian manifold, and a further elastic deformation back to Euclidean space, we do not make use of the initial Euclidean space in our analysis, nor do we detail any specific mechanism driving the anelastic deformation. The Euclidean reference configuration appears in the initial presentation of these universal solutions in nonlinear elasticity, therefore we use it as a comparison when displaying the elastic stress generated by these deformations, but there is no need for it to be the initial state of our undefined anelastic process. Provided our anelastic evolution arrives at the material manifolds derived here, the remainder of the deformation can be accomplished elastically. Therefore, the anelastic deformation can in principle map from any configuration; there is no need for a Euclidean reference. Even if the reference is Euclidean (as we tend to model physical space as Euclidean), there is no need for the coordinates used to be the typical Cartesian, cylindrical polar, or spherical polar coordinates used in the elastic case; they can be curvilinear coordinates, dramatically broadening the range of anelastic deformations to which our results are applicable.
Finally, we remark that neither the classical Ericksen's problem, nor the anelastic Ericksen's problem presented here has been proved to be fully solved and there may still be universal solutions unaccounted for. However, our conjecture, based on the correspondence between solution families and their symmetry groups, is that both classifications are actually complete. An additional work demonstrating that the strain fields of these universal solutions must by necessity be symmetric with respect to two one-dimensional Lie subgroups of SE(3) would prove that this classification is complete.

A Methods and Tools
In the course of our analysis, we shall use a few tools that are infrequently used in nonlinear elasticity.
Here we present a brief summary of these tools, and provide references to further sources for interested readers.

A.1 Algebraic tools
Some techniques from elimination theory [Cox et al., 1992] will be used in our analysis of universal solutions. Chiefly among these is the method of resultants. We will not use resultants in their full generality, but rather will only need to compute the resultant of two quadratic polynomials, and as such will only provide the details necessary to substantiate our usage. Consider two polynomials The resultant of p 1 and p 2 is a multivariate polynomial in {a 0 , ..., a k , b 0 , ..., b r } that vanishes if and only if there exists a common solution to the equations p 1 (x) = 0, and p 2 (x) = 0. This will be useful to us because we will consider multivariate polynomials recursively as single variable polynomials, with coefficients in some extended field, and the method of resultants gives us a way of reducing the number of necessary equations that must be satisfied, while reducing the number of variables we must consider.
As an example, consider two quadratic polynomials p 1 (x) = a 1 x 2 + b 1 x + c 1 , and p 2 (x) = a 2 x 2 + b 2 x + c 2 . We seek a condition on {a 1 , b 1 , c 1 , a 2 , b 2 , c 2 } such that there existsx satisfying p 1 (x) = p 2 (x) = 0. Taking the linear combination a 2 p 1 (x) − a 1 p 2 (x) = 0 gives us a linear condition onx, namely Then taking the combination b 1 p 2 (x)−b 2 p 1 (x) = 0 gives the condition (a 2 b 1 − a 1 b 2 )x 2 +b 1 c 2 −b 2 c 1 = 0, which under the above linear restriction becomes Cross multiplying and subtracting the two linear equations (254) and (255) gives a necessary condition on the coefficients of p 1 and p 2 for there to be a common solution, namely The left hand side is precisely the resultant of two quadratic polynomials, and so we denote The vanishing of the resultant is a necessary condition for the existence of a common root of the two quadratic equations in x.
Additionally, we will repeatedly use the fact that if a polynomial p(x) vanishes on an open set U , i.e., p(x) = 0 , ∀x ∈ U , then, by the fundamental theorem of algebra, all its coefficients vanish identically, i.e., p(x) is the zero polynomial.

A.2 Group action on manifolds
Symmetry will play an important role in our construction. As group theory lies at the heart of any discussion of symmetry, we present some definitions from the theory of Lie groups, and reference the reader to Gorbatsevich et al. [2013] for a full treatment. We recall that a semi-direct product is a generalization of a direct product, where only one factor must be a normal subgroup of the result. For instance, T(n), the group of translations of Euclidean space, is a normal subgroup of the special Euclidean group SE(n), while the group of rotations, SO(n) is not, hence SE(n) = SO(n) × T(n), but rather SE(n) = SO(n) T(n). The special Euclidean group, denoted SE(n), consists of all orientation preserving global isometries of Euclidean space, and is a semi-direct product of SO(n), and T(n). Therefore, an element of SE(n) can be identified with a tuple (Q|c) consisting of an element Q of SO(n), and an element c of T(n). The defining feature of SE(n) being how it acts on E n , we must now express this action, and hence the natural group operation of SE(n) on E n in terms of (Q|c).
The action of a group (G, ) on a manifold M, informally, is a map ρ : G × M → M that preserves the group structure of G. Denoting this action for g ∈ G and x ∈ M as ρ (g, x) = g • x, this demands (m 2 m 1 ) • x = m 2 • (m 1 • x) for all m 1 , m 2 ∈ G, and all x ∈ M. Additionally, denoting the identity of G as e, we demand ρ (e, x) = e • x = x, ∀x ∈ M. Concisely, a group G acts on an object M via a homomorphism ρ : G → Aut (M).

Example A.1
Treating E n as a vector space, i.e., fixing an origin, the action of SE(n) on E n in terms of the tuple (Q|c) sends the point x ∈ E n to the point Qx + c. One can easily verify that this is an isometry. The action of the element (Q 1 |c 1 ) followed by the action of the element (Q 2 |c 2 ) is then x → Q 2 Q 1 x+Q 2 c 1 +c 2 , hence, in this representation, the product of the special Euclidean group, , takes the form (Q 2 |c 2 ) (Q 1 |c 1 ) = (Q 2 Q 1 |Q 2 c 1 + c 2 ). We will eventually see that the right Cauchy-Green stretch tensor for each of the known families of universal solutions is preserved under the prolonged action of some Lie subgroup of SE(3).
The action of a group on a manifold can then be prolonged to its tangent bundle. This prolonged action can be determined by fixing an arbitrary group elementg, and considering the action of this element as a map ρ(g) : M → M. The existence of inverse elements in G guarantees that this map is invertible, since ρ(g) −1 = ρ(g −1 ). Provided M has a smooth structure, and ρ(g) is a smooth function of x, this map can then be differentiated to obtain the corresponding induced tangent map, i.e., the push forward map, which then determines the action ofg on the tangent bundle of M. The invertibility of ρ(g) then provides a group action on the cotangent bundle via the pull-back induced by the inverse map. Notice that, generally, we consider the prolonged action of a group on bundles over M, not merely on individual tangent spaces, because underlying points in the base space are not generally fixed, i.e., the base space and the total space are transformed together. Additionally, even if certain points in the base space are fixed points under the action of a group element, this does not guarantee that the tangent spaces at these points are similarly preserved. For example, in E 3 , a rigid rotation preserves the position of the points on its axis, but rotates the tangent spaces at those points.

B Explicit Calculations of the Anomalous Solutions
For families containing anomalous solution branches, we have the necessary (but not sufficient) condition that the invariants of b are constant, with det b = 1 for incompressibility. For each of these families, we have four linear differential equations, one linear algebraic equation, and two nonlinear algebraic equations for the six unknown functions comprising the components of M AB . We will use the linear equations to solve for five of these unknown functions in terms of the sixth, and then characterize the common solutions to the remaining two equations to determine the final component.

Family 1
For this family, we consider the case where the invariants of b are constant. This gives the equation which, when applying the constant invariant condition and forcing this to hold for all energy functions requires The first of these has two nonzero components after substituting the form of M AB . One of these components yields the differential equation 3M 12 (r) +rM 12 (r) = 0. This equation is readily integrated to obtain M 12 (r) = α 1 r 2 + α 2 .
We also have the constant trace condition on b, which becomes We can express this system of equations, linear in M 22 (r), M 23 (r), and M 33 (r), as the matrix equation which is invertible, because the determinant of the matrix on the left hand side is since M 12 (r) and M 13 (r) cannot simultaneously vanish. We invert these equations to obtain expressions for these components of the inverse metric in terms of M 11 (r), r, and various constants.
The factors explicitly shown are identically nonzero, since A and B are nonzero for the deformation to be invertible, r > 0, and the other factor only vanishes if both M 12 (r) and M 13 (r) vanish, in which case we are no longer on the anomalous solution branch. Therefore, we take (...) = 0. This factor is massive, being approximately 8000 terms, so it it far too large to list here, but enough information has been provided to compute it explicitly if the reader desires. We next take its coefficients to vanish independently, and factor each coefficient. The shortest of these factors is which can be satisfied in one of two ways. Either both µ 2 and γ 1 are zero, or ν = B 2 µ2 γ1 is an eigenvalue of b. In the first case, after simplification of the other coefficients, we obtain another equation which implies either α 2 = β 2 = 0, or ν = β2 A 2 B 2 α2 is an eigenvalue of b. Taking α 2 = β 2 = 0, we obtain another similar eigenvalue equation that implies γ 2 = µ 1 = 0 or ν = B 2 µ1 γ2 is an eigenvalue of b.
Backing up a branch, we can take ν = B 2 µ1 γ2 as an eigenvalue of b. We then perform the substitutions I 1 = e 1 + 1 e2 and I 2 = e1 e2 + e 2 with ν = 1 e2 , which expresses the invariants of b in terms of the elementary symmetric polynomials in the other two eigenvalues. This reveals an equation with β 1 − A 2 B 2 α1−ACγ2 e2 as a factor. If this factor is zero, we satisfy all of the necessary equations. If this factor is not zero, we have either e 3 2 − e 2 e 1 + 1 = 0, or α 1 (e 2 β 1 + ACγ 2 ) = 0. In the later case, plugging in α 1 = 0 we obtain , which corresponds to the vanishing of the other factor. Likewise, if we take e 2 β 1 +ACγ 2 = 0, we obtain α 1 = 0 as a condition. Both of these together however imply that β 1 − A 2 B 2 α1−ACγ2 e2 = 0, which is a contradiction.
If e 3 2 − e 1 e 2 + 1 = 0, this implies that λ 1 = 1 λ 2 2 or λ 2 = 1 λ 2 1 . In either case we can express the remaining equation in terms of only one remaining eigenvalue. This equation has one factor that we know is nonzero because it corresponds to β 1 = A 2 B 2 α1−ACγ2 e2 , which would yield a contradiction. So we take the remaining factor to vanish. This factor is quadratic in α 1 . Taking the discriminant of this equation in α 1 , we obtain ∆ α1 = −4A 6 B 6 γ 4 2 λ 4 a , where λ a is the repeated eigenvalue. This discriminant must be non-negative in order for the factor to vanish with real values of α 1 . However, this discriminant is identically non-positive, which means it must be zero. However, the only way for this to happen would be for γ 2 = 0, which is a contradiction.
Having exhausted the options corresponding to α 2 = β 2 = 0, we consider ν = β2 A 2 B 2 α2 as an eigenvalue of b, and perform the substitutions on the invariants to express the invariants in terms of elementary symmetric polynomials in λ 1 and λ 2 . Doing this gives five remaining polynomial equations, which are still rather long and complicated. Ordering these equations by their length, and taking the second shortest one, we note that this equation is quadratic in e 1 . Taking the discriminant of this equation in e 1 , and demanding it be non-negative, we obtain This quantity is identically non-positive, and so the only way we can have solutions with real values for e 1 is if this quantity is zero. There are four factors that can possibly be zero: α 1 , β 1 + AB 2 Cµ 1 , A 2 B 2 α 1 − e 2 β 1 − AB 2 Ce 2 µ 1 , and AB 2 Cα 2 β 1 − β 1 γ 2 + A 2 B 4 α 1 µ 1 + A 2 B 4 C 2 α 2 µ 1 − AB 2 Cγ 2 µ 1 . First consider α 1 = 0. Inserting this, many of the remaining equations have a factor β 1 + AB 2 Cµ 1 . If this factor vanishes, the remaining equations both contain the factor AB 2 Cα 2 − γ 2 + B 2 e 2 µ 1 . The vanishing of this factor satisfies all of the equations. If this factor does not vanish, we can take the resultant of the remaining factors in e 1 , and obtain A 2 B 6 α 2 2 AB 2 Cα 2 − γ 2 2 µ 2 1 AB 2 Cα 2 − γ 2 + B 2 e 2 µ 1 2 = 0.

e2
. This second option satisfies the remaining equation, so we then consider α 2 = 0. With this, we obtain β 1 = A 2 B 2 α1−ACγ2 e2 , which is a special case of the previous option.
We then consider e 3 2 − e 1 e 2 + 1 = 0, which demands a repeated eigenvalue. This then demands that this repeated eigenvalue λ a = 1, which means that all eigenvalues are the same. This allows us to solve for β 1 and µ 1 as This satisfies the remaining conditions, and so our analysis is complete, having exhausted all possible branches of solutions.
If α 1 = 0, we solve this expression for β 2 and obtain the necessary and sufficient condition with If α 1 = 0, we can solve for µ 2 to obtain This requires Alternatively, if Cα 1 + Aγ 1 = 0, we obtain γ 1 = −Cα1 A , which implies α 1 µ 1 = 0. If we assume α 1 = 0, and insert this relation into the equations, we obtain µ 2 1 + (Aβ 1 − Cµ 1 ) 2 = 0, so we can freely take µ 1 = 0. With this, we obtain an eigenvalue equation that demands that either ν = − β1 B 2 α1 is an eigenvalue of b, or that β 1 = α 1 = 0. In the first case, we perform the usual substitutions and take discriminants in e 1 and demand non-negativity, which yields or not, in which case In the first case, we have the sufficient condition or not, in which case e 3 2 − e 1 e 2 + 1 = 0. With this, we then obtain the condition λ 1 = λ 2 = 1, and then requiring discriminants in β 2 to be non-negative, Cα 2 + Aγ 2 = 0. With this, we take γ 2 = − C A α 2 , which requires β 2 = −B 2 α 2 , which is a special case of the previous solution.
Taking γ 2 = 0 shows that in this case the sufficient condition is also necessary, and taking λ 1 = 1 λ 2 2 reveals that the sufficient condition is necessary in all cases.
Finally, we take α 2 = β 2 = 0. This reveals only one remaining equation E 6 µ 3 2 −I 1 γ 2 E 4 µ 2 2 +I 2 γ 2 2 E 2 µ 2 − γ 3 2 = 0, which, since γ 2 = 0, requires the eigenvalue ν = E 2 µ2 γ2 , and hence µ 2 = γ2 E 2 e2 . The analysis tree is then yields U 4 . Of course, these families are not simply symmetric with respect to an arbitrary choice of generators of these natures; the translational generator in U 3 is orthogonal to the plane of rotation determined by the rotational generator of U 3 , and both of the rotational generators for U 4 fix a common point. It is then natural to ask if there are other universal solutions that are likewise equivariant with respect to a similar subgroup, but without these specific generator choices.
We first want to choose our coordinates in such a way that to simplify our calculations. We seek to align our coordinate frame with the axial vector of the skew symmetric submatrix Ω. The axial vector of Ω lies in the null space of Ω, which is spanned by the vector δ, , ζ T , unless Ω = 0, in which case we do not have to do anything at this step. These two options are exhaustive, since the eigenvalues of Ω are {0, ± −δ 2 − 2 − ζ 2 }.
Provided Ω = 0, we can choose a Cartesian coordinate system such that e 3 is the normalized axial vector: We do this by considering any rotation mapping the normalized axial vector to e 3 . Denoting such a rotation R, we change coordinates by computing When we apply this coordinate transformation, our chosen element of the Lie algebra takes the form where ω = δ 2 + 2 + ζ 2 , and the α, β, and γ here have been relabeled, being independent linear combinations depending on R of the old α, β, and γ, which were arbitrary to begin with. Next, we seek to apply a coordinate translation to simplify the translation portion of our chosen element. To do this, we seek to identify the fixed points of this action. The velocity of points under the action of the one-parameter subalgebra generated by this element is given by Hence, if ω = 0, we can choose a coordinate translation that sets the point −β/ω α/ω 0 T to be the Here, u = γ, but we shall explicitly use u and ω to emphasize that we have expressed this element of se(3) in a coaxial coordinate system. In the case where Ω = 0, we simply choose our coordinate rotation so that our translation vector is aligned with e 3 , which sets our chosen Lie algebra element to the form above with ω = 0.

C.2.1 2-Dimensional Subalgebras
Obviously, provided that the generators we select are linearly independent, they span a two-dimensional vector space, hence all subalgebras containing them are at least 2-dimensional. In order for us to identify 2-dimensional subalgebras, we simply need to establish necessary and sufficient conditions for the two generators and their bracket to be linearly dependent. We select a coordinate system that is coaxial with one of our generators, and hence have and select another arbitrary generator, Taking the Lie bracket of these two elements, we obtain We then require the Lie bracket of our generators to be within their span, i.e., we seek all solutions to the equations which explicitly become a 2 α − u − βω = 0, a 2 β + uδ + αω = 0, a 1 u + a 2 γ = 0, a 2 δ − ω = 0, a 2 + δω = 0, a 2 ζ + a 1 ω = 0.
Since ω = 0, we require α = β = 0. This leaves us with v 1 as initially specified and which means that v 1 and v 2 generate independent screw motions about the same axis, corresponding to the symmetry of family U 3 .

C.2.2 3-Dimensional Subalgebras
Defining v 3 = [v 1 , v 2 ], and provided v 3 = 0, v 1 , v 2 , and v 3 span a three-dimensional vector space. The span of these three vectors must be closed under the Lie bracket, hence we require and We know that if v 3 = 0 then v 1 and v 2 generate a 2-dimensional subalgebra, hence we can freely assume v 3 = 0. First, we recognize that if ω = 0, both v 1 and v 3 are pure translations. They are linearly independent provided δ = 0 and = 0, in which case v 3 = 0. Hence, the bracket of a pure translation with a non-coaxial rotation yields another translation that is linearly independent of the original translation. Hence, t(2) is contained in such a Lie subalgebra, and hence all universal solutions that are symmetric with respect to the subgroups corresponding to these subalgebras are contained in U 2 . If the rotation is coaxial with the translation, then the bracket vanishes and we are reduced to the already solved 2-dimensional case; hence from now on, we can safely assume ω = 0.
Taking the linear combination δ a 2 + a 3 δω − ω 2 = 0 − a 2 δ − a 3 ω − δω 2 = 0 = a 3 δ 2 + 2 ω = 0, coupled with the condition ω = 0 yields either a 3 = 0 or both δ = 0 and = 0. If δ = = 0, v 3 is a pure translation that is orthogonal to the axis of v 1 , hence, taking [v 1 , v 3 ] generates another pure translation orthogonal to the axis of v 1 and that of v 3 , hence we capture the symmetry t(2) as a subgroup of our symmetry group, and hence this case is captured in family U 2 . If either = 0 or δ = 0, we have a 3 = 0, which upon substitution yields δ a 2 − ω 2 = 0 and a 2 − ω 2 = 0. These equations together imply a 2 = ω 2 . Substituting this new relation into our equations, two of our equations reduce to − 2uδω = 0, −2u ω = 0, which together imply that u = 0, since δ and cannot simultaneously vanish and ω = 0; hence v 1 must be a pure rotation, not simply a screw motion. With this, our first equation becomes γω 2 = 0, hence γ = 0 as well. When we insert this relation into our equations, we obtain 2 (αδ + β ) ω = 0, which implies that the inner product of v 2 's axial vector with its translation vector is zero. This implies that v 2 is also a pure rotation, since this inner product is unchanged under coordinate transformations. This can be seen by noting that the velocity field u induced by the action of an element of se(3) is given by  Taking the inner product of this with the embedding of the axial vector δ ζ 1 T yields αδ + β + γζ, which, not depending on position, is an invariant of the velocity field. Since the velocity field is coordinate independent, we know that this invariant will be preserved under coordinate changes. When we express our generator in a coordinate system aligned with its axis, this invariant becomes ωu, which vanishes if either our generator is a pure translation or a pure rotation. In our analysis, we have for v 2 , αδ + β = 0 together with γ = 0, hence we know v 2 is either a pure translation or a pure rotation. We know that v 2 is not a pure translation since either δ or is nonzero. Additionally, we know that the axial vectors of v 1 and v 2 are linearly independent, since either δ or is nonzero. Summing up our progress thus far, we have shown that both v 1 and v 2 must be pure rotations. In fact, their axes of rotation intersect, hence they generate so (3), the symmetry present in U 4 , indicating that our classification captures all 3-dimensional cases. To see this, note that we have aligned our coordinates so that the axis of rotation for v 1 is the z axis. We seek to show that the axis of rotation for v 2 intersects the z axis.
First notice that for rotations about the origin, the velocity field generated is of the form where ω is the axial vector of Ω, and ∧ is the standard cross product. This implies that the velocity vector at a point is orthogonal to the plane spanned by the axial vector of the rotation, and the position vector X. Since (446) assumes we have chosen our origin such that the axis of rotation passes through the origin, this plane is equivalently the plane containing the axis of rotation and the point X. Therefore, for the generator v 2 , we can examine the velocity generated at the origin, and recognize that it lies entirely in the x, y plane. If this velocity is nonzero, we know that the plane passing through the origin that is orthogonal to this translation contains v 2 's axis of rotation. This plane also contains the z axis, since all planes passing through the origin that are orthogonal to a nonzero vector in the x, y plane contain the z axis. Therefore the axes of rotation for v 1 and v 2 are coplanar. We have already established that they are not parallel, since the axial vectors for v 1 and v 2 are linearly independent, hence they must intersect at some point. If the velocity generated by v 2 at the origin is zero, then the axis of rotation of v 2 passes through the origin, and hence not only intersects the z axis, but intersects it at the origin. We have therefore shown that all three-dimensional Lie subalgebras of se(3) that are generated by two linearly independent generators either contain t(2) as a subalgebra, or are so(3), the algebra associated with the set of rotations about a fixed point, and hence universal solutions that are equivariant with respect to the associated Lie groups are contained in one of our discovered families.

C.2.3 4+ Dimensional Subalgebras
Without loss of generality, we assume v 1 , v 2 , v 3 , and v 4 = [v 1 , v 3 ] are linearly independent, since the other choice would be v 4 = [v 2 , v 3 ], which would be equivalent. Specifically, we denote V 2 = Span (v 1 , v 2 ), and V 3 = Span (v 1 , v 2 , [v 1 , v 2 ]). Provided that v 1 , v 2 , and [v 1 , v 2 ] are linearly independent, we can write V 3 = V 2 ⊕ Span ([v 1 , v 2 ]). It suffices to take the fourth linearly independent element to be of the form since for all u, w ∈ V 2 , [u, w] ∈ V 3 , and for all u, w ∈ Span ([v 1 , v 2 ]), [u, w] = 0. Since v 1 and v 2 are arbitrary, we can choose this fourth linearly independent element to be [v 1 , v 3 ]. Doing this, we have Notice that the axial vectors of v 1 , v 3 , and v 4 are [0, 0, ω] T , [− ω, δω, 0] T , and −δω 2 , − ω 2 , 0 T respectively. These vectors are mutually orthogonal, hence provided ω = 0 and that = 0 or δ = 0, these span R 3 , and hence the rotational components of these three generators can be used to reduce any fourth linearly independent generator to a pure translation. As shown earlier, taking the bracket of a pure translation with any other linearly independent element of se(3) generates a 2-dimensional subalgebra: either t(2) or so(2) × t(1). Therefore, all subalgebras of dimension four or higher contain one of these two-dimensional subalgebras, hence universal solutions that are symmetric with respect to such a 4-dimensional subalgebra will be contained in either U 2 or U 3 .