Analytic approach to bifurcation cascades in a class of generalized Hénon-Heiles potentials

Analytic approach to bifurcation cascades in a class of generalized Hénon-Heiles potentials

Sergey N. Fedotkin,¹Alexander G. Magner,^1,

*

and Matthias Brack²

1Institute for Nuclear Research, 03680 Prospekt Nauki 47, Kiev, Ukraine

2Institut für Theoretische Physik, Universität Regensburg, D-93040 Regensburg, Germany 共Received 14 March 2008; published 27 June 2008兲

We investigate the bifurcation cascades of a linear librational orbit in a generalized class of Hénon-Heiles potentials. The stability traces of the new orbits created at its bifurcations are found numerically to intersect linearly at the saddle energy共e= 1兲, forming what we term the “Hénon-Heiles fans.” In the limit close to the saddle energy共e→1兲, where the dynamics is nearly chaotic, we derive analytical asymptotic expressions for the stability traces of both types of orbits and confirm the numerically determined properties of the generalized Hénon-Heiles fans. As a bonus of our results, we obtain analytical approximations for the bifurcation energies e_nwhich become asymptotically exact fore_n→1.

DOI:10.1103/PhysRevE.77.066219 PACS number共s兲: 05.45.⫺a

I. INTRODUCTION

The approximation of the exact density of states of a quantum system in terms of classical periodic orbits via semiclassical trace formulas is a fascinating subject which has triggered a lot of research 共see 关1,2兴 and the literature quoted therein兲. It presents a nice illustration of the corre- spondence between classical and quantum mechanics, be- sides allowing one to approximately determine quantum shell structures in terms of classical mechanics 共see 关2兴 for applications in various fields of physics兲. In Hamiltonian systems that are classically neither regular nor purely chaotic, this semiclassical theory is enriched—but also complicated—by the many facets of nonlinear dynamics.

One of them is the bifurcation of periodic orbits when they undergo changes of stability关3兴.

An essential ingredient to determine the stability of a pe- riodic orbit is its so-called stability matrix M_⬜, appearing in the amplitudes of Gutzwiller’s trace formula 关4兴, which is determined from the linearized equations of motion around the periodic orbit. The analytical calculation of M_⬜for non- integrable systems with mixed dynamics is in general not possible; the only nontrivial example is, to our knowledge, that of a two-dimensional quartic oscillator 关5兴.

In this paper we investigate the stability matrix M_⬜of the simplest orbit in a class of two-dimensional potentials which are a generalization of the famous Hénon-Heiles 共HH兲 po- tential 关6兴that has become a textbook example of a system with mixed classical dynamics. For small energies the mo- tion is dominated by a harmonic-oscillator part and is quasi- regular; at energies close to and above the saddles 共e= 1兲, over which a particle can escape, the motion is quasichaotic 共see, e.g.,关1,2,6,7兴and the literature quoted therein兲. At all energies below the saddle, there exists a straight-line librat- ing orbitA oscillating toward one of the saddles. This orbit undergoes an infinite sequence of stability oscillations and hence a cascade of bifurcations, which can be understood as the main mechanism of the transition from regular motion to chaos关8–10兴. The stability traces of the new orbitsRandL

generated at the bifurcations are found numerically 关9兴 to intersect linearly at the saddle energy 共e= 1兲, forming what have been termed the “Hénon-Heiles fans”关11兴.

In the present paper we present analytical calculations of the stability traces of both the A orbit and the new orbitsR and L bifurcating from it. The results are obtained in the limit close to the saddle共e→1兲and hence are asymptotically valid as the bifurcation energies e_n approach the saddle en- ergye= 1. They confirm analytically the numerical properties of the Hénon-Heiles fans also in the generalized HH poten- tial. As a bonus, we obtain analytical expressions for the bifurcation energies en, which are mathematically valid as- ymptotically for e_n→1, i.e., for n→⬁, and practically for nⱖ7 within five digits.

In Sec. II we present the generalized Hénon-Heiles sys- tem and discuss its shortest orbits, the bifurcation cascade of the linear A orbit and, in particular, the properties of the Hénon-Heiles fans. In Sec. III we present the basic ideas of our analytical approach and the essential results, while the technical details of our calculations are given in the Appen- dixes A and B. In Sec. IV we present an alternative pertur- bative approach for evaluating the stability traces, with the details given in Appendix C, and compare its results with those of the nonperturbative calculations.

II. BIFURCATION CASCADES IN THE HÉNON-HEILES SYSTEM

A. The generalized Hénon-Heiles Hamiltonian In this paper we investigate the following family of Hamiltonians:

H_GHH=1

2共px2+p_y²兲+1

2共x²+y²兲+␣

冉

^{−13y3+␥x2y}

冊

^{, 共1兲}

where ␥ⱖ0 is a parameter specifying specific members of the family, and ␣⬎0 is a chaoticity parameter that can be scaled away with the energy as shown below. For␥^{= 1, the} Hamiltonian 共1兲 reduces to the standard Hénon-Heiles Hamiltonian 关6兴; we therefore call 共1兲 here the generalized Hénon-Heiles 共GHH兲 Hamiltonian. The HH system with

␥^{= 1 hasC}3vsymmetry: it is invariant under rotations around

*magner@kinr.kiev.ua

the origin by 2␲/3 and 4␲/3, and under reflections at three symmetry lines with the angles⫾␲/6 and␲/2 with respect to the x axis. It exhibits three saddles at energy E_sad

= 1/6␣², the equipotential lines at E=E_sadforming an equi- lateral triangle. For␥⫽1, theC_3vsymmetry is lost and only the reflection symmetry at theyaxis remains; there are, how- ever, still three saddles over which the particle can escape.

For ␥= 0 the system becomes separable and has only one saddle on the yaxis共cf.关12,13兴兲.

After multiplying the Hamiltonian共1兲by a factor 6␣^2and introducing the scaled variablesx

⬘

^,y

⬘

^,eby

x

⬘

⁼␣^{x, y}

⬘

⁼␣^{y, e= 6}␣^2E=E/Esad, 共2兲 the scaled Hamiltonian becomes independent of␣, and for a given␥there is only one parameterethat regulates the clas- sical dynamics. For simplicity of notation, we omit in the following the primes of the scaled coordinatesx,ybut keep using the scaled energye.

