# Elastic Effects in Disordered Nematic Networks

## Abstract

Elastic effects in a model of disordered nematic elastomers are numerically investigated in two dimensions. Networks crosslinked in the isotropic phase exhibit unusual soft mechanical response against stretching. It arises from gradual alignment of orientationally correlated regions that are elongated along the director. A sharp crossover to a macroscopically aligned state is obtained on further stretching. The effect of random internal stress is also discussed.

Nematic elastomers and gels exhibit rich mechanical effects due to elasticity-orientation coupling . While a considerable number of theoretical studies has been directed to homogeneous systems, nematic elastomers are often in a highly non-uniform polydomain state, in which the correlation length for the director orientation is typically of micron scales. Polydomain networks show unusual non-linear elastic response against stretching, often with an extremely low stress over a large interval of strain . As the strain is increased, the directors gradually rotate toward the direction of stretching until a macroscopically aligned state is attained. This structural change is called the polydomain-monodomain (P-M) transition. Attempting to describe the presumably equilibrium polydomain textures, Terentjev and coworkers  proposed a random-field model analogous to those for random anisotropy magnets. They argued that crosslinkers of anisotropic shapes act as sources of quenched disorder. On the other hand, the mechanical response is not yet well understood. It is known that elasticity-mediated long-range interactions among spatial inhomogeneities are crucial in systems such as metallic alloys  and gels . For polydomain networks, the role of elastic interactions among orientationally correlated regions (“domains”) is yet to be clarified. In this Rapid Communication, we numerically investigate the mechanical response and the domain structure of model nematic networks incorporating both rubber elasticity and quenched random anisotropy. Unusual soft response is obtained and is explained in terms of the elastic interaction. We briefly discuss the effect of random internal stress as another kind of quenched disorder which can destroy long-range orientational order .

The total free energy of our model system is of the form $`F=F_{el}+F_R+F_F`$, where $`F_{el}`$, $`F_R`$, and $`F_F`$ are, respectively, the rubber-elastic, random disorder, and Frank contributions. We assume networks brought deep into the nematic phase after crosslinking in the isotropic phase, and apply the affine-deformation theory of nematic rubber elasticity due to Warner et.al.. Then $`F_{el}`$ is written in terms of the symmetric deformation tensor $`W_{ij}=(r_i/r_k^0)(r_j/r_k^0)`$, where $`r_i^0`$ and $`r_i`$ are the Cartesian coordinates of the material point at the moment of crosslinking and after deformation, respectively. Summation over repeated indices is implied throughout the paper unless otherwise stated. It is convenient to rewrite the original form  of $`F_{el}`$ using the tensor $`Q_{ij}=n_in_j\delta _{ij}/d`$, where $`𝒏`$ is the director, to obtain 

$`F_{el}={\displaystyle \frac{\mu }{2}}{\displaystyle 𝑑𝒓(W_{ii}\alpha Q_{ij}W_{ij})}.`$ (1)   
The dimensionless coupling constant $`\alpha (>0)`$ is determined by chain anisotropy and does not exceed $`d/(d1)`$. The modulus $`\mu `$ is given by $`k_BT`$ multiplied by the crosslink number density and a numerical prefactor ($`1`$) which is weakly dependent on the temperature. We consider the incompressible limit and impose the constraint $`det𝖶=1`$. The disorder free energy is assumed in the form given in , and is rewritten as

$`F_R={\displaystyle 𝑑𝒓P_{ij}Q_{ij}},`$ (2)   
where $`P_{ij}`$ is a symmetric, traceless, Gaussian random tensor with vanishing quenched average ($`P_{ij}(𝒓)=0`$) and with variance

$`P_{ij}(𝒒)P_{kl}(𝒒)=U\left(\delta _{ik}\delta _{jl}+\delta _{il}\delta _{jk}{\displaystyle \frac{2}{d}}\delta _{ij}\delta _{kl}\right).`$ (3)   
For the Frank free energy we assume the form

