Architected Elastomer Networks for Optimal Electromechanical Response

Dielectric elastomers (DEs) that couple deformation and electrostatics have the potential for use in soft sensors and actuators with applications ranging from robotic, biomedical, energy, aerospace and automotive technologies. However, currently available DEs are limited by weak electromechanical coupling and require large electric fields for significant actuation. In this work, a statistical mechanics-based model of DE chains is applied to elucidate the role of a polymer network architecture in the performance of the bulk material. Given a polymer network composed of chains that are cross-linked, the paper examines the role of cross-link density, orientational density of chains, and other network parameters in determining the material properties of interest including elastic modulus, electrical susceptibility, and the electromechanical coupling. From this analysis, a practical strategy is presented to increase the deformation and usable work derived from (anisotropic) dielectric elastomer actuators by as much as $75-100\%$.

The key advantages of dielectric elastomers (DEs) are that they are generally inexpensive, lightweight, easily shaped, pliable, and can undergo significant deformations [BC04]. However, despite these advantages, DEs are limited by weak electromechanical coupling so that large voltages are often required to achieve meaningful actuation [BC04]. A better understanding of these materials -in particular, how mesoscale polymer network characteristics affect the bulk material response -would not only enable microstructural design with currently available DEs, but can also potentially lead to the development of advanced DEs with stronger electromechanical coupling, as discussed in this paper.
Manufacturing technologies are rapidly improving, thereby enabling the synthesis of materials with complex and hierarchical patterning, e.g. [ABB + 17, FAK + 19]. Thus, it is the goal of this work to: (1.) develop a model, based in statistical mechanics, to connect the molecular-and macromolecular-properties of a given DE to its performance as an actuator and/or energy harvester, and to (2.) use the multiscale model to develop insights and predictions into how novel DEs with enhanced electromechanical coupling may be designed and manufactured, particularly by tailoring the polymer network architecture. Broadly, our strategy and goal is similar to work in nonlinear elasticity-based homogenization (e.g. [LP14,SC14,CG11,GC13]), in that we aim to tailor the microstructure to achieve desired or optimal responses, except that we focus on smaller lengthscales and use physical models that are appropriate to these scales.
The quintessential example of a DE actuator (DEA) is a thin DE film sandwiched between two compliant electrodes-a soft parallel plate capacitor. When a voltage difference is applied across the electrodes, a positive charge density accumulates on one of the electrodes and an equal and opposite charge density accumulates on the other. The electrodes are attracted to each through simple Coulombic forces, causing the DE film to compress in the thickness direction and, because DEs are roughly incompressible, the expands in the plane of the electrodes [PKK00,WM07,Kof08,KZSK12].
Although the Coulombic attraction is an important factor, it has been pointed out (both theoretically and experimentally) that if the susceptibility of the DE film is a function of deformation, then an additional stress develops in the DE film [ZS08,Suo10,Cd16]. In this work, we aim to design novel DEs that (1.) have maximal susceptibility (in their typical operating conditions) in order to increase the capacitance (and, therefore, performance) of its corresponding thin film DEA and (2.) use this additional stress in order to increase the electromechanical coupling of DEAs.
The electromechanical modeling of DEs can be grouped into two categories: (1.) continuum based approaches where the general form of the energy density (of either a polymer chain or polymer network) as a function of mechanical and electrical loads is inferred by assuming a form of the equation, which usually depends on electroelastic invariants, that leads to behavior observed in experiments (see for example [Tou56, DO06, ZS08, Tou56, DO14, HCB13, ZDDP17, LS18, KLS19, LLS15, DDLS19, Sd13, FG08, RYBS19]) and (2.) statistical mechanics based approaches that build from molecular-scale responses up to the levels of a single chain and eventually the continuum level response.
Statistical mechanics has a rich history in modeling and understanding rubber elasticity, e.g. [Tre75,KG42,Wei12,Flo44,Hil86]; however, its use in the modeling of DEs is much less developed. The work by deBotton and coworkers [CDd16,Cd16] appears to be the first in this direction 1 , and an explicit approximate expression for the response of a single electro-active polymer chain was provided in [GD20]. This explicit approximation provides the starting point for the analysis in the current paper.
Structure of the paper. In Section 2, we provide a summary of the statistical mechanical model of a single electro-responsive polymer chain from [GD20], and the related issues of averaging over chains to obtain the continuum response. In Section 3, we identify certain network properties as design parameters, motivated heuristically by the ease of being able to tune these parameters experimentally. Subsequent sections (4,5,6) examine the role of network properties on the elastic modulus, the susceptibility, and the electromechanical coupling modulus respectively. Finally, Section 7 provides an application of these ideas to the setting of a simple DEA geometry. Appendix D provides a summary of the key notation that is used in the paper.
carried by the i-th monomer, denoted µ i ; and (4.) the electric field, -∇φ(x), which, for now, is a general function of position x. All of these can be varied independently.
We make the following assumptions in our model at this stage. (1.) Following standard practice, we have implicitly assumed above that only the orientation of the monomers is relevant and that there is no stretching. That is, stretching of monomers costs energy that is much larger than kT , and hence the chain is assumed to be inextensible. (2.) Again following standard practice, we assume for simplicity that the bending energy -i.e., energy associated with the change in orientation of the monomers along the chain -is much less than kT and can hence be neglected. (3.) While dipole effects are due to electronic and nuclear motion and hence can have interactions between monomers [BJDM17], we assume for simplicity that the dipole induced in a monomer can be modeled without regard to the environment; equivalently, we assume that the electrical energy of the monomers can be decomposed additively.
Under these assumptions, the potential energy of a microstate can be written: as the local electrical field at the location of the monomer, and where we have chosen to work in Gaussian units (i.e. unit system in which ǫ 0 is unity). The first term in the summation is the energyũ required to separate charges to form a dipole µ i , and will be discussed further in Section 2.A.1. The second term in the summation is the energy of interaction between the local electric field at the monomer location and the induced dipole. The volume integral is the electrical field energy.

2.A.1. Dipole Response of a Single Monomer to an Electric Field
Following the classical Born-Oppenheimer approximation, we assume that electrons reach their ground state configuration -under electric field -rapidly compared to the timescale of the thermal motion of the atoms [YM16]. While we do not consider quantum effects explicitly, the Born-Oppenheimer approximation justifies neglecting thermal effects in modeling the dipole response of a single monomer to an electric field. In other words, we are assuming that the first excited state of the electrons has energy much larger than kT , and the system is always in the ground state with respect to the electron configuration.
An important implication of this assumption is that the dipole moment of an individual monomer, µ i , is uniquely determined -through energy minimization -givenn i , E i , and x i . To find an expression for µ i in terms of the other quantities, we differentiate (2.1) with respect to µ i to obtain the polarization response ∂ũ/∂µ i = E i at the ground state.
Here we model the monomer response through the choice: where χ is the tensorial polarizability of the monomer. Following [CDd16], we model the tensor χ as transversely isotropic: χ(n) = χ n ⊗n + χ ⊥ (I −n ⊗n). The non-negative material constants χ and χ ⊥ are measures of the susceptibility along the monomer direction and transverse to the monomer direction respectively. We refer to monomers with χ > χ ⊥ as uniaxial, and monomers with χ < χ ⊥ as transverse isotropic (TI). For this choice of χ, we have the ground state energy, u =ũ − µ i · E i , as: If χ > χ ⊥ , i.e. the monomer is uniaxial, then the monomer has minimum energy whenn is parallel or anti-parallel to E i , and maximum energy whenn lies in the plane orthogonal to E i . If χ < χ ⊥ , i.e. the monomer is TI, the situation is reversed: the minimum energy state is whenn lies in the plane orthogonal to E i , and the maximum energy state is whenn is parallel or anti-parallel to E i .

2.A.2. Multiscale Structure of the Electrical Field Energy, and Consequent Nonlocal-to-Local Decoupling
Under the assumptions in Section 2.A.1, the potential energy from (2.1) reduces to: The energy posed in (2.4) has a highly nonlocal structure [MD14,JK90,YD11,Liu13]. Physically, this is due to the fact that we need to solve the electrostatics equation to find the field at every monomer location in Ω. We can see this by examining the ground state with respect to the electric potential. Taking the variation φ → φ + ψ and requiring this to be 0 for all variations ψ, we find that: ⇒ div ∇φ = divp, subject to boundary conditions. (2.5) Here, δ x i is the Dirac mass located at x i , andp(x) is the dipole moment of the chain, treated here as a field through the use of Dirac masses that represent the point dipoles carried by the monomers, to be interpreted in the sense of distributions.
This final equation shows the nature of the nonlocal problem: to evaluate the energy in (2.4), we need to solve a boundary value problem. Therefore, statistical mechanical averaging over the energy will need to average over all fields φ that are consistent with the specified average electric field ensemble. This is challenging and would require methods of statistical field theory.
Instead, we begin by assuming that the energy in (2.4) has a separation of scales such that the energy in forming a dipole by separating charges is of order kT , while the stored energy of the field in vacuum is of much higher order. This allows us to first minimize with respect to the field energy to find the electric field, and then use this electric field as a constraint in performing the statistical averaging.
The potential energy, with the separation of scales explicitly highlighted, is: We first minimize the field energy, of order much greater than kT , by setting the first variation to zero to obtain div ∇φ = 0. We find that −∇φ = E 0 , where E 0 is a fixed (average) electric field specified by the ensemble (see [GD20] for a more detailed discussion).
The final form of the reduced potential energy that we will use for statistical averaging reads simply: For simplicity, moving forward, we will drop the subscript zero and let it be understood that E is a locally average electric field at a chain or a material point.
We highlight two additional important simplifications that are a consequence of the assumption of separation of energy scales. First, the degrees of freedom required to describe a microstate are vastly reduced from our general starting point. In particular, the spatial position plays no role, as the electric field is uniform in space; the point dipole is completely specified given E and the orientation of the monomer; and the electric field is also completely specified −∇φ = E. Therefore, the monomer orientationsn i are the only remaining degrees of freedom over which to conduct statistical averaging. Second, the simplification from having to solve a nonlocal boundary value problem to determine ∇φ makes it possible to perform a Legendre transform to go between the free energy written as a function of E, which is relatively simpler to evaluate using statistical mechanics, and the free energy written as a function of p, which is more convenient for applications since it has a minimum rather than a saddle-point structure [GD20]. We will see later on that this has implications regarding Legendre transforms of the continuum-scale free energy density as well.

