Ridge energy for thin nematic polymer networks

Abstract

Minimizing the elastic free energy of a thin sheet of nematic polymer network among smooth isometric immersions is the strategy purported by the mainstream theory. In this paper, we broaden the class of admissible spontaneous deformations: we consider ridged isometric immersions, which can cause a sharp ridge in the immersed surfaces. We propose a model to compute the extra energy distributed along such ridges. This energy comes from bending; it is shown under what circumstances it scales quadratically with the sheet’s thickness, falling just in between stretching and bending energies. We put our theory to the test by studying the spontaneous deformation of a disk on which a radial hedgehog was imprinted at the time of crosslinking. We predict the number of folds that develop in terms of the degree of order induced in the material by external agents (such as heat and illumination).

Graphic Abstract

1 Introduction

Nematic elastomers are rubber materials with a fluid-like component constituted by elongated, rod-like molecules appended to the crosslinked polymer strands that form the background matrix. The fluid component is ordered as nematic liquid crystals can be, which makes these solid materials very susceptible to external stimuli, such as heat, light, and environmental humidity. The prompt response to these stimuli, so characteristic of liquid crystals, once transferred to the solid matrix, makes it possible to do work and change the shape of bodies with no direct contact. The possible technological applications of these materials are boundless (see, for example, the papers [1,2,3,4,5,6,7,8,9], and above all the review [10]), but a number of theoretical challenges remain open [11]; this paper is concerned with one of them.

The order established in the material by the mutual interaction of nematic molecules is described by a scalar-order parameter, representing the degree of molecular alignment, and a director, representing the average direction of alignment. Actually, there are two sets of these order parameters, namely, the pair $(s_0,\varvec{m})$ for the reference configuration of the rubber matrix, which here will be taken to be the configuration where the crosslinking takes place, and the pair $(s,\varvec{n})$ for the current (deformed) configuration, the one the rubber matrix takes on in response to an applied stimulus (more details are given in Sect. 2). The director $\varvec{n}$ can be tied to the deformation of the body in several ways, the spectrum going from complete independence to complete enslaving. Following the terminology introduced in [10], we call nematic polymer networks the nematic elastomers in which the crosslinking in so tight that the nematic director remains enslaved to the deformation;¹ these are the specific nematic elastomers treated here. The reason for this choice will soon become clear.

The most striking manifestation of the ability of nematic polymer networks to perform changes in shape is perhaps achieved when they are thin sheets. We represent one such sheet as a slab ${\mathsf {S}}$ of thickness 2h extending in the reference configuration on both sides of a flat surface S. The director $\varvec{m}$ is blueprinted on S (in its own plane), uniformly reproduced across the thickness, with a given scalar-order parameter $s_0$. External stimuli may act on the degree of order, changing $s_0$ into s, in a programmable way. The system is thus carried out of equilibrium and a deformation ensues, for the free energy to attain a minimum under the changed circumstances.

An elastic free-energy density, $f_e$, is available for bulk materials in three-space dimensions since the pioneering work [17] (a comprehensive introduction to the subject is offered by the landmark book [18]); it is delivered by the “trace formula,” derived from assuming an anisotropic Gaussian distribution for the polymer chains that constitute the rubber matrix.² This formula features both the deformation $\varvec{f}$ of the three-dimensional body ${\mathscr {B}}$ occupied by the material and measures of anisotropy in both reference and current configurations of ${\mathscr {B}}$ (see Sect. 2.) For a sufficiently thin slab ${\mathsf {S}}$, however, one’s desire is to reduce $f_e$ to a function of the mapping $\varvec{y}$ that only changes the flat reference mid-surface S into a curved surface ${\mathscr {S}}$ in the current configuration.

In a nematic polymer network, for which $f_e$ eventually depends only on $\varvec{f}$, such a dimension reduction was performed in [21] by revisiting (and extending) a standard method of the theory of plates, known as the Kirchhoff–Love hypothesis [22]. As expected, this method delivers a surface elastic energy with two components, a stretching energy $f_s$ scaling like h, and a bending energy $f_b$ scaling like $h^3$; $f_s$ depends only on the two-dimensional stretching (or metric) tensor ${\mathbf {C}}:=(\nabla \varvec{y})^{\mathsf {T}}(\nabla \varvec{y})$, while $f_b$ also depends on the invariant measures of curvature of ${\mathscr {S}}$ and the relative orientation of $\varvec{n}$ in the frame of principal directions of curvature. Not only do $f_s$ and $f_b$ scale differently with h, they are also basically different things. By Gauss’ theorema egregium [23, p. 139], the Gaussian curvature K of ${\mathscr {S}}$ is fully determined by the metric tensor ${\mathbf {C}}$, thus deserving the name of intrinsic curvature. As a consequence, $f_s$ depends only on the intrinsic curvature, whereas $f_b$ also depends on extrinsic measures of curvature, relating on how ${\mathscr {S}}$ is embedded in three-dimensional space. Studying the complete equilibrium problem, where $f_b$ is treated on the same footing as $f_s$ has proven so far difficult. A number of strategies have been devised to circumvent the energy coupling.

For moderately curved surfaces ${\mathscr {S}}$ and sufficiently thin slabs ${\mathsf {S}}$, for which $f_b$ can be neglected relative to $f_s$, the energy minimizing shapes are isometric immersions³ of the metric tensor ${\mathbf {C}}_0$ that minimizes $f_s$. The search for such immersions corresponding to a variety of imprinted $\varvec{m}$ fields has been the subject of a vast, elegant literature (see, among others, [15, 24,25,26,27,28,29,30]) This may seem to solve the direct morphic mechanics problem for nematic polymer networks, namely, how to identify the shapes produced by a certain imprinted director field $\varvec{m}$. More difficult (and less visited), but affordable is the inverse problem of assigning $\varvec{m}$ so as to produce a desired shape upon stimulation [31].

As reassuring as this picture may appear, things are unfortunately more complicated than they look like: there are at least two conflicting, as it were, unresolved issues. A smooth isometric immersion with prescribed metric tensor ${\mathbf {C}}_0$ may altogether fail to exist in the large. On the other hand, if we renounce the smoothness requirement for the immersion, the number of admissible solutions may easily become too large.⁴

A remedy for the first issue was proposed by the theory of geometric elasticity [32, 33]. If the target metric corresponding to ${\mathbf {C}}_0$ is geometrically incompatible with a smooth immersion, this theory proposes to replace it with the one that minimizes an appropriate $L^2$-distance from it. It is a viable approximation, if you do not wish to renounce regularity.

A remedy for the second issue would be provided by a selection criterion that single out one shape out of many, preferably on energetic grounds. Here the essential question is: what extra energy should be attached to a singular shape? This is the avenue taken here. We allow ${\mathscr {S}}$ to have ridges, that is, lines along which the outer unit normal $\varvec{\nu }$ suffers a discontinuity. As for the extra energy cost to be associated with a ridge, we extract it from the bending energy density $f_b$. We conceive a ridge as a limiting tight fold, for which we justify an expression for a ridge linear density $f_r$, which depends (in a symmetric way) on the traces $\varvec{\nu }_1$ and $\varvec{\nu }_2$ of the unit normal $\varvec{\nu }$ on both sides of the ridge. It turns out that in our theory $f_r$ scales like $h^2$, just in between $f_s$ and $f_b$, so that $f_r$ becomes the effective substitute for $f_b$. This justifies an approximation alternative to geometric elasticity: finding piecewise isometric $C^2$-immersions that minimize the total ridge energy.

The vicarious nature of our theory is to be stressed from the start. The real (still unresolved) challenge is minimizing the total elastic energy of a thin nematic polymer network, with both stretching and bending contributions. Failing to do so, we find it expedient to replace a distributed bending energy with a concentrated one, which is simpler than the former (and scales differently too), but is not ad hoc. Were the same replacement adopted for the Euler elastica, one would find a similarly viable theory.

The paper is organized as follows. In Sect. 2, we recall both stretching and bending energies for nematic polymer networks, as they emerged from the dimension reduction of the bulk energy density delivered by the trace formula. Section 3 plays the role of an intermezzo in our development: there we show how to destructure the classical elastica in a chain of rigid rods connected by articulated joints (edges) encapsulating the bending energy of the parent body. We shall see how this simplified model is capable of capturing the known qualitative behavior of the elastica, thus paving the way to our ridge construction. In Sect. 4, we construct the ridge energy as a limit of the bending energy entrapped in a folded sheet. Section 5 is concerned with the general equations that govern piecewise $C^2$-immersions with discontinuities of the unit normal field $\varvec{\nu }$ concentrated along smooth curves; such ridged isometric immersions are the shapes competing for a minimum in our theory. In Sect. 6, we consider some special symmetric ridged isometric immersions that mimic the folds generated in a disk when the imprinted director field $\varvec{m}$ is the radial hedgehog; we compute the ridge energy that acts as an obstruction to the proliferation of folds and we determine their optimal number. Section 7 is where we draw our conclusions and comment on other possible uses of our theory. The paper is closed by an appendix, where we illustrate a geometric construction apt to produce the analytic solution proposed in Sect. 6 for the ridged isometric immersion of a hedgehog.

2 Stretching and bending energies

In this section, we recall the outcomes of the dimension reduction method applied in [21] to the trace formula of the neo-classical theory for nematic elastomers (for which we refer the reader to Chap. 6 of [18]). Two director fields feature in this theory; these are $\varvec{m}$, defined in the reference configuration ${\mathscr {B}}$ of the body, and $\varvec{n}$, defined in the current configuration $\varvec{f}({\mathscr {B}})$ obtained from ${\mathscr {B}}$ through the deformation $\varvec{f}$. ${\mathscr {B}}$ is a region in three-dimensional Euclidean space ${\mathscr {E}}$ and $\varvec{f}:{\mathscr {B}}\rightarrow {\mathscr {E}}$ is a diffeomorphism of ${\mathscr {B}}$.

