## Abstract

Compression of a stiff film on a soft substrate may lead to surface wrinkling when the compressive strain reaches a critical value. Further compression may cause a wrinkling–folding transition, and the sinusoidal wrinkling mode can then give way to a period-doubling bifurcation. The onset of the primary bifurcation has been well understood, but a quantitative understanding of the secondary bifurcation remains elusive. Our theoretical analysis of the branching of surface patterns reveals that the wrinkling–folding transition depends on the wrinkling strain and the prestrain in the substrate. A characteristic strain in the substrate is adopted to determine the correlation among the critical strain of the period-doubling mode, the wrinkling strain and the prestrain in an explicit form. A careful examination of the total potential energy of the system reveals that beyond the critical strain of period-doubling, the sinusoidal wrinkling mode has a higher potential energy in comparison with the period-doubling mode. The critical strain of the period-doubling mode strongly depends on the deformation state of the hyperelastic solid, indicating that the nonlinear deformation behaviour of the substrate plays a key role here. The results reported here on the one hand provide a quantitative understanding of the wrinkling–folding transition observed in natural and synthetic material systems and on the other hand pave the way to control the wrinkling mode transition by regulating the strain state in the substrate.

## 1. Introduction

Surface instability is frequently observed in squeezed soft materials [1–6]. The physics of this nonlinear phenomenon could be harnessed to create tuneable patterns, to fabricate functional surfaces and to understand the correlation between growth and morphogenesis in some natural systems. Therefore, it has received considerable attention of many scientists and engineers [7–14]. Biot [15] made a pioneering effort to understand surface instability of a neo-Hookean solid subject to plane strain compression. His linear perturbation analysis leads to a critical compressive (nominal) strain of 45.6% for the onset of surface instability. Recent analysis [11] based on Koiter's elastic stability theory reveals that surface wrinkling of a neo-Hookean solid is highly unstable because of the nonlinear interaction among the multiple modes associated with the critical compressive state. Concomitantly, wrinkling is sensitive to exceedingly small initial imperfections and spontaneously collapses into a local crease. The nonlinear analysis of Hohlfeld & Mahadevan [7,10] and Hong *et al*. [8] indicates that a crease/sulcus represents a subcritical instability mode, which is energetically favourable when the nominal compressive strain is beyond 35.4%. The surface instability of a stiff film resting on a soft substrate (bilayer system) is usually believed to be fundamentally different from that of a hyperelastic solid. In the former, surface wrinkles are rather stable, especially when the modular ratio of the film to substrate is large, and could be reliably observed over a wide loading range [14]. This paper is concerned with an interesting experimental phenomenon observed by Brau *et al*. [5] in Poly(dimethylsiloxane) (PDMS) film/substrate systems and illustrated by our experiments (figure 1). The experiments show that the sinusoidal wrinkling mode will transform to the period-doubling mode when the overall compressive strain reaches around 19%. This phenomenon was recently confirmed by nonlinear finite-element simulations [11]. Moreover, the nonlinear analysis [11] reveals that the critical condition for the onset of the wrinkling pattern transition largely depends on the nonlinear deformation behaviour of the substrate. Precompression favours the occurrence of folding while pre-stretch will retard the wrinkling mode transition; over pre-stretch may lead to a mountain ridge mode [12]. Hutchinson [14] investigated the role of substrate nonlinearity on the wrinkling of thin films bonded to a compliant substrate within the initial post-bifurcation range when wrinkling first emerges. His insightful analysis revealed that when the elastic modulus of the film is about five times greater than that of the substrate, the wrinkling bifurcation will be stable, whereas for a relatively soft film, the wrinkles can be unstable. However, the initial post-bifurcation expansion does not capture the unusual wrinkling modes observed at compressive strains well above the bifurcation strain such as period-doubling, folding [5,11,16] and ridging [12,17,18]. These previous studies help us to understand the physics behind the wrinkling pattern evolution in hyperelastic bilayer systems. However, a quantitative understanding of the wrinkling mode transition from sinusoidal to period-doubling in bilayer systems with the presence of substrate pre-deformation has thus far remained elusive.

