Vol.0 (20xx) No.0, 000–000

11institutetext: School of Physics and Materials Science, Guangzhou
University, Guangzhou 510006, China<EMAIL_ADDRESS>

# On the nature of long period radio pulsar GPM J1839$-$10: death line and
pulse width

H. Tong

(Received 2009 month day; accepted 2009 month day)

###### Abstract

Recently another long period radio pulsar GPM J1839$-$10 is reported, similar
to GLEAM-X J162759.5$-$523504.3. Previously, the energy budget and rotational
evolution of long period radio pulsars had been considered. This time, the
death line and pulse width for neutron star and white dwarf pulsars are
investigated. The pulse width is included as the second criterion for neutron
star and white dwarfs pulsars. It is found that: (1) PSR J0250+5854 and PSR
J0901$-$4046 etc should be normal radio pulsars. They have narrow pulse width
and they lie near the radio emission death line. (2) The two long period radio
pulsars GLEAM-X J162759.5$-$523504.3 and GPM J1839$-$10 is unlikely to be
normal radio pulsars. Their possible pulse width is relatively large. And they
lie far below the fiducial death line on the $P-\dot{P}$ diagram. (3) GLEAM-X
J162759.5$-$523504.3 and GPM J1839$-$10 may be magnetars or white dwarf radio
pulsars. At present, there are many parameters and uncertainties in both of
these two possibilities.

###### keywords:

stars: magnetar – pulsars: general – pulsars: individual (GPM J1839$-$10)

##  1 Introduction

Pulsars are rotating magnetized neutron stars (Lyne & Graham-Smith 2012).
Their rotation period can range from milliseconds to that of several seconds
(Young et al. 1999). The period of magnetars can be as large as 12 seconds
(Olausen & Kaspi 2014). Recently, there are new discoveries toward the long
period end of pulsars. PSR J0250+5854 is the slowest radio pulsar at that
time, with a period of $23.5\ \rm s$ (Tan et al. 2018). Later more radio
pulsars with periods larger than 10 seconds are discovered (Han et al. 2021).
The rotational evolution of these long period radio pulsars may evolve physics
in the magnetosphere or magnetic field decay (Kou et al. 2019). The central
compact object inside supernova remnant RCW 103 showed magentar activities and
was confirmed to be a magnetar with a period about 6.6 hours (D’Ai et al.
2016; Rea et al. 2016). Its rotational evolution may require a magnetar spin-
down by a fallback disk (Tong et al. 2016). However, the magnetar inside RCW
103 is radio quiet. It is not sure why it is radio quiet, e.g. beaming or no
radio emission at all.

The long period radio pulsar (dubbed as LPRP) is push forward by the discovery
of GLEAM-X J162759.5$-$523504.3 (GLEAM-X J1627 for short, Hurley-Walker et al.
2022). GLEAM-X J1627 has a period of $1091\ \rm s$ (about 18 minutes) and it
is a transient in radio. The existence of LPRPs is further strengthened by PSR
J0901$-$4046 with a period of $76\ \rm s$ (Caleb et al. 2022). However, it is
not certain whether these LPRPs are neutron stars or white dwarfs. For such
slow rotators, it is more natural that they are white dwarfs (Katz 2022; Loeb
& Maoz 2022). At the same time, similar to the case of RCW 103 magnetar, these
LPRPs may also be magnetar+fallback disk systems (Ronchi et al. 2022; Tong
2023). At present, their nature remains unclear (Rea et al. 2022).

A source similar to GLEAM-X J1627 is recently reported, GPM J1839$-$10
(Hurley-Walker et al. 2023). GPM J1839$-$10 has a pulsation period of $1318\
\rm s$ (about 21 minutes). Compared with GLEAM-X J1627, GPM J1839$-$10 has (1)
a more stringent period derivative upper limit $\dot{P}<3.6\times 10^{-13}\
\rm s\ s^{-1}$ (Hurley-Walker et al. 2023), which can place more constraint on
the pulsar death line. (2) its radio emission has lasted for more than 30
years. This can enable later more monitoring.

For the energy budge and rotational evolution of GPM J1839$-$10, similar
calculations can be applied to it as that for GLEAM-X J1627 (Katz 2022; Loeb &
Maoz 2022; Ronchi et al. 2022; Tong 2023). Since their periods are similar (21
minutes verse 18 minutes), the same conclusion may also be applied to GPM
J1839$-$10, i.e., magnetar+fallback disks or white dwarf pulsars. We will not
repeat the calculations here. This time, we focus on (1) the death line for
neutron star and white dwarf pulsars, (2) pulse width of LPRPs. From these two
aspects, we want to discuss the nature of GPM J1839$-$10, which is the most
recent example of LPRPs.

##  2 Model calculations

###  2.1 Death line for neutron star and white dwarf pulsars

The death line for radio emission of normal pulsars and magnetars has been
discussed in Section 2.5 in Tong (2023). The potential drop at a specific
angle across the polar cap is just a classical electrodynamics exercise. For
normal radio pulsars with a dipole magnetic field, the maximum acceleration
potential across the polar cap is (Ruderman & Sutherland 1975; Zhou et al.
2017):

$\Phi_{\rm max}=\frac{B_{p}R^{3}\Omega^{2}}{2c^{2}}\equiv 10^{12}\ \rm V,$ (1)

where $B_{p}$ is the surface magnetic field at the pole region (which is two
times the equatorial value, Lyne & Graham-Smith 2012), $R$ is the stellar
radius, $\Omega$ is the star’s angular velocity. When the maximum acceleration
potential equals $10^{12}\ \rm V$, it is defined as the radio emission death
line. Below the death line, the star is not expected to have radio emission.
The value of $10^{12}\ \rm V$ is just a fiducial value. A rough estimation of
the physics evolved is that (Ruderman & Sutherland 1975): an electron
accelerated in such a potential attains a Lorentz factor about $\gamma\sim
10^{6}$. This electron may emit curvature photons of energy

$h\nu=h\frac{3\gamma^{3}c}{4\pi\rho}\geq 1\ \rm MeV.$ (2)

These curvature photons may be converted to electron-position pairs in strong
magnetic fields. If the electron Lorentz factor (i.e., acceleration potential)
is lower, the curvature photons may have energy less than $1\ \rm MeV$. The
subsequent pair production process can not continue. And this may result in
cease of the radio emission (i.e., radio emission death line).

For magnetar+fallback disk system, the magnetosphere may be modified by (1)
the fallback disk if the disk is still active, (2) the magnetar’s twisted
magnetic fields. If the death line is modified by the fallback disk, the
corotation radius will replace the light cylinder radius as the maximum radial
extent of the field lines. The corresponding maximum acceleration potential
and death line is presented in Eq. (10) in Tong (2023):

$\Phi_{\rm
max,disk}=\frac{B_{p}R^{3}\Omega^{2}}{2c^{2}}\frac{R_{lc}}{R_{co}}\equiv
10^{12}\ \rm V,$ (3)

where $R_{lc}=Pc/(2\pi)$ is the neutron star’s light cylinder radius,
$R_{co}=(GM/4\pi^{2})^{1/3}P^{2/3}$ is the corotation radius. If the death
line is modified by the twist of the field lines, a twisted field line will
result in a larger polar cap and a larger potential drop. The maximum
acceleration potential and death line for a twisted magnetic field is
presented in Eq. (12) in Tong (2023):

$\Phi_{\rm
max,twist}=\frac{B_{p}R^{3}\Omega^{2}}{2c^{2}}\left(\frac{R_{lc}}{R}\right)^{1-n}\equiv
10^{12}\ \rm V,$ (4)

where $n$ is a parameter characterizing the twist of field lines. $n=1$
corresponds to the dipole case (for $n=1$, the equation returns to the dipole
case Eq.(1)). $0<n<1$ corresponds to the twisted dipole case.

For white dwarf pulsars, it is possible that they have similar magnetospheric
precess to that of neutron star pulsars (Zhang & Gil 2005; Katz et al. 2022).
However, there are several changes in the definition of death line (Eq.1) for
white dwarf pulsars.

  1. 1.

For a typical neutron star, the radius is usually set to be $10\ \rm km$. For
a typical white dwarf (figure 5.17 in Camenzind 2007), it has a radius of one
precent the solar radius $0.01\ R_{\odot}$ (the corresponding white dwarf mass
is about $0.8\ M_{\odot}$).

  2. 2.

The torque of a rotating magnetized object can be approximated as magnetic
dipole braking (Xu & Qiao 2001). The star’s magnetic field can be obtained
from the period and period-derivative measurement, which is a crude estimate
of the star’s true magnetic field. Assuming a perpendicular rotator, the
magnetic field is

$B=\sqrt{\frac{3Ic^{3}}{8\pi^{2}R^{6}}P\dot{P}},$ (5)

where $B$ the equatorial magnetic field at the surface (it is two times
smaller than the polar magnetic field), $I$ is the star’s moment of inertia.
For typical neutron stars, with $I\approx 10^{45}\ \rm g\ cm^{2}$, $R\approx
10^{6}\ \rm cm$, it is the commonly cited formula for radio pulsars:

$B=3.2\times 10^{19}\sqrt{P\dot{P}}\ \rm G.$ (6)

For white dwarfs, assuming a typical mass of $0.8\ M_{\odot}$ and radius of
$0.01\ R_{\odot}$, the white dwarf’s moment of inertia is about $3\times
10^{50}\ \rm g\ cm^{2}$. Therefore, the characteristic magnetic field for a
white dwarf pulsar is:

$B=5.2\times 10^{13}\sqrt{P\dot{P}}\ \rm G.$ (7)

This formula will be employed when drawing the death line on the $P-{\dot{P}}$
diagram.

  3. 3.

From Eq. (2), a minimum Lorentz factor (i.e. acceleration potential) is
require to generate curvature photons with energy higher than $1\ \rm MeV$. In
the case of white dwarfs, the stellar radius is larger. The curvature radius
of the magnetic field line is also larger, which is order of $\sqrt{rR_{lc}}$
($r$ is the emission height, Xu & Qiao 2001). Therefore, a higher acceleration
potential may be required, e.g. as high as $10^{13}\ \rm V$ in Eq. (1). The
exact value depends on the detailed modeling of the white dwarf’s
magnetosphere (as in the case of neutron star pulsars).

The death line for normal radio pulsars, magnetars and magnetar+fallback disk
systems is shown in figure 1, along with GPM J1839$-$10 and other LPRPs.
Figure 1 is updated from figure 2 in Tong (2023). The death line for white
dwarf pulsars is so different from that of neutron star case that it is shown
separately in figure 2. $\Phi_{\rm max}=10^{12}\ \rm V$ and $\Phi_{\rm
max}=10^{13}\ \rm V$ are shown respectively. The characteristic magnetic field
for white dwarf pulsars is also shown, for $B=10^{8}\ \rm G$ and $B=10^{9}\
\rm G$. Most of the presently observed pulsating white dwarfs have magnetic
field smaller than $10^{9}\ \rm G$ (Zhang & Gil 2005; Katz 2022; Marsh et al.
2016; Pelisoli et al. 2023).

Figure 1: Definition of death line and distribution of long period radio
pulsars (red circles) on the $P$-$\dot{P}$ diagram. The fiducial pulsar death
line, the death line for a twisted magnetic field ($n=0.8$), and the death
line modified by the fallback disk are also shown. Updated from figure 2 in
Tong (2023).
Figure 2: Definition of death line and contour of constant magnetic field for
white dwarf pulsars. Two death lines are shown, with $\Phi_{\rm max}=10^{13}\
\rm V$ (upper one) and $\Phi_{\rm max}=10^{12}\ \rm V$ (lower one)
respectively. Two contours of constant magnetic field are shown, with
$B=10^{9}\ \rm G$ and $B=10^{8}\ \rm G$ respectively. The limiting Keplerian
period for a typical white dwarf pulsar is also shown.

###  2.2 Pulse width

GPM J1839$-$10 has a period of $1318\ \rm s$. The single pulse vary in a pulse
window $400\ \rm s$ (Hurley-Walker et al. 2023). Similar things also happen in
normal radio pulsars and radio emitting magnetars (Levin et al. 2012; Yan et
al. 2015; Huang et al. 2021). At present, the integrated pulse profile of GMP
J1839$-$10 is not available. From previous experiences in pulsars and
magnetars, the pulse window may be an estimate of the integrated pulse width.
Then the pulse width of GPM J1839$-$10 is $PW\approx 400/1318=30\%$ of the
pulse phase (in this case, the pulse width may also be called the duty cycle,
Tan et al. 2018). If this is the pulse width of GPM J1839$-$10, it can also
constrain the nature of the source.