The directors $\varvec{m}$ and $\varvec{n}$ represent the average alignment of the elongated molecules appended to the rubber polymeric matrix in the reference and current configurations. They are properly defined through the tensorial measures of anisotropy that characterize the end-to-end Gaussian distribution of polymer strands. These are the polymer step tensors ${\mathbf {L}}_{\varvec{m}}$ and ${\mathbf {L}}_{\varvec{n}}$, in the reference and current configurations, respectively, which, following [34] and [35], we write as

$$\begin{aligned} {\mathbf {L}}_{\varvec{m}}:={\mathfrak {a}}_0({\mathbf {I}}+s_0\varvec{m}\otimes \varvec{m}) \end{aligned}$$

(1a)

and

$$\begin{aligned} {\mathbf {L}}_{\varvec{n}}:={\mathfrak {a}}({\mathbf {I}}+s\varvec{n}\otimes \varvec{n}). \end{aligned}$$

(1b)

Here ${\mathbf {I}}$ is the identity (in three-dimensional space), ${\mathfrak {a}}_0$ and ${\mathfrak {a}}$ are fixed positive parameters (representing the persistence lengths perpendicular to $\varvec{m}$ and $\varvec{n}$, respectively), $s_0$ and s are nematic scalar-order parameters, which can be expressed as $s_0=r_0-1$ and $s=r-1$ in terms of the ratios $r_0$ and r of the parallel (along $\varvec{m}$ and $\varvec{n}$) and perpendicular (across $\varvec{m}$ and $\varvec{n}$) step chain lengths in the reference and current configurations, respectively.

The neo-classical theory of nematic elastomers expresses the elastic free-energy density $f_e$ (per unit volume in the reference configuration) as

$$\begin{aligned} f_e:=\frac{1}{2}k{\text {tr}}({\mathbf {F}}^{\mathsf {T}}{\mathbf {L}}_{\varvec{n}}^{-1}{\mathbf {F}}{\mathbf {L}}_{\varvec{m}}), \end{aligned}$$

(2)

where ${\mathbf {F}}:=\nabla \varvec{f}$ is the deformation gradient and $k>0$ is an elastic modulus (which scales linearly with both absolute temperature and number density of polymer chains). This is usually called the trace formula.

In nematic elastomers, $\varvec{n}$ and ${\mathbf {F}}$ are fully independent. In contrast, in nematic polymer networks, $\varvec{n}$ is enslaved to ${\mathbf {F}}$. In these materials, with which we are concerned in this paper, the director field $\varvec{m}$ is blueprinted in the elastic matrix [36] and conveyed by the deformation into $\varvec{n}$, which is thus delivered by

$$\begin{aligned} \varvec{n}=\frac{{\mathbf {F}}\varvec{m}}{|{\mathbf {F}}\varvec{m}|}. \end{aligned}$$

(3)

In general, elastomers are incompressible, and so ${\mathbf {F}}$ must satisfy

$$\begin{aligned} \det {\mathbf {F}}=1. \end{aligned}$$

(4)

Both (3) and (4) will be enforced as constraints on all admissible deformations $\varvec{f}$ of ${\mathscr {B}}$.

With $\varvec{m}$ (and $s_0$) imprinted in the reference configuration at the time of crosslinking and $\varvec{n}$ enslaved to the deformation, the only residual freedom lies with s, which can be changed by either thermal or optical stimuli. For example, by heating the sample above the crosslinking temperature, we reduce the nematic order of the chains, so that $s<s_0$; this in turn induces a spontaneous deformation so as to minimize the total elastic free energy. Thus, s can be regarded as the activation parameter of our theory, driven by external stimuli. For definiteness, we shall assume that both $s_0$ and s range in the interval $(-1,1)$.

It is shown in [21] that by use of (1) and (3) $f_e$ can be given the following form

$$\begin{aligned} f_e=\frac{1}{2}k\frac{{\mathfrak {a}}_0}{{\mathfrak {a}}}F({\mathbf {C}}_{\varvec{f}}), \end{aligned}$$

(5)

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

$$\begin{aligned} F({\mathbf {C}}_{\varvec{f}})={\text {tr}}{\mathbf {C}}_{\varvec{f}}+\frac{s_0}{s+1}\varvec{m}\cdot {\mathbf {C}}_{\varvec{f}}\varvec{m}-\frac{s}{s+1}\frac{\varvec{m}\cdot {\mathbf {C}}_{\varvec{f}}^2\varvec{m}}{\varvec{m}\cdot {\mathbf {C}}_{\varvec{f}}\varvec{m}}.\nonumber \\ \end{aligned}$$

(6)

The properties of this function will illuminate the role of s as activation parameter.

As a consequence of (4), ${\mathbf {C}}_{\varvec{f}}$ is also subject to the constraint

$$\begin{aligned} \det {\mathbf {C}}_{\varvec{f}}=1. \end{aligned}$$

(7)

The tensors ${\mathbf {C}}_{\varvec{f}}$ that make $F({\mathbf {C}}_{\varvec{f}})$ stationary subject to (7) are solutions to the equation

$$\begin{aligned} \frac{\partial F}{\partial {\mathbf {C}}_{\varvec{f}}}=\lambda \frac{\partial }{\partial {\mathbf {C}}_{\varvec{f}}}(\det {\mathbf {C}}_{\varvec{f}}), \end{aligned}$$

(8)

where $\lambda$ is a Lagrange multiplier. It is not difficult to see that this equation reduces to

$$\begin{aligned}&{\mathbf {I}}+\frac{1}{s+1}\left( s_0+s\frac{\varvec{m}\cdot {\mathbf {C}}_{\varvec{f}}^2\varvec{m}}{(\varvec{m}\cdot {\mathbf {C}}_{\varvec{f}}\varvec{m})^2}\right) \varvec{m}\otimes \varvec{m}\nonumber \\&\quad -\frac{s}{s+1}\frac{1}{\varvec{m}\cdot {\mathbf {C}}_{\varvec{f}}\varvec{m}}({\mathbf {C}}_{\varvec{f}}\varvec{m}\otimes \varvec{m}+\varvec{m}\otimes {\mathbf {C}}_{\varvec{f}}\varvec{m})\nonumber \\&\quad =\lambda {\mathbf {C}}_{\varvec{f}}^{-1}. \end{aligned}$$

(9)

It follows from (9) that changing $\varvec{m}$ into ${\mathbf {Q}}\varvec{m}$, for any orthogonal tensor ${\mathbf {Q}}$, transforms a solution ${\mathbf {C}}_{\varvec{f}}$ into ${\mathbf {Q}}{\mathbf {C}}_{\varvec{f}}{\mathbf {Q}}^{\mathsf {T}}$, which makes any solution ${\mathbf {C}}_{\varvec{f}}$ of (9) an isotropic tensor-symmetric-valued function of $\varvec{m}$. By the representation theorem of such functions [37] and (7), we know that ${\mathbf {C}}_{\varvec{f}}$ must have the form

$$\begin{aligned} {\mathbf {C}}_{\varvec{f}}=\lambda _f^2\varvec{m}\otimes \varvec{m}+\frac{1}{\lambda _f}({\mathbf {I}}-\varvec{m}\otimes \varvec{m}), \end{aligned}$$

(10)

for some $\lambda _f\in {\mathbb {R}}^+$. Making use of (10) in (9), we readily conclude that

$$\begin{aligned} \lambda =\frac{1}{\lambda _f}\quad \text {and}\quad \lambda _f=\root 3 \of {\frac{s+1}{s_0+1}}. \end{aligned}$$

(11)

By expressing $F({\mathbf {C}}_{\varvec{f}})$ in terms of $\lambda _f$ with the aid of (10), it is easy to show that for $\lambda _f$ as in (11) this function attains its unique minimum.

Thus, when $s<s_0$, the spontaneous deformation induced in the material would be a contraction along $\varvec{m}$, accompanied by a dilation in the plane orthogonal to $\varvec{m}$, to preserve the volume.⁵ Of course, it remains to be seen whether, for an assigned $\varvec{m}$, a deformation with a metric that minimizes F locally is indeed geometrically compatible in the large; differently put, whether there is an isometric immersion in three-space dimensions of the desired target metric ${\mathbf {C}}_{\varvec{f}}$ as in (10).

Here we are interested in thin sheets and in the appropriate dimension reduction of $F({\mathbf {C}}_{\varvec{f}})$ to be attributed to the mid-surface S of the slab ${\mathsf {S}}$ of thickness 2h. Formally, S is a flat region in the $(x_1,x_2)$ plane of a fixed Cartesian frame $(\varvec{e}_1,\varvec{e}_2,\varvec{e}_3)$ and ${\mathsf {S}}$ is the set in three-space defined as ${\mathsf {S}}:=\{(\varvec{x},x_3)\in S\times [-h,h] \}$. The mapping $\varvec{y}:S\rightarrow {\mathscr {E}}$ describes the deformation of S into the surface ${\mathscr {S}}=\varvec{y}(S)$ in the deformed slab $\varvec{f}({\mathsf {S}})$; we shall assume that $\varvec{y}$ is of class $C^2$ and that $\varvec{m}$ is a two-dimensional field imprinted on S, so that $\varvec{m}\cdot \varvec{e}_3\equiv 0$ (see Fig. 1).⁶

Fig. 1

A flat surface S in the $(x_1,x_2)$ plane of a fixed Cartesian frame $(\varvec{e}_1,\varvec{e}_2,\varvec{e}_3)$ is deformed by the mapping $\varvec{y}$ into a smooth surface ${\mathscr {S}}$ embedded in three-dimensional Euclidean space ${\mathscr {E}}$. The blueprinted orientation is denoted by $\varvec{m}$ in the reference configuration and by $\varvec{n}$ in the current one; $\varvec{e}_3$ is the outer unit normal to S, while $\varvec{\nu }$ is the outer unit normal to ${\mathscr {S}}$; correspondingly, $\varvec{m}_\perp :=\varvec{e}_3\times \varvec{m}$ and $\varvec{n}_\perp :=\varvec{\nu }\times \varvec{n}$