2.B. Statistical Mechanics of a Single Electro-responsive Polymer Chain
Given the potential energy of a microstate, we derive the free energy in [GD20]. To this end, we derive a mean-field theory and determine that the density of monomers oriented in the directionn is given by where κ = E 2 ∆χ/2kT , and the unknowns, C and τ , are determined by enforcing the constraints where S 2 denotes the surface of the unit sphere, n is the number of monomers in the chain, r is the endto-end vector, and b is the length of a single monomer. In terms of the monomer density function, one can show that the free energy is approximately: There are still two remaining difficulties in solving for C and τ : (1.) the integrals in (2.9) are difficult to evaluate and (2.) the resulting systems of equations are nonlinear. However, there are two limits in which a closed-form solution can be obtained. Let γ = r/nb denote the absolute chain stretch, where r = |r|. In the limit of γ → 0, we have that |τ | → 0. Thus, one can obtain a solution in the limit of small stretch by using a Taylor expansion on the τ ·n term in the exponential of the monomer density function (ρ given by (2.8)) up to linear order. Similarly, ρ takes a simpler form when γ → 1. Physically, this is because when the chain approaches the fully stretched limit, the distribution of possible monomer orientations are narrowly centered aroundr := r/r; and in this case, all of the possible orientations are of approximately the same energy. Hence the energy term, −κ Ê ·n 2 , in the exponential of ρ can be neglected. Neglecting the energy term, one obtains the Kuhn and Grün [KG42] solution. Let F sτ denote the approximation of the free energy using the small stretch solution and similarly, let F KG denote the approximation near the fully stretched limit (i.e. the result of substituting the Kuhn and Grün density into (2.10)). Then, using what is known about the limiting behavior, we construct a free energy approximation: where κ ⊥ = E 2 χ ⊥ /2kT , and L −1 is the inverse Langevin function. Notice that (2.11) recovers the exact solution when κ = 0 and is exact in the limits of zero stretch and full stretch. This approximation has been shown to agree well with numerical experiments for a large variety of general chain conditions (e.g. stretches, orientations with respect to the electric field, |κ|, etc.) [GD20].
Having obtained an approximation of the free energy, one can obtain the chain polarization, p, by differentiating the free energy with respect to the electric field, which is equivalent to integrating the dipole moment over the monomer density function; i.e.
The significance of the relationship between the partial derivative, ∂F ∂E and the chain polarization is that it establishes the fact that F is the Legendre transform of A = A (r, p), that is, the Helmholtz free energy of the chain. The derivation of this relationship is given in Appendix A.
We note an important feature of the free energy-stretch relationship of electro-responsive chains: numerical experiments suggest that the F/kT vs γ curve is convex and its minimum is at zero stretch, implying that a chain will not stretch when subject to electric field. Therefore, the additional electrostriction that occurs in dielectric elastomers when the permittivity is deformation-dependent cannot be explained by the notion that individual chains spontaneously stretch or contract within the network due the applied field 2 . Physically, this feature of the F/kT -γ relationship can be understood as a consequence of: (1.) the electrostatic energy of the monomers being quadratic in E ·n (see (2.3)) and (2.) the assumption that monomer-monomer interactions are negligible. This means that if a monomer's orientation is reversed (i.e.n → −n), its energy, and hence, its contribution to the Boltzmann factor is the same. And since there is no energy penalty associated with large or small bond angles between neighboring monomers, the chain is free to fold back on itself. So in terms of the Boltzmann factor, a longer end-to-end vector is never any more favorable than a shorter end-to-end vector. However, in terms of entropy, the shorter end-to-end vector is more favorable. For these reasons, the free energy versus stretch relationship for an electro-responsive polymer chain is expected to be convex with its minimum at zero stretch.

2.C. Orientational Averaging Over a Polymer Chain Network
In order to relate the continuum scale deformation to individual chains in the network, we use the affine deformation assumption [Tre75] and the full network model [WVDG93,Bea03].
At each material point in the reference (stress-free) state, it is assumed that chains have their most probable 2 By additional electrostriction, we mean the contribution to electromechanical coupling that is not due to the Coulomb attraction between the electrodes. length, b √ n, and are randomly oriented. For instance, consider a material that is isotropic due to the chain orientations being uniformly distributed. That is, the probability density function of finding a chain with end-to-end vector,r, is 4π .
By the affine deformation assumption, each chain gets mapped from the reference configuration to the current configuration by the deformation gradient, F , i.e. r = Fr. Finally, the free energy density is taken to be the product of the average chain free energy and the number of chains per unit volume, N : where N is the number of chains per unit volume and · r denotes an average over the distribution of chains.
Similarly, we can obtain the (continuum-scale) polarization, i.e. dipole moments per unit volume, P by: Appendix C discusses an approximate closed-form expression to evaluate the integration in (2.12) for a commonly-used class of P. In general, numerical integration of (2.12) can be challenging to perform accurately [Ver15].
Let Ω 0 denote the body in the reference configuration and Ω = Υ (Ω 0 ) the body in the current configuration, where Υ is the deformation mapping. The position of a material point in the reference configuration, X, is mapped to a position in the current configuration by the deformation mapping, i.e. x = Υ (X); and the deformation gradient is given by F = Grad Υ .
Then given some boundary conditions, B, and some external work, W, the equilibrium configuration of the body is given by the minimization of the Gibbs free energy such that the boundary conditions are satisfied; that is: Thus, the constitutive response is encoded in the form of the Helmholtz free energy density function, W, which is given by whereP = JP is the pullback of the (continuum-scale) polarization 3 , and W * is the Legendre transform (in theP slot) of the Helmholtz free energy density. We will denote W * as the closed-dielectric free energy; more precisely, it is the Gibbs free energy when the dielectric material can be considered as a closed system in a constant electric field. Physically, the significance of the closed-dielectric free energy is that it has a minimum principle at constant temperature and constant electric field. See [GD20] Section 3.3 for further discussion of this issue.
Finally, recall that because of the long-range nature of dipole-dipole interactions, the relationship between E andP is nonlocal; they are related by the electrostatic equation. Therefore, in general, the Legendre transform cannot be taken in terms of the free energy density, rather this must be done at the level of the total free energy, which is a formidable challenge as well as highly problem-specific. However, when the monomer-monomer dipolar interactions are neglected, the nonlocal relationship becomes local (recall Section 2.A.2 of this paper, or see Section 2.2 of [GD20]) and the transform taken in (2.14) is indeed valid.

2.E. Electromechanical Coupling Due to Chain Torque
We highlight some important physical implications of the model outlined above. We notice that an applied electric field applies a moment to individual polymer chains, due to the fact that the chain carries an effective dipole. The importance of this is two-fold: (1.) it contributes to the stress beyond that imposed by the Coulombic attraction of the electrodes, and (2.) it provides a means to engineer the orientational distribution of chains so as to obtain a desired network architecture.
(1.) represents an opportunity to design for enhanced electromechanical coupling. Specifically, one could potentially tune the network architecture such that the deformation is increased due to contributions from (1.). To illustrate (1.), we consider a typical DEA, but with the stress caused by the electrodes counteracted. This counteraction could be achieved, for instance, by applying a traction, t  , to the outside surfaces of the electrodes that is equal and opposite to the Coulomb attraction and by ensuring the voltage difference is adjusted to keep the electric field constant when the distance between the electrodes changes. This is shown in fig. 2.1 (inset).
In this case, the change in energy stored in the electric field (or equivalently, the work of the battery) and the work of the applied traction will cancel each other out; as a result, the equilibrium configuration will be the one that minimizes the free energy of the DE film. We neglect the fringe fields and assume that the electric field is uniform in space, implying that the electric field is then aligned along the thickness direction. Further, assume homogeneity, isotropy, and incompressibility; the stretch in the thickness direction is denoted λ E , implying that the stretch in any in-plane direction is λ − 1 2 E . Evaluating the free energy for this class of deformations using (2.14) with uniform orientational distribution P = 1/4π, we then minimize to find the energy-minimizing stretch λ E . FIG. 2.1 shows λ E as a function of κ. We see that the model predicts a spontaneous deformation of the film, even in the absence of pressure from the top and bottom electrodes. When the chain is composed of TI monomers (κ > 0), the film shrinks in the thickness direction (i.e. λ E < 1). Alternatively, when the chain is composed of uniaxial monomers (κ < 0), the film elongates in the thickness direction (λ E > 1).

Design Parameters
In this section, we identify design parameters which describe the polymer network architecture. That the design parameters constitute a description of the network is supported by the affine deformation assumptionwhich states that chain end-to-end vectors in the reference configuration are mapped to the current configuration by the deformation gradient-and the negligibly interacting assumption-which allows us to decompose complex chain architectures in the network into separate linear parts.
The key variables that we aim to optimize over are related to the network connectivity and geometry, and the orientational distribution of the chains. These aspects are relatively feasible to control using applied electromagnetic fields or stresses [ABB + 17]. In contrast, we treat the dielectric response, specified here by χ , χ ⊥ of individual monomers as a fixed and given, as modifying these characteristics require molecularscale chemical methods that are much less developed in terms of being able to achieve target properties.
The following is a list of the design variables: mass density, i.e. N 0 n where N 0 = JN is the number of chains per unit volume in the reference configuration; fraction of loose-end monomers, α; density of cross-links; and orientational distribution of chains, P; which will be introduced shortly.