For normal radio pulsars, the colatitude of the last open field line at
emission height $r$ is:

$\theta_{\rm open}=\sin^{-1}\left(\frac{r}{R_{lc}}\right)^{1/2}.$ (8)

The emission beam radius is $3/2$ times the angle $\theta_{\rm open}$
(Eq.(17.9) in Lyne & Graham-Smith 2012):

$\rho_{\rm beam}=\frac{3}{2}\sin^{-1}\left(\frac{r}{R_{lc}}\right)^{1/2}.$ (9)

The emission height of normal radio pulsars is always less than $100$ times
the neutron star radius (Johnston et al. 2023). Therefore, the emission height
in Eq.(9) can be set as $100R$, where $R$ is the stellar radius. The observed
pulse width depends also on the inclination angle $\alpha$ (angle between
rotation axis and line of sight), and impact angle $\beta$ (closest approach
between the magnetic axis and the line of sight). The impact angle may always
be a small quantity, in this case the observed pulse width is related to the
emission beam radius as (Eq.(15.2) in Lyne & Graham-Smith 2012):

$W=2\rho_{\rm beam}\frac{1}{\sin\alpha}.$ (10)

In terms of pulse phase, the observed pulse width is:

$PW=\frac{W}{2\pi}=\frac{3}{2\pi}\sin^{-1}\left(\frac{r}{R_{lc}}\right)^{1/2}\frac{1}{\sin\alpha}.$
(11)

A small inclination angle $\alpha$ can result in a large pulse width.

For a magnetar+fallback disk system, the calculation is similar to the
calculations for death line. If the magnetosphere is regulated by the fallback
disk, the corotation radius replaces the role of light cylinder radius in Eq.
(11). For a twisted magnetic field, the colatitude of the last open field line
will be larger (Eq. 11 in Tong 2023). The pulse width in units of pulse phase
is:

$PW=\frac{3}{2\pi}\sin^{-1}\left(\frac{r}{R_{lc}}\right)^{n/2}\frac{1}{\sin\alpha}.$
(12)

For $n=1$, the above equation returns to the dipole case (Eq.(11)). For a
white dwarf pulsar, the expression for the pulse width is the same as Eq.
(11), except that the stellar radius should be the white dwarf radius.

The theoretical pulse width for normal pulsars, magnetar+fallback disk
systems, and white dwarf pulsars is shown in figure 3. Magnetars and white
dwarfs pulsars have wider pulse width compared with that of normal radio
pulsars. This is especially true for LPRPs.

Figure 3: Theoretical pulse width as a function of period. The pulse width is
in units of pulse phase. From bottom to top are: normal radio pulsars (black),
magnetars with twisted magnetic field (dashed blue, for $n=0.5$),
magnetar+fallback disk systems (solid blue), and white dwarf pulsars (green).
Since the inclination angle is unknown, the plotted pulse width is actually
$PW\times\sin\alpha$. For a small $\alpha$, the actual pulse width can be
larger.

###  2.3 On the nature of GPM J1839$-$10

From the death line on the $P-\dot{P}$ diagram, the LPRP GPM J1839$-$10 can
not be normal radio pulsars. A twisted magnetic field or the presence of a
fallback disk can help to lower the position of the death line on the
$P-\dot{P}$ diagram. The quantitative result depends on the parameters
involved, e.g. the twist parameter $n$. Figure 1 shows the death line for a
twisted magnetic field with $n=0.8$. For a more twisted magnetic field (i.e.
$n=0.5$), the position of the death line is lower on the $P-\dot{P}$ diagram.
However, the typical untwisting timescale may be smaller, which may have
difficulties in explaining why GPM J1839$-$10 can have radio emission lasting
more than 30 years. A magnetar+fallback disk system may be consistent with the
position of GPM J1839$-$10 on the $P-\dot{P}$ diagram. However, the generation
of radio emission in the presence of an active fallback disk is uncertain,
although there are such possibilities (see Section 2.5 in Tong 2023 for
discussions).

The death line for white dwarf pulsars is also consistent with the position of
GPM J1839$-$10 on the $P-\dot{P}$ diagram. However, there may be three
constrains for white dwarf pulsars: (1) The maximum acceleration potential for
white dwarf pulsars may be higher, e.g., as high as $10^{13}\ \rm V$. (2)
White dwarf pulsars generally have magnetic field lower than $10^{9}\ \rm G$
(Zhang & Gil 2005; Katz 2022; Marsh et al. 2016; Pelisoli et al. 2023). (3)
For a white dwarf with mass $0.8M_{\odot}$ and radius $0.01R_{\odot}$, the
limiting Keplerian period is:

$P_{K}=2\pi\sqrt{\frac{R^{3}}{GM}}\approx 11\ \rm s.$ (13)

These three constrains will limit the existence of white dwarf pulsars to a
small triangle on the $P-\dot{P}$ diagram, see Figure 2. Such a triangle of
parameter space may explain why there are so few white dwarf radio pulsars
(LPRPs are just candidates).

The pulse with of GPM J1839$-$10 is unknown at present. If its pulse window
represents its pulse width, then it will have a pulse width of $30\%$. From
the theoretical pulse width, normal radio pulsars may have difficulties in
explaining the pulse width of GPM J1839$-$10\. The solution in the normal
radio pulsar case is that they should have an extremely large emission height.
For a twisted magnetic field or magnetar+fallback disk system, the pulse width
is larger, typically several precent. A small inclination angle or a slightly
higher emission height may explain the possible observed pulse width of GPM
J1839$-$10\. White dwarf pulsars can naturally have larger pulse width.

By combining the death line and pulse width requirement, it is unlikely that
GPM J1839$-$10 is a normal radio pulsar. It is possible that GPM J1839$-$10 is
a magnetar (including magnetar+fallback disk system) or a white dwarf pulsar.
For these two possibilities, there are many parameters at present. This
conclusion is consistent with previous calculations for GLEAM-X J1627, based
on energy budget and rotational evolution (Katz 2022; Loeb & Maoz 2022; Ronchi
et al. 2022; Tong 2023). Population synthesis of neutron star and white dwarf
pulsars also get similar conclusions (Rea et al. 2023).

##  3 Discussion

Comparison with other LPRPs. All the three LPRPs PSR J2144-3933 (Young et al.
1999), PSR J0250+5854 (Tan et al. 2018) and PSR J0901$-$4046 (Caleb et al.
2022) have a very narrow pulse width, typically less than $1\%$. Assuming a
maximum emission height of $100R$ (Johnston et al. 2023), the theoretical
upper limit on pulse width is (Eq.(11)):

$PW<7\%P^{-1/2}\frac{1}{\sin\alpha}.$ (14)

Therefore, PSR J0250+5854 (with a period of $23.5\ \rm s$) and PSR
J0901$-$4046 (with a period of $76\ \rm s$) etc is consistent with a normal
radio pulsar origin, from the pulse width point of view. However, both of
GLEAM-X J1627 and GPM J1839$-$10 (Hurley-Walker et al. 2022, 2023) showed
possible signatures of a large pulse width. Therefore, we propose that the
pulse with of LPRPs may be taken as the second criterion in identifying their
nature, in addition to their relative position to the death line on the
$P-\dot{P}$ diagram. If a future LPRP has narrow pulse width consistent with
that of normal radio pulsars, it may be viewed as an extreme radio pulsar. All
we have to do is to consider the corresponding magnetospheric physics
involved, e.g. the definition of the death line etc. If a LPRP has a large
pulse width, then we must consider the possibility of magnetars or white dwarf
pulsars.

Comparison with other radio emitting magnetars. Radio emitting magnetars
generally have a large pulse width compared with normal radio pulsars (Levin
et al. 2012; Yan et al. 2015; Huang et al. 2021). Especially, the single pulse
width of magnetars are generally very narrow. They may vary randomly in the
pulse window, therefore resulting in a wide integrated pulse width (Levin et
al. 2012; Yan et al. 2015). We think this may also be the case of GPM
J1839$-$10: narrow single pulse and a wide pulse window. From the pulse width
point of view, the two LPRPs GLEAM-X J1627 and GPM J1839$-$10 are similar to
radio-emitting magnetars.

Maximum period of LPRPs. From Figure 1 and 2 (or Figure 1 in Rea et al. 2023),
the existence of maximum magnetic field in combination with the death line
implies there is a maximum period for radio-emitting neutron stars and white
dwarfs. For neutron stars, the maximum magnetic field may be about $10^{16}\
\rm G$. The definition of death line is rather uncertain (which will result in
a death valley). However, the possible maximum period for radio emission may
be around $10^{4}\ \rm s$ or so. In this respect, the magnetar inside RCW 103
(with a period of 6.6 hours) is not expected to have radio emissions. For
white dwarfs, assuming of a maximum magnetic field of $10^{9}\ \rm G$, the
expected maximum period may be about several thousand seconds. Future LPRPs
with longer periods may help to unveil their nature, i.e., neutron stars or
white dwarfs.

Comparison with previous works. Observationally, GPM J1839$-$10 is similar to
that of GLEAM-X J1627 (Hurley-Walker et al. 2023). It confirms the existence
of LPRPs as sub-population of radio pulsars, which may deserved future more
searchings. Theoretically, the modeling for GPM J1839$-$10 may seem similar to
that of GLEAM-X J1627. Compared with previous theoretical works (Ronchi et al.
2022; especially Tong 2023), the present work has two outstandings: (1) The
death line is extended to the case of magnetar’s twisted magnetic field and
magnetar+fallback disk system in previous works. Here, the pulsar death line
is extended to the white dwarf pulsar case, considering possible physics
involved during the definition of the death line. (2) The pulse width is newly
included in this work. We propose that it may be taken as the second criterion
concerning the nature of LPRPs, in addition to that of timing parameters
(e.g., period, period-derivative, death line etc). By combining the rotational
evolution, energy budget, death line and pulse with, something is already
clear for LPRPs: (1) For PSR J0901$-$4046 (Caleb et al. 2022) etc with period
less than $76\ \rm s$, they can be understood in the normal radio pulsar
domain. (2) For the two LPRPs GLEAM-X J1627 and GPM J1839$-$10 (Hurley-Walker
et al. 2022, 2023), they are unlikeky to be normal radio pulsars. Whether they
are magnetars or white dwarf pulsars remain inconclusive at present.

##  4 Conclusion

By considering the four aspects together (rotational evolution and energy
budget in previous works, death line and pulse width in this work), it is
found that: (1) PSR J0250+5854 (Tan et al. 2018) and PSR J0901$-$4046 (Caleb
et al. 2022) should be normal radio pulsars. They have narrow pulse width and
they lie near the radio emission death line. Further investigations of their
magnetospheric physics is required. (2) The two LPRPs GLEAM-X J1627 and GPM
J1839$-$10 (Hurley-Walker et al. 2022, 2023) is unlikely to be normal radio
pulsars. Their possible pulse width is relatively large. And they lie far
below the fiducial death line on the $P-\dot{P}$ diagram. (3) GLEAM-X J1627
and GPM J1839$-$10 may be (a) magnetars with twisted magnetic field or
magnetar+fallack disk systems, or (b) white dwarf radio pulsars. At present,
there are many uncertainties in both of these two possibilities. More
multiwave observations are required in order to tell whether they are
magnetars or white dwarfs.

## acknowledgments

H. Tong would like to thank Dr. Huang Zhi-Peng and Yan Zhen for discussions on
pulse width and emission height. This work is supported by National SKA
Program of China (No. 2020SKA0120300) and NSFC (12133004).