For␥= 1, the three saddles are at the scaled energye= 1;

one of them is positioned atx= 0,y= 1. For␥⫽1, the saddle with energye= 1 persists at the same position, while the two other saddles lie at different energies and are positioned sym- metrically to the y axis. For a more detailed description of the topology of the potential共1兲 共and an even larger class of generalized HH potentials兲 and its shortest periodic orbits, we refer to关14兴. The shortest periodic orbits of the standard HH system 共␥^{= 1}兲 have been extensively discussed in the literature关8–10,15兴, and their use in semiclassical trace for- mulas for the quantum density of states of the HH system was investigated in关12,16–18兴.

B. The motion along theAorbit

As mentioned above, we use henceforth the symbolsx,y for the scaled coordinates 共corresponding to ␣^{= 1}兲, along with the scaled energyegiven in共2兲. The equations of mo- tion for the Hamiltonian共1兲are then

x¨共t兲+关1 + 2␥y共t兲兴x共t兲= 0, 共3兲 y¨共t兲+y共t兲−y²共t兲+␥^x2共t兲= 0. 共4兲 In the present work we focus on the linear orbit that librates along the yaxis, here called theAorbit. It goes through the origin 共x,y兲=共0 , 0兲 and toward the saddle at 共x,y兲=共0 , 1兲, which it, however, reaches only asymptotically for e→1 with a period TA→⬁. Since this orbit has xA共t兲=x˙A共t兲= 0 at all times t, its equation of motion is

y¨A共t兲+yA共t兲−y_A²共t兲= 0, 共5兲 which can be solved analytically关12兴. We give here the re- sult in the most general form, relevant for our subsequent development, where the initial point along theyaxis is given as y₀=yA共t= 0兲. The solution is then

yA共t兲=y₁+共y2−y₁兲sn²共z,␬兲, 共6兲 z=a_␬t+F共␸^,␬兲. 共7兲 Here sn共z,␬兲is a Jacobi elliptic function关19兴with argument z; its modulus␬ and the constanta_␬are given by

␬⁼

冑

^yy²₃⁻−^yy¹₁, a_␬=

冑

共y₃−y₁兲/6, 共8兲 in terms of the three real solutions of the equation e= 3y²

− 2y³ⱕ1 given by

y₁= 1/2 − cos共␲/3 −␾/3兲, y₂= 1/2 − cos共␲/3 +␾/3兲, y₃= 1/2 + cos共␾/3兲, 共9兲 with cos␾= 1 − 2e. The functionF共␸^,␬兲in共6兲is the incom- plete elliptic integral of the first kind with modulus ␬^{, the} argument ␸being determined by the initial condition:

␸^{= arcsin}

冑

^yy⁰₂⁻−^yy¹₁. 共10兲 y₁ and y₂ are the lower and upper turning points, respec- tively, of the A orbit along the y axis. The period and the action of the 共primitive兲A orbit are given by

T_A共e兲= 2

a_␬K共␬兲, S_A共e兲=12a_␬

5␣^2关E共␬兲+c_␬K共␬兲兴, 共11兲 withc_␬= −2共y3−y₂兲共2y3−y₂−y₁兲/9, in terms of the complete elliptic integrals of first and second kinds,K共␬兲andE共␬兲 共we use the notation of关19兴兲.

Note that in the limit e→1 we havey₂→1,y₃→1, and

␬→1, so thatK共␬兲andTAdiverge共whileSAremains finite兲.

The Aorbit then is no longer periodic 共and may be called a

“homoclinic orbit”关3兴兲. ExpandingTAarounde= 1, one finds the asymptotic form关9兴

TA共e兲 ⬇˜T

A共e兲= ln

冉

^{1 −432e}

冊

^{共e→1兲. 共12兲}

C. The bifurcation cascade of theAorbit in the standard HH potential

While approaching the saddle ase→1, theAorbit under- goes an infinite cascade of pitchfork bifurcations, giving rise to a sequence of new orbits R₅, L₆, R₇, L₈, . . . . This sce- nario, which has some similarities to the Feigenbaum sce- nario关20兴, was discussed extensively in关9兴, and the analyti- cal forms of the newly created R and L orbits in terms of periodic Lamé functions were discussed in关10兴.

In Fig.1 we show the traces of the stability matrix M_⬜, defined in共14兲below, of theAorbit and the orbits bifurcated from it, plotted versus energye. Whenever trM_⬜= 2, a bifur- cation occurs. We see the successive bifurcations at increas- ing energiesen; upon repeated zooming of the upper end of the energy scale neare= 1共from bottom to top兲, the pattern repeats itself in a self-similar manner. The bifurcation ener- gies e_n form a geometrically progressing series 共see 关9,10兴 for details兲cumulating at the saddle energy共e= 1兲 such that e₅共R5兲⬍e₆共L6兲⬍e₇共R7兲⬍¯⬍1, where the parentheses contain the names of the new orbits created at the pitchfork bifurcations. These are alternatively ofR共rotations兲and ofL type共librations兲.共The subscripts in the orbit names indicate the Maslov indices appearing in the semiclassical trace for- mulas; the index of theAorbit increases by one unit at each

bifurcation.兲 Due to the discrete symmetries of the system, all these pitchfork bifurcations are isochronous and hence not generic共cf. 关14兴.兲.

In Fig.2, we show again trM_⬜—in the following briefly termed the “stability traces”—of the same orbits, but this time plotted versus their respective periodsT. On this scale, trM_⬜A共TA兲 共shown by the heavy line兲 is numerically found 关8兴 for large T_A to vary as a sine function; its period ⌬T

= 3.6276 was shown in 关9兴 to be given analytically by ⌬T

= 2␲/

冑

3. The exact calculation of the function trM_⬜A共TA兲 is, however, not trivial at all. It is one of the objects of our present investigations 共see Sec. III B兲.

D. The Hénon-Heiles fans

An interesting property of the stability traces of theRand Lorbits created at the bifurcations, which has been observed numerically 关9兴 and termed the Hénon-Heiles fan structure 关11兴, is emphasized in Fig.3. Here we plot the stability traces of the primitiveAorbit and the first three primitive pairs ofR andLorbits versus the scaled energye. We note two promi- nent features共which can also be recognized in Fig.1兲.