3.A. Mass Density
Since we take χ and χ ⊥ to be fixed and, for many applications, we would prefer a high polarization susceptibility, we may be tempted to pack as many monomers per unit volume, N 0 n, into the polymer network as possible. There are, however, some subtleties associated with this strategy. For one, since there is a mass associated with each monomer, this would serve to increase the mass density of the material; however, it is likely that we do not want the density of the material to exceed a certain threshold for many of the applications of interest (e.g. wearable sensors, soft robotics, etc.). Similarly, packing more monomers into the network may also affect the compliance of the material. Therefore we will be interested in the product N 0 n as a tunable parameter; and, specifically, its influence on, and interplay between, the aforementioned bulk properties.
Here H A , H B , are the Hamiltonians of systems A, B respectively; and H A↔B , is a correction term due to interactions between system A and system B (see [TM11] section 7.3.3, for example). In this work, we assume that the chains of the network are in weak interaction with each other such that their interactions are negligible. Because H A↔B is negligible, we neglect it completely. As we will see, in the negligible interacting regime, one can characterize many different chain and network architectures in a general and straightforward way.

3.C. Individual Chain Geometry: Backbone Monomers, Loose-end Monomers, and Linear Chains
A key design parameter is the architecture of the chains that make up the cross-linked network. Specifically, we mean the geometrical features of the chain, the types of monomers the chain consists of, and the various sequences that the monomer types are arranged in. FIG. 3.1 shows a few examples of polymer chain architectures. We divide the monomers into two groups: (1.) n b monomers that contribute to the backbone of the chain; and (2.) n l monomers that are in the loose ends of the network 4 . By backbone of the chain, we mean all segments that span from cross-link to cross-link; by loose ends, we mean the non-backbone segments. This distinction allows us to consider monomers which interact through the end-to-end vector constraint separately from those that do not need to satisfy such a constraint. To make this clear, consider the example architectures in fig. 3.1. If there are cross-links at each of the stars in fig. 3.1.a, then we say this is a linear chain without any loose-end monomers (i.e. n b = 11, n l = 0). In contrast, if a cross-link is missing at either star, then the chain only contributes to the average number of loose-end monomers in the network, (i.e. n b = 0, n l = 11). Similarly, in the branched polymer given in (b), if there are cross-links at the stars but not the squares, then our decomposition says that this is a single chain with n b = 11, n l = 8. In the last example, the star polymer, (c), if there are only cross-links at the stars, this results in a single chain with n b = 22 and n l = 33; whereas cross-links at both the stars and squares results in five linear chains, each with n b = 11.
Having decomposed complex architectures into linear parts, we make the assumption that loose end monomers are in negligible interaction with the backbone monomers and that separate chains are in negligible interaction with each other 5 . This reduces all chain architectures to collections of negligibly interacting linear chains with some fraction of loose ends, α = n l /n. This is justified as follows. First, in regards to the electrical energy of the system, we have assumed that all dipole-dipole interactions will be neglected. Second, for the entropic contributions, monomers within a chain interact through the enforcement of the end-to-end vector constraint; however, loose end monomers, do not take part in this constraint, and obviously monomers in separate chains do not interact through such a constraint. The negligible-interaction assumption could be violated if chains are "too short" such that a significant portion of monomers in the chain are in the vicinity of another chain; if long-range interactions are significant; or if excluded volume effects are important.
To simplify further, we assume that the chain architectures are uniform enough to model their behaviour in terms of average numbers of backbone monomers and loose-end monomers per chain 6 . In summary, we describe the network through the number of linear chains per unit volume, N , and the (average) number of backbone monomers and loose-end monomers per chain. Thus, N n b is the number of backbone chains per unit volume; in other words, the number of monomers per unit volume which are related to continuum scale deformation through F and chain averaging. Whereas, N n l is the number of loose-end chains per unit volume; in other words, the number of monomers which only contribute to the polarization of the network and do not contribute at all to the elasticity of the network.

3.D. Density of Cross-links
The junction point at which chains in the network are joined together are called (effective) cross-links. Specifically, given our decomposition of complex chain geometries into chains with linear backbones and loose-ends, we say that there are cross-links at the beginning and end of each linear backbone. For a given mass density then, N 0 is proportional to, and n is inversely proportional to, the number of cross-links per unit volume. Therefore the density of cross-links in the network is both a parameter which is known to be controllable [Tre75] and will affect the electroelasticity of the architected network (by (2.11) and (2.15), for example).

3.E. Chain Probability Distribution Function
The orientation of polymer chains with respect to the local electric field influences the electroelastic response of the network. As a result of both the importance of orientation and our decomposition of complex chain architectures into (multiple) linear chains, we use the orientational distribution of these linear chains within the network as a design variable. More specifically, let P = P (r) denote a probability density function which describes the fraction of chains with reference end-to-end vectorr. The design space of P is the space of all probability measures on R 3 ; and hence, an infinite dimensional space. However, we reduce the dimensionality of the problem as follows.
First, we assume that the length of all chains depends only on n and b; specifically, we assume |r| = b √ n, which is the expectation of the length of a random walk of n steps with step length b. This choice is motivated by the fact that, prior to cross-linking, chains are free to move in a random manner. While this choice is primarily motivated by classical polymer theory, it will also be discussed below when we discuss residual stresses.
Second, having reduced the support of P to the surface of a sphere of radius b √ n, we next consider the symmetries inherent in the physical problem. For instance, in regards to chain statistical mechanics: F (r, ...) = F (−r, ...); or, in other words, although chains have well defined ends (cross-linking points), either end may be identified as the start and finish. Therefore P need only be defined on a half sphere.
Third, we further reduce the dimensionality by considering the setting of a thin film dielectric elastomer actuator with the voltage applied across the thickness of the film. Physically, the direction across the film thickness is significant. We therefore restrict our attention to functions P = P (θ) which satisfy where θ is the polar angle with respect to the direction of the electric field,Ê. This choice is consistent with the symmetry of the DEA and with the r → −r symmetry. For materials with the above symmetries and the BVP associated with a thin film DEA, we expect a homogeneous deformation of the form Fourth, we obtain a drastic dimension reduction by using the following as an ansatz of P: where C is a normalization constant such that P = 1. Notice that this ansatz is consistent with a transversely isotropic material.
Some examples of this ansatz are shown in fig. 3.2. The reasons for this choice of ansatz are as follows: (1) the properties of Gaussian distributions are well understood; (2) ease of evaluating various integrations, e.g. (2.15); (3) the symmetry given by (3.1) is satisfied; (4) it reduces our search space from infinitedimensional to a two-parameter space spanned by (θ 0 , σ); and (5) it is consistent with heuristic ideas of modeling manufacturing tolerances; for instance, tolerances could be modeled by placing lower bounds on σ.
In regards to manufacturing for a design P: we envision controlling P either by cross-linking while under an applied electric field (i.e. taking advantage of chain torque) and/or applied stresses; or by an advanced type 3D printing [ABB + 17].

3.F. Reference States to Define Material Response
Consider a uniform distribution P = 1/4π. In the absence of external field and external mechanical load, assuming incompressibility gives from symmetry that there is no deformation at equilibrium, i.e. F eq = I.
Interestingly, this is not the case for a general choice of P, i.e. the material is not at equilibrium for F = I in the absence of external loads. This can be interpreted as a residual stress introduced during manufacturing. We use the relaxed Helmholtz energy minimizing deformation, denoted F ⋆ , to define linearized material response properties that we will aim to optimize.
It is for this reason that we remarked earlier that assuming fixed |r| = b √ n would be further justified. This justification lies in the fact that it is apparent that chain lengths and orientations cannot be controlled independently of each other-they are related through the relaxation of the residual stresses. Indeed, after F ⋆ the chain lengths in the network will no longer be homogeneous, but will depend onr (i.e. |F ⋆r | = const.).
Considering the Taylor expansion of the free energy density about the relaxed configuration F = F ⋆ ,P = 0 : where ǫ 1 = |F − F ⋆ | and ǫ 2 = |P |. We assume W is convex inP and its minimum isP = 0, justified empirically for general DEs. Therefore, the linear terms in (3.3) vanish. The constitutive response of the These quantities correspond to the stiffness tensor, the inverse of the polarization susceptibility tensor (i.e. X −1 )-which is a measure of the stiffness of the atomic bonding of charges bound to the DE-and the cross modulus tensor, respectively. The magnitude of the cross modulus signifies the intrinsic electromechanical coupling of the material; that is, the electromechanical coupling irrespective of external loads or Coulombic attraction of external electrodes. Thus, we will be interested in how our design variables affect these three quantities. However, we will not, strictly speaking, seek to optimize all of these properties. In particular, stability and positive definiteness of the energy near the equilibrium state bounds the cross modulus from above by, roughly, the geometric mean of the susceptibility and the stiffness. Therefore, a practical strategy is to keep the stiffness from getting too high or too low, while optimizing the susceptibility and the cross modulus.

3.G. Incompressibility
As a final note to this section, we make a few remarks regarding the incompressibility of our hypothetical, proposed, anisotropic electroresponsive elastomer. While most elastomers are incompressible to a very good approximation and are often modeled as such, this is empirically motivated and taken as an assumption in both statistical mechanics and continuum mechanics models. Since we are proposing to design an anisotropic material -which is potentially very different from existing elastomers -it is possible that incompressibility will not be a good approximation for our designed materials. However, for simplicity, we will assume here that the proposed anisotropic electro-responsive elastomers are incompressible. Formally, this means that the deformation gradient must satisfy J = 1 where J := det F . As a consequence, N = N 0 and P =P .