## References

  * [1] Caleb, M., Heywood, I., Rajwade, K., et al. 2022, NatAs, 6, 828
  * [2] Camenzind, M., 2007, Compact objects in astrophysics, Springer, Berlin
  * [3] D’Ai, A., Evans, P. A., Burrows, P. A., et al. 2016, MNRAS, 463, 2394
  * [4] Han, J. L., Wang, C., Wang, P. F., et al., 2021, RAA, 21, 107
  * [5] Huang, Z. P., Yan, Z., Shen, Z. Q., et al., 2021, MNRAS, 505, 1311
  * [6] Hurley-Walker, N., Zhang, X., Bahramian, A., et al. 2022, Nature, 601, 526
  * [7] Hurley-Walker, N., Rea, N., McSweeney, S. J., et al., 2023, Nature, 619, 487
  * [8] Johnston, S., Kramer, M., Karastergiou, A., et al., 2023, MNRAS, 520, 4801
  * [9] Katz, J. I. 2022, Ap&SS, 367, 108
  * [10] Kou, F. F., Tong, H., Xu, R. X., & Zhou, X. 2019, ApJ, 876, 131
  * [11] Levin, L., Bailes, M., Bates, S. D., et al. 2012, MNRAS, 422, 2489
  * [12] Loeb, A., & Maoz, D. 2022, RNAAS, 6, 27
  * [13] Lyne, A. G., & Graham-Smith, F. 2012, Pulsar astronomy (4th ed.), Cambridge University Press, Cambridge
  * [14] Marsh, T. R., Gansicke, B. T., Hummerich, S., et al., 2016, Nature, 537, 374
  * [15] Olausen, S. A., & Kaspi, V. M. 2014, ApJS, 212, 6
  * [16] Pelisoli, I., Marsh, T. R., Buckley, D. A. H., et al., 2023, arXiv:2306.09272 (NatAs in press)
  * [17] Rea, N., Borehese, A., Esposito, P., et al. 2016, ApJ, 828, L13
  * [18] Rea, N., Coti Zelati, F., Dehman, C., et al. 2022, ApJ, 940, 72
  * [19] Rea, N., Hurley-Walker, N., Pardo-Araujo, D., et al., arXiv:2307.10351
  * [20] Ronchi, M., Rea, N., Graber, V., Hurley-Walker, N. 2022, ApJ, 934, 184
  * [21] Ruderman, M. A., & Sutherland, P. G. 1975, ApJ, 196, 51
  * [22] Tan, C. M., Bassa, C. G., Cooper, S., et al. 2018, ApJ, 866, 54
  * [23] Tong, H., Wang, W., Liu, X. W., & Xu, R. X. 2016, ApJ, 833, 265
  * [24] Tong, H., 2023, ApJ, 943, 3
  * [25] Xu, R. X., Qiao, G. J., 2001, ApJL, 561, L85
  * [26] Yan, Z., Shen, Z. Q., Wu, X. J., et al., 2015, ApJ, 814, 5
  * [27] Young, M. D., Manchester, R. N., Johnston S., 1999, Nature, 400, 848
  * [28] Zhang, B., Gil, J., 2005, ApJL, 631, L143
  * [29] Zhou, X., Tong, H., Zhu, C., & Wang, N. 2017, MNRAS, 472, 2403


# A Molecular Theory for Liquid Crystal Elastomers: Nematic Ordering, Shape
Deformation and Mechanical Response

Luofu Liu  Department of Chemical and Biomolecular Engineering, University of
California Berkeley, Berkeley, California 94720, United States  Rui Wang
<EMAIL_ADDRESS>Department of Chemical and Biomolecular Engineering,
University of California Berkeley, Berkeley, California 94720, United States

###### Abstract

Modeling liquid crystal elastomers (LCEs) at the molecular level is crucial
for the predictable design of energy-conversion and stimuli-responsive
materials. Here, we develop a self-consistent field theory for LCEs which
captures the coupling between nematic ordering, backbone alignment and network
deformation. Molecular features such as density of elastic strands, strength
and architecture of local chemical hinge, and LC grafting density are
systematically included. Crosslinking suppresses nematic ordering as a result
of the elastic energy stored during network deformation. Higher work capacity
can be achieved by less crosslinked LCEs. The spontaneous shape change of end-
on side-chain LCEs can be either elongation or contraction depending on the
competition between the local and global couplings. Adjusting LC grafting
density is found to be an effective way to fine-tune the deformation mode. We
elucidate a universal scaling relationship between the transition temperature
and the shear modulus as $T_{\rm NI,0}-T_{\rm NI}\sim{\mu^{*}}^{\frac{4}{5}}$.
Furthermore, we predict that the first-order nematic phase transition can be
degraded into a continuous manner upon applied stress. Coupled with nematic
ordering, the mechanical response of LCEs significantly deviates from
classical rubber elasticity. A plateau in the stress-deformation curve appears
accompanied by the nematic phase transition. Our theoretical predictions are
in good agreement with the experimental results reported in the literature.

Materials Sciences Division, Lawrence Berkeley National Lab, Berkeley,
California 94720, United States

##  1 Introduction

There is an ever-increasing societal demand for efficient energy storage and
conversion. Liquid crystal elastomers (LCEs) are crosslinked liquid crystal
polymers (LCPs), where spontaneous shape change can be triggered by the
ordering of LC moieties 1. This ordering forces the polymer backbone into a
stretched/compressed conformation which relaxes into the random coil when the
ordering disappears. Such reversible shape change not only enables the energy
conversion from thermal/optical/electrical/chemical forms to elastic energy,
but also makes LCEs promising stimuli-responsive smart materials 2 in the
applications like artificial muscles 3, 4, 5, 6, actuators 7, 8, 9, 10,
sensors 7, 11, 12, micro-pumps 13, and adhesives 9, 14, etc.

The phase and mechanical behaviors of LCEs strongly depend on their molecular
structure and interactions. Compared to the uncrosslinked LCPs, it has been
observed that crosslinking suppresses the nematic ordering by lowering the
transition temperature and reducing the order parameter 15. The first-order
feature of the phase transition can even degrade to a continuous one 16, 15,
17, 18. Recent experiments show a significant dependence of the transition
temperature, actuation strain, and work capacity on the molecular details of
the LCEs, including the LC content, crosslinker functionality and crosslinking
density 19, 20, 21, 22, 23. While mesogens can be incorporated into the
polymer backbone as main-chain LCEs, connecting mesogens to the backbone as a
pendant group through a flexible spacer (side-chain LCEs) facilitates the
versatile design of materials: the network backbone, mesogen, and spacer can
be chosen separately 24, 25. It is interesting that the deformation mode of
side-chain LCEs (SCLCEs) is highly related to the local architecture of the
attachment. The spontaneous deformation of side-on SCLCEs is always
elongation, in contrast, that of end-on SCLCEs can be either elongation or
contraction 26, 27, 28. The response of LCEs to applied stress is more
complicated. It depends on the geometric relationship between the force
direction and the nematic director 29, 30, 31, 32, 33.

Modeling LCEs at the molecular level is crucial for the predictable design of
LCE-based materials. The pioneering theory to describe the phase and
mechanical behaviors of LCEs is developed by Warner, Terentjev and co-workers,
29, 34, 35, 36, 37, 38, 39 which combines the nematic energy in the Landau-de
Gennes form 40 and the neo-classical elastic energy. Such a combination
provides an initial attempt to model the coupling between the nematic ordering
and network deformation. However, it relies on phenomenological parameters
that lack clear physical significance in the Landau energy expansion, limiting
its capability to account for molecular features. An alternative approach is
to replace the Landau-de Gennes form of nematic energy with the Maier-Saupe
form.41, 42, 43, 44, 45, 46, 47, 48, 49 This enables the applications to
systems that have multiple nematic order parameters, e.g., LCEs consisting of
mixture of mesogens and networks swelling in nematic solvents.47, 48, 49 The
Maier-Saupe form also facilitates the description of the response under
various external stimuli, such as light 44, 45, 46, electric and magnetic
fields 42, 47. However, all the existing theories fail to distinguish the
alignment of polymer segments from mesogen orientation, which greatly hinders
their application to SCLCEs. This problem is particularly severe for the end-
on architecture where the local geometry prefers a perpendicular alignment
between mesogens and the backbone. Furthermore, it is difficult for the
existing theories to model complex LCEs with spatial inhomogeneity due to
local defects or microdomains.50 In a previous work, Wang and Wang developed a
molecular model to describe the global and local coupling effects on the phase
behavior and chain conformation in side-chain liquid crystal polymers.51
However, how to systematically incorporate this model into the crosslinked
polymers and properly account for network elasticity remains a great
challenge.

In this work, we develop a self-consistent field theory for LCEs, which is
built upon the previous work by Wang and Wang on liquid crystal polymers.51
The theory systematically captures the coupling between nematic ordering,
backbone alignment, and network deformation. The couplings between backbone
and mesogen are treated separately by a global term associated with nematic
field, and a local term originating from chemical linkage. The structure
information at the molecular level, including crosslinking density, mesogen
grafting density, as well as the geometry and strength of mesogen attachment,
are explicitly incorporated. The phase behavior and shape deformation both in
the spontaneous case and under external stress are investigated in a unified
framework.

##  2 Theory

As illustrated in Fig. 1, the LCE is modeled by a network formed by end-
linking $n$ liquid crystal polymers. Here, we consider a general case of
SCLCEs. Each LCP between two adjacent crosslinkers is described by the freely-
jointed chain model with $N$ repeating sections. Each section contains $m$
polymer backbone segments and one mesogen, with a Kuhn length $b$. The number
of segments for one entire LCP is thus $Nm$. Without loss of generality, the
mesogen is assumed to be grafted on the first segment of the repeating
section. $m$ hence represents the mesogen grafting density where $m=1$ is the
fully grafted case. We use $\bm{p}_{\alpha,1}$, $\bm{p}_{\alpha,2}$, …,
$\bm{p}_{\alpha,Nm}$ and $\bm{q}_{\alpha,1}$, $\bm{q}_{\alpha,2}$, …,
$\bm{q}_{\alpha,(N-1)m+1}$ to denote the unit orientational vectors of the
backbone segment and mesogen in the $\alpha$-th chain. The model recovers
MCLCE in the absence of the grafted mesogens (i.e., all the $\bm{q}$ vectors).

Figure 1: Schematic of modeling liquid crystal elastomers (LCEs). The
elastomer is a polymer network described by the neo-classical model. Each
elastic strand between two crosslinks is modeled as a freely-jointed chain.
Global nematic interactions between mesogens ($\bm{q}$) and backbone segments
($\bm{p}$) are denoted by red arrows. The chemical hinge connecting the
mesogen and its attached backbone segment is considered as the local coupling
(green arrows).

We start by defining the following microscopic operators to describe the local
density and molecular orientation of the backbone segment and mesogen. The
density operators are given by

$\displaystyle\hat{\phi}_{\rm
p}(\bm{r})=v_{\rm{p}}\sum_{\alpha=1}^{n}\sum_{j=1}^{N}\delta(\bm{r}-\bm{r}_{\alpha,(j-1)m+1})$
(1) $\displaystyle\hat{\phi}_{\rm
p^{*}}(\bm{r})=v_{\rm{p}}\sum_{\alpha=1}^{n}\sum_{j=1}^{N}\sum_{k=2}^{m}\delta(\bm{r}-\bm{r}_{\alpha,(j-1)m+k})$
$\displaystyle\hat{\phi}_{\rm
q}(\bm{r})=v_{\rm{q}}\sum_{\alpha=1}^{n}\sum_{j=1}^{N}\delta(\bm{r}-\bm{r}_{\alpha,(j-1)m+1})$

where $\rm p$, $\rm p^{*}$, $\rm q$ refer to the LC-grafted backbone segment,
non-grafted segment, and mesogen, respectively. $v_{\rm p}$ and $v_{\rm q}$
are the occupied volume of backbone segment and mesogen. it should be noted
that the volume of the spacer connecting the mesogen with the backbone is not
accounted for in this work. As discussed later, the spacer will be introduced
as a local hinge effect. The position vector of the $k$-th segment in the
$j$-th repeating section on the $\alpha$-th chain is given by
$\bm{r}_{\alpha,(j-1)m+k}=\bm{R}_{\alpha,0}+b\sum_{k=1}^{(j-1)m+k-1}\bm{p}_{\alpha,k}+(b/2)\bm{p}_{\alpha,(j-1)m+k}$
with $\bm{R}_{\alpha,0}$ the position of one of the chain-end. The nematic
operators are given by