This paper is organized as follows. A theoretical analysis is first performed in §2 on the branching of the wrinkling mode to characterize the correlation between the wrinkling transition strain and the parameters of the system. Computational studies are carried out in §3 to explore the deformation behaviour of the substrate in the vicinity of the second bifurcation point. A characteristic (nominal) strain in the substrate is identified at the incipient occurrence of the period-doubling, which corresponds to a critical point beyond which the period-doubling mode is energetically favourable in comparison with the sinusoidal wrinkling mode. In §4, explicit solutions are given to predict the critical overall strain for the onset of the period-doubling mode in the cases with and without substrate pre-deformation, which are validated by numerical simulations and experiments. Section 5 gives some concluding remarks.

## 2. Theoretical analysis of wrinkling pattern transition

### (a) Model

To determine the critical compressive strain for the onset of a period-doubling mode, we first perform a theoretical analysis on the wrinkling mode branching. For the system illustrated in figure 2, we explore a displacement-controlled loading procedure. The upper surface of the thin film is traction-free. Here, *x* is coordinate of the initial configuration and ℓ and *h* are the length and thickness of the film, respectively, and ℓ is assumed to be much larger than the wrinkling wavelength of the system. We assume that the film follows the generalized Hooke's law and the substrate is an incompressible neo-Hookean material. The overall compressive strain in the film in the fundamental state is *δ*. In the plane strain case, the resultant membrane stresses in the film in the fundamental state can be obtained by
*E* and *ν* being Young's modulus and Poisson's ratio, respectively.

The nonlinear von Kármán plate equations [19–21] are used to model the film in the wrinkling state (figure 2)
*w*(*x*) is the additional deflection of the film from the fundamental state. *F* is the stress function and
*R*_{0} is a constant term in the expansion of *R*. If the expansion is up to the linear term, *R*_{0} is zero, ensuring that the total vertical resistant force is zero. However, *R*_{0} may be non-zero for a nonlinear foundation, but it can be seen from subsequent analysis that *R*_{0} has no contribution to the potential energy of the system. The linear stiffness coefficient of the foundation is given by *K*=*cμ*_{s} [11,14] with *μ*_{s} being the initial shear modulus of the incompressible neo-Hookean substrate. The wavenumber of the wrinkling state is *p* and *λ*_{0} is the pre-stretch ratio of the substrate. The constant *K*_{2} is proportional to the initial shear modulus and depends on the deformation state of the substrate as described in the sequel. The stress function *F*(*x*,*y*) gives the additional stresses from the fundamental state by
*α*,*β*,*γ*=1,2 and *δ*_{αβ}=*et al*. [21] showed that at low strains, hyperelastic film/substrate systems also satisfy the Föppl–von Kármán equations. Our finite-element simulations in the sequel also demonstrate that the assumption of the film as linear elastic or hyperelastic basically has no effect on the wrinkling mode evolution when the modular ratio is relatively large.

Linearizing equation (2.2*a*) about the pre-buckling state gives the following buckling equation:
*A* refers to the amplitudes of

Substituting equation (2.6) into (2.5) gives
*p* gives

### (b) Wrinkling mode branching

Now consider the sinusoidal wrinkling mode with the additional deflection expressed by
*b*_{0} being the amplitude of wrinkling mode.

Equation (2.2) admits the stress function in the following form:
*a*_{1} and *a*_{2} represent the change of average membrane stresses along the *x*- and *y*−directions, respectively.

Invoking the following compatibility conditions [9,22]:
*et al*. [9], and the last term in equation (2.13) refers to the elastic energy of the foundation.

Minimizing the potential energy with respect to *b*_{0}, ∂*P*_{I}[*w*,*F*]/∂*b*_{0}=0, the following relationship is obtained
*K*_{2} has no influence on the potential energy for this wrinkling mode when we take the resistant force of the foundation in the form of equation (2.3) and omit higher order terms.