The (two-dimensional) deformation gradient has the following general representation,

$$\begin{aligned} \nabla \varvec{y}=\varvec{a}\otimes \varvec{m}+\varvec{b}\otimes \varvec{m}_\perp , \end{aligned}$$

(12)

where $\varvec{m}_\perp :=\varvec{e}_3\times \varvec{m}$. In (12), $\varvec{a}$ and $\varvec{b}$ are vector fields defined on S; they live in ${\mathscr {V}}$, the translation space of ${\mathscr {E}}$, and are everywhere tangent to ${\mathscr {S}}$. It follows from (12) that the two-dimensional stretching tensor ${\mathbf {C}}$ is represented as

$$\begin{aligned} {\mathbf {C}}= & {} (\nabla \varvec{y}^{\mathsf {T}})(\nabla \varvec{y})=a^2\varvec{m}\otimes \varvec{m}+\varvec{a}\cdot \varvec{b}(\varvec{m}\otimes \varvec{m}_\perp +\varvec{m}_\perp \otimes \varvec{m})\nonumber \\&+b^2\varvec{m}_\perp \otimes \varvec{m}_\perp , \end{aligned}$$

(13)

where $a^2:=\varvec{a}\cdot \varvec{a}$ and $b^2:=\varvec{b}\cdot \varvec{b}$. We shall require that S is inextensible, which amounts to the constraint $|\varvec{a}\times \varvec{b}|=1$. Thus, since $\det {\mathbf {C}}=a^2b^2-(\varvec{a}\cdot \varvec{b})^2=|\varvec{a}\times \varvec{b}|^2$, we shall require that

$$\begin{aligned} \det {\mathbf {C}}=1. \end{aligned}$$

(14)

Under this constraint, the outer unit normal $\varvec{\nu }$ to ${\mathscr {S}}$ will be delivered by

$$\begin{aligned} \varvec{\nu }=\varvec{a}\times \varvec{b}. \end{aligned}$$

(15)

Applying (3) to the present setting, we obtain that

$$\begin{aligned} \varvec{n}=\frac{(\nabla \varvec{y})\varvec{m}}{|(\nabla \varvec{y})\varvec{m}|}, \end{aligned}$$

(16)

and so we may write $\varvec{a}=a\varvec{n}$ and define $\varvec{n}_\perp :=\varvec{\nu }\times \varvec{n}$, so that the frame $(\varvec{n},\varvec{n}_\perp ,\varvec{\nu })$ is oriented as $(\varvec{m},\varvec{m}_\perp ,\varvec{e}_3)$ (see Fig. 1).

In [21], we extend the classical Kirchhoff–Love hypothesis [22] to obtain a dimension reduction of $F({\mathbf {C}}_{\varvec{f}})$ in (6), that is, a method that converts $f_e$ in (5) into a surface energy-density (to be integrated over S). As standard in the theory of plates, such a surface energy is delivered by a polynomial in odd powers of h, conventionally truncated so as to retain the first two relevant ones, the first and the third power. The former is the stretching energy $f_s$, accounting for the work done to alter distances and angles in S, while the latter is the bending energy $f_b$, accounting for the work done to fold S. Thus, dropping the scaling constant $\frac{1}{2}k\frac{{\mathfrak {a}}_0}{{\mathfrak {a}}}$, which has the physical dimensions of an energy per unit volume, we can write

$$\begin{aligned} f_e=f_s+f_b+O(h^5), \end{aligned}$$

(17)

where (to within an inessential additive constant)

$$\begin{aligned} f_s= & {} \frac{2h}{s+1}\left( {\text {tr}}{\mathbf {C}}+s_0\varvec{m}\cdot {\mathbf {C}}\varvec{m}+\frac{s}{\varvec{m}\cdot {\mathbf {C}}\varvec{m}}\right) , \end{aligned}$$

(18a)

$$\begin{aligned} f_b= & {} \frac{2h^3}{3}\left\{ 2(8H^2-K){+}\frac{1}{s+1}\left[ \left( \frac{3s}{a^2}{-}a^2s_0{-}{\text {tr}}{\mathbf {C}}\right) K\right. \right. \nonumber \\&\left. \left. -\frac{4s}{a^2}(2H-\kappa _n)\kappa _n\right] \right\} . \end{aligned}$$

(18b)

Here H and K are the mean and Gaussian curvatures of ${\mathscr {S}}$, defined as

$$\begin{aligned} H:=\frac{1}{2}{\text {tr}}(\nabla \!_\mathrm {s}\varvec{\nu })\quad \text {and}\quad K:=\det (\nabla \!_\mathrm {s}\varvec{\nu })\end{aligned}$$

(19)

in terms of the (two-dimensional) curvature tensor $\nabla \!_\mathrm {s}\varvec{\nu }$, and

$$\begin{aligned} \kappa _n:=\varvec{n}\cdot (\nabla \!_\mathrm {s}\varvec{\nu })\varvec{n}. \end{aligned}$$

(20)

The (scaled) total elastic free energy then reduces to the functional

$$\begin{aligned} {\mathscr {F}}[\varvec{y}]:=\int _S(f_s+f_b)\mathrm {d}A, \end{aligned}$$

(21)

where A is the area measure.

A perturbation approach to the minimization of ${\mathscr {F}}$ is justified when the length scale associated with the average radius of curvature of ${\mathscr {S}}$ is large compared to h, which is the smallest length in the system: then $f_s$ and $f_b$ are well scale-separated and the latter can be viewed as a higher-order correction to the former. In this approach, it is justified to ask what stretching tensor ${\mathbf {C}}_0$ would minimize $f_s$, the leading term in ${\mathscr {F}}$. The answer is easily obtained [21], that is,

$$\begin{aligned} {\mathbf {C}}_0=\lambda _1^2\varvec{m}\otimes \varvec{m}+\lambda _2^2\varvec{m}_\perp \otimes \varvec{m}_\perp , \end{aligned}$$

(22)

where

$$\begin{aligned} \lambda _1:=\root 4 \of {\frac{s+1}{s_0+1}}\quad \text {and}\quad \lambda _2=\frac{1}{\lambda }_1. \end{aligned}$$

(23)

A deformation $\varvec{y}$ for which (22) is valid is an isometric immersion; it minimizes the (leading) stretching energy. The problem is then whether such immersions do exist and how many they are.⁷ This is when the bending energy comes to play. If there are no isometric immersions, it means that $f_s$ must be blended with $f_b$ and more elaborate minimizing shapes ${\mathscr {S}}$ must be sought for, presumably exhibiting regions where the average radius of curvature is not much larger than h. On the other hand, if there are many isometric immersions, we may hope to use the bending energy as a selection criterion, choosing the isometric immersion with the least bending energy.

Both scenarios, however, are overoptimistic. The first, because minimizing the blended energy is not an easy task, also numerically, as the functional ${\mathscr {F}}$ depends on the second as well as the first gradient of $\varvec{y}$. The second, because the two-step minimization, which unleashes $f_b$ over the minimizers of $f_s$, may actually turn out to be rather disappointing; for example, only spheres are allowed among surfaces ${\mathscr {S}}$ with positive K, if one insists in minimizing $f_b$ uniformly [21].

So far we have thought of isometric immersions as smooth mappings. The regularity issue now becomes relevant and tackling it opens up new perspectives. As shown in Sect. 5, one can easily incur in many an isometric immersion by relaxing the requirement that $\varvec{y}$ be $C^2$. We shall consider mappings $\varvec{y}$ that are piecewise $C^1$, with $\nabla \varvec{y}$ allowed to jump across one or several ridges, which will be assumed to be smooth curves of class $C^1$. One such mapping is a piecewise isometric immersion if $(\nabla \varvec{y})^{\mathsf {T}}(\nabla \varvec{y})\equiv {\mathbf {C}}_0$ on the whole domain S, despite the discontinuities of $\nabla \varvec{y}$ across ridges.

The issue with such ridged immersions is that they effectively encapsulate a bending energy in the ridges across which the outer unit normal to ${\mathscr {S}}$ jumps abruptly. In Sect. 4, by regarding each of these ridges as a tight fold with continuous principal curvatures extending over a fixed length $\varepsilon \sim h$, we shall extract out of $f_b$ an elastic ridge energy-density (per unit length) $f_r$. This energy, which scales like $h^2$, will replace $f_b$ in the purpose of mitigating the multiplicity of ridged isometric immersions. It will form the basis of our (simplified) model for nematic polymer networks. Before all this, to motivate better our moderately unritual approach to isometric immersions, we pause briefly and apply a similar approach to a well established elastic problem, that of the elastica.⁸

3 Intermezzo: disembodied elastica

As is well known, the elastica is a one-dimensional continuum body represented by an inextensible curve in space endowed with bending stiffness. We may think of this as the mid-line of a thin three-dimensional body, extended in one direction much more than in the other two. In its simplest incarnation, the elastic energy stored in the elastica is

$$\begin{aligned} {\mathscr {E}}_b=\frac{1}{2}B\int _0^L\kappa ^2\mathrm {d}\xi , \end{aligned}$$

(24)

where $B>0$ is the bending modulus, L is the length of the midline, $\xi$ the arc-length parameter and $\kappa$ the curvature.

Fig. 2

Sketchy justification of the reduction of ${\mathscr {E}}_b$ to ${\mathscr {E}}_e$. The fibers of an elastica are bent over a fixed length $\varepsilon$ along a circle of radius R. We extract from them the relevant portion of bending energy ${\mathscr {E}}_b$ and attribute it to an edge where the unit tangent $\varvec{t}$ suffers an abrupt disalignment, from $\varvec{t}_1$ to $\varvec{t}_2$, measured by the angle $\Delta \theta =\varepsilon /R$

Now, instead of distributing the bending energy all along the midline, we concentrate it in a finite number of places, N. In each of these, as shown in Fig. 2, we imagine to explore the three-dimensional body at the length scale of the diameter 2h of its (circular) cross section. The midline is bent by an angle $\Delta \theta$ along a circular arc of fixed length $\varepsilon$, bearing an elementary edge energy

