## 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 [12–14]. 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.

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\u2208Z>0,M\u2208Z>0$).

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.

### 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.

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.

*m*th hoop of the waterbomb tube. Because of rotational symmetry, the configuration of the

*m*th hoop can be represented by using parameter

*x*

_{m},

*y*

_{m}for each congruent module and the common dihedral angles between isosceles trapezoids 2

*ρ*

_{m}(Fig. 4, right). In order that the

*m*th hoop to be a closed cylinder, we have the following constraint (see Appendix A for derivation).

*ρ*

_{m}can be derived as the function of

*x*

_{m}and

*y*

_{m}, the configuration of each hoop can be uniquely represented by parameters

*x*

_{m},

*y*

_{m}. However, note that not all values of (

*x*

_{m},

*y*

_{m}) represent valid states as cos

*ρ*

_{m}may become nonreal number for certain given pair of values.

*x*

_{m},

*y*

_{m}and

*x*

_{m+1},

*y*

_{m+1}that need to be satisfied. The constraints result in the significant property of waterbomb tube that the state of

*m*th 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

*V*

_{1},

*V*

_{2},

*V*

_{3}are fixed by constructing three spheres

*S*

_{1},

*S*

_{2},

*S*

_{3}of radii

*XV*

_{1},

*XV*

_{2},

*XV*

_{3}centering at

*V*

_{1},

*V*

_{2},

*V*

_{3}, 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

*V*

_{1},

*V*

_{2}, and

*V*

_{3}. Thus, MV-assignments of creases

*XV*

_{3}and

*X*′

*V*

_{3}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

*m*th module. Note that the existence of a solution is not guaranteed because intersections of three spheres may be empty.

*m*th hoop and

*m*+ 1th hoop can be formulated by following recurrence relation:

*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

*x*

_{0},

*y*

_{0}(0 <

*x*

_{0},

*y*

_{0}< 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 *x*_{m}, *y*_{m}. Also, the system (3) is *autonomous system* because the function *f*, *g* is independent of indices *m*.

*reversibility*. Because the module of waterbomb tessellation is point symmetric about the central vertex (see Fig. 2, left), the following equation holds $(m=1,\u2026,M\u22121)$:

*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 (*x*_{0}, *y*_{0}) 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\u2208Z}$ is called *orbit* of (*x*_{0}, *y*_{0}) 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.

### 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.

*w*,

*w*)(0 <

*w*< sin

*α*) represents this uniform state, this type of solution can be represented or defined by following equation:

*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**.

#### 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 *y*_{m} = *x*_{m}. The symmetry of the orbits is the result of the reversibility of the system. Generally, the orbit of reversible systems initiated from **x**_{0} = (*x*_{0}, *y*_{0}), that is defined as the set ${x=Fm(x0)|m\u2208Z}$, 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 **x**_{0} are invariant under **G**, are actually symmetric about the graph of *y*_{m} = *x*_{m}.

#### 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*, (*x*_{m}, *y*_{m}) deviate from the region (0, sin *α*) × (0, sin*α*), that is, there is no state of modules corresponding to parameter (*x*_{m}, *y*_{m}). 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 (*x*_{8}, *y*_{8}) is (0.689573, 0.724059), which *y*_{8} 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 *x*_{0} = *y*_{0}, 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.

*linear stability analysis*. The linear stability analysis is the stability analysis on the following linearized system at the cylinder solutions

**w**= (

*w*,

*w*):

*D*

**F**(

**w**) is the Jacobian matrix of the system at the cylinder solutions defined by

*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 (*w*_{1}, *w*_{1}) and (*w*_{2}, *w*_{2}), respectively. The linear stability of (*w*_{1}, *w*_{1}) and (*w*_{2}, *w*_{2}) are unstable saddle and neutrally stable center, respectively.

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 (*w*_{2}, *w*_{2}). 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 $\Theta \u2261arcsin(w/(sin\alpha ))\u2208(0,90\u2218)$, 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).

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 $\alpha =45\u2218,N=8,\beta \u2208(0,\beta 1*)\u22c3(\beta 2*,135\u2218)$ cannot be wave-like surfaces. Thus, whether waterbomb tube can be wave-like surface or not changes at the bifurcation values.

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 $\beta \u2208(0,\beta 1*)\u22c3(\beta *,65\u2218)\u22c3(\beta 2*,115\u2218)$. Hence, bifurcation values are closely related to the possible forms of waterbomb tubes as in the case of *α* = 45°.

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.