$`F_F={\displaystyle \frac{K}{2}}{\displaystyle 𝑑𝒓(𝒏)^2}.`$ (4)   
Here we treat the two-dimensional case for numerical and analytical advantages. Then, in the absence of elasticity, our model reduces to the random-anisotropy 2D XY model by regarding the unit vector $`𝒎=(2Q_{xx},2Q_{xy})=(\mathrm{cos}2\theta ,\mathrm{sin}2\theta )`$ as the spin variable, where $`\theta `$ is the director orientation defined by $`𝒏=(\mathrm{cos}\theta ,\mathrm{sin}\theta )`$. We consider deformations of the form $`r_i=\lambda _ir_i^0+u_i`$ (no summation) where $`\lambda _x=\lambda `$ and $`\lambda _y=1/\lambda `$ express the average deformation, and $`𝒖=𝒖(𝒓)`$ denotes the internal displacement. Cooling into the nematic phase tends to induce spontaneous elongation along the director . If the directors are uniformly aligned along the $`x`$-axis, the elastic free energy (1) is minimized at $`\lambda =\lambda _m`$ and $`𝒖=\text{0}`$ with

$`\lambda _m=\left({\displaystyle \frac{1+\alpha /2}{1\alpha /2}}\right)^{1/4}.`$ (5)   
However, if the random field is sufficiently strong, the system under no external stress energetically favors a macroscopically isotropic state with $`\lambda =1`$. This can be identified with the polydomain state. Our questions concern how domains spontaneously deform and to what degree the elastic free energy is reduced in this highly non-uniform state.

The mechanical response was numerically simulated by varying the macroscopic strain $`\lambda `$ and minimizing the free energy with respect to $`𝒏`$ and $`𝒖`$ for each value of $`\lambda `$. We solved the Langevin equation for the director,

$`{\displaystyle \frac{𝒏}{t}}=(𝑰𝒏𝒏)\left(L{\displaystyle \frac{\delta F}{\delta 𝒏}}+𝜼\right),`$ (6)   
where $`𝜼`$ is an uncorrelated Gaussian thermal noise introduced to facilitate structural evolution. Without the noise the minimization process would stop at one of the local minima close to the initial configuration. After approaching the global minimum we turned off the noise as explained below. The displacement $`𝒖`$ was determined by solving the non-linear equation $`\delta (F_{el}+F_v)/\delta 𝒖=0`$ with a relaxation method, where $`F_v`$ is an artificial free energy functional of $`𝒖`$, which penalizes volume change. In this way the local volume was kept constant with errors below $`1\%`$ throughout the runs. Periodic boundary conditions were imposed on $`𝒏`$ and $`𝒖`$. The simulation was performed on a $`128^2`$ square lattice with the grid size $`\mathrm{\Delta }x=1`$. We set $`K=4`$ and $`U=1`$ for all the runs, whereas $`\mu `$ and $`\alpha `$ were varied for different runs. In each run the external strain $`\lambda `$ was slowly increased after an initial equilibration stage at $`\lambda =1`$. Occasionally, we stopped the increase of $`\lambda `$ and turned off the thermal noise for an interval of time. Thus a single run consisted of alternating periods of annealing (with increasing strain) and quenching. In each quench period we computed the spatially-averaged free energy density $`f=f_{el}+f_R+f_F`$ and the orientation $`S=2Q_{xx}=\mathrm{cos}2\theta `$. This procedure enabled us to approximately minimize the free energy at numerous values of $`\lambda `$ in reasonable computational time. For a further check, we then decreased $`\lambda `$ back from a large value in a similar manner. A small hysteresis was obtained but it does not affect the description below.

In Fig.1 we show the strain-stress and strain-orientation relations for several values of $`\alpha `$ with $`\mu \alpha ^2=4`$ fixed. In both plots we can see a sharp crossover around $`\lambda =\lambda _m(\alpha )`$. Below $`\lambda _m`$ the total stress $`f/\lambda `$ is vanishingly small and slightly positive. The average orientation increases almost linearly in the same region. The free energy densities are plotted in Fig.2. The elastic free energy has a slightly negative slope in the region $`\lambda <\lambda _m`$, while the disorder free energy has a positive slope and makes the total stress slightly positive. We chose the parameters so that the Frank contribution is much smaller than $`\mu \alpha ^2`$, which is considered to be the case in typical experiments. We also studied a few cases with stronger or weaker elastic effects. For larger values of $`\mu \alpha ^2`$ the shapes of the strain-elastic stress and strain-orientation curves were almost unchanged. For cases with $`\mu \alpha ^2\stackrel{<}{}\mathrm{\hspace{0.33em}0.2}`$ these two curves exhibited less sharp crossovers.