$$\begin{aligned} \Delta E:=\frac{1}{2}\frac{B}{\varepsilon }(\Delta \theta )^2, \end{aligned}$$

(25)

which easily follows from (24). Next, we replace the elastica by a chain of N rigid rods, each of length L/N, connected to the adjoining ones by an edge, to which we assign the energy $\Delta E$ in (25), with $\varepsilon$ fixed and $\Delta \theta$ expressing the mismatch of the concurring rods (see Fig. 2). In a way, here we go backwards along Hencky’s route [38]. We extract from the energy of Euler’s elastica the energy that Hencky’s model attributes to an articulated system of rigid rods connected through torsional springs.⁹ Admittedly, this is a cruder model, which is however amenable to a simple numerical study revealing the main qualitative features of classical solutions, as we now proceed to show by example.

First, we write the total energy of the chain as

$$\begin{aligned} {\mathscr {E}}_e:=\frac{1}{2}\frac{B}{h}\sum _{i=1}^{N-1}\arccos ^2(\varvec{t}_{i+1}\cdot \varvec{t}_i), \end{aligned}$$

(26)

where we have set $\varepsilon =h$, which represents the smallest length scale in the model. Here $\varvec{t}_i$ is the unit tangent vector along the i-th rod in the chain. In particular, we want to study the equilibrium problem of a chain whose first and last rods are clamped one on top of the other at a distance $a<L$ (see Fig. 3), so that, in a Cartesian frame $(\varvec{e}_x,\varvec{e}_y,\varvec{e}_z)$, $\varvec{t}_1=\varvec{t}_N=\varvec{e}_z$.

Fig. 3

A gallery of absolute energy minimizers for ${\mathscr {E}}_e$ in (26) subject to the constraints (29) are computed for $N=10$ and decreasing values of the (normalized) end-to-end separation $\alpha$ of a disembodied clamped elastica. Red balls represent the clamped ends, while blue balls represent the articulated edges where two adjacent rods concur. The buckling transition to a twisted (out-of-plane) configuration can be placed in the interval $0.6<\alpha <0.7$. The system is seen from different perspectives in different panels to help the reader visualizing the twisted configurations. The Python code (with adjustable parameters $\alpha$ and N) that produced these stable equilibrium configurations can be found (and run) in [44]. All pictures were drawn using VESTA software [45]

Letting $\varvec{t}_i$ be represented in spherical coordinates as

$$\begin{aligned} \varvec{t}_i=\sin \theta _i\cos \phi _i\varvec{e}_x+\sin \theta _i\sin \phi _i\varvec{e}_y+\cos \theta _i\varvec{e}_z \end{aligned}$$

(27)

with $\theta _i\in [0,\pi ]$ and $\phi _i\in [0,2\pi )$, we easily give ${\mathscr {E}}_e$ in (26) the following form