## 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.

*α*,

*β*,

*N*) = (45°, 45°, 6). The function

*f*under this crease pattern is given as follows:

*w*(0 <

*w*< sin 45°):

*w*

_{1}be the solution with a larger radius, and

*w*

_{2}be another one with a smaller radius.

*D*

**F**(

*x*,

*y*) at cylinder solutions (

*x*,

*y*) = (

*w*

_{i},

*w*

_{i}) (

*i*= 1, 2). Note that ∂

_{x}

*g*(

*w*

_{i},

*w*

_{i}) and ∂

_{y}

*g*(

*w*

_{i},

*w*

_{i}) can be computed from symbolic forms of ∂

_{x}

*f*(

*w*

_{i},

*w*

_{i}), ∂

_{y}

*f*(

*w*

_{i},

*w*

_{i}) using chain rule:

*D*

**F**(

*w*

_{i},

*w*

_{i}) (

*i*= 1, 2) are computed as

*x*,

*y*) = (

*w*

_{1},

*w*

_{1}) are real, the absolute values of which are greater than 1 and smaller than 1, respectively; therefore, the type of cylinder solution (

*w*

_{1},

*w*

_{1}) as the fixed point of the system is saddle. Hence, both linear and nonlinear stabilities of (

*w*

_{1},

*w*

_{1}) are unstable, so there are no cyclic solutions around (

*w*

_{1},

*w*

_{1}).

CYLINDER SOLUTION | EIGENVALUE | TYPE | STABILITY | |
---|---|---|---|---|

$w1=1410\u221233+\u22121+43$ | $|\lambda 1|=|4.69\u2026|>1$ $|\lambda 2|=|0.212\u2026|<1$ $(\lambda 1,\lambda 2\u2208R)$ | Saddle | Unstable | |

$w2=1410\u221233\u2212\u22121+43$ | $|\lambda |=|0.855\u2026\xb1i0.517\u2026|=1$$(\lambda \u2208C)$ | Center | Neutrally Stable |

CYLINDER SOLUTION | EIGENVALUE | TYPE | STABILITY | |
---|---|---|---|---|

$w1=1410\u221233+\u22121+43$ | $|\lambda 1|=|4.69\u2026|>1$ $|\lambda 2|=|0.212\u2026|<1$ $(\lambda 1,\lambda 2\u2208R)$ | Saddle | Unstable | |

$w2=1410\u221233\u2212\u22121+43$ | $|\lambda |=|0.855\u2026\xb1i0.517\u2026|=1$$(\lambda \u2208C)$ | Center | Neutrally Stable |

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

Now, to prove the existence of the nested closed symmetric quasi-periodic orbits around (*w*_{2}, *w*_{2}) 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 (*w*_{2}, *w*_{2}) is an elliptic symmetric fixed point. Also, (*w*_{2}, *w*_{2}) is nonresonant because of the following reason. The eigenvalues of *D***F**(*w*_{2}, *w*_{2}) 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 (*w*_{2}, *w*_{2}) are not algebraic integer, which implies that the eigenvalues of (*w*_{2}, *w*_{2}) are not the roots of unity. Thus, we can apply KAM theorem to the system at (*w*_{2}, *w*_{2}), 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

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.

*l*

_{1},

*l*

_{2}in Fig. 14 can be written in the following form:

*ρ*

_{m}. The notation of vertices is defined as Fig. 15. Because of the

*N*-fold symmetry of waterbomb tube, vertices $Am,01,Am,11,\u2026,Am,N\u221211$ form regular

*N*-sided polygon. Based on this, we derive Eq. (6). Let

*θ*

_{m,n},

*M*

_{m,n},

*P*

_{m,n},

