First we consider the steady states and their bifurcations using the evolutionary solution methods of Sec. II A and II B. Numerical experiments were carried out at a succession of closely spaced parameter values, with the initial condition always taken as the steady state obtained with the previous parameter values. Since the set of all initial conditions is a function space, it is not possible to try all initial conditions, and there could be other attractors in addition to those found.

A. Overview of bifurcations

The main regions of different steady states are shown in Fig. 1 as a function of the two parameters δ and L. In some regions, two different types of steady states coexist at the same parameter values.

In regions I and III, to the left and above bifurcation arc P, the equilibrium at θ = Sin^-1δ is stable. Crossing arc P this equilibrium undergoes a supercritical Hopf bifurcation to a type (B) steady state, a limit cycle of the first kind, which is stable in the areas of II and IV near arc P.

Moving across region IV from arc P toward arc R, this type (B) solution period doubles repeatedly and becomes a chaotic steady state containing only unstable periodic motions of type (B). As arc R is crossed, this chaotic attractor loses stability in a blue sky catastrophe or boundary crisis. Observation of this bifurcation is a new result of this study, and was not indicated on the bifurcation diagram of Ref. 3.

In the areas of regions III and IV slightly above and to the right of arc Q there is a steady state of type (C), a limit cycle of the second kind. As arc Q is crossed moving into region I or II, this type (C) solution loses stability as it touches the unstable equilibrium point at θ = π - Sin^-1δ; that is, the limit cycle of the second kind becomes a homoclinic connection and vanishes in a blue sky catastrophe. The resulting transient jump settles to the stable equilibrium type (A) solution in region I, or to the limit cycle type (B) in region II.

The bifurcation arcs Q and R approach each other for small values of δ, and it becomes very difficult to locate the bifurcations because the steady states become very complicated.

B. Steady states and bifurcations for δ = 0.5

In order to illustrate the variety of steady states most fully, we chose to study in detail the behavior at fixed δ = 0.5 over a range of values of the delay L. Figure 2 summarizes the results of this study in the form of a schematic bifurcation diagram. The lower branch of this diagram follows the development of a limit cycle of the first kind from birth via Hopf bifurcation crossing arc P in Fig. 1, to ultimate demise in a global bifurcation, a chaotic blue sky catastrophe or boundary crisis crossing arc R. The upper branch of the bifurcation schematic corresponds to a limit cycle of the second kind which appears suddenly in the phase portrait as a result of a global homoclinic connection at bifurcation arc Q, and extends to L = 5 where a chaotic attractor has developed.

The schematic bifurcation diagram of Fig. 2 is labeled with lower case letters a through x which correspond to attractor portraits in Fig. 3 for various values of L. These portraits are shown projected onto the unwrapped cylindrical identification space with coordinates θ and dθ/dt. The sequence of portraits corresponds to the result of increasing the parameter L from 1.5 to 5 in small steps, and subsequently decreasing L back to 2 in order to capture all coexisting multiple attractors. There is also a sequence of portraits v, w, x with L decreasing from 3.06 beginning on the branch of steady states of the first kind, to capture additional coexisting attractors. We note that the parameter values presented in Fig. 3 are only a sample of the simulations actually performed, which included many more numerous and smaller increments and decrements of L; we shall return to this point later. The letters corresponding to individual portraits in Fig. 3 are marked on the schematic Figure 2 beneath the relevant attractor branch when they correspond to gradually increasing L, and are marked above the attractor branch when they correspond to decreasing L.

Confirmation of this entire bifurcation scenario is given in numerically computed bifurcation diagrams shown in Fig. 4.

C. Unstable limit cycles

Figure 5 shows unstable limit cycles obtained by the Harmonic balance method, together with two coexisting stable limit cycles of the second kind, for three different values of the parameter L. The unstable limit cycle approaches one stable limit cycle near the fold at L ≅ 3.407, and it approaches the other stable limit cycle as L nears the other fold at L ≅ 3.646. This confirms that both bifurcations involve dynamic saddle-node collisions of the same unstable solution branch.

In addition an unstable rotation 1 limit cycle of the second kind has been found corresponding to the rotation 2 steady state in Fig. 3(l); this is the unstable branch left by the supercritical flip bifurcation between (s) and (l). Also a rotation 1 unstable limit cycle of the second kind has been found by the Harmonic balance method at L = 4.3 where the steady state is chaotic; this is evidently the simplest of the unstable periodic solutions in the interior of the chaotic attractor. These rotation 1 unstable solutions are shown in Fig. 6.

D. Further evidence concerning bifurcations

Concerning the coexistence of two limit cycles of the first kind shown in Fig. 3 as cases (f) and (x), we surmise that between them lies an unstable limit cycle of the first kind which has a saddle-node collision with the stable (x) branch at L ≅ 3.041 and with the stable (f), (g) branch at L ≅ 3.051. We have obtained such unstable solutions, and they are illustrated in Fig. 7 as shown in Fig. 5. In Fig. 7(c), one limit cycle, which is stable in Fig. 7(b), becomes unstable due to flip bifurcation between cases (x) and (w). Also unstable limit cycles of the first kind resulting from flip bifurcations have been obtained corresponding to case (e), as shown in Fig. 8 (a). Corresponding to case (w) and (h), unstable limit cycles resulting from another flip bifurcation are shown in Fig. 8 (b) and (c). Unstable solutions corresponding to case (v) or (i), where the steady state is chaotic, have also been obtained, although they are not illustrated here.

The bifurcation at arc Q which causes loss of stability of the type (C) solution and a transient jump to type (B) has been investigated by continuing the stable solution of case (u), decreasing L to 2.229. The result is a marginally stable type (C) limit cycle shown in Fig. 9. Clearly this trajectory is essentially a rotation 1 homoclinic connection at an unstable equilibrium of saddle-focus type with θ coordinate π - Sin^-1δ. The bifurcation at Q is thus a global bifurcation, in which the period of the type (C) solution becomes infinite. A trajectory with this structure was already observed by Jelonek and Cowan in Ref. 6, but they did not mention an explanation in terms of a homoclinic connection.

Finally we consider the bifurcation at arc R. An enlarged view of the chaotic attractor of the first kind, case (i), is shown in Fig. 10(a). This picture clearly suggests that the chaotic attractor is close to touching an unstable equilibrium point of saddle-focus type with θ coordinate π - Sin^-1δ, essentially the same saddle implicated in the bifurcation at Q, but involving the opposite branch of the one-dimensional unstable manifold or outset. Figure 10(a) bears some resemblance to a chaotic blue sky catastrophe occurring in a system of three first-order ordinary differential equations introduced by Rossler in Ref. 7; see also Ref. 5, p. 238 and Fig. 12.2.

Figure 10(b) shows a transient trajectory for L = 3.068, at which value the chaotic attractor of the first kind has become unstable: the trajectory shown spends some time near the attractor which existed in case (i), and finally crosses the separator (presumed to be the inset or stable manifold of a saddle-focus point) and exits the right edge of the figure, beginning a jump to the branch of type (C) solutions.