Elastic Modulus
We first consider the elastic modulus in the purely mechanical setting. Let E (x) = 0 for all x ∈ Ω. Then U = 0 and consequently: The affine deformation assumption gives that r = Fr, and therefore the Helmholtz free energy density is given by: wherer can be written as b √ n (cos φ sin θ, sin φ sin θ, cos θ). This integral is difficult to evaluate in closedform, but simplifies significantly if we use a Taylor expansion of L −1 (γ) in powers of γ: For the remainder of this section, we neglect the higher-order terms in (4.2). For now, we consider F of the form diag (λ 1 , λ 2 , λ 3 ); thus, First, let P = 1/4π. It is easy to show that we obtain the neo-Hookean energy density: where G iso := N kT is the shear modulus that is predicted by the Gaussian chain approximation in classical rubber elasticity [Tre75]. This result is to be expected since we took a Taylor expansion of the Langevin chain statistics (4.1) about zero stretch and because the chosen form of P is isotropic. Clearly then, for isotropic elastomers, the stiffness-for a constant mass density, N 0 n-increases with the density of crosslinks; that is, increases with N 0 (recall: N 0 = JN ) (this is well known; see, for example, [Tre75]). Moreover, the slope of the inverse Langevin function is a monotonically increasing function of its argument. Its argument in (4.1) is O n −1 , so increasing the density of cross-links for fixed mass density-which effectively lowers n-also increases the stiffness through higher order terms in (4.2). Physically, this is because the higher order terms account for the finite extensibility of the chain and, as n decreases, so does the maximum length of a chain. Similarly, the stiffness increases with the fraction of loose end monomers (for fixed n) because, as α increases, the maximum length of the chain, n b b, decreases and finite extensibility effects are more relevant for shorter stretches [Tre75, KG42] 7 .
Next we consider a transversely isotropic elastomer. In this case, we take P to be given by the ansatz (3.2). We make the definition: where k ∈ N, µ ∈ 0, π 2 , and σ ∈ [0, ∞). Then the average square relative stretch and the Helmholtz free energy density are: We can approximate these expressions in the limits of σ ≪ 1 and σ ≫ 1 by using (C1) and (C2), respectively.
As mentioned previously, there is no systematic theory which connects the molecular structure of elastomers to their effective Poisson's ratio; instead, incompressibility is generally taken as an assumption, and we use this assumption here.
We consider loading the body with displacement-control along the axis of symmetry, and free on the transverse faces. With incompressibility, the symmetry of the system gives the deformation gradient of the form where the 3-direction is aligned along the axis of symmetry.
The slope of the stress-strain curve gives the tangent elastic modulus Y in the direction of the axis of symmetry. Similarly, the tangent elastic modulus in the plane orthogonal to the axis of symmetry is denoted by Y ⊥ 8 .
In the limit of σ ≪ 1, our model predicts: which shows a strain hardening in compression and a strain softening in tension-except when θ 0 = 0 and σ = 0. Interestingly, (θ 0 , σ) = (0, 0) recovers the isotropic elastic modulus; as does any (θ 0 , σ) at λ = 1. At first glance, this seems physically unreasonable, as one would expect changing the distribution of chain directions in the network to affect the stiffness in various directions. However, recall from our discussion in Section 3.E that λ = 1 is no longer the equilibrium configuration in the absence of external loads. Indeed, let λ ⋆ be such that ∂W ∂λ λ=λ ⋆ = 0; then, Equation (4.6) is shown in FIG. 4.1 as a function of θ 0 for σ = 0, π/32 and π/16. Let θ iso be the polar angle for which F ⋆ = I when σ = 0; one can show that θ iso = arctan √ 2. When θ 0 < θ iso , we have λ ⋆ < 1. This is because there is simultaneously a higher density of chains that are oriented closer to the direction of stretch (i.e. the axis of symmetry), and a lower density of chains oriented towards the orthogonal to the direction of stretch. The elasticity of polymer chains is due entirely to entropy in the absence of electric field in our model; consequently, polymer chains are in tension, at any finite temperature, if their end-to-end vector is finite. Put differently, the maximum chain entropy is given by a vanishing end-to-end vector; thus, chains always want to contract if we neglecting excluded volume effects. The additional incompressibility assumption accounts empirically for the excluded volume effects and prevents the collapse of the network to zero volume. 8 The property Y ⊥ corresponds to the slope of the stress-strain curve if a sample were stretched biaxially in the plane orthogonal to the axis of symmetry and traction-free on the lateral faces. It can be calculated from the deformation considered here by When incompressibility is enforced, for chains to contract in one direction requires chains to elongate in other directions. This trade-off between stretch and elongation in different directions gives, for isotropic networks, that λ ⋆ = 1. When θ 0 < θ iso , we obtain λ ⋆ < 1 and thus there is a contraction along the axis of symmetry, because there are more chains oriented toward this direction and hence a net increase in entropy can be gained from some λ ⋆ < 1 compared to the decrease in entropy of fewer chains that must elongate in other directions. Conversely, when θ 0 > θ iso , there are fewer chains along the axis of symmetry, and a net increase in entropy can be gained by contracting orthogonal to the axis despite the decrease in entropy to the fewer chains along the symmetry axis that must elongate. In the limit of σ = 0 when all the chains are aligned precisely along a single direction, this leads to singularities at (θ 0 = 0, σ = 0) and (θ 0 = π/2, σ = 0) as can be seen in FIG. 4.1.
We return our attention to the elastic modulus. In particular, we are interested in Y ⋆ := Y (λ = λ ⋆ ); that is, the stiffness of the material at its equilibrium state when no loads are applied. Using (4.6) in (4.5), we obtain: Similarly, As expected, the maximum Y ⋆ occurs at (θ 0 , σ) = (0, 0)-as this is the case which has the maximum amount of chains oriented in the direction of the axis of symmetry for a given N 0 . Interestingly, although, as previously mentioned, there is a singularity at (θ 0 , σ) = (0, 0) such that λ ⋆ = 0, the elastic modulus is finite. The maximum Y ⋆ is given by 9N kT , which is the number of spatial dimensions times the elastic modulus for an isotropic network as expected, Y ⋆ is a minimum and, more specifically, vanishes at (θ 0 , σ) = (π/2, 0). Similarly, Y ⋆ ⊥ is a maximum at (θ 0 , σ) = (π/2, 0) and vanishes at its minimum (θ 0 , σ) = (0, 0). This analysis establishes upper bounds 9 for Y ⋆ and Y ⋆ ⊥ as 9N kT and 9/4N kT , respectively; and shows that theoretically, zero stiffness can be achieved. However, as previously mentioned, there are many applications, such as soft robotics, when it will be desirable to control the DE stiffness within some range while optimizing over other properties. In this case, (4.7) and (4.8) can be used as design tools in the limit when σ ≪ 1; that is, when directionality in the network is highly controlled.
The approximations in (4.9) share the same physical character as the σ ≪ 1 approximations in the sense that, when σ is large enough, Y ⋆ ,∞ is at its maximum (on the interval [0, π/2]) at θ 0 = 0 and minimum at θ 0 = π/2 (this can also be seen in FIG. 4.3). The situation is reversed for Y ⋆ ⊥,∞ ; that is, its maximum is at θ 0 = π/2 and its minimum is at θ 0 = 0. Again, this is due to the fact that there is a greater stiffness in the directions of higher chain densities. Lastly, note that the approximations in (4.9) predict a nonphysical singularity on the interval θ 0 ∈ [0, π/2] when σ is not large enough (σ π). This can result in, incorrectly, predicting negative and/or diverging Y or Y ⊥ . Thus, one should take care that σ is large enough when using (4.9) for design of transversely isotropic elastomers.  4.3. The dimensionless elastic moduli in the directions of the axis of symmetry, Y ⋆ , and orthogonal to the axis of symmetry, Y ⋆ ⊥ , as a function of the upper Gaussian center, θ 0 , for finite σ: 4π/3 and 3π/2. The dashed, black line represents the elastic modulus for an isotropic network, Y iso = 3G iso = 3N kT .

Susceptibility
Using the negligibly interacting assumption, we have that W * = N F r . However, (2.10) is challenging to evaluate in closed-form. Proceeding as in the case of the elastic modulus, we use a Taylor expansion of the inverse Langevin function about zero stretch to find γ/L −1 (γ) = 1 3 − 1 5 This provides an approximation for the free energy: The four terms on the right side of (5.1) can be understood as follows.
1. the first term is the closed-dielectric free energy density of a collection of (N n) monomers that are kinematically free; that is, monomers that are not constrained to satisfy a given end-to-end vector; 2. a correction to the first term that is related to the average magnitude of chain stretch; 3. a further correction that takes into account the average amount that the chain is stretched parallel to the direction of the electric field; 4. higher order terms related to the finite extensibility of the chain. Note that the first term is invariant when changing the cross-linking density at fixed mass density.
The second and third terms contribute to electromechanical coupling, while only the third term captures the effect of chain torque. Interestingly, if one neglects O γ 4 n terms in (5.1), then it is invariant under α. Therefore changing the fraction of loose end monomers can only have an effect on higher order terms, terms which are related to the finite extensibility of the chain.