*Q*

_{m,n}be the base angle $\u2220Am,n1Bm,n1Bm,n2=\u2220Am,n2Bm,n1Bm,n1$, a midpoint of $Am,n1Bm,n1$, a foot of perpendicular from point

*M*

_{m,n}to edge $Am,n1Am,n2$, and a foot of perpendicular from point

*M*

_{m,n}to edge $Bm,n1Bm,n2$, respectively. Here, following equations hold:

### 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 *m*th hoop, i.e., *x*_{m} and *y*_{m}, are given. First, we represent some quantities such as distance between vertices or angles by crease-pattern parameters and *x*_{m} and *y*_{m}. Next, we derive the function *f* and *g* by introducing the coordinate system and representing *x*_{m+1} and *y*_{m+1} by the coordinates of some vertices of waterbomb tube.

#### Some Necessary Quantities

First, we represent some quantities using (*α*, *β*, *N*, *M*) and (*x*_{m}, *y*_{m}) based on the preservation of edge length. The notation of vertices is defined in Fig. 16, and for simplicity, we omit subscripts of *x*_{m} and *y*_{m}.

*θ*

^{1},

*θ*

^{2}be the base angle $\u2220B1A1A2=\u2220B2B1A1$, $\u2220A2B2B1=\u2220A1A2B2$, respectively (0 <

*θ*

^{i}<

*π*). Their cosine and sine can be represented as follows:

*h*: = |

*OO*′|, where

*O*′ corresponds to the circumcenter of isosceles trapezoid

*A*

^{1}

*B*

^{1}

*B*

^{2}

*A*

^{2}because point

*O*is equidistant from

*A*

^{i},

*B*

^{i}(

*i*= 1, 2); therefore, if we represent the length of diagonals of

*A*

^{1}

*B*

^{1}

*B*

^{2}

*A*

^{2}, radius of circumcircle of

*A*

^{1}

*B*

^{1}

*B*

^{2}

*A*

^{2}as

*d*,

*r*, respectively,

*d*,

*r*,

*h*can be expressed as follows:

*i*= 1, 2), the following equations hold:

#### Derivation of *x*_{m+1} and *y*_{m+1}

*x*

_{m+1}and

*y*

_{m+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

*x*

_{m}and

*y*

_{m}. Then, we obtain

*x*

_{m+1}and

*y*

_{m+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

*m*th 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

*V*

_{x},

*V*

_{y}, and

*V*

_{z}, respectively, and

**R**

_{y}(

*η*) and

**t**denote the rotation matrix of angle

*η*with

*Y*-axis as the rotation axis and translation vector defined as

**t**≡ [0, 0,

*x*

_{m}/tan (

*π*/

*N*)]

^{T}. Here, we can represent the coordinates of each vertex as follows using Eqs. (A1), (A3):

**R**

_{y}(

*η*) we have

*m*th hoop by rotating the corresponding vertex by 2

*nπ*/

*N*around the

*X*-axis. For example,

*f*(Fig. 6).

*g*by finding

*y*

_{m+1}(Fig. 6). To do so, we first find the coordinates of vertex

*O*

_{m+1,1}, and then the coordinates of vertex $Bm+1,12$. From Fig. 6, vertex

*O*

_{m+1,1}is located at a distance 1, 1,

*l*

_{2}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

*O*

_{m+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\u2212|PCm,0|2=1\u2212xm+12$ and

*l*

_{2}, respectively. Using

**x**

_{1}: = [

*P*

_{x},

*P*

_{y},

*P*

_{z}]

^{T}, $x2:=[Cm+1,11xCm+1,11y,Cm+1,11z]T,$$R1:=1\u2212xm+12,R2:=l2$, and

**n**: =

**R**

_{x}(

*π*/

*N*)[0, 1, 0]

^{T}which is the normal vector of $P$, we can represent the intersection of the two circles as follows:

**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 *O*_{m+1,0}, *O*_{m+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\u2212(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\alpha ,x2=[Qx,Qy,Qz]T,R2=$$1\u2212(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\u22c5[0,cos\pi N,sin\pi 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