共i兲The functions trM_⬜R,L共e兲 are approximately linear up to共and even beyond兲the barrier energye= 1.

共ii兲The curves trM_⬜R,L共e兲intersect at e= 1 in one point each for allRandLtype orbits with Maslov indices greater than 8, positioned at the values 2⫾dwithd= 6.183⫾0.001, thus forming two fans emanating from these points. The un- certainty in the parameter dcomes from the numerical diffi- culty of finding periodic orbits 共which was done using a Newton-Raphson iteration procedure兲 close to bifurcations;

our result for d was obtained for Rn and Ln⬘ orbits with 9 ⱕn,n

⬘

^ⱕ13, evaluated ate= 1. The upper limitn= 13 is due to the numerical problems only; we expect that the same value d= 6.183⫾0.001 holds also for all highern.

We found exactly the same types of HH fan for the gen- eralized HH systems given by the Hamiltonian 共1兲 for the bifurcation cascade of theAorbit along theyaxis, whereby the slopes of the fans and hence the value ofddepend on the parameter ␥. The GHH fans can be described, for large enoughn, by the empirical formula

trM_⬜^共emp兲_R,L共e兲= 2⫿cRL共␥兲共e−en兲

共1 −e_n兲 共eⱖen兲, 共13兲 where the negative and positive signs belong to theR andL type orbits, respectively. At e= 1 the curves trM_⬜R,L共e兲 in- tersect linearly at the two values trM_⬜R,L共1兲= 2⫿cRL共␥兲, so that the parameterd given above for the standard HH poten- tial isd=cRL共1兲. The numerical values forcRL共␥兲are shown by crosses in Fig. 5below.

The main goal of our paper is to find analytical support for these numerical findings. In Sec. III we will, indeed, con-

0.99998 0.99999 1.0

e

-4 -2 0 2 4

trM

0.99925 0.9995 0.99975 1.0

e

-4 -2 0 2 4

trM

0.97 0.98 0.99 1.0

e

-4 -2 0 2 4

trM

A₆

A₈

A₁₀

R₅

R₇

R₉

R₅

R₅, R₇

L₆

L₈

L₁₀

R₇

R₉

R₁₁

A₇

A₉

A₁₁

A₅

A₇

A₉

<zoom>

<zoom>

FIG. 1. Trace of stability matrix M_⬜of orbitA and the orbits created at successive pitchfork bifurcations in the standard HH sys- tem 共␥= 1兲, plotted vs the scaled energy e. From bottom to top:

successively zoomed energy scale neare= 1共after关9兴兲.

6 8 10 12 14 16 18 20 22 24

T

-4 -2 0 2 4 6 8

trM(T)

B₄ L₆ L₈ L₁₀ L₁₂ L₁₄

C₃ R₅ R₇ R₉ R₁₁ R₁₃

T=3.6276

A₅ A₇ A₉ A₁₁ A₁₃

A₆ A₈ A₁₀ A₁₂ A₁₄

FIG. 2. Trace of the stability matrixM_⬜of the orbitsA共heavy line兲,B, andC, and the orbitsR_2m−1,L_2m共mⱖ3兲created at succes- sive pitchfork bifurcations of orbitAin the standard HH potential, plotted vs their individual periodsT.⌬Tis the asymptotic period of the curve trM_⬜A共T_A兲for largeT_A共after关9兴兲.

0.99925 0.9995 0.99975 1.0

e

-4 0 4 8

trM

R₅ R₇

R₉ R₁₁ L₆

A

L₈ L₁₀ L₁₂

FIG. 3. The Hénon-Heiles fans. Trace of stability matrix of primitive Aorbit共solid line兲and the first four pairs ofR共dashed兲 and L orbits共dash-dotted兲 in the standard HH system, plotted vs scaled energye; the latter forming two fans for theRandLorbits.

firm the empirical formula共13兲analytically in the asymptotic limit e→1.

III. ASYMPTOTIC EVALUATION OF STABILITY TRACES In this section we derive analytic expressions for trM_⬜共e兲 of theA,R, andLorbits in the GHH system, which are valid in the asymptotic limite→1, i.e., close to the barrier. Before presenting them in Secs. III B and III C, we recall the defi- nitions of the stability matrix trM_⬜and of the monodromy matrixM of which it is a submatrix.

A. Monodromy and stability matrices 1. Stability matrix and the Hill equation

The analytical calculation of the stability matrixM_⬜of a periodic orbit in a nonintegrable system is in general a diffi- cult task. We recall that the stability matrix is obtained from a linearization of the equations of motion and defined by

␦␰_⬜共T兲=M_⬜␦␰_⬜共0兲, 共14兲 where␦␰⬜共t兲is the共2N− 2兲-dimensional phase-space vector of infinitesimally small variations transverse to the given periodic orbit 共N being the number of independent degrees of freedom兲, and T is the period of the orbit. For 共N= 2兲-dimensional systems, we may choose ␰⬜共t兲=共q,p兲 where q is the coordinate and p the canonical momentum transverse to the orbit in the plane of its motion.共q,p兲 then form a “natural” canonical pair of Poincaré variables, nor- malized such that共q,p兲=共0 , 0兲is the fixed point of the peri- odic orbit on the projected Poincaré surface of section共PSS兲.

For two-dimensional Hamiltonians of the form “kinetic + potential energy” H=T+V 共and particles with mass m= 1, so that p=q˙兲, the Newtonian form of the linearized equation of motion for q共t兲 becomes the Hill equation 共see the text- book关21兴for an explicit discussion兲

q¨共t兲+Vqq共t兲q共t兲= 0, 共15兲 whereVqq共t兲is the second partial derivative of the potential V with respect to q, taken along the periodic orbit, and the two-dimensional stability matrix is given by

冉

^qq˙共T兲^共T兲

冊

^=M⬜

冉

^qq˙共0兲^共0兲

冊

^{. 共16兲}

For isolated periodic orbits, solutions of 共15兲 withq共t兲⫽0 are in general not periodic. However, when the orbit under- goes a bifurcation, 共15兲 has at least one periodic solution which describes the transverse motion of the new orbit cre- ated at the bifurcation; the criterion for the bifurcation to occur is trM_⬜= + 2 共cf.关21兴兲.