$\displaystyle\hat{\bm{S}}_{\rm
p}(\bm{r})=v_{\rm{p}}\sum_{\alpha=1}^{n}\sum_{j=1}^{N}\delta(\bm{r}-\bm{r}_{\alpha,(j-1)m+1})\left[\bm{p}_{\alpha,(j-1)m+1}\otimes\bm{p}_{\alpha,(j-1)m+1}-\frac{1}{3}\bm{I}\right]$
(2) $\displaystyle\hat{\bm{S}}_{\rm
p^{*}}(\bm{r})=v_{\rm{p}}\sum_{\alpha=1}^{n}\sum_{j=1}^{N}\sum_{k=2}^{m}\delta(\bm{r}-\bm{r}_{\alpha,(j-1)m+k})\left[\bm{p}_{\alpha,(j-1)m+k}\otimes\bm{p}_{\alpha,(j-1)m+k}-\frac{1}{3}\bm{I}\right]$
$\displaystyle\hat{\bm{S}}_{\rm
q}(\bm{r})=v_{\rm{q}}\sum_{\alpha=1}^{n}\sum_{j=1}^{N}\delta(\bm{r}-\bm{r}_{\alpha,(j-1)m+1})\left[\bm{q}_{\alpha,(j-1)m+1}\otimes\bm{q}_{\alpha,(j-1)m+1}-\frac{1}{3}\bm{I}\right]$

where $\bm{I}$ stands for the identity tensor.

Following the treatment by Wang and Wang, the coupling between the chain
backbone and LC side groups is properly accounted for by two types of terms:
the global nematic interaction and the local chemical hinge effect.51 The
Hamiltonian of the system consists of nematic, hinge, and elastic
contributions:

$\mathcal{H}=\mathcal{H}_{\rm nem}+\mathcal{H}_{\rm hinge}+\mathcal{H}_{\rm
ela}$ (3)

The nematic interaction is described by two-body pseudopotentials between
selected pairs of mesogen, grafted backbone segment and non-grafted backbone
segment in the Maier-Saupe type,51, 52, 53 which yields

$\displaystyle\mathcal{H}_{\rm nem}$ $\displaystyle=-\frac{U_{\rm
pp}}{2}\int{\rm d}\bm{r}\hat{\bm{S}}_{\rm p}(\bm{r}):\hat{\bm{S}}_{\rm
p}(\bm{r})-\frac{U_{\rm p^{*}p^{*}}}{2}\int{\rm d}\bm{r}\hat{\bm{S}}_{\rm
p^{*}}(\bm{r}):\hat{\bm{S}}_{\rm p^{*}}(\bm{r})-\frac{U_{\rm qq}}{2}\int{\rm
d}\bm{r}\hat{\bm{S}}_{\rm q}(\bm{r}):\hat{\bm{S}}_{\rm q}(\bm{r})$ (4)
$\displaystyle-U_{\rm pp^{*}}\int{\rm d}\bm{r}\hat{\bm{S}}_{\rm
p}(\bm{r}):\hat{\bm{S}}_{\rm p^{*}}(\bm{r})-U_{\rm pq}\int{\rm
d}\bm{r}\hat{\bm{S}}_{\rm p}(\bm{r}):\hat{\bm{S}}_{\rm q}(\bm{r})-U_{\rm
p^{*}q}\int{\rm d}\bm{r}\hat{\bm{S}}_{\rm p^{*}}(\bm{r}):\hat{\bm{S}}_{\rm
q}(\bm{r})$

where $U_{\kappa\kappa^{\prime}}$ is the nematic interaction parameter between
species $\kappa$ and $\kappa^{\prime}$ ($\kappa,\kappa^{\prime}={\rm p},{\rm
p^{*}},{\rm q}$).

The hinge effect accounts for the local coupling between the LC side group and
its attached chain segment as a result of chemical connection from the spacer.
For simplicity, we assume that the local coupling only exists between the
mesogen and the backbone segment it directly connects to. Its effect on the
adjacent backbones is neglected in this work. The hinge effect is modeled by a
simple quadratic form as

$\mathcal{H}_{\rm
hinge}=-\sum_{\alpha=1}^{n}\sum_{j=1}^{N}u_{C}(\bm{p}_{\alpha,(j-1)m+1}\cdot\bm{q}_{\alpha,(j-1)m+1})^{2}$
(5)

where $u_{C}$ is the coupling parameter. The magnitude of $u_{C}$
characterizes the strength of the local coupling, which is determined by
factors like spacer length and stiffness. The sign of $u_{C}$ reflects the
geometry of the alignment between the mesogen ($\bm{q}_{\alpha,(j-1)m+1}$) and
the backbone ($\bm{p}_{\alpha,(j-1)m+1}$). $u_{C}>0$ represents the case of
side-on SCLCEs where the mesogen prefers to align parallel to its attached
backbone. $u_{C}<0$ represents the case of end-on SCLCEs where the
perpendicular alignment is energetically favorable.

We consider the LCE as a uniformly crosslinked ideal network where topological
defects such as dangling ends and cyclic loops are ignored. The elastic energy
is described by the neo-classical model 1, given by

$\mathcal{H}_{\rm ela}=\frac{1}{2}\mu\int{\rm d}\bm{r}{\rm
Tr}[\bm{l}_{0}\bm{\lambda}^{\rm
T}(\bm{r})\hat{\bm{l}}^{-1}(\bm{r})\bm{\lambda}(\bm{r})]$ (6)

where $\mu$ is the shear modulus proportional to the density of the elastic
strands. $\bm{l}_{0}$ and $\hat{\bm{l}}$ are the step length tensor for the
as-prepared state and deformed state, respectively. $\bm{\lambda}$ is the
deformation gradient tensor. $\hat{\bm{l}}$ can be related to the orientations
of backbone segments as

$\hat{\bm{l}}(\bm{r})=b\left\\{\bm{I}+\frac{1}{m}\left[\frac{3\hat{\bm{S}}_{\rm
p}(\bm{r})}{\hat{\phi}_{\rm
p}(\bm{r})}\right]+\frac{m-1}{m}\left[\frac{3\hat{\bm{S}}_{\rm
p^{*}}(\bm{r})}{\hat{\phi}_{\rm p^{*}}(\bm{r})}\right]\right\\}$ (7)

The prefactors $3/\hat{\phi}_{\kappa}$ are introduced for the consistency with
the standard definition of order parameter tensors.1 The terms $1/m$ and
$(m-1)/m$ weight the contributions from grafted and non-grafted backbone
segments to the step length based on their number fractions.

For LCEs under an applied nominal stress $\bm{\sigma}_{\rm n}$, the partition
function of the system is then calculated by integrating all configurations of
segments and mesogens as well as different deformations:

$\displaystyle Z$ $\displaystyle=\int{\rm
D}\bm{\lambda}\exp\left[\beta\int{\rm d}\bm{r}\bm{\sigma}_{\rm
n}:\bm{\lambda}(\bm{r})\right]\prod_{\alpha=1}^{n}\left(\int{\rm
d}\bm{R}_{\alpha,0}\prod_{j=1}^{N}\prod_{k=1}^{m}\int{\rm
d}\bm{p}_{\alpha,(j-1)m+k}\cdot\prod_{j=1}^{N}\int{\rm
d}\bm{q}_{\alpha,(j-1)m+1}\right)$ (8)
$\displaystyle\prod_{\bm{r}}\delta\left[\hat{\phi}_{\rm
p}(\bm{r})+\hat{\phi}_{\rm p^{*}}(\bm{r})+\hat{\phi}_{\rm
q}(\bm{r})-1\right]\cdot\prod_{\bm{r}}\delta\left[\det\bm{\lambda}(\bm{r})-1\right]\exp(-\beta\mathcal{H})$

The functional integral over $\bm{\lambda}$ results from the transform of
ensemble from constant deformation to constant applied stress. The first
$\delta$-functional enforces the local incompressibility, whereas the second
one originates from the volume conservation of incompressible materials under
deformation.

We follow the standard self-consistent field approach which involves (1)
identity transformation to decouple the many-body interactions into a single
chain interacting with fluctuating fields and (2) replacing the functional
integral with its saddle-point value 54. The detailed derivation is provided
in Sec. I of Supporting Information. The theory is general for LCEs with
spatial inhomogeneity in both composition and nematic ordering. For a simple
case of homogeneous networks, the self-consistent field theory ultimately
leads to the computation of the single-chain partition function in spatially
homogeneous fields $\bm{M}_{\kappa}$ as

$\displaystyle Q$ $\displaystyle=\prod_{j=1}^{N}\prod_{k=1}^{m}\int{\rm
d}\bm{p}_{(j-1)m+k}\cdot\prod_{j=1}^{N}\int{\rm
d}\bm{q}_{(i-1)m+1}\exp\left[\sum_{j=1}^{N}\beta
u_{C}(\bm{p}_{(j-1)m+1}\cdot\bm{q}_{(j-1)m+1})^{2}\right]$ (9)
$\displaystyle\exp\bigg{\\{}\sum_{j=1}^{N}\left[\beta v_{\rm p}\bm{M}_{\rm
p}:\left(\bm{p}_{(j-1)m+1}\otimes\bm{p}_{(j-1)m+1}-\frac{1}{3}\bm{I}\right)\right]$
$\displaystyle+\sum_{j=1}^{N}\sum_{k=2}^{m}\left[\beta v_{\rm p}\bm{M}_{\rm
p^{*}}:\left(\bm{p}_{(j-1)m+k}\otimes\bm{p}_{(j-1)m+k}-\frac{1}{3}\bm{I}\right)\right]$
$\displaystyle+\sum_{j=1}^{N}\left[\beta v_{\rm q}\bm{M}_{\rm
q}:\left(\bm{q}_{(j-1)m+1}\otimes\bm{q}_{(j-1)m+1}-\frac{1}{3}\bm{I}\right)\right]\bigg{\\}}$

where

$\displaystyle\bm{M}_{\rm p}=U_{\rm pp}\bm{S}_{\rm p}+U_{\rm
pp^{*}}\bm{S}_{\rm p^{*}}+U_{\rm pq}\bm{S}_{\rm q}+\frac{3\mu}{2m\phi_{\rm
p}}b\bm{l}^{-1}\bm{\lambda}\bm{l}_{0}\bm{\lambda}^{\rm T}\bm{l}^{-1}$ (10a)
$\displaystyle\bm{M}_{\rm p^{*}}=U_{\rm pp^{*}}\bm{S}_{\rm p}+U_{\rm
p^{*}p^{*}}\bm{S}_{\rm p^{*}}+U_{\rm p^{*}q}\bm{S}_{\rm
q}+\frac{3(m-1)\mu}{2m\phi_{\rm
p^{*}}}b\bm{l}^{-1}\bm{\lambda}\bm{l}_{0}\bm{\lambda}^{\rm T}\bm{l}^{-1}$
(10b) $\displaystyle\bm{M}_{\rm q}=U_{\rm pq}\bm{S}_{\rm p}+U_{\rm
p^{*}q}\bm{S}_{\rm p^{*}}+U_{\rm qq}\bm{S}_{\rm q}$ (10c)

$\bm{S}_{\kappa}$ is the tensorial order parameter corresponding to the
nematic operators. The equilibrium values of $\bm{S}_{\kappa}$ and
$\bm{\lambda}$ can be obtained from the saddle-point approximation, which are
given by

$\displaystyle\bm{S}_{\rm p}=\frac{nv_{\rm
p}}{V}\left\langle\sum_{j=1}^{N}\left(\bm{p}_{(j-1)m+1}\otimes\bm{p}_{(j-1)m+1}-\frac{1}{3}\bm{I}\right)\right\rangle$
(11a) $\displaystyle\bm{S}_{\rm p^{*}}=\frac{nv_{\rm
p}}{V}\left\langle\sum_{j=1}^{N}\sum_{k=2}^{m}\left(\bm{p}_{(j-1)m+k}\otimes\bm{p}_{(j-1)m+k}-\frac{1}{3}\bm{I}\right)\right\rangle$
(11b) $\displaystyle\bm{S}_{\rm q}=\frac{nv_{\rm
q}}{V}\left\langle\sum_{j=1}^{N}\left(\bm{q}_{(j-1)m+1}\otimes\bm{q}_{(j-1)m+1}-\frac{1}{3}\bm{I}\right)\right\rangle$
(11c) $\displaystyle\bm{l}_{0}\bm{\lambda}^{\rm T}\bm{l}^{-1}+\bm{l}_{0}^{\rm
T}\bm{\lambda}^{\rm T}\bm{l}^{-\rm
T}-\gamma\bm{\lambda}^{-1}-2\mu^{-1}\bm{\sigma}_{\rm n}=\bm{0}$ (11d)