Next, we consider the following period-doubling mode:
*K*_{2}=0, *ΔP* is always positive, indicating that the sinusoidal wrinkling mode has a lower potential energy and the system is unlikely to evolve into the period-doubling mode. Our finite-element simulations have also demonstrated that the system will maintain the sinusoidal mode if the substrate is assumed to be a linear foundation. This demonstrates that the nonlinear behaviour of the substrate deformation plays a fundamental role in the occurrence of the period-doubling mode in film/substrate bilayer systems. For a nonlinear foundation, *ΔP*=0 determines the critical overall compressive strain *δ*_{2} in the film for the transition from sinusoidal wrinkling to a period-doubling pattern, which gives
*et al*. [5] have performed an insightful theoretical analysis on the occurrence of the period-doubling mode. They considered an inextensible film with *δ*_{c}=0 but did not explore the effect of substrate pre-stretch, i.e. *λ*_{0}=1. In this case, equation (2.23) degenerates to
*et al*. and for the case of *λ*_{0}=1, equation (2.24) gives
*et al*. [5]. Equation (2.23) indicates that *δ*_{2} depends on *λ*_{0} and *δ*_{c} as well. Provided that *K*_{2} is known, the function *f*(*δ*_{2},*λ*_{0},*δ*_{c}) can be determined and then *δ*_{2} can be solved from equation (2.23). However, it is by no means trivial to determine *K*_{2} analytically. Without considering the substrate compression, Brau *et al*. [5] have obtained an analytical solution for *K*_{2}, which is zero for an incompressible substrate; with this value, *δ*_{2} predicted by equation (2.25) will be infinite and inconsistent with nonlinear finite-element simulations [11]. Very recently, Hutchinson [14] made an effort to calculate *K*_{2}. He also obtained *K*_{2}=0 in the absence of substrate compressive deformation. His analysis [14] revealed that *K*_{2} strongly depends on the pre-deformation of the substrate and the boundary condition at the upper surface, for which an analytical solution is difficult to obtain if realistic interfacial conditions are taken into account. Bearing this in mind, a nonlinear finite-element analysis is carried out in this paper to explore the deformation behaviour of the substrate at the onset of the period-doubling mode. We will derive an explicit solution to predict *δ*_{2} without solving *K*_{2} directly.

## 3. Nonlinear finite-element analysis

### (a) Computational model

Nonlinear finite-element analysis is carried out to investigate the plane strain problem (*λ*_{3f}=*λ*_{3s}=1) of a neo-Hookean bilayer system at compressive strain levels well beyond the bifurcation strain. We concentrate on the deformation behaviour of the substrate when the wrinkling mode evolves from the sinusoidal pattern to the period-doubled pattern. Plane strain finite-element simulations were performed via the commercial software ABAQUS [25]. The incompressible neo-Hookean material model is adopted for both the film and the substrate, and the anticipated wavelength is much greater than the element size. The hybrid element (CPE8RH in ABAQUS) was used which is applicable for simulations of incompressible materials.

In the post-buckling analysis, two schemes were used to introduce the stress-free geometric imperfections to the system to initiate growth of the post-bifurcation modes [11]. With the absence of substrate pre-stretch, a linear perturbation was first accomplished using the ‘BUCKLE’ function in ABAQUS. The critical eigenmode scaled by a very small factor (e.g. 0.005 h) was introduced as a geometric imperfection into the mesh. With the presence of substrate pre-stretch, finite-element simulations were first run by specifying a sinusoidal surface deflection at the upper surface. The computed displacement field was then introduced as a stress-free geometric imperfection into the mesh. Displacement-controlled loading is employed with *u*_{1} (independent of *x*_{2}) and zero shear traction specified on the vertical sides of the model. The nominal compressive overall strain imposed on the system after the film is attached to the substrate, *δ*, is defined in terms of the film stretch, *λ*_{1f}, calculated from the relative difference between *u*_{1} on the two sides of the model. The displacement *u*_{2} and the shear traction on the bottom surface of the model are taken to be zero. The computational model has a width of the order of 10 wavelengths of the sinusoidal wrinkling mode. The depth of the substrate is taken to be more than 10 times the sinusoidal wavelength and, thus, sufficiently deep such that the substrate can be regarded as infinite. In order to capture the occurrence of the period-doubling mode, a pseudo-dynamic algorithm was adopted in the loading paths. Both loading and unloading procedures were explored to show the effect of damping on the deformation behaviour of the substrate at the bifurcation point. The results indicate that the unloading procedure permits the identification of the critical overall strain for the onset of the period-doubling mode more accurately as shown in detail below.

### (b) A characteristic strain at the onset of period-doubling mode