For particular systems, the Hill equation共15兲may become a differential equation with known periodic solutions. For the GHH systems under investigation here, the Hill equation for theAorbit directed along theyaxis is given by共3兲, withy共t兲 replaced by yA共t兲 in 共6兲, and becomes the Lamé equation 共see, e.g., 关22兴兲 whose periodic solutions are the periodic Lamé functions共see 关10兴for the details兲. However, the ele- ments ofM_⬜in共16兲can in general not be found analytically.

One of the rare exceptions is that of the coupled two- dimensional quartic oscillator for which Yoshida关5兴derived an analytical expression for trM_⬜as a function of the chao- ticity parameter共cf.关23兴兲.

Magnus and Winkler关21兴have given an iteration scheme for the computation of trM_⬜ for periodic orbits in smooth Hamiltonians. We have tried their method for the Aorbit in the HH system, but we found 关24兴 that its convergence is too slow for computing trM_⬜A共e兲with a sufficient accuracy that would allow us to deduce the properties of the HH fans.

However, in the limit e→1, it is possible to use an asymptotic expansion of the function sn appearing inyA共t兲of 共6兲, which allows us to compute trM_⬜A共e兲 analytically, as discussed in Sec. III B.

2. Matrizant and monodromy matrix

For curved periodic orbits—such as the R and L type orbits bifurcating from the A orbit in the GHH systems—

which usually can only be found numerically, the phase- space variables␰_⬜transverse to the orbit used in the defini- tion 共14兲 of the stability matrix M_⬜ cannot be constructed analytically. Instead, one must in general use Cartesian coor- dinates and resort to the full monodromy matrix M defined below. ForN= 2, one first linearizes the equations of motion to find the matrizant X共t兲 which propagates small perturba- tions of the full phase-space vector ␰共t兲defined by

␰共t兲=兵x共t兲,y共t兲,x˙共t兲,y˙共t兲其 共17兲 from their initial values at t= 0 to a finite timet:

␦␰共t兲=X共t兲␦␰共0兲. 共18兲 For a Hamiltonian of the form H共x,y,x˙,y˙兲=₂¹共x˙²+y˙²兲 +V共x,y兲, the differential equation for X共t兲is

d

dtX共t兲=

冉

⁻U共t兲^{0 I}0²

冊

^{X共t兲 共19兲}

with the initial conditions

X共0兲=I₄, 共20兲

whereI₂andI₄are the two- and four-dimensional unit ma- trices andU共t兲is the two-dimensional Hessian matrix of the potential, taken along the periodic orbit 共PO兲:

Uij共t兲= ⳵^2V

⳵^xi⳵^xj

兵x共t兲,y共t兲其PO 共xi,xj=x,y兲. 共21兲 After共19兲is solved, the monodromy matrixM of the given periodic orbit with periodT is defined by

M=X共T兲. 共22兲

In an autonomous system,M has always two unit eigenval- ues corresponding to small initial variations along the peri- odic orbit and transverse to the energy shell. After a trans- formation to an “intrinsic” coordinate system, in which one of the coordinates is always in the directionr_储共with momen- tum p储=r˙储兲of the periodic orbit 关4兴,M can be brought into the form

M=

冢

^M0⬜

^冉

^{1 . . .0. . .1}

^冊 冣

^{, 共23兲}

where the dotes denote arbitrary nonzero real numbers and M_⬜ is the stability matrix. The diagonal elements in the lower right block of共23兲then correspond to

Mr_储r_储=Mr˙_储r˙_储= 1. 共24兲 The transformation to such an intrinsic coordinate system is quite nontrivial关25兴and not unique. For curved orbits it can in general only be found numerically and is therefore not suitable for analytical calculations. For the curved R and L orbits of our system, we therefore have to resort to the full monodromy matrix M 共22兲via the solution of共19兲. For the evaluation of their stability traces, we only need the diagonal elements of M and can then use the obvious relation trM_⬜

= trM− 2.

B. Asymptotic evaluation of the stability trace trM_A(e) fore\1

In the limit e→1, where the modulus ␬ ^{defined in} 共8兲 goes to unity, we may approximateyA共t兲by the leading term in the expansion of the function sn共z,␬兲 around ␬^{= 1} 共see 关19兴兲:

sn共z,␬兲 ⬇tanh共z兲 共␬→1兲. 共25兲 Since the function tanh共z兲 is not periodic, we have to ap- proximateyA共t兲in two portions. Takingt₂as the time where the orbit passes through its maximum aty₂, i.e.,

y_A共t₂兲=y₂ ⇐⇒ t₂=关K共␬兲−F共␸^,␬兲兴/a_␬, 共26兲 we define the asymptotic expression for theAorbit over one period by

˜y_A共t兲=⌰共t2−t兲Y1共t兲+⌰共t−t₂兲Y2共t兲, 0ⱕtⱕTA, 共27兲 where the functionsY₁共t兲andY₂共t兲are given by

Y₁共t兲=y₁+共y2−y₁兲tanh²共z兲,

Y₂共t兲=y₁+共y2−y₁兲tanh²关z− 2K共␬兲兴, 共28兲 withzgiven in共7兲. Although the function共27兲is not analytic at t=t₂, it suffices to find an asymptotic expression for trM_⬜A共e兲valid fore→1.

The details of our calculation are given in Appendix A.

The analytical asymptotic result for trM_⬜A共e兲 is given in 共A26兲in terms of associated Legendre functions. In the limit e→1, the energy dependence of trM_⬜A共e兲 goes only through the periodTA共e兲:

trM_⬜A共e兲 ⬇trM_⬜^共as兲_A共e兲= trM_⬜^共as兲_A„T_A共e兲,␥… 共e→1兲, 共29兲 where trM_⬜^共as_A^兲共TA,␥兲 is a universal function given by

trM_⬜^共as兲_A共TA,␥兲= + 2兩F˜

A共␥兲兩cos关

冑

1 + 2␥TA−⌽˜

A共␥兲兴. 共30兲 The phase function⌽˜