5.A. Definition of the Susceptibility
In Section 3.F we defined the inverse susceptibility tensor: where the derivative is evaluated at the reference state. However, (5.1) gives us W * ; that is, the Legendre transform of W in the polarization slot. Further, although it is clear that (where by E P we mean the electric field as a function of polarization), it is not straightforward to write out W explicitly. This is because, while we have derived the polarization as a function of the electric field, we are unable to invert this function, for general F andP , to obtain E = E P . Therefore, we instead look for a correspondence between derivatives of W and of W * . The desired identity is as follows 10 : This is shown in Appendix B.

5.B. Tangent and Secant Susceptibility
We return our attention to the first term in (5.1)-the term that corresponds to a dielectric that consists of monomers that are unconstrained. We call the result of taking −∂/∂E of this term the free polarization.
The ratio of free chain polarization to its maximum theoretical value is shown in FIG. 5.1 as a function of the applied electric field. We distinguish between the secant susceptibility, X sec , and the tangent susceptibility X . The secant susceptibility is defined such that P = X sec E. Thus, fig. 5.1 shows a measure of the secant susceptibility of free monomers. Clearly X sec = X = − ∂ 2 W * ∂E∂E in general; in fact, X = X sec + ∂X sec ∂E E. They are equivalent when the dielectric is linear.
All of the dielectric nonlinearity in our DE materials is encoded in the function w * f , which is a per monomer contribution to the closed-dielectric free energy free energy for an unconstrained polymer chain (see (5.2)). In fact, the closed-dielectric free energy of an unconstrained polymer chain is nkT w * f − κ ⊥ . Taking The chains consist of TI (|∆χ| = χ ⊥ ) and uniaxial (|∆χ| = χ ) monomers. In the limit of small electric field or large temperature (i.e. |κ| → 0), the nondimensional (lab) susceptibility per monomer is 2/3 for TI monomers and 1/3 for uniaxial monomers. In the limit of large electric field or small temperature, the chain polarization approaches its theoretical maximum: nχ M E.
The first regime corresponds to |κ| → 0; then X sec /n → χ + 2χ ⊥ /3. In this limit (i.e. E → 0 and/or kT → ∞), the pdf of monomer orientations is uniformly distributed over the unit sphere because the electrostatic energy is vanishingly small as compared to the thermal energy. The factors of 1/3 for χ and 2/3 for χ ⊥ correspond to the dimensionality of each dipole susceptibility: the monomers are in three dimension space while χ is the dipole susceptibility along a line (i.e. spann) and χ ⊥ is the dipole susceptibility in the plane orthogonal ton. In addition, note that X sec does not vanish as kT → ∞. This is because χ µ is quadratic inn. Even thoughn is uniformly distributed, the polarization does not cancel.
There are two main takeaways that are relevant to this discussion: (1.) the dielectric nonlinearity of our DE materials arises because the average alignment of the monomers is determined by a balancing of electrostatic potential energy and the entropy, and because the dipole susceptibility of a monomer depends on its alignment with respect to the applied electric field; and (2.) despite this nonlinearity, we have X = X sec in the limits of |κ| → 0 and |κ| → ∞. Moving forward, we will consider both X and X sec ; as they will both prove useful. Specifically, X is particularly relevant to material stability while X sec is relevant to the capacitance-and hence, electromechanical coupling-for a given thin film DEA geometry. Moreover, for convenience, we make the definitions: We now consider what we have learned about free polarization in the context of the design of anisotropic dielectric elastomers. When used for a DEA, it is desirable that the secant susceptibility of the DE be as large is possible. This is because the susceptibility serves to increase the capacitance of the DEA. The increased capacitance means a greater accumulation of charge on the electrodes for a given voltage difference-which leads to a greater Coulomb attraction between the electrodes and hence a greater electromechanical coupling of the DEA. Now, for illustrative purposes, consider a single chain. If the chain is allowed to contract to zero stretch, then all the terms in (5.1) vanish except the first one-which is consistent with our previous discussion. In this case, the susceptibility is isotropic and is such that X sec = X sec I, X sec ∈ N n χ + 2χ ⊥ /3, max χ , χ ⊥ . Recall (FIG. 5.1) that the secant susceptibility is maximized at low temperature or large electric field (i.e. |κ| ≫ 1). This is an issue, from a practical standpoint, because we would prefer not to need to operate our DEA under these conditions. Both temperature control as well as large applied field are challenging and cumbersome to impose. Now consider this same chain, but in the reference, load-free state of an isotropic network such that γ r = 1. Then the second and third terms-the electromechanical terms-are O (1) while the first is of O (n). This means, for a dielectric elastomer that consists of "long chains" (i.e. n ≫ 1, n 1000), that the electromechanical terms are negligible when the elastomer is in the load-free state (and in many cases, they are even negligible at large macroscopic deformations). Thus, if we are to significantly improve on the small |κ| secant susceptibility in the load-free state, then one may first want to increase the density of cross-links of the elastomer (i.e. increase N 0 /n)-while keeping in mind that the stiffness also scales with the density of cross-links. Before moving on, it is also worth recalling the discussion on the "residual stresses" of our hypothetical anisotropic elastomers. Since the load-free state has some initial deformation to relieve the "residual stresses", γ r = 1 in general. At first glance, one may consider this as an opportunity to increase X h ⋆ . However, it is easy to show that γ 2 r ⋆ r ≤ 1 and that equality only holds for an isotropic network. The key is that, in the absence of electrical loads, the relaxed state should maximize the entropy of the elastomer and, by our approximation, the chain entropy as proportional to −γ 2 r ; thus, for entropy to not decrease with respect to the manufactured state (i.e. −γ 2 r r = 1), we require γ 2 r ⋆ r ≤ 1. Similarly, the lower bound on γ 2 r Ê ·r 2 ⋆ r is clearly zero since it is strictly nonnegative and it vanishes when the elastomer is manufactured such that all of the chains are orthogonal to the eventual direction of the applied electric field. The upper bound, however, is much less clear and warrants further investigation. For the purposes of this investigation, we split the susceptibility into two contributions: one associated with free monomers, X free , and a correction term due to the electromechanical coupling of the material,X ; that is: Taking derivatives of (5.1) with respect to E and taking the appropriate limits, we obtain:

5.C. Design for susceptibility
We again consider transversely isotropic materials and use (3.2) as our ansatz for P. This gives us: (5.5) In the limit of σ ≪ 1, we have the approximation: We could then derive, using (5.5) and (4.3), W * for transversely isotropic DEs. However, we turn our attention directly toX h . In the limit of σ ≪ 1, we have the approximation: where again, ⋆ denotes a quantity evaluated at λ = λ ⋆ . Interestingly, the electromechanical correction term vanishes at λ ⋆ for all θ 0 and σ when σ ≪ 1; and hence, at least in the limit of the load-free state, the DE effectively behaves as a collection of free monomers.
Similarly, using (C2) in (4.3) and (5.5), then subsequently (5.3) and (4.9), we arrive at an approximation of X h ⋆ in the limit of weak directional control (i.e. σ ≫ 1). This approximation is: where f was first defined in (4.9). The approximation is shown in FIG. 5.2. Notice that, while the electromechanical correction (at λ = λ ⋆ ) does not vanish, X h ,∞ ⋆ /N ∆χ is small compared to 1 in the limit of σ ≫ 1. Since the free monomer susceptibility is O (N n), this correction term is negligible; and, what is more, is that it vanishes as σ → ∞.
The electromechanical correction to the susceptibility in the load-free state, X h ⋆ , vanishes in the limit of σ ≪ 1 and is negligible in the limit of σ ≫ 1. It would seem then that there is little hope for increasing the initial susceptibility of our design DEs. However, it is clear from (5.3) thatX h is deformation dependent.
Thus, to maximize X h ,0 ⋆ , one must find a way to alter the polarization properties while simultaneously maintaining a (nearly) mechanically isotropic network.
Correction to the free monomer susceptibility in the load-free state in the limit of σ ≫ 1. The correction to the free monomer susceptibility is much less in this limit than the σ ≪ 1 limit, and it vanishes as σ → ∞.

