## Abstract

Folded surfaces of origami tessellations have attracted much attention because they often exhibit nontrivial behaviors. It is known that cylindrical folded surfaces of waterbomb tessellation called waterbomb tube can transform into peculiar wave-like surfaces, but the theoretical reason why wave-like surfaces arise has been unclear. In this paper, we provide a kinematic model of waterbomb tube by parameterizing the geometry of a module of waterbomb tessellation and derive a recurrence relation between the modules. Through the visualization of the configurations of waterbomb tubes under the proposed kinematic model, we classify solutions into three classes: cylinder solution, wave-like solution, and finite solution. Through the stability and bifurcation analyses of the dynamical system, we investigate how the behavior of waterbomb tube changes when the crease pattern is changed. Furthermore, we prove the existence of a wave-like solution around one of the cylinder solutions.

## 1 Introduction

Origami tessellations are origami obtained by tiling translational copies of a modular crease pattern. Origami tessellations can be folded from flat sheets of paper and transform into various shapes. Even though origami tessellations are corrugated polyhedral surfaces, they sometimes approximate smooth surfaces of various curvatures, including synclastic and anticlastic surfaces [1]. Such macroscopic surfaces exhibit nontrivial behaviors we cannot expect from the periodicity of its crease pattern and recently attracted much attention of scientists and engineers. For example, Schenk and Guest [2] showed that Miura-ori forms anticlastic surfaces while egg-box pattern forms synclastic surfaces through numerical modeling using bar-and-hinge model. Nasser et al. [3,4] analyzed the approximated surface of Miura-ori and egg-box pattern based on the knowledge of differential geometry.

In this paper, we focus on folded surfaces of waterbomb tessellation, specifically, cylindrical folded surfaces called waterbomb tube. Waterbomb tube is known as origami work “Namako” [5] or more commonly, the name “Magic Ball.” Based on the property of waterbomb tube, various engineering applications are attempted, such as origami-stent graft [6], earthworm-like robot [7], soft gripper [8], and robot wheel [9] using the property that the radius of the tube can vary. The dynamics of the stent graft and the wheel inspired by waterbomb tube has been analyzed [10,11]. Also, the kinematics and structural deformation of waterbomb tube has been studied [1214]. One of the most interesting phenomena reported on waterbomb tube is that it can form wave-like surfaces [15,16] like in Fig. 1. This is the unique phenomenon not yet observed in the folded surfaces of other tessellations such as Miura-ori and egg-box pattern; however, the theoretical reason why wave-like surfaces arise has been unclear.

Fig. 1
Fig. 1
Close modal

Here, our objective is to know “why” waterbomb tube produces wave-like surfaces, i.e., to clarify the mathematics behind the behavior. In this paper, we model waterbomb tube as the sequence of modules and focus on its recurrence behavior instead of simultaneously solving the network of 6R spherical linkages as in [12,13,15,17]. We first provide the kinematic model of each module of waterbomb tube and the relation with adjacent modules to obtain the recurrence relation dominating the folded states of modules based on the rigid origami model with symmetry assumption (Sec. 2). Then, through computation and visualization of the folded states of waterbomb tube using the recurrence relation, we observe that the solutions fall into three types: cylinder solution, wave-like solution, and finite solution (Sec. 3). Furthermore, by applying the stability and bifurcation analyses of the dynamical system of waterbomb, we illustrate how the behavior of the system, i.e., whether waterbomb tube can be wave-like surface or not, changes when the crease-pattern parameters are changed (Sec. 4). Finally, we prove the existence of wave-like solution by applying theorems of discrete dynamical systems (Sec. 5).

## 2 Model

### 2.1 Definition.

First, we introduce parameters of the crease pattern of waterbomb tessellation. We can obtain the crease pattern of waterbomb tessellation by tiling translational copies of a unit module. The entire pattern is controlled by four parameters (Fig. 2): a unit module shown in Fig. 2, left, is controlled by two parameters α, β (α ∈ (0, 90°), β ∈ (0, 180° − α)), and the repeating number of modules in column and row directions is given by two integers N and M ($N∈Z>0,M∈Z>0$).

Fig. 2
Fig. 2
Close modal

We consider a waterbomb tube as the rigid folded state of the crease pattern of waterbomb tessellation, such that the left and right sides in Fig. 2 are connected to form a cylindrical form. In our model, we assume N-fold symmetry about an axis and mirror-symmetry about N planes passing through the axis as in the existing research [13] (see Fig. 3). However, unlike the previous research [13], we do not assume mirror symmetry with respect to a plane perpendicular to the axis. Here, the behavior of waterbomb tube is governed by three parameters (α, β, N), because M specifies the finite subset interval of infinitely continuing waterbomb tube.

Fig. 3
Fig. 3
Close modal

### 2.2 Kinematics of Module.

First, we consider the degrees of freedom (DOF) of each module. Because the kinematic DOF of n-valent vertex in general is n − 3 [18,19], the DOF of waterbomb module (without symmetry) is 3. With the assumed mirror symmetry, the DOF drops by one more degree. Therefore, the module has 6 − 3 − 1 = 2 DOF.

Here, we focus on the module with correct mountain and valley (MV) assignment and pop-up pop-down assignment, where we assume that the center vertex of the module is popped-down, i.e., the sum of fold angles (positive for valley and negative for mountain) are positive, and other vertices are popped-down. The definition of pop-up and pop-down follows that of [20], and it is known that rigid origami cannot continuously transform from popped-up and popped-down state without being completely flat.

To represent the folded states of modules, we consider the isosceles trapezoid formed by two pairs of mirror reflected vertices of the module (Figure 4, left). The parallel bottom and top edges are aligned in the hoop direction, while the hypotenuses are along the longitudinal direction. The length of hypotenuses of the trapezoid is fixed at 2cosα, while the half lengths of bottom and top edges, denoted by x and y, change in (0, sinα) by folding the module. For each given pair of x and y in (0, sinα), there exists eight different states of a module (see Fig. 5). However, only one of the folded state has the correct MV and pop-up/pop-down assignment for the following reason, thus we may safely use x and y as the parameters to uniquely represent the folded state of the module.

Fig. 4
Fig. 4
Close modal
Fig. 5
Fig. 5
Close modal