In order to understand the origin of the soft response it is useful to examine the structure of the polydomain state at $`\lambda =1`$, for which an analytical treatment is possible in the weak coupling case $`\alpha 1`$. We expand $`\mathrm{\Delta }F_{el}=F_{el}[𝒖]F_{el}[0]`$ with respect to $`𝒖`$ to obtain

$`\mathrm{\Delta }F_{el}=\mu {\displaystyle 𝑑𝒓\left[\frac{1}{4}\left(\frac{u_i}{r_j}+\frac{u_j}{r_i}\right)^2\alpha Q_{ij}\frac{u_i}{r_j}\right]}.`$ (7)   
Eliminating the elastic field using the conditions of mechanical equilibrium $`\delta \mathrm{\Delta }F_{el}/\delta 𝒖=0`$ and incompressibility $`𝒖=0`$, we have a non-local elastic interaction among orientational inhomogeneities. We define new variables $`Q_1(𝒓)`$ and $`Q_2(𝒓)`$ through their Fourier transforms,

$`Q_1(𝒒)`$ $`=`$ $`\mathrm{sin}(2\phi )Q_{xx}(𝒒)\mathrm{cos}(2\phi )Q_{xy}(𝒒),`$ (8)   
$`Q_2(𝒒)`$ $`=`$ $`\mathrm{cos}(2\phi )Q_{xx}(𝒒)+\mathrm{sin}(2\phi )Q_{xy}(𝒒),`$ (9)   
where $`\phi `$ is the azimuthal angle of the wave-vector $`𝒒=q(\mathrm{cos}\phi ,\mathrm{sin}\phi )`$. Then the average free energy density reads 

$`f_{el}|_{\lambda =1}=\mu \left(1{\displaystyle \frac{\alpha ^2}{2}}Q_1^2\right)`$ (10)   
to order $`\alpha ^2`$. Note that $`Q_1`$ and $`Q_2`$ satisfy $`Q_1^2+Q_2^2=Q_{xx}^2+Q_{xy}^2=1/4`$. We have $`Q_1^2=Q_2^2=1/8`$ in the absence of the elastic coupling. In its presence, asymmetry $`Q_1^2>Q_2^2`$ arises to reduce the elastic free energy (10). In the elasticity-dominated limit where $`\mu \alpha ^2`$ is much larger than the disorder and the Frank free energy densities, we expect $`Q_1^21/4`$, $`Q_2^20`$, and $`f_{el}|_{\lambda =1}\mu (1\alpha ^2/8)`$. Indeed these are numerically confirmed as shown in the next paragraph. On the other hand, the elastic free energy density under the uniform deformation with $`\lambda =\lambda _m`$ is also given by $`\mu (1\alpha ^2/8)`$ to order $`\alpha ^2`$. Thus, in the above limit, the P-M transition accompanies only a small change of order $`\alpha ^3`$ in the elastic free energy. To see how each domain is deformed at $`\lambda =1`$, we consider the local elastic stress which is given as $`\sigma _{ij}=\mu (_iu_j+_ju_i\alpha Q_{ij})`$ from (7). After some calculation, its variance in the mechanical equilibrium is obtained as

$`\sigma _{ij}^{}{}_{}{}^{2}=2\mu ^2\alpha ^2Q_2^2.`$ (11)   
In the elasticity-dominated limit the variance of the quantity $`\mu ^1\sigma _{ij}=_iu_j+_ju_i\alpha Q_{ij}`$ vanishes due to the factor $`Q_2^2`$ in (11), which means that each part of the system is elongated by $`1+\alpha /4(\lambda _m)`$ times along the local director. This, together with the numerical result on the mechanical response, supports the following simple picture : In the polydomain state each domain is uniaxially elongated by $`\lambda _m`$ times along the local director, and thus the elastic free energy is equal to that for the monodomain state at $`\lambda =\lambda _m`$ (Fig.3). The P-M transition in the region $`1<\lambda <\lambda _m`$ proceeds via rotation of domains and does not change the elastic free energy.

