Journal of the Mechanics and Physics of Solids Forceless folding of thin annular strips

Thin strips or sheets with in-plane curvature have a natural tendency to adopt highly symmetric shapes when forced into closed structures and to spontaneously fold into compact multicovered configurations under feed-in of more length or change of intrinsic curvature. This disposition is exploited in nature as well as in the design of everyday items such as foldable containers. We formulate boundary-value problems (for an ODE) for symmetric equilibrium solutions of unstretchable circular annular strips and present sequences of numerical solutions that mimic different folding modes. Because of the high-order symmetry, closed solutions cannot have an internal force, i.e., the strips are forceless. We consider both wide and narrow (strictly zero-width) strips. Narrow strips cannot have inflections, but wide strips can be either inflectional or non-inflectional. Inflectional solutions are found to feature stress localisations, with divergent strain energy density, on the edge of the strip at inflections of the surface. ‘Regular’ folding gives these singularities on the inside of the annulus, while ‘inverted’ folding gives them predominantly on the outside of the annulus. No new inflections are created in the folding process as more length is inserted. We end with a discussion of an intriguing apparent connection with a deep result on the topology of curves on surfaces.


Introduction
Spontaneous folding of two-dimensional structures into complex three-dimensional shapes is found in nature (e.g., in flowers, leaves and insect wings (Faber et al., 2018)) and is exploited in the design of deployable structures used in space (Miura, 1993) or medical (Kuribayashi et al., 2006) applications. Similarly, the age-old art of paper folding has inspired the design of origami mechanisms (Hanna et al., 2014;Saito et al., 2016) and self-folding electronics (Miyashita et al., 2014;Sundaram et al., 2017), while controlled folding is of interest from a robotics perspective (Balkcom and Mason, 2008). Spontaneous assembly is also exploited in some everyday objects such as foldable baskets, so-called '2-second' pop-up tents and other light-weight foldable trekking items (Mouthuy et al., 2012).
During the folding and manipulation of thin sheets stress localisations may occur (Witten, 2007). When designing foldable or deployable structures it may be important to be able to predict how and where these occur. For instance, twisting a thin rectangular strip causes a repetitive pattern of alternating points of high stress on the two long edges of the strip (Korte et al., 2011). The topology of the classical Möbius strip gives an example of a surface with only one such point of stress localisation van der Heijden, 2007, 2015). In both cases the high stresses are associated with inflection points of the centreline of the strip, i.e., points of vanishing centreline curvature.
Thin sheets under most conditions deform predominantly by bending, without stretching, and such sheets have therefore been studied within the theory of inextensible plates. Absence of stretching means that an intrinsically flat strip within this theory is described by a developable surface. This geometrical structure can be used to reduce the strain energy of the strip to its centreline, giving the one-dimensional Wunderlich functional (Wunderlich, 1962). The equilibrium (Euler-Lagrange) equations for this functional are thus ordinary differential equations (ODEs)  and therefore much easier to analyse than the usual partial differential equations (PDEs) of plate theory.
In the reduced ODE formulation the stress localisations are described by logarithmic singularities of the strain energy density. These lead to singularities of the equilibrium equations that are therefore generally hard to solve. Narrow limits have therefore been considered (Sadowsky, 1931), but these are not accurate near inflection points.
Numerical solutions for the full Wunderlich model are given in Starostin and van der Heijden (2015) and Korte et al. (2011) for closed and open strips, respectively. The developable Möbius strip solution in Starostin and van der Heijden (2015) has recently been validated experimentally (Kumar et al., 2021). The equations have been extended to annular strips (Dias and Audoly, 2015) and developable shapes of diametrically folded circular discs with a hole are computed in Yu et al. (2022) and Yu (2022). Here we compute equilibrium shapes of annular strips and investigate their spontaneous folding and stress localisation.
Many of the foldable structures mentioned above consist of naturally curved strips. One folding mechanism is illustrated by cutting a ring along its radius and inserting extra material. This leads to frustration, i.e., inconsistency between the length of the new ring and its intrinsic shape. The structure, although developable, can no longer adopt a planar shape and is forced to deform out of the plane into a wavy saddle shape (as in Fig. 1), which may occur in various modes. The same effect is observed if instead of increasing the length the intrinsic curvature is increased. Such saddle shapes have for instance been created in annular bilayer nanostructures by means of heat-induced differential stress (Cho et al., 2010) and also in temperature-responsive hydrogel ribbons (Bae et al., 2014). We call such shapes overcurved. By contrast, if material is cut out of a ring, or its curvature is decreased, then the strip may be pulled into a conical shape (see Fig. 1). We call this shape undercurved. As we shall see later, the cone is not the only possible equilibrium shape in this case.
Shapes as in Fig. 1 have high-order symmetry. Because of this symmetry these closed structures (just as closed rods of similar shapes) cannot have an internal force; they are forceless. Knowing in advance that the force is zero allows us to significantly simplify the equilibrium equations and to make analytical progress that does not seem possible for general developable shapes. This progress is even stronger in the narrow limit, which is governed by a strain energy functional analogous to the Sadowsky functional (Sadowsky, 1931) for rectangular strips. Forceless (non-inflectional) solutions for this functional can be obtained explicitly (as is also the case for the Sadowsky functional Starostin and van der Heijden, 2018).
Thus here we consider forceless equilibria of circular annular strips closed without a twist. We consider both non-inflectional and inflectional shapes but are particularly interested in the latter as they have stress localisations at inflectional generators of the surface. Where the centreline crosses such a generator, it must have vanishing normal curvature. Thus we extend the work in Starostin and van der Heijden (2015) on developable surfaces with straight centreline in their planar development to the non-geodesic case of developable surfaces with non-straight centreline in their planar development.
The paper is structured as follows. Section 2 presents the geometry of developable surfaces for deformed annular strips. The mean curvature is computed, which is required for the (bending) strain energy. In Section 3 the Wunderlich-like reduction of this strain energy is performed. This leads to a second-order variational problem on a space curve for which the equilibrium equations are written down. Our treatment in these two sections is self-contained because our subsequent analysis requires details (e.g., the Hamiltonian) not touched on in Dias and Audoly (2015). The derivation of the equations (initially a differential-algebraic system of equations that we however immediately turn into a system of ordinary differential equations) is different from that in Dias and Audoly (2015), following more closely our treatment of rectangular strips (Starostin and van der Heijden, 2015). We also identify simple analytical solutions (circular cones) of the equations. Then in Section 4 we formulate boundary-value problems for (two-sided) strip configurations with high-order symmetry and present numerical solutions (all forceless), both inflectional and non-inflectional, the latter showing stress localisation. In Section 5 we consider the special Sadowsky-like case of narrow strips for which a reduction to a planar system is possible. Section 6 ends the study with a discussion of our results.