A共␥兲is defined through Eqs.共A28兲and 共A30兲, and the amplitude function兩F˜_A共␥兲兩is given explicitly in 共A31兲. We recall that ␥ is the potential parameter of the GHH potential 共1兲with␥= 1 for the standard HH potential.

For this case, the result 共30兲becomes

trM_⬜^共as兲_A共T_A,1兲= 2.680 439 76 cos共

冑

3T_A+ 1.567 826 96兲, 共31兲 where the numerical constants have been calculated for ␥

= 1. The period of the cosine function in共31兲 was correctly shown in 关9兴to be 2␲/

冑

3, but the phase ⌽˜

A共␥= 1兲 and the amplitude 2兩F˜

A共␥= 1兲兩were only obtained numerically. The asymptotic relation共29兲had already been observed numeri- cally in关8,9兴.

The result共31兲is shown in Fig.4 by the dotted line and compared to the exact numerical result from 关9兴, shown by the solid line. We see that the agreement becomes nearly perfect for T_Aⲏ10.5, corresponding to eⲏe₆. The asymptotic result 共29兲 and共30兲allows us to give analytical expressions for the bifurcation energiesenin the asymptotic limit en→1. The pitchfork bifurcations of theA orbit occur when trM_⬜A= + 2. We therefore define approximate bifurca- tion energiese_n^ⴱ by

trM_⬜^共as兲_A„T_A共e_n^ⴱ兲,␥…= + 2. 共32兲 Using the asymptotic form ofTA共e兲in共12兲and共30兲, we can give the solutions of共32兲in the following formulas:

e_2k−1^ⴱ ⬇1 − 432 exp„−兵⌽˜

A共␥兲− arccos关1/兩F˜

A共␥兲兩兴 + 2␲k其/

冑

1 + 2␥… 共R兲,

e_2k^ⴱ ⬇1 − 432 exp„−兵⌽˜

A共␥兲+ arccos关1/兩F˜

A共␥兲兩兴 + 2␲^k其/

冑

1 + 2␥… 共L兲, 共33兲 wherek= 3 , 4 , 5 , . . ., and the odd numbersn= 2k− 1 refer to theRtype and evenn= 2kto theLtype bifurcations, respec- tively. For e_n^ⴱ sufficiently close to 1, i.e., for large enoughn, the above values should reproduce the numerically obtained

“exact” valuesen.

7.5 10.0 12.5 15.0 17.5

T_A

-3 -2 -1 0 1 2 3 4

trMA(TA)

analytical, asymptotic numerical

FIG. 4. Stability discriminant trM_⬜Aof theAorbit in the HH potential共␥= 1兲, plotted vs periodT_A. Solid line:numerical result 共as in Fig. 2, after关9兴兲. Dotted line: analytical asymptotic result trM_⬜A^共as兲共T_A, 1兲given in共31兲.

This is demonstrated for ␥= 1 in Table I. In the second column we give the resulting values ofe_n^ⴱwith 5ⱕnⱕ16 for the standard HH system, and in the third column we repro- duce their numerical values en obtained in 关10兴 as roots of the equation trM_⬜A共en兲= + 2. As we see, the asymptotic re- sults e_n^ⴱ approach the numerical values e_n very well already starting fromn= 7, as could be expected from Fig.4. In view of the numerical difficulties in determining the en from a search of periodic orbits 共see the remarks after Fig. 3兲, the agreement is very satisfactory for allnⱖ7.

This is in itself a remarkable result, because we are not aware of any analytical results for bifurcation energies 共or bifurcation values of any chaoticity parameter兲 in noninte- grable Hamiltonian systems, except for the coupled two- dimensional quartic oscillator 共see 关10,23兴兲. In the present case, the bifurcation energiese_n can be related to the eigen- values of the Lamé equation. These can, in principle, be given by infinite continued fractions关26兴, but their determi- nation is hereby only possible numerically by iteration, which becomes even less accurate than the numerical solu- tion of trM_⬜A共e_n兲= + 2 as done in关10兴. The analytical ex- pressions共33兲therefore represent an important achievement of this paper.

C. Asymptotic evaluation of trM_R,L(e) fore\1 For the stability traces of theR andL orbits we need, as mentioned in Sec. III A 2 above, to know the diagonal ele- ments of the full monodromy matrix M, i.e., the elements X_ii共t=T兲withi=x,y,x˙,y˙. Since Eqs. 共19兲couple all 16 ele- ments of X共t兲, this is still a considerable task. It can, how- ever, be simplified considerably in the asymptotic limit e

→1. First, we can make use of the “frozenymotion approxi- mation”关in short, the frozen approximation共FA兲兴introduced in Refs. 关9,10兴. It exploits the fact that, near the bifurcation energies en at which the R and L orbits are created, their motion in theydirection is close to that of the bifurcatingA

orbit and, for increasing energye, changes only very little. It can be shown 共cf. 关23兴 and Sec. IV below兲 that this may correspond to the first order in a perturbative expansion in the parametere−en, valid to leading order in the small quan- tity 1 −e_n. Second, we can exploit some symmetry relations between the elements of M if the initial point att= 0 for the calculation ofX共t兲is chosen as the upper turning point in the direction of theA orbit, i.e., its maximum along the y axis.

These symmetry relations are derived in Sec. B1; their main consequence is that we only need to calculate the 4⫻4 sub- matrix of Xij共t兲 with spatial indices i,j=x,y, and that we have the asymptotic equality trM_⬜R,L⬇2Myy fore→1; see 共B23兲. As shown below, these symmetry relations can be used also beyond the FA, and only in order to simplify them will some properties of the FA be exploited in our further derivations.

With these approximations, the calculation of trM_⬜R,L共e兲 proceeds similarly to that of trM_⬜A共e兲discussed in the pre- vious section; its details are presented in Appendix B. The analytical result is given in共B46兲in terms of associated Leg- endre functions. After their expansion in the asymptotic limit e→1 we obtain the result

trM_⬜^共as兲_R,L共e兲= 2⫿cRL共␥兲共e−en兲

共1 −en兲 共eⱖen→1兲, 共34兲 where the − and + signs belong to theRandLorbits, respec- tively. The slope function cRL共␥兲is found analytically to be

cRL共␥兲= 4

冑

1 + 2␥

sinh关2␲

冑