where $\langle...\rangle$ is the ensemble average evaluated based on the
Boltzmann factor in the single-chain partition function in Eq. 9. $\gamma$ in
Eq. 11d comes from Lagrangian multiplier as a result of volume conservation.
In principle, Eqs. 11 can describe LCEs with the nematic ordering and applied
stress in any direction. In this work, we focus on the case of uniaxial order
and deformation along the $z$-axis. Both $\bm{\lambda}$ and $\bm{\sigma}_{\rm
n}$ are thus diagonal: $\lambda_{zz}=\lambda$,
$\lambda_{xx}=\lambda_{yy}=1/\sqrt{\lambda}$, and $\sigma_{{\rm
n},zz}=\sigma_{\rm n}$, $\sigma_{{\rm n},xx}=\sigma_{{\rm n},yy}=0$. The
tensorial order parameters are diagonal and traceless, with their scaler forms
defined by normalizing the $zz$-component by the volume fractions:

$\displaystyle s_{\rm p}=\frac{3}{2}\frac{\bm{S}_{{\rm p},zz}}{\phi_{\rm
p}}=\frac{1}{N}\left\langle\sum_{j=1}^{N}P_{2}\left(\cos\theta_{{\rm
p},(j-1)m+1}\right)\right\rangle$ (12a) $\displaystyle s_{\rm
p^{*}}=\frac{3}{2}\frac{\bm{S}_{{\rm p^{*}},zz}}{\phi_{\rm
p^{*}}}=\frac{1}{N(m-1)}\left\langle\sum_{j=1}^{N}\sum_{k=2}^{m}P_{2}\left(\cos\theta_{{\rm
p^{*}},(j-1)m+k}\right)\right\rangle$ (12b) $\displaystyle s_{\rm
q}=\frac{3}{2}\frac{\bm{S}_{{\rm q},zz}}{\phi_{\rm
q}}=\frac{1}{N}\left\langle\sum_{j=1}^{N}P_{2}\left(\cos\theta_{{\rm
q},(j-1)m+1}\right)\right\rangle$ (12c)

where $P_{2}(x)=3x^{2}/2-1/2$ is the second Legendre polynomial, and
$\theta_{\kappa,i}$ is the angle between the $i$-th segment (of species
$\kappa$) and the nematic director. The scaler form of Eq. 11d in uniaxial
deformation becomes

$\sigma_{\rm n}=\mu\left(\lambda\frac{1+2\bar{s}_{\rm p,0}}{1+2\bar{s}_{\rm
p}}-\frac{1}{\lambda^{2}}\frac{1-\bar{s}_{\rm p,0}}{1-\bar{s}_{\rm p}}\right)$
(13)

where $\bar{s}_{\rm p}=(1/m)s_{\rm p}+(1-1/m)s_{\rm p^{*}}$ is the average
nematic order parameter of the backbone segments, and $\bar{s}_{\rm p,0}$ is
that at the as-prepared state. It should be noted that Eq. 13 recovers the
classical stress-deformation relationship for an isotropic rubber,
$\sigma_{\rm n}=\mu(\lambda-1/\lambda^{2})$, if nematic ordering is absent
($\bar{s}_{\rm p}=\bar{s}_{\rm p,0}=0$). Furthermore, Eq. 13 reduces to
$\lambda^{3}=[(1-\bar{s}_{\rm p,0})(1+2\bar{s}_{\rm p})]/[(1+2\bar{s}_{\rm
p,0})(1-\bar{s}_{\rm p})]$ for the force-free case, which is consistent with
the equation for spontaneous deformation used in the literature 1, 37, 36.

The freely-jointed chain model together with the self-coupling form of the
hinge effect allows a further decoupling of the single-chain partition
function $Q$ down to the level of individual segments:
$Q=(Q_{1}Q_{2}^{m-1})^{N}$. $Q_{1}$ and $Q_{2}$ are the partition functions of
the LC-grafted segment and the non-grafted segment respectively, which are
given by

$\displaystyle Q_{1}$ $\displaystyle=\int{\rm d}\bm{p}\int{\rm
d}\bm{q}\exp\left[\beta u_{C}(\bm{p}\cdot\bm{q})^{2}\right]$ (14)
$\displaystyle\exp\left\\{\frac{m+r}{1+r}\left[\beta u_{\rm pp}\phi_{\rm
p}^{2}s_{\rm p}+\beta u_{\rm pp^{*}}\phi_{\rm p}\phi_{\rm p^{*}}s_{\rm
p^{*}}+\beta u_{\rm pq}\phi_{\rm p}\phi_{\rm q}s_{\rm
q}+E\right]P_{2}(\cos\theta_{\rm p})\right\\}$
$\displaystyle\exp\left\\{\frac{m+r}{1+r}\left[\beta u_{\rm pq}\phi_{\rm
p}\phi_{\rm q}s_{\rm p}+\beta u_{\rm p^{*}q}\phi_{\rm p^{*}}\phi_{\rm q}s_{\rm
p^{*}}+\beta u_{\rm qq}\phi_{\rm q}^{2}s_{\rm q}\right]P_{2}(\cos\theta_{\rm
q})\right\\}$

and

$Q_{2}=\int{\rm d}\bm{p}^{*}\exp\left\\{\frac{m+r}{1+r}\left[\frac{\beta
u_{\rm pp^{*}}\phi_{\rm p}\phi_{\rm p^{*}}s_{\rm p}+\beta u_{\rm
p^{*}p^{*}}\phi_{\rm p^{*}}^{2}s_{\rm p^{*}}+\beta u_{\rm p^{*}q}\phi_{\rm
p^{*}}\phi_{\rm q}s_{\rm q}}{m-1}+E\right]P_{2}(\cos\theta_{\rm
p^{*}})\right\\}$ (15)

$E$ in Eqs. 14 and 15 is the elastic contribution:

$E=\frac{\mu^{*}}{m}\left(\frac{\lambda^{2}(1+2\bar{s}_{\rm
p,0})}{(1+2\bar{s}_{\rm p})^{2}}-\frac{1-\bar{s}_{\rm
p,0}}{\lambda(1-\bar{s}_{\rm p})^{2}}\right)$ (16)

with $r=v_{\rm q}/v_{\rm p}$ the volume ratio between mesogen and backbone
segment. $u_{\kappa\kappa^{\prime}}=(2/3)U_{\kappa\kappa^{\prime}}(v_{\rm
p}+v_{\rm q})$ is the rescaled nematic interaction parameter.
$\mu^{*}=\beta\mu(v_{\rm p}+v_{\rm q})$ is the dimensionless shear modulus
which is only proportional to the density of elastic strands but independent
of temperature. The elastic effect is systematically incorporated in the
single-segment partition functions and thus implicitly alters the chain
conformation and nematic ordering through Eq. 12. Furthermore, the three
volume fractions can be expressed by two structure parameters $m$ and $r$ as
$\phi_{\rm p}=1/(m+r)$, $\phi_{\rm q}=r/(m+r)$ and $\phi_{\rm
p^{*}}=(m-1)/(m+r)$. It is worth noting that Eqs. 14 and 15 present the
molecular nature of our theory: partition functions are determined by
molecular interaction parameters ($u_{\kappa\kappa^{\prime}}$ and $u_{C}$) as
well as molecular structural parameters ($m$, $r$ and $\mu^{*}$).

The order parameters in Eq. 12 can also be evaluated at the segment level as

$s_{\kappa}=\langle P_{2}(\cos\theta_{\kappa})\rangle,\ \ \kappa={\rm
p,p^{*},q}$ (17)

where the average for ${\rm p}$ and ${\rm q}$ is taken based on $Q_{1}$, and
that for ${\rm p^{*}}$ is based on $Q_{2}$. By solving the self-consistent
equations 13 and 17 iteratively, we can obtain the equilibrium value of the
nematic ordering, segment alignment, and network deformation under a given
applied uniaxial stress (see Sec. II of Supporting Information for numerical
details). The free energy of one repeating section is given by

$\displaystyle\beta F_{\rm S}$
$\displaystyle=\frac{m+r}{1+r}\bigg{(}\frac{1}{2}\beta u_{\rm pp}\phi_{\rm
p}^{2}s_{\rm p}^{2}+\frac{1}{2}\beta u_{\rm p^{*}p^{*}}\phi_{\rm
p^{*}}^{2}s_{\rm p^{*}}^{2}+\frac{1}{2}\beta u_{\rm qq}\phi_{\rm q}^{2}s_{\rm
q}^{2}+\beta u_{\rm pp^{*}}\phi_{\rm p}\phi_{\rm p^{*}}s_{\rm p}s_{\rm p^{*}}$
(18) $\displaystyle+\beta u_{\rm pq}\phi_{\rm p}\phi_{\rm q}s_{\rm p}s_{\rm
q}+\beta u_{\rm p^{*}q}\phi_{\rm p^{*}}\phi_{\rm q}s_{\rm p^{*}}s_{\rm
q}\bigg{)}-\ln Q_{1}-(m-1)\ln Q_{2}$
$\displaystyle+\frac{m+r}{1+r}\frac{\mu^{*}}{2}\left[\lambda^{2}(1+2\bar{s}_{\rm
p,0})\frac{1+4\bar{s}_{\rm p}}{(1+2\bar{s}_{\rm
p})^{2}}+\frac{2}{\lambda}(1-\bar{s}_{\rm p,0})\frac{1-2\bar{s}_{\rm
p}}{(1-\bar{s}_{\rm p})^{2}}-\frac{2\sigma_{\rm n}}{\mu}\lambda\right]$

It should be noted that Eqs. 13 and 17 can also obtained via minimization of
$F_{\rm S}$ with respect to $s_{\kappa}$ and $\lambda$.

Although the equations are presented in detail for ideal homogeneous networks
under uniaxial stress along the director, our self-consistent field theory is
universal for LCEs with spatial inhomogeneity and for any form of stress. The
theory can easily incorporate topological defects in real networks like
dangling ends and cyclic loops.55, 56, 57 The theory can also be generalized
to the swelling of LCEs in either isotropic or nematic solvents,58, 47, 48, 49
and the response of LCEs to other external stimuli such as optical, electric
and magnetic.44, 45, 46, 42, 47

##  3 Results and Discussion

As discussed in the above section, the macroscopically measurable properties,
nematic ordering ($s_{\rm q}$) and network deformation ($\lambda$) under
stress $\sigma_{\rm n}$ can be calculated based on the molecular information
of interactions ($u_{\kappa\kappa^{\prime}}$ and $u_{C}$) and structures ($m$,
$r$ and $\mu^{*}$). The LCE theory recovers the main-chain case by ignoring
the mesogen terms, i.e., setting $r=0$ and $u_{C}=0$. The case of spontaneous
deformation is naturally recovered when $\sigma_{\rm n}=0$. For simplicity, we
here consider networks prepared at the isotropic state such that $s_{\rm
p,0}=0$.

###  3.1 Effects of the Network Structure on Nematic Ordering and Spontaneous
Deformation

In this subsection, we focus on the interplay between the nematic ordering and
spontaneous shape change in the absence of applied stress. The molecular
nature of our theory facilitates the systematic examination of different
aspects of the network structure, including the density of elastic strands,
the strength and architecture of chemical hinge, and the LC grafting density.
Figure 2 plots the nematic order parameter of mesogen $s_{\rm q}$ and network
deformation $\lambda$ as a function of temperature for different densities of
elastic strands in side-on SCLCEs. The density of elastic strands is
represented by the dimensionless shear modulus $\mu^{*}$, where $\mu^{*}=0$
stands for the special case of uncrosslinked LCP. As $\mu^{*}$ increases, Fig.
2a shows that the nematic-isotropic transition shifts to a lower temperature,
while the order parameter in the nematic phase is also reduced. This trend
predicted by the theory is consistent with the experimental observation by
Disch et al.15 The suppression of nematic ordering by crosslinking originates
from the requirement to store extra energy in the elastic form. For SCLCEs,
this suppression is realized through the coupling between the backbone and
mesogen, by means of both the global nematic field and the local chemical
hinge. The coupling effect critical for SCLCEs cannot be captured by existing
theories.