Basic geometrical description
We consider an inextensible plate of uniform thickness that in its unstressed state has the form of a plane strip bounded by two concentric circular arcs. By the developability property, the strip, however deformed, can be reconstructed from an arbitrary reference curve on the surface of the strip. We take for this reference curve the natural and convenient choice of the centreline of the strip, i.e., the curve that is equidistant from both circular arcs in the intrinsic geometry of the strip. We denote this curve by ( ) ∈ R 3 , where ∈ [0, ] is arclength along the curve and is its length. We assume that ( ) is a regular curve of differentiability class 3 . The strip may be arbitrarily long so that when developed, it may have multiple overlappings, or it may be too short to complete a single full turn (Fig. 1).
Developable surfaces are special cases of ruled surfaces. We describe the deformed strip therefore by the general parametrisation of a ruled surface: where = ′ is the unit tangent to the centreline and is a transverse unit vector (perpendicular to ). Here and in what follows a prime denotes differentiation with respect to . The straight lines = const are the generators of the surface. They make an angle  = arctan(1∕ ) with the positive tangent direction of the curve ( ) (see Fig. 2). The interval bounds − and + are the coordinates (along − ) of the points on the long edges of the strip where the generator intersects in the planar development of the strip. They are therefore functions of . They will be defined in Section 3.1. At this point is an arbitrary function; we will later derive a condition on such that (1) describes a developable surface. We require that the short edges = 0 and = of the strip are generators so that closure of the strip is equivalent to the condition (0) = ( ) (we only consider two-sided surfaces). We need to avoid the case = ∞, which corresponds to = 0, i.e., a generator that is tangent to the centreline. Deformations in which such a tangency happens can be studied but would need a different reference curve to be chosen.
Let ∶= × be the surface unit normal. The three vectors then define the orthonormal Darboux frame { ( ), ( ), ( )} at the centreline (see Fig. 2). After choosing a coordinate system we may identify orientations of the Darboux frame with elements of the group of orthogonal 3 × 3 matrices: This defines a skew-symmetric 3 × 3 matrix in the Lie algebra so(3) as follows: where we have introduced the 'hat' isomorphism between skew-symmetric matriceŝ= (3) and axial (or rotation) vectors = ( 1 , 2 , 3 ) ⊺ in R 3 . 1 By definition of the Darboux frame, we have 1 = , 2 = , 3 = , which are the geodesic torsion, the geodesic curvature and the normal curvature, respectively.
We compute the first fundamental form of the surface coordinate patch given by Eq. (1): to the surface is defined as ( , ) = × ‖ × ‖ (for = 0 it becomes the above ( )). We also need the second fundamental form defined by We can now compute the shape operator (or Weingarten map) , i.e., the linear operator on the tangent plane defined by ( ) = − , where is a unit tangent vector to the surface (Spivak, 1999). ( ) is the gradient of the unit normal to the surface in the tangent direction . It therefore encodes information about the curvature of the surface. In fact, the eigenvalues of the shape operator at each point are the principal curvatures at this point and the eigenvectors are the principal directions. We have ( ) = − , ( ) = − , and hence in the basis ( , ), [1 − ( ′ + (1 + 2 ))] 2 + 2 (1 + 2 )( − ) 2 , We then compute the Gaussian curvature where 1 and 2 are the two principal curvatures of the surface. Now we set = ∕ . It is easy to check that this choice makes the Gaussian curvature vanish. Thus the ruled surface described by Eq. (1) is actually developable. The developability condition is therefore the same as for the geodesic case of a straight centreline (Starostin and van der Heijden, 2015), but with the curvature replaced by the normal curvature and the torsion replaced by the geodesic torsion. For a developable surface the surface normal is constant along generators, i.e., ( , ) = ( ). The area element becomes while the coefficients of the shape operator reduce to For the mean curvature we then find Because the surface described by Eq. (1) is henceforth taken to be developable, the deformed configuration of the strip can be obtained from a planar reference configuration by bending only, without stretching, provided both reference and deformed configurations have bounding end generators with the same angle (or, with the same (0) and ( )) (for closed configurations this proviso is superfluous). The deformation implicit in Eq. (1) is therefore isometric. We verify this property in Appendix A by first constructing the reference configuration and then calculating the relevant strain tensor.

Closed strips
For a closed annular strip of finite radius = −1 and centreline length we define the relative length excess parameter ∶= ∕(2 ) − 1, > −1. If < 0 then we say that the closed strip is undercurved (relative to the single-covered annulus of length 2 ∕ ) and if > 0 then the strip is called overcurved. The fact that ∈ Z ≥0 does not imply that the closed strip is planar, but all closed planar solutions have ∈ Z ≥0 (of these solutions only the one with = 0 is non-self-intersecting).

Edge of regression
The asymptotic completion of a developable strip is defined as the surface obtained by extending all generators to infinity in both directions, i.e., taking in all of R in Eq. (1). From Eq. (4) we see that the mean curvature becomes singular if the parameter in Eq. (1) equals 1∕[ ′ + (1+ 2 )]. If ′ = − (1+ 2 ) then there is no singularity on the extended generator in the asymptotic completion. We call points on the centreline where ′ = − (1 + 2 ), 'cylindrical'. At such points the mean curvature is independent of , so the principal curvatures are constant along the local generator. The meaning of the cylindrical point condition becomes transparent when we rewrite it in terms of the angle between the generator and the tangent vector: ′ = . At a cylindrical point the generators must remain parallel. As the tangent vector turns with the angular rate , the generator should rotate with respect to the tangent at the same rate.
Away from cylindrical points we can define the curve which is called the edge of regression. When [ ′ ( ) + (1 + 2 ( ))] changes sign the edge of regression jumps within the asymptotic completion from one side of the centreline to the other. The strip cannot be wider than the critical value of , i.e., we require that ∈ [ − , + ] and either 1∕[ ′ + (1 + 2 )] ≤ − or + ≤ 1∕[ ′ + (1 + 2 )]. The developable surface is then the envelope of tangents to the edge of regression (i.e., of the generators of the surface) and is therefore also the tangent developable. This envelope meets the edge of regression in two sheets that form cusps in the normal plane to (Naokawa, 2013). The edge of regression may have its own singularities. By differentiation we find

At points where
the tangent to the edge of regression is discontinuous. This corresponds to a cusp point on the curve and a swallow tail singularity of the asymptotic completion of the strip. We call isolated points of the centreline where Eq. (5) holds 'conical'.