Figure 3 shows the nominal compressive strain of the centre of the peak in the substrate along the loading direction in both loading and unloading paths. Electronic supplementary material, Video S1, shows the evolution of the wrinkling morphology of the film–substrate system for the case without substrate pre-stretch. The strain in figure 3 is the maximum compressive strain in the substrate when the wrinkling mode of the system is sinusoidal. The location of this strain is shown in figure 4. A positive value in figure 3 refers to a compressive strain. In theory, the loading and unloading paths should coincide with each other, since the system is hyperelastic. However, the occurrence of a period-doubling mode may not be captured when the pseudo-dynamic algorithm is not adopted in our simulations. Therefore, the pseudo-dynamic algorithm is used in the loading paths, but not in the unloading paths (*β* is the damping factor). It can be seen that the two paths are not consistent with each other around the bifurcation point (figure 3) because of the effects of damping. The difference of the two paths is deceased when *β* becomes smaller. However, figure 3 shows that two unloading paths are consistent for the two different damping factors *β* and the loading path will basically coincide with the unloading path when *β* tends to zero. Thus, the unloading procedure permits us to identify the occurrence of period-doubling more accurately. The results of the unloading path show that when the nominal strain reaches a maximum value (0.381, point *d* in figure 3), the strain in the centre of the peak is released and strain localization occurs in the substrate (figure 3*c*). Point *d* in figure 3 is regarded as the bifurcation point for the onset of the period-doubling mode. The nominal compressive strain at point *d* is defined as the characteristic strain *ε*_{c}, and *ε*_{c}≈0.381. The overall compressive strain at point *d* is defined as the critical overall strain *δ*_{2}, which is approximately 0.185. To understand the physics underpinning this characteristic strain, we examine the total potential energy of the system during the wrinkling pattern evolution.

In finite-element simulations, without perturbation the system may keep the sinusoidal wrinkling mode even when the overall compressive strain is greater than *δ*_{2}, which permits us to evaluate the total potential energy *P*_{s} corresponding to the sinusoidal wrinkling mode. With the presence of perturbation for instance by introducing the damping factor, the system will evolve from the sinusoidal mode to the period-doubling wrinkling mode; in this case, the total potential energy *P*_{a} can be obtained from finite-element results. Here, *P*_{a} corresponds to the unloading path. Figure 5 shows the potential energy difference *ΔP**=*P*_{a}−*P*_{s} between the wrinkling mode that actually occurs and the sinusoidal wrinkling mode. It can be seen that when the nominal compressive strain at the point shown in figures 4 and 5 reaches the value of 0.381, the sinusoidal wrinkling mode gives way to the period-doubling mode. Beyond this critical point, the period-doubling mode is energetically favourable in comparison with the sinusoidal mode. The corresponding overall compressive strain at the bifurcation point is approximately 0.185.

Our analysis shows that the characteristic strain *ε*_{c} is basically independent of the modular ratio of the film to the substrate and insensitive to the geometric imperfection introduced in the system, as shown in figures 6 and 7. The parameter *η* in figure 7 represents the ratio of the imperfection amplitude to the film thickness. The strain *ε*_{c} is not far from the crease strain identified by Hohlfeld & Mahadevan [7,10] and Hong *et al*. [8]. But the connection between these two strains remains unclear. We also explored the case with substrate pre-stretch, i.e. *λ*_{0}≠1. Our results show that *ε*_{c} basically linearly depends on *λ*_{0} (figure 8). A linear function with *ε*_{c}=0.381+0.15(*λ*_{0}−1) is used to fit the results.

The characteristic strain identified here and the analysis on the total potential energy of the system may not only help us understand the physics behind the occurrence of the period-doubling mode and localized folding in bilayer systems but also provide a means to determine the critical overall compressive strain at the onset of the wrinkling mode transition. This issue will be further addressed below.

## 4. Explicit solutions to predict the critical overall strain for the onset of period-doubling

In this section, we show that the characteristic strain identified in §3 permits the determination of *δ*_{2} in an explicit form. Our analysis reveals that *δ*_{2} depends on *λ*_{0} and *δ*_{c}, i.e.
*λ*_{0}=1,
*f*_{1} and *δ*_{2} in equation (4.2). For the sinusoidal wrinkling mode, the maximum nominal strain *δ*_{c} and the overall compressive strain *δ* without substrate pre-stretch (figure 4).