1 + 2␥兴cosh

冉

^␲2

冑

48␥^{− 1}

冊

^{. 共35兲}

Equation共34兲has exactly the functional structure of the em- pirical GHH fan formula 共13兲. Mathematically, it holds as- ymptotically in the limiten→1 to leading order in the small parameter

冑

1 −e_n. We emphasize that this result confirms also the numerical finding that, for large enough n 共practi- cally, for n⬎8兲, the functions trM_⬜R,L共e兲 are linear in e fromen up to at leaste= 1.

In Fig.5we show by crosses the values ofc_RL共␥兲, evalu- ated numerically from the stability of the R and L orbits at e= 1, as a function of ␥. The solid line shows the analy- tical result 共35兲. In the lower part of the figure, we show the region of small ␥. The curve cRL共␥兲 goes through zero with a finite slope which can easily be found by Taylor expanding 共35兲 after the replacement cosh共␲

冑

48␥^{− 1}/2兲

→cos共␲

冑

1 − 48␥/2兲. The slope at␥= 0 becomes c_RL

⬘

共0兲=

冏

^dd␥^cRL共␥兲

冏

␥=0= 48␲

sinh共2␲兲= 0.563 209 42.

共36兲 This value is found analytically 关11兴 from a semiclassical perturbative approach, in which the term␥^x2yof the Hamil- tonian共1兲is treated as a perturbation. Using the perturbative trace formula given by Creagh关27兴 one can extract the sta- bilities of the R and L orbits which in this approach are created from the destruction of rational tori 共see 关11兴 for details兲. To first order in the perturbation, one obtains exactly the correct linear approximation to cRL共␥兲, with the slope TABLE I. Bifurcation energies in the standard HH potential

共␥= 1兲.e_n^ⴱare asymptotic values, calculated from the analytical ex- pressions共33兲up to 15 digits withMATHEMATICA.e_nare numerical values taken from关10兴.

n e_n^ⴱ e_n

5 0.969455224681049 0.9693090904

6 0.986829936340510 0.9867092353

7 0.999188121903970 0.9991878410

8 0.999649940584051 0.9996498

9 0.999978420334217 0.999978390

10 0.999990695444011 0.9999906955

11 0.999999426413919 0.999999424

12 0.999999752685521 0.9999997525

13 0.999999984754120 0.99999998475

14 0.999999993426398 0.99999999343

15 0.999999999594766 0.9999999996046

16 0.999999999825274 0.9999999998249

共36兲, shown in the lower part of Fig. 5 by the dotted line 关28兴.

The theoretical value ofcRL共1兲= 6.181 997 17 agrees very well with the valued= 6.183⫾0.001 that was found from the numerical stabilities of theRnandLn⬘orbits in the standard HH potential共␥= 1兲 for 9ⱕn,n

⬘

^ⱕ13, evaluated ate= 1.

Our result共34兲obeys a known “slope theorem” for pitch- fork bifurcations 关14,29,30兴. It states that the slope of trM_⬜共e兲of the new orbits createdat the bifurcation point en

equals minus twice that of the parent orbit. Specifically, in the present system, it says

d

detrM_⬜R,L共en兲= − 2d

detrM_⬜A共en兲. 共37兲 We can easily obtain the slopes of trM_⬜A共e兲 at e=en from the asymptotic result for trM_⬜^共as兲_A共e兲given in共30兲. By its Tay- lor expansion around the asymptotic bifurcation energy e_n^ⴱ given by共33兲, we find up to first order in e−e_n^ⴱ

trM_⬜^共as兲_A共e兲= 2⫾cA共␥兲共e−e_n^ⴱ兲

共1 −e_n^ⴱ兲+O

冉

共1 −^共e−e^e_n^ⴱnⴱ兲^兲3/22

冊

^{. 共38兲}

The alternating sign of the linear term is + for the R and − for theL type orbit bifurcations and thus opposite to that in 共34兲. The slope functioncA共␥兲is found to be

c_A共␥兲=

冏

^dTdA

trM_⬜^共as兲_A共T_A,␥兲

冏

T_A=T˜ A共e

nⴱ兲

= 2

冑

1 + 2␥

冑

_兩F˜

A共␥兲兩²− 1

= 2

冑

1 + 2␥

sinh共2␲

冑

1 + 2␥兲cosh

冉

^␲2

^冑

^{48␥− 1}

冊

^{, 共39兲}

where兩F˜_A共␥兲兩is given in共A31兲andT˜_A共e兲in共12兲. Note that cA共␥兲 does not depend on the bifurcation energy e_n^ⴱ since trM_⬜^共as兲_A共TA,␥兲 is a periodic function ofTA. Comparing Eqs.

共35兲and共39兲, we see thatcRL共␥兲= 2cA共␥兲so that the theorem 共37兲is, indeed, satisfied with the correct sign.

IV. PERTURBATIVE EVALUATION OF trM_R,L(e) NEARe= 1

Here we present an iterative perturbative approach for the calculation of the stability trace of the new orbits created at the bifurcation energiesenof theAorbit for the GHH Hamil- tonian共1兲, takingRorbits as the example. This approach can be useful for Hamiltonians for which we do not find symme- try properties like those given in 共B17兲 and 共B18兲, which allow for a nonperturbative calculation of the stability traces.

As the small perturbation parameter we introduce the available energy above the bifurcation point,

␧=e−en, 共40兲 which is always positive. Thexandycoordinates of the new orbits, labeledx_POandy_PO, and the relevant elements of their monodromy matrices, all as functions of time t, can be ex- panded in powers of the small perturbation parameter␧:

y_PO共t兲=yA共t兲+␧y_PO^共1兲共t兲+ ¯,

x_PO共t兲=u_PO关x_PO^共0兲共t兲+␧x_PO^共1兲共t兲+¯兴, 共41兲 X_ii共t兲=X_ii^共0兲共t兲+␧X_ii^共1兲共t兲+ ¯,

Xij共t兲=u_PO关Xij共0兲共t兲+␧Xij共1兲共t兲+¯兴 共i⫽j兲, 共42兲 wherei,j=x,y. The superscripts共m兲indicate in an obvious manner the power ␧^mat which the corresponding terms ap- pear at the mth order of the expansion. The normalization constantsu_POof x_PO共t兲are given by