Energy functional
We assume that the annular plate is planar in its relaxed state and remains developable throughout deformation. Its width is 2 = const and < 1, = const > 0. The length of the strip is arbitrary and may exceed 2 ∕ . Technologically, if it exceeds this length then the strip cannot be cut out of a sheet of film or paper, so other means of obtaining such models must be used. A kind of growing material may for instance be considered. For a sufficiently thin plate the stretching energy may be neglected (Cerda and Mahadevan, 2005) and we are left with the bending energy, which is where = 2 ℎ 3 ∕[3(1 − 2 )] is the flexural rigidity, with 2ℎ the thickness of the plate, Poisson's ratio and Young's modulus. By substituting the mean curvature from Eq. (4), the area element from Eq. (3) and performing integration first in the radial direction and then along the centreline we get To find the limits − and + , we solve the equation for the edges of the strip in its planar development (see Fig. 2): ( − 1∕ ) 2 + 2 = (1∕ ± ) 2 , where and are the coordinates of the orthogonal reference frame with origin at the current point of the centreline, the coordinate along the current tangent direction and the coordinate in the perpendicular direction, positive in the direction of the centroid of the annulus; so the generator is described by = − . There exists a critical value * = 1− √ (2− ) such that for | | > * the generator never crosses the inner edge. Restricting ourselves first to the regular case by assuming that | | ≤ * , we find the limits Note that in the conical limit → 1, * = 0. Generally, when = 0, + = , − = − . Now consider briefly the case when there exists a generator that is tangent to the inner edge. Such a generator crosses the centreline twice with = ± * at points 0 and 1 , say. Assuming | ( )| > * for ∈ [ 0 , 1 ] we conclude that there must be a point of discontinuity ( ) = ∞ within this interval where the generator becomes tangent to the centreline (and where the width of the strip must vanish). We have excluded such cases from our consideration by requiring finite everywhere, i.e., the generators must always cross the centreline transversely. The other possibility is that | ( )| < * for ∈ [ 0 , 1 ] and the generators cross the inner edge all in the same point (a singular point of the edge of regression). In this case the single generator tangent to the inner edge can simply be viewed as two (aligned) generators: one at point 0 with ( 0 ) = * and the other at 1 with ( 1 ) = − * (or vice versa). The segment of the strip between 0 and 1 must remain planar.
The -integration in Eq. (7) can be carried out analytically (as in the case of straight parallel edges Wunderlich, 1962; Starostin and van der Heijden, 2015), and we arrive at The argument of the logarithm is 0 or ∞ when generators intersect on the outer or inner edge, respectively, corresponding to divergence of the mean curvature , and hence the bending energy density, on the edge of the strip. Note that equilibrium shapes do not depend on the material properties of the strip: is a simple factor. Also note that in the limit of a narrow strip, → 0, we have → 1 and no derivative enters the integrand in Eq. (8). Thus, in this case we obtain a functional analogous to Sadowsky's (1930Sadowsky's ( , 1931, in which the curvature and torsion of the centreline are replaced by the normal curvature and the geodesic torsion, respectively. The other limiting case is obtained by pushing the geodesic curvature to zero, → 0, i.e., the radius of the annular strip tends to infinity and the strip, when developed onto the plane, becomes a band with straight edges. Then → and → and the bending energy reduces to Eq. (8) with ℎ( , , ′ ) = 2 ( 1 + 2 ) 2 ( ′ ), which is the Wunderlich functional (Wunderlich, 1962;van der Heijden, 2007, 2015). For solutions of closed strips or strips with fixed end points we impose a constant end-to-end distance constraint by adding the following integral expression to the bending energy: where we have used inextensibility of the centreline, i.e., = ′ , and is a (constant) Lagrange multiplier (with the physical meaning of an internal force). Finally, we have to impose the constraint that the centreline is circular in the intrinsic geometry, i.e. its geodesic curvature is a constant: 2 = = const. We enforce this constraint by adding the integral = ∫ 0 2 2 d , where 2 = 2 ( ) is another (local) Lagrange multiplier.
In conclusion, equilibrium shapes of a thin, inextensible and intrinsically planar strip are given by stationary points of the functional where  is the Lagrangian of the problem and = ∕( ) (all force and moments in the rest of the paper are thus normalised by ). This represents a 1D variational problem on a curve in R 3 cast in Euclidean invariant form. Although we are here interested in closed-strip solutions, the functional  can also be used to study open-strip solutions. If for such open strips the ends are free to move relative to each other then represents the work done by any applied end force . However, since we use the parametrisation of Eq. (1) in the -integration to obtain the one-dimensional integral (8), the short edges of the strip must be generators and therefore straight. This is quite natural for fixed ends (held by a straight clamp covering the entire edge of the strip) but constitutes a restriction on deformations considered for strips with free ends.

Equilibrium equations
The derivation of the equilibrium equations follows closely the approach in Starostin and van der Heijden (2015) to which we refer for details. Alternative methods for deriving the equations can be found in van der Heijden (2009), Hornung (2010) and Dias and Audoly (2015).
In terms of the original rotation matrix the Lagrangian  may be expressed as  =  ( , ′ , ′′ ), where is considered a parameter. This form shows that the variational problem is second-order, i.e., defined on the second tangent bundle of the symmetry group SO(3). This is a consequence of using developability to reduce the variational problem for a (inextensible) two-dimensional elastic body (a strip) to one for a one-dimensional elastic body (the strip's centreline), giving equilibrium equations in the form of ODEs rather than PDEs. We can then apply Euler-Poincaré reduction to the tangent space (in the presence of advected parameters) and obtain the symmetry-reduced functional ( , ′ , ) =  ( , ′ , ′′ ), where is the force in the body frame (Gay-Balmaz et al., 2012). In fact, this is the natural form of the functional  already derived in (11). We can therefore immediately write down the (Euler-Lagrange) equilibrium equations in the form of Euler-Poincaré equations. They consist of Starostin and van der Heijden (2015): (a) balance equations for the components of the internal force = ( 1 , 2 , 3 ) ⊺ and moment = ( 1 , 2 , 3 ) ⊺ expressed in the Darboux frame (Starostin and van der Heijden, 2009;Hornung, 2010) and (b) the 'constitutive' equations It follows directly from Eqs. (12) and (13) that | | 2 and ⋅ are first integrals, i.e., are independent of .

Hamiltonian
Legendre transformation of the second-order reduced Lagrangian gives the symmetry-reduced Hamiltonian For a uniform strip with ℎ not explicitly depending on arclength ,  is a conserved quantity, as can be verified directly by differentiating the right-hand side of Eq. (19) with respect to and using Eqs. (12), (13), (15) and (16) to show that  ′ = 0. This conserved quantity will be useful in reducing the equilibrium equations in various special cases (Sections 3.8, 3.9 and 5.1). It is also useful for measuring the accuracy of a given numerical solution.
With the help of Eq. (15) and the obvious property of ℎ: ℎ = 2ℎ, we can rewrite Eq. (19) as and after substitution of ℎ from Eq. (10) in the explicit form In the limit → 0 the above simplifies to (cf. Starostin and van der Heijden (2015)). The first term in Eq. (21) may be broken up into two so that we may write it as Intriguingly, each of the two new terms may be viewed as the bending energy density (in units of ) at one of the edges of the strip (compare Eq. (7)). For narrow annular strips,

Symmetries
The equilibrium equations (Eqs. (12), (13), (17), (18)) are invariant under the non-reversing involution and the reversing involution Note that in both cases remains the same (does not change sign). The following non-reversing involution in which changes sign will also be useful:

Kinematics equations
Reconstruction of the centreline of the strip requires solving for the tangent and integrating this to get . We choose a parametrisation of the Darboux frame { , , } in terms of three Euler angles , and (Love, 1927) (not to be confused with the angles defined in Fig. 2 − sin sin + cos cos cos cos sin + cos sin cos − sin cos − sin cos − cos cos sin cos cos − cos sin sin sin sin The Euler angles are related to the Darboux vector by the kinematics equations ′ = ( sin − cos ) csc , Note that the angle should stay away from 0 and . To guarantee this we must choose the axis of the laboratory reference frame such that the tangent never aligns with the ± axis (alternatively, we can use a parametrisation of the Darboux frame in terms of Euler parameters Starostin and van der Heijden, 2014).

Special solution (circular cone)
We begin the exploration of equilibrium solutions by testing our equations for the existence of the conical solution in Fig. 1. For a conical solution all centreline points are 'conical', i.e., Eq. (5) is satisfied for all . It immediately follows that the only uniform solution = const. is ≡ 0, i.e., the generator is everywhere perpendicular to the centreline.
For the analysis of this case we shall need The first constitutive equation Eq. (15) gives us 3 = 2 0 , while the second Eq. (16) implies 1 = −2 ′ 1 . From the moment balance equations (the 4th to 6th Eq. (26)) we find 2 = 2 1 (we assume non-vanishing ), while for 1 we use the expression of the Hamiltonian:  0 = 2 0 + 2 + 1 . With these , = 1, 2, 3, the first force balance equation in Eq. (26) is identically satisfied, while the second and third give, respectively, The first equation is a reminder that conical deformations of elastic sheets under various conditions are governed by the elastica equation in which the ordinary curvature is replaced by the normal curvature (Cerda and Mahadevan, 2005; Starostin and van der Heijden, 2015).
The above system immediately implies that the only solution is = const. The centreline is thus a circle of radius 1∕ √ 2 + 2 and the length of the closed single-covered strip is = 2 ∕ √ 2 + 2 . Therefore, the uniform conical solution exists only for an undercurved annular strip. The surface of the strip belongs to a circular cone. All the generators intersect in a single point, the apex of the cone, and make a constant right angle with the centreline (Fig. 3).
For constant the forces and moments reduce to 1 = 2 = 3 = 0 (i.e., the solution is forceless) and 1 = 0, 3 = 2 (which means that the constant moment vector is directed along the axis of the cone). Explicitly, 2 = 1 log 1+ 1− and 3 = log 1+ 1− . In the narrow limit ( → 0) this becomes 2 = 2 , 3 = 2 . Let the centreline lie in the coordinate plane; to have this we set ≡ ∕2. Then the kinematics equations Eq. (24) imply = − arctan( ∕ ) = const and ′ = ± √ 2 + 2 . The latter can be integrated to find the angle as a linear function of arclength: = ± √ 2 + 2 + 0 . Further integration of Eq. (25) provides explicit expressions for the coordinates of the circular centreline. We finally comment that, as we shall see in the following, conical solutions are not the only equilibria of undercurved closed annular strips.

Perturbation of conical solutions and linearisation
We first expand Eqs. (15), (16) to express the moments 1 and 3 as Our aim here is to study non-inflectional shapes that appear as small perturbations of multicovered circular conical solutions. Thus we assume | ( )| ≪ 1, | ′ ( )| ≪ 1. Then we can write = 0 + 1 ′ + ( 2 , ′2 ) with 0 and 1 defined in Eqs. (27), (28), respectively. For the derivatives, we have We only consider forceless solutions. From Eq. (20) we then obtain  − 2 = 2 ( 0 + 2 1 ′ ) + ( 2 , ′2 ), from which we find the remaining second component of the moment Now let the normal curvature be the sum of a constant and a small variable term: ( ) = 0 + 1 ( ), 0 = const, | 1 ( )| ≪ 1. Assuming that 0 ≠ 0, we introduce the normalised variable term ( ) ∶= 1 ( )∕ 0 . Keeping only constant and linear terms in , 1 and their derivatives we then have We substitute the above into the first moment balance equation Eq. (13) and equate the zero-and first-order terms. The first results in an expression for the value of the Hamiltonian  = (2 2 + 2 0 ) 0 , while the second, together with the other two equations Eq. (13), gives us three linear homogeneous differential equations for ( ) and ( ), only two of which are independent because of the conservation of the moment vector. These equations may be given the following form: The system Eqs. (31), (32) has an oscillatory solution sin( ( − 0 )), = , with The period of the solution equals = 2 ∕ , and for mode number where is the length of the centreline closed after ( ∈ Z + ) turns, so we can write = 2 √ 2 + 2 0 = 2 √ . Thus, we arrive at where ∶= ( ∕ ) 2 . We can now replace in the above equation and write 2 = which can be rewritten as the equation where the coefficients are = ( − 1)( 0 − 1), = 0 ( − 2), = − 0 2 . For = 1, Eq. (33) is linear and the root is negative, therefore we discard this case. For ≠ 1, Eq. (33) is quadratic with roots This means that ( ) has no positive roots in this case. For > 1, > 0, hence − < 0 and must be discarded, while (1) = 1− − 0 (1+ ) < 0 implies + > 1, which is consistent with its definition. We conclude that the number of turns must always exceed the mode number . The Hamiltonian is given by = 0 2 (1 + + ).
We will use these linearised solutions to compute non-inflectional annulus configurations in Section 4.3.

Symmetry
Inspection of Fig. 1, or simple experimentation with plastic or paper annular strips, reveals that moderately overcurved closed strips possess antiprismatic symmetry, or -symmetry in the Schönflies notation of point symmetry groups in 3D (Müller, 2006). An object with -symmetry has an -fold central axis of rotational symmetry (i.e., rotation by an angle 2 ∕ leaves the object invariant), side axes of 2-fold symmetry (i.e., -rotation symmetry) perpendicular to the central axis (together giving symmetry) and in addition mirror-symmetry planes containing the central axis and passing between 2-fold axes (see Fig. 4). 1 -symmetry (equal to 2ℎ -symmetry) is degenerate in having a continuum of central axes of symmetry under rotations by 2 ∕ = 2 : any line in the plane of mirror symmetry that intersects the side axis of -rotation symmetry is a central axis (see, e.g., Figs. 5 and 8).
This high-order symmetry has consequences for the forces and moments that can exist in a deformed strip having this symmetry. Specifically, the force must be perpendicular to a mirror-symmetry plane cutting the structure, while the moment, as an axial vector, must lie in the plane. Furthermore, at a point of the structure through which passes an axis of rotational symmetry there can neither be a force nor a moment component along that axis. Since for ≥ 1 there is always a rotational side axis that does not belong to a mirror plane, it follows immediately that the force vector must vanish. But the force vector is constant in space, so we conclude that the force must vanish identically, = . This in turn means that the moment vector is conserved, i.e., = . The condition that it must belong to every mirror plane then implies that the moment is directed along the central symmetry axis (or, in case of = 1, along one of the central axes).
We note that the symmetry group is also characteristic of equilibrium shapes of excess cones (e-cones) (see Müller et al. (2008)), which may be considered as limiting cases of an annulus when its hole shrinks to a point, i.e., when → 1.
In the absence of the side axes of 2-fold symmetry, -symmetry reduces to -symmetry ( being a subgroup of ). Since for ≥ 2 the force vector must then be perpendicular to at least two non-parallel planes, this lower-order symmetry still implies a zero force for ≥ 2. However, 1 -symmetry does not imply a zero force, as the (single) axis of rotational symmetry lies in the (single) mirror plane. A non-zero force must then be normal to this plane (and hence to this axis).

Boundary-value problem
Motivated by the above symmetry observations we focus on -symmetric shapes for both undercurved ( < 2 ∕ ) and overcurved strips ( > 2 ∕ , the length of the closed centreline). We assemble a closed strip from 2 pairs of fragments (building blocks), see Fig. 4. Each pair consists of a fragment, of length ∕(4 ), and its mirror image glued together. We formulate a boundary-value problem for the arclength interval [0, ∕(4 )].
Note that the 3 ( ∕(4 )), ( ∕(4 )) and ( ∕(4 )) conditions together enforce the constitutive relation (15) (taking the explicit form (30)) that was differentiated in deriving our system of Eqs. (26). The conditions (34)-(39) place the solution in the fixed-point set of the reversing involution at = 0 and in the fixed-point set of the reversing involution at = ∕(4 ). This makes the plane through 0 ∶= (0) and the central axis a plane of mirror symmetry and the line through 1 ∶= ( ∕(4 )) and perpendicular to the central axis a side axis of 2-fold symmetry (see Fig. 4). Consequently, we have an inflection of the surface at 1 , where goes through zero. At both ends of the elementary fragment the generator is orthogonal to the centreline ( (0) = ( ∕(4 )) = 0). At = 0 it therefore lies in a symmetry plane, while at = ∕(4 ) it coincides with a side axis of rotational symmetry. The six kinematics conditions (40)-(45) are arbitrary (but for the requirement to avoid Euler-angle singularities) and fix the position and orientation of the strip in space.
Consider the pair consisting of the fragment from = 0 to = ∕(4 ) and its mirror image in the plane passing through 0 and with normal vector 0 . The vector̃1 ∶= 1 −2( 0 ⋅ 1 ) 0 , where 1 = ( ∕(4 )), is the mirror image of 1 with respect to this plane. Our pair of fragments havẽ1 and 1 at their ends. Let the angle between them be . Then we can write cos = 1 − 2( 0 ⋅ 1 ) 2 . For a closed strip we require 2 = 2 (for the planar multicovered annulus, is simply the number of turns around the central axis). This gives us the final boundary condition The two signs correspond to the two symmetrically related solutions, with opposite signs of , under the involution . We are taking > 0 so we have to accept both signs in Eq. (46), which will give different solutions. This final condition (46) guarantees that concatenation of 2 copies of pairs each made of the elementary piece and its mirror image, produces a closed strip (with 2 inflections). By substituting the Darboux frame vectors from Eq. (23), this condition may be expressed in terms of Euler angles.
Let 0 ∶= (0). Then the generator lines at both ends 0 + 0 and ∶= 1 + 1 cross the central axis for some , ∈ R, so that we get the vectorial equation where = ‖ ‖ is a unit vector along the central axis (and the line of the moment) and ∈ R is a third unknown. By solving the above equation we find , , ; in particular, The central symmetry axis is now defined by computing one point that lies on it, for instance . Eq. (47) has no solution if 0 ⋅ 1 = 0 (i.e., if the rotational axis is parallel to the mirror-symmetry plane) unless that axis belongs to the plane. If it does not belong to the plane then the assembled shape lacks -symmetry and the concatenation results in an open periodic configuration. If it does belong to the plane then the strip has 1 symmetry. In both cases the solution may have a non-vanishing force directed normal to the symmetry plane.
The boundary-value problem is solved numerically by the continuation code AUTO (Doedel et al., 2007). There are numerical difficulties solving this problem as a result of the logarithmic singularity in ( , ′ ), corresponding to intersection of generators on the edge of the strip. The issues are the same as in the case of rectangular strips and we refer to Starostin and van der Heijden (2015) for further discussion of our numerical approach. Singularities are found to occur when = 0 and as a result of the divergence of all coefficients , in Eqs. (17) and (18) diverge then as well. On the other hand, 1 never approaches zero. The only measure we take in our computations to circumvent singularities is that we allow to stay away from zero at the right end point of the integration interval.
( ∕(4 )) can therefore be interpreted as a regularising parameter. We try to make it as small as possible, typically reaching values between 10 −2 and 10 −7 , depending on the solution at hand. With this approach, singularities can only be accommodated at end points of the integration interval. This is not a limitation for the symmetric solutions we are here interested in, which never develop a singularity in the interior of the interval under the various parameter continuations. The generator also never becomes tangent to the inner edge of the strip.

1 -symmetric
We start by presenting results for varying length (or ) in the case = 1, taking = 0.5, = 0.5. Fig. 5 shows a sequence of four undercurved ( < 0) 1 -symmetric solutions. The left column shows the shape in space coloured according to the elastic bending energy density (from violet for low to red for high bending). The middle column shows the corresponding development in the plane, while the right column gives the normal curvature and geodesic torsion for a quarter of the centreline. In this figure, and in subsequent figures, red curves outside the strip show the edge of regression.
We see that even for very small length deficit (top row), the solution has two inflections, at = ∕4 and = 3 ∕4. At these points the edge of regression touches the outer edge of the strip, corresponding to logarithmic singularities of ( , ′ ) and hence of the bending energy density. Then also = 0, so the generator is perpendicular to the centreline (both generators belonging to the same straight line in space). Out of the singularities radiate flat triangular regions. The same spontaneous bending features were found in one-sided (Starostin and van der Heijden, 2015) and twisted open (Korte et al., 2011) inextensible strips and are more generally seen in crumpled membranes (Witten, 2009).
As the length of the strip is reduced the solution approaches a folded state with sharp creases connecting four vertices (bottom row in Fig. 5). The edge of regression comes very close to the inner edge of the strip at two points, but the cusp point remains off the strip until the limiting planar multi-covered state is reached. This limiting folded state has the shape of a double-covered rhombus and has length = 8 arctan( ) (see Appendix B for details). For = 0.5, = 0.5, this gives a length deficit of −2 ∕ = −8.64671 ( = −0.68808). For comparison, the limiting state of a rectangular strip with 4 symmetry has lim →0 = 8 and the final configuration is a double-covered square with sides 2 √ 2 so that each edge of the strip makes up one of the diagonals of the square while covering it twice (Starostin and van der Heijden, 2015). In another limiting case, → 1∕ implies → 2 ∕ and the strip becomes a flat disc with no hole. The variation of the normalised magnitude of the moment under varying length deficit for these solutions (for = 1.0) is shown by the dashed curve in Fig. 6. The moment diverges as the length shrinks to zero. A flat rhombus may be folded in a different way from what is illustrated in Fig. 5. If the outer and inner edges are swapped we obtain the shape shown at the bottom row of Fig. 7. The length of the strip in both flat limits is given by the same expression = 8 arctan( ) (see Appendix B). If we now gradually increase the length we trace another branch of solutions carrying the 'inverted' equilibria displayed in Fig. 7. The strip of length = 2 ∕ of the unperturbed plane annulus (i.e., = 0) is now significantly warped. As the strip becomes longer it starts to intersect itself. Continuing with the length increase we eventually approach the flat triple-covered annulus for the limiting length = 3 × 2 ∕ (i.e., = 2). Dependence of the normalised moment on the relative length deficit/excess for this branch of strip solutions is given by the dotted = 1 curve in Fig. 6.
We also include here in Fig. 8 a similar equilibrium with the same symmetries but for a rectangular strip, with geodesic centreline (i.e., = 0). The simple closed two-sided solution for such a strip is a circular cylinder, but this 'inverted' state is also an equilibrium and can indeed be folded from a paper strip (suggesting the solution is stable). The strip folds into a flat square with diagonal 4 when its aspect ratio tends to the critical value = 8 (see Eq. (61)).

2 -symmetric
Undercurved 2 -symmetric solutions are shown in Fig. 9. Singularities occur on the outer edge of the strip (close inspection reveals tiny gaps between cusps of the edge of regression and the inner edge of the strip). There is a critical length deficit at which the strip touches itself (bottom row in the figure). By symmetry, a pair of contacts occur simultaneously, both on the central axis. The condition for this contact can be obtained from Eq. (48), which can be rewritten as The inner edge of the strip touches itself when = − . Overcurved 2 -symmetric solutions at increasing length are shown in Fig. 10. The top row displays the familiar saddle shape obtained for small excess length. We see that singularities now occur on the inner edge of the strip. On the second row we see the solution at first self-contact (of the outer edges of the strip), when Eq. (49) is satisfied for = . Longer strips cross themselves and the limiting state is a flat triple-covered annulus of length = 3 × 2 ∕ . Similar spontaneous folding of a single-covered planar equilibrium into a triple-covered planar equilibrium has been observed in growing thin isotropic intrinsically curved rods (Moulton et al., 2013) and in slender overcurved structures (Goto et al., 1992;Mouthuy et al., 2012;Audoly and Seffen, 2015). The behaviour is confirmed in the bifurcation diagram in Fig. 6, which also shows a discontinuity in the limiting value of between under-and overcurved folding at = 0.  -symmetric folding diagram for the first four folding modes, = 1 (red), 2 (blue), 3 (green) and 4 (violet): normalised moment | | against relative excess length for = 1.0, = 0.5. Solid curves are for overcurved 'regular' folding ( > 0) with solutions going planar at = 0 (see, e.g., Fig. 5). Dashed curves are for undercurved 'regular' folding ( < 0). Dotted curves are for 'inverted' folding, passing smoothly through = 0, where the solution is not planar. At large , 'regular' solutions fold into closed (2 − 1)-covered planar annuli having = 2 − 2. The vertical asymptotes correspond to flat folded configurations with given by in Eq. (60) ( = 1). Fig. 11 at fixed = 1. All three shapes have two pairs of contacting points on their outer edge. The graphs in the right column show that the curvature of the centreline (equal to √ 2 + 2 ) increases with the width, while the deviation of the generators from the direction orthogonal to the centreline (measured by = ∕ ) becomes smaller. Thus we may say that the surface is getting closer to conical as is increased. Another characteristic feature that is seen in these graphs is that as the strip becomes narrower, the curvature develops a knee in the approach to the inflection. This was also observed for the curvature of the Möbius strip van der Heijden, 2007, 2015).

Dependence of the shape at first contact on the width of the strip is demonstrated by the series of solutions in
Similar to the 1 -symmetric case, there exists another branch of 'inverted' 2 -symmetric solutions, also included in Fig. 6. In the large-(i.e., large-) limit it starts with a flat 5-times-covered annulus of length = 5 × 2 ∕ . Solutions on this branch are shown in Fig. 12, the top solution being close to the 5-times-covered limit. As the length is decreased the shape becomes free of self-intersections and the strip turns only once around the central symmetry axis (see the bottom row of the figure, which is for length equal to 2 ∕ , the length of the unperturbed planar annulus). For < 0 the strip eventually collapses into a compact state at the value given by Eq. (60) for = 2 (similar to the sequence in Fig. 7 for 1 -symmetry, but now self-intersecting). Fig. 13 shows three undercurved 3 -symmetric solutions. The top solution is for a tiny length deficit. It has six inflections with logarithmic singularities (stress localisations) on the outer edge of the strip and six alternating near-singularities (cusps of the edge of regression just off the strip) on the inner edge. The ridges connecting the singularities on the outer edge and the points of high curvature on the inner edge form a star-like polygon that divides the strip into twelve flat triangular domains. As the strip shrinks, six pairwise contacts (three at the top and three at the bottom) occur simultaneously on the inner edge at some critical length deficit (closely approximated by the bottom shape in the figure). Unlike the case for = 2, the contacts are located off the central symmetry axis.

3 -symmetric
Overcurved 3 -symmetric solutions of increasing length are shown in Fig. 14. For small solutions have a monkey-saddle-like shape that becomes more wavy as the length increases. Six singularities lie on the inner edge of the strip, where the surface has inflections. At some critical length six pairs of points on the outer edges touch each other. Further increase of length leads to self-intersections and the strip eventually collapses into a 5-covered annulus of length = 5 × 2 ∕ .
Similar to the case of 2 -symmetry, there exists another branch of 'inverted' solutions with the same 3 -symmetry (see Fig. 6). Fig. 15 presents three strips on this branch. The top row shows a shape that is a perturbation of the 7-covered flat annulus with selfintersections for a length slightly shorter than 7 × 2 ∕ . Further contraction leads to a 3 + 3-lobes shape without self-intersections. An example is shown on the bottom row of Fig. 15, which is for length exactly equal to the length 2 ∕ of the undeformed annulus. For < 0 the strip eventually collapses into a compact state at the value given by Eq. (60) for = 3 and = 1 (similar to the sequence in Fig. 7 for 1 -symmetry). Fig. 16 presents three overcurved solutions with 4 symmetry. The top row configuration, which is already quite far from the flat unperturbed annulus, has an articulated wavy look resembling a ruff. Further increase of length causes multiple self-intersections and ends in a 7-covered flat annulus. Again, there also exists a branch of 'inverted' solutions (included in Fig. 6), but its further discussion is omitted.

Non-inflectional solutions
For the inflectional solutions we constructed solutions from fragments of 1∕(4 )th of the strip. This was possible because of the -symmetry but was also necessary because of the singularities (stress localisations) at inflections of the surface (where = 0) that prevent numerical integration through such points unless special measures are taken. Here we compute non-inflectional solutions for which we solve our 15th-order system (26) over the full arclength domain [0, ] subject to the following 15 boundary conditions consisting of: 5 conditions at = 0, Note that the 3 (0) and (0) conditions together enforce the constitutive relation (15). Solutions with mode numbers ( , ) can be computed by starting from the multi-covered circles calculated in Section 3.9 and then varying the length (or any other parameter). Not having inflection points, which used to coincide with side axes of 2-fold rotational symmetry, we now do not have these side symmetry axes and the symmetry group is reduced from to . In particular, the perturbed conical solutions constructed in Section 3.9 are -symmetric. Figs. 17 and 18 illustrate how the shape of an annular 2 -symmetric strip of unit width and radius 1∕ = 10 with ( , ) = (3, 2) varies as its length (and hence ) decreases. The top figure in Fig. 17 shows a slightly perturbed triplecovered solution; the bottom corresponds to a shorter strip ( < 0). As the length is reduced further the self-intersections disappear and the saddle-like strip approaches a flat folded state, as illustrated in Fig. 18. Fig. 19 similarly shows an ( , ) = (9, 2) solution for two different lengths. Neither of these solutions features singularities of the strain energy density ℎ (i.e., stress concentrations) on the edge of the strip, as confirmed by the edge of regression (in red) staying away from the edge of the strip, although such singularities are approached in the flat folded limit (see Fig. 18). Note that the self-intersections of these solutions are not isolated events. They persist under parameter variations and suggest that the centreline of the strip lies on an imaginary surface. The same phenomenon is observed for forceless narrow rectangular strips, which lie on a sphere (Starostin and van der Heijden, 2018). It is not clear on what ( -dependent) surface these annular strips lie.
Note that this system is semi-decoupled: the first equation contains only .
The first derivatives vanish at the critical points, all of which are located on the -axis ( ′ = 0), one in the origin = 0 and a pair at = ± ⋆ , ⋆ = √ 1 √  2 − 1 (provided  > 2 2 ). For the critical point at the origin, = 16( − 2 2 ), ′ = 0, ′ ′ = 32. This critical point is therefore a (stable) centre if  > 2 2 and a (unstable) saddle otherwise. For = ⋆ , ′ = 0, we find = −64 2 ⋆ < 0, ′ = 0, ′ ′ = 32(1 + 2 ⋆ ) 2 > 0 and we conclude that these critical points are (unstable) saddles, when they exist. We therefore have a subcritical pitchfork bifurcation under variation of  (i.e., the value of the Hamiltonian), with an unstable critical point becoming stable with the simultaneous birth of a pair of unstable critical points. The critical points ± ⋆ move away from the origin and tend to infinity as approaches zero, i.e., in the geodesic limit (indeed, for the rectangular strip there is only a (non-bifurcating) critical point at the origin Starostin and van der Heijden, 2018). Fig. 20 shows the phase portrait for parameters  = 0.5, = 0.1.  The colouring is according to the magnitude of the normal curvature, which vanishes at the upper separatrix (thick black line) and is imaginary in the region above the separatrix (coloured grey, except for the small white subregions where the moment is imaginary as well). The saddles are located at ± ⋆ = ±2. Substituting = ⋆ , ′ = 0 into Eq. (53) gives for the saddles and separatrices (heteroclinic orbits connecting the saddles) 2 = 2 ⋆ =  2 2 . The heteroclinic orbits may be found explicitly as follows. Setting 2 =  2 2 in Eq. (53) we obtain and hence For  > 2 2 , evaluation of the integral yields

This gives an implicit equation for as a function of arclength .
A similar calculation for  < 2 2 yields while for  = 2 2 we get the simpler result The heteroclinic solution, with corresponding pince-nez shape of the narrow annular strip, is displayed in Fig. 21.
To find the normal curvature for the heteroclinic solution, we first rewrite Eq. (54) as Here the '+' sign corresponds to the separatrix with ′ > 0 for | | < ⋆ or the separatrix with ′ < 0 for | | > ⋆ , while the '−' sign corresponds to the other separatrices (reflected about the axis). An expression for can then be obtained from Eq. (51) after Heteroclinic pince-nez-shaped solution for the narrow annular strip with  = 0.5, = 0.1 (four views). The shape is drawn for half-width = 0.9 and coloured to show the elastic bending energy density (violet for low, red for high bending). The moment vector is drawn as a black arrow. The bottom row shows graphs of ( ) and the normal curvature (solid) and geodesic torsion (dashed) of the centreline.
elimination of the derivative with help of Eq. (56). The result for the '−'-curve is The right-hand side of the above is non-negative only for | | ≤ ⋆ . The '+'-curve has identically-zero . It separates physical solutions from unphysical ones (see Fig. 20). The rings in Fig. 21, which correspond to the saddles in Fig. 20, have = 0 and therefore have radius 1∕ .

Linearisation
Linearisation of Eq. (50) yields We again assume that the normal curvature is the sum of a constant and a small variable term: ( ) = 0 (1 + ( )), 0 = const, 0 ≠ 0. Then, keeping only the zero-and first-order terms in Eq. (51), we find  = 2 2 + 2 0 and = − , we find = 2 √ 2 − 2 2 , 2 0 = 2 2 2 −2 2 2 ,  = 2( 2 − 2 ) 2 −2 2 2 and 2 = 4 2 2 −2 2 2 . We conclude that the number of turns must exceed ⌈ √ 2 ⌉. Fig. 22 shows a bifurcation diagram with ( , ) modes bifurcating from flat multi-covered configurations for = 3, 5 and 7 (having, respectively, 1, 2 and 3 bifurcating branches). To obtain this diagram we solve the system of Eqs. (50), (51), coupled to the kinematics equations (24), (25), imposing the boundary condition (0) = 1 and the closure conditions (0) = ( ) and (0) = ( ), while treating the value of the Hamiltonian, , as a free parameter. Bifurcation points are given by = √ 2 − 2 2 −1 and = 4 , in agreement with the calculation in Section 5.2. All bifurcating modes are non-inflectional and have -symmetry. All solutions remain closed along the bifurcating branches in Fig. 22 except for the = 1 modes, which have (0) ≠ ( ) (conversely, if we impose (0) = ( ) then the tangent is non-closed). We recall from Section 4.1 that 1 -symmetry is exceptional: non-closure for = 1 is possible because in 1 -symmetry the axis of rotational symmetry lies in the mirror plane. We note that for each the branch for the highest possible ends at = −1, i.e., = 0, with finite . A selection of closed solutions is shown in Figs. 23-26 with phase-space orbits indicated. Strip configurations are drawn with a small width and coloured to illustrate how they twist and bend in space. As in the case of the wide non-inflectional solutions in Section 4.3 (and the case of narrow rectangular strips Starostin and van der Heijden, 2018), the configurations have persistent self-intersections, suggesting that narrow strips lie on some imaginary surface. Solutions on orbits near the fixed point at the origin are close to the flat multi-covered limit, while solutions on orbits near the heteroclinic orbit demonstrate the pince-nez nature of these orbits, but in 'quantised' form, showing two-cycle solutions for = 2 (Fig. 25) and (lotus-flower-like) three-cycle solutions for = 3 (Fig. 26). Looking back at Fig. 19 we also note that the pince-nez nature of non-inflectional solutions is retained by wide strips.

Discussion
Over-and undercurved thin annular strips or sheets forced into closed structures have a natural tendency to adopt highly symmetric shapes and to spontaneously fold into compact multi-covered configurations under feed-in of more length (or varying intrinsic curvature). We have formulated and solved boundary-value problems for -symmetric equilibrium solutions of unstretchable annular strips. Continuation of the length parameter allows us to mimic the folding process for different solution modes governed by two integer numbers. Because of the high-order symmetry, closed solutions cannot have an internal force, i.e., the strips are forceless. We have considered both wide and narrow (strictly zero-width) strips. Narrow strips cannot have inflections, but wide strips can be either inflectional or non-inflectional.
Inflectional solutions are found to feature stress localisations on the edge of the strip at generators of surface inflection and corresponding to logarithmic singularities of the strain energy density. Numerics suggests that all inflections induce such infinite stress localisation and all such stress localisations are associated with inflections. It is not clear why this should be so. This behaviour is similar to what was found for the geodesic case of rectangular strips (  Starostin and van der Heijden, 2015), but note that in that case inflections correspond to inflection points of the (straight) centreline of the strip (where the standard curvature is zero), while here, in the non-geodesic, annular case, inflections correspond to centreline points with zero normal curvature (note that the centreline in this case cannot have an inflection point because 2 = 2 + 2 ≥ 2 > 0). No new inflections are created in the folding process as more length is inserted.
From our explorations of forceless folding the following picture emerges (see Fig. 6). For each mode number ≥ 1 there exist (at least) two paths of -symmetric inflectional solutions under variation of , representing 'regular' and 'inverted' spontaneous -mode folding under length feed-in. The path of 'regular' folding goes through the unstressed planar annulus (at = 0) and has both undercurved ( < 0, length deficit) and overcurved ( > 0, length excess) solutions. Undercurved shapes have 2 conical singularities (stress localisations) on the outside of the annulus, while overcurved shapes have them on the inside of the annulus.
For overcurved solutions the general rule appears to be that for every > 1 there is a continuous path of closed solutions from the planar annulus with = 0 to the (2 − 1)-covered planar annulus with = 2 − 2. For = 1 the path collapses into a single point. The path of 'inverted' folding has both undercurved and overcurved solutions, but does not go through the single-covered E.L. Starostin and G.H.M. van der Heijden Fig. 27. Three views of a periodic open-strip solution with ∞ -symmetry having zero force ( = 0.5, = 0.5) and a paper model. unstressed planar annulus. It connects the (2 + 1)-covered planar annulus having = 2 with a compact, planar, self-intersecting configuration having = (4 ∕ ) arctan( tan( ∕(4 ))) − 1. The singularities are on the inside of the annulus. The two types of folding, 'regular' and 'inverted', are related by the involution in which changes sign. One can therefore obtain an 'inverted' solution from a 'regular' one (or vice versa) by parameter continuation in to its opposite value (provided this homotopy exists). This property has proved useful in getting our complete set of solutions up to = 4. Note that the rectangular strip, having = 0, does not have this splitting into two different types of folding. In this sense, annular strips can be viewed as symmetry-broken.
The flat states at the ends of these folding paths (the end points of the bifurcation curves in Fig. 6) have a normal curvature that is infinitely peaked at the sharp folds and is zero elsewhere (see the approximate states at the bottom of Figs. 5 and 7). These solutions have the constant moment vector perpendicular to their flat surface. The only non-zero moment component of these flat states is therefore 2 , the Lagrange multiplier enforcing the geodesic constraint 2 = . Note that for forceless solutions it follows indeed immediately from the equilibrium equations (26) that when = 0, 2 is constant. In the limit → ∞, -symmetry degenerates into ∞ -symmetry (equal to ∞ℎ -symmetry). In that case the 2-fold rotation axes are parallel to the mirror planes. We can compute solutions with this symmetry by replacing the closure condition (46) by (0) ⋅ ( ∕4) = 0. Generically, the rotation axes will not lie in the mirror planes. The solution is then necessarily open and periodic (and the symmetry group is a translation group, not a point group). Fig. 27 shows three views of such a solution in the form of a 3D curved elastica shape (also easily produced by hand). This solution has inherited zero force, but because rotation axes and mirror planes are parallel, ∞ -symmetric solutions need not be forceless. Indeed we can now apply a force to close the solution. At the point of closure the rotation axes join the mirror planes (so we have point symmetry again) and the result is a solution with 1ℎ -symmetry (this 1ℎ being the only point group among and ℎ allowing for a non-zero force). Three solutions with this 1ℎ -symmetry (but for different values of ) are shown in Fig. 28. Generally, an object with ℎsymmetry has -symmetry with in addition a plane of mirror-symmetry perpendicular to the central symmetry axis and, as a consequence, also perpendicular mirror planes each containing the -fold central axis and one of the 2-fold side axes. For = 1 the 2-fold side axis lies in the symmetry plane, so there can be a non-zero force, which must then be perpendicular to the symmetry plane (and the 2-fold axis). The (constant) force vector is now along the central axis. Fig. 28 shows that for a relatively short strip the shape is an annular figure-of-eight (for a rectangular strip it would be a cylindrical figure-of-eight with centreline given by the planar figure-of-eight Euler elastica). The 2-fold symmetry axis coincides with the inflectional generator and the line of self-intersection. As the strip is made longer (larger ), its centreline becomes more warped with higher geodesic torsion. When the length approaches = 2 × 2 , the strip flattens to form a double-covered flat annulus with circular centreline (right image in Fig. 28). We have seen that non-inflectional forceless solutions pulled out of a multi-covered flat state have persistent self-intersections (see Figs. 17,19,24,25,26) that suggest that they lie on a surface (that may vary with the parameter that does the 'pulling', for instance ). This potentially has consequences for the topology of the centreline of such a closed strip. The classical four-vertex theorem for planar curves states that the curvature of a simple, closed, smooth plane curve has at least four local extrema (DeTurck et al., 2007). An extreme point of the curvature is called a vertex. Extensions of the four-vertex theorem to space curves exist but since space curves are much less constrained, extra conditions are usually imposed. For instance, the curve can be assumed to lie on a surface. Various results and conjectures exist for closed curves that lie on the boundary of a convex body. Importantly, in the context of space curves a vertex is defined as a point of vanishing torsion. The strongest result to date appears to be the theorem proved recently by Ghomi (2017). It states that under mild conditions a curve without inflections that bounds a convex surface has torsion that is either identically zero or changes sign four times.
Our results suggest that a version of the four-vertex theorem is relevant for the solutions of our boundary-value problems. Note that the various 'modes' bifurcating from the centre point under integer conditions (Section 5.2) have an increasing number of zeros of and hence of the geodesic torsion . We can compute the (Frenet) torsion of the centreline as follows, using the fact that the tangent is common to the Frenet frame ( , , ) of tangent, principal normal, and binormal, and our Darboux frame ( , , ). By writing out the triple product ⋅ (̇×̈) for both frames and using the corresponding (generalised) Frenet-Serret equations, we find = + ′ 2 + 2 .
We may therefore expect also to have an increasing number of zeros as increases. For the = 1 mode, however, we expect , and , to have only two zeros. (Note that zeros of occur at intersections of the phase-plane orbit in Fig. 20 with the vertical axis, where = 0 and ′ = 0, the latter indicated by a 'colour extremum'). If the centreline lies on a surface, as it seems to do, then the four-vertex theorem would appear to give restrictions on the type of surface we can have. Interestingly, we find that = 1 solutions (and generally multi-covered ( , ) solutions with ∕ = an integer), unlike solutions for higher , immediately open up as they bifurcate from the flat state (see Fig. 22) at the critical value of given in Section 5.2, as allowed by the exceptional nature of 1 -symmetry. By opening up these solutions of course avoid violation of any four-vertex theorem. However, by applying a force it is possible to close them. What happens when doing this, though, is that the torsion develops initially a tiny wiggle near the interior zero that creates two more zeros giving a total of four! (see Fig. 29, which actually shows wide-strip solutions obtained by additional continuation in from 0 to 0.5 immediately after bifurcation). Of course we know very little about any surface the strip might lie on, but it is interesting that the strip seems to 'feel' the constraint posed by a four-vertex theorem for space curves and shows behaviour that appears to avoid violating it.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
No data was used for the research described in the article.

Appendix A. Calculation of the strain tensor
By calculating the relevant strain tensor we verify that the deformation implicit in the parametrisation (1) is isometric, i.e., does not involve stretching of the material of the strip. The calculation follows a similar calculation for the geodesic case in van der Heijden and Starostin (2022), which disproved the erroneous claim of non-isometry made in Chen and Fried (2016).
We first need to find the parametrisation of the reference configuration of the strip in the plane. Referring to Fig. 2 we can write for this reference configuration = cos 1 + sin 2 .
Here ( 1 , 2 ) is a Cartesian coordinate system with origin at the centroid of the annulus. We use as coordinates of an arbitrary point the height and the arclength (measured from the 1 axis) of the centreline point where the generator through intersects. Thus we consider = ( , ).
For the deformed configuration = ( , ), given in Eq.

Appendix B. Geometry of flat folded states
Here we derive geometrical properties of the flat folded configurations of -or -symmetric annular strips. Consider the building block of the strip bounded by two generators aligned with 0 = (0) and 1 = ( ) (recall that = 0 at both ends of the fragment, which implies alignment of the generators with the vector ) (see Fig. 30). The generator at = 0 crosses the inner and outer edges of the strip at points ′′ and , respectively. The generator at = = ∕(4 ) crosses the inner and outer edges of the strip at points and ′ , respectively. The fold line passes diagonally through points and so that the generator ′′ maps to the line ′ . The preceding neighbouring piece is a mirror image of our folded building block with respect to the plane orthogonal to the strip and passing through ′ . The following neighbouring piece is a congruent copy of our fragment glued along the line ′ . Let be perpendicular to ′ and be perpendicular to ′ . We denote = ∠ and = ∠ . For the closed strip the angle between (0) and 1 must be 2 − 2 , but the same angle equals 2 − ( + ) in the folded state, hence we have + = 2 . Now consider the triangle , where is the centre of the annulus, so that | | = − and | | = + and the angle ∠ = ∕ . As ∠ = − , ∠ = − . Applying the sine rule to , we get (1 + ) sin = (1 − ) sin , which allows us to find The length of the closed centreline is = 4 = 4 ( − ) = 4 ( 2 − 2 ) , from which we finally obtain In the limit of a straight centreline this gives lim →0 = 8 tan ( 4 ) .
Also note that in the limit → ∞, Eqs. (60) and (61) both give lim →∞ = 2 , independent of . Thus, folded strips within our family can be immersed in R 3 if their ratio ∕(2 ) is larger than . For comparison, a smooth developable Möbius strip can be immersed in R 3 only if this ratio ∕(2 ) is larger than ∕2 (Halpern and Weaver, 1977).