When *ε*_{c}, *δ* will have the value of *δ*_{2}, i.e.
*g* is known, we can determine the function *f*_{1} in equation (4.2). Based on the nonlinear finite-element simulations in §2, we have determined the function *g*. This function gives the correlation between *δ* and *δ*_{c} as
*δ*=*δ*_{c}. 2/*hp*(*δ*−1/2 *δ*^{2}−*δ*_{c}+ 1/2 δ_{c}^{2})^{1/2} can well describe the evolution of the wrinkling amplitude given by finite-element simulations. For illustration, the function *g* determined from a finite-element example is shown in figure 9.

From equations (4.3) and (4.4), *δ*_{2} is derived in the explicit form:
*δ*_{2} on *δ*_{c} is rather weak when the modular ratio *μ*_{f}/*μ*_{s} is large. For instance, when *μ*_{f}/*μ*_{s}>50, equation (4.5) gives *δ*_{2}≈0.185, as shown in figure 10. Figure 11 shows that *δ*_{2} is insensitive to the geometric imperfections introduced in the post-buckling analysis. The solution given by equation (4.5) is consistent with our experimental results (figure 1) and finite-element simulations performed by Cao & Hutchinson [11].

In the case of *λ*_{0}≠1, our theoretical analysis shows that *δ*_{2}=*f*_{1}(*λ*_{0},*δ*_{c)}, and the prestrain in the substrate is *δ*_{0}=1−*λ*_{0}. In the presence of pre-stretch, the maximum (nominal) compressive strain in the substrate based on the finite-element results is given by
*λ*_{0}=1. The comparison of equation (4.6) with the finite-element simulations shown in figure 9 indicates that equation (4.6) is a very good approximation. When *ε*_{c}=0.381+0.15(*λ*_{0}−1). We have compared the predicted *δ*_{2} given by equation (4.7) with our finite-element simulations. The results are plotted in figure 12 for two different modular ratios and various pre-stretch ratios. Equation (4.7) matches the finite-element results well. It is also seen that the pre-stretch has pronounced effects on *δ*_{2}. With the increase of *λ*_{0} from 0.8 to 1.3, *δ*_{2} increases by a factor of five.

To validate further our analysis and equation (4.7), experiments have been performed. PDMS bulk was chosen as the soft substrate, and the stiff film was generated by subjecting the PDMS to a UV-ozone treatment for 10 min. Experimental details are given in appendix A. The surface wrinkling morphology of the PDMS bulk with occurrence of period-doubling is shown in figure 13 for *λ*_{0}=0.83, 0.9, 1, 1.1, respectively. It is a challenging issue to accurately determine *δ*_{2} in our experiments; therefore, it is difficult to quantitatively validate equation (4.7). However, figure 13 indeed shows that the onset of period-doubling strongly depends on the pre-stretch ratio and the experimental results basically support our theoretical and numerical results. Equations (4.5) and (4.7) are valid for the systems with relatively large modular ratio *μ*_{f}/*μ*_{s}; when the modular ratio is small, say smaller than 10, some localized buckling modes may occur instead of period-doubling [11], which is beyond the scope of this study. Very recently, Auguste *et al*. [26] studied the post-wrinkling bifurcations in films on pre-stretched substrate through a combination of simulations and experiments. They found the interesting ‘chaotic’ wrinkling when *λ*_{0}=0.7. In our experiments, the minimum pre-stretch ratio *λ*_{0}=0.83 at which stable period-doubling was observed.

The solution in equation (4.7) allows us to predict the critical condition for the onset of the period-doubling mode, and thus it provides a guideline for controlling surface wrinkling patterns by tuning the stress state in the substrate. Understanding the morphogenesis and origin of shapes has long been a central goal of developmental biology [27–29]. Growth is responsible for the morphogenesis of biological organs and tissues, which consists of a series of orchestrated steps [30]. Besides genetic and chemical effects, mechanical environments play a significant role in regulating pattern formation. Differential growth or atrophy of tissues and organs could induce residual stresses, which are important in the morphogenesis and physiological functions of soft tissues and organs [28,31–33]. Recently, Ben Amar & Jia [34] have studied the surface wrinkling and the zigzag mode of a hyperelastic bilayer soft tissue. In their study, differential growth of the soft tissue elicits residual stresses which lead to surface instability. In this paper, we have studied the development of a period-doubling mode in the plane strain compression of film/substrate bilayer systems with focus on the effects of prestresses. Equation (4.7) clearly demonstrates that prestresses in the substrate affect the evolution of wrinkling patterns and therefore it would be interesting to take this factor into account in understanding the correlation between the growth and morphogenesis of biological organs and tissues.

It should be pointed out that this study is limited to the system with a flat surface. In this sense, although figure 11 shows that *δ*_{2} is insensitive to geometric imperfections introduced in the post-buckling analysis, the geometric imperfections should be smaller than the film thickness in order to use equations (4.5) or (4.7) as appropriate. To illustrate this point further, we explore a case where geometric imperfections are introduced based on a buckled deformation corresponding to the sinusoidal mode. In the initial configuration, the buckled amplitude of film is much greater than its thickness. Electronic supplementary material, Video S2, shows the deformed configurations of the system at different overall compressive strains. In this case, equation (4.5) is no longer valid to predict the critical overall compressive strain *δ*_{2} for the onset of period-doubling. However, it is interesting to find that period-doubling still occurs when the local nominal strain in the substrate at the centre of the peak (figure 4) reaches the characteristic strain *ε*_{c}. Both this example and the solution in equation (4.5) indicate that the period-doubling mode employs the sinusoidal wrinkling mode as the base state and is almost independent of the overall strain in the fundamental state given in figure 2. Besides the wavelength, the sinusoidal wrinkling mode is primarily characterized by its amplitude. This may help us to understand why the critical condition for period-doubling is given in terms of the characteristic maximum nominal strain in the substrate.

## 5. Concluding remarks

Understanding the unusual wrinkling modes at compressive strains well above the bifurcation strain such as period-doubling and folding is a challenging issue [5,11,16,34]. In this paper, we have studied the occurrence of a period-doubling mode in the plane strain compression of film/substrate bilayer systems. Our analysis shows that the overall compressive strain *δ*_{2} in the film for the onset of the period-doubling mode depends on the substrate pre-stretch ratio and the critical overall strain for the onset of wrinkling. When the modular ratio of the film to the substrate is large, say greater than 50, *δ*_{2} only depends on the substrate pre-stretch and is insensitive to the wrinkling strain. Nonlinear finite-element analysis (compressing and unfolding the system) reveals a characteristic strain of 0.381 (with no substrate pre-stretch) for the onset of period-doubling. When maximum compressive strain reaches the value of *ε*_{c}, the period-doubling mode occurs and the potential energy of system in this case is lower than that of the sinusoidal wrinkling mode. This indicates that the period-doubling mode is indeed energetically favourable and likely to occur in practical systems. Finally, explicit solutions based on the characteristic strain have been formulated to predict the overall compressive strain for the onset of the period-doubling in the cases with or without substrate pre-stretch. The solutions proposed here match both experimental and finite-element results well.

## Funding statement

Supports from the National Natural Science Foundation of China (grant nos 11172155, 11432008), Tsinghua University and 973 Program of MOST (2010CB631005) are acknowledged.

## Acknowledgements

The authors thank the valuable and insightful comments and remarks of Prof. J. W. Hutchinson (Harvard University) on this study.

## Appendix A. Experimental details

**(a) Materials**

PDMS bulk was prepared by mixing the base and curing agent (Sylgard 184, Dow Corning) in a ratio of 40:1 by weight. The mixture was then cured at 60^{°} for 3 h, after being degassed for 1 or 2 h to remove excess bubbles. The PDMS bulk was cut into cuboids of size 20×15×7 *mm* in the experiment. The initial shear modulus of the PDMS bulk was measured as being around 0.043 MPa through performing indentation tests by using the Bose ElectroForce3100 test instrument.

**(b) Surface treatments of Poly(dimethylsiloxane)**

The PDMS (with or without pre-stretch) was subjected to a UV-ozone treatment for 10 min. A thin, stiff film was then formed onto the PDMS surface. Wrinkles formed when the film/substrate system was compressed. The tensile/compressive tests were performed using a locally made uniaxial tensile instrument. UV-ozone conditioning was provided by a UV-ozone cleaner. The modulus of the film induced by the UVO treatment is of the order 10 GPa [35].

**(c) Examination of surface morphology**

The PDMS surface was observed using a laser scanning microscope (VK-X100 K, Keyence, Japan). We observed the surface morphology while applying compression to the treated PDMS bulk. The compressive strain *δ* in figure 13 corresponded to the situation where the period-doubling mode was clearly observed.

- Received September 14, 2014.
- Accepted November 6, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.