uR=

冑

^e−3^en, uL=

冑

3共1 + 2^e−e␥^ny2兲. 共43兲 Note that they both are proportional to

冑

␧, so thatx_PO共t兲goes to zero in the limit e→e_n. The solution of Eqs.共B21兲with the initial conditions共B22兲for the stability trace trM_⬜R,L共e兲 using the perturbative expansions 共41兲and共42兲is presented in Appendix C for the case of the Rtype orbits. The calcu- lation for the L type orbits is completely analogous. The asymptotic result for trM_⬜R共e兲is given in共C10兲.

We now compare the nonperturbative result 共B42兲 and the perturbative approximation 共C10兲for the stability traces trM_⬜R,L共e兲with numerical results. Figure 6 shows by solid lines the asymptotic analytical results 共B42兲 for the case

0.0 0.5 1.0 1.5 2.0

0 10 20 30

c

RL

( )

numerical analytical

0.0 0.05 0.1 0.15 0.2

0.0 0.1 0.2 0.3

c

RL

( )

numerical 1st order pert.

analytical

FIG. 5. Upper panel:slope parameterc_RLof the HH fans plotted vs the potential parameter␥. Crosses:numerical values. Solid line:

the functionc_RL共␥兲given in共35兲. Lower panel: excerpt for small values of ␥. The dotted line gives the linear approximation to c_RL共␥兲, with the slope given in 共36兲, as found in a semiclassical perturbative approach关11兴.

␥= 1. They form the HH fans with their linear energy depen- dence of trM_⬜R,L共e兲 arounde= 1, intersecting at the values trM_⬜R,L共1兲− 2 =⫿cRL共1兲 with cRL共1兲⬇6.182 for the R and L type orbits, respectively. As seen from this figure, they become approximately symmetric with respect to the line trM_⬜= + 2, starting from n= 9 in good agreement with the numerical results 关10兴. Note that the linear dependence of trM_⬜R,L共e兲 共B42兲is obtained up to terms of relative order

冑

1 −en. The perturbative result for theRtype orbits共C10兲is shown by the dotted lines, already in good agreement with the analytical result共B42兲fornⱖ9.

For further comparison with numerical results, we define the slope parameter

dn=兩trM_⬜R,L共e= 1兲− 2兩, 共44兲 evaluating trM_⬜R,Lat the barrier共e= 1兲for a given orbit Rn

or Ln created at the bifurcation energyen. As shown in Sec.

III C and Appendix B 2, this parameter tends to the asymptotic limitcRL共␥兲, given in共35兲, for n→⬁.

Table II shows the slope parameter d_n 共44兲 for 7ⱕn ⱕ20, evaluated for ␥= 1 in various approximations; in the left part forRtype orbits共oddn兲and in the right part forL

type orbits共evenn兲.d_n^anin columns 3 and 7 are the nonper- turbative analytical results from共B42兲,d_n^sain column 2 rep- resents the perturbative semianalytical result共C10兲for theR orbits, andd_n^numin columns 5 and 9 are the numerical results 关10兴. Columns 4 and 8 contain d_n^num*obtained numerically from solving the equations of motion 共3兲 and 共4兲 for the periodic orbits with the FA initial conditions共B12兲and共B19兲 at the top turning point 共B1兲, and Eqs. 共19兲 att=T for the monodromy matrix elements. This approximation is in good agreement with the full numerical results for large enoughn, the better the largern, as seen from comparison of the fourth and fifth共and the last two兲columns in TableII. The bifurca- tion energies fornⱖ12 were taken analytically from TableI.

For smaller n, they were obtained by numerically solving the equation trM_⬜A共e_n兲= 2 with a precision better than 兩trM_⬜A共en兲− 2兩ⱗ10⁻⁹. As seen from this table, one has good agreement of the asymptotic behavior ofdn of the perturba- tive d_n^saand even better of the analytical resultsd_n^anas com- pared with these numerical calculations. It should be noted also that the slope parameter 共44兲 of the perturbative ap- proach 共C10兲 within the FA 关see 共C1兲兴, even without the correction共C4兲to the periodic orbit yA共t兲, is in rather good agreement with the numerical results presented in Table II, especially for asymptotically largen, with a precision better than 5%. However, the second correction in共C10兲above the FA essentially improves the slope parameter 共44兲 in this asymptotic region. As noted above, the asymptotic values of the perturbative d_n^sa and the nonperturbatived_n^an, as well as the numerical FA result for d_n^num*, all converge sufficiently rapidly to the asymptotic analytical number cRL共1兲

= 6.181 997 17 given by Eq.共35兲, in line with the analytical convergence found above from共B46兲.

Figure 7 shows good agreement between the analytical 共B42兲, semianalytical 共C10兲, and numerical solution of the GHH equations 共3兲and共4兲 for classical periodic orbits and 共19兲for the monodromy matrix with FA initial conditions for L₁₂ andR₁₃as examples. Both these curves agree very well with the asymptotic analytical slopes cRL共␥兲 within a rather wide interval of ␥ even for not too large n of the orbits mentioned above. This comparison is improved with increas- ing n, the better the largern, which gives a numerical con- firmation of the analytical convergence of trM_⬜R,L共e,␥兲 共B42兲 to the asymptotic cRL共␥兲 共35兲 at the barrier e= 1 for

0.99998 0.99999 1

-6 -4 -2 0 2 4 6

trM -2

e

L

R

L L L

R

R R

8

10 12 14

9

7

11 13

⊥

FIG. 6. Stability traces trM_⬜R,L− 2 as functions of the energye at␥= 1. Solid lines show the analytical expression共B42兲forR_nand L_n orbits with n= 7 − 14. Dotted lines are the perturbative results 共C10兲for a fewRorbits as examples.

TABLE II. The slope parametersd_n^saof the semianalytical共C10兲andd_n^anof the analytical共B42兲expres- sions vs the numerical valuesd_n^num*for solving共3兲,共4兲, and共19兲within the FA for the initial conditions at the top point共B1兲of the periodicR_n共left兲andL_n共right兲orbits, andd_n^numare the exact full numerical results关10兴 共␥= 1 in all cases兲.