5.D. Hybrid networks: maximizing the operating susceptibility
It was clear, particularly in FIG. 5.3, that X h ⋆ could be made greater if somehow the initial, load-free deformation, λ ⋆ , could be as close to unity as possible. For this reason, we propose the following (rough) manufacturing process and design algorithm for optimizing X h ⋆ . First, imagine that we have two types of polymer chains that are compatible with each other in the sense that they can be cross-linked together in a network. Further, imagine that the monomers for the one type of chain are such that ∆χ > 0 and the other are such that ∆χ < 0. Call these type A and type B, respectively. Now, if an electric field is applied in a constant direction just prior to and during cross-linking, then the density of A chains oriented orthogonal to the electric field direction will increase, as will the density of B chains oriented parallel to the electric field. Let σ A and σ B be the design standard deviations for chains of type A and type B, respectively. Then, clearly, σ A and σ B would depend on the magnitude of the applied electric field; and, σ A and σ B would likely not be able to be controlled independently of each other-or rather, the envisioned manufacturing process would need to be specialized further in order to control the two independently of each other. In this scenario, the monomer susceptibilities χ A ⊥ , χ A , χ B ⊥ , and χ B are given; or at best, selected from a catalog of possible monomer types. Then the design variables consist of 0 is the number of chains of type A per unit volume, n A is the number of monomers per chain in chains of type A, θ 0 is the center of the Gaussian in the upper half of the unit sphere for chains of type A, σ A is the standard deviations of the pdf for chains of type A, and the remaining quantities are the same but for chains of type B. The closed-dielectric free energy density for this hybrid DE is: where The design space for this problem is of a higher dimension than those that we have considered thus far. Therefore, we proceed by using some of the intuition that we have gained thus far. By FIG. 5.2 and FIG. 5.3, we reason that an optimal X h ⋆ will result from taking θ A 0 = π/2, θ B 0 = 0, and σ A = σ B = 0. Now let N m 0 := N A 0 n A + N B 0 n B be the number of monomers per unit volume in the reference configuration and N m := J −1 N m 0 be the number of monomers per unit volume in the current configuration. The free monomer susceptibility scales with the number of monomers per unit volume; thus, assume that N m 0 is already taken as large as possible or desired given some range of acceptable DE mass density. Similarly,X h scales with N 0 , but so does the stiffness. Thus, one should take N 0 as large as possible while still keeping the stiffness within some desired range. Let considering FIG. 5.3, we want to pick Ξ A ∈ [0, 1] such that λ ⋆ = 1. In other words, we require that Ξ A satisfies: which sets Ξ A = 2/3.
Next, we need to pick n A and n B . We define the ratios ξ A := N 0 n A /N m 0 and ξ B := N 0 n B /N m 0 . Since, by assumption, the total number of monomers is already given, by definition, we require: Consequently: ξ A ∈ 0, 1/Ξ A and ξ B = 1 − Ξ A ξ A / 1 − Ξ A . Because the only contribution to X h that depends on N m is the free monomer susceptibility, if χ A + 2χ A ⊥ > χ B + 2χ B ⊥ then we try to reach the limit ξ A → 1/Ξ A . Where as, if χ A + 2χ A ⊥ < χ B + 2χ B ⊥ then we try to reach the limit ξ A → 0. Lastly, if χ A + 2χ A ⊥ = χ B + 2χ B ⊥ then, in our truncated theory, X h free is invariant with respect to ξ A (in this case, it is likely best to let n A = n B = N m /N 0 ). Let M = arg max χ A + 2χ A ⊥ , χ B + 2χ B ⊥ . If the limiting process is successfully carried out such that X h free is maximized, then: (5.6) In summary, the proposed design process for maximizing X h is as follows: 1. Choose a preferred mass density, N m 0 . 2. Let θ A 0 = π/2 and θ B 0 = 0. Try to approach σ → 0. 3. Let Ξ A = 2/3 and, consequently, Ξ B = 1/3. 4. Maximize N 0 /N m 0 without exceeding the desired stiffness threshold(s). Note, when Ξ A = 2/3, Ξ B = 1/3: The target ξ A should be determined by: However, approaching either the upper limit of ξ A , which would result in chains of type B to be "short" (i.e. n B small) or the lower limit of ξ A , which would result in chains of type A to be "short", may affect the electromechanical response of the DE. This is because, as chains become shorter, higher order terms in (5.1) become more relevant at even moderate deformations. Specifically, these higher order terms would cause monomers to be constrained toward the direction of chain stretch more quickly as the DE deforms. This will lead to an increase in both strain hardening and electromechanical coupling. An additional consideration is the effect of chain length on the validity of the negligibly interacting assumption; it could be that, as chain lengths become shorter, interactions between chains become more important. Thus, either limit should be approached iteratively until it is determined how closely the limit can be approached without affecting the desired stiffness properties.
We have shown that X h free ⋆ can theoretically be improved upon by deliberately designing and manufacturing the network architecture of a dielectric elastomer. However, practically speaking, it can be seen from (5.6) that, unless the network has a high density of cross-links and therefore consists of "short chains" (i.e. N m 0 /N 0 ≤ 10) the increase is modest at best (e.g. 1% for N m 0 /N 0 100); and a high density of cross-links may not be desirable because the stiffness per mass density scales linearly with the density of cross-links. In addition, the theoretical improvement may be lost entirely once manufacturing error inevitably occurs.
While this analysis shows that the susceptibility can only be marginally improved at best, it does provide useful insights: (1.) it provides theoretical limitations of what can be achieved through the design of dielectric elastomer network architectures, and (2.) although the gains are modest for X h free ⋆ , we will see below that there can be significant increases in the electromechanical coupling. Indeed, recall that X h free is also a function of deformation; so that larger increases in susceptibility may be realized at deformations such that λ = λ ⋆ . This effect will be investigated in the next subsection.