Figure 2: Effect of density of elastic strands on the mesogen ordering and
spontaneous shape change under the force-free condition ($\sigma_{\rm n}=0$).
(a) The nematic order parameter of mesogen $s_{\rm q}$ and (b) network
deformation $\lambda$ are plotted as a function of rescaled temperature for
various dimensionless shear moduli $\mu^{*}$. The temperature is in a rescaled
form as $T^{*}=k_{\rm B}T/u_{\rm qq}$. $m=1$, $r=1$, $u_{\rm pp}=0.07u_{\rm
qq}$, $u_{\rm pq}=0.2u_{\rm qq}$, and $u_{C}=0.25u_{\rm qq}$.

The isotropic-to-nematic transition is accompanied by the spontaneous shape
change of the elastomer as shown in Fig. 2b. The deformation becomes smaller
as $\mu^{*}$ increases, because the network with a higher modulus has a larger
resistance to deform. Our results suggest that less crosslinked LCE is
preferred to achieve higher work capacity after returning to its isotropic
state due to larger spontaneous shape change. The effect of density of elastic
strands on the LCE work capacity predicted by our theory is also in good
agreement with the experimental findings by Saed et al.21

It has been observed in experiments that how the mesogen is attached to the
backbone in LCEs has a profound effect on both the phase behavior and shape
change.26, 28 While the spontaneous shape change of side-on SCLCEs is always
elongation, that of the end-on SCLCEs can be either elongation or contraction
depending on the chemical hinge and other structural parameters. The chemical
hinge effect is modeled by the local coupling term in our theory, where the
coupling parameter $u_{C}$ characterizes the geometry and the strength of the
hinge. $u_{C}>0$ stands for the side-on attachment whereas $u_{C}<0$ stands
for the end-on case. Figure 3 depicts the local coupling effect on the mesogen
ordering and network deformation. As shown in Fig. 3a, the nematic ordering of
mesogen is weakened with a corresponding reduction of the transition
temperature, when the hinge geometry changes from side-on to end-on. This is
because the preferable perpendicular alignment of the mesogen and backbone in
the end-on geometry leads to a disturbance to the global nematic field.
Furthermore, Fig. 3b shows that the elongational shape change of the side-on
SCLCE also becomes less pronounced as $u_{C}$ decreases and eventually turns
to contraction (manifested by $\lambda<1$) when $u_{C}$ gets more negative. At
$u_{C}=-0.08u_{\rm qq}$, it is interesting to note that the end-on SCLCE does
not exhibit any spontaneous shape change in the entire temperature range even
though nematic ordering occurs at low temperatures. For this critical $u_{C}$,
the local hinge effect exactly offsets the global nematic field, which leads
to effective decoupling of the backbone segments from mesogens.

Figure 3: Effect of the local chemical hinge on the mesogen ordering and
spontaneous shape change. (a) $s_{\rm q}$ and (b) $\lambda$ are plotted as a
function of $T^{*}$ for different values of local coupling parameter $u_{C}$.
$m=1$, $r=1$, $u_{\rm pp}=0.07u_{\rm qq}$, $u_{\rm pq}=0.2u_{\rm qq}$, and
$\mu^{*}=0.5$.

The underlying mechanism that end-on SCLCEs can undergo either elongation or
contraction is a result of the competition between the global and local
couplings. While the global nematic coupling prefers the parallel alignment of
both mesogen and backbone along the director, the local coupling restrained by
the chemical hinge favors a perpendicular alignment. For a small negative
value of $u_{C}$ (loosely hinged end-on), the global coupling dominates,
leading to elongation. On the other hand, the local coupling becomes more
important as $u_{C}$ gets more negative (tightly hinged end-on), which
enforces a contractive deformation. Our theory can well capture the trend
found by Xu et al. that the spontaneous elongation of LCEs changes to
contraction by gradually replacing side-on mesogens with end-on ones.28 Our
theory can also explain the fact that LCE with the same end-on mesogen but
different network structures can undergo either elongation or contraction.26

Figure 4: The reduction of nematic-isotropic transition temperature as a
function of dimensionless shear modulus for MCLCE and SCLCEs with different
hinge effects. $T_{\rm NI,0}^{*}$ is the rescaled transition temperature for
uncrosslinked LCPs. The results are presented in the log-log plot, leading to
a universal exponent $4/5$. For SCLCEs, $m=1$, $r=1$, $u_{\rm pp}=0.07u_{\rm
qq}$, $u_{\rm pq}=0.2u_{\rm qq}$.

Figure 4 summarizes the crosslinking-induced reduction of the transition
temperature $T_{\rm NI}$ from that of the uncrosslinked LCP ($T_{\rm NI,0}$).
The results of both MCLCE and SCLCEs with different hinge effects are
presented. The reduction of $T_{\rm NI}$ is most significant in MCLCE because
mesogens are directly incorporated into the backbone and hence have the
largest impact on network deformation. For SCLCEs, The reduction of $T_{\rm
NI}$ shows a non-monotonic dependence on the local coupling parameter $u_{C}$
due to the competition between the global and local coupling effects. In the
limiting case of loosely hinged end-on SCLCE (e.g. $u_{C}/u_{\rm qq}=-0.05$),
the two coupling effects almost cancel each other, resulting in comparably
small deformation (see Fig. 3b) and negligible reduction of $T_{\rm NI}$.
Although the reduction of $T_{\rm NI}$ is different for different LCE
structures, it follows a universal scaling relationship with shear modulus as

$T_{\rm NI,0}-T_{\rm NI}({\mu^{*}})\sim{\mu^{*}}^{\frac{4}{5}}$ (19)

Using the Landau-de Gennes form of nematic energy, Warner and Gelling
predicted a different linear scaling as $T_{\rm NI,0}-T_{\rm
NI}({\mu^{*}})\sim{\mu^{*}}$.37 We think the $1/5$ difference in the scaling
exponent comes from the fact that the Landau-de Gennes treatment only includes
the leading-order terms of the nematic order parameter. Furthermore, the
scaling relationship predicted by our theory facilitates access to LCEs with a
targeted transition temperature. The complete information of transition
temperature can be obtained by only measuring that of a few samples.

###  3.2 Multiple Deformation Modes by Adjusting LC Grafting Density

The competition between the global and local couplings in SCLCEs and its
resulting shape change can be modulated by varying LC grafting density. In our
freely jointed chain model, LC grafting density is represented by the
parameter $m$ which controls the interval between two adjacent mesogens in the
same chain. Figure 5a shows the effect of LC grafting density on the
spontaneous shape change of end-on SCLCE ($u_{C}=-0.25u_{\rm qq}$). Both the
fully grafted case ($m=1$) and the sparsely grafted case ($m=3$) exhibit a
single deformation mode: the former is always contraction due to the dominant
local coupling effect, whereas the latter is always elongation controlled by
the global coupling. For the case of intermediate LC grafting density ($m=2$),
it is interesting to note that the LCE has dual deformation modes. The LCE
contracts as temperature decreases just below $T_{\rm NI}$, which is followed
by a subsequent re-elongation as temperature further decreases. Such non-
monotonic shape change sheds light on potential applications of LCEs as smart
actuators with multiple responses to external stimuli and shape memory
materials with diverse shape patterns.9

The multiple deformation modes modulated by LC grafting density originate from
the fact that the coupling effects experienced by the LC-grafted (${\rm p}$)
and non-grafted (${\rm p^{*}}$) backbone segments are different. While the
global coupling associated with the nematic field affects both types of
backbone segments, the local coupling due to the chemical hinge only affects
the directly grafted ones. As shown in Fig 5b, non-grafted segments are
aligned parallel to the director ($s_{\rm p^{*}}>0$), whereas the LC-grafted
ones prefer a perpendicular alignment ($s_{\rm p}<0$) due to the tight hinge.
Compared to $s_{\rm p}$, $s_{\rm p^{*}}$ is more sensitive to the decrease of
temperature, particularly showing a more rapid change at low temperatures.
This is because the parallel alignment of $\rm p^{*}$ concentrates only in one
direction; however, the alignment of $\rm p$ is uniformly distributed in the
entire plane perpendicular to the director. The overall effect is a non-
monotonic change in the average backbone alignment $\bar{s}_{\rm p}$ from
negative to positive, leading to the dual deformation modes as indicated by
Eq. 13.

Figure 5: Effect of LC grafting density on the spontaneous shape change for
end-on SCLCE ($u_{C}=-0.25u_{\rm qq}$). LC grafting density is represented by
$m$, the interval between two adjacent mesogens in the same chain. (a)
Deformation $\lambda$ as a function of temperature for different $m$.
Temperature is rescaled relative to the transition temperature $T_{\rm NI}$
for each $m$. (b) Order parameters of backbone alignment for LC-grafted
($s_{\rm p}$) and non-grafted ($s_{\rm p^{*}}$) backbone segments, as well as
the average value ($\bar{s}_{\rm p}$) for the case of $m=2$. $r=1$, $u_{\rm
pp}=u_{\rm p^{*}p^{*}}=u_{\rm pp^{*}}=0.07u_{\rm qq}$, $u_{\rm pq}=u_{\rm
p^{*}q}=0.2u_{\rm qq}$, and $\mu^{*}=0.5$.

Recent experiments showed a significant impact of liquid crystal loading on
the phase transition temperature and actuation strain.20, 22, 23 By replacing
part of mesogens on the backbone of MCLCE with non-mesogenic molecules, Barens
et al. found a linear reduction of the transition temperature with respect to
the LC fraction.20 To validate our theory, we compare the theoretical
predictions with the experimental results of Barens et al.20 As shown in Fig.
6, our theory quantitatively captures the experimental data using
$\mu^{*}=0.3$. The reduction of $T_{\rm NI}$ as $m$ increases is due to the
dilution of mesogens by non-mesogenic molecules.

Figure 6: Effect of the LC fraction on the nematic-isotropic transition
temperature for MCLCE. The LC fraction is represented by $m$, where $m=1$
denotes the LC fully loaded case without any non-mesogenic segment. Our
theoretical predictions using $\mu^{*}=0.3$ are compared to experimental data
adopted from ref. [20].

###  3.3 Deformation and Critical Behavior under Applied Stress

LCE-based materials such as artificial muscles, often work under the condition
of applied external stress. Such applied stress deforms the elastic network,
changes the alignment of backbone segments, and hence alters the nematic
ordering of SCLCEs via the local and global coupling effects. The nematic
ordering in turn affects the deformation, leading to the deviation from
classical rubber elasticity. Our self-consistent field theory enables a
systematic inclusion of these effects in a unified framework. Figure 7 shows
the nematic order parameter of mesogens $s_{\rm q}$ and deformation $\lambda$
for side-on SCLCEs ($u_{C}=0.25u_{\rm qq}$) under different nominal stress
$\sigma_{\rm n}$. As shown in Fig. 7a, nematic ordering is enhanced by the
elongational stress ($\sigma_{\rm n}>0$) along the director, manifested by the
increase of both transition temperature and order parameter. In contrast, the
nematic ordering is suppressed by the compressive stress ($\sigma_{\rm n}<0$).
Compared to the force-free counterpart, the mesogens in SCLCE show pre-
ordering at temperatures higher than the transition point, which is known as
the paranematic phase.59, 18 Furthermore, the two order parameters belonging
to the coexistent phases at the transition point approach each other as
$\sigma_{\rm n}$ increases. At the critical stress, they eventually merge into
a single point. Beyond that point, the transition from the paranematic to
nematic phase occurs in a continuous manner, i.e., supercritical behavior. The
trend of $\lambda$ shown in Fig. 7b follows a similar behavior as $s_{\rm q}$.

Figure 7: Effect of applied nominal stress along the director on (a) mesogen
nematic ordering $s_{\rm q}$ and (b) network deformation $\lambda$ for side-on
SCLCEs ($u_{C}=0.25u_{\rm qq}$). $\sigma_{\rm n}^{*}=\sigma_{\rm n}/(2/3U_{\rm
qq})$ is the dimensionless nominal stress. The dotted lines denote the binodal
envelope. $m=1$, $r=1$, $u_{\rm pp}=0.07u_{\rm qq}$, $u_{\rm pq}=0.2u_{\rm
qq}$, and $\mu^{*}=0.5$.

