2021

[2]Epifanio G. Virga These authors contributed equally to this work.

1]Laboratory for Computation and Visualization in Mathematics and Mechanics,
Institute of Mathematics, École Polytechnique Fédérale de Lausanne, Lausanne,
CH-1015, Switzerland

[2]Dipartimento di Matematica, Università di Pavia, via Ferrata 5, Pavia,
I-27100, Italy

# Bending and Stretching in a Narrow Ribbon of Nematic Polymer Networks

Harmeet Singh<EMAIL_ADDRESS><EMAIL_ADDRESS>[  *

###### Abstract

We study the spontaneous out-of-plane bending of a planar untwisted ribbon
composed of nematic polymer networks activated by a change in temperature. Our
theory accounts for both stretching and bending energies, which compete to
establish equilibrium. We show that when equilibrium is attained these energy
components obey a _complementarity_ relation: one is maximum where the other
is minimum. Moreover, we identify a _bleaching_ regime: for sufficiently large
values of an activation parameter (which measures the mismatch between the
degrees of order in polymer organization in the reference and current
configurations), the ribbon’s deformation is essentially independent of its
thickness.

###### keywords:

Nematic polymer networks, Nematic glasses, Soft matter elasticity, Nematic
elastomers, Ribbon theory, Thermally activated elastic materials

###### pacs:

[

MSC Classification]Primary 74, Secondary 74B20, 74K10, 74K35, 76A15

##  1 Introduction

A nematic elastomer differs from a standard elastomer in that it hosts
_nematogenic_ rod-like molecules within its polymer strands. If the
temperature is sufficiently low, the nematogenic components develop a tendency
to be aligned alike (the very signature of a nematic liquid crystal), thus
introducing a degree of anisotropy in the spatial organization of polymer
strands. This suffices to alter the elastic response of these materials. Since
the degree of nematic order depends on temperature, a change in the latter can
induce a spontaneous deformation of the body.

The nematic order is described by a scalar, here denoted $Q$, and a unit
vector $\bm{n}$ (or _director_), where the former has a bearing on the spatial
organization of polymer strands and the latter represents the average
orientation of the nematogenic constituents.

Nematic elastomers can be envisioned as a fluid-solid mixture, where the
nematic component is fluid and the cross-linked polymer matrix is solid
corbett:photomechanics . Generally, fluid and solid components enjoy a certain
degree of mutual independence. A full spectrum of materials exist, whose
mechanical behavior can be classified in accordance with the degree of freedom
allowed to $\bm{n}$ within the polymer matrix white:photomechanical_collection
; ware:programmed .

_Nematic polymer networks_ (NPNs) are nematic elastomers of a special kind:
their cross-linking is so tight that the nematic director is completely
_enslaved_ to the matrix deformation. Only these materials will form the
object of our study.

We shall assume that a NPN is cross-linked in the nematic phase, which will
constitute the stress-free, reference configuration of the bodies considered
here. Deformations will be induced by a change in temperature, which drives
the reference configuration out of equilibrium.111A change in temperature
affects the degree of orientational order among nematogenic moulecules; its
effect is similar to that driven at constant temperature by the shape change
induced by light on photoresponsive nematic molecules (not involved here).
See, for example, sonnet:model  for a discussion on photoactivable NPNs. A
scalar activation parameter will describe the mismatch between nematic order
in the reference and current configurations; it will act as a control
parameter.

We are interested in the slimiest of bodies, a _ribbon_. Our continuum theory
is based on the three-dimensional theory of nematic elastomers revolving
around the celebrated _trace formula_ for the volume elastic energy density,
which has long been studied blandon:deformation ; warner:theory ;
warner:elasticity ; warner:nematic  (as also illustrated in Sect. 6.1 of
warner:liquid  under the suggestive heading of _neo-classical_ theory).

For the readers interested in broadening their background on nematic
elastomers who are intimidated by the vast available literature, a recommended
reference is the influential book by Warner and Terentjev warner:liquid ; the
theoretical literature that preceded and prepared for it blandon:deformation ;
warner:soft ; terentjev:orientation ; verwey:soft ; verwey:multistage ;
verwey:elastic  is also of interest. General continuum theories have also been
proposed anderson:continuum ; zhang:continuum ; mihai:nematic , some also very
recently. Applications are plenty; a collection can be found in a book
white:photomechanical_book  and a recent special issue korley:introduction .
Finally, some valuable guidance can be gained from the reviews
mahimwalla:azobenzene ; ube:photomobile ; white:photomechanical ; ula:liquid ;
pang:photodeformable ; kuenstler:light ; warner:topographic .

This paper is an outgrowth of a previous study singh:model , where we derived
our theory for NPN ribbons from a reduction procedure from three to two space
dimensions proposed in ozenda:blend  with the intent of enclosing both
stretching and bending contents in the surface elastic energy of thin NPN
sheets. To make our development self-contained, we recall in Sect. 2 the
fundamentals of our theory; the articulation in subsections will aid the
experienced reader to gather the essential information. In Sect. 3, we specify
the kinematics of the simple problem that we consider here. Section 4, which
is the heart of the paper, is where we state and solve the variational problem
for the optimal shape of a thermally activated ribbon. A plethora of activated
equilibria are described in Sect. 5 for a whole range of activation and
geometric parameters. Finally, in Sect. 6, we collect two qualitative
conclusions of our study: one identifies a regular pattern in the relative
distribution of stretching and bending energies at equilibrium, the other
suggests the existence of an activation regime where the ribbon’s deformation
is insensitive to its thickness.

##  2 Preliminaries

To make our paper self-contained, we summarize in this section some results
already obtained in ozenda:blend  and singh:model . The reader familiar with
this material may jump ahead to the following sections.

Our theory is based on the _trace formula_ for the elastic volume energy
density (in the reference configuration) of nematic elastomers in the form put
forward by finkelmann:elastic ,

$f_{\mathrm{e}}=\frac{1}{2}k\bigg{\\{}\operatorname{tr}(\mathbf{F}^{\mathsf{T}}\mathbf{L}^{-1}\mathbf{F}\mathbf{L}_{0})+\ln\left(\frac{\det\mathbf{L}}{\det\mathbf{L}_{0}}\right)\bigg{\\}},$
(1)

where $\mathbf{F}$ is the deformation gradient and the _shear_ elastic modulus
$k$ is given by

$k=n_{\mathrm{s}}k_{B}T,$ (2)

in terms of the number density of polymer strands $n_{\mathrm{s}}$, the
Boltzmann constant $k_{B}$, and the absolute current temperature $T$. In (1),
$\mathbf{L}_{0}$ and $\mathbf{L}$ are the _step-length_ tensors of the
material in the reference and current configurations, respectively. They are
generally different because different are the arrangements of polymer chains
in the two configurations, if a change in temperature has occurred, altering
the orientational order of the nematogenic molecules that comprise the polymer
strands. Following verwey:elastic ; nguyen:theory , we represent them as

${\bf
L}_{0}=A_{0}(\mathbf{I}+S_{0}{\bm{n}}_{0}\otimes{\bm{n}}_{0})\quad\text{and}\quad{\bf
L}=A(\mathbf{I}+S{\bm{n}}\otimes{\bm{n}}),$ (3)

where $\mathbf{I}$ is the identity (in three-dimensional space) and
${\bm{n}}_{0}$ and $\bm{n}$ are the nematic directors in reference and current
configurations. Moreover,

$A_{0}=\ell_{\mathrm{m}}(1-Q_{0}),\quad S_{0}=\frac{3Q_{0}}{1-Q_{0}}$ (4a) and
their twins, $A=\ell_{\mathrm{m}}(1-Q),\quad S=\frac{3Q}{1-Q},$ (4b)

are expressions derived within the classical statistical mechanics model for a
polymer chain thought of as consisting in freely jointed nematogenic monomers,
each of length $\ell_{\mathrm{m}}$ (see Chapt. 6 of warner:liquid  and also
corbett:nonlinear ; corbett:polarization ). In (4), $Q_{0}$ (and $Q$)
represents the Maier-Saupe nematic scalar order parameter describing the
degree of orientation of monomers in the reference (and current) configuration
of the material.222This theory is an extension of the classical Gaussian
rubber elasticity. An exposition of statistical theories for rubber can be
found in the reference book treloar:non-linear_third . Adaptation of the
simplest realization of these theories deam:theory  to the case where the
distribution of monomers in a polymer chain is anisotropic delivers (1).

$Q$ (and, correspondingly, $Q_{0}$) is defined as

$Q:={\left\langle P_{2}(\bm{n}\cdot{\bm{u}})\right\rangle},$ (5)

where $P_{2}$ is the second Legendre polynomial, ${\bm{u}}\in\mathbb{S}^{2}$
is the unit vector along a single monomer, and the brackets
${\left\langle\cdots\right\rangle}$ denote ensemble average. $Q$ ranges in the
interval $[-\frac{1}{2},1]$, whose end-values represent the limiting cases of
${\bm{u}}$’s distributed isotropically in the plane orthogonal to $\bm{n}$ and
${\bm{u}}$’s aligned with $\bm{n}$, respectively.

Here, we shall assume that $Q$ (and $Q_{0}$) ranges in $(0,1)$. In this
interval, by (4b), $S$ is a monotonically increasing function of $Q$, which
also increases when the temperature $T$ is decreased below the isotropic-to-
nematic transition. Thus, a decrease in temperature results in an increase in
$S$. This explains why the mechanism of thermal activation can be effectively
described as having $S\neq S_{0}$. A temperature increase makes $S<S_{0}$,
whereas a temperature decrease makes $S>S_{0}$.

When, as in our case, the scalar order parameters $S_{0}$ and $S$ are
prescribed, the second term in (1) is not affected by the deformation and can
be omitted, thus reducing (1) to the following _bare_ trace formula of
blandon:deformation  (also discussed in warner:new ), which depends only on
$\mathbf{F}$ and $\bm{n}$ and will be adopted from now on,

$f_{\mathrm{e}}=\frac{1}{2}k\operatorname{tr}({\bf F}^{\mathsf{T}}{\bf
L}^{-1}{\bf F}{\bf L}_{0}).$ (6)

The degree of cross-linking in the material is responsible for the mobility of
the nematic director relative to the network matrix ware:programmed . A
_nematic polymer network_ (NPN) represents the tightest end of the cross-
linking spectrum, where $\bm{n}$ is _enslaved_ to the deformation.333The name
_nematic polymer network_ was proposed in white:programmable . Other authors
prefer to call these materials _nematic glasses_ (see, for example,
modes:disclination ). For these materials, which are the only ones considered
here, we assume that

$\bm{n}=\frac{\mathbf{F}\bm{n}_{0}}{|\mathbf{F}\bm{n}_{0}|},$ (7)

which says that the nematic director $\bm{n}_{0}$ imprinted in the polymer
network at the cross-linking time is _conveyed_ by the solid matrix of the
body. We further assume that the material is _incompressible_ ,444Even if, as
discussed in white:photomechanical_collection , not all NPNs are strictly so.
so that $\mathbf{F}$ is subject to the constraint

$\det\mathbf{F}=1.$ (8)

By use of (3), (2), (4b), and (7), (6) can be given the following form
ozenda:blend ,

$f_{\mathrm{e}}=\frac{1}{2}e_{0}E(\mathbf{C}_{\bm{f}}),$ (9)

where

$e_{0}:=k\frac{A_{0}}{A}=n_{\mathrm{s}}k_{B}T\frac{3+S}{3+S_{0}},$ (10)

$\mathbf{C}_{\bm{f}}:={\bf F}^{\mathsf{T}}{\bf F}$ is the right Cauchy-Green
tensor associated with the deformation $\bm{f}$ and

$E(\mathbf{C}_{\bm{f}}):=\operatorname{tr}\mathbf{C}_{\bm{f}}+\frac{S_{0}}{S+1}\bm{n}_{0}\cdot\mathbf{C}_{\bm{f}}\bm{n}_{0}-\frac{S}{S+1}\frac{\bm{n}_{0}\cdot\mathbf{C}_{\bm{f}}^{2}\bm{n}_{0}}{\bm{n}_{0}\cdot\mathbf{C}_{\bm{f}}\bm{n}_{0}}.$
(11)

To apply our theory to the spontaneous deformation of a ribbon out of its
plane (which is the theme of the following sections), we must first learn how
to reduce the volume free-energy density in (9) to a surface free-energy
density to be attributed to a thin sheet, thought of as a two-dimensional
body. This task was accomplished in ozenda:blend ; in the following
subsection, we recall the main results of this study.

###  2.1 Dimension reduction

The method adopted in ozenda:blend  is classical and consists in performing an
expansion of $f_{\mathrm{e}}$ retaining up to the cubic terms in the sheet’s
thickness, thus identifying _stretching_ and _bending_ contents in the surface
energy by the power in the thickness they scale with.

We identify the undeformed body with a slab $\mathsf{S}$ of thickness $2h$ and
midsurface $\mathscr{S}_{0}$ in the $\\{\bm{e}_{3},\bm{e}_{1}\\}$ plane of a
Cartesian frame. We further assume that $\bm{n}_{0}$ is imprinted in
$\mathsf{S}$ so that it does not depend on the $x_{2}$ (out-of-plane)
coordinate and $\bm{n}_{0}\cdot\bm{e}_{2}\equiv 0$. Moreover, we represent the
three-dimensional deformation $\bm{f}$ as

$\bm{f}(\bm{x},x_{2})=\bm{y}(\bm{x})+\phi(\bm{x},x_{2})\bm{\nu},$ (12)

where $\bm{x}$ varies in $\mathscr{S}_{0}$, $x_{2}$ ranges in the interval
$[-h,h]$, $\bm{\nu}$ is the normal to the midsurface
$\mathscr{S}=\bm{y}(\mathscr{S}_{0})$ of the deformed slab
$\bm{f}(\mathsf{S})$ (see Fig. 1),

Figure 1: Schematic representation of the deformation of a thin sheet.
$\mathscr{S}_{0}$ is the planar midsurface of the slab $\mathsf{S}$ of
thickness $2h$ representing the reference configuration of the body; it lies
in the plane $\\{\bm{e}_{3},\bm{e}_{1}\\}$ of Cartesian frame
$\\{\bm{e}_{1},\bm{e}_{2},\bm{e}_{3}\\}$. $\mathscr{S}$ is the image under the
mapping $\bm{y}$ of $\mathscr{S}_{0}$; it is an oriented surface with unit
normal $\bm{\nu}$.

and $\phi(\bm{x},x_{2})$ is a function to be determined,555In the classical
theory of plates, the _Kirchhoff-Love_ hypothesis stipulates that $\phi\equiv
x_{2}$ (see, for example, podio:exact , for a modern treatment). In
ozenda:kirchhoff , the Kirchhoff-Love hypothesis was reformulated in the more
general form adopted here and criteria were suggested to identify the function
$\phi$, none of which delivered exactly the original Kirchhoff-Love form.
enjoying the property

$\phi(\bm{x},0)=0.$ (13)

As shown in ozenda:blend , the constraint of incompressibility (8) determines
$\phi$ in the form

$\phi(\bm{x},x_{2})=x_{2}-Hx_{2}^{2}+\frac{1}{3}(6H^{2}-K)x_{2}^{3}+O(x_{2}^{4}),$
(14)

where $H$ and $K$ are the _mean_ and _Gaussian_ curvature of $\mathscr{S}$,
respectively, defined as

$H(\bm{y}(\bm{x})):=\frac{1}{2}\operatorname{tr}(\nabla\\!_{\mathrm{s}}\bm{\nu})\quad\text{and}\quad
K(\bm{y}(\bm{x})):=\det(\nabla\\!_{\mathrm{s}}\bm{\nu}),$ (15)

in terms of the two-dimensional curvature tensor
$\nabla\\!_{\mathrm{s}}\bm{\nu}$ at the point $\bm{y}(\bm{x})$ on
$\mathscr{S}$.

Moreover, under assumption (7) and for $h$ sufficiently small, the nematic
director $\bm{n}$ in the present configuration $\bm{f}(\mathsf{S})$ is
delivered by (see sonnet:model )

$\bm{n}(\bm{f}(\bm{x},x_{2}))=\frac{(\mathbf{I}+\phi\nabla\\!_{\mathrm{s}}\bm{\nu})(\nabla\bm{y})\bm{n}_{0}}{\lvert(\mathbf{I}+\phi\nabla\\!_{\mathrm{s}}\bm{\nu})(\nabla\bm{y})\bm{n}_{0}\rvert}.$
(16)

Since $\bm{n}_{0}\cdot\bm{e}_{2}\equiv 0$,
$(\nabla\bm{y})\bm{n}_{0}\cdot\bm{\nu}\equiv 0$, and
$\nabla\\!_{\mathrm{s}}\bm{\nu}$ is a symmetric tensor mapping the local
tangent plane to $\mathscr{S}$ into itself, (16) shows that
$\bm{n}\cdot\bm{\nu}=0$ everywhere within $\bm{f}(\mathsf{S})$, but while
$\bm{n}_{0}$ is _uniform_ across the thickness of $\mathsf{S}$, $\bm{n}$ is
_not_ so across the thickness of $\bm{f}(\mathsf{S})$. However, by (16) and
(13), on $\mathscr{S}$ $\bm{n}$ is given by

$\bm{n}(\bm{f}(\bm{x},0))=\frac{(\nabla\bm{y})\bm{n}_{0}}{\lvert(\nabla\bm{y})\bm{n}_{0}\rvert},$
(17)

and so it appears to be conveyed by the deformation of $\mathscr{S}_{0}$, in
accordance with the three-dimensional constraint (7). Accordingly, to mimic
(8), we shall assume that

$\det(\nabla\bm{y})=1,$ (18)

which makes $\mathscr{S}_{0}$ an _inextensible_ surface.

Laborious calculations, building upon (14), (17), and (18), established in
ozenda:blend  that the surface elastic energy density for the deformation of
$\mathscr{S}_{0}$ into $\mathscr{S}$ is given by

$\widetilde{f}_{\mathrm{e}}=\frac{1}{2}e_{0}h(f_{\mathrm{s}}+f_{\mathrm{b}})+O(h^{5}),$
(19)

where

$\displaystyle f_{\mathrm{s}}$ $\displaystyle=$
$\displaystyle\frac{2}{{\color[rgb]{0,0,0}S}+1}\left(\operatorname{tr}\mathbf{C}+S_{0}\bm{n}_{0}\cdot\mathbf{C}\bm{n}_{0}+\frac{S}{\bm{n}_{0}\cdot\mathbf{C}\bm{n}_{0}}\right),$
(20a) $\displaystyle f_{\mathrm{b}}$ $\displaystyle=$
$\displaystyle\frac{2h^{2}}{3}\bigg{\\{}2(8H^{2}-K)$ (20b) $\displaystyle+$
$\displaystyle\frac{1}{S+1}\left[\left(\frac{3S}{a_{0}^{2}}-a_{0}^{2}S_{0}-\operatorname{tr}\mathbf{C}\right)K-\frac{4S}{a_{0}^{2}}(2H-\kappa_{n})\kappa_{n}\right]\bigg{\\}}$

are the (dimensionless) stretching and bending free-energy _contents_ ,
respectively. In (20), $\mathbf{C}:=(\nabla\bm{y})^{\mathsf{T}}(\nabla\bm{y})$
is the right Cauchy-Green tensor associated with the deformation $\bm{y}$,
$a_{0}^{2}:=\bm{n}_{0}\cdot\mathbf{C}\bm{n}_{0}$, and

$\kappa_{n}:=\bm{n}\cdot(\nabla\\!_{\mathrm{s}}\bm{\nu})\bm{n},$ (21)

which couples $\bm{n}$ with the geometry of $\mathscr{S}$.

As shown in detail in Appendix A, in the absence of activation, that is, for
$S=S_{0}$, $f_{\mathrm{s}}$ attains its minimum subject to (18) for
$\mathbf{C}=\mathbf{I}_{2}$, where $\mathbf{I}_{2}$ is the identity in two-
dimensional space. Thus, in the membrane approximation, the surface elastic
energy of an _inactive_ sheet is minimized by isometries.

Finally, dropping in (19) the higher order terms in $h$, we write the total
elastic free-energy functional in the form

$\mathscr{F}[\bm{y}]:=\frac{1}{2}e_{0}h\int_{\mathscr{S}_{0}}(f_{\mathrm{s}}+f_{\mathrm{b}})dA,$
(22)

where now $A$ denotes the area measure and $e_{0}h$ is a characteristic
surface energy density precisely as $e_{0}$ in (10) is a characteristic volume
energy density.

###  2.2 General theory of NPN ribbons

We reproduce here only the essential aspects of the theory of NPN ribbons,
obtained through a dimension reduction procedure developed in singh:model .
The procedure reduces the surface energy density of a sheet of NPNs, derived
in ozenda:blend  and summarized in Sect. 2, to a line energy density valid for
ribbons, whose reference and deformed configurations are assumed to be ruled
surfaces.

Consider a ribbon made of NPNs, whose stress-free natural configuration
$\mathscr{S}_{0}$, and deformed configuration $\mathscr{S}$ are depicted in
Fig. 2.

Figure 2: (a) The reference configuration $\mathscr{S}_{0}$ of a planar
ribbon, ruled by the directrix ${\bm{r}}_{0}(s)$ and the generators in the
direction of the unit vector ${\bm{q}}_{0}(s)$. An orthonormal set of
directors $\\{{\bm{e}}_{1}(s),{\bm{e}}_{2},{\bm{e}}_{3}(s)\\}$ are adapted to
the centreline such that ${\bm{e}}_{3}(s)$ is tangential to ${\bm{r}}_{0}(s)$,
and ${\bm{e}}_{2}$ is pointing out of the plane of the paper. The angle
between ${\bm{q}}_{0}$ and ${\bm{e}}_{1}$ is denoted by $\alpha_{0}$. (b) A
current, non-planar, configuration $\mathscr{S}$ of the ribbon with the
directrix denoted by ${\bm{r}}(s)$, and the direction of the generators given
by ${\bm{q}}(s)$. An orthonormal triad
$\\{{\bm{d}}_{1}(s),{\bm{d}}_{2}(s),{\bm{d}}_{3}(s)\\}$ is attached to the
centreline with ${\bm{d}}_{3}(s)$ oriented along the tangent vector, and
${\bm{d}}_{1}(s)$ in the tangent plane at $\mathscr{S}$. The angle between
${\bm{d}}_{1}(s)$ and ${\bm{q}}_{0}(s)$ is shown by $\alpha$.

We should think of them as follows. $\mathscr{S}_{0}$ is a stress-free
configuration of the inactive ribbon (when $S=S_{0}$). Once activated, so that
$S\neq S_{0}$, $\mathscr{S}_{0}$ is no longer stress-free and the ribbon falls
out of equilibrium; a deformation then restores equilibrium by morphing
$\mathscr{S}_{0}$ into $\mathscr{S}$.

We choose a material line ${\bm{r}}_{0}\in\mathscr{S}_{0}$ in the natural
configuration, endowed with a triad of orthonormal vectors
$\\{{\bm{e}}_{1},{\bm{e}}_{2},{\bm{e}}_{3}\\}$, oriented as shown in Fig. 2.
Here $\bm{e}_{2}$ is independent of $s$, while $\\{\bm{e}_{3},\bm{e}_{1}\\}$
is a movable frame in the plane orthogonal to $\bm{e}_{2}$. Letting
$\bm{\omega}$ denote the Darboux vector associated with the frame
$\\{{\bm{e}}_{1},{\bm{e}}_{2},{\bm{e}}_{3}\\}$, so that
$\partial_{s}\bm{e}_{i}=\bm{\omega}\times\bm{e}_{i}$, we easily see that
$\bm{\omega}=\omega_{2}\bm{e}_{2}$. We refer to ${\bm{r}}_{0}$ as the
directrix, or the centerline, of the ribbon. We parameterize the natural
configuration as a ruled surface with the following representation,

${\bm{x}}(s,t)={\bm{r}}_{0}(s)+t{\bm{q}}_{0}(s)\,,$ (23)

where ${\bm{q}}_{0}(s)$ is a unit vector given by
${\bm{q}}_{0}(s)=\cos\alpha_{0}(s){\bm{e}}_{1}+\sin\alpha_{0}(s){\bm{e}}_{3}$
with $-\frac{\pi}{2}\leq\alpha_{0}\leq\frac{\pi}{2}$. Here $s$ and $t$ are
material coordinates along the centerline ${\bm{r}}_{0}$ and ${\bm{q}}_{0}$,
respectively, and $\alpha_{0}$ is a smooth function of $s$ such that
${\bm{q}}_{0}(0)$ and ${\bm{q}}_{0}(L)$ coincide, respectively, with the short
edges of the ribbon at $s=0$ and $s=L$. Formally, in the $(s,t)$ coordinates,
the ribbon is represented by the set $\\{(s,t):0\leqq s\leqq L,t_{-}(s)\leqq
t\leqq t_{+}(s)\\}$, where $t_{\pm}(s)$ are functions describing the two long
sides of the ribbon. We will identify the unit vector ${\bm{q}}_{0}$ with the
imprinted nematic director ${\bm{n}}_{0}$, thereby imparting material
character to the former,

${\bm{q}}_{0}={\bm{n}}_{0}\,.$ (24)

When the ribbon deforms, the centerline ${\bm{r}}_{0}(s)$ and the unit vectors
${\bm{q}}_{0}(s)$ are convected to ${\bm{r}}(s)\in\mathscr{S}$ and
${\bm{q}}(s)\in\mathscr{S}$, where $\mathscr{S}$ denotes a deformed
configuration of the ribbon. Similar to ${\bm{r}}_{0}(s)$, the deformed
centerline ${\bm{r}}(s)$ is endowed with an ordered orthonormal director frame
$\\{{\bm{d}}_{1}(s),{\bm{d}}_{2}(s),{\bm{d}}_{3}(s)\\}$, where
${\bm{d}}_{3}(s)$ is constrained to lie along the tangent (see Fig. 2). The
kinematics of the centerline are captured by the following two relations,

$\partial_{s}{\bm{r}}=v_{3}{\bm{d}}_{3}\quad\text{and}\quad\partial_{s}{\bm{d}}_{i}={\bm{u}}\times{\bm{d}}_{i}\,,$
(25)

where $v_{3}(s)$ is identified as the stretch of the centerline, and
${\bm{u}}(s)$ is the Darboux vector associated with the frame
$\\{{\bm{d}}_{1},{\bm{d}}_{2},{\bm{d}}_{3}\\}$. The director components in the
expansion ${\bm{u}}=u_{1}{\bm{d}}_{1}+u_{2}{\bm{d}}_{2}+u_{3}{\bm{d}}_{3}$
represent the bending strains about the corresponding directors.

The deformed configuration $\mathscr{S}$ of the ribbon is parameterized as

${\bm{y}}(s,t)={\bm{r}}(s)+\psi(s,t){\bm{q}}(s)\,,$ (26)

where ${\bm{q}}(s)$ is a unit vector given by

${\bm{q}}(s)=\cos\alpha(s){\bm{d}}_{1}(s)+\sin\alpha(s){\bm{d}}_{3}(s)\quad\text{with}-\frac{\pi}{2}\leq\alpha\leq\frac{\pi}{2}$
(27)

and $\psi(s,t)$ is a smooth function. Letting ${\bm{q}}_{0}^{\perp}$ be the
unit vector in the $\\{\bm{e}_{3},\bm{e}_{1}\\}$ plane defined by

${\bm{q}}_{0}^{\perp}:=\bm{e}_{2}\times{\bm{q}}_{0},$ (28)

we represent the deformation gradient $\nabla\bm{y}$ as

$\nabla{\bm{y}}=\bm{a}\otimes{\bm{q}}_{0}+{\bm{b}}\otimes{\bm{q}}_{0}^{\perp}\,,$
(29)

where $\bm{a}$ and ${\bm{b}}$ are the images of ${\bm{q}}_{0}$ and
${\bm{q}}_{0}^{\perp}$ in the current configuration. In the present setting,
(18) amounts to the condition

$|\bm{a}\times{\bm{b}}|=1.$ (30)

It was proved in singh:model  that (30) is equivalently represented by the
following relations,

$\displaystyle\psi$ $\displaystyle=a(s)t\quad\text{with}\quad
a(s)\geq\frac{\cos\alpha_{0}}{v_{3}}>0\,,$ (31a) $\displaystyle\alpha$
$\displaystyle=\arccos\left(\frac{\cos\alpha_{0}}{v_{3}a}\right)\,,$ (31b)
$\displaystyle u_{2}$
$\displaystyle=\partial_{s}\alpha-\frac{\partial_{s}\alpha_{0}}{a^{2}}\,,$
(31c) $\displaystyle u_{3}$ $\displaystyle=u_{1}\tan\alpha\,,$ (31d)

where $a$ and $v_{3}$ are independent deformation measures.

Equations (31) have important consequences (see again singh:model ), namely,

$\partial_{s}{\bm{q}}=-\frac{\partial_{s}\alpha_{0}-\omega_{2}}{a^{2}}{\bm{q}}^{\perp}\,,$
(32)

where ${\bm{q}}^{\perp}:={\bm{d}}_{2}\times{\bm{q}}$, and

$\displaystyle\bm{a}$ $\displaystyle=a{\bm{q}},$ (33a) $\displaystyle{\bm{b}}$
$\displaystyle=-\left[\frac{v_{3}\sin\alpha-a\sin\alpha_{0}+t\partial_{s}a}{\cos\alpha_{0}+t(\partial_{s}\alpha_{0}-\omega_{2})}\right]{\bm{q}}+\frac{1}{a}{\bm{q}}^{\perp}\,.$
(33b)

The latter two equations express $\bm{a}$ and ${\bm{b}}$ in the orthonormal
frame $\\{{\bm{q}},{\bm{q}}^{\perp}\\}$ that spans the local tangent plane at
$\bm{y}(\bm{x})$ to the current configuration of the ribbon. Moreover, in
light of (33a), (29), and (24), (17) implies that

${\bm{q}}=\bm{n},$ (34)

which confers a material character to the rulings of the surface representing
via (26) the current configuration of the ribbon.

It also follows from the inextensibility constraint (30) that the unit normal
$\bm{\nu}$ to $\mathscr{S}$ can be written as

${\bm{\nu}}(s,t)={\bm{q}}(s)\times{\bm{q}}^{\perp}(s)={\bm{d}}_{2}(s)\,,$ (35)

from which we derived in singh:model  the following formula for the curvature
tensor,

$\nabla\\!_{\mathrm{s}}{\bm{\nu}}=\sigma{\bm{q}}^{\perp}\otimes{\bm{q}}^{\perp}\,,$
(36)

where

$\sigma=\frac{a^{2}u_{1}v_{3}}{\cos\alpha_{0}\left[\cos\alpha_{0}+t(\partial_{s}\alpha_{0}-\omega_{2})\right]},$
(37)

so that, in particular, the Gaussian curvature $K$ of $\mathscr{S}$ vanishes
identically, whereas the mean curvature $H=\sigma/2$.

In singh:model , using the kinematics described above along with relations
(31), we obtained a dimension reduction of the sheet elastic energy in (22) to
an energy defined over the centerline of a generic ribbon with a planar
natural configuration. The resulting energy functional is given by

$\mathscr{F}[\bm{y}]=\frac{1}{2}e_{0}h\int_{0}^{L}\\!f\,ds,$ (38)

where $e_{0}$ is the characteristic volume energy density introduced in (10),

$f=F_{0}\ln\left(\frac{\cos\alpha_{0}+t_{+}M}{\cos\alpha_{0}+t_{-}M}\right)+F_{1}(t_{+}-t_{-})+F_{2}(t_{+}^{2}-t_{-}^{2})\,,$
(39)

and

$\displaystyle F_{0}$
$\displaystyle=\left[\frac{8a^{4}h^{2}u_{1}^{2}v_{3}^{2}}{3M\cos\alpha_{0}^{2}}+\frac{2\left(VM-\partial_{s}a\cos\alpha_{0}\right)^{2}}{(1+S)M^{3}}\right]\,,$
(40a) $\displaystyle F_{1}$
$\displaystyle=2\left[\frac{\cos\alpha_{0}}{a^{2}}+\frac{1+S_{0}}{1+S}a^{2}\cos\alpha_{0}+\frac{\partial_{s}a(2VM-\partial_{s}a\cos\alpha_{0})}{(1+S)M^{2}}\right]\,,$
(40b) $\displaystyle F_{2}$
$\displaystyle=\left[\frac{M}{a^{2}}+\frac{1+S_{0}}{1+S}a^{2}M+\frac{(\partial_{s}a)^{2}}{(1+S)M}\right]\,,$
(40c) $\displaystyle V$ $\displaystyle=v_{3}\sin\alpha-a\sin\alpha_{0}\,,$
(40d) $\displaystyle M$ $\displaystyle=\partial_{s}\alpha_{0}-\omega_{2}\,.$
(40e)

##  3 A Narrow Rectangular Ribbon

Our focus in this article is on a simple ribbon with a rectangular geometry,
as shown in Fig. 3.

Figure 3: Schematic of natural and deformed configurations of a narrow NPN
ribbon of length $L$ and width $2w$. The width of the ribbon has been
exaggerated for presentation purposes.

We assume the centerline ${\bm{r}}_{0}$ in the reference configuration to be a
straight line such that it lies along ${\bm{e}}_{3}$, so that
$\omega_{2}\equiv 0$. We further impose a condition on $\alpha_{0}$ such that,

$\alpha_{0}(0)=\alpha_{0}(L)=0\,.$ (41)

which ensures that the shorter edge of the ribbon is along ${\bm{e}}_{1}$, and
that the rulings (23) cover the entire material surface of the natural
configuration. As a consequence of this choice, we have,

$t_{\pm}=\pm\frac{w}{\cos\alpha_{0}}\,.$ (42)

To obtain an energy density for a _narrow_ ribbon, we consider the leading
order term in $w$ in the expansion of (39), under which the energy functional
takes the following form,

$F[\alpha,v_{3}]:=\int_{0}^{1}\left(f_{1}h+f_{3}h^{3}\right)ds\,.$ (43)

Here $F$ is renormalized by the scaling energy $2wL^{2}e_{0}$. The total
length $L$ of the centerline in the reference configuration is taken to be the
length scale, and consequently, all lengths from this point onward will be
stated in units of $L$. Furthermore, $f_{1}$ and $f_{3}$ are given by,

$\displaystyle f_{1}$
$\displaystyle=\frac{v_{3}^{2}\cos^{2}\alpha}{\cos^{2}\alpha_{0}}+\frac{1+S_{0}}{1+S}\frac{\cos^{2}\alpha_{0}}{v_{3}^{2}\cos^{2}\alpha}+\frac{(v_{3}^{2}\cos\alpha\sin\alpha-\cos\alpha_{0}\sin\alpha_{0})^{2}}{(1+S)v_{3}^{2}\cos^{2}\alpha_{0}\cos^{2}\alpha}\,,$
(44a) $\displaystyle\quad f_{3}$
$\displaystyle=\frac{4}{3}\frac{u_{1}^{2}}{v_{3}^{2}\cos^{4}\alpha}\,.$ (44b)

Although $a$ has a transparent geometric meaning, defining $F$ in (43) we
preferred to express it in terms of two other (independent) measures of
deformation, namely, $v_{3}$ and the angle $\alpha$ that $\bm{n}$ makes with
${\bm{d}}_{1}$. They are related to $a$ through the equation

$a=\frac{\cos\alpha_{0}}{v_{3}\cos\alpha},$ (45)

which makes (31b) identically satisfied. Similarly (31c) can be reduced to the
following

$u_{2}=\partial_{s}\alpha-\partial_{s}\alpha_{0}\left(\frac{v_{3}\cos\alpha}{\cos\alpha_{0}}\right)^{2}.$
(46)

Our objective here is to minimize $F$ in (43) under constraints (46) and (31d)
with an appropriate set of boundary conditions leading to out-of-plane bending
of the ribbon. $S_{0}$ is the scalar order parameter frozen in the reference
configuration at the time of cross-linking, while $S$ is the current scalar
order parameter induced by a change in temperature.666It should be recalled
that they are related through (4) to the Maier-Saupe nematic scalar order
parameters $Q_{0}$ and $Q$ in the corresponding configurations of the
material. Since both $S$ and temperature $T$ are held fixed as the spontaneous
deformation unfolds, the minimizers of $F$ are not affected by its being
scaled to a quantity, $e_{0}$, depending on both $S$ and $T$.

##  4 The Variational Problem

We seek solutions with out-of-plane deformation of a narrow rectangular NPN
ribbon on which the following boundary conditions are prescribed,

$\displaystyle{\bm{r}}(1)-{\bm{r}}(0)$ $\displaystyle={\bm{e}}_{3}\,,$ (47a)
$\displaystyle{\bm{d}}_{1}(1)-{\bm{d}}_{1}(0)$ $\displaystyle={\bm{0}}\,,$
(47b)

i.e., we fix the end-to-end position vector, as well as the relative rotation
of ${\bm{d}}_{1}$ between the two ends. Using (25), constraints (47) can be
given the following integral forms,

$\displaystyle\int_{0}^{1}\\!\partial_{s}{\bm{r}}ds=\int_{0}^{1}\\!v_{3}{\bm{d}}_{3}\,ds={\bm{e}}_{3}\,,$
(48a)
$\displaystyle\int_{0}^{1}\\!\partial_{s}{\bm{d}}_{1}\,ds=\int_{0}^{1}(u_{3}{\bm{d}}_{2}-u_{2}{\bm{d}}_{3})\,ds={\bm{0}}\,.$
(48b)

We restrict our interest to the cases where the centerline deforms in the
$\\{{\bm{e}}_{2},{\bm{e}}_{3}\\}$ plane only, meaning that the deformed
centerline can be accorded the following representation,

$\displaystyle{\bm{r}}=x_{2}{\bm{e}}_{2}+x_{3}{\bm{e}}_{3}\,.$ (49)

In line with this assumption of planar deformations, we further restrict the
set of possible deformations with the following conditions on the strains,

$\displaystyle u_{2}=u_{3}=0\,,$ (50)

which makes (31d) satisfied for any $\alpha$. Boundary condition (48b) is
identically satisfied under this assumption, as the only non-zero component
left of ${\bm{u}}$ is $u_{1}$. Equations (50), (31c), (31d), and (41) together
imply the following for non-zero $u_{1}$,

$\displaystyle\alpha(s)=\alpha_{0}(s)=0\,,$ (51)

Next, we parameterize ${\bm{d}}_{3}$ and ${\bm{d}}_{2}$ using an angle
$\vartheta$ as follows,

$\displaystyle{\bm{d}}_{2}$
$\displaystyle=\cos\vartheta{\bm{e}}_{2}-\sin\vartheta{\bm{e}}_{3}\,,$ (52a)
$\displaystyle{\bm{d}}_{3}$
$\displaystyle=\sin\vartheta{\bm{e}}_{2}+\cos\vartheta{\bm{e}}_{3}\,.$ (52b)

It is then easy to deduce using the above and (25)2 the following
parametrization for the only non-zero component of the strain vector,

$u_{1}=\vartheta^{\prime}\,,$ (53)

where a prime ′ denotes differentiation with respect to $s$.

Using (51) and (53), the functions $f_{1}$ and $f_{3}$ in (44) reduce to,

$f_{1}=v_{3}^{2}+\frac{\nu^{2}}{v_{3}^{2}}\quad\text{and}\quad
f_{3}=\frac{4}{3}\left(\frac{\vartheta^{\prime}}{v_{3}}\right)^{2}\,,$ (54)

where

$\nu:=\sqrt{\frac{S_{0}+1}{S+1}}$ (55)

is the effective activation parameter.777Note that $\nu$ defined here is the
reciprocal of $\mu$ defined in singh:model . The convenience of such a choice
will become clear from the appearance of the bifurcation diagrams in Fig. 4.
With this definition of $\nu$, the centerline of the ribbon would be expected
to expand when $\nu>1$. The function $f_{1}$, representing the stretching
energy density in the ribbon, attains its minimum at $v^{*}_{3}=\sqrt{\nu}$
(see also Appendix A). The coupling between bending and stretching in the
ribbon under consideration is clearly delineated in the expression for
$f_{3}$, where the stretch $v_{3}$ appears in the denominator. This functional
dependence of $f_{3}$ on $v_{3}$ effectively reduces the bending rigidity of
the ribbon along the length in a variable fashion.

The energy functional (43) augmented with the constraint (48a) can then be
written as follows in the unknown functions $\vartheta(s)$ and $v_{3}(s)$,

$\displaystyle\widetilde{F}[\vartheta,v_{3}]:=\int_{0}^{1}\left(f_{1}h+f_{3}h^{3}+\Lambda_{2}v_{3}\sin\vartheta+\Lambda_{3}(v_{3}\cos\vartheta-1)\right)ds\,,$
(56)

where $\Lambda_{i}={\bm{\Lambda}}\cdot{\bm{e}}_{i}$, $i\in\\{2,3\\}$, are the
Cartesian components of a constant vector multiplier $\bm{\Lambda}$ employed
to enforce (48a).

Using standard variational procedure, we compute the variation in
$\widetilde{F}[\vartheta,v_{3}]$ due to
$\vartheta\rightarrow\vartheta+\delta\vartheta$ as follows,

$\begin{split}\delta\widetilde{F}(\vartheta,v_{3})[\delta\vartheta]&=\left(\frac{8h^{3}}{3v_{3}^{2}}\vartheta^{\prime}\delta\vartheta\right)\bigg{|}_{0}^{1}\\\
&-\int_{0}^{1}\left\\{\left(\frac{8h^{3}}{3v_{3}^{2}}\vartheta^{\prime}\right)^{\prime}-\Lambda_{2}v_{3}\cos\vartheta+\Lambda_{3}v_{3}\sin\vartheta\right\\}\delta\vartheta\,ds\,,\end{split}$
(57)

delivering the Euler-Lagrange (EL) equation for $\vartheta$ along with the
corresponding boundary conditions,

$\displaystyle\left(\frac{8h^{3}}{3v_{3}^{2}}\vartheta^{\prime}\right)^{\prime}-\Lambda_{2}v_{3}\cos\vartheta+\Lambda_{3}v_{3}\sin\vartheta$
$\displaystyle=0\,,$ (58a)
$\displaystyle\left(\frac{8h^{3}}{3v_{3}^{2}}\vartheta^{\prime}\right)_{s=0}=\left(\frac{8h^{3}}{3v_{3}^{2}}\vartheta^{\prime}\right)_{s=1}$
$\displaystyle=0\,.$ (58b)

Similarly, the EL equation corresponding to $v_{3}$ is given by a
straightforward partial derivative in $v_{3}$ of the integrand of (56),

$\displaystyle\frac{8h^{3}\vartheta^{\prime
2}}{3v_{3}^{3}}+2h\left(\frac{\nu^{2}}{v_{3}^{3}}-v_{3}\right)-\Lambda_{2}\sin\vartheta-\Lambda_{3}\cos\vartheta=0\,.$
(59)

There are no boundary conditions associated with $v_{3}$.

###  4.1 Equilibrium equations

Here we obtain a complete system of first order ODEs representing our
equilibrium problem. The state space of our system is described by the vector
$\\{\vartheta,u_{1},v_{3},x_{2},x_{3},\Lambda_{2},\Lambda_{3}\\}$, with $7$
fields. The evolution of $\vartheta,u_{1},x_{2},x_{3}$ with $s$ can be easily
obtained using (53), (58a), (25)1, (49), and (52b). The two multipliers
$\Lambda_{2}$ and $\Lambda_{3}$ are constant, so their evolution is trivial.
The only field whose evolution needs to be determined is $v_{3}$, as the
corresponding EL equation (59) is algebraic in $v_{3}$. To that end, we
differentiate (59) and use it to resolve for $u^{\prime}_{1}$ and
$v_{3}^{\prime}$ using (58a) (and (53) to eliminate $\vartheta^{\prime}$). As
a result, we arrive at the following set of first order ODEs,

$\displaystyle\vartheta^{\prime}$ $\displaystyle=u_{1}\,,$ (60a)
$\displaystyle u^{\prime}_{1}$
$\displaystyle=\frac{-3v_{3}^{3}(\Lambda_{2}\cos\vartheta-\Lambda_{3}\sin\vartheta)(9\nu^{2}+4h^{2}u_{1}^{2}+3v_{3}^{4})}{8h^{3}\left[4h^{2}u^{2}_{1}-3(3\nu^{2}+v_{3}^{4})\right]}\,,$
(60b) $\displaystyle v_{3}^{\prime}$
$\displaystyle=\frac{-3u_{1}v_{3}^{4}(\Lambda_{2}\cos\vartheta-\Lambda_{3}\sin\vartheta)}{8h^{3}u_{1}^{2}-6h(3\nu^{2}+v_{3}^{4})}\,.$
(60c) $\displaystyle x_{2}^{\prime}$ $\displaystyle=v_{3}\sin\vartheta\,,$
(60d) $\displaystyle x_{3}^{\prime}$ $\displaystyle=v_{3}\cos\vartheta\,.$
(60e) $\displaystyle\Lambda_{2}^{\prime}$ $\displaystyle=0\,,$ (60f)
$\displaystyle\Lambda_{3}^{\prime}$ $\displaystyle=0\,.$ (60g)

These 7 ODEs need to be complemented with $7$ boundary conditions. Two of
these are given by (58b), and the rest are the following 5,

$\begin{split}x_{2}(0)=0\,,\quad x_{2}(1)=0\,,\quad x_{3}(0)=0\,,\quad
x_{3}(1)=1\,,\\\
2h\left(\frac{\nu^{2}}{v_{3}(0)^{3}}-v_{3}(0)\right)-\Lambda_{2}(0)\sin\vartheta(0)-\Lambda_{3}(0)\cos\vartheta(0)=0\,,\end{split}$
(61)

where the last condition ensures that the boundary condition on $v_{3}(0)$ is
consistent with the parent algebraic equation (59) from which (60c) was
derived. Our system is now complete, and can be viewed as a two point boundary
value problem with $\nu$ as a control parameter.

Physical intuition suggests that for $\nu=1$, that is, when $S=S_{0}$, no
spontaneous deformation of the ribbon should take place, as the material is
not activated. In the absence of deformation, the ribbon remains in the
straight configuration, where $\vartheta\equiv 0$, $u_{1}\equiv 0$, and
$v_{3}\equiv 1$, which solve (60) with $\Lambda_{2}=\Lambda_{3}=0$ and
$x_{2}=0$, $x_{3}=s$. It is indeed a simple matter to show that the
configuration

$\begin{split}\vartheta=0,\quad u_{1}=0,\quad
v_{3}=1,\quad\Lambda_{2}=0,\quad\Lambda_{3}=2h(\nu^{2}-1),\\\ x_{2}=0,\quad
x_{3}=s\,,\end{split}$ (62)

is a solution of (60) for all $h>0$ and $\nu>0$. This is the _trivial_
solution of our equilibrium problem and we would like to obtain solutions to
(60) _emanating_ from it.

The system of governing equations (60) is invariant under the transformations
$\vartheta\mapsto-\vartheta,\,u_{1}\mapsto-u_{1},\,v_{3}\mapsto
v_{3},\,\Lambda_{2}\mapsto-\Lambda_{2},\,\Lambda_{3}\mapsto\Lambda_{3},\,x_{2}\mapsto-
x_{2}\,\,x_{3}\mapsto x_{3}$. Thus, for every equilibrium configuration
$\left(x_{3}(s),x_{2}(s)\right)$, its mirror image
$\left(x_{3}(s),-x_{2}(s)\right)$ is also a valid equilibrium configuration.
All equilibrium configurations come in mirror-symmetric pairs.

###  4.2 Bifurcation analysis

Consider the following solution to (60) which differs from the trivial
solution (62) by an infinitesimal amount,

$\begin{split}\vartheta=\varepsilon{\tilde{\vartheta}},\
u_{1}=\varepsilon{\tilde{u}}_{1},\ v_{3}=1+\varepsilon{\tilde{v}}_{3},\
\Lambda_{3}=2h(\nu^{2}-1)+\varepsilon{\tilde{\Lambda}}_{3},\
\Lambda_{2}=\varepsilon{\tilde{\Lambda}}_{2},\\\
x_{2}=\varepsilon{\tilde{x}}_{2},\
x_{3}=s+\varepsilon{\tilde{x}}_{3}\,,\end{split}$ (63)

where $\varepsilon$ is a small non-zero perturbation parameter. Using (58b)
and (61), we deduce the following boundary conditions on the perturbations,

$\displaystyle{\tilde{\vartheta}}^{\prime}(0)={\tilde{\vartheta}}^{\prime}(1)=0,\quad{\tilde{x}}_{2}(0)={\tilde{x}}_{2}(1)=0,\quad{\tilde{x}}_{3}(0)={\tilde{x}}_{3}(1)=0\,,$
(64)

Substituting (64) into (60), and equating equal powers of $\varepsilon$, we
obtain,

$\displaystyle{\tilde{\vartheta}}^{\prime}$ $\displaystyle={\tilde{u}}_{1}\,,$
(65a) $\displaystyle{\tilde{u}}^{\prime}_{1}$
$\displaystyle=\frac{3}{8h^{3}}\left[{\tilde{\Lambda}}_{2}-2h(\nu^{2}-1){\tilde{\vartheta}}\right]\,,$
(65b) $\displaystyle{\tilde{v}}^{\prime}_{3}$ $\displaystyle=0\,,$ (65c)
$\displaystyle{\tilde{x}}^{\prime}_{2}$ $\displaystyle={\tilde{\vartheta}}\,,$
(65d) $\displaystyle{\tilde{x}}^{\prime}_{3}$
$\displaystyle={\tilde{v}}_{3}\,,$ (65e)
$\displaystyle{\tilde{\Lambda}}_{2}^{\prime}$ $\displaystyle=0\,,$ (65f)
$\displaystyle{\tilde{\Lambda}}^{\prime}_{3}$
$\displaystyle=-2h(3\nu^{2}+1){\tilde{v}}_{3}\,.$ (65g)

The boundary conditions on ${\tilde{x}}_{3}$ in (64) imply that
${\tilde{v}}_{3}$ vanishes (from (65c) and (65e)), and so does
${\tilde{\Lambda}}_{3}$ (from (65g)). Using (65a) and (65b) we arrive at the
following second-order linear ODE for ${\tilde{\vartheta}}$,

$\displaystyle{\tilde{\vartheta}}^{\prime\prime}+\frac{3}{4h^{2}}\left(\nu^{2}-1\right){\tilde{\vartheta}}=\frac{3{\tilde{\Lambda}}_{2}}{8h^{3}}\,.$
(66)

For $\nu<1$, the above ODE has no solutions that satisfy the boundary
conditions on ${\tilde{\vartheta}}^{\prime}$ stated in (64). The implication
being that all such configurations are likely to be locally _stable_ , as no
equilibrium solution emanates from it. This was to be expected, as $\nu<1$
implies that $S>S_{0}$, and so the activated configuration is more ordered
than the reference configuration, entailing an extension of the ribbon along
${\bm{d}}_{1}$ (see Fig. 3), and by the surface inextensibility constraint, a
contraction along ${\bm{d}}_{3}$, which is incompatible with boundary
condition (47a).888Here, we effectively regard $\nu$ as a _continuation_
parameter, which drives new equilibrium solutions out of the trivial one.

On the other hand, for $\nu\geq 1$, which is the only case that interests us
here, equation (66) admits non-trivial solutions consistent with the boundary
conditions in (64) when $\nu=\nu_{n}$, where

$\nu_{n}:=\sqrt{1+\frac{4n^{2}\pi^{2}h^{2}}{3}}\quad\text{with}\quad
n\in\mathbb{N}\,.$ (67)

Consequently, ${\tilde{x}}_{2}$ is the _jet_ (eigenfunction) generated by
${\tilde{x}}_{2}=\sin(n\pi s)$. Expression (67) indicates that for a given
branch, identified by an integer value of $n$, the separation between two
adjacent values of $\nu_{n}$ must increase with increasing value of $h$. In
the next section, we quantify this effect, illuminating the interplay of two
_tendencies_ in the mechanical behavior of an activated ribbon, which for
brevity we call “membrane-like” and “plate-like”. The former tendency is
naively expected to prevail for small values of $h$, where the stretching
energy is dominant and the bending energy would be regarded as a perturbation;
the latter, on the contrary, is expected to play a role for larger values of
$h$. We shall see how the activation parameter can be dysfunctional towards
such naive expectations.

##  5 Activated Equilibria

The governing set of equations (60), accompanied by the boundary conditions
(58b) and (61), were solved numerically using a parameter continuation method
doede:numerical , implemented in AUTO-07P auto07p . The bifurcation parameter
was chosen to be the effective activation parameter $\nu$ in (55), which was
varied from $1.0$ to $2.0$, and the trivial solution (62) was given as the
base solution from which the continuation was carried out.

We resolve the degeneracy induced by mirror symmetry by selecting a single
member in each symmetric solution pair. Each equilibrium profile shown below
should ideally be accompanied by its mirror image.

(a)

(b)

(c)

(d)

Figure 4: Plots in panels (4(a)) to (4(d)) show bifurcation diagrams of the
ribbon with increasing thickness. The curves are generated numerically using
AUTO-07P, a continuation and bifurcation software for ODEs, while the circles
denote the analytical predictions of the bifurcation points made using
equation (67). The thick curve in all the plots corresponds to the first mode
of deformation with one crest (see Fig. 5(a)). The spacing between bifurcation
points increases with thickness, as predicted by equation (67).

Our main objective behind performing the forthcoming simulations is to
highlight the interplay between bending and stretching of the ribbon. As we
saw in the expressions for the bending and stretching energies in (54), the
two mechanisms of deformation are coupled via the stretch $v_{3}$ appearing in
the denominator in the expression for $f_{3}$. The following plots reveal and
quantify this non-trivial coupling between the two competing modes of
deformation.

We begin by tracking the evolution of the bifurcation diagrams of the ribbon
as its thickness is increased. The bifurcation plots, i.e. the $L_{2}$ norm of
the function $x_{2}(s)$ vs $\nu$, generated from the numerical simulations for
increasing values of the thickness $h$ are shown in Fig. 4. The bifurcation
points $(\nu_{n},0)$ as predicted by the analytical expression (67) are shown
in the plots using circles, and are in good agreement with the numerical
solutions. Figures 4(a) to 4(d) show the distances between the bifurcation
points increasing with thickness, confirming our assertion made at the end of
Sect. 4.2. This increasing gap between bifurcation points is an indication of
the gradual transformation of the ribbon’s behavior from membrane-like to
plate-like.

(a)

(b)

Figure 5: Plots in panels (5(a)) and (5(b)), respectively, correspond to the
first and the second modes of deformation of a ribbon of thickness $h=0.08$
(in units of $L$), for different values of parameter $\nu$. The total energy
$F$ in (43) of each configuration is shown in the legend below each panel. For
a given value of $\nu$, the energy of the first mode is very close to the
energy of the second mode, however, the latter remains higher than the former
in all cases, and both are lower than the energy $h(1+\nu^{2})$ of the trivial
solution (62).

(a)

(b)

Figure 6: Panel (6(a)): Configurations of a thin NPN ribbon with different
values of $\nu$ and thickness $h$, where each color represents a fixed value
of $\nu$. In all cases, increasing thickness results in decreasing the maximum
displacement of the ribbon. However, this effect appears to be more pronounced
for smaller values of $\nu$. Panel (6(b)): Maximum displacement of the first
mode of deformation of an NPN ribbon plotted against its normalized thickness
for various values of $\nu$.

A series of configurations of the ribbon centerline for different values of
the effective activation parameter $\nu$ are shown in Fig. 5. The
configurations in Fig. 5(a) lie on the first bifurcating branch (the thickened
curve) in Fig. 4(b), while configurations in Fig. 5(b) lie on the second
bifurcating branch (adjacent to the thickened curve) in Fig. 4(b). The total
energy of each configuration is stated in a legend below the two plots in Fig.
5. A curious observation is that for a given value of $\nu$ the energies of
the first and the second mode of deformation are very close, although the
latter is always higher than the former. As expected, the maximum displacement
is larger for the first mode than for the second mode (a rough estimate is
that the former is twice as efficient as the latter). For the first mode, the
maximum displacement also dramatically increases with activation, ranging from
roughly $20\%$ to above $40\%$ of the ribbon’s length, as $\nu$ ranges from
$1.2$ to $2$.

Next we explore the effect of thickness on the maximum displacement of the
ribbon in the first mode. We plot various configurations in Fig. 6(a) for
different values of $\nu$, ranging from $1.2$ to $2.0$, and thickness $h$
ranging from $0.01$ to $0.18$ (in $L$ units). We find that increasing the
thickness of the ribbon reduces its maximum displacement for any given value
of $\nu$. However, as was perhaps to be expected, this effect is much more
pronounced for smaller values of $\nu$, as can be easily observed from Fig
6(a). So much so that for $\nu=2$ the deformed shapes of the ribbon
corresponding to different thicknesses are nearly indiscernible from one
another. We may say that high values of the activation parameter _obliterate_
the distinction between “membrane-like” and “plate-like” behavior of the
ribbon. The dependence of displacement on thickness is quantified in Fig.
6(b), which shows the maximum displacement of the ribbon centerline plotted
against its thickness for several values of the effective activation
parameter.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 7: Several plots depicting the interplay of bending energy
($f_{3}h^{3}$) and stretching energy ($f_{1}h$) of a rectangular NPN ribbon
for different values of $h$ and $\nu$ plotted against $s$. In most panels, the
two energies complement one another such that the local maximum of one
corresponds to the local minimum of the other. This behavior slightly changes
in panels (7(i)), (7(j)), (7(k)), and (7(l)), in that a new local maximum of
the stretching energy forms in the center, where a local maximum of the
bending energy lies.

Finally, we track the distribution of the bending energy vs the stretching
energy along the arc-length coordinate in Fig. 7. The gray curve represents
the bending energy, while the blue curve represents the stretching energy
along the arc-length coordinate. We observe that the stretching energy is
always quantitatively dominant (by at least one order of magnitude) over the
bending energy. However, there appears to be a sort of _complementarity_
relation between the two. Boundary conditions dictate that the bending energy
vanishes at the ribbon’s end-points, whereas the stretching energy can be
freely distributed at equilibrium. Now, although dominant in magnitude, the
latter chooses to be concentrated where the bending energy must be absent and
to be minimum where the bending energy is maximum. Thus, the stretching energy
is primarily concentrated at the extreme ends of the ribbon, while the bending
energy is concentrated in the middle. The stretching energy attains a local
extremum at the center of the ribbon in all cases. This extremum is the
absolute minimum for most of the plots shown. However, for plots in panels
(7(i)), (7(j)), (7(k)), and (7(l)) this extremum morphs curiously into a local
maximum.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 8: Several plots showing the distribution of the curvature $-u_{1}>0$
and the stretch $v_{3}$ of a rectangular NPN ribbon for different values of
$h$ and $\nu$ against $s$. The location of the maxima and minima of the
curvature and the stretch on the centerline coincides for all cases. The
values of $v_{3}$ hover close to the local minimum at $v_{3}^{*}=\sqrt{\nu}$
of the stretching energy function $f_{1}$ in (54).

On similar lines, we plot in Fig. 8 the stretch $v_{3}$ and the curvature
$-u_{1}$ for various values of $\nu$ and $h$. We observe that the curvature
and the stretch both attain their respective maxima at the center of the
ribbon, with their minima located at the two ends. The coincidence of maximum
curvature and maximum stretch was perhaps to be expected, as the bending
rigidity in (54) decreases with increasing $v_{3}$. We also observe that the
variation in the stretch along the length is relatively much smaller, and
$v_{3}$ does not differ much from the value $v_{3}^{*}=\sqrt{\nu}$, shown by
the dashed line on each plot, at which the function $f_{1}$ in (54) attains
its minimum. The stretching measure $v_{3}$ is moderate (corresponding to a
length increment of nearly $20\%$) for small $\nu$ and significant
(corresponding to a length increment of nearly $40\%$) for large $\nu$.

Figure 9: The values of the average stretch $\bar{v}_{3}$, as defined by (68),
are compared with the minimum $v_{3}^{\ast}$ of the stretching energy for the
different values of $h$ and $\nu$ considered in our analysis. The gray
datapoint “$\color[rgb]{.5,.5,.5}{\square}$” has been added for $\nu=1.2$ at
$\bar{v}_{3}=1$, corresponding to the trivial solution in (62), which, as
shown in Fig. 4(d), is the only energy minimizer for that value of $\nu$.

We see from Fig. 9 that, for given $\nu$, an increase in $h$ results into a
decrease in the average stretch defined by

$\bar{v}_{3}:=\int_{0}^{1}v_{3}\,ds.$ (68)

Comparing this latter with $v^{\ast}_{3}$, we see that $\bar{v}_{3}\leq
v^{\ast}_{3}$ for all values of $h$ and $\nu$ that we explored. Moreover, for
given $h$, $\bar{v}_{3}$ and $v^{\ast}_{3}$ get closer as $\nu$ increases and,
for given $\nu$, they become farther apart as $h$ increases (see Fig. 9).
Another telling feature is revealed by Fig. 9: for the smallest value of $h$,
$\bar{v}_{3}$ falls nicely over $v_{3}^{\ast}$, marking the expected
prevalence of a membrane-like behavior, which is soon lost as $h$ is
increased.

Both Fig. 9 and the plots in Figs. 8(j), 8(k), and 8(l) show clearly that, for
large values of $\nu$, $h$ has little effect on the deformation of the ribbon,
in complete accordance with Fig. 6. We shall designate this as the _bleaching_
regime. Our analysis justifies approximating $v_{3}$ with $v^{\ast}_{3}$ in
the bleaching regime: there the membrane-like behavior prevails _irrespective_
of the value of $h$ and the bending energy acts as a _selection criterion_
that singles out a deformation out of myriads with the same minimum stretching
energy. This, however, is no longer the case for smaller values of the
activation parameter (say, near $\nu=1$), where the independent measures of
deformation $(v_{3},u_{1})$ are closely interlaced and the mechanical behavior
of the ribbon cannot be classified neither as “pure-membrane” nor as “pure-
plate”.

##  6 Conclusions

In this paper, which is the ideal companion of singh:model , we studied the
out-of-plane deformations of a narrow ribbon consisting of a nematic polymer
network, a nematic elastomer where the degree of cross-linking is sufficiently
high to justify the enslavement of the nematic director to the polymer matrix
deformation.

The material is activated by a change in temperature: this produces a
mismatch, measured by a control parameter $\nu>1$, between the degrees of
order in the polymer chains organization in the reference and current
configurations. This drives a spontaneous deformation of the ribbon. In our
theory, both stretching and bending energies are contemplated, scaling with
different powers of the ribbon’s thickness $h$.

We deliberately kept simple the deformation pattern, choosing boundary
conditions compatible with the absence of twist, as our primary aim was
exploring the interplay between stretching and bending components of the
elastic energy while both $\nu$ and $h$ were independently varied.

Two major conclusions were reached here. First, stretching and bending
energies appear to obey a _complementarity_ relation: one is maximum where the
other is minimum. Second, our findings point to the existence of a _bleaching_
regime: for $\nu$ sufficiently large, the deformation of the ribbon is
essentially independent of $h$. This latter feature may be relevant to
applications where displacements induced by activation are desired to be large
in not too thin ribbons.

Clearly, envisioning twist in the admissible class of deformations might have
a bearing on our present conclusions, but not—we believe—to the point of
obliterating them. It would be desirable to address this issue in the future.

Acknowledgments H.S. is grateful to S.P.C. Dhanakoti and T. Yu for helpful
discussions regarding AUTO-07P. He also acknowledges partial support by Swiss
National Science Foundation Grant $200020\\_182184$ to J.H. Maddocks.

## Declarations

### Competing interests

The authors have no conflicts of interest to declare that are relevant to the
content of this article.

##  Appendix A Minimizing stretch

Here we show that the stretching content $f_{\mathrm{s}}(\mathbf{C})$ in (20a)
has precisely one stationary point $\mathbf{C}_{0}$ subject to the
inextensibility constraint

$\det\mathbf{C}=1$ (69)

generated by (18).

Since

$\frac{\partial}{\partial\mathbf{C}}\det\mathbf{C}=(\det\mathbf{C})(\mathbf{C}^{\mathsf{T}})^{-1},$
(70)

see for example (gurtin:mechanics, , p. 24), by (69) and the symmetry of
$\mathbf{C}$, it follows from (20a) that the stationary condition

$\frac{\partial
f_{\mathrm{s}}}{\partial\mathbf{C}}=\lambda\frac{\partial}{\partial\mathbf{C}}\det\mathbf{C}$
(71)

is equivalently written as

$\mathbf{I}_{2}+\left(S_{0}-\frac{S}{(\bm{n}_{0}\cdot\mathbf{C}\bm{n}_{0})^{2}}\right)\bm{n}_{0}\otimes\bm{n}_{0}=\lambda\mathbf{C}^{-1},$
(72)

where $\mathbf{I}_{2}$ is the two-dimensional identity and $\lambda$ is the
Lagrange multiplier associated with (69).

Acting with $\mathbf{C}$ on both sides of (72), we arrive at

$\mathbf{C}+\left(S_{0}-\frac{S}{(\bm{n}_{0}\cdot\mathbf{C}\bm{n}_{0})^{2}}\right)\mathbf{C}\bm{n}_{0}\otimes\bm{n}_{0}=\lambda\mathbf{I}_{2},$
(73)

which easily implies that a unit vector $\bm{n}_{0}^{\perp}$ perpendicular to
$\bm{n}_{0}$ in the $\\{\bm{e}_{3},\bm{e}_{1}\\}$ plan is an eigenvector of
$\mathbf{C}$ with eigenvalue $\lambda$, while $\bm{n}_{0}$ is an eigenvector
with eigenvalue $1/\lambda$, by (69). Taking the inner product with
$\bm{n}_{0}\otimes\bm{n}_{0}$ of both sides of (73), we readily arrive at

$\lambda=\nu,$ (74)

where $\nu$ is as in (55), and we write the minimizer $\mathbf{C}_{0}$ of
$f_{\mathrm{s}}$ as999That the _only_ stationary point $\mathbf{C}_{0}$
actually minimizes $f_{\mathrm{s}}$ follows easily from being $f_{\mathrm{s}}$
bounded below (but not above).

$\mathbf{C}_{0}=\frac{1}{\nu}\bm{n}_{0}\otimes\bm{n}_{0}+\nu\bm{n}_{0}^{\perp}\otimes\bm{n}_{0}^{\perp}.$
(75)

Clearly, $\mathbf{C}_{0}$ reduces to $\mathbf{I}_{2}$ for $\nu=1$ (that is,
for $S=S_{0}$).

## References

  *   * (1) Corbett, D., Modes, C.D., Warner, M.: Photomechanics: Bend, curl, topography, and topology. In: White, T.J. (ed.) Photomechanical Materials, Composites, and Systems. Wireless Transduction of Light Into Work, pp. 79–116. John Wiley & Sons, Hoboken, NJ (2017)
  * (2) White, T.J.: Photomechanical effects in liquid-crystalline polymer networks and elastomers. In: White, T.J. (ed.) Photomechanical Materials, Composites, and Systems. Wireless Transduction of Light Into Work, pp. 153–177. John Wiley & Sons, Hoboken, NJ (2017)
  * (3) Ware, T.H., White, T.J.: Programmed liquid crystal elastomers with tunable actuation strain. Polym. Chem. 6, 4835–4844 (2015). https://doi.org/10.1039/C5PY00640F
  * (4) Sonnet, A.M., Virga, E.G.: Model for photoresponsive nematic elastomers. J. Elast. (2022). https://doi.org/10.1007/s10659-022-09959-4
  * (5) Bladon, P., Terentjev, E.M., Warner, M.: Deformation-induced orientational transitions in liquid crystals elastomer. J. Phys. II France 4(1), 75–91 (1994). https://doi.org/10.1051/jp2:1994100
  * (6) Warner, M., Gelling, K.P., Vilgis, T.A.: Theory of nematic networks. J. Chem. Phys. 88(6), 4008–4013 (1988). https://doi.org/10.1063/1.453852
  * (7) Warner, M., Wang, X.J.: Elasticity and phase behavior of nematic elastomers. Macromolecules 24(17), 4932–4941 (1991). https://doi.org/10.1021/ma00017a033
  * (8) Warner, M., Mostajeran, C.: Nematic director fields and topographies of solid shells of revolution. Proc. R. Soc. London A 474(2210), 20170566 (2018). https://doi.org/10.1098/rspa.2017.0566
  * (9) Warner, M., Terentjev, E.M.: Liquid Crystal Elastomers. International Series of Monographs on Physics, vol. 120. Oxford University Press, New York (2003)
  * (10) Warner, M., Bladon, P., Terentjev, E.M.: “Soft elasticity”—deformation without resistance in liquid crystal elastomers. J. Phys. II France 4(1), 93–102 (1994). https://doi.org/10.1051/jp2:1994116
  * (11) Terentjev, E.M., Warner, M., Bladon, P.: Orientation of nematic elastomers and gels by electric fields. J. Phys. II France 4(4), 667–676 (1994). https://doi.org/10.1051/jp2:1994154
  * (12) Verwey, G.C., Warner, M.: Soft rubber elasticity. Macromolecules 28(12), 4303–4306 (1995). https://doi.org/10.1021/ma00116a036
  * (13) Verwey, G.C., Warner, M.: Multistage crosslinking of nematic networks. Macromolecules 28(12), 4299–4302 (1995). https://doi.org/10.1021/MA00116A035
  * (14) Verwey, G.C., Warner, M., Terentjev, E.M.: Elastic instability and stripe domains in liquid crystalline elastomers. J. Phys. II France 6(9), 1273–1290 (1996). https://doi.org/10.1051/jp2:1996130
  * (15) Anderson, D.R., Carlson, D.E., Fried, E.: A continuum-mechanical theory for nematic elastomers. J. Elast. 56, 33–58 (1999). https://doi.org/10.1023/A:1007647913363
  * (16) Zhang, Y., Xuan, C., Jiang, Y., Huo, Y.: Continuum mechanical modeling of liquid crystal elastomers as dissipative ordered solids. J. Mech. Phys. Solids 126, 285–303 (2019). https://doi.org/10.1016/j.jmps.2019.02.018
  * (17) Mihai, L.A., Wang, H., Guilleminot, J., Goriely, A.: Nematic liquid crystalline elastomers are aeolotropic materials. Proc. R. Soc. London A 477(2253), 20210259 (2021). https://doi.org/10.1098/rspa.2021.0259
  * (18) White, T.J. (ed.): Photomechanical Materials, Composites, and Systems: Wireless Transduction of Light Into Work. John Wiley & Sons, Hoboken, New Jersey (2017)
  * (19) Korley, L.T.J., Ware, T.H.: Introduction to special topic: Programmable liquid crystal elastomers. J. Appl. Phys. 130(22), 220401 (2021). https://doi.org/10.1063/5.0078455
  * (20) Mahimwalla, Z., Yager, K.G., Mamiya, J.-i., Shishido, A., Priimagi, A., Barrett, C.J.: Azobenzene photomechanics: prospects and potential applications. Polym. Bull. 69, 967–1006 (2012). https://doi.org/10.1007/s00289-012-0792-0
  * (21) Ube, T., Ikeda, T.: Photomobile polymer materials with crosslinked liquid-crystalline structures: Molecular design, fabrication, and functions. Angew. Chem. Int. Ed. 53(39), 10290–10299 (2014). https://doi.org/10.1002/anie.201400513
  * (22) White, T.J.: Photomechanical effects in liquid crystalline polymer networks and elastomers. J. Polym. Sci. Part B: Polym. Phys. 56(9), 695–705 (2018). https://doi.org/%****␣Bending_Ribbon_SN_rev_galleys.tex␣Line␣1500␣****10.1002/polb.24576
  * (23) Ula, S.W., Traugutt, N.A., Volpe, R.H., Patel, R.R., Yu, K., Yakacki, C.M.: Liquid crystal elastomers: an introduction and review of emerging technologies. Liquid Cryst. Rev. 6(1), 78–107 (2018). https://doi.org/10.1080/21680396.2018.1530155
  * (24) Pang, X., Lv, J.-a., Zhu, C., Qin, L., Yu, Y.: Photodeformable azobenzene-containing real polymers and soft actuators. Adv. Mater. 31(52), 1904224 (2019). https://doi.org/10.1002/adma.201904224
  * (25) Kuenstler, A.S., Hayward, R.C.: Light-induced shape morphing of thin films. Curr. Opin. Colloid & Interface Sci. 40, 70–86 (2019). https://doi.org/10.1016/j.cocis.2019.01.009
  * (26) Warner, M.: Topographic mechanics and applications of liquid crystalline solids. Annu. Rev. Condens. Matter Phys. 11(1), 125–145 (2020). https://doi.org/10.1146/annurev-conmatphys-031119-050738
  * (27) Singh, H., Virga, E.G.: A ribbon model for nematic polymer networks. J. Elast. (2022). https://doi.org/10.1007/s10659-022-09900-9
  * (28) Ozenda, O., Sonnet, A.M., Virga, E.G.: A blend of stretching and bending in nematic polymer networks. Soft Matter 16, 8877–8892 (2020). https://doi.org/10.1039/D0SM00642D
  * (29) Finkelmann, H., Greve, A., Warner, M.: The elastic anisotropy of nematic elastomers. Eur. Phys. J. E 5, 281–293 (2001). https://doi.org/10.1007/s101890170060
  * (30) Nguyen, T.-S., Selinger, J.V.: Theory of liquid crystal elastomers and polymer networks. Eur. Phys. J. E 40, 76 (2017). https://doi.org/10.1140/epje/i2017-11569-5
  * (31) Corbett, D., Warner, M.: Nonlinear photoresponse of disordered elastomers. Phys. Rev. Lett. 96, 237802 (2006). https://doi.org/10.1103/PhysRevLett.96.237802
  * (32) Corbett, D., Warner, M.: Polarization dependence of optically driven polydomain elastomer mechanics. Phys. Rev. E 78, 061701 (2008). https://doi.org/10.1103/PhysRevE.78.061701
  * (33) Treloar, L.R.G.: The Physics of Rubber Elasticity, 3rd edn. Oxford Classic Texts in the Physical Sciences. Oxford University Press, Oxford (2005)
  * (34) Deam, R.T., Edwards, S.F.: The theory of rubber elasticity. Philos. Trans. R. Soc. London A 280, 317–353 (1976). https://doi.org/10.1098/rsta.1976.0001
  * (35) Warner, M.: New elastic behaviour arising from the unusual constitutive relation of nematic solids. J. Mech. Phys. Solids 47, 1355–1377 (1999). https://doi.org/10.1016/S0022-5096(98)00100-8
  * (36) White, T.J., Broer, D.J.: Programmable and adaptive mechanics with liquid crystal polymer networks and elastomers. Nature Mater. 14(2210), 1087–1098 (2015). https://doi.org/10.1038/nmat4433
  * (37) Modes, C.D., Bhattacharya, K., Warner, M.: Disclination-mediated thermo-optical response in nematic glass sheets. Phys. Rev. E 81, 060701 (2010). https://doi.org/10.1103/PhysRevE.81.060701
  * (38) Podio-Guidugli, P.: An exact derivation of the thin plate equation. J. Elast. 22, 121–133 (1989). https://doi.org/10.1007/BF00041107
  * (39) Ozenda, O., Virga, E.G.: On the Kirchhoff-Love hypothesis (revised and vindicated). J. Elast. 143, 359–384 (2021). https://doi.org/10.1007/s10659-021-09819-7
  * (40) Doedel, E., Keller, H.B., Kernevez, J.P.: Numerical analysis and control of bifurcation problems (II): Bifurcation in infinite dimensions. Int. J. Bif. Chaos 1, 745–772 (1991). https://doi.org/10.1142/S0218127491000555
  * (41) Doedel, E.J., Fairgrieve, T.F., Sandstede, B., Champneys, A.R., Kuznetsov, Y.A., Wang, X.: AUTO-07P: Continuation and bifurcation software for ordinary differential equations (2007). http://indy.cs.concordia.ca/auto/
  * (42) Gurtin, M.E., Fried, E., Anand, L.: The Mechanics and Thermodynamics of Continua. Cambridge University Press, Cambridge (2010)