$$\begin{aligned} {\mathscr {E}}_e= & {} \frac{1}{2}\frac{B}{h}\sum _{i=1}^{N-1}\arccos ^2\left[ (\sin \theta _{i+1}\sin \theta _i\cos (\phi _{i+1}-\phi _i)\right. \nonumber \\&\left. +\cos \theta _{i+1}\cos \theta _i\right] , \end{aligned}$$

(28)

subject to the constraints

$$\begin{aligned}&\frac{1}{N}\sum _{i=1}^N\cos \theta _i=\frac{a}{L}=:\alpha ,\nonumber \\&\sum _{i=1}^N\sin \theta _i\cos \phi _i=0,\quad \sum _{i=1}^N\sin \theta _i\sin \phi _i=0, \end{aligned}$$

(29)

which (with $\theta _1=\theta _N=0$) ensure that the boundary conditions are met.

This chain of articulated rods, which pays an elastic disalignment cost at the edges, is our disembodied elastica. We minimized numerically ${\mathscr {E}}_e$ in (28) subject to (29) for decreasing values of $0<\alpha <1$; we used a stochastic gradient descent method for which a Python code is available in [44] (and can be run with parameters $\alpha$ and N of the user’s choice). The outcomes of our computations are shown in Fig. 3 for $N=10$. Upon decreasing $\alpha$, we first find the absolute minimum of ${\mathscr {E}}_e$ in a planar configuration (one of the infinitely many possible was selected with a bias for $\phi =0$), but as soon as $\alpha$ becomes smaller than a critical value $\alpha _\mathrm {c}$ the absolute minimum of ${\mathscr {E}}_e$ is attained on either of two (equally energetic) out-of-plane configurations exhibiting a spontaneous twist, which persists and grows upon further reducing $\alpha$. Such a buckling instability, which injects chirality into the system, was already found in [46] by a bifurcation analysis of the equilibrium equations for the parent elastica (illustrated, in particular, in their Figs. 3 and 4); in our parameterization, their bifurcation point lies at $\alpha _\mathrm {c}\doteq 0.63$ (corresponding to the case where the torsional stiffness vanishes, as is implicit in (24)). Our estimate ($0.6<\alpha _\mathrm {c}<0.7$) is clearly approximate, but the qualitative agreement between disembodied and full fleshed elasticae is undeniable.

Reassured by this agreement, achieved with just a small number of rods, in the following section, we shall propose a similar simplified representation for the bending energy of nematic polymer networks.¹⁰ What here are edges, will there be ridges.

4 Ridge energy

In this section, we describe how we envision what at large scale is a ridge on ${\mathscr {S}}$: this is generated by what at short scale is a sharp bend of ${\mathsf {S}}$; we shall derive from $f_b$ in (18b) the energy that can be associated with it. We start from the large-scale perspective. Let C be a smooth (plane) curve on S (say, of class $C^1$) splitting S in two sides, $S_1$ and $S_2$ and let $\varvec{e}$ be a unit tangent vector to C (see Fig. 4).

Fig. 4

The reference surface S is split by a smooth curve in two sides, $S_1$ and $S_2$, which a deformation $\varvec{y}$, continuous through C, but with discontinuous gradient, maps into the sides ${\mathscr {S}}_1$ and ${\mathscr {S}}_2$ of the ridge ${\mathscr {C}}$. The unit tangent vector $\varvec{t}$ to ${\mathscr {C}}$ is related by (31) to the unit tangent vector $\varvec{e}$ to C

A deformation $\varvec{y}:S\rightarrow {\mathscr {E}}$, continuous across C but with discontinuous gradient $\nabla \varvec{y}$, must obey the following kinematic compatibility condition,

$$\begin{aligned} {\llbracket {\nabla \varvec{y}}}\rrbracket \varvec{e}=\varvec{0}, \end{aligned}$$

(30)

where the jump ${\llbracket {(\cdot )}}\rrbracket :=(\cdot )_2-(\cdot )_1$ is taken on the two sides $S_2$ and $S_1$ of C. The deformation $\varvec{y}$ transforms C into a smooth curve ${\mathscr {C}}$ on ${\mathscr {S}}$ along which the outer unit normal $\varvec{\nu }$ is discontinuous; we call ${\mathscr {C}}$ a ridge of ${\mathscr {S}}$. We shall denote by $\varvec{\nu }_1$ and $\varvec{\nu }_2$ the traces of $\varvec{\nu }$ taken on the sides ${\mathscr {S}}_1$ and ${\mathscr {S}}_2$, connecting along ${\mathscr {C}}$ and corresponding to the sides $S_1$ and $S_2$ of C, respectively.¹¹

We designate by $\varvec{\gamma }$ the parameterization of ${\mathscr {C}}$ in the arc-length $\ell$ and correspondingly we call $\varvec{t}(\ell )=\varvec{\gamma }'(\ell )$ its unit tangent vector; $\varvec{t}$ is related to $\varvec{e}$ through the equation

$$\begin{aligned} \varvec{t}=\frac{(\nabla \varvec{y})_i\varvec{e}}{|(\nabla \varvec{y})_i\varvec{e}|}, \end{aligned}$$

(31)

where by (30) i can take either value $i=1,2$, according to the side upon which ${\mathscr {C}}$ is approached.

The mapping $\varvec{y}$ is a ridged isometric immersion if

$$\begin{aligned} {\llbracket {(\nabla \varvec{y})^{\mathsf {T}}(\nabla \varvec{y})}}\rrbracket =\varvec{0}, \end{aligned}$$

(32)

meaning that the stretching tensor ${\mathbf {C}}$ is continuous across C.¹² The field $\varvec{m}$ is taken to be continuous across C, but by (16) $\varvec{n}$ generally fails to be continuous across ${\mathscr {C}}$. However, because of the identities

$$\begin{aligned}&{\llbracket {(\nabla \varvec{y})\varvec{m}\cdot (\nabla \varvec{y})\varvec{m}}}\rrbracket =0,\quad {\llbracket {(\nabla \varvec{y})\varvec{e}\cdot (\nabla \varvec{y})\varvec{e}}}\rrbracket =0,\nonumber \\&\quad \text {and}\quad {\llbracket {(\nabla \varvec{y})\varvec{e}\cdot (\nabla \varvec{y})\varvec{m}}}\rrbracket =0, \end{aligned}$$

(33)

which follow from (32) and the continuity of $\varvec{m}$, we arrive at

$$\begin{aligned} {\llbracket {\varvec{n}\cdot \varvec{t}}}\rrbracket =0, \end{aligned}$$

(34)

so that the projection of $\varvec{n}$ along the ridge must be the same on both its sides.

Now, we turn to the short-scale perspective. Much in tune with the geometric construction employed in disembodying the elastica in the preceding section, we imagine that a ridge ${\mathscr {C}}$ results from a sharp bend in the mid-surface ${\mathscr {S}}$ of $\varvec{f}({\mathsf {S}})$ extending over a fixed length $\varepsilon$ comparable with the shortest length scale h in the model. To establish a clear connection between such a short-scale bent surface and the large-scale ridge ${\mathscr {C}}$, we need digress slightly.

Consider a tube surface ${\mathscr {T}}$ (see, for example, pp. 649–650 of [47]) generated by the motion of a sphere of (possibly variable) radius R whose center travels along a curve ${\mathscr {C}}^*$, see Fig. 5a.

Fig. 5

The tube construction for a single ridge

Such a surface can be parameterized as follows,

$$\begin{aligned}&\varvec{p}(\xi ,\theta )=\varvec{\gamma }^*(\xi )+R(-\cos \theta \varvec{n^*}(\xi )+\sin \theta \varvec{b}^*(\xi ))\nonumber \\&\text {for}\quad 0\leqq \xi \leqq L^*\quad \text {and}\quad 0\leqq \theta \leqq 2\pi , \end{aligned}$$

(35)

where $\varvec{\gamma }^*$ is the parameterization of ${\mathscr {C}}^*$ in the arc-length $\xi$, and $(\varvec{t}^*,\varvec{n^*},\varvec{b}^*)$ designates its Frenet–Serret frame.

The situation we envision is shown in Fig. 5b. The tube surface ${\mathscr {T}}$ connects two components, ${\mathscr {S}}_1$ and ${\mathscr {S}}_2$, of the mid-surface ${\mathscr {S}}$ over a length $\varepsilon$, which is related to the (finite) angle $\Delta \theta =|\theta _2-\theta _1|$ delimited by the connecting normals, $\varvec{\nu }_1$ and $\varvec{\nu }_2$, through

$$\begin{aligned} \varepsilon =R\Delta \theta . \end{aligned}$$

(36)

At this scale, we identify the curve ${\mathscr {C}}$ that will serve as ridge at a coarser scale by taking the intersection on the $(\varvec{n^*},\varvec{b}^*)$ plane between the lines orthogonal to $\varvec{\nu }_1$ and $\varvec{\nu }_2$, as shown in Fig. 5b. In this construction, the radius R of the traveling sphere, as well as the angles $\theta _1$ and $\theta _2$ can be taken as (smooth) functions of the arc-length parameter $\xi$; for the derivatives of these functions, we shall assume that both

$$\begin{aligned} \varepsilon |\theta _1'|\ll 1\quad \text {and}\quad \varepsilon |\theta _2'|\ll 1. \end{aligned}$$

(37)

In the parameters $(\xi ,\theta )$, the outer unit normal $\varvec{\nu }$ to ${\mathscr {T}}$ reads as

$$\begin{aligned} \varvec{\nu }=-\cos \theta \varvec{n^*}(\xi )+\sin \theta \varvec{b}^*(\xi ), \end{aligned}$$

(38)

and we can represent ${\mathscr {C}}$ as (see Fig. 5b)

$$\begin{aligned} \varvec{\gamma }=\varvec{\gamma }^*+\frac{R}{\cos \left( \frac{\Delta \theta }{2}\right) }(-\cos {\overline{\theta }}\varvec{n^*}+\sin {\overline{\theta }}\varvec{b}^*), \end{aligned}$$

(39)

where we have set ${\overline{\theta }}:=\frac{1}{2}(\theta _1+\theta _2)$. It is now a simple matter to show that, if in addition to (37) we also assume that the curvature $\kappa ^*$ and torsion $\tau ^*$ of ${\mathscr {C}}^*$ are such that

$$\begin{aligned} \varepsilon \kappa ^*\ll 1\quad \text {and}\quad \varepsilon |\tau ^*|\ll 1, \end{aligned}$$

(40)

the curves ${\mathscr {C}}$ and ${\mathscr {C}}^*$ are nearly parallel and the curvature $\kappa$ of ${\mathscr {C}}$ can be identified with the curvature $\kappa ^*$ of ${\mathscr {C}}^*$: they differ by terms vanishing with $\varepsilon$ as do the arc-length parameters $\ell$ and $\xi$.¹³ With this identification, it is not difficult to see that the curvature tensor of ${\mathscr {T}}$ is given by (see also [47, p. 650])

$$\begin{aligned} \nabla \!_\mathrm {s}\varvec{\nu }=\frac{\kappa \cos \theta }{1+R\kappa \cos \theta }\varvec{t}\otimes \varvec{t}+\frac{1}{R}\varvec{t}_\perp \otimes \varvec{t}_\perp , \end{aligned}$$

(41)

where $\varvec{t}_\perp :\varvec{\nu }\times \varvec{t}$, so that

$$\begin{aligned} H= & {} \frac{1}{2}\left( \frac{1}{R}+\frac{\kappa \cos \theta }{1+R\kappa \cos \theta }\right) \quad \text {and}\nonumber \\ K= & {} \frac{1}{R}\frac{\kappa \cos \theta }{1+R\kappa \cos \theta }. \end{aligned}$$

(42)

The area element is correspondingly delivered by

$$\begin{aligned} \mathrm {d}A=(1+R\kappa \cos \theta )\mathrm {d}\ell \mathrm {d}\sigma , \end{aligned}$$

(43)

where $\mathrm {d}\sigma =R\mathrm {d}\theta$.

The aim of this construction is to extract from $f_b$ the bending energy concentrated in a jump of $\varvec{\nu }$ that is to be assigned to ${\mathscr {C}}$ as an energy distributed over its length. To this end, we recall that, at the leading order in $\varepsilon \kappa$, $f_b$ in (18b) can be estimated from (41) as

$$\begin{aligned}&H=\frac{\Delta \theta }{2\varepsilon } +O(\varepsilon \kappa ),\quad K=\frac{\Delta \theta }{\varepsilon }\kappa \cos \theta +O(\varepsilon \kappa ),\nonumber \\&\text {so that}\quad \frac{K}{H^2}=O(\varepsilon \kappa ), \end{aligned}$$

(44a)

$$\begin{aligned}&\text {and}\quad (2H-\kappa _n)\kappa _n=\left( \frac{\Delta \theta }{\varepsilon }\right) ^2(\varvec{n}\cdot \varvec{t})^2(\varvec{n}\cdot \varvec{t}_\perp )^2+O(\varepsilon \kappa ). \end{aligned}$$

(44b)

Making use of (44) in (18b), we arrive at

$$\begin{aligned} f_b= & {} \frac{8}{3}\frac{h^3}{\varepsilon ^2}\left( 1-\frac{s}{s+1}\frac{1}{a^2}(\varvec{n}\cdot \varvec{t})^2(\varvec{n}\cdot \varvec{t}_\perp )^2\right) (\Delta \theta )^2\nonumber \\&+O(\varepsilon \kappa ), \end{aligned}$$

(45)

where we recall that $a^2=\varvec{m}\cdot {\mathbf {C}}\varvec{m}$. Integrating this in the tube delimited by $\theta _1$ and $\theta _2$, since by (34) and the constraint $\varvec{n}\cdot \varvec{n}=1$ both $(\varvec{n}\cdot \varvec{t})^2$ and $(\varvec{n}\cdot \varvec{t}_\perp )^2$ are continuous across ${\mathscr {C}}$ and can be taken as independent of $\sigma$ over the tube, by (43) we estimate a single ridge energy ${\mathscr {F}}_r$ as

$$\begin{aligned} {\mathscr {F}}_r:= & {} \int _0^L\int _0^\varepsilon f_b\mathrm {d}\ell \mathrm {d}\sigma \nonumber \\= & {} \frac{8}{3}h^2\int _0^L \left( 1{-}\frac{s}{s{+}1}\frac{1}{a^2}(\varvec{n}\cdot \varvec{t})^2 (\varvec{n}\cdot \varvec{t}_\perp )^2\right) (\Delta \theta )^2\mathrm {d}\ell \nonumber \\&+\ O(h^2\kappa ), \end{aligned}$$

(46)

where, on physical grounds, we have taken $\varepsilon =h$, for h is the smallest length scale meaningful in our model. This justifies reducing ${\mathscr {F}}_r$ to a line integral along ${\mathscr {C}}$ with density (per unit length)

$$\begin{aligned} f_r:{=}\frac{8}{3}h^2\arccos ^2(\varvec{\nu }_1\cdot \varvec{\nu }_2)\left( 1{-}\frac{s}{s{+}1}\frac{1}{a^2}(\varvec{n}\cdot \varvec{t})^2(\varvec{n}\cdot \varvec{t}_\perp )^2\right) , \end{aligned}$$

(47)

which is the main outcome of our tube construction.

If, at a length scale larger than h, ${\mathscr {S}}_1$ and ${\mathscr {S}}_2$ are isometrically immersed and meet at the ridge ${\mathscr {C}}$, this latter is endowed with the extra energy ${\mathscr {F}}_r$ in (46). In case of multiple ridges ${\mathscr {C}}_j$, we shall simply assume that

$$\begin{aligned} {\mathscr {F}}_r:=\sum _{j=1}^N\int _{{\mathscr {C}}_j}f_r\mathrm {d}\ell , \end{aligned}$$

(48)

where N is the total number of ridges present on ${\mathscr {S}}$. This is the total ridge energy that we shall assign here to a ridged isometric immersion.

Two comments are in order. First, by direct inspection of (47), it is evident that for $s>0$ the ridge energy-density would promote an alignment of $\varvec{n}$ at $\frac{\pi }{4}$ with the ridge (on both adjoining sides), whereas for $s<0$ it would equally promote an alignment either parallel or orthogonal to the ridge. Second, and more importantly, since ${\mathscr {F}}_r$ scales like $h^2$, it dominates over the bending energy distributed over the smooth components ${\mathscr {S}}_k$ of ${\mathscr {S}}$ isometrically immersed in three-dimensional space. Thus, in our theory ${\mathscr {F}}_r$ becomes the effective substitute for the bending energy.¹⁴

In the following section, we shall write the equations that describe a ridged immersion in a special representation. An example of ${\mathscr {F}}_r$ will be computed explicitly in Sect. 6.

5 Representing ridged isometric immersions

Away from possible point defects (the only ones allowed here), the director field $\varvec{m}$ imprinted on S and its orthogonal companion $\varvec{m}_\perp$ are assumed to have continuous gradients, which can be represented as

$$\begin{aligned} \nabla \varvec{m}=\varvec{m}_\perp \otimes \varvec{c}\quad \text {and}\quad \nabla \varvec{m}_\perp =-\varvec{m}\otimes \varvec{c}\end{aligned}$$

(49)

in terms of the planar connector field $\varvec{c}$ [21]. We shall use the frame $(\varvec{m},\varvec{m}_\perp ,\varvec{e}_3)$ to represent a deformation $\varvec{y}$ of S,

$$\begin{aligned} \varvec{y}=y_1\varvec{m}+y_2\varvec{m}_\perp +y_3\varvec{e}_3, \end{aligned}$$

(50)

where $y_i$ are smooth scalar fields on S, so that by (49)

$$\begin{aligned} \nabla \varvec{y}= & {} y_1\varvec{m}_\perp \otimes \varvec{c}+\varvec{m}\otimes \nabla y_1-y_2\varvec{m}\otimes \varvec{c}+\varvec{m}_\perp \otimes \nabla y_2\nonumber \\&+\varvec{e}_3\otimes \nabla y_3. \end{aligned}$$

(51)

Letting, similarly, $\varvec{c}=c_1\varvec{m}+c_2\varvec{m}_\perp$, we easily see that the vectors $\varvec{a}$ and $\varvec{b}$ in (12) can be given the representation

$$\begin{aligned} \varvec{a}= & {} (\nabla \varvec{y})\varvec{m}=(y_{1,1}-c_1y_2)\varvec{m}+(y_{2,1}+c_1y_1)\varvec{m}_\perp +y_{3,1}\varvec{e}_3,\nonumber \\ \end{aligned}$$

(52a)

$$\begin{aligned} \varvec{b}= & {} (\nabla \varvec{y})\varvec{m}_\perp =(y_{1,2}-c_2y_2)\varvec{m}+(y_{2,2}+c_2y_1)\varvec{m}_\perp +y_{3,2}\varvec{e}_3,\nonumber \\ \end{aligned}$$

(52b)

where we have used the expressions

$$\begin{aligned} \nabla y_1= & {} y_{1,1}\varvec{m}+y_{1,2}\varvec{m}_\perp +y_{1,3}\varvec{e}_3, \end{aligned}$$

(53a)

$$\begin{aligned} \nabla y_2= & {} y_{2,1}\varvec{m}+y_{2,2}\varvec{m}_\perp +y_{2,3}\varvec{e}_3, \end{aligned}$$

(53b)

$$\begin{aligned} \nabla y_3= & {} y_{3,1}\varvec{m}+y_{3,2}\varvec{m}_\perp +y_{3,3}\varvec{e}_3. \end{aligned}$$

(53c)

Now, also in view of (13), we see that requiring $\varvec{y}$ in (50) to be an isometric immersion satisfying (22) reduces to enforcing the following conditions

$$\begin{aligned} a^2=\lambda _1^2,\quad b^2=\lambda ^2_2,\quad \varvec{a}\cdot \varvec{b}=0. \end{aligned}$$

(54)

These, with the aid of (52), read explicitly as

$$\begin{aligned}&(y_{1,1}-c_1y_2)^2+(y_{2,1}+c_1y_1)^2+y_{3,1}^2=\lambda _1^2, \qquad \end{aligned}$$

(55a)

$$\begin{aligned}&(y_{1,2}-c_2y_2)^2+(y_{2,2}+c_2y_1)^2+y_{3,2}^2=\lambda _2^2, \qquad \end{aligned}$$

(55b)

$$\begin{aligned}&(y_{1,1}{-}c_1y_2)(y_{1,2}{-}c_2y_2)\nonumber \\&\quad +\,(y_{2,1}{+}c_1y_1)(y_{2,2}{+}c_2y_1){+}y_{3,1}y_{3,2}=0, \end{aligned}$$

(55c)

which constitute a nonlinear system of PDEs for the unknown functions $y_1$, $y_2$, and $y_3$. As a consequence of Gauss’ theorema egregium, an isometric immersion characterized by (22) has Gaussian curvature dictated by $\varvec{m}$ through the equation [21, 26]

$$\begin{aligned} K=\left( \lambda _1^2-\lambda _2^2\right) (c_2^2-c_1^2+c_{12}), \end{aligned}$$

(56)

where we have set $c_{12}=\varvec{m}\cdot (\nabla \varvec{c})\varvec{m}_\perp$.

Moreover, for a ridged isometry, Eq. (55) must be supplemented by the form appropriate to this setting of the jump condition in (30). Since both $\varvec{m}$ and $\varvec{m}_\perp$ are continuous across any plane curve C (with unit tangent $\varvec{e}$), by (12), (30) becomes

$$\begin{aligned} (\varvec{m}\cdot \varvec{e}){\llbracket {\varvec{a}}}\rrbracket +(\varvec{m}_\perp \cdot \varvec{e}){\llbracket {\varvec{b}}}\rrbracket =\varvec{0}. \end{aligned}$$

(57)

Letting $\varvec{e}=\cos \chi \varvec{m}+\sin \chi \varvec{m}_\perp$, since both $\varvec{y}$ and $\varvec{c}$ are continuous across C, (57) reduces to the three scalar equations

$$\begin{aligned}&{\llbracket {y_{1,1}}}\rrbracket \cos \chi +{\llbracket {y_{1,2}}}\rrbracket \sin \chi =0,\nonumber \\&{\llbracket {y_{2,1}}}\rrbracket \cos \chi +{\llbracket {y_{2,2}}}\rrbracket \sin \chi =0,\nonumber \\&{\llbracket {y_{3,1}}}\rrbracket \cos \chi +{\llbracket {y_{3,2}}}\rrbracket \sin \chi =0. \end{aligned}$$

(58)

While Eqs. (55) hold on the whole of S, despite the jumps that the gradients $\nabla y_1$, $\nabla y_2$, and $\nabla y_3$ may suffer across the curves $C_j$ that $\varvec{y}$ transforms into the ridges ${\mathscr {C}}_j$, Eqs. (58) are valid only along such curves.

In the following section, we shall find solutions to (55) and (58) in a special case. It will be expedient to compute on a ridge the inner product $\varvec{\nu }_1\cdot \varvec{\nu }_2$, which features in the expression for $f_r$ in (47). To this end, we first recall (15) and remark that

$$\begin{aligned} \varvec{\nu }_2= & {} (\varvec{a}_1{+}{\llbracket {\varvec{a}}}\rrbracket ){\times }(\varvec{b}_1+{\llbracket {\varvec{b}}}\rrbracket )\nonumber \\= & {} \varvec{\nu }_1+\varvec{a}_1\times {\llbracket {\varvec{b}}}\rrbracket +{\llbracket {\varvec{a}}}\rrbracket \times \varvec{b}_1, \end{aligned}$$

(59)

where use has also been made of (57). Since $\varvec{a}_1=\lambda _1\varvec{n}_1$ and $\varvec{b}_1=\lambda _2\varvec{\nu }_1\times \varvec{n}_1$, by (54) it follows from (59) that

$$\begin{aligned} \varvec{\nu }_1\cdot \varvec{\nu }_2=\frac{1}{\lambda _1^2}\varvec{a}_1\cdot \varvec{a}_2+\lambda _1^2\varvec{b}_1\cdot \varvec{b}_2-1, \end{aligned}$$

(60)

where we also employed (23). With $\varvec{a}_i$ and $\varvec{b}_i$ given by (52), we easily revert (60) into an expression featuring the traces of the gradient components $y_{i,j}$ on the two sides of the ridge under consideration.

6 Ridged cones

It is time now to put our theory to the test. In this section, we shall consider a classical example, already treated within the traditional theory [24], that of a disk S of radius R upon which the planar radial hedgehog $\varvec{m}$ has been imprinted. In polar coordinates $(\varrho ,\vartheta )$, with associated orthonormal frame $(\varvec{e}_\varrho ,\varvec{e}_\vartheta )$, $\varvec{m}=\varvec{e}_\varrho$ and $\varvec{m}_\perp =\varvec{e}_\vartheta$. It is an easy exercise to check with the aid of (49) that then $c_1=0$, $c_2=\frac{1}{\varrho }$, and $c_{12}=-\frac{1}{\varrho ^2}$, so that by (56) $K=0$, independently of the prescribed principal stretches $\lambda _1$ and $\lambda _2$.

We shall use the representation (50) for $\varvec{y}$ with

$$\begin{aligned}&y_1=f\sin \psi \cos \varphi ,\quad y_2=f\sin \psi \sin \varphi ,\nonumber \\&y_3=f\cos \psi , \end{aligned}$$

(61)

where f, $\psi$, and $\varphi$ are assumed to be picewise $C^2$-functions of $(\varrho ,\vartheta )$. Thus $\psi$ and $\varphi$ represent the polar and azimuthal angles of $\varvec{y}$ in the movable frame $(\varvec{e}_\varrho ,\varvec{e}_\vartheta ,\varvec{e}_\varphi )$, while f measures radial dilation (or contraction). Some labor is required to see that with this choice the isometry conditions (55) become

$$\begin{aligned}&f_{,\varrho }^2+f^2(\psi _{,\varrho }^2+\varphi _{,\varrho }^2\sin ^2\psi )=\lambda _1^2, \end{aligned}$$

(62a)

$$\begin{aligned}&\frac{1}{\varrho ^2}\left[ f_{,\vartheta }^2+f^2\left( \psi _{,\vartheta }^2+(1+\varphi _{,\vartheta })^2\sin \psi ^2\right) \right] =\lambda _2^2, \end{aligned}$$

(62b)

$$\begin{aligned}&\frac{1}{\varrho }\left\{ f_{,\varrho }f_{,\vartheta }+f^2\left[ \psi _{,\varrho }\psi _{,\vartheta }+\varphi _{,\varrho }(1+\varphi _{,\vartheta })\sin ^2\psi \right] \right\} =0,\nonumber \\ \end{aligned}$$

(62c)

where commas denote partial derivatives in the variables $(\varrho ,\vartheta )$.

Fig. 6

Geometric construction of a ridged immersion of the unit disk S

Fig. 7

On the left: Functions $\psi$ (red) and $\varphi$ (blue) for $\mu =2.7$ and different values of n. On the right: Corresponding ridged immersions of the unit disk S

Fig. 8

When the geometric construction of a ridged immersion of S becomes singular. On the left: Functions $\psi$ (red) and $\varphi$ (blue) for $\mu =n=3$. On the right: Corresponding ridged immersion of S with all ridges collapsed on the vertical axis

We shall look for solutions of (62) under the simplifying assumption that f is a positive function of $\varrho$ only and both $\psi$ and $\varphi$ are functions of $\vartheta$ only. The deformed surface ${\mathscr {S}}$ is thus conical, in the same class employed to represent crumpled sheets of paper in [49, 50].¹⁵ Further requiring that the centre of the disk S is held fixed in a spontaneous deformation, we see that (62a) has the unique (positive) solution

$$\begin{aligned} f(\varrho )=\lambda _1\varrho , \end{aligned}$$

(63a)

and that (62c) is identically satisfied, while (62b) reduces to

$$\begin{aligned} \psi '^2+(1+\varphi ')^2\sin ^2\psi =\mu ^2, \end{aligned}$$