It has been observed in experiments that the first-order feature of the
nematic phase transition in uncrosslinked LCP can degrade into a continuous
one upon crosslinking.1, 15, 17, 18 Combining the results in Fig. 2 and Fig.
7, we demonstrate that the crosslinking alone will not change the first-order
nature of the phase transition. Instead, we suggest that this degradation can
be explained by the residual stress retained in the elastomer when stretching
a polydomain LCE to prepare a monodomain LCE at the second stage of the widely
adopted two-step crosslinking procedure. 15, 1, 25

Figure 8: Effect of local chemical hinge on the critical behavior of SCLCE.
The critical nominal stress $\sigma_{\rm n,crit}^{*}$ and critical transition
temperature $T_{\rm crit}^{*}$ are plotted as a function of the local coupling
parameter $u_{C}$. The vertical dashed line denotes the discontinuity of
$T_{\rm crit}^{*}$ and the singularity of $\sigma_{\rm n,crit}^{*}$ at the
particular $u_{C}=-0.08u_{\rm qq}$. $m=1$, $r=1$, $u_{\rm pp}=0.07u_{\rm qq}$,
$u_{\rm pq}=0.2u_{\rm qq}$, and $\mu^{*}=0.5$.

To predictably control the ordering and deformation of LCEs at various working
conditions, it is important to understand the relationship between the
critical point and the network structures. Figure 8 summarizes the effect of
local chemical hinge on the critical applied nominal stress $\sigma_{\rm
n,crit}$ and the corresponding critical transition temperature $T_{\rm crit}$.
$T_{\rm crit}$ increases with $u_{C}$ but shows a discontinuity at the
critical hinge value, i.e. $u_{C}=-0.08u_{\rm qq}$. On the other hand,
$\sigma_{\rm n,crit}$ changes non-monotonically, and even exhibits a
singularity at the critical $u_{C}$. For tightly hinged side-on SCLCE (a very
positive $u_{C}$), only a small elongational stress is sufficient to bring the
system to the supercritical regime. Similar small compressive stress is needed
for the case of tightly hinged end-on SCLCE (a very negative $u_{C}$). As the
hinge is loosened (i.e., the magnitude of $u_{C}$ decreases), the magnitude of
the stress diverges towards the opposite limit near the singularity. At
$u_{C}=-0.08u_{\rm qq}$, the global coupling effect is exactly canceled by the
local coupling, such that mesogens are effectively decoupled from backbone
segments. Therefore, no matter how large the network is deformed by the strong
applied stress, it does not have any effect on the nematic ordering of
mesogens, and the first-order transition will never degrade into a continuous
one in this critical case.

###  3.4 Mechanical Response to Applied Stress

Understanding the stress-strain relationship plays is crucial for the
application of rubbery materials. In LCEs, the nematic ordering causes
deviation in the mechanical response from the classical rubber elasticity. The
problem becomes more complicated due to the interplay between mechanical
deformation and nematic phase transition: the response is different between
the nematic and isotropic phases.

Figure 9a plots the nominal stress as a function of deformation for side-on
SCLCEs ($u_{C}=0.25u_{\rm qq}$) with various LC grafting densities at the same
fixed temperature. The behavior of isotropic rubber $\sigma_{\rm
n}=\mu(\lambda-1/\lambda^{2})$ is also presented as a reference. The
mechanical responses can be classified into three different types depending on
the relation of the working temperature in comparison with $T_{\rm NI}$ and
$T_{\rm crit}$ of the LCE. For the first case as the sparsely grafted LCE
($m=3$), the working temperature chosen is higher than $T_{\rm crit}$ (red
point in Fig. 9b), the LCE is always in the supercritical state during the
deformation. The $\sigma_{\rm n}$-$\lambda$ curve is similar to that of the
classical rubber but with a lower modulus. The LCE is easier to be deformed
even in its paranematic state, because deformation is promoted by the
alignment of backbone segments. For the second case as the LCE with
intermediate grafting density ($m=2$), the working temperature is located
between $T_{\rm NI}$ and $T_{\rm crit}$ (blue point in Fig. 9b). In the regime
of small deformation, the LCE remains in its paranematic state and is softer
than the case of $m=3$. The stronger nematic ordering due to the increase of
LC loading enhances backbone alignment as a result of the coupling effects. As
$\sigma_{\rm n}$ becomes large enough, the stress triggers the discontinuous
transition to the nematic phase, leading to a plateau in the $\sigma_{\rm
n}$-$\lambda$ curve. After the transition, the deformation becomes easier as
the elastomer is stretched in the nematic state. For the third case as the
fully grafted LCE ($m=1$), the working temperature is lower than $T_{\rm NI}$
(green point in Fig. 9b). The elastomer has a spontaneous deformation prior to
the applied stress. The initial plateau at $\sigma_{\rm n}=0$ is followed by
the stretch of the nematic elastomer.

Figure 9: Mechanical response of SCLCEs ($u_{C}=0.25u_{\rm qq}$) with
different LC grafting density at the same fixed temperature. LC grafting
density is represented by $m$, the interval between two adjacent mesogens in
the same chain. (a) Rescaled nominal stress $\sigma_{\rm n}^{*}$ as a function
of network deformation $\lambda$ for different $m$ at the working temperature
$T^{*}=0.048$. (b) Rescaled nematic-isotropic transition temperature $T_{\rm
NI}^{*}$ and critical temperature $T_{\rm c}^{*}$ as a function of $m$. The
horizontal dashed line denotes the working temperature $T^{*}=0.048$. The red,
blue, and green points locate the case of $m=1$, $m=2$, and $m=3$,
respectively. $r=1$, $u_{\rm pp}=u_{\rm p^{*}p^{*}}=u_{\rm pp^{*}}=0.07u_{\rm
qq}$, $u_{\rm pq}=u_{\rm p^{*}q}=0.2u_{\rm qq}$, and $\mu^{*}=0.5$.

Our results elucidate that the pattern of the mechanical response can be well
controlled by the rational design of the LCE structure. In the current work,
we predict that the plateau behavior during stretching originates from the
discontinuous deformation accompanied by the first-order nematic phase
transition. We note that similar plateau signature in the stress-strain curve
of LCE-based materials can also be found in other cases. One example is the
soft elasticity with respect to the director rotation if the stress is applied
perpendicularly to the original nematic director.32 Another example is the re-
alignment of local directors when stretching a polydomain LCE into a
monodomain one.29, 20

##  4 Conclusions

In this work, we develop a self-consistent field theory which provides a
unified description of both main-chain and side-chain LCEs at the molecular
level. Molecular features, such as density of elastic strands, the strength
and architecture of chemical hinge, and LC grafting density are systematically
included. Our theory fully captures the coupling between nematic ordering,
backbone alignment and network deformation. The phase behavior and shape
deformation both without and with the applied external stress are described in
a single framework. We examine the effects of molecular structures on the
nematic ordering and spontaneous shape change under the force-free condition.
Compared to the uncrosslinked LCP, crosslinking is found to suppress the
nematic ordering by reducing both the nematic order parameter and transition
temperature as a result of the elastic energy stored during network
deformation. Our results suggest that less crosslinked LCE achieves higher
work capacity. While side-on SCLCEs always prefer elongational deformation,
our theory captures the experimental observation that the shape change of end-
on SCLCEs can be either elongation or contraction, determined by the
competition between the global coupling associated with the nematic field and
the local coupling due to the chemical hinge. We elucidate a scaling
relationship between the nematic-isotropic transition temperature and the
shear modulus as $T_{\rm NI,0}-T_{\rm
NI}({\mu^{*}})\sim{\mu^{*}}^{\frac{4}{5}}$, which is universal for both MCLCE
and SCLCEs with different hinge architectures. Furthermore, we show that the
pattern of shape change can be well controlled by adjusting the LC grafting
density. For end-on SCLCE with intermediate LC grafting density, the network
deformation exhibits dual modes: contraction at higher temperatures and
elongation at lower temperatures. The dependence of the transition temperature
on LC grafting density predicted by our theory is in quantitive agreement with
experimental data reported in previous literature.

The phase behavior and mechanical response of LCEs under applied stress depend
on the direction of the force and the nematic director. Nematic ordering is
enhanced for side-on SCLCEs under elongational stress. The first-order
transition eventually degrades to a continuous manner beyond the critical
stress. We suggest that the continuous transition for LCEs observed in
experiments can be explained by the residual stress retained in the elastomer
when they are prepared via a two-step crosslinking procedure. At a critical
strength of the local hinge, the local and global couplings cancel each other,
leading to a discontinuity in the critical temperature and a singularity in
the critical stress. Furthermore, the mechanical response of the LCEs deviates
significantly from the classical rubber elasticity. The stress-deformation
curve shows different patterns depending on the working temperature relative
to the transition temperature and critical temperature. A plateau appears
accompanied by the first-order nematic phase transition.

Although the current work focuses on a homogeneous network with uniaxial
deformation along the director, our theory is universal for LCEs with spatial
inhomogeneity and under any form of stress. Our theory can be easily
generalized to more complex systems, for example, heteropolymer networks
consisting of multi-types of mesogens and LCEs swollen in isotropic or nematic
solvents. Our theory can also incorporate the effect of other types of
external stimuli, such as light, electric and magnetic fields. Modeling the
phase and mechanical behaviors of LCEs at the molecular level is crucial for
the establishment of their structure-property relationship. Our work is thus
an important step towards the predictable design of LCE-based materials.

R. W. acknowledges the support from the University of California, Berkeley.
This research used the computational resources provided by the Kenneth S.
Pitzer Center for Theoretical Chemistry.

Detail derivation of self-consistent field theory for liquid crystal
elastomers; numerical details for solving the self-consistent equations. This
material is available free of charge via the Internet at http://pubs.acs.org/.