5.E. Mechanically induced susceptibility
The susceptibility of dielectric elastomers is deformation-dependent, because deformation can cause chains in the network to rotate (thereby changing the average monomer orientation in each of the chains) and cause chains to stretch (thereby increasing the concentration of monomers oriented toward the direction of stretch for a given chain). Using the results developed thus far, we consider a few examples.
First, as a baseline, we consider an isotropic dielectric elastomer such that P = 1/4π. Using (4.3), (5.5), and (5. Notice that when ∆χ > 0,X h uni increases (relative to the reference configuration) when the DE is compressed in the direction of the axis of symmetry (λ < 1) and decreases when stretched; when ∆χ < 0, vice versa.
For the hybrid network, where θ A 0 = π/2, θ B 0 = 0, σ A = σ B = 0, Ξ A = 2/3, Ξ B = 1/3, we have: Interestingly, since both coefficients are strictly nonnegative and physically we require λ > 0, it is the case thatX h hybrid is semi-convex in λ for the domain of admissible λ. Further, it is convex when |∆χ A | = 0 and |∆χ B | = 0. This is significant because it means that when either compressing (i.e. λ < 1) or stretching (i.e. λ > 1) -even though initially there may be a drop inX h -eventuallyX h will begin increasing again and do so monotonically. The electromechanical increase of X h is bidirectional for the hybrid network. If |∆χ A | = 2|∆χ B |, thenX h has its minimum at λ = 1 so that any deformation increasesX h .
Also, importantly, (5.7) and (5.8) show that the mechanically induced susceptibility of the hybrid network is greater than that of the uniform network for all admissible deformations, λ.
Next, in contrast to either a uniformly distributed network or the hybrid network, we consider a network that consists of a single type of monomer and has been manufactured in the limit of high control, i.e. σ ≪ 1. In this case: using (C1), (4.3), (5.5), and (5.3): (∆χ) (1 + 2 cos (2θ 0 )) e −4σ 2 − 3 λ −1 + 2 (1 + 2 cos (2θ 0 )) e −4σ 2 + 2 λ 2 (5.9) When θ 0 = 0 the coefficient of λ −1 vanishes and when θ 0 = π/2 the coefficient of λ 2 vanishes. This is expected because the sign and coefficient of λ −1 and λ 2 determine the effect of compression and stretching, respectively, onX h . The above equation could be used as a design tool for anisotropic elastomers that consist of a single monomer type. For further physical insight, we let σ = 0. Then (5.9) simplifies to: Importantly, the theoretical factor for mechanically induced susceptibility can be larger in this case than it is for the uniform or hybrid networks. That is because, in this case, we can orient all of the chains in their preferred electromechanical susceptibility direction instead of a fraction. The maximum factors occur at θ 0 = 0 and θ 0 = π/2, as expected. At θ 0 = 0, the λ −1 factor (i.e. compression factor) vanishes and the λ 2 factor (i.e. expansion factor) is maximized-its value being 2N/5. And at θ 0 = π/2, the λ 2 factor vanishes and the λ −1 factor is maximized-its value being N/5. These factors are 3/2 and 2 times larger than the hybrid network factors, respectively. However, recall that Y ⊥ vanishes as θ 0 → 0 and Y vanishes as θ 0 → π/2. Thus, while a greater mechanically-induced susceptibility can be achieved (over the uniform and hybrid networks), there is a trade-off in terms of stiffness and mechanical stability. The implications of this in terms of the electromechanical coupling, operation, and failure of DEAs will be explored in the next section.

Cross modulus, Electrostriction, and Material Stability
We now turn our attention to the electromechanical or cross modulus, ∂ 2 W ∂P ∂F , which provides a measure of the electromechanical coupling of the material. Among other things, we will show here that it is related to both the stress that develops in the DE when loaded by a constant electric field and as well as the stability of the material.
To derive the cross modulus of our DEs, we begin with the relation: Let K denote the cross modulus tensor. Then, taking derivatives of both sides of the above equation with respect to the deformation gradient, we obtain: which can be directly related to derivatives of W * , as outlined in Section 5. In order to highlight some important properties of the cross modulus, we will again consider the setting of a thin film as described in Section 2.E, but now using transversely isotropic materials. By symmetry, we expect F = diag 1/ √ λ, 1/ √ λ, λ in a Cartesian basis with the 3 direction aligned along the thickness direction. Further, the symmetries of the material, deformation, and applied electric field provide a simplified form of (6.1) along the axis of symmetry Note that since K ∝ E, we do not see a frozen in polarization (i.e. the material does not spontaneously become an electret) for any relaxation stretch, λ ⋆ , as expected.
For brevity, we restrict our attention to the limit of |σ| ≪ 1. The cross modulus as a function of θ 0 is shown in FIG. 6.1 and FIG. 6.2 κ = 1 and κ = −1, respectively. For both figures, n = 100. Interestingly, the absolute value of the cross modulus does not have its maximum, in either case, at θ 0 = 0, π/2, or θ iso ; and it does not decay as the distribution spreads out (σ > 0). Instead, the cross modulus vanishes for σ = 0 at θ 0 = 0 and π/2.
The significance of the cross modulus is subtle and is perhaps best understood through the lens of the experiment shown in FIG. 2.1: a constant electric field, E 0 , is applied to a dielectric elastomer with (effectively) traction free boundary conditions. Since the deformation gradient and electric field are both homogeneous, the Gibbs free energy minimization ((2.13)) can be carried out pointwise. Furthermore, the polarization is work-conjugate to the applied electric field; thus, the Gibbs free energy density is given by: W λ,P −P E 0 or, equivalently, W * . Now consider an anisotropic DE that has been manufactured and allowed to relax to λ = λ ⋆ . Then, after relaxation, the electric field is applied. The stress that develops in the DE as a result of the electric field (denote by Σ) is given by: There is a clear connection between the cross modulus and the stress as they are both proportional to change in the secant susceptibility with respect to deformation. After our discussion in Section 5.E, there should be no surprise that this contribution to stress exists. The existence of a mechanically induced susceptibility implies a stress induced by changing susceptibility, and vice versa.

2/3
Notice that the stress is also proportional to κ and its sign depends on whether the chains consist of TI monomers or uniaxial monomers. The stress attains its maximum absolute value at θ 0 = 1 2 arccos 1 3 which, to a good approximation, is also where K ⋆ attains its maximum absolute value.
There are a few important takeaways here that will be relevant to the design of DEs for dielectric elastomer actuators in Section 7. First, there are two apparent contributions to the electromechanical coupling in DEAs: (1.) the Coulomb attraction between the top and bottom electrodes of the DEA and (2.) the contribution to the stress in the DE due to the change in the secant susceptibility with respect to deformation. As discussed previously, the former can be increased by increasing the susceptibility of the DE. The latter, as we have shown in this section, can be increased by increasing the cross modulus of the material. Since the susceptibility and cross modulus are optimal for different network architectures, there will be a competition between the two when optimizing the network architecture for a DEA.
The second takeaway is this: the amount of electrostriction of the anisotropic DEs under consideration, for a fixed E 0 , is independent of the chosen network architecture; that is, if we make the decomposition λ = λ ′ λ ⋆ where λ ⋆ is the initial relaxation of the DE and λ ′ is the deformation that occurs as a result of the electrically induced stress, then λ ′ is invariant with respect to P, N , n, and α. This is indeed surprising, but we can begin to explain it as follows: 1. There are two consequences of the linearized theory that result in invariance with respect to α and n: (a) As before, α is not accounted for in the linearized free energy density because it is has a higher order effect related to the finite extensibility of the chains in the network. (b) Although the change in susceptibility with respect to deformation does depend on n (since it depends on the absolute stretch of the chain), this is also a higher order effect related to finite extensibility. 2. The invariance with respect to n occurs because both Y and Σ are linearly proportional to N kT . 3. The invariance with respect to P is surprising since there is clearly an influence of the parameters (θ 0 , σ) on the cross modulus and the stress. However, note that Y also depends on (θ 0 , σ); and, importantly, K and Y both depend on the deformation (see (6.1) and (4.5), respectively). This interplay between θ 0 , σ, λ, Y and K results in λ ′ being invariant with respect to θ 0 and σ.
For completeness, we also provide the formula for λ ′ : where w * f was defined earlier in (5.2).

6.A. Material stability
In addition to the electromechanical coupling of our anisotropic, design DEs, we will also be interested in designing against failure; here, we define failure as loss of stability of the free energy, not considering fracture, breakdown &c. For this purpose, consider again the Taylor expansion of the Helmholtz free energy density given in (3.3). However, here we take the expansion about a general equilibrium configuration, λ eq ,P eq . Then, for perturbations about the equilibrium configuration, ∆λ, ∆P , the change in the free energy density is: where |ǫ| = max |∆λ|, |∆P | ≪ 1 is a small parameter and the superscript, eq , on the material properties are there to emphasize that each of the second partial derivatives are evaluated at ∆λ, ∆P . Moving forward, we drop the superscript as it will be clear what is meant from the context. Next, define the Hessian of the free energy density: For a material to be stable (at a given equilibrium configuration), we require that the energy increase for all allowable perturbations 11 . Here, we neglect the role of constraints, thus, a sufficient condition for stability is that H be positive definite.
As a first step toward understanding the properties of H, let us consider its diagonal: Y and X −1 . Directly, 11 As pointed out in [Eri98], although we are actually interested in changes in the free energy (as opposed to the free energy density), it often suffices to consider stability in the context of the free energy density. Here we do not consider global instabilities, such as wrinkling and buckling, which are often geometry and/or boundary condition dependent, and instead focus on local failure modes. For dielectric elastomers, the principal local failure mode is characterized by a dramatic thinning [ZS07,ZDDP17]. Typical global failure modes include wrinkling, buckling, and dielectric breakdown. See [ZDDP17,YZS17] for a more comprehensive analysis and discussion. Note that, when failure is studied and observed in dielectric elastomers, the material is generally subjected to a Maxwell stress. Importantly, we show later in this section that, for both standard isotropic and the proposed anisotropic materials alike, a local failure mode will not occur in the absence of a Maxwell stress.
It is not difficult to see that the quantity enclosed in the braces is nonnegative and 45 + 22κ − 30w * f is strictly positive. Thus, the stiffness is never negative in this case. Similarly, although the formula for X is tediously long, we can easily reason that it too is nonnegative. The quantity X is the result of averaging across scales, starting at the monomer dipole susceptibility. Since the monomer dipole susceptibility is nonnegative, the subsequent averages that we obtain must also be nonnegative. Further, X −1 > 0, since by similar reasoning there must be an upper bound on X (this upper bound is, of course, N nχ M where, as before, χ M = max χ , χ ⊥ ). Now, ordinarily H ≥ 0 would only be a necessary, and not a sufficient, condition for the positive semi-definiteness (positive definiteness when H > 0) of H. However, since the trace is also nonnegative, the above condition is also sufficient. Thus, if H > 0 then stability is guaranteed; otherwise, the material may be unstable. We forego the nuances of exactly when the design DE is unstable, for simplicity, and focus on how H depends on different network architectures and loading conditions.
It is reasonable to then ask: does the stability picture change as |κ| increases? FIG. 6.4 and FIG. 6.5 show H as a function of |κ| for architectures consisting of TI and uniaxial chains, respectively. For each figure, n = 100, σ = 0, and θ 0 is varied by series: π/12, π/6, π/4, π/3, and 5π/12. FIG. 6.4 shows that, for architectures consisting of TI chains, H decreases slightly before increasing with |κ|. Importantly, it does not approach zero in any case. For architectures consisting of uniaxial chains (FIG. 6.5), H decreases with |κ| before leveling off asymptotically to some number greater than zero. Therefore, we see, these anisotropic DEs do not become unstable when loaded by a constant electric field.
We conclude this section by explaining why we have not seen instability in these materials. The reason becomes apparent when we consider how Y , K and X scale. Specifically: And, consequently,  The θ 0 considered are π/12, π/6, π/4, π/3, and 5π/12. H decreases slightly before increasing with |κ|; and importantly, it does not approach zero. Therefore, material instability does not occur.
Since in many cases n ∼ 100-10000, the strictly nonpositive contributions to H are negligible compared to the strictly nonnegative. The main reason why this scaling occurs is because X scales with N n but The θ 0 considered are π/12, π/6, π/4, π/3, and 5π/12. H decreases with |κ| before leveling off asymptotically to some positive value; and importantly, it does not approach zero. Therefore, material instability does not occur.
In the next section, we explore how the performance of a dielectric elastomer actuator changes with the network architecture. The significance of the results developed in this section will be two-fold. First, as mentioned in Section 6, the electrostrictive stress that develops in the DE film due to its polarization response will effect the performance of the DE as an actuator; and, importantly, the electrostrictive stress is not optimal for the same architectures in which the polarization susceptibility is optimal. Secondly, we will also consider the failure of DEAs due to instability. When we do, we will be justified in dropping the Hessian-based approach and instead only consider the sign of ∂ 2 G/∂λ 2 -as the sign of ∂ 2 G/∂λ 2 captures the instability associated with the Coulomb attraction overcoming the stiffness of the film.

Application to Dielectric Elastomer Actuators
In this section, we explore the effect of our design parameters on the deformation and usable work obtained from a dielectric elastomer actuator (DEA) as a function of its electrical input. There are two main goals associated with this design: (1.) we would like to maximize the deformation and/or usable work that results from a fixed electrical input; and (2.) we would typically like to maximize the deformation and/or usable work that the DEA can produce before its failure.
Let the dimensions of the DEA in the reference configuration be l × l × l 3 , where l 3 ≪ l is the thickness of the DE film. Following [YZS17] and [ZS07], the Gibbs free energy of a DEA undergoing homogeneous deformation, F = diag λ −1/2 , λ −1/2 , λ , is: where the first term in the parentheses is the free energy density of the DE film and the second term in the parentheses is the energy density of the electric field; and the last term in (7.1) is the work of the battery. However, (7.1) is in terms of W instead of W * . Recalling that W = W * + P · E and using Gauss's law: In deriving (7.2), we have assumed that E = − (∆ϕ) /λl 3 inside the DE film, which is a good approximation when l 3 ≪ l and the DE film has homogeneous material properties.
As mentioned previously, there is an initial relaxation, in general, after manufacturing an anisotropic dielectric elastomer. To account for this, we propose the following process: 1. The DE film is manufactured; it relaxes to some λ = λ ⋆ . In anticipation of this relaxation, its initial, pre-relaxation thickness is l 3 /λ ⋆ . 2. A voltage difference, ∆ϕ, is applied across the electrodes. 3. The DE film deforms further by some amount λ ′ such that λ = λ ′ λ ⋆ .

7.A. Deformation
Let E :=Ẽ nX h free /kT . This is a nondimensional measure of the electric field, similar to the dimensionless electric field used in [ZHS07,ZDDP17]. Because we would like to emphasize optimization here, we will focus on three types of networks: (1) a uniform (isotropic) network as a baseline, (2) network designs in the limit of σ = 0, and (3) the aforeproposed hybrid network. Similarly, because (for unitype networks) elastomers that consist of TI monomers will have an greater electromechanically induced susceptibility than elastomers consisting of uniaxial monomers, we will consider elastomers such that ∆χ = χ ⊥ for the uniform and σ = 0 networks. First, FIG. 7.2 shows λ ′ vs θ 0 contours of constant E for networks of type (1) and (2). Note that the end of each data series (vertical line to the x-axis) represents failure of the DEA for the corresponding network architecture and applied electric field. It can be seen in FIG. 7.2 that, for each of the contours, θ 0 = θ iso has the least amount of deformation. In fact, the contours are symmetric about θ iso . However, while the deformation increases as θ 0 has a larger deviation from θ iso , the DEA also fails at smaller E. Thus, there is a trade-off between the maximizing the deformation for a fixed E and maximizing the operating E before failure-and, consequently, the maximum deformation possible. Put differently, there is a trade-off between electromechanical efficiency and stability. Note that it can be seen from numerical examples that θ 0 = θ iso has an equivalent electromechanical response to the uniform network. Thus, θ 0 = θ iso represents our isotropic baseline.
Next, we consider the hybrid network within this context. FIG. 7.3 compares the DEA deformation of the hybrid network to the unitype networks with σ = 0. For a given E, the hybrid network deforms more than the unitype isotropic network (i.e. θ 0 = θ iso ). However, for E low enough, there is always some θ 0 , σ = 0 for which the unitype network has a larger induced deformation than the hybrid network.

7.B. Usable Work
Another performance metric for our DEAs is the amount of usable work that can be derived from the system. We will use −∆G as a measure of the usable work, where FIG. 7.4 shows −∆G/N 0 kT vs θ 0 contours of constant E for networks of type (1) and (2). Similar to the case of deformation, it can be seen that, for a given E, the usable work can be maximized by picking the θ 0 such that |θ iso − θ 0 | is maximized without an instability occurring.
Next, consider the usable work of the hybrid network compared to the unitype σ = 0 networks. FIG. 7.5 shows −∆G/N 0 kT as a function of E for σ = 0; θ 0 = π/6, π/4, θ iso , π/3 and the hybrid network. Again, it can be seen that for a fixed E, a unitype network with properly chosen θ 0 can outperform the hybrid network. However, the hybrid network combines both a higher electromechanical coupling (than the isotropic network) and maintains its stability at larger E. In particular, the hybrid network shows ≈ 75% increase in usable work over the isotropic network for general E. While σ = 0, θ 0 = θ iso can endure a larger E before failure, the maximum usable work before failure of the hybrid network is still greater than that of the isotropic network. In summary the unitype, σ = 0 networks can be optimized for a specific electrical load more so than the hybrid network. However, the hybrid network is more stable and would be much more  (1) and (2). The end of each data series (vertical line to the x-axis) represents failure of the DEA for the corresponding network architecture and applied electric field. For each of the contours, θ 0 = θ iso has the least amount of deformation and the contours are symmetric about θ iso . However, while the deformation increases as θ 0 has a larger deviation from θ iso , the DEA also fails at smaller E. Note that θ 0 = θ iso has an equivalent electromechanical response to the uniform network-so θ 0 = θ iso represents our isotropic baseline. Notice: for a given E, the hybrid network deforms more than the unitype isotropic network (i.e. θ 0 = θ iso ). However, for E low enough, there is always some θ 0 , σ = 0 for which the unitype network has a large induced deformation than the hybrid network.  (1) and (2). For each of the contours, θ 0 = θ iso has the least amount of usable work and the contours are symmetric about θ iso . However, while the deformation increases as θ 0 has a larger deviation from θ iso , the DEA also fails at smaller E. Note that θ 0 = θ iso has an equivalent electromechanical response to the uniform network-so θ 0 = θ iso represents our isotropic baseline.
preferable than the unitype network when used in an application with a wider range of operating E.

Summary of Findings
In this work, we utilized a multiscale modelling approach to investigate how the microstructure of a dielectric elastomer affects its material properties and performance as an actuator. The multiscale model lead to important insights such as the phenomenon of electrostatic chain torque, which not only contributes to the electromechanical coupling of the material but could be used as a possible mechanism for aligning chains during the DE manufacturing process. Towards material design, we proposed the following design parameters: the distribution of chain end-to-end vectors in the network, P; the mass density (i.e. N 0 n); the density of cross-links; and the fraction of loose-end monomers, α. Further, we restricted our attention to forms of P that would result in a transversely isotropic material. It was found that substantial gains in the deformation and usable work, for fixed electrical input, can be obtained by designing and manufacturing an anisotropic DE material.
Our key findings include the following: 1. To a good approximation, the electromechanical coupling of the network was invariant with respect to α. 2. The tangent stiffness modulus in the direction of the axis of symmetry, Y ⋆ , is a maximum for (θ 0 , σ) = (0, 0) (Y ⋆ = 9N kT and vanishes as θ 0 → π/2, σ → 0. Similarly, Y ⋆ ⊥ , is a maximum for (θ 0 , σ) = (π/2, 0) (Y ⋆ ⊥ = 9/4N kT and vanishes as θ 0 → 0, σ → 0. These results are intuitive as they suggest that the elastomer is more stiff in directions in which more chains are aligned. In general, the stiffness, at fixed density, is proportional to the density of cross-links in the network. 3. The total susceptibility scales as O (nN ) while (to a good approximation) the influence of the network architecture on the susceptibility only scales as O (N ). This means that unless the cross-linking While the most usable work can be derived from a unitype, σ = 0 network with properly chosen θ 0 for a given E, the hybrid network combines both stability and an enhanced usable work. In fact, the hybrid network performs ≈ 75% better than the isotropic network.
density is high-which would reduce compliance-the potential increase in susceptibility from tuning the network architecture would likely only be modest. 4. Surprisingly, the true electrostriction of a transversely isotropic DE is invariant with respect to its network architecture. 5. Despite the gains in susceptibility that result from adjusting the network architecture being somewhat modest, its influence (as well as the influence of some of the other aforementioned material properties) on the deformation and usable work of a dielectric elastomer actuator is nonlinear. In fact, the performance of a DEA can be significantly increased by properly designing the elastomer network. 6. For single chain type networks, there is a trade-off between the maximum amount of deformation and usable work that may be extracted from a DEA before failure and the efficiency of the DEA (i.e. mechanical output per electrical input). 7. A so-called "hybrid network" was proposed which consisted of uniaxial chains oriented withÊ and −Ê and TI chains oriented orthogonal toÊ. The hybrid network preserved the isotropic stiffness and electromechanical stability of the DE, while increasing its usable work output by ≈ 75%.
Outlook. Here we have chosen to focus on optimizing the microstructure of dielectric elastomers for a general type of actuator; for simplicity, we chose an ansatz for P with only two tunable parameters and assumed the material was homogeneous, However, one could imagine relaxing this and using homogenization, topology optimization, machine learning, and other tools to design a more complex, possibility heterogeneous, microstructure to achieve novel behavior such as bending actuation and shape morphing, or for optimizing the microstructure for specific applications in wearable electronics, energy harvesting, soft robotics, prosthetics, etc. In general, going beyond isotropic materials and exploiting emerging manufacturing techniques to introduce ordering into the material can lead to significant performance increases, following broadly along the lines of work in nonlinear homogenization [CG11, GC13, SC14].

A. Polarization as Energy Conjugate to Electric Field
We show that a consequence of assuming (2.8) as the form of the monomer density function and enforcing the constraints given in (2.9), we arrive at the relationship: p = − ∂F ∂E . This was first shown in [GD20] and is provided here again for completeness.
Taking derivatives of both sides of (2.9) with respect to E, we obtain: We are able to interchange the operations of derivation and integration because of the smoothness of the integrands; and in the last equalities we use the fact that neither the number of the monomers in the chain nor the end-to-end vector constraint depend on E.
Now, we obtain the desired result by taking derivatives of both sides of (2.10): It is obvious then that this is not true of all chains with polarizable monomers such that u includes a dipoleelectric field interaction term (i.e. −µ · E). Instead, this follows from the fact that the dipole susceptibility tensor, χ µ , is symmetric since it derives from a potential; and also from the fact that the energy to separate bound charges in the monomer is harmonic such that u bond = 1 2 µ · χ −1 µ µ.
As a corollary, we show that This follows as a consequence of − ∂F ∂E = p and the assumption that chains in the network are in negligible interaction (see Section 3.B). By the negligibly interacting assumption, W * = N F r , where by r we mean the average over the chain pdf. Hence, taking derivatives of both sides: (Recall: incompressibility implies that J = 1.) This result, along with the statistical mechanical derivation of W * , establishes W * (F , E) as the Legendre transform of W F ,P , the Helmholtz free energy density.

B. Derivatives of W and W *
Since we consider the physical regime in which monomer-monomer interactions are negligible (and thereby the nonlocal E-P relationship simplifies into a local one), and since we assume the material is incompressible, i.e. J = 1, the following hold: 1. W F ,P = W * F , E P +P · E P , 2. ∂W ∂P = E,
Assume E P invertible such thatP =P (E), ∂P ∂E is invertible, and W, W * , E P andP (E) are all smooth functions.
By assumption, ∂W ∂P = E P . Taking derivatives of both sides: Similarly, taking derivatives of both sides of ∂W * ∂E = −P : From the two above equations we can conclude that The physical implication of the above result is that: − ∂ 2 W * ∂E∂E = X .

D. List of Symbols
r Chain end-to-end vector (current configuration) r Chain end-to-end vector (reference configuration) n Number of monomers per chain b Length of a single monomer γ Chain stretch, |r|/nb N Density of chains per unit volume (current configuration) N 0 Density of chains per unit volume (reference configuration), N 0 = JN n b Number of backbone monomers n l Number of loose-end monomers α Fraction of loose-end monomers, n l /n P Probability density function which describes the fraction of chains with reference end-to-end vectorr θ 0 , σ Mean of pdf