Next we present numerical results on the polydomain structure at $`\lambda =1`$, which was studied through the correlation function $`G(r)=2Q_{ij}(𝒓)Q_{ij}(0)`$ and the degree of structural asymmetry $`A=Q_1^2Q_2^2`$. To accelerate computation of the elastic field we assumed a weak coupling $`\alpha =0.1`$ and solved $`\delta \mathrm{\Delta }F_{el}/\delta 𝒖=0`$ under the constraint $`𝒖=0`$ using FFT, instead of the relaxation method above. The amplitude of the thermal noise was set constant in an initial stage and then gradually reduced to zero at a constant rate. The correlation function is computed for the final state and averaged over $`20`$ independent runs for each set of parameters. Runs were sufficiently long to insure that the initial configurations with uniform and random orientations give indistinguishable results for $`G(r)`$. Shown in Figs.4 and 5 are the correlation function and the correlation length $`R`$ defined by $`G(R)/G(0)=1/2`$. The elastic coupling increases the correlation length without qualitatively affecting the form of the correlation function. We could not deduce a quantitative decay law for $`G(r)`$ from the relatively small number of samples, but the decay was slightly faster than exponential near the origin. For the non-elastic case the same feature was obtained in the Monte Carlo simulation by Gingras and Huse  in the presence of thermal noise, while Yu et.al.  obtained exponential decay using free boundary conditions. Another important factor affecting $`G(r)`$ is the disorder strength. More systematic study of the decay law is left to future work. In Fig.5 the degree of asymmetry $`A`$ is also shown. With increasing the magnitude of the elastic interaction it approaches to the upper limit $`1/4`$ as expected.

Finally we discuss the effect of random internal stress arising from microscopic heterogeneities in the network structure, which are intrinsic to gels . We restrict our discussion to the case $`\lambda =1`$ with small internal deformations. In the expansion of the elastic free energy with respect to $`𝒖`$ there will arise an additional term,

$`\mathrm{\Delta }F_{el,R}={\displaystyle 𝑑𝒓R_{ij}\frac{u_i}{r_j}},`$ (12)   
where $`R_{ij}`$ is the Gaussian random stress with $`R_{ij}(𝒓)=0`$ and

$`R_{ij}(𝒒)R_{kl}(𝒒)=V_1\delta _{ij}\delta _{kl}+V_2(\delta _{ik}\delta _{jl}+\delta _{il}\delta _{jk}).`$ (13)   
Eliminating the elastic field from $`\mathrm{\Delta }F_{el}+\mathrm{\Delta }F_{el,R}`$ we have a new interaction free energy $`\alpha 𝑑𝒓R_1^SQ_1`$, where $`R_1^S`$ is defined using the shear component $`R_{ij}^S=R_{ij}R_{kk}\delta _{ij}/d`$ by an equation parallel to (8) as

$`R_1^S(𝒒)`$ $`=`$ $`\mathrm{sin}(2\phi )R_{xx}^S(𝒒)\mathrm{cos}(2\phi )R_{xy}^S(𝒒).`$ (14)   
Treating this interaction as a weak perturbation as in , we can see that it renders the equilibrium correlation length finite even in the absence of the disorder free energy (3). We mention that Golubović and Lubensky  discussed another mechanism of long-range-orientational-order breaking due to random stress. Their argument is based on the observation that the amplitude of thermal fluctuations around a uniformly aligned state diverges. Its relevance to the present case of nematic networks is limited in that their free energy does not explicitly include the orientational degree of freedom.

To summarize, we have numerically obtained a soft mechanical reponse during the P-M transition. It originates from structural self-organization of domains due to the long-range elastic interaction, and should be distinguished from the soft elasticity  of uniformly oriented networks. The elastic contribution to the stress is slightly negative in the transition region. We have found a positive disorder contribution to the stress. The elastic interaction is found to increase the correlation length. We have demonstrated that random internal stress acts as a random field on the director. Further experimental and theoretical studies are necessary to examine its relevance to real polydomain textures.

The author gratefully acknowledges Prof. A. Onuki for helpful discussions and a critical reading of the manuscript, and Dr. E. M. Terentjev for valuable comments on our related work.