## References

  * Warner and Terentjev 2007 Warner, M.; Terentjev, E. M. _Liquid crystal elastomers_ ; Oxford University Press, 2007; Vol. 120
  * Stuart et al. 2010 Stuart, M. A.; Huck, W. T.; Genzer, J.; Müller, M.; Ober, C.; Stamm, M.; Sukhorukov, G. B.; Szleifer, I.; Tsukruk, V. V.; Urban, M.; Winnik, F.; Zauscher, S.; Luzinov, I.; Minko, S. Emerging applications of stimuli-responsive polymer materials. _Nat. Mater._ 2010, _9_ , 101–113
  * De Gennes et al. 1997 De Gennes, P.-G.; Hébert, M.; Kant, R. Artificial muscles based on nematic gels. _Macromol. Symp._ 1997, _113_ , 39–49
  * Li et al. 2006 Li, M. H.; Keller, P.; Antonietti, M.; Lacey, D.; Meyer, R. B. Artificial muscles based on liquid crystal elastomers. _Phil. Trans. R. Soc. A_ 2006, _364_ , 2763–2777
  * Tian et al. 2018 Tian, H.; Wang, Z.; Chen, Y.; Shao, J.; Gao, T.; Cai, S. Polydopamine-Coated Main-Chain Liquid Crystal Elastomer as Optically Driven Artificial Muscle. _ACS Appl. Mater. Interfaces_ 2018, _10_ , 8307–8316
  * Chen et al. 2020 Chen, Y.; Chen, C.; Rehman, H. U.; Zheng, X.; Li, H.; Liu, H.; Hedenqvist, M. S. Shape-memory polymeric artificial muscles: Mechanisms, applications and challenges. _Molecules_ 2020, _25_ , 4246
  * Ohm et al. 2010 Ohm, C.; Brehmer, M.; Zentel, R. Liquid crystalline elastomers as actuators and sensors. _Adv. Mater._ 2010, _22_ , 3366–3387
  * Kularatne et al. 2017 Kularatne, R. S.; Kim, H.; Boothby, J. M.; Ware, T. H. Liquid crystal elastomer actuators: Synthesis, alignment, and applications. _J. Polym. Sci., Part B: Polym. Phys._ 2017, _55_ , 395–411
  * Ko and Javey 2017 Ko, H.; Javey, A. Smart Actuators and Adhesives for Reconfigurable Matter. _Acc. Chem. Res._ 2017, _50_ , 691–702
  * Guin et al. 2018 Guin, T.; Settle, M. J.; Kowalski, B. A.; Auguste, A. D.; Beblo, R. V.; Reich, G. W.; White, T. J. Layered liquid crystal elastomer actuators. _Nat. Commun._ 2018, _9_ , 2531
  * Shafiq et al. 2020 Shafiq, Y.; Henricks, J.; Ambulo, C. P.; Ware, T. H.; Georgakopoulos, S. V. A Passive RFID Temperature Sensing Antenna with Liquid Crystal Elastomer Switching. _IEEE Access_ 2020, _8_ , 24443–24456
  * Mistry et al. 2020 Mistry, D.; Nikkhou, M.; Raistrick, T.; Hussain, M.; Jull, E. I.; Baker, D. L.; Gleeson, H. F. Isotropic Liquid Crystal Elastomers as Exceptional Photoelastic Strain Sensors. _Macromolecules_ 2020, _53_ , 3709–3718
  * Fleischmann et al. 2012 Fleischmann, E. K.; Liang, H. L.; Kapernaum, N.; Giesselmann, F.; Lagerwall, J.; Zentel, R. One-piece micropumps from liquid crystalline core-shell particles. _Nat. Commun._ 2012, _3_ , 1178
  * Guo et al. 2023 Guo, H.; Saed, M. O.; Terentjev, E. M. Mechanism of Pressure-Sensitive Adhesion in Nematic Elastomers. _Macromolecules_ 2023, _56_ , 6247–6255
  * Disch et al. 1994 Disch, S.; Schmidt, C.; Finkelmann, H. Nematic elastomers beyond the critical point. _Macromol. Rapid Commun._ 1994, _15_ , 303–310
  * Tajbakhsh and Terentjev 2001 Tajbakhsh, A. R.; Terentjev, E. M. Spontaneous thermal expansion of nematic elastomers. _Eur. Phys. J. E_ 2001, _6_ , 181–188
  * Cordoyiannis et al. 2007 Cordoyiannis, G.; Lebar, A.; Zalar, B.; Žumer, S.; Finkelmann, H.; Kutnjak, Z. Criticality controlled by cross-linking density in liquid single-crystal elastomers. _Phys. Rev. Lett._ 2007, _99_ , 197801
  * Cordoyiannis et al. 2009 Cordoyiannis, G.; Lebar, A.; Rožič, B.; Zalar, B.; Kutnjak, Z.; Žumer, S.; Brömmel, F.; Krause, S.; Finkelmann, H. Controlling the critical behavior of paranematic to nematic transition in main-chain liquid single-crystal elastomers. _Macromolecules_ 2009, _42_ , 2069–2073
  * Donato et al. 2023 Donato, S.; Martella, D.; de Luna, M. S.; Arecchi, G.; Querceto, S.; Ferrantini, C.; Sacconi, L.; Brient, P. L.; Chatard, C.; Graillot, A.; Wiersma, D. S.; Parmeggiani, C. The Role of Crosslinker Molecular Structure on Mechanical and Light-Actuation Properties in Liquid Crystalline Networks. _Macromol. Rapid Commun._ 2023, _44_ , 2200958
  * Barnes et al. 2022 Barnes, M.; Cetinkaya, S.; Ajnsztajn, A.; Verduzco, R. Understanding the effect of liquid crystal content on the phase behavior and mechanical properties of liquid crystal elastomers. _Soft Matter_ 2022, _18_ , 5074–5081
  * Saed et al. 2017 Saed, M. O.; Torbati, A. H.; Starr, C. A.; Visvanathan, R.; Clark, N. A.; Yakacki, C. M. Thiol-acrylate main-chain liquid-crystalline elastomers with tunable thermomechanical properties and actuation strain. _J. Polym. Sci., Part B: Polym. Phys._ 2017, _55_ , 157–168
  * Saed et al. 2017 Saed, M. O.; Volpe, R. H.; Traugutt, N. A.; Visvanathan, R.; Clark, N. A.; Yakacki, C. M. High strain actuation liquid crystal elastomers via modulation of mesophase structure. _Soft Matter_ 2017, _13_ , 7537–7547
  * Saed et al. 2019 Saed, M. O.; Ambulo, C. P.; Kim, H.; De, R.; Raval, V.; Searles, K.; Siddiqui, D. A.; Cue, J. M. O.; Stefan, M. C.; Shankar, M. R.; Ware, T. H. Molecularly-Engineered, 4D-Printed Liquid Crystal Elastomer Actuators. _Adv. Funct. Mater._ 2019, _29_ , 1806412
  * Rogez et al. 2018 Rogez, D.; Krause, S.; Martinoty, P. Main-chain liquid-crystal elastomers: Versus side-chain liquid-crystal elastomers: Similarities and differences in their mechanical properties. _Soft Matter_ 2018, _14_ , 6449–6462
  * Herbert et al. 2022 Herbert, K. M.; Fowler, H. E.; McCracken, J. M.; Schlafmann, K. R.; Koch, J. A.; White, T. J. Synthesis and alignment of liquid crystalline elastomers. _Nat. Rev. Mater._ 2022, _7_ , 23–38
  * Greve and Finkelmann 2001 Greve, A.; Finkelmann, H. Nematic elastomers: The dependence of phase transformation and orientation processes on crosslinking topology. _Macromol. Chem. Phys._ 2001, _202_ , 2926–2946
  * Davis 1993 Davis, F. J. Liquid-crystalline Elastomers. _J. Mater. Chem._ 1993, _3_ , 551–562
  * Xu et al. 2021 Xu, Y.; Dupont, R. L.; Yao, Y.; Zhang, M.; Fang, J. C.; Wang, X. Random Liquid Crystalline Copolymers Consisting of Prolate and Oblate Liquid Crystal Monomers. _Macromolecules_ 2021, _54_ , 5376–5387
  * Schätzle et al. 1989 Schätzle, J.; Kaufhold, W.; Finkelmam, H. Nematic elastomers: The influence of external mechanical stress on the liquid-crystalline phase behavior. _Macromol. Chem. Phys._ 1989, _190_ , 3269–3284
  * Bladon et al. 1993 Bladon, P.; Terentjev, E. M.; Warner, M. Transitions and instabilities in liquid crystal elastomers. _Phys. Rev. E_ 1993, _47_ , R3838–R3840
  * Mihai and Goriely 2021 Mihai, L. A.; Goriely, A. Instabilities in liquid crystal elastomers. _MRS Bull._ 2021, _46_ , 784–794
  * Warner et al. 1994 Warner, M.; Bladon, P.; Terentjev, E. “Soft elasticity”—deformation without resistance in liquid crystal elastomers. _J. Phys. II_ 1994, _4_ , 93–102
  * Ware et al. 2016 Ware, T. H.; Biggins, J. S.; Shick, A. F.; Warner, M.; White, T. J. Localized soft elasticity in liquid crystal elastomers. _Nat. Commun._ 2016, _7_ , 10781
  * Warner and Terentjev 1996 Warner, M.; Terentjev, E. M. Nematic elastomers—a new state of matter? _Prog. Polym. Sci._ 1996, _21_ , 853–891
  * Finkelmann et al. 2001 Finkelmann, H.; Greve, A.; Warner, M. The elastic anisotropy of nematic elastomers. _Eur. J. Phys. E_ 2001, _5_ , 281–293
  * Warner and Wang 1991 Warner, M.; Wang, X. J. Elasticity and phase behavior of nematic elastomers. _Macromolecules_ 1991, _24_ , 4932–4941
  * Warner et al. 1988 Warner, M.; Gelling, K. P.; Vilgis, T. A. Theory of nematic networks. _J. Chem. Phys._ 1988, _88_ , 4008–4013
  * Zhang et al. 2019 Zhang, Y.; Xuan, C.; Jiang, Y.; Huo, Y. Continuum mechanical modeling of liquid crystal elastomers as dissipative ordered solids. _J. Mech. Phys. Solids_ 2019, _126_ , 285–303
  * Zhao and Meng 2023 Zhao, J.; Meng, F. Modeling Viscoelasticity and Dynamic Nematic Order of Exchangeable Liquid Crystal Elastomers. _Phys. Rev. Lett._ 2023, _131_ , 068101
  * De Gennes and Prost 1993 De Gennes, P.-G.; Prost, J. _The physics of liquid crystals_ ; Oxford University Press, 1993
  * Liarte et al. 2011 Liarte, D. B.; Salinas, S. R.; Yokoi, C. S. Elastic Maier-Saupe-Zwanzig model and some properties of nematic elastomers. _Phys. Rev. E_ 2011, _84_ , 011124
  * Warsono et al. 2015 Warsono; Yusuf, Y.; Nurwantoro, P.; Abraha, K. Theoretical Studies of the Effects of Magnetic Field on the Phase Transition of Swollen Liquid Crystal Elastomers. _Adv. Mater. Res._ 2015, _1123_ , 46–54
  * Pasini et al. 2005 Pasini, P.; Skačej, G.; Zannoni, C. A microscopic lattice model for liquid crystal elastomers. _Chem. Phys. Lett._ 2005, _413_ , 463–467
  * Corbett and Warner 2006 Corbett, D.; Warner, M. Nonlinear photoresponse of disordered elastomers. _Phys. Rev. Lett._ 2006, _96_ , 237802
  * Corbett and Warner 2008 Corbett, D.; Warner, M. Polarization dependence of optically driven polydomain elastomer mechanics. _Phys. Rev. E_ 2008, _78_ , 061701
  * Bai and Bhattacharya 2020 Bai, R.; Bhattacharya, K. Photomechanical coupling in photoactive nematic elastomers. _J. Mech. Phys. Solids_ 2020, _144_ , 104115
  * Matsuyama and Kato 2001 Matsuyama, A.; Kato, T. Volume phase transitions of nematic gels under an external field. _J. Chem. Phys._ 2001, _114_ , 3817–3822
  * Urayama et al. 2005 Urayama, K.; Arai, Y. O.; Takigawa, T. Anisotropic swelling and phase behavior of monodomain nematic networks in nematogenic solvents. _Macromolecules_ 2005, _38_ , 5721–5728
  * Cheewaruangroj and Terentjev 2015 Cheewaruangroj, N.; Terentjev, E. M. Shape instability on swelling of a stretched nematic elastomer filament. _Phys. Rev. E_ 2015, _92_ , 042502
  * Wang et al. 2016 Wang, X.; Miller, D. S.; Bukusoglu, E.; Pablo, J. J. D.; Abbott, N. L. Topological defects in liquid crystals as templates for molecular self-assembly. _Nat. Mater._ 2016, _15_ , 106–112
  * Wang and Wang 2010 Wang, R.; Wang, Z.-G. Theory of side-Chain liquid crystal polymers: Bulk behavior and chain conformation. _Macromolecules_ 2010, _43_ , 10096–10106
  * Zhuang and Wang 2012 Zhuang, B.; Wang, Z.-G. Anomalous concentration effects on phase behavior and nematic order in mixtures of side-chain liquid crystal polymers and low-molecular-weight liquid crystals. _Macromolecules_ 2012, _45_ , 6220–6229
  * Maier and Saupe 1959 Maier, W.; Saupe, A. Eine einfache molekular-statistische Theorie der nematischen kristallinflüssigen Phase. Teil l1. _Z. Naturforsch._ 1959, _14_ , 882–889
  * Fredrickson 2006 Fredrickson, G. _The equilibrium theory of inhomogeneous polymers_ ; Oxford University Press, 2006
  * Rubinstein and Colby 2003 Rubinstein, M.; Colby, R. H. _Polymer Physics_ ; Oxford University Press, 2003
  * Wang et al. 2016 Wang, R.; Alexander-Katz, A.; Johnson, J. A.; Olsen, B. D. Universal Cyclic Topology in Polymer Networks. _Phys. Rev. Lett._ 2016, _116_ , 188302
  * Zhong et al. 2016 Zhong, M.; Wang, R.; Kawamoto, K.; Olsen, B. D.; Johnson, J. A. Quantifying the impact of molecular defects on polymer network elasticity. _Science_ 2016, _353_ , 1264–1268
  * Urayama et al. 2005 Urayama, K.; Arai, Y. O.; Takigawa, T. Volume phase transition of monodomain nematic polymer networks in isotropic solvents accompanied by anisotropic shape variation. _Macromolecules_ 2005, _38_ , 3469–3474
  * Lebar et al. 2005 Lebar, A.; Kutnjak, Z.; Žumer, S.; Finkelmann, H.; Sánchez-Ferrer, A.; Zalar, B. Evidence of supercritical behavior in liquid single crystal elastomers. _Phys. Rev. Lett._ 2005, _94_ , 197801