n d_n^sa d_n^an d_n^num* d_n^num n d_n^an d_n^num* d_n^num

7 4.7476 5.5863 6.1688 6.1801 8 5.7796 6.2661 6.1803

9 5.9234 6.0901 6.1800 6.1819 10 6.1209 6.1951 6.1820

11 6.1391 6.1685 6.1817 6.1820 12 6.1731 6.1841 6.1897

13 6.1750 6.1801 6.1819 6.1837 14 6.1807 6.1823

15 6.1808 6.1817 6.1820 16 6.1818 6.1821

17 6.1818 6.1820 6.1820 18 6.1820 6.1820

19 6.1820 6.1820 6.1820 20 6.1820 6.1820

any ␥. For larger ␥, one needs larger n in order to obtain convergence of all the compared curves.

V. SUMMARY AND CONCLUSIONS

In this paper we have investigated the bifurcation cas- cades of the linear Aorbit in a class of generalized Hénon- Heiles potentials. We were able to derive analytical expres- sions for the stability traces trM_⬜A共e兲 of the A orbit and trM_⬜R,L共e兲of theRandLorbits bifurcating from it as func- tions of the energy, which are asymptotically valid for ener- gies close to the saddle at e= 1, i.e., in the limit where the bifurcations energiesenapproach the saddle:en→1. Our re- sults confirm analytically the empirical numerical properties of the Hénon-Heiles fans that are formed by the asymptoti- cally linear intersection of the functions trM_⬜R,L共e兲ate= 1, as given in Eq.共34兲. We found good agreement of our alter- native nonperturbative and perturbative asymptotic results for trM_⬜R,L共e兲 with the numerical results. As a bonus, we have also obtained asymptotically exact expressions for the bifurcation energies e_n of the A orbit in the GHH system, given in Eq.共33兲. Our results can be interpreted in the sense that the nonintegrable, chaotic GHH Hamiltonian becomes approximately integrable locally at the barrier, i.e., for e= 1.

Both our approaches may be useful also for more general Hamiltonians, for semiclassical calculations of the Gutzwiller trace formula for the level density 关4兴, and ex- tended to bifurcation cascades with the help of suitable nor- mal forms and corresponding uniform approximations 关3,29兴. A normal form with uniform approximation for two successive pitchfork bifurcations has been derived and suc- cessfully applied to the HH system in 关12兴. In future re- search, we hope to generalize the normal form theory to infinitely dense bifurcation sequences with the help of the results of 关13,12兴 and the theory of Fedoryuk 关31,32兴. Hereby the HH fan phenomenon for the stability traces might be useful.

ACKNOWLEDGMENTS

S.N.F. and A.G.M. acknowledge the hospitality at Regensburg University during several visits and financial support by the Deutsche Forschungsgemeinschaft 共DFG兲 through the graduate college 638 “Nonlinearity and Non- equilibrium in Condensed Matter.”

APPENDIX A: ASYMPTOTIC EVALUATION OF trM_A(e) FORe\1

To obtain the stability matrix M_⬜A for the orbit A, we have to solve the linearized equation of motion共15兲for small perturbations around the orbit in the perpendicular direction.

Since the A orbit moves along the y axis, we haveq=x,p

=x˙ and 共15兲 becomes 共3兲, which is already linear in x. We thus find M_⬜A from the nonperiodic solutions of 共3兲 with small initial values x₀=x共t= 0兲, x˙₀=x˙共t= 0兲. Let us denote these solutions byx共t;x₀,x˙₀兲. The elements ofM_⬜A共we omit the subscript “⬜A” for simplicity兲are then given by

M_qq= lim

x₀→0

x共T_A;x₀,0兲

x₀ , M_qp= lim

x˙₀→0

x共T_A;0,x˙₀兲 x˙₀ , 共A1兲

Mpq= lim

x₀→0

x˙共TA;x₀,0兲

x₀ , Mpp= lim

x˙₀→0

x˙共TA;0,x˙₀兲 x˙₀ . 共A2兲 We could not find exact analytical solutions of 共3兲using the exact function yA共t兲 共6兲 for the A orbit, for which 共3兲 be- comes the Lamé equation. Only at the bifurcation energiesen

is one of its solutions a periodic Lamé function which has known expansions 关22兴. For the nonperiodic solutions, no expansions could be found in the literature. We can, how- ever, solve 共3兲 if instead of the exact yA共t兲 we use the ap- proximation˜y_A共t兲given in共27兲, which becomes exact in the asymptotic limite→1, and for which 共3兲can be reduced to the Legendre equation as shown below. We proceed sepa- rately for the two time intervals 0ⱕtⱕt₂ andt₂ⱕtⱕTA, as specified after共25兲.

共a兲0ⱕtⱕt₂. Solve the equation

x¨₁共t兲+关1 + 2␥^Y1共t兲兴x₁共t兲= 0, 共A3兲 with the initial conditions

x₁共0兲=x₀= 0, x˙₁共0兲=x˙₀→0, 共A4兲 and obtainx₁共t2兲.

共b兲t₂ⱕtⱕTA. Solve the equation

x¨₂共t兲+关1 + 2␥^Y2共t兲兴x2共t兲= 0, 共A5兲 with the initial conditions

x₂共t2兲=x₁共t2兲, x˙₂共t2兲=x˙₁共t2兲, 共A6兲 and obtainx₂共TA兲.

To do so, we transform Eqs. 共A3兲 and 共A5兲 by defining the following variables:

z₁=z, z₂=z− 2K共␬兲. 共A7兲 Then, Eqs.共A3兲and共A5兲can be written compactly as

0 0.5 1 1.5 2 2.5

-60 -40 -20 0 20 40 60

trM -2

⊥

γ

R₁₃ L₁₂

FIG. 7. Stability traces trM_⬜R,L− 2 as functions of␥forRand L orbits, respectively, evaluated at the barrier energy e= 1. Solid lines show the analytical expression 共B42兲 for the orbitsR₁₃and L₁₂; dashed lines the asymptotic results共35兲for⫿c_RL共␥兲; dots the perturbative results to 共C10兲; and crosses the numerical results

⫿d_n^num*with the FA initial conditions as in TableII. The bifurcation energiese_n共␥兲are obtained analytically through Eqs.共33兲.