(63b)

where use has also been made of (63a) and we have set

$$\begin{aligned} \mu :=\frac{\lambda _2}{\lambda _1}=\frac{1}{\lambda _1^2}. \end{aligned}$$

(63c)

Hereafter in this section, a prime $'$ will denote differentiation with respect to $\vartheta$.¹⁶

We now focus on solving (63b) subject to the periodic boundary conditions

$$\begin{aligned} \psi (0)=\psi (2\pi ),\quad \varphi (0)=\varphi (2\pi ). \end{aligned}$$

(64)

For a given function $\psi$ satisfying the first equality in (64) and such that $\sin \psi \ne 0$ and $\psi '^2\leqq \mu ^2$, (63b) is solved by integrating

$$\begin{aligned} \varphi '=-1+\frac{\sqrt{\mu ^2-\psi '^2}}{\sin \psi }, \end{aligned}$$

(65)

provided that the following condition is met,

$$\begin{aligned} I:=\int _0^{2\pi }\frac{\sqrt{\mu ^2-\psi '^2}}{\sin \psi }\mathrm {d}\vartheta =2\pi , \end{aligned}$$

(66)

for the second equality in (64) to be valid too. Clearly, for $\mu \leqq 1$, (66) is solved by $\psi \equiv \arcsin \mu$, which corresponds to a circular cone when $\mu <1$.

It is for $\mu >1$ where our discrete ridge model comes to play. We look for ridged isometric immersions representable in the class (61). First, we see how to write the jump conditions (58) in the present context. Here $\chi =0$, as jumps may only occur along radii of the disk S. Moreover, by (63a), Eq. (61) implies that ${\llbracket {y_{1,1}}}\rrbracket ={\llbracket {y_{2,1}}}\rrbracket ={\llbracket {y_{3,1}}}\rrbracket =0$, independently of $\psi$, so that (58) is identically satisfied. We shall thus look for solutions $\psi$ of (63b) that are piecewise of class $C^1$. Since (63b) must be valid on the whole of S, a jump in $\psi '$ is admissible only if

$$\begin{aligned} {\llbracket {\psi '^2}}\rrbracket =0, \end{aligned}$$

(67)

while, by (65), no jump in $\varphi '$ is allowed.

It is a simple matter to show that for a conical surface described by (61) the vectors $\varvec{a}$ and $\varvec{b}$ in (52) can be given the following expressions (also by use of (63a)),

$$\begin{aligned} \varvec{a}= & {} \lambda _1\left( \sin \psi \cos \varphi \varvec{e}_\varrho +\sin \psi \sin \varphi \varvec{e}_\vartheta +\cos \psi \varvec{e}_3\right) , \end{aligned}$$

(68a)

$$\begin{aligned} \varvec{b}= & {} \lambda _1\{[\psi '\cos \psi \cos \varphi -(1+\varphi ')\sin \psi \sin \varphi ]\varvec{e}_\varrho \nonumber \\+ & {} [\psi '\cos \psi \sin \varphi \,{+}\,(1\,{+}\,\varphi ')\sin \psi \cos \varphi ]\varvec{e}_\vartheta \nonumber \\- & {} \psi '\sin \psi \varvec{e}_3\}. \end{aligned}$$

(68b)

It follows from these equations and (60) that on a ridge

$$\begin{aligned} \varvec{\nu }_1\cdot \varvec{\nu }_2=1-2\frac{\psi '^2}{\mu ^2}, \end{aligned}$$

(69)

where use has also been made of both (63b) and (67). Moreover, since $\varvec{t}_\perp \cdot \varvec{n}=0$, by (47) the ridged energy density (scaled to $\frac{8}{3}h^2$) simply reduces to $f_r=\arccos ^2(\varvec{\nu }_1\cdot \varvec{\nu }_2)$.

We must still enforce (64) and (66). We do so by means of a geometric construction, which is fully substantiated in “Appendix A”. We take full advantage of our unconventional approach that here renounces to approximate the smooth deformed surface ${\mathscr {S}}$ with a large (albeit finite) number of ridges. Granting (up–down) symmetry to ${\mathscr {S}}$, we shall be contented with capturing the simplest feature: how many (up ad down) folds it possesses. Thus, here the unknown number of folds will just be the number of admitted ridges. Minimizing the total ridge energy will provide the optimal number of folds.

By rescaling lengths, we can assume that S is the unit disk, stretched circumferentially by $\mu$ into a symmetric ridged immersion through the following steps (see Fig. 6). (1) For a given $\mu >1$, choose an integer $n\geqq \mu$ and take a circular sector of S of amplitude $\alpha :=\frac{\pi }{2n}$ (the reference sector). (2) Stretch it so that its amplitude becomes $\mu \alpha$. (3) Rotate the stretched sector around one of its edges by an appropriately chosen angle $\beta$. (4) Reflect the rotated sector across the vertical plane containing the lifted edge. (5) Rotate by $\pi$ the reflected sector around its edge lying on S. (6) Repeat $n-1$ times the preceding steps.

This construction works, making sure that the generated surface closes on itself, if $\beta$ is such that the reflection plane in step (4) cuts S where lied the edge of the reference sector (see Fig. 6d), that is, if (see “Appendix A”).

$$\begin{aligned} \beta := {\left\{ \begin{array}{ll} \arccos \frac{\tan \alpha }{\tan (\mu \alpha )} &{}\text {if}\quad \mu <n,\\ \frac{\pi }{2} &{}\text {if}\quad \mu =n. \end{array}\right. } \end{aligned}$$

(70)

A continuous, piecewise differentiable $4\alpha$-periodic function $\psi$ is obtained in (A5) by extending over ${\mathbb {R}}$ the function $\psi ^*$ defined on $[0,\alpha ]$ by

$$\begin{aligned} \psi ^*(\vartheta ):=\arccos (\sin \beta \sin (\mu \vartheta )). \end{aligned}$$

(71)

Correspondingly, the function $\varphi ^*$ associated with $\psi ^*$ is given in $[0,\alpha ]$ by (see again “Appendix A”)

$$\begin{aligned} \varphi ^*(\vartheta ):=-\vartheta +\arccos \frac{\cos (\mu \vartheta )}{\sin \psi ^*(\vartheta )}. \end{aligned}$$

(72)

Its $C^1$-extension $\varphi$ over ${\mathbb {R}}$ compatible with (65) is recorded in (A6). It is a boring, but simple exercise to check that functions $\psi ^*$ and $\varphi ^*$ satisfy (63b) identically, and so do their extensions $\psi$ and $\varphi$.

Figure 7 shows examples of the functions $\psi$ and $\varphi$ so generated alongside with the conical surface produced by the corresponding ridged immersions of the unit disk S.

Fig. 9

The function $F_r$ in (73) is plotted against $\mu$ in the intervals [1, n], for several values of n (blue graphs). The lower envelope (red graph) is the plot of $F_r(\lceil \mu \rceil ,\mu )$

It should be noted that whenever $\mu =n$ this construction becomes singular, though it is still applicable. In such a case, all ridges are degenerate and lie on the vertical axis $\varvec{e}_3$ of the disk S; all faces of ${\mathscr {S}}$ are vertical as well and $\varphi$ becomes discontinuous (as shown in Fig. 8 for $\mu =3$).

The above construction produces $N=2n$ ridges, at each of which $\psi '^2=\psi '^2(\alpha )$, so that by (69) the total (dimensionless) ridge energy $F_r$ is (see (48) and (A7))

$$\begin{aligned} F_r(n,\mu )=2n\arccos ^2\left( 1{-}\frac{2}{\sin ^2\left( \frac{\mu \pi }{2n}\right) }\left[ \cos ^2\left( \frac{\pi }{2n}\right) {-}\cos ^2\left( \frac{\mu \pi }{2n}\right) \right] \right) .\nonumber \\ \end{aligned}$$

(73)

Plots of $F_r$ against $\mu$ for several values of n are depicted in Fig. 9. They show that, for a given $\mu$, $F_r$ is minimized for $n=\lceil \mu \rceil$, when the number of ridges is the least possible (as was perhaps to be expected).¹⁷

7 Conclusions

Common wisdom has it that in sufficiently thin sheets of nematic polymer networks, as in all elastic materials for that matter (see, for example, [51, p. 396] or [52, p. 404]), the bending energy (which scales as the cube of the thickness) may be neglected relative to the stretching energy (which scales linearly in the thickness). When activated, a nematic polymer network suffers a spontaneous deformation that attempts to transfer on the current shape the metric tensor that minimizes the stretching energy, which (with a slight abuse of language) we called an isometric immersion, for short. Such an immersion would generally depend on the nematic director imprinted on the sheet at the time of crosslinking, and, as is well known, it may fail to exist.

We started from relaxing the requirement of smoothness for an isometric immersion, thus removing a possible obstacle to its existence. We allowed for ridges in the immersed surfaces ${\mathscr {S}}$ representing deformed sheets; these are lines where the normal to ${\mathscr {S}}$ suffers a jump. Clearly, ridges do not come for free (nothing does). If they did, we would be overwhelmed with a superabundance of shapes, for which we would lack a selecting energy criterion (as all would have the same stretching energy).

We thought of ridges as concentrations of bending energy; we put forward a model to compute the energy they bear distributed along their length. To accomplish this task, we employed a formula for the bending energy recently derived from the “trace formula” valid in three-space dimensions [21].

We showed that the ridge energy density (per unit length) $f_r$ scales quadratically with the sheet’s thickness, and so it represents a contribution intermediate between stretching and bending energies. The formula we obtained for $f_r$ not only depends (symmetrically) on the normals to the adjoining sides (as was perhaps to be expected), but also on the orientation of the nematic director relative to the tangent to the ridge.

We applied our theory to the case where a planar hedgehog is imprinted on a flat disk at the time of crosslinking. We studied the total elastic energy, including the new ridge energy, in a class of conical deformations not new in the literature. In the regime where the radii of the reference disk shrink and the circumferences expand, by identifying the ridges of the discrete model with the folds of the continuum model, we used the total ridge energy as a selection criterion to determine the optimal number of folds.

The ridged cones that we found as energy minimizers are neither the developable cones of [48, 49, 53,54,55,56] nor the excess cones of [50], as they do not share the degree of smoothness that both the latter and the former have in common. But we trust that, under similar circumstances, they have in common the same number of folds.

Of course, we cannot expect that a rubber-like material will spontaneously take on sharp ridges when activated by a change in its internal material organization. Our model, as applied here, has more the flavor of a vicarious theory, where a distributed bending energy is replaced by one concentrated at a number of places. As we have shown, it is, nevertheless, predictive, at least of the expected number of folds. Since $F_r$ in (73) is minimized for $n=\lceil \mu \rceil$ and, correspondingly, the optimal number of ridges/folds is $N=2\lceil \mu \rceil$, it follows from (63c) and (23) that

$$\begin{aligned} N=2\left\lceil \sqrt{\frac{s_0+1}{s+1}}\right\rceil \geqq 4\quad \text {for}\quad s<s_0, \end{aligned}$$

(74)

where s is the activation parameter of the theory.

We expect that (74) would reproduce the number of folds predicted by an elastic theory based on the full blown energy, where stretching and bending components are blended together and compete on different length scales.

Symmetry and identification of folds with ridges played a role in deriving (74). The ridge energy obtained in this paper, like the notion itself of ridged immersions, is susceptible of further applications, if we relax the ridge/fold identification and take the more traditional approach of considering our (simplified) discrete model as an approximation of the (more difficult) continuum model (with bending energy). This approximation, which is expected to improve upon increasing the number of ridges, should be justified by a convergence assessment. For a finite (but large) number of ridges unrestricted by symmetry requirements, it would provide a good test for the number of folds predicted by (74). More generally, an appropriate decomposition of the reference surface S could be ridge-immersed in a triangulation of the deformed surface ${\mathscr {S}}$ to determine the optimal shape of the activated film. These extensions are presently being studied.

Appendix A: Geometric construction

In this appendix, we provide further details on the geometric construction employed in Sect. 6 to produce a ridged isometric immersion of the unit disk S. In particular, our objective will be to justify both the formula for $\beta$ in (70) and that for $\psi ^*$ in (71).

Take a circular sector of S with amplitude $\alpha$ (the reference sector) and stretch it uniformly to obtain a sector with amplitude $\mu \alpha$. Then rotate the stretched sector around $\varvec{e}_{1}$ by $\beta$ and reflect the rotated sector across the vertical plane passing through the unit vector $\varvec{e}(\alpha ):=\cos \alpha \varvec{e}_1+\sin \alpha \varvec{e}_2$. Finally, rotate by $\pi$ the reflected sector around $\varvec{e}(2\alpha )$. Thus, we have ridged-immersed a sector with amplitude $4\alpha$. By replicating $n-1$ times this immersion, we finally obtain the entire surface ${\mathscr {S}}$. Figure 10 shows the geometric details of the first immersed sector (with amplitude $\alpha$).

Fig. 10

Geometric details for the first immersed sector

More precisely, for all $\vartheta \in [0,\alpha ]$ we define the mapping

$$\begin{aligned} \varvec{e}(\vartheta )\mapsto \varvec{E}(\vartheta ):={\mathbf {R}}_{\varvec{e}_{1},\beta }\varvec{e}(\mu \vartheta ), \end{aligned}$$

(A1)

where ${\mathbf {R}}_{\varvec{e}_{1},\beta }$ is the rotation around $\varvec{e}_{1}$ by the angle $\beta$. Explicitly, $\varvec{E}(\vartheta )$ is given by

$$\begin{aligned} \varvec{E}(\vartheta )= & {} \left( {\mathbf {I}}+\sin \beta {\mathbf {W}}(\varvec{e}_1)-(1-\cos \beta ){\mathbf {P}}(\varvec{e}_1)\right) \varvec{e}(\vartheta )\nonumber \\= & {} \cos (\mu \vartheta )\varvec{e}_{1}+\sin (\mu \vartheta )\cos \beta \varvec{e}_{2}\nonumber \\&+\sin (\mu \vartheta )\sin \beta \varvec{e}_{3}, \end{aligned}$$

(A2)

where ${\mathbf {P}}(\varvec{e}_1)$ denotes the projection onto the plane orthogonal to $\varvec{e}_1$ and ${\mathbf {W}}(\varvec{e}_1)$ is the skew-symmetric tensor associated with $\varvec{e}_1$. It follows from (61) and (A2) that for all $\vartheta \in [0,\alpha ]$ the polar angle of the first immersed sector is given by (71) in the main text.

The appropriate value of $\beta$ is determined by requiring that the plane spanned by $\varvec{e}(\alpha )$ and $\varvec{E}(\alpha )$ is the vertical plane used to reflect the first immersed sector in step 4 of our construction (see Fig. 6d). Formally, this geometric requirement demands that ${\mathbf {P}}(\varvec{e}_3)\varvec{E}(\alpha ) = \delta \varvec{e}(\alpha )$, for some $\delta \in [0,1]$, where ${\mathbf {P}}(\varvec{e}_{3})$ is the projection onto the plane orthogonal to $\varvec{e}_{3}$, that is,

$$\begin{aligned}&\cos (\mu \alpha )\varvec{e}_{1}{+}\sin (\mu \alpha )\cos \beta \varvec{e}_{2}\nonumber \\&\quad = {\left\{ \begin{array}{ll} \delta (\cos \alpha \varvec{e}_{1}{+}\sin \alpha \varvec{e}_2) &{}\text {if}\quad \mu \alpha <\frac{\pi }{2},\\ \varvec{0}&{}\text {if}\quad \mu \alpha =\frac{\pi }{2}. \end{array}\right. } \end{aligned}$$

(A3)

Solving (A3), we easily find that $\beta$ is given by (70), while

$$\begin{aligned} \delta =\frac{\cos \mu \alpha }{\cos \alpha }. \end{aligned}$$

(A4)

Since, by (61), $\varvec{E}(\vartheta )\cdot \varvec{e}_1=\sin \psi \cos (\varphi +\vartheta )$ (see Fig. 10c), from (A2) we also retrieve (72).

Finally, we note that the function $\psi ^*$ can be extended continuously to the entire ${\mathbb {R}}$ by making it periodic with period $4\alpha$ and oscillating symmetrically about $\frac{\pi }{2}$. The extended function $\psi$ is represented as follows in one period,

$$\begin{aligned} \psi (\vartheta ) = \left\{ \begin{array}{ll} \psi ^*(\vartheta ) &{}\quad \text {if}\quad \vartheta \in [0,\alpha ],\\ \psi ^*(2\alpha -\vartheta ) &{}\quad \text {if}\quad \vartheta \in [\alpha ,2\alpha ],\\ \pi -\psi ^*(\vartheta -2\alpha )&{}\quad \text {if}\quad \vartheta \in [2\alpha ,3\alpha ],\\ \pi -\psi ^*(4\alpha -\vartheta )&{}\quad \text {if}\quad \vartheta \in [3\alpha ,4\alpha ]. \end{array}\right. \end{aligned}$$

(A5)

Similarly, the periodic $C^1$-extension of $\varphi ^*$ compatible with (65) is given by

$$\begin{aligned} \varphi (\vartheta ) = \left\{ \begin{array}{ll} \varphi ^*(\vartheta )&{}\quad \text {if}\quad \vartheta \in [0,\alpha ],\\ -\varphi ^*(2\alpha -\vartheta ),&{}\quad \text {if}\quad \vartheta \in [\alpha ,2\alpha ],\\ \varphi ^*(\vartheta -2\alpha )&{}\quad \text {if}\quad \vartheta \in [2\alpha ,3\alpha ],\\ -\varphi ^*(4\alpha -\vartheta ),&{}\quad \text {if}\quad \vartheta \in [3\alpha ,4\alpha ]. \end{array}\right. \end{aligned}$$

(A6)

Since both $\sin \psi$ and $\psi '^2$ are $2\alpha$-periodic functions, by (65) so is also $\varphi$. Moreover, it is easy to check that

$$\begin{aligned}&\psi '^2((2k+1)\alpha )=\frac{\mu ^2}{\sin ^2\alpha }(\cos ^2\alpha -\cos ^2(\mu \alpha ))\nonumber \\&\text {and}\quad \varphi (k\alpha )=0\quad \text {for all}\ k\in {\mathbb {N}}, \end{aligned}$$

(A7)

as illustrated in Fig. 7.