To show the uniqueness, first notice that the top and bottom modules illustrated in the same column in Fig. 5 are the mirror reflection of each other through the plane containing the isosceles trapezoid, so the top states are popped-down and the bottom states are popped-up. Therefore, the folded state needs to be one of the four top states. Within the top row we compare the first and second left configurations. Vertex P′ is the reflection of vertex P through the plane defined by vertices O, A, and B, so crease OP is mountain and OP′ is valley, thus we choose vertex P (counterclockwise side of cycle OAB) is selected as the correct side. In a similar manner, vertex Q is chosen over Q′. Thus, the top-left configuration with P and Q is chosen. Finally, we can show that the top-left configuration has indeed correct MV-assignment also for other creases; this can be verified by checking that P, D, and C are on the counterclockwise side of OAB and that P is on the symmetry plane between C and D. Therefore, only the top-left configuration has the correct MV-assignment.

### 2.3 Kinematics of Entire Tube.

Next, we consider representing the folded states of an entire waterbomb tube. First, we take out the mth hoop of the waterbomb tube. Because of rotational symmetry, the configuration of the mth hoop can be represented by using parameter xm, ym for each congruent module and the common dihedral angles between isosceles trapezoids 2ρm (Fig. 4, right). In order that the mth hoop to be a closed cylinder, we have the following constraint (see Appendix  A for derivation).
$cosρm=−(xm−ym)2sec2α+4sin2πN4−(xm−ym)2sec2α$
(1)
Because 2ρm can be derived as the function of xm and ym, the configuration of each hoop can be uniquely represented by parameters xm, ym. However, note that not all values of (xm, ym) represent valid states as cosρm may become nonreal number for certain given pair of values.
Therefore, we can represent folded states of waterbomb tube by the following sequence (Figure 4, right):
$(xm)m=0,1,…,M−1,(ym)m=0,1,…,M−1$
(2)
In order that two consecutive hoops are compatible, there are constraints between xm, ym and xm+1, ym+1 that need to be satisfied. The constraints result in the significant property of waterbomb tube that the state of mth hoop uniquely determines that of m + 1th hoop. This is confirmed by the fact that any vertex of m + 1th modules can be solved by sequentially applying three-sphere-intersections three times (see Fig. 6). In each step of three-sphere-intersection, position-unknown vertex X incident to three position-known vertices V1, V2, V3 are fixed by constructing three spheres S1, S2, S3 of radii XV1, XV2, XV3 centering at V1, V2, V3, respectively, and taking their intersection. In each step, there are at most two candidates for X, denoted by X and X′. They are the mirror reflection of each other through the plane defined by V1, V2, and V3. Thus, MV-assignments of creases XV3 and XV3 are the inversions of each other, and we choose the one consistent with MV-assignments shown in Fig. 2 (for more details, see Appendix  B). Therefore, the m + 1th module with consistent MV-assignment is uniquely determined from the mth module. Note that the existence of a solution is not guaranteed because intersections of three spheres may be empty.
Fig. 6
Fig. 6
Close modal
When the solution exists, the relationship between mth hoop and m + 1th hoop can be formulated by following recurrence relation:
${xm+1=f(xm,ym;α,β,N)ym+1=g(xm,ym;α,β,N)$
(3)
where function f and g are complex nonlinear function which parameters are α, β, and N, but can be analytically derived (see Appendix  B). Consequently, waterbomb tube can be interpreted as a 2-DOF mechanism for an arbitrary M, whose entire shape is determined if initial values x0, y0 (0 < x0, y0 < sinα) are given.

### 2.4 Kinematics Interpreted As Dynamical Systems.

Hereafter, we clarify the mathematics behind wave-like surfaces by interpreting the recurrence relation (3) as a discrete dynamical system. Recurrence relation (3) is categorized as a nonlinear discrete dynamical system in two dimension, as functions f, g are nonlinear functions of two variables xm, ym. Also, the system (3) is autonomous system because the function f, g is independent of indices m.

An important property of the system is reversibility. Because the module of waterbomb tessellation is point symmetric about the central vertex (see Fig. 2, left), the following equation holds $(m=1,…,M−1)$:
${xm−1=g(ym,xm)ym−1=f(ym,xm)$
(4)
From the above formula, if we define the map $F:(x,y)↦(f(x,y),g(x,y))$ and linear involution $G:(x,y)↦(y,x)$, the following one holds:
$G=F∘G∘F$
(5)
Like Eq. (5), the system where exist some involution reversing time direction called reversible dynamical system [21], and involution G is the symmetry of the system.

## 3 Visualization of Configuration

To understand the 2-DOF kinematics of waterbomb tube, we computed the sequences (2) under different pairs of initial values (x0, y0) and visualized their configurations as the sequence of points in x, y-plane, i.e., the phase space. The sequence ${(xm,ym)=Fm(x0,y0)|m∈Z}$ is called orbit of (x0, y0) and the plot of them is called phase diagram. Here, we fix the crease-pattern parameters (α, β, N, M) = (45°, 34.6154°, 24, 100) and show the phase diagram under this parameter in Fig. 7. Under this parameter, we can observe all possible types of solutions we explain below.

Fig. 7
Fig. 7
Close modal

### 3.1 Three Types of Solutions.

It can be observed that the characteristics of the solution change depending on the given initial values: the solutions can be classified into three types; cylinder solution where it forms a constant radius cylinder, wave-like solution where the corresponding waterbomb tube form wave-like surface, and finite solution where the system fails to obtain a solution after some iterations. Note that the behavior of the system (3) can change under different crease-pattern parameters as we discuss in Sec. 4 (i.e., we cannot observe some types of solutions in different parameters). However, we can still classify solutions under different parameters into these three types.

#### 3.1 Cylinder Solution.

The top and bottom of waterbomb tube in Fig. 7, left, forms a constant-radius cylinder that corresponds to fixed two points shown in black on the phase space. In this type of solution folded states of modules composing waterbomb tube are fixed (i.e. $(x0,y0)=(x1,y1)=…=(xM−1,yM−1)$), and corresponding isosceles trapezoids become rectangles. Hence, if (w, w)(0 < w < sin α) represents this uniform state, this type of solution can be represented or defined by following equation:
$f(w,w)=g(w,w)=w$
(6)
Hereafter, we call (w, w) satisfying Eq. (6)cylinder solution of crease pattern under parameters (α, β, N, M). Remarkably, cylinder solutions defined by Eq. (6) are called symmetric fixed points of the system (3) which is invariant under the map F and G.

Cylinder solutions shown in black in Fig. 7 are numerically calculated by solving Eq. (6) under crease-pattern parameter (α, β, N, M) = (45°, 34.6154°, 24, 100), which two solutions exist and each have different radius.

#### 3.1 Wave-Like Solution.

The second case corresponds to third to fifth from the top of waterbomb tube shown in Fig. 7, left. The wave-like folded state corresponds to the sequence of points on the phase space moving in a clockwise direction. We call such solution a wave-like solution. In a wave-like solution, the point continues rotating around the same closed curve without divergence or convergence. Note that the points along the closed curve does not coincide, thus the folded state is not periodic. However, the overall folded shape seems to approximate a smooth periodic wave-like surface. These sequences and corresponding cycles are nested around one of two cylinder solutions (smaller one).

In addition, these plots of nested closed curves seemed to be symmetric about the graph of ym = xm. The symmetry of the orbits is the result of the reversibility of the system. Generally, the orbit of reversible systems initiated from x0 = (x0, y0), that is defined as the set ${x=Fm(x0)|m∈Z}$, is invariant under the symmetry G when the orbit has a point which is invariant under G [21]. For this reason, the orbits shown in Fig. 7, which initial terms x0 are invariant under G, are actually symmetric about the graph of ym = xm.

#### 3.1 Finite Solution.

The second waterbomb tube from the top in Fig. 7, left, corresponds to the third type, where its points are plotted just a little outside of the above-mentioned concentric plots. The reason plots stop at the middle is that, at some index m, (xm, ym) deviate from the region (0, sin α) × (0, sinα), that is, there is no state of modules corresponding to parameter (xm, ym). In other words, finite solution appears in the case that the intersection of three spheres become empty at some step, when computing vertices of modules as shown in Fig. 6. Specifically, in the sequence corresponding to second waterbomb tube in Fig. 7, the numerical value of ninth term (x8, y8) is (0.689573, 0.724059), which y8 is greater than sin α = sin 45° ≈ 0.707107. So, these solution terminates at some point, and only a finite portion of the paper can be folded along this solution.

### 3.2 Kinematics.

From this visualization, we can observe a single disk region in the phase space as the set of initial values yielding solutions for any m (the gray region in Fig. 7). The region is the union of all wave-like solutions and the smaller cylinder solution. As the region is the configuration space of the mechanism with m → ∞, there exists a 2-DOF rigid folding motion. The motion can be represented by the “amplitude” and “phase” of wave-like surfaces. As the initial configuration gets closer to the cylinder solution, the amplitude of wave shapes gets smaller. We can change the phase by rotating the initial configuration along the closed curve (Supporting Movie).1 In the example state shown in Fig. 7, each initial value is taken along x0 = y0, so the left-most module forms the “valley” of the wave.

### 3.3 Outline of Subsequent Analysis.

From Fig. 7, we found that the system behaves differently around the different cylinder solutions. In Sec. 4, we perform a stability analysis to classify the symmetric fixed points, i.e., cylinder solutions, by the behavior of the system around them. We also investigate how the cylinder solutions emerge and disappear and how the stability changes when the crease-pattern parameters are changed.

We conjecture that quasi-periodic solutions exist around a neutrally stable cylinder solution, which explains the nested closed solutions representing wave-like solutions. In Sec. 5, we give a proof for the conjecture for a particular crease pattern.

## 4 Stability and Bifurcation Analysis

In this section, we investigate how the behavior of the system around cylinder solutions changes when the crease-pattern parameters are changed. For this, we numerically analyze the stability of the cylinder solutions and visualize the changes of the stability using the bifurcation diagram. Figure 8 shows the concept of the analysis in this section. From the bifurcation diagram, we found that there are some critical values of crease-pattern parameters where the stability or the behavior of the system changes. Note that the result of this section is obtained numerically (not symbolically) using Mathematica based on the explicit form of system (3).

### 4.1 Stability Analysis.

First, we introduce the concept of stability and classify cylinder solutions. In order to know the stability of the cylinder solutions, we perform the linear stability analysis. The linear stability analysis is the stability analysis on the following linearized system at the cylinder solutions w = (w, w):
$xm+1−w≈DF(w)(xm−w)$
(7)
where DF(w) is the Jacobian matrix of the system at the cylinder solutions defined by
$DF(x,y)=[∂xf(x,y)∂yf(x,y)∂xg(x,y)∂yg(x,y)]$
(8)
The linearized behavior of the nonlinear system around the fixed point can be characterized by using the eigenvalues of the Jacobian matrix at that point. In a typical case of the system (3), we found that the eigenvalues of the Jacobian matrix at the symmetric fixed point (w, w) can be classified into two types. In the first type, the eigenvalues are complex conjugate whose magnitude equals to 1 at w. Such a point is called elliptic fixed point, where the linear system (7) has concentric closed elliptical orbits around the origin. The linear stability around elliptic fixed point is neutrally stable, meaning that the orbit around the fixed point keeps some distance from that point called center. In the second type, the eigenvalues are real, and one absolute value is less than 1 and the other is greater than 1. Such a point is called a hyperbolic fixed point. The linear stability at a hyperbolic fixed point is unstable, and the fixed point is called a saddle, i.e., meaning that it has both stable and unstable directions.2

In Fig. 8, the top-left figure shows the two different cylinder configurations with larger and smaller radius under (α, β, N) = (45°, 40°, 8) represented by the numerically calculated cylinder solutions (w1, w1) and (w2, w2), respectively. The linear stability of (w1, w1) and (w2, w2) are unstable saddle and neutrally stable center, respectively.

Fig. 8
Fig. 8
Close modal

A fixed point classified as an unstable saddle point in the linearized system (7) is also an unstable saddle point in the original nonlinear system (see Sec. 5). However, a fixed point classified as a neutrally stable fixed point in the linearized system is not necessarily neutrally stable in the original system. Therefore, the linear stability analysis alone does not guarantee the existence of nested closed solutions. Nevertheless, from the phase diagram in Fig. 8, we found that there are nested closed plots around (w2, w2). This leads us to conjecture that if the linearized system (7) has a neutrally stable symmetric fixed point, it is also a neutrally stable symmetric fixed point in the original system, around which there are quasi-periodic solutions. We show the conjecture is correct for a particular crease pattern in Sec. 5.

### 4.2 Bifurcation Analysis.

Next, we analyze how the linear stability of cylinder solutions changes when β changes under fixed α and N. For the visualization of the result, we use a bifurcation diagram. The bottom of Fig. 8 shows the bifurcation diagram for α = 45°, N = 8.

To create bifurcation diagram for given α and N, we move the remaining parameter β in the range (0, 180° − α), and compute the cylinder solutions (w, w) of the crease pattern by solving Eq. (6) at each value. Here, we transform parameter w to $Θ≡arcsin(w/(sinα))∈(0,90∘)$, the half of the dihedral angle formed by the two isosceles triangles (Fig. 8), so that the plot range is independent of α. Then, we perform the linear stability analysis and plot the cylinder solutions in β–Θ planes by solid and dotted line if it is neutrally stable and unstable, respectively. In Fig. 8, the vertical line β = 40° intersects with the bifurcation diagram at two different points, (40°, Θ1) on the dotted line and (40°, Θ2) on the solid line, shown in red and blue, respectively. The red and the blue points correspond to the unstable cylinder with a larger radius and the neutrally stable cylinder with a smaller radius, respectively.

### 4.3 Result.

Now, we describe how the linear stability of the cylinder solutions depends on the parameters using the bifurcation diagram. First, we observe the bifurcation diagram for fixed N = 8. By varying α, we found the critical value α* ≈ 60.4° at which the appearance of the bifurcation diagram changes. So, we use the representative cases α = 45° < α* and α = 65° > α* for further observing the bifurcation diagram for these cases. We found some critical values of β where the number of cylinder solutions or their linear stability change, i.e., whether waterbomb tube can be wave-like surface changes. We also explain the differences between the two cases. Finally, we show that we can apply the observations of the case of N = 8 for the other N.

Here, we fix N = 8. Under fixed N, we can create the 3D bifurcation diagram by varying the value α in range (0, 90°) and arranging the 2D bifurcation diagrams for each α (right in Fig. 9). From the 3D diagram, we found the critical value α* ≈ 60.4° where the behavior of the bifurcation diagram changes. The behaviors of the bifurcation diagram for α < α* are represented by (α, N) = (45°, N = 8) (blue bifurcation diagram on the top-left in Fig. 9), while the bifurcation diagram for α > α* is represented by (α, N) = (65°, N = 8) (red bifurcation diagram on the bottom-left in Fig. 9).

Fig. 9
Fig. 9
Close modal

First, we describe the behavior of the representative case of (α, N) = (45°, 8) described in the blue bifurcation diagram on the top-left in Fig. 9. As we change the remaining variable β, there are some critical value of β denoted by β*1, β*2, where the linear stability of a cylinder solution changes. As we increase β from 0, at β = β*1 ≈ 35.8°, the graph of y = f(w, w) and y = w touches and the fixed point that is saddle and center generated (saddle-node bifurcation). The branch of the fixed point that is saddle extends from the bottom of the graph at β = 45°. Then, at β = β*2 ≈ 46.1°, saddle-node bifurcation occurs again where the center and the saddle coincide and disappear. Thus, the system (3) has no cylinder solutions under β ∈ (0, β*1), two cylinder solutions (center and saddle) under β ∈ (β*1, 45°), three cylinder solutions (center and two saddle) under β ∈ (45°, β*2), and only one cylinder solution (saddle) under β ∈ (β*2, 135°). Here, Fig. 10 shows the part of the bifurcation diagram and the phase diagrams for β = 35°, 40°, 45.5°, and 50° belonging to (0, β*1), (β*1, 45°), (45°, β*2), and (β*2, 135°), respectively. Under β = 40° ∈ (β*1, 45°) and β = 45.5° ∈ (45°, β2°), there are nested cyclic plots around the center, which suggests the existence of a wave-like configurations. On the other hand, there are no cyclic plots in the other two diagrams, i.e., all solutions are finite solutions; therefore, waterbomb tube under $α=45∘,N=8,β∈(0,β1*)⋃(β2*,135∘)$ cannot be wave-like surfaces. Thus, whether waterbomb tube can be wave-like surface or not changes at the bifurcation values.

Fig. 10
Fig. 10
Close modal

Next, for the case of (α, N) = (65°, 8), the red diagram in Fig. 9 shows that there are three critical values for β: β*1, β*2, and β*. As we increase β from 0, saddle-node bifurcation occurs at β = β*1 ≈ 27.8° in the same way as the α = 45°; however, as β is further increased, the linear stability changes at β = β* ≈ 52.3°, which is not observed in the case of α < α*. Furthermore, the branch of the center fixed point extends from the bottom of the graph at β = 65°. Finally, at β = β*2 ≈ 65.5°, saddle-node bifurcation occurs again where the saddle and the center coincide and disappear. From the above, in the case of α = 65°, the system (3) has zero, two (center and saddle), two (two saddles), three (two saddles and center), and one (saddle) cylinder solutions when β ∈ (0, β*1), (β*1, β*), (β*, 65°), (65°, β*2), and (β*2, 115°), respectively. Figure 11 shows the part of the bifurcation diagram and the phase diagrams for β = 25°, 40°, 49°, 57°, 65.1°, and 70° belonging to (0, β*1), (β*1, β*), (β*1, β*), (β*, 65°), (65°, β*2), and (β*2, 115°), respectively. In the top-right, middle-left, and bottom-left phase diagrams, there are the plots of wave-like solutions around the center; therefore, waterbomb tube has wave-like configurations. However, the wave-like configurations placed in each diagram self-intersect at some parts in the same way as the configuration in the bottom on the left side of Fig. 11, or this is caused by too small radius as for β = 65.1°, which means these wave-like configurations are not realizable. As for the middle-left diagram, we found the unstable nonsymmetric fixed points (x, y) ≈ (0.698488, 0.375731), (x, y) ≈ (0.375731, 0.698488) satisfying f(x, y) = g(x, y) = (x, y), which is not observed in the case of α = 45°, but the configurations of the module and the tube corresponding to the nonsymmetric fixed points are self-intersecting complicatedly as shown in Fig. 11. The other three diagrams, top-left, middle-right, and bottom-right on the right side of Fig. 11, having no wave-like solutions, which suggests that waterbomb tube cannot be wave-like surfaces if $β∈(0,β1*)⋃(β*,65∘)⋃(β2*,115∘)$. Hence, bifurcation values are closely related to the possible forms of waterbomb tubes as in the case of α = 45°.

Fig. 11
Fig. 11
Close modal

For cases other than N = 8, the top, middle, and bottom figure in Fig. 12 shows the 3D bifurcation diagrams for N = 6, N = 30, and N = 100, respectively. Comparing 3D diagrams under the different values of N in Figs. 9 and 12, the shape of the diagrams and the critical values α* depends on the values of N. However, the pattern formed by the dotted and solid lines is the same for the four different values of N; therefore, we can apply the features described for N = 8, i.e, the existence of critical values such as α*, β*1, β2*, β* and their relation with the behavior of waterbomb tube for N = 6, 30, and 100.

Fig. 12
Fig. 12
Close modal

## 5 Proof of Wave-Like Solution

In this section, we claim that there are nested wave-like solution around the neutrally stable cylinder solution and they are quasi-periodic orbits of the system; and we give the proof of existence for fixed crease-pattern parameters (α, β, N) = (45°, 45°, 6). We also show the system behaves differently around the saddle cylinder solution. For the proof of quasi-periodic orbits, we use the Kolmogorov–Arnold–Moser (KAM) theorem that guarantees the equivalence of the behaviors in the nonlinear and linearized systems under some conditions as described later.

Figure 13 visualizes the solutions under (α, β, N) = (45°, 45°, 6). The function f under this crease pattern is given as follows:
$f(x,y)=14(y2−1)((x−y)2−2)×(x(3((2y2−1)(x2+y2−1)+1)+2y(2y2−1)(2(x−y)2−1))−y(2y((2y2−1)(2(x−y)2−1)−2(2(x−y)2−1)(x2+y2−1))+3((2y2−1)(x2+y2−1)−3)+33y2)+2((2y2−1)(2(x−y)2−1)−(2(x−y)2−1)(x2+y2−1))+(−3)x2y)$
(9)
We derived the symbolic forms of cylinder solutions, Jacobian matrix, and its eigenvalues for linear stability analysis by using Mathematica. We obtain two cylinder solutions analytically as the following by solving Eq. (6) for w (0 < w < sin 45°):
$w=1410−33±−1+43$
(10)
Hereafter, let w1 be the solution with a larger radius, and w2 be another one with a smaller radius.
Fig. 13
Fig. 13
Close modal
Next, we derive the Jacobian matrix DF(x, y) at cylinder solutions (x, y) = (wi, wi) (i = 1, 2). Note that ∂xg(wi, wi) and ∂yg(wi, wi) can be computed from symbolic forms of ∂xf(wi, wi), ∂yf(wi, wi) using chain rule:
$∂g∂x(wi,wi)=−∂f∂y(wi,wi)$
(11)
$∂g∂y(wi,wi)=1−(∂f∂y(wi,wi))2∂f∂x(wi,wi)$
(12)
The elements of the Jacobian matrices DF(wi, wi) (i = 1, 2) are computed as
$(DF(w1,w1))1,1=26243+3843+609+47(DF(w1,w1))1,2=18(−153−51123−8845+28)(DF(w1,w1))2,1=18(153+51123−8845−28)(DF(w1,w1))2,2=1104(5753+403884323−69922045−688)$
(13)
and
$(DF(w2,w2))1,1=18(−33+963−165+8)(DF(w2,w2))1,2=18(−153+51123−8845+28)$
(14)
$(DF(w2,w2))2,1=18(153−51123−8845−28)(DF(w2,w2))2,2=1104(5753−403884323−69922045−688)$
(14)
Finally, we consider eigenvalues of the derived Jacobian matrices (13), (14) and the linear stability of corresponding cylinder solutions. Table 1 shows the numerical values of the eigenvalues and the linear stability of the cylinder solutions. The eigenvalues for both cylinder solutions are given as the roots of the following polynomial:
$13w8+292w7+161868w6−1083460w5+1704206w4−1083460w3+161868w2+292w+13$
(15)
The eigenvalues for larger (x, y) = (w1, w1) are real, the absolute values of which are greater than 1 and smaller than 1, respectively; therefore, the type of cylinder solution (w1, w1) as the fixed point of the system is saddle. Hence, both linear and nonlinear stabilities of (w1, w1) are unstable, so there are no cyclic solutions around (w1, w1).
Table 1

Eigenvalues and stability at two cylinder solutions

CYLINDER SOLUTIONEIGENVALUETYPESTABILITY
$w1=1410−33+−1+43$$|λ1|=|4.69…|>1$
$|λ2|=|0.212…|<1$
$(λ1,λ2∈R)$
$w2=1410−33−−1+43$$|λ|=|0.855…±i0.517…|=1$
$(λ∈C)$
CenterNeutrally
Stable
CYLINDER SOLUTIONEIGENVALUETYPESTABILITY
$w1=1410−33+−1+43$$|λ1|=|4.69…|>1$
$|λ2|=|0.212…|<1$
$(λ1,λ2∈R)$
$w2=1410−33−−1+43$$|λ|=|0.855…±i0.517…|=1$
$(λ∈C)$
CenterNeutrally
Stable

On the other hand, at the smaller cylinder solution (x, y) = (w2, w2) with a smaller radius, DF(w2, w2) has complex conjugate eigenvalues whose magnitudes are exactly 1. Hence, the type of the cylinder solution (w2, w2) is center and linear stability of (w2, w2) is neutrally stable. Nevertheless, because (w2, w2) is an elliptic fixed point, it is not yet guaranteed that the original nonlinear system (3) has nested closed orbits around (w2, w2) as we mentioned in Sec. 4.1.

Now, to prove the existence of the nested closed symmetric quasi-periodic orbits around (w2, w2) in the original system, we use the one derived form of KAM theorem. The KAM theorem guarantees the existence of such orbits called KAM curve around a fixed point of the reversible dynamical system if the fixed point is elliptic symmetric fixed point and nonresonant, that is, the eigenvalues are not being roots of unity [21,22]. We already know that the system is reversible and (w2, w2) is an elliptic symmetric fixed point. Also, (w2, w2) is nonresonant because of the following reason. The eigenvalues of DF(w2, w2) are the roots of the irreducible polynomial (15). Then, some coefficients of the minimal polynomial, obtained as the monic form of (15), are not integers (e.g., coefficients for seventh-order term is 292/13), so the eigenvalues of (w2, w2) are not algebraic integer, which implies that the eigenvalues of (w2, w2) are not the roots of unity. Thus, we can apply KAM theorem to the system at (w2, w2), and the existence of wave-like solution around the point is proved.

## 6 Conclusion

In this paper, we revealed a part of the mathematical structure behind wave-like surfaces of waterbomb tube. First, we have described the kinematic model representing the configuration of each module of waterbomb tube and its entire shape to derive recurrence relation of two variables dominating its kinematics. Then, based on the observation of the plots on the phase space corresponding to the folded states, we classified solutions into three types to identify the wave-like solutions surrounding one of two cylinder solutions. We applied the knowledge of discrete dynamical systems to the recurrence relation, and investigated how the stability of the symmetric fixed point of the system, i.e., the behavior around cylinder solutions, changes when we change the crease-pattern parameters. Then, we observed that around the neutrally stable cylinder solution, there exists nested wave-like solutions, and found critical values where the behavior of the system changes. Finally, we gave the proof of the existence of quasi-periodic solutions, i.e., wave-like solutions, around the neutrally stable cylinder solution.

Note that we gave proof of the existence of quasi-periodic solutions only for fixed crease-pattern parameters; characterizing the existence of such solutions is a future work of this paper. Also, the physical interpretation, i.e., the relation between the mathematical properties of solutions and the mechanical properties are still left to be investigated.

More generally, our proposed approach based on the geometry of a module and applying the knowledge of dynamical system to the recurrence relation can be useful to explain behaviors of various origami tessellations. The analysis of tessellations consisting of modules with higher or lower DOF is of our interest.

## Footnotes

2

Here, stable or unstable refers to the stability of the dynamical systems as m increases. It does not refer to mechanical stability.

## Acknowledgment

This work is supported by JST PRESTO (Grant No. JPMJPR1927).

### Appendix A: Derivation of Equation Imposing Cylinder Constraint

Here, we derive Eq. (6) imposing cylinder constraint. Assuming that values of crease-pattern parameters are given, we represent some distance of vertices or angle using (α, β, N, M) based on the preservation of edge length.

We derive some quantities for considering the cylinder constraint. First, edge length l1, l2 in Fig. 14 can be written in the following form:
$l1=sinαsin(α+β),l2=cosα−sinαtan(α+β)$
Next, we consider the dihedral angle ρm. The notation of vertices is defined as Fig. 15. Because of the N-fold symmetry of waterbomb tube, vertices $Am,01,Am,11,…,Am,N−11$ form regular N-sided polygon. Based on this, we derive Eq. (6). Let θm,n, Mm,n, Pm,n, Qm,n be the base angle $∠Am,n1Bm,n1Bm,n2=∠Am,n2Bm,n1Bm,n1$, a midpoint of $Am,n1Bm,n1$, a foot of perpendicular from point Mm,n to edge $Am,n1Am,n2$, and a foot of perpendicular from point Mm,n to edge $Bm,n1Bm,n2$, respectively. Here, following equations hold:
$cosθm,n=xm−ym2cosα|Mm,nPm,n|=|Mm,nQm,n|=xmsinθm,n$
Because $∠Mm,n+1Am,n1Mm,n=π(N−2)/N$, by using cosine theorem in $△Mm,n+1Am,n1Mm,n$,
$|Mm,nMm,n+1|=2xm2−2xm2cosπ(N−2)N$
Thus, by using cosine theorem in $△Mm,n+1Pm,nMm,n$,
$cos2ρm=cos∠Mm,n+1Pm,nMm,n=2|Mm,nPm,n|2−|Mm,nMm,n+1|22|Mm,nPm,n|2$
By using half angle formula, we can get the following equation equivalent to Eq. (6):
$cosρm=−(xm−ym)2sec2α+4sin2πN4−(xm−ym)2sec2α$
Fig. 14
Fig. 14
Close modal
Fig. 15
Fig. 15
Close modal

### Appendix B: Derivation of Function f

Here, we derive the function f and g assuming that values of crease-pattern parameters and the state of modules belonging to mth hoop, i.e., xm and ym, are given. First, we represent some quantities such as distance between vertices or angles by crease-pattern parameters and xm and ym. Next, we derive the function f and g by introducing the coordinate system and representing xm+1 and ym+1 by the coordinates of some vertices of waterbomb tube.

#### Some Necessary Quantities

First, we represent some quantities using (α, β, N, M) and (xm, ym) based on the preservation of edge length. The notation of vertices is defined in Fig. 16, and for simplicity, we omit subscripts of xm and ym.

Fig. 16
Fig. 16
Close modal
Let θ1, θ2 be the base angle $∠B1A1A2=∠B2B1A1$, $∠A2B2B1=∠A1A2B2$, respectively (0 < θi < π). Their cosine and sine can be represented as follows:
$cosθ1=x−y|A1A2|=x−y2cosα,sinθ1=1−cos2θ1$
(A1)
$cosθ2=y−x|A1A2|=y−x2cosα,sinθ2=1−cos2θ2$
(A2)
Next, we calculate h : = |OO′|, where O′ corresponds to the circumcenter of isosceles trapezoid A1B1B2A2 because point O is equidistant from Ai, Bi (i = 1, 2); therefore, if we represent the length of diagonals of A1B1B2A2, radius of circumcircle of A1B1B2A2 as d, r, respectively, d, r, h can be expressed as follows:
$d=|A2B1|=|B1A1|2+|A1A2|2−2|B1A1‖A1A2|cosθ1=2xy+cos2αr=|A1O′|=|A2B1|2sin∠B1A1A2=d2sinθ1h=|OO′|=|A1O′|2−|A1O′|2=1−r2$
On the other hand, if $ψi=∠OMiCi$, $ζi=∠OMiO′$, $ξi=ψi−ζi(|ξi|=∠RiMiCi)$ (i = 1, 2), the following equations hold:
$cosψi=|MiCi|2+|OMi|2−|CiO|22|MiCi‖OMi|={(l12−x2)+(1−x2)−l222l12−x21−x2(i=1)(l12−y2)+(1−y2)−l222l12−y21−y2(i=2)sinψi=1−cos2ψisinζi=|OO′||OMi|={h1−x2(i=1)h1−y2(i=2)cosζi=1−cos2ζi$
$cosξi=cosψicosζi+sinψisinζisinξi=sinψicosζi−cosψisinζi$
(A3)

#### Derivation of xm+1 and ym+1

Next, we derive xm+1 and ym+1, i.e., the functions f and g. Here, we take the coordinate system shown in Fig. 17 and represent the coordinates of each vertex as the functions of xm and ym. Then, we obtain xm+1 and ym+1, i.e., the functions f and g, by using the coordinates of vertices. In Figs. 17 and 18, the coordinate system is taken so that the rotation axis of symmetry of the waterbomb tube is the X-axis and one of the modules belonging to the mth hoop which vertices are denoted as □m,0 is mirror-symmetric for the XZ-plane. In the following, the coordinate components of vertex V are denoted as Vx, Vy, and Vz, respectively, and Ry(η) and t denote the rotation matrix of angle η with Y-axis as the rotation axis and translation vector defined as t ≡ [0, 0, xm/tan (π/N)]T. Here, we can represent the coordinates of each vertex as follows using Eqs. (A1), (A3):
$[Am,01x,Am,01y,Am,01z]T=Ry(η)[0,−xm,0]T+t[Bm,01x,Bm,01y,Bm,01z]T=Ry(η)[0,xm,0]T+t$
$[Am,02x,Am,02y,Am,02z]T=Ry(η)[2cosαsinθm,01,−ym,0]T+t[Bm,02x,Bm,02y,Bm,02z]T=Ry(η)[2cosαsinθm,01,ym,0]T+t[Om,0x,Om,0y,Om,0z]T=Ry(η)[rm,02−xm2,0,−hm,0]T+t[Cm,01x,Cm,01y,Cm,01z]T=Ry(η)[l12−xm2cosξm,01,0,l12−xm2sinξm,01]T+t[Cm,02x,Cm,02y,Cm,02z]T=Ry(η)[l12−ym2cosξm,02,0,l12−ym2sinξm,02]T+t$
Note that for the elements of the matrix Ry(η) we have
$sinη=−Am,02z−Am,01z2cosαsinθm,01,cosη=1−sin2η$
Using the coordinates of the vertices, we can obtain the coordinate of the vertices of the other modules belonging to the mth hoop by rotating the corresponding vertex by 2/N around the X-axis. For example,
$[Cm,12x,Cm,12y,Cm,12z]T=Rx(2π/N)[Cm,02x,Cm,02y,Cm,02z]T$
Here, $xm+1=[Cm,02x,Cm,02y,Cm,02z]T⋅[0,cosπN,sinπN]T$, which gives the function f (Fig. 6).
Fig. 17
Fig. 17
Close modal
Fig. 18
Fig. 18
Close modal
Next, we derive the function g by finding ym+1 (Fig. 6). To do so, we first find the coordinates of vertex Om+1,1, and then the coordinates of vertex $Bm+1,12$. From Fig. 6, vertex Om+1,1 is located at a distance 1, 1, l2 from the vertices $Cm,02,Cm,12,Cm+1,12$, respectively. Vertices $Cm,02$ and $Cm,12$ are mirror-symmetric across plane $P$ obtained by rotating the XZ-plane by π/N around X-axis. Hence, vertex Om+1,1 is contained in the intersection of the two circles on $P$ centered at P and $Cm+1,11$ which radius are $|POm+1,11|=1−|PCm,0|2=1−xm+12$ and l2, respectively. Using x1 : = [Px, Py, Pz]T, $x2:=[Cm+1,11xCm+1,11y,Cm+1,11z]T,$$R1:=1−xm+12,R2:=l2$, and n : = Rx(π/N)[0, 1, 0]T which is the normal vector of $P$, we can represent the intersection of the two circles as follows:
$x1+R(n,±arccos(R12+‖x1−x2‖2−R222R1‖x1−x2‖))R1(x2−x1)‖x2−x1‖$
(A4)
Here, $R(v,θ)(v∈R3,θ∈R)$ is the rotation matrix in which rotation axis and rotation angle are v and θ, respectively. The signs +, − of the angle of rotation results crease $Cm+1,11Om+1,1$ to mountain-folded and valley-folded, respectively; therefore, in this case, + sign is consistent with the prescribed MV-assignment.

Next, we derive the coordinates of vertex $Bm+1,12$. Vertex $Bm+1,12$ is located at a distance 1, 1, 2 cos α from vertices $Om+1,0,Om+1,1,Cm,02$, respectively. Here, since the vertices Om+1,0, Om+1,1 are the mirror-reflection with respect to XZ-plane, the vertex $Bm+1,12$ is contained in the intersection of a circle of radius $1−(Om+1,1y)2$ centered at the vertex Q on the XZ-plane and a circle of radius 2cosα centered at the vertex $Cm,02$. We can obtain the intersection of the circles by Eq. (A4) where $x1=[Cm,02x,Cm,02y,Cm,02z]T,R1=2cosα,x2=[Qx,Qy,Qz]T,R2=$$1−(Om+1,1y)2,n=[0,1,0]T$. In this case, the positive sign is consistent with the prescribed MV-assignment of $Cm,02Bm+1,12$. Therefore, we obtain the coordinates of vertex $Bm+12$. Since $ym+1=[Bm+12x,Bm+12y,Bm+12z]T⋅[0,cosπN,sinπN]T$, the function g is obtained.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

## References

1.
Callens
,
S. J.
, and
,
A. A.
,
2018
, “
From Flat Sheets to Curved Geometries: Origami and Kirigami Approaches
,”
Mater. Today
,
21
(
3
), pp.
241
264
.
2.
Schenk
,
M.
, and
Guest
,
S. D.
,
2011
, “
Origami Folding: A Structural Engineering Approach
,”
The 5th International Meeting on Origami in Science, Mathematics and Education (5OSME)
,
Singapore
,
July
.
3.
Nassar
,
H.
,
Lebée
,
A.
, and
Monasse
,
L.
,
2017
, “
Curvature, Metric and Parametrization of Origami Tessellations: Theory and Application to the Eggbox Pattern
,”
Proc. R. Soc. A: Math. Phys. Eng. Sci.
,
473
(
2197
), p.
20160705
.
4.
Nassar
,
H.
,
Lebée
,
A.
, and
Monasse
,
L.
,
2018
, “
Fitting Surfaces With the Miura Tessellation
,”
The 7th International Meeting on Origami in Science, Mathematics and Education (7OSME)
,
Oxford, UK
,
Sept. 5–7
.
5.
Fujimoto
,
S.
,
1976
,
Sōzō sei wo kaihatsu suru rittai Origami
,
Shin-Shashoku-Shuppan
,
Osaka, Japan
.
6.
Kuribayashi
,
K.
,
Tsuchiya
,
K.
,
You
,
Z.
,
Tomus
,
D.
,
Umemoto
,
M.
,
Ito
,
T.
, and
Sasaki
,
M.
,
2006
, “
Self-deployable Origami Stent Grafts As a Biomedical Application of Ni-rich Tini Shape Memory Alloy Foil
,”
Mater. Sci. Eng. A
,
419
(
1–2
), pp.
131
137
.
7.
Fang
,
H.
,
Zhang
,
Y.
, and
Wang
,
K.
,
2017
, “
Origami-Based Earthworm-Like Locomotion Robots
,”
Bioinspiration Biomimetics
,
12
(
6
), p.
065003
.
8.
Li
,
S.
,
Stampfli
,
J. J.
,
Xu
,
H. J.
,
Malkin
,
E.
,
Diaz
,
E. V.
,
Rus
,
D.
, and
Wood
,
R. J.
,
2019
, “
A Vacuum-Driven Origami ‘Magic-Ball’ Soft Gripper
,”
2019 International Conference on Robotics and Automation (ICRA)
,
,
May 20–24
, IEEE, pp.
7401
7408
.
9.
Lee
,
D.-Y.
,
Kim
,
J.-K.
,
Sohn
,
C.-Y.
,
Heo
,
J.-M.
, and
Cho
,
K.-J.
,
2021
, “
,”
Sci. Robot.
,
6
(
53
).
10.
Rodrigues
,
G. V.
,
Fonseca
,
L. M.
,
Savi
,
M. A.
, and
Paiva
,
A.
,
2017
, “
Nonlinear Dynamics of an Adaptive Origami-Stent System
,”
Int. J. Mech. Sci.
,
133
, pp.
303
318
.
11.
Fonseca
,
L. M.
,
Rodrigues
,
G. V.
,
Savi
,
M. A.
, and
Paiva
,
A.
,
2019
, “
Nonlinear Dynamics of an Origami Wheel With Shape Memory Alloy Actuators
,”
Chaos Solitons Fractals
,
122
, pp.
245
261
.
12.
Feng
,
H.
,
Ma
,
J.
,
Chen
,
Y.
, and
You
,
Z.
,
2018
, “
Twist of Tubular Mechanical Metamaterials Based on Waterbomb Origami
,”
Sci. Rep.
,
8
(
1
), pp.
1
13
.
13.
Ma
,
J.
,
Feng
,
H.
,
Chen
,
Y.
,
Hou
,
D.
, and
You
,
Z.
,
2020
, “
Folding of Tubular Waterbomb
,”
Research
,
2020
, p.
1735081
.
14.
Fonseca
,
L. M.
, and
Savi
,
M. A.
,
2021
, “
On the Symmetries of the Origami Waterbomb Pattern: Kinematics and Mechanical Investigations
,”
Meccanica
,
56
(
10
), pp.
2575
2598
.
15.
Feng
,
H.
,
2018
,
Kinematics of Spatial Linkages and its Applications to Rigid Origami
,
PhD thesis
,
Université Clermont Auvergne
,
Clermont-Ferrand
.
16.
,
T.
,
Ma
,
J.
,
Feng
,
H.
,
Hou
,
D.
,
Gattas
,
J. M.
,
Chen
,
Y.
, and
You
,
Z.
,
2020
, “
Programmable Stiffness and Shape Modulation in Origami Materials: Emergence of a Distant Actuation Feature
,”
Appl. Mater. Today
,
19
, p.
100537
.
17.
Chen
,
Y.
,
Feng
,
H.
,
Ma
,
J.
,
Peng
,
R.
, and
You
,
Z.
,
2016
, “
Symmetric Waterbomb Origami
,”
Proc. R. Soc. A: Math. Phys. Eng. Sci.
,
472
(
2190
), p.
20150846
.
18.
Kawasaki
,
T.
,
1994
, “
R (γ)= 1
,”
The 2nd International Meeting of Origami Science and Scientific Origami (2OSME)
,
Otsu, Japan
,
November–December
, Vol. 3, pp.
31
40
.
19.
Belcastro
,
S.-M.
, and
Hull
,
T. C.
,
2001
, “
A Mathematical Model for Non-flat Origami
,”
The 3rd International Meeting of Origami Mathematics, Science, and Education (3OSME)
,
Monterey, CA
,
Mar. 9–11
, pp.
39
51
.
20.
Abel
,
Z.
,
Cantarella
,
J.
,
Demaine
,
E. D.
,
Eppstein
,
D.
,
Hull
,
T. C.
,
Ku
,
J. S.
,
Lang
,
R. J.
, and
Tachi
,
T.
,
2016
, “
Rigid Origami Vertices: Conditions and Forcing Sets
,”
J. Comput. Geom.
,
7
(
1
), pp.
171
184
.
21.
Roberts
,
J. A.
, and
Quispel
,
G.
,
1992
, “
Chaos and Time-Reversal Symmetry. Order and Chaos in Reversible Dynamical Systems
,”
Phys. Rep.
,
216
(
2–3
), pp.
63
177
.
22.
Sevryuk
,
M. B.
,
1986
,
Reversible Systems
,
Springer
,
Berliln/Heidelberg
.