The effect of parallel static and microwave electric fields on excited hydrogen atoms

Motivated by recent experiments we analyse the classical dynamics of a hydrogen atom in parallel static and microwave electric fields. Using an appropriate representation and averaging approximations we show that resonant ionisation is controlled by a separatrix, and provide necessary conditions for a dynamical resonance to affect the ionisation probability. The position of the dynamical resonance is computed using a high-order perturbation series, and estimate its radius of convergence. We show that the position of the dynamical resonance does not coincide precisely with the ionisation maxima, and that the field switch-on time can dramatically affect the ionisation signal which, for long switch times, reflects the shape of an incipient homoclinic. Similarly, the resonance ionisation time can reflect the time-scale of the separatrix motion, which is therefore longer than conventional static field Stark ionisation. We explain why these effects should be observed in the quantum dynamics. PACs: 32.80.Rm, 33.40.+f, 34.10.+x, 05.45.Ac, 05.45.Mt


Introduction
A strong electromagnetic fields can perturb an atom in many unexpected and complicated ways that are difficult to understand. If the atom is initially in an excited state usually a large number of unperturbed bound states are coupled, making the numerical solution of Schrödinger's equation difficult. Moreover, the corresponding classical dynamics is normally partially chaotic -meaning that there are both unstable and stable orbits close to the initial unperturbed torus -and the wave function mimics this behaviour, see for instance Leopold andRichards 1994 andRichards 1996b, thus making the interpretation of numerical solutions difficult.
The investigation of the effects of strong periodic fields on an excited atom dates from the original experiments of Bayfield and Koch (1974) which showed that a relatively weak field could produce a multiphoton transition into the continuum, contrary to the received wisdom of quantal perturbation theory. The subsequent history of this interesting problem is told by Koch (1990). In 1974 conventional quantal theory required high-order perturbation theory to describe the 80 photon jumps to the continuum of this early experiment. This was, and remains, an impossible calculation, but theory was rescued by Delone et al (1978) who proposed that the classical ionisation mechanism involved diffusion of the electron through atomic states highly perturbed by the field. A year later Meerson et al (1979) proposed a different classical diffusion approximation. In the same year Leopold and Percival (1979) used a classical Monte-Carlo trajectory method to estimate classical ionisation probabilities and obtained qualitative agreement with experiment. In 1985 the experiment was repeated with better control of all important parameters and the comparison between these results and a classical Monte-Carlo simulation, van Leeuwen et al (1985), showed remarkable agreement.
Since then our understanding of the dynamics of this type of system and the relationship between classical and quantal solutions in different parameter regimes has developed. For instance we now know that the scaled frequency -the ratio between the driving frequency and the Kepler frequency of the initial unperturbed motion -is one of the most significant parameters and that there are six separate scaled frequency regions in which the dynamics has quite different characteristics, see for instance Koch and van Leeuwen (1995).
Linearly polarised fields with low scaled-frequencies were considered in Richards et al (1989) and there it was found that at particular scaled frequencies quantal effects were important due to resonances between two adiabatic states, see also Dando and Richards (1993). Low frequency elliptically polarised fields, Bellermann et al (1996Bellermann et al ( , 1997 and Koch and Bellermann (2000), showed the existence of complicated resonance structures that, on the other hand, can be explained using classical dynamics, Richards (1997). When the scaled frequency is close to unity the main classical resonance island plays a dominant and similar role in both the classical and quantal dynamics, except at certain frequencies where scars produce differences, Leopold and Richards (1994) and Richards (1996b). At higher scaled frequencies classical dynamics fails, as predicted by Casati et al (1984), and demonstrated by Galvez et al (1988) for linearly polarised fields and, for elliptically polarised fields, by Wilson (2003). These different behaviours have been the subject of several reviews Jensen et al (1991), Koch (1990Koch ( , 1995, Koch and van Leeuwen (1995).
In this paper we consider the effect of strong, parallel static and oscillatory fields on an excited hydrogen atom, the work being motivated by recent experiments of Professor Koch's group. These two fields affect the system in complicated ways. Roughly, the static field splits the hydrogen degeneracy introducing another frequency to the system with which the external field can resonate. Classically this means that the Kepler ellipse moves periodically and this motion can resonate with the driving field. Such resonances can enhance the ionisation probability so a significant part of this paper is devoted to understanding the dynamics of these resonances and the mechanism causing enhanced ionisation. Preliminary experimental results and some comparisons with classical Monte-Carlo calculations are given in Galvez et al (2000); a more detailed discussion of the experimental method with more extensive results and comparisons with classical calculations is provided in Galvez et al (2004), subsequently referred to as I.
In these experiments excited hydrogen atoms are subjected to strong parallel static and periodic fields, so the external force on the electron is (F s +F µ cos Ωt)ẑ: the accurate theoretical description of this system present a challenge. The main observed effect of this field combination is to produce an ionisation probability that, for fixed F µ and Ω, rises steadily as F s is increased, but which is punctuated by a series of approximately equally spaced sharp local maxima, see for instance Galvez et al (2000, figure 1) and figure 1 of this paper. Each local maximum is produced by a resonance between the periodic part of the perturbation,ẑF µ cos Ωt, and the mean rotational motion of the Kepler ellipse induced by the static part,ẑF s . These peaks are relatively narrow but their shapes, widths and heights vary significantly with the system parameters. Moreover, they disappear at certain values of F µ and F s , for reason that are sketched in Galvez et al (2000); similar explanations are provided by Oks and Uzer (2000) and Ostrovsky and Horsdal-Pedersen (2003).
Since publication of the original Galvez et al (2000) paper we have made detailed comparisons between the ionisation probabilities of the experiment and classical Monte-Carlo simulations, which will be published in Galvez et al (2004). These comparisons show remarkable agreement in many instances, but also some differences, some of which may be attributed to unquantifiable differences between the assumed substate distribution. In experiment and calculations we observe systematic differences between both the positions of the resonances and of their disappearances and the predictions of the simple theories. In an attempt to understand these differences more detailed numerical investigations of the classical dynamics were undertaken and these calculations show that the dynamics underlying the apparently simple, averaged results seen in figure 1 is, in fact, very complicated: even an accurate theoretical prediction of positions of the local maxima in the ionisation probability is fraught with difficulties in classical and quantum dynamics. Further, it transpires that the nature of the resonance that causes these maxima is unusual in that the classical resonant island moves as the field is switched on: in some circumstances this means that the classical ionisation probabilities reflect the development of a homoclinic tangle, see figure 18 and 19. The size of the quantum numbers needed for the quantum dynamics to mimic this behaviour is not known. This paper has two main aims. First, to provide a theoretical understanding of the classical dynamics and to isolate those features that determine the nature of these resonances, for example their positions and temporal development. Second we need to derive an approximate Hamiltonian that leads to a numerically tractable Schrödinger equation. Both aims are achieved by using a representation in which coupling between basis states is relatively small.
In section 2 notation is defined and we present some numerical results illustrating the theoretical problems that need to be solved. In section 3 the approximations are developed; this is algebraically complicated but, in essence, it is simply a standard use of perturbation theory and averaging approximations. However, for two reasons care is need with this theory. First, the convergence of the perturbation expansion needs to be understood because the fields used are strong and, it transpires, sometimes beyond the radius of convergence of the series. Second, it is shown that particular observable effects are produced by subsets of terms and it is necessary to determine their origin because only some subsets can be computed to high order -this means that some effects, such as the resonance positions, can be computed to high order, but others cannot. This analysis is necessary partly to extend the earlier approaches beyond first order and partly to show how ionisation occurs, a feature missing from all earlier theories of this system. Furthermore, from this analysis emerges an approximate Hamiltonian that should facilitate a quantal calculation.
In section 4 the mechanism connecting the dynamical resonances, described in section 3, with ionisation is described. A significant result of this analysis is that the positions of the dynamical resonance is not precisely the same as the positions of the local maxima seen in ionisation curves; although small, an estimate of the difference seems beyond current theory. In this section we determine some necessary conditions for a dynamical resonance to affect ionisation and also obtain accurate estimates for the position at which the resonances observed by Galvez et al (2000) disappear, see also Schultz (2003) and Schultz et al (2004). We also discuss resonance widths and show that for the current experiments, which involve an average over substates, the width is not due to the variations of the resonance position with substate (which is relatively small), but is due mainly to non-adiabatic dynamical effects, that are difficult to quantify.
In section 5 we analyse the time required for a resonance to develop and show this to be relatively long. Further, we observed that the nature of the resonance island means that the classical ionisation probability can be significantly affected by the field envelope. In particular in some circumstances the ionisation probability can reflect an incipient homoclinic tangle that develops as the initial state moves slowly through a separatrix.

Notation
The Hamiltonian for a hydrogen atom in parallel static and microwave electric fields, as in the experiments of Galvez et al (2000Galvez et al ( , 2004 (the latter henceforth being referred to as I), is derived in Leopold and Richards (1991) and, provided the field envelope λ(t) changes sufficiently slowly, is given by where µ is the atomic reduced mass, e the electron charge and λ(t) is the envelope function describing the passage of the atom through the cavity. This Hamiltonian has azimuthal symmetry so the z-component of angular momentum, I m , is conserved. For the particular experiments described in I, λ(t) has the 16-113-16 configuration, meaning that it rises monotonically from zero to unity in 16 field periods, remains constant for 113 periods and then decreases monotonically to zero in 16 periods. In all calculations reported here the initial rise over N a field periods is taken to be and the decrease as the field is switched off has the same shape. Classical ionisation probabilities are normally insensitive to small changes in λ(t); exceptions to this rule are discussed in section 6. For excited atoms it is convenient to use units defined by the initial unperturbed Kepler ellipse which has semi-major axis a = I 2 0 /(µe 2 ) and frequency ω K = µe 4 /I 3 0 , I 0 = n 0h , where I 0 and n 0 are the initial values of principal action and quantum number respectively: scaled units are convenient because the magnitude of most scaled parameters that produce similar physical effects change little with n 0 . The scaled frequency and field amplitude are defined by The scaled time is t 0 = ω K t and a scaled action I is I/I 0 . In the current experiment Ω = 8.105 GHz, so Ω 0 = (0.010722n 0 ) 3 . In the following we use the symbols F s and F µ for both the scaled and actual field amplitudes and t for scaled and actual time: this misuse of notation avoids a clutter of subscripts but should not cause confusion because scaled quantities are dimensionless.

Some numerical results
Before describing the theory we show the results of a few classical calculations in order to provide the reader with an idea of the features that a theory needs to describe. For the present calculations a Monte-Carlo method, as described in Abrines and Percival (1966), is used in which N initial conditions are chosen from a microcanonical ensemble: if M of these orbits ionise the estimate of the ionisation probability is P i = M/N . Without stratification the standard deviation of this estimate is (Hammersley and Handscomb, 1964) meaning that there is a 68% and 95% probability of the true result being in the ranges (P i −σ, P i +σ) and (P i −2σ, P i +2σ) respectively. In the present calculations the sample of initial conditions is stratified by dividing the range of each variable into equally probable intervals and choosing, from a microcanonical distribution, one point at random in each sub-interval, see Abrines and Percival (1966). Stratification reduces the statistical errors and sample calculations suggest that with this form of stratification the true value of σ is about half the above estimate. For the numerical integration of Hamilton's equations the problem associated with the Coulomb singularity is circumvented by using the regularisation method described in (Rath and Richards 1988): numerical integration was performed with the NAG routine D02CAF. In the first figure is shown a 'typical' classical ionisation curve in which F µ and Ω 0 are fixed and the variation of the ionisation probability, P i , with the static field F s is shown: here a microcanonical distribution of substates is used. For this illustration we choose Ω 0 = 0.0980, (n 0 = 43) and F µ = 0.10: the field envelope was 16-80-16 and 1296 orbits, for each value of F s , were used. The broad features are clear: P i = 0 for F s < 0.013 and P i = 1 for large F s , typically F s > 0.2. As F s increases between these values the steady increase in P i (F s ) is punctuated only by sharp local maxima at almost equal intervals in F s , at F s = 0.0316, 0.0620, 0.0916 and 0.1203; a close inspection of the data shows another small maximum at F s ≃ 0.148, marked by the arrow. The small amplitude undulations seen in this ionisation curve are assumed to be caused by the statistical errors mentioned above, and provide a visual estimate of the magnitude of these errors. The variation in the scalelength of these oscillations in F s is due to variations in F s -interval used in the calculations which was smallest near the local maxima of P i . We denote the positions of these maxima in P i (F s ) by F (j) s . Theory associates these maxima with a resonance between the driving field and the precession of the atomic Kepler-ellipse and we denote the position of these resonances by F ≃ Ω 0 j/3 (in scaled units), giving 0.0327, 0.0653, 0.0980 and 0.131, so the relative differences between Ω 0 j/3 and F (j) s are 3.5%, 5.3%, 7% and 9% respectively. The more accurate theory developed in section 4.2, gives F (j) s = 0.0316, 0.0623, 0.0904 and 0.112, for j = 1-4, with relative differences of 0.5%, 1.3% and 7%, respectively, for j = 2, 3 and 4, with the values for j = 1 agreeing to three figure accuracy.
The averaged ionisation probabilities disguises a richer and more complicated behaviour. In the next three figures are shown ionisation probabilities from a given substate, I 2 = I m = 0.2 -these scaled actions are defined below, equation 8 and 9. Here Ω 0 = 0.0528, (n 0 = 35), and 1600 orbits, for each F s , were used. The arrows point to  At the lowest microwave field, F µ = 0.13, a number of local maxima are seen: those labelled j = 1-4 can be associated with a dynamical resonance, and most are visible for F µ = 0.14 and 0.15. For F µ = 0.13 there are also a number of other local maxima, at F s = 0.0430, 0.0535 and 0.0705 the origin of which is not known, but in section 5 results are presented which suggest that the maximum at F s = 0.0430 is a non-integer resonance, equivalent to j = 2 2 3 ; In this example P i = 0 for F s < 0.016 and P i = 1 for F s > 0.076, with an underlying steady increase in P i (F s ) for F s > 0.045.
For F µ = 0.14 the four labelled maxima of the previous figure persist but have shifted slightly and the j = 3 maxima has split: we have confirmed that the latter effect is not due to statistical sampling errors. In other calculations with the shorter rise time of four field periods there is no split, suggesting that it is caused by the field envelope: this and other effects of the field envelope are discussed in the section 6. The j = 1 maximum near F s = 0.0166 is lower and the theory developed below shows that it disappears completely for F µ ≃ 0.147, F s ≃ 0.0161. A new feature in this graph is the large value of P i at F s = 0, with P i falling rapidly to zero as F s increases to 0.007. For F s > 0.035 the underlying trend in P i is a steady increase to unity, but all other structure is not understood.
For F µ = 0.150 the three labelled maxima persist, but have shifted and broadened. Now P i is large for F s < 0.013 and F s > 0.04 and there is a new maximum at F s = 0.0192, the existence and magnitude of which depends upon the switch-on time, see figure 3: we have confirmed that, at these low frequencies, P i is affected insignificantly by the switch-off time, provided the total interaction time is sufficiently long.
The numerical details of the resonance positions seen in these figures are given in table 1. Here the values of F (j) s are computed using a grid ∆F s = 0.0001 and defining F (j) s to be the value of F s at which P i (F s ) has a local maximum. If P i > 0.99 for a range of fields, the range is quoted and sometimes there is more than one clear maximum. For fixed F µ and increasing j there are clear systematic differences between these two values. It is not known what causes these differences; the breakdown of the perturbation expansion used to compute F (j) s maybe significant, but also the discussions in sections 4.2 and 6, show that the relation between the values of F  The regular but not equally spaced series of four minima beyond the arrows, in figure 2, have, at present, no dynamical explanation. Numerical evidence, however suggests that some of this structure is caused by the field switch, with shorter switch times removing most of the structure and longer times producing more, see figure 3: the duration of the centre part of the envelope is irrelevant provided it is long enough, see section 5. In the following two figures we show ionisation curves for the same parameters used in the right hand panel of figure 2, but with different envelopes.  (2) s = 0.0315. For the longer switch, 40-50-40, the above two maxima persist, at the same field values, but now the ionisation curve shows a great deal of other structure, some of which is sensitive to the switch time. We discuss other effects of the field switch in section 6.

Theory
In the situations of interest here the field frequency Ω is small by comparison with the Kepler frequency, so the scaled frequency Ω 0 , equation 2, is small; typically it varies between 0.04 and 0.11 for n 0 between 32 and 45 for the 8.105 GHz cavity. Hence the variation of the field, F (t), is slow by comparison with the electron orbital motion and an averaging approximation may be used to remove one degree of freedom, but first it is necessary to choose the correct representation.
The fields encountered here are sufficiently large to couple together many states of the field free atom and to ionise some of these states, so useful theoretical descriptions of the experimentally observed signal must include a mechanism for ionisation and, ideally, should use a representation in which coupling between bound states is relatively small.
In order to understand the magnitude of this problem, and to motivate the following analysis, we consider the static field ionisation of the one-dimensional atom with the Hamiltonian where F is constant and here λ increases monotonically from zero to unity over a time long compared to a Kepler period. We denote the quasi-bound states of H 1 when λ = 1 by |n F . If initially the atom is in the state |n 0 0 there is no classical ionisation provided F < F crit , where, in scaled units, F crit = 2 10 /(3π) 4 ≃ 0.1298, Richards (1987), and complete ionisation if F > F crit . In quantum mechanics tunnelling decreases these thresholds by an amount that depends upon the interaction time and the initial principal quantum number. For an excited 1d hydrogen atom initially in state |n F the probability of remaining bound at time t can be deduced from the relations derived in Richards (1987) and behaves approximately as, where for n > 20, c(n) ≃ 1/(1 + 1.65/n + 173.1/n 2 − 249.5/n 3 ) is derived by a numerical fit to the theory. This probability approaches a step function as n → ∞. For t = 2πn 3 , that is one Kepler period, the probability P b decreases from 0.9 to 0.1 as F increases over an interval γF crit where γ = 0.2, 0.1, 0.04 and 0.02 for n = 5, 10, 30 and 50, respectively. For F < F crit the state |n 0 F can be approximated by a linear combination of the field free states; the matrix elements n 0 F | n 0 provide some idea of how many unperturbed states are required to accurately describe the wave function in the presence of a strong field. This matrix element is estimated by Richards et al (1989) where it is shown to be significant for n 0 < n < m, where m ≃ n 0 (2F ) −1/4 (= 1.5n 0 for F = 0.1). Thus any realistic approximation using a basis of unperturbed states requires about 2n 0 states in a 1d system and about n 2 0 states for each of the n values of the azimuthal quantum number, m. The method used by Robicheaux et al (2002) avoids these problems, but as n 0 increases the number of grid points increases and the computational time increases commensurately.
The theory presented here minimises coupling between basis states by describing the motion in a basis that diagonalises the static-field, or Stark, Hamiltonian If |n F is a bound eigenstate of H S with energy E n (F ) then a basis that may be used to describe the bound motion of the time-dependent Hamiltonian is obtained simply by replacing the constant F by F (t); in this basis coupling between states is caused only by the rate of change F (t) not by its magnitude. For the examples of interest here it will be seen that this coupling is, in scaled units, O(Ω 0 F µ ), which is typically an order of magnitude smaller than the coupling between the unperturbed states. This method was first used in the present context by Richards (1987) and Richards et al (1989), where it was applied to the one-dimensional hydrogen atom and shown to explain important features of the three-dimensional experimental results. Because this approximation uses a bound-state basis the continuum needs to be introduced as an extra approximation, described later. In addition, the eigenfunctions r | n F are not conveniently represented by simple functions, consequently approximations to these are necessary. Finally, since canonical transformations are easier to handle than their corresponding quantal unitary transformations, it is easier to develop this approximation using classical dynamics and to quantise the resulting Hamiltonian, rather than tackle the quantum mechanics directly.

The classical Stark effect
The first goal is to find a suitable approximation to the generating function, S(I, r, F ), for the canonical transformation to the angle-action variables of H S ; we also need expressions for H S and ∂S/∂F in terms of these variables. This is a relatively routine, but complicated, calculation because it is necessary to expand to high orders in F . The main result of these calculations is the adiabatic Hamiltonian, defined by equation 18 below, which forms the basis of further approximations.
The theory starts with the Coulomb-Stark Hamiltonian, equation 3, in which the force on the electron is static and in the negative z-direction. This Hamiltonian is separable in the parabolic coordinates; following Born (1960, section 35) we use the coordinates x = ξη cos φ, y = ξη sin φ, z = 1 2 (ξ 2 − η 2 ), ξ ≥ 0, η ≥ 0, sometimes named squared parabolic coordinates, giving In the following we assume F ≥ 0 and E < 0. The Hamilton-Jacobi equation is

∂S ∂ξ
the general solution of which defines the generating function S(I, r, F ) for the canonical transformation to the required angle-action variables. This equation is separable so where α 1 and α 2 are the dimensionless separation constants that satisfy α 1 + α 2 = 2, with α 1 > 0 and α 2 > 0. Motion in the ξ-direction can be bound or unbound whereas motion in the η-direction is always bound. The bound motion is restricted to the regions 0 ≤ ξ 1 ≤ ξ ≤ ξ 2 and 0 ≤ η 1 ≤ η ≤ η 2 where ξ 1 and η 1 are zero only if I m = 0. The action variables are defined by the integrals and satisfy the relation I 1 + I 2 + |I m | = I n and 0 ≤ I 1, 2 ≤ I n − |I m |. They are related to the usual quantum numbers n 1 and n 2 by I k = (n k + 1/2)h with n 1 + n 2 + |m| + 1 = n, 0 ≤ n k ≤ n − |m| − 1.
These equations relate (I 1 , I 2 ) to (E, α 1 ) and may be inverted to give E and α 1 in terms of (I 1 , I 2 ). However, for F = 0 the integrals cannot be evaluated in closed form. One method of inverting these equations is to invert the series obtained by expanding as a power series in F . A method of performing these calculations is outlined in the appendix; the resulting algebra is complicated and performed using Maple. For reasons that will soon become apparent, we have computed these series to O(F 17 ), but here quote lower order expansions. The resulting perturbation series for the energy is where and I e = I 2 − I 1 ; this last action variable is related to the electric quantum number, I e = n eh , though in some text, for instance Bethe and Salpeter (1957), the electric quantum number is defined to be n 1 − n 2 , because coordinates are chosen so the force on the electron is in positive z-direction; the value of I e is the projection of the Runge-Lenze vector along Oz. Up to O(F 5 ) the above expansion agrees with the series given in Damburg and Kolosov (1983, page 45), see also Silverstone (1978), ash → 0. Note that the odd components E 2k+1 (I) have a term linear in I e ; it will be shown that these components determine the resonance position. The separation constant is, to O(F 3 ), The angle-variables corresponding to the action variables defined in equations 8 and 9 are defined by θ k = ∂S/∂I k , k = 1, 2. It is shown in the appendix that the two angle variables can be expressed in terms of relations where the two auxiliary angles, (ψ, χ) are defined by the relations and where (P k (x), Q k (x)) are odd, 2π-periodic functions with zero mean value; expressions for these, accurate up to O(F ), are given in the appendix, equation 50. However, for reasons discussed below, we require these functions only in the limit F = 0 and then we have where σ k = I k (I k + I m )/I n , k = 1, 2.

Dynamic Stark effect
When the field amplitude F varies with time the function, S(I, r, F (t)), generates a time-dependent canonical transformation and the Hamiltonian becomes The first term of this is the Stark Hamiltonian, equation 10; the second term is more difficult to find, but it is important because only this term mixes states. It is shown in the appendix that the function ∂S/∂F can be expressed as a Fourier series of the following form where the angles (ψ, χ) are defined in equation 13 and where the coefficients (A k , B k ) are functions of the action variables and F . It is important to note that there is no term independent of both ψ and χ. In our applications F (t) = λ(t)(F s + F µ cos Ωt) andλ/λ ≪ Ω, soḞ ≃ −λΩF µ sin Ωt and since, in scaled units Ω, F s and F µ are numerically similar and small, a second order approximation is obtained by evaluating the derivative ∂S/∂F at F = 0, which considerably simplifies the analysis: the following result is derived in the appendix, It is more convenient to use new angle-action variables, so, when F = 0 equation 14 gives 2φ e = θ 2 − θ 1 = χ − ψ. This gives an approximation we name the adiabatic Hamiltonian, where E(I, F ) is the Stark energy given in equation 10. The angles ψ and χ are not conjugate to the action variables and to develop further approximations it is necessary to express all quantities involving ψ and χ in terms of θ 1 and θ 2 , using equations 14. This is most easily achieved by expressing each function as a multiple Fourier series, The adiabatic Hamiltonian 18 is useful because the coupling term is O(Ω 0 F µ ) rather than O(F ) as in the original Hamiltonian, so the resulting Schrödinger equation may be solved using a far smaller basis.
Hamilton's equations in the original representation are singular at r = 0 and this is dealt with using a regularisation method, see Szebehely (1967) for a general introduction and Rath and Richards (1988) for an application to the perturbed hydrogen atom. The equivalent singularity in the adiabatic Hamiltonian occurs when I m = 0 and σ 1 + σ 2 = 1 and this also needs to be removed. For instance the equation for ψ isψ = [ψ, K] and the right hand side is proportional to 1/J where J = I n (1 − σ 1 cos ψ − σ 2 cos χ), which can be zero. A method of avoiding numerical problems when J is small is to define a new time, τ , by the equation dt/dτ = J, to give the equations of motion where κ = I 4 n /(2µ 2 e 6 ) and I n = I 1 + I 2 + I m . In figure 5 below we compare ionisation probabilities computed using these equations and the original Hamiltonian, but first it is necessary to re-introduce ionisation into this approximation.

Ionisation
The adiabatic Hamiltonian, equation 18, does not allow for ionisation because angleaction variables exist only for bound orbits. Ionisation therefore has to be included as an extra approximation which is described here.
For static fields each classical state, or torus, labelled by the actions (I n , I e , I m ), has a critical field F crit such that it exists only if 0 ≤ F < F crit : the approximation to F crit given by Banks and Leopold (1978) is used here. Note that if I m = 0 bound orbits exist for all F , Howard (1995), but those orbits that exist for large F are so special that they do not affect the current problem. Adiabatic invariance suggests that this behaviour persists for sufficiently slowly varying fields, that is small scaled frequencies.
The extreme values of F crit occur when the atom is aligned along the field, I m = 0: in scaled units min(F crit ) ≃ 0.13 for (I e , I m ) = (I n , 0), and max(F crit ) ≃ 0.38 for (I e , I m ) = (−I n , 0). The variation of F crit with I e = I 2 − I 1 , for various values of I m , is shown in the following figure. I e Figure 4 The scaled classical critical field F crit as a function of the scaled action Ie = I 2 − I 1 , for various Im.
In the adiabatic limit I n and I e are almost constant, so we may define the timedependent critical field F crit (t) = F crit (I e (t), I m ) and assume that ionisation occurs at the time when the actual field F (t) exceeds this. If F (t) changes sign, as happens if F µ > F s , the quantisation axis changes direction and the relevant critical field is given by F crit (−I e , −I m ) (F < 0) so the combined ionisation criterion is: This approximation is accurate when the field changes very little during one Kepler period, that is Ω 0 ≪ 1, and then this criteria may be used to include ionisation in the adiabatic equations of motion, equations 20. A useful guide to the behaviour of the system is obtained by putting I n (t) and I e (t) equal to their initial values. This gives two boundaries beyond which P i = 1: (0) is the initial value of this action variable.
For the ionisation curves shown in figure 5, I e (0) = −0.6 and I m = 0 giving F + crit = 0.219 and F − crit = 0.142. Thus, since F µ = 0.15 the condition II gives P i (F s ) = 1 if F s < 0.008 or F s > 0.069: these boundaries are shown by the arrows in the figures and are consistent with the calculations. Similar boundary-arrows are included in figures 2 and 3.
Non-adiabatic dynamics affect this simple picture in two ways. First, they blur and slightly shift the boundaries. Second, and more important, isolated resonances produce large changes in I e (t) and can enhance ionisation at particular combinations of (F s , F µ ), other than those defined by I and II above. These dynamical effects produce the peaks labelled j = 1-4 seen in figure 5.
In the following figure we compare values of P i computed using the original Hamiltonian 1, the solid lines, and the adiabatic equations 20, the dashed lines.
Here These figures show broad agreement between the two calculations, but there are three marked differences; consider the left hand panel.
• The maxima j = 2, 3 and 4 are at different values of F s . This is because only terms upto O(F 2 ) were included in the expansion of E(I), equation 10, and is not an inherent inaccuracy of the adiabatic approximation.
• The main difference is the shift in the shoulder near F s = 0.01.
Apart from these differences the agreement between the two approximations is good. The same remarks apply to the right hand panel but now we see that the new maxima at F s ≃ 0.039 and some of the structure at F s ≃ 0.07 are reproduced in the adiabatic approximation. The adiabatic approximation also has a local minimum at F s = 0.0025, not present in the 'exact' probabilities: however, similar behaviour is seen in the exact results for other parameters values, see for instance the right hand panel of figure 3.
These results, and other comparisons that cannot be shown here, suggest that the adiabatic Hamiltonian provides a good approximation to the true dynamics. This is important because, for the principal quantum numbers used in current experiments, the numerical solution of the Schrödinger equation derived from this Hamiltonian is a feasible computational task, unlike that derived from the exact Hamiltonian using either an unperturbed or a static-Stark basis. In the quantal approximation ionisation is included by adding a complex part to the energies, which may be computed semiclassically, as in Leopold and Richards (1991) and Sauer et al (1992).

Averaged equations of motion
The adiabatic Hamiltonian 18 can be simplified further by noting that the two natural frequencies of the motion are quite different, so |ω e | ≪ ω n ≃ ω K ; this means that the orbital elements of the Kepler-ellipse change relatively little during one Kepler period. Hence the first averaged approximation is obtained by averaging over φ n , which is most easily achieved by ignoring all terms containing φ n = (θ 1 + θ 2 )/2 in the Fourier series 19. This gives sin ψ = − 1 2 σ 2 sin 2φ e , sin χ = 1 2 σ 1 sin 2φ e and sin kψ = sin kχ = 0, k ≥ 2.
Substituting these mean values into the adiabatic Hamiltonian 18 gives the mean motion Hamiltonian where In quantum mechanics this approximation corresponds to ignoring all transitions between states of different n-manifolds of the adiabatic basis.
Numerical integration of the equations of motion derived from 22, is not straightforward because of the square root singularity of B(I e ) at I e = ±(I n − I m ). It is best accomplished by introducing the three-dimensional vector Z = (B(I e ) cos 2φ e , B(I e ) sin 2φ e , I e ), the components of which satisfy the commutation relations [Z i , Z j ] = 2ǫ ijk Z k and |Z| 2 =constant so the vector Z moves on the surface of a sphere. The equations of motion,Ż k = [Z k , H], become where κ = I 3 n /(2µ 2 e 6 ). These equations are trivially solved numerically.

The Resonance Hamiltonian
The mean motion Hamiltonian 22 needs to be rearranged in order to extract a clear picture of the dynamics. Observe that the odd terms in the series for the Stark energy, equation 10, contain components that are linear in I e and that these terms produce a slow secular change in φ e which physically corresponds to a rotation of the angular momentum vector about the Runge-lenze vector. We shall see that the mean part of this motion determines the position of the resonances and the oscillatory part causes these resonances to disappear at certain field ratios. Denoting these linear terms by E L (I e , F ) and on setting µ = e = 1, we have, to The remaining part of the Hamiltonian,Ḟ ∂S/∂F , has a different form, equation 16, and does not give rise to factors like I e F k ; this is important for the following analysis.
Using E L only in Hamilton's equations we find that φ e (t) = φ e (0)− 3I n g(t)/2 where (24) The series in this integral has a finite radius of convergence, F rc (I m ), so it is important that max(F ) < F rc . We have computed this series to O(F 17 ), and used these nine terms to estimate F rc . By extrapolating the ratios of coefficients using Richardson's extrapolation we estimate that F rc ∼ 0.17, 0.19 and 0.21 for I m = 0, 0.8 and 1, respectively. Using Padé approximants we obtain F rc ∼ 0.18, 0.20 and 0.22, respectively. This provides a rough guide to the range of fields for which the following theory is valid.
Since the field amplitude, F (t), is periodic in t, the function g(t) can be written in the form g(t) = gt +g(t), whereg(t) is periodic with zero mean and g is the mean of E L over a field period. With F (t) = F s + F µ cos Ωt this becomes The periodic functiong(t) can be expressed as the Fourier series The dominant harmonic isg 1 : bothg 2 andg 3 are 0(F 2 ) and all higher harmonics are O(F 4 ) and may be neglected. Resonances in the dynamics occur when the angular frequency, 3I n g/2, resonates with the field frequency: the magnitude of the effect of any resonance depends upon the periodic component of g(t), and principally upong 1 .
In order to see this we change to a moving reference frame, in which φ e (t) is approximately stationary, by defining a new angle ψ e = φ e + 3I n g(t)/2 using the generating function F 2 (p, φ e ) = (φ e + 3I n g(t)/2) p, where (ψ e , p) are the new conjugate variables. Since g(t) is, by definition, independent of I e we have p = I e . In these variables the Hamiltonian 22 becomes ) . (27) No further approximation has been made in deriving this Hamiltonian from K m defined in equation 22. By definition the curly brackets contains terms quadratic and higher in I e which are independent of ψ e ; the leading term is and sinceḞ = F µ Ω sin Ωt we see that the terms of K m are O(F 2 ) and O(F Ω), and since Ω ∼ F ∼ 0.1 (in scaled units) the mean motion of (ψ e , I e ) is slow by comparison with the field oscillations: hence the relatively high frequency oscillations of E − E L do not qualitatively affect the motion -a fact that has been confirmed numericallyand hence we may replace E − E L by its mean over a field period. Retaining only the dominant quadratic term gives sin Ωt sin (2ψ e − 3I n g(t)) .
The second term is, for most values of Ω, F µ and F s , an oscillatory function of time with small mean value: in these circumstances it has little effect and may be ignored so K m ∼ I 2 e , giving I e (t) ∼ constant with ψ e approximately proportional to t. However, for any given (Ω, F µ ) there are particular resonant values of F s for which the longtime average of the second term is proportional to sin 2ψ e and then the nature of the resonant motion is qualitatively different. Near these values of F s , I e (t) can vary over a large portion of its accessible range and in some circumstances this leads to enhanced ionisation.
The functiong(t) is periodic and odd so we may write where the Fourier coefficients, J k , depend upong s , s ≥ 1. The coefficient J k is dominated byg 1 = 1 + O(F 2 ), butg 2 andg 3 are also O(F 2 ), and to this order where z k = 3g k I n F µ /Ω. Using equation 28 the mean-motion Hamiltonian becomes, where the functions A(I e ) and B(I e ) are defined after equation 22 and J j = (J j−1 − J j+1 )/2, so that, to the lowest order, J j = J ′ j (3F µ /Ω 0 ). If |ν j | is small the jth term of the sum changes more slowly than all other terms, which may therefore be averaged over to give the resonance Hamiltonian, We define the jth dynamical resonance to be at the static field, F s , where ν j = 0, that is where coupling between I e -states is largest: the equation for F (j) s is 3g(F s , F µ , I m )I n = jΩ or to lowest order, in scaled units, (31) We show below that near this value of F s ionisation may be enhanced, but in section 6, it is shown that there is no clearly defined, precise relation between F (j) s and F (j) s , the position the maximum in the ionisation probability seen in figures 1 and 2; the two fields are close but the difference can be larger than the resonance width, see section 4.2, in particuar table 2.
The position of the jth resonance, F s , is, to a first approximation, independent of the substate quantum numbers; if one substate is ionised by this resonance others will be similarly affected so the effect of the resonance is not significantly changed by an average over substates.
With the translation θ = ψ e − ν j t/2 + jπ/2, the resonance Hamiltonian may be converted to the conservative system, where The Hamiltonian K R shows that the jth resonance disappears when J j = 0: to the lowest order this gives, in scaled units where j ′ jk is the kth positive zero of J ′ j (x). This critical value of F µ was first derived in a linear quantal approximation, Galvez et al (2000), and later by Oks and Uzer (2000) using a Floquet approximation and by Ostrovsky and Horsdal-Pedersen (2003) using a linear approximation. Recent experiments, Schultz (2003) and Shultz et al (2004), and comparisons with classical calculations suggest that this estimate can be inaccurate by up to 10%. Later we show that the present theory can be used to improve upon these estimates.
The derivation of the resonance position 31 and the resonance Hamiltonian 32 involves a series of approximations. Before proceeding it is useful to list these in order to estimate the effect of ignored terms.
1. We have used the Stark angle-action variables, defined by H S , equation 3. For the action variables we use a series representation in F . For the angle variables we use the F = 0 limit because these variables appear in the Hamiltonian K, equation 15, only in the term which is O(Ω 0 F µ ).
2. The term ∂S/∂F , equation 16, has zero mean value when averaged over (ψ, χ) and because it is multiplied by the factor Ω 0 F µ we may approximate it by its value at F = 0.
3. The mean motion Hamiltonian, K m equation 22, is derived by averaging over φ n = (θ 1 + θ 2 )/2 which replaces sin kψ and sin kχ by Fourier series in sin 2φ e . The approximations described in points 1 and 2 truncate this series at the first term and approximate its coefficient to second order.
The inclusion of higher-order terms in the mean-motion Hamiltonian 22 introduces corrections O(ΩF µ F ) to the factor ΩF µ A(I e )B(I e ) and adds further terms corresponding to the harmonics sin 2pφ e , p = 2, 3, · · ·. Crucially this means that the estimate, ν j = 0, of the j resonance position is not affected by the approximations made: that is, the position of the dynamical resonance is determined solely by the parts of the Stark Hamiltonian, E(I), equation 10, which are linear in I e .
A better estimate of F (j) s is therefore obtained using the series 25, which has been evaluated, using computer assisted algebra, to O(F 17 ). In section 4.2 we use this to obtain better estimates of F

Qualitative discussion
Here we show how the dynamical resonance described in the previous section can enhance ionisation. The connection is qualitative, but explains many featues of the ionisation probablity.
The Hamiltonian K R , equation 32, is similar to that of a vertical pendulum, but there are two significant differences. First, I e is confined to the region |I e | ≤ I n − I m with natural boundaries at I e = ±(I n − I m ), where B(I e ) = 0, see equation 22. Second, the coefficient of cos 2θ depends upon I e . The fixed points of K R are at the roots of ∂K R /∂I e = ∂K R /∂θ = 0, and for this analysis it is convenient to replace J j by | J j |; when J j < 0 this represents a physically unimportant translation in θ. The equation ∂K R /∂θ = 0 gives, for 0 ≤ θ < π, θ = 0 and π/2. At θ = π/2 the equation ∂K R /∂I e = 0 has a single root near I e = α j and this fixed point is a centre. At θ = 0 there are generally three roots: a saddle near I e = α j , but there are two others near I e = ±(I n − I m ), because of the square-root singularity in A(I e ). If the (θ, I e ) phase plane is projected onto a sphere with latitude ψ, so I e = (I n − I m ) cos ψ, it is seen that there are phase curves with centres close to, but not at, the poles and which enclose the poles: in the Cartesian coordinate system (θ, I e ) this produces the two extra fixed points. The physically significant fixed points are near I e = α j and these exist only if |α j | < I n − I m , approximately.
By plotting the contours of K R for fixed F µ and Ω and with F s increasing so α j increases from below −(I n − I m ) to above I n − I m we see how the resonance develops and why, in certain circumstances, ionisation is enhanced if F s ≃ F (j) s . The following five figures show the contours of K R , near the j = 1 resonance for I m = 0.2, Ω 0 = 0.0528 and F µ = 0.13, corresponding to figure 2. For these graphs we use K R withg 1 and g given by equations 26 and 25 (to O(F 5 )) respectively; these give ν 1 = 0 when F (1) s = 0.0168. In each figure F s and F µ are fixed so, according to the adiabatic ionisation criterion, there is a critical value of I e , above which orbits ionise, given by the solution of F crit (I e , I m ) = F s + F µ : for the parameters used here the critical value of I e changes from I e = 0.50 (F s = 0.0152) to I e = 0.42 (F s = 0.0183). The maximum of F crit (I e , I m ) is at there is no ionisation, even at a resonance. The upper solid horizontal line in each figure is at the value of I e at which F crit (I e , 0.2) = F s + F µ , so orbits straying above this line will ionise. The lower solid line is I e = −0.4 is taken, for illustrative purposes, to be the initial state. Note that in this case F µ − F s < 0.13, so the other adiabatic boundary defined in equation 21 does not lead to ionisation. As F s increases from 0.0152 to 0.0183, through F (1) s = 0.0168, the centre of the resonance island moves upwards. The physical effect of this is understood by considering a field suddenly switched on, with the initial value of φ e uniform in (0, 2π) and the initial value of I e , to be −0.04.
• Figure 6: F s = 0.0152 (α 1 = −0.75). The adiabatic condition shows that orbits for which I e (t) > 0.50 ionise; at this field there is no resonance island and no ionisation from the initial state.
• Figure 7: F s = 0.0163 (α 1 = −0.25). The resonance island exits; it intersects the initial state, but does not overlap the ionising region, so there is no ionisation. In practice the demarcation between ionising and non-ionising regions is less sharp because the averaging approximations used to derive this simple picture replaces unstable manifolds by separatrixes.
• Figure 8: F s = F (1) s = 0.0168 (α 1 = 0). The resonance island is in the centre of the phase space. The ionisation criterion is practically the same as in figure 6 so orbits with I e (t) > 0.50 ionise. Now, however, the resonance island can transport initial states to the ionising region, I e > 0.46. Note that not all orbits trapped in the resonance island will ionise, but only those near the separatrix. We shall see in section 5 how this affects the ionisation times.
• Figure 9: F s = 0.0173 (α 1 = 0.25). The centre of the island is now at I e = 0.2 and its separatrix just dips below the initial state, so few orbits ionise. In these circumstances it is shown in section 6 that the ionisation probability can be affected significantly by the way the field is switched on.
• Figure 10: F s = 0.0183 (α 1 = 0.75). As for figure 6 the island no longer exists and there is no ionisation from the initial state.
This qualitative description of the ionisation process suggests that for a microcanonical distribution of substates the background ionisation increases as F s increases across a resonance, because I c e decreases, as shown in figure 14. The centre of the jth resonance island is at approximately see equation 32, so it exists only for F s in the interval In this field range a high proportion of initial values of I e (0) may lead to ionisation: outside this interval the resonance does not exist. This qualitative description of the ionisation mechanism shows that for a system initially in a given I m -substate there are several conditions necessary for a resonance to enhance the ionisation probability.
1) The field amplitudes must be sufficiently large that there is ionisation for some value of I e , for the given I m .
2) The field amplitudes must not be so large that P i = 1.
3) If F crit (I c e , I m ) = F s + F µ , then the island width must exceed I c e − I e (0), otherwise the initial state cannot be transported to an ionising state. There is, of course, a similar relation for the boundary defined by F crit (−I c e , −I m ) = F µ − F s when F µ > F s .
The first two of these conditions define a region in the (F s , F µ )-plane in which the resonance may enhance ionisation. The boundaries of this region depend upon I m and are the complement of the region defined by the two conditions P i = 0 and P i = 1 for all I e . Using the adiabatic assumption, P i = 1 if F s + F µ > max(F crit ) and P i = 0 if F s + F µ < min(F crit ) and F µ − F s < min(F crit ), if F s < F µ , see equation 21. In the case I m = 0.1 this region is shown by the shaded area in the following figure: outside this region the resonance can have no effect. If |I m | > 0.1 the equivalent region lies inside the area shown. Inside the shaded region a resonance affects the ionisation probability only if condition 3 above is satisfied. Below, but near, the upper boundary P i ∼ 1, so resonance peaks are barely noticeable, see for instance the j = 5 resonance in figure 1 at F s ≃ 0.148. The break down of adiabatic invariance broadens these boundaries but in a manner difficult to estimate, although the effect increases as Ω 0 increases; tunnelling also affects these boundaries.
There are three parameters of the resonant island that affect P i (F s ). These are most easily estimated by setting F s = F (j) s , so α j = 0, and I m = 0 as well as using the lowest-order estimates of all variables: these values are used in the remainder of this section.
The first parameter is the island area, A j , which determines how the classical resonance affects the quantum dynamics. An estimate for this is, For other values of I m and α j the form of the resonance Hamiltonian shows that A j is proportional to |J ′ j |. For the parameters of figure 1, Ω 0 = 0.098, F µ = 0.1 this gives A j /2π = (0.80, 0.05, 0.35, 0.25)I n for j = 1, 2, 3 and 4. In this case n 0 = 43 so the approximate number of states associated with these islands are 34, 2, 15 and 11, respectively. Thus if n 0 is decreased, all other scaled variables remaining the same, we should expect, in quantum dynamics, the j = 2 resonance to become less prominent than the other resonances. The second parameter is the island width, ∆I e , that is the maximum distance between the two branches of the separatrix. It is difficult to derive a simple estimate of ∆I e ; here we simply note that it is proprtional to |J ′ j |, defined above. A necessary condition for enhanced ionisation is that ∆I e is larger that I c e − I e (0); otherwise transport to ionising regions does not occur. Notice that this condition is independent of the principal quantum number, n, unlike that discussed above.
The third important classical parameter is the period of the mean motion inside the island; as we shall see, this determines how rapidly a resonance develops, section 5, and how it is affected by the field envelope, section 6. The frequency, ω j , of the motion inside the resonance island is approximated by expanding the resonance Hamiltonian, equation 32 about its centre and near the island centre we obtain, in scaled units This estimate gives the largest value of the frequency in the island; for motion nearer the separatrix ω j is smaller. In this section we have shown how a dynamical resonance can enhance the ionisation probability and have derived some approximate necessary conditions. The analysis uses the resonance Hamiltonian, K R , equation 32, derived using two stages of averaging. Moreover, an important part of this Hamiltonian is the factor J j , which, for small Ω 0 , oscillates between its maximum and minimum values for relatively small changes in F µ . Hence, whilst K R provides a good qualitative description, for any fixed (Ω 0 , F µ ) the details may be wrong; for instance the field at which a resonance disappears is given inaccurately by this approximation if Ω 0 is small.

Resonance positions
In this section we examine ionisation from a particular substate and compare theoretical predictions with exact numerical calculations. We choose the low frequency Ω 0 = 0.011414, (n 0 = 21) (to minimise non-adiabatic effects), fix I m = 0.2, use the initial condition I e = −0.4 (so there is no average over substates) and put F µ = 0.13.
Since I e is an approximate constant of the motion and F crit (−0.4, 0.2) = 0.1984, if the dynamics were adiabatic we should expect complete ionisation when F s exceeds F crit − F µ = 0.0684 and no ionisation for smaller static fields. At a resonance I e (t) varies over part of its accessible range and since min(F crit ) = 0.1357 we might see the effect of the jth resonance if F (j) s ≥ 0.0057, provided the size of the the resonance island is sufficiently large. This simple analysis suggests that the resonances 2 ≤ j ≤ 18 could be seen via ionisation: in practice, for reasons to be discussed later the j = 2, · · · , 6, 11 and 14 resonances are not observed.
In figure 12 we show the classical ionisation probabilities for the envelope 16-50-16 in which the j =7-10, 12, 13 and 15-19 resonances are clearly visible, but the j = 6, 11 and 14, marked by the arrows, are missing: other calculations show that the j = 5 resonance is also missing and theory suggests that the j = 2-4 resonance islands are too narrow to affect the ionisation probability, that is ∆I e < I c e − I e (0), as discussed in the previous section. In s , the static field at which P i is largest, are computed using a grid 10 −5 in F s , and F (j) s is taken to be at the maximum value of P i on this grid. The value of F (j) s is the root of g(F s , F µ ) = jΩ 0 /3: below the double lines F (j) s + F µ > F rc , the approximate radius of convergence of the series for g, see equation 24. The function g(F s , F µ ) is known only via its series expansion, equation 25, which has been computed to O(F 17 ), that is the first nine terms. From the discussion after equation 24, since max(F ) = F s + F µ , we expect any theory based on the series representation of g to be valid for those resonances satisfying F (j) s + F µ < 0.18 (I m = 0) and 0.22 (I m = 1).
It transpires that if F is near the upper boundary the values of F (j) s are sensitive to the number of terms in the series for g and extrapolation is necessary to estimate converged values. Here we consider two methods of extrapolating and give reasons which suggest that the Padé approximant is more reliable. All the following results are obtained by substituting F µ = 0.13 and I m = 0.2 into the series for g and then manipulating the resultant power series in F s : for completeness we give this series: g(F s ) F s = 1.057 + 0.08759x + 0.1079x 2 + 0.1073x 3 + 0.05903x 4 + 0.01588x 5 +0.001861x 6 + 7.700 × 10 −5 x 7 + 6.509 × 10 −7 x 8 , x = (10F s ) 2 . (36) For the j = 7 resonance the lowest order approximation gives F (7) s ≃ 7Ω 0 /3 = 0.02663, which is about 7% too large. Eight other estimates can be obtained by truncating the series for g(F s )/F s at F 2k s , k = 1, 2, · · · , 8: these are 0.02566, 0.02536, 0.02522, s , using p = 1, 2, · · · , M terms of the series 36, then using the method of Richardson we assume that F p = F s . Another approach is to form a Padé approximant of g(F s )/F s using the expression 36, treated as an eighth degree polynomial in F 2 s . The coefficients, and hence the positions of the poles, of these approximants depend upon I m and F µ : we find that the position of the pole nearest the origin is relatively insensitive to I m , but changes significantly with F µ . Therefore for F µ < 0.09 we use a [2/2] approximant (in F 2 s ) and for F µ ≥ 0.09 we use a [3/3] approximant.
For the case considered here, F µ = 0.13, I m = 0.2, the relevant Padé is s . In this case the Padé approximant seems to provide a more reliable method of extrapolating the truncated series for g(F s ). We note that this resonance field is on the edge of the validity of the series expansion, so any estimate of F (15) s based on the series may not be accurate; in these circumstances, however, the Padé approximant is more likely to provide an accurate estimate of the exact function. In all future estimates of F (j) s we therefore use the Padé approximant. We now turn our attention to the resonances missing from figure 12. For this analysis we use the resonance Hamiltonian used to plot the contours in figures 6-10 which, for the reasons discussed at the end of the previous section, provides only a qualitative description.
Using the simple approximations for g andg 1 , equations 25 and 26, we find than the width of the resonance island is proportional to An overview of the widths of the j = 2-15 resonances is given by the graph of this function with j taking all real values. This graph is shown next, with integer values of j being marked by the circles. decrease with increasing j; for j = 2, I c e = 0.75 and for j = 5 I c e = 0.43, but for these cases the maximum possible size of the resonance island is smaller than I c e − I e (0). For j = 6, 11 and 14 the resonance island is seen from figure 13 to be very small so these resonances do not affect the ionisation probability.
For j = 7-10 the simple Hamiltonian with contours shown in figures 6-10 suggests that the resonance island is slightly too small for enhanced ionisation. But non-adiabatic effects, the approximations used -recall that two averaging approximations have been used to derive this Hamiltonian -and the field envelope will broaden these boundaries. The j = 11 and 14 resonance islands are predicted to be too small to promote ionisation, whereas this simple approximation predicts enhanced ionisation for all other j.

Resonance disappearance
The jth resonance has no effect on the dynamics at those values of F µ where J j = 0, see equation 32. With g(F s , F µ ) approximated by a Padé approximant, as in equation 36, the equation g = jΩ 0 /3 provides an expression for F  Observe that for j = 6, 11 and 14, F (j,k) µ ≃ 0.13 and that for j = 8 it is close to 0.13 and at this resonance max(P i ) is relatively small.
In the experimental results reported by Galvez et al (2000) it was shown that resonances disappear at certain field values (F (j) s , F (j,k) µ ), given approximately by equation 33. Since then two sets of more accurate measurements have been made. First, disappearances in the 8.105 Ghz cavity, with scaled frequencies mostly in the range 0.0731 ≤ Ω 0 ≤ 0.136, are reported by Galvez et al (2004). Second, Schlultz (2003) has reported results for a cavity with frequency 3.5539 GHz, and scaled frequencies in the range 0.035 < Ω 0 < 0.16, see also Schultz et al (2004). In those papers experimental values of (F For the 29 cases considered by Schultz (2003) the relative difference between the experimental and computed values of F (j) s is less than 1% for 22 cases and between 1% and 2% in 5 cases: for F (j,k) µ the corresponding percentages are 11 and 16 respectively. For a comparison with theory we compute (F ). Of the 32 cases considered the relative difference between F (j) s and F (j) s is less than 1% for 18 cases and between 1% and 2% in 10 cases: for F (j,k) µ the corresponding percentages are 24 and 6 respectively.
Some typical comparisons between the Monte Carlo calculations and the theoretical values, for various scaled frequencies, are shown in table 4.  with I m . From the series for g(F ), equation 25, we see that F (j) s decreases as I m increases, but that the difference between the largest and smallest values is approximately 11jΩ 0 (9F 2 µ + (jΩ 0 ) 2 )/432: for the data in table 2 this gives 4.4j × 10 −5 approximately. The results obtained using the Padé approximant for g are given in table 5. Columns 2 to 4 show the mean F (j) s , averaged over I m , and the minimum and maximum values of F (j) s , for F µ = 0.13 and Ω 0 = 0.0114. This data shows that, for this low frequency, the above simple estimate of the spread is reasonable for j ≤ 10, and that it is comparable with the difference between the F (j) s and F (j) s , given in columns 2 and 3 of table 2 and to the resonance width, column 5. Note also that as j increases, and F (j) s + F µ tends towards the radius of convergence of the series for g, the difference decreases: it is not known if this effect is real or due to the approximations used.
Although the position of the dynamical resonance is very weakly dependent upon I m its effect on P i (F s ) can depend strongly upon I m , because I c e -that is the solution of F crit (I c e , I m ) = F µ + F (j) s -depends upon I m : for instance with j = 10, I c e varies from 0.04 (I m = 0) to 0.3 (I m = 1).  Table of the theoretical resonance positions for the parameters defined in table 2

Resonance widths
The shapes of the resonances seen in the experimental data and the classical simulations are complicated. In particular the resonances are not normally symmetrical about F (j) s , with details depending upon F µ and j. Here we show show the classical widths can be defined and computed and discuss some of the reason for the shapes observed. A detailed comparison between the classical simulations and the experimental data is provided in the companion paper, Galvez et al (2004), and there it is shown that the experimental results display the same complexities.
In order to estimate the widths it is necessary to isolate the resonances from the background. This is achieved by noting that on either side of the resonance P i (F s ) increases approximately linearly, figure 14 below. The background may be eliminated by subtracting these straight line segments from P i (F s ) to give an adjusted probability that is approximately zero on both sides of the resonance, as shown on the right of figure 14. The only complication with this procedure is that the straight line segments have different gradients, so we form a new fit to the background, P fit i = m(F s )F s +c(F s ), where the gradient, m(F s ), and constant, c(F s ), change smoothly between the values either side of the resonance, (m 1 , c 1 ) and (m 2 , c 2 ) respectively. If the straight line segments are on the intervals F 1 ≤ F s ≤ F 2 and F 3 ≤ F s ≤ F 4 , with F 2 < F 3 , (chosen by eye) we set with a similar fit for c(F ). In the following two figures we show how this process works for the j = 2 resonance with n 0 = 47 (Ω 0 = 0.1278) and F µ = 0.1, which is a fairly typical example: in these calculations a microcanonical distribution of initial states with 1296 orbits, for each value of F s , is used and the envelope is 16-113-16. The left hand panel of figure 14 shows that either side of a resonance P i (F s ) is approximately linear, but with different gradients. The difference between these gradients changes with F µ and, of course, is zero when the resonance disappears. The right hand panel shows the graph of P i (F s ) − P fit i (F s ) which highlights the resonance shape; this is clearly asymmetrical about the maximum. The graph shown is typical though the degree of asymmetry changes with j and F µ . The position of the resonance, the two gradients m 1 and m 2 and the width of the adjusted ionisation probability all provide tests for any theory. Comparisons of these parameters obtained from classical calculations and experiment are given in Galvez et al (2004).
In this particular example F (2) s = 0.0802: the width at half-height is about ∆F s = 0.0018, with the base being nearly 4 times wider, 0.007. The calculated position of this dynamical resonance varies from F (2) s = 0.0797 (I m = 1) to 0.0800 (I m = 0). The spread in F (2) s due to the variations in I m is therefore about 0.0003, about 1 6 of the 1 2 -height width seen in figure 14: this difference is fairly typical.
We now consider some of the factors determining the resonance shapes, and show how these are partly determined by a combination of substate averaging and non-adiabatic effects. Consider ionisation from a given substate: for illustrative purposes choose n 0 = 47 (Ω 0 = 0.1278), I m = 0, F µ = 0.1 and use a 4-50-4 envelope with initial conditions I e (0) = −0.9, −0.8, −0.6 and −0.4. The values of F crit − F µ are depicted by the arrows in each of the four graphs and adiabatic invariance suggests that P i = 1 to the right of these arrows,   This data also shows that as I e increases and the adiabatic boundary encroaches upon each resonance it broadens, acquires an asymmetry and eventually disappears. In other cases, when F µ > F s , an ionisation boundary can also encroach from the left, figure 5, and this will also change the resonance shape. This example, which is typical, suggests that the width and shape of the microcanonical averaged resonances is caused mainly by the effect of the separatrix between bound and free motion, which distorts nearby resonances, rather than the variation in the resonance position with I m . It is therefore difficult to provide theoretical estimates of the resonance width and shape.

Resonance time-scales
The classical adiabatic ionisation mechanism, described in section 3.1, suggests that, in the absence of a resonance, ionisation occurs when F = λ(t)(F s + F µ ) reaches a critical value defined by the condition F = F crit (I e , I m ): at this time ionisation from a particular orbit occurs within a Kepler period. This behaviour has been checked numerically when F s = 0 (Rath, 1990) and the results of the present calculations, where F s > 0, show the same behaviour. With increasing scaled frequency the dynamics becomes less adiabatic although this behaviour persists, albeit with the boundaries becoming blurred; this behaviour is seen clearly in figure 2 of Richards (1996a).
The classical ionisation mechanism at a resonance is different because here for an orbit to be ionised it must first be transported into a region of large I e , that is smaller F crit , by motion round the resonance island. Hence the rate of ionisation will depend upon the period of this motion; at the island centre this is given approximately by equation 35, which shows the period to be O(F −3/2 Ω −1/2 0 ). Thus we should expect the time dependence of the ionisation probability on and off resonance to be quite different.
These predictions can be checked by computing the time at which P i (F s , t) reaches a given proportion of its final value. In the following figure we show two graphs which allow comparison of this ionisation time with the ionisation probability. The upper graph is the ionisation probability, P i (F s ), for Ω 0 = 0.0528 (n 0 = 35), F µ = 0.13 starting in the initial state (I e , I m ) = (−0.4, 0.2), using a 16-50-16 envelope and 1600 orbits, which is the same as in figure 2: the lower graph is the time, T h , at which P i (F s , t) reaches half its final value: with T h is measured in units of the field period: the horizontal line is at T = 16T f , the time when the field amplitude reaches its maximum. There are several features of this comparison worthy of note.
(i) For F s > 0.06 ionisation occurs close to end of the switch-on time, T = 16T f .
Since F crit (−0.4, 0.2) = 0.198 the adiabatic condition suggests, that away from resonances and these initial conditions, bound states exist only for F s < 0.07.
(ii) At the j = 1, 2 and 3 resonances, ionisation occurs some time after the field has reached its maximum amplitude, with the longest delay occurring at the edges of the resonance and the shortest near the maximum in P i . This is consistent with the description given in section 4, where it is shown that close to the resonance edge transport is near the separatrix where motion is slowest. The formula 35 for ω j , gives, for these parameters, the 1/2-period near the resonance island centre of about 14T f , which is consistent with the lower graph of figure 17.
(iii) The local maximum in T h at F s ≃ F (2) s cannot, at present, be explained.
(iv) The ionisation time near the local maximum in P i at F s = 0.0428 has the same shape as those near the j = 1 and 2 resonances, but ionisation clearly takes longer, suggesting that this structure is due to a higher-order resonance. Linear interpolation between the j = 1-3 resonances suggests this could be the j = 2 2 3 resonance. A similar calculation suggests that local maximum in T h at F s ≃ 0.0537 could be the j = 3 1 3 resonance. It is not easy to see what produces these non-integer resonances. Second-order perturbation theory applied to the mean motion Hamiltonian 22 does not appear to give 1/3 resonances; this suggests that higher harmonics of φ e are required and these occur, at this level of averaging, only if higher order terms of θ 1 (ψ, χ) and θ 2 (ψ, χ), equations 12 and 49, are included.
(v) For F s ≥ 0.05 the boundary at F crit −F µ = 0.068 is begining to affect the dynamics and seems to be interfering with the j = 3 resonance.

Envelope effects
In the previous section it was shown how resonance islands affect ionisation times. Here we examine the effect of the envelope switch-on time, T a = 2πN a /Ω, on a particular resonance. At this point it is useful to recall that the dynamical resonances discussed here are unusual because each exist only for a narrow range of F s and within this interval the resonance island moves from the lower to the upper edge of phase space, see figures 6-10. Moreover, the motion inside a resonance island is very slow, equation 35.
s , for some j, then for most of the switch-on period the resonance island does not exist. But for some time close to T a the island develops at the bottom edge of phase space and as t → T a it moves up through phase space and through the initial phase line. As this happens the line is distorted, with the amount of distortion depending upon the relative values of dλ/dt and the frequency of the motion in the island. For short switch times the initial phase line I e (0) evolves into a nearby line at t = T a , as shown in the left panel of figure 20 below. For relatively long switch times there may be enough time for the initial phase line to develop an incipient homoclinic tangle and become quite complicated, as seen in the right panel of figure 20. The examples considered next show how changes in the switch-time, T a , can dramatically affect the ionisation probability.
The demonstration of this effect is in two parts. First we show some exact numerical results illustrating how the ionisation probability changes with T a . In this example we chose the low frequency, Ω 0 = 0.011414, (n 0 = 21), I m = 0.2, F µ = 0.13 and examine P i (F s ) in the vicinity of the j = 7 resonance for the envelope N a -50-N a , for N a = 1-40. From table 2 we see that when N a = 16, P i (F s ) has its maximum at F s = 0.024950 and that the resonance width is ∆F s ≃ 8 × 10 −5 .
In the following two figures we show how P i (F s ) changes with N a : for these calculations we used 900 orbits and a grid δF s = 2 × 10 −6 , so there are 50 data points for each unit of f = (F s − 0.024925)10 4 .
The first figure shows ionisation probabilities for N a = 1-17. The variation of P i with N a shows a surprising amount of variation; in particular we note that for N a = 15 the ionisation probability is zero across the resonance.
Also observe that F (j) s changes by ∆F ∼ 5.2 × 10 −5 for N a in this range, and that this is comparable to the resonance width. This explains why an unambiguous relation between F In the next two figures are shown the ionisation probabilities for 23 ≤ N a ≤ 40. With these longer switch times more structure is seen. For instance, with N a = 31 and 36, P i (F s ) has two local maxima and for both N a = 31 and 37 the probability has a long, low plateau after the maximum. A qualitative explanation of these feature is given next. The behaviour depicted in figures 18 and 19 can be understood qualitatively using a combination of the mean-motion Hamiltonian, equation 39 below, and the resonance Hamiltonian. We assume that initially the system is in a given I e -state with its conjugate variable uniformly distributed in (0, π). This initial phase line, C 0 , evolves during the switch-on period in the mean-motion Hamiltonian: in scaled units, with I n = 1, this is where g(t) = t 0 dt λ(t) (F s + F µ cos Ω 0 t). For t ≥ T a (and before the switch-off time), λ = 1 and this Hamiltonian simplifies to sin Ω 0 T a . Near the jth resonance this can be approximated by the resonance Hamiltonian, where α j is defined after equation 32 and 2θ e = 2ψ e − (3F s − jΩ 0 )t − 3g a , (t ≥ T a ). During the period 0 ≤ t ≤ T a the initial phase curve evolves according to the Hamiltonian 39 into the line C a . Hence by plotting the line C a and the contours of K j (the dashed lines), we obtain a qualitative picture showing how P i can be affected by the field switch.
In the following three figures we show the separatrix of K j and the line C a for F s = 0.0263, F µ = 0.14, Ω 0 = 0.0114, I m = 0.2 and the initial state I e (0) = −0.3, when N a = 1, 3 and 35: the field values are slightly different from those used to generate figures 18 and 19 because an approximate Hamiltonian is used. For N a = 1 the line C a is close to the initial line, I e (0) = −0.3. Slightly less than half of C a lies inside the separatrix and these orbits will be transported to regions of larger I e and some will ionise depending upon the value of I c e . For N a = 3 the curve C a is more distorted; a smaller proportion of orbits lie inside the separatrix, but all of these are close to it and will all be transported to larger I e than the equivalent points of the previous example.
For N a = 35, C a has developed a complicated shape due to the motion inside the island. In this example a significant proportion of the orbits inside the separatrix are close to the horizontal line through the island centre, so will not be transported to regions of large I e . Clearly this structure is very sensitive to changes in N a and this sensitivity will be reflected in P i .
These figures provide a qualitative explanation for the complications seen in figures 18 and 19. In particular they show why there is no simple, precise relation between F (j) s and F (j) s ; they also show that the dynamics underlying the apparently simple resonances seen in figure 1 is very complicated.

Conclusions
In this paper we examine the behaviour of a classical hydrogen atom in parallel static and microwave fields, with frequencies that are low by comparison to the unperturbed orbital frequency. There are three main reason why the classical atom is considered.
First, it is not necessary to make dynamical approximations in order to numerically integrate the classical equations of motion. The errors of the estimated ionisation probability are determined mainly by the Monte-Carlo sampling errors and with modern computers these can be made acceptably small. This is in contrast to quantal calculations in which unquantifiable approximations have to made in order to solve the corresponding equations of motion.
Second, within the framework of classical dynamics there is a range of easily applied approximations that help provide understanding of observed phenomena. The corresponding approximations are not so easy to apply to either Schrödinger's of Heisenberg's equations of motion.
Finally, using techniques of analytical dynamics it is possible to construct an approximate Hamiltonian, which provides a fairly accurate approximation to the exact classical dynamics, see figure 5, and which may be used as a basis for feasible quantal calculations.
The main effects of interest here are the resonances between the microwave field and the Stark frequency induced by the static field. These resonances were first observed by Galvez et al (2000) and this paper also presents the first theory to describe these resonances qualitatively. It was shown how these resonances are responsible for an enhanced ionisation signal over a narrow range of static field strengths, for fixed microwave field amplitude and frequency. Additionally these signals disappear at particular combinations of the two field amplitudes and recent experiments, Schultz (2003), Galvez et al (2004) and Schultz et al (2004), have extended the measurements of these "disappearance fields".
Since the first observations of these resonances three theoretical papers describing the phenomena from different perspectives have been published. Oks and Uzer (2000) use a Floquet analysis to derive zero-order estimates of F (j) s and F (j,k) µ . Robicheaux et al (2002) have solved Schrödinger's equation for this problem using a split-operator method and have made comparisons with the classical and experimental ionisation probabilities for n 0 = 39 across the j = 1 resonance with F µ = 0.144. These calculations suggest that the classical and quantal values of F (1) s are very close, but that quantal value of P i is smaller than the classical value for F s > F (1) s . In addition the time-dependence of the ionisation probability is described for three values of F s and this appears to contradict the results summarised in figure 17, though it is difficult to make comparisons between substate averaged and unaveraged data. Ostrovsky and Horsdal-Pedersen (2003) use an energy shell subspace with a time-dependent electric field and weak, perpendicular magnetic field, with the aim of understanding oscillations seen in the experimental result, subsequently attributed to another cause, see Wilson et al (2004). This analysis is based on the same type of averaging approximation that leads to the Hamiltonian 22, see also Born (1960, section 38), and inevitably gives zero-order estimates of F (j) s and F (j,k) µ .
In this paper we have described the classical dynamics of this system in more detail: in particular more accurate values of F (j) s and F (j,k) µ are determined and the properties of the classical resonances are described in some detail. We now list these features and discuss the probable consequences to the quantum mechanics. can depend on the field envelope.
2) We have isolated the terms in the Hamiltonian that give rise to the dynamical resonance. This allows the computation of the dynamical resonance position using high-order perturbation theory; where this series converges we obtain improved estimates of the resonance position.
This analysis is essentially the same as the corresponding quantal theory which we therefore expect to give the same result. Further, this suggests that the discrepancy noted in (1) above will also occur in an accurate quantal calculation.
3) Using a classical approximation, based on two stages of averaging, we have derived a number of conditions necessary for the dynamical resonance to affect the ionisation probability. These depend upon properties of the classical resonance island, the most significant being the island width.
A dynamical resonance affects the ionisation probability only if it is wide enough to bridge the gap in phase space between the initial state and those states that ionise, see section 4.2. Because the ratio of these two actions is independent of the initial principal number, we expect a similar story in the quantal description, though quantal effects will inevitably blur these boundaries. This suggests that there may be cases where a resonance not seen in the classical ionisation probabilities will be visible in the quantal probabilities.
Besides the island width, its area also plays a role in quantum mechanics; this area is proportional to the initial principal quantum number, n 0 , so we expect the resonances seen here and in current experiments to change, and possibly disappear, as n 0 decreases.
4) The resonances in the ionisation probability, not averaged over substates, are generally very sharp, see for instance figure 12 and table 2. The full-width of an isolated resonance is generally smaller than its theoretical shift produced by changing I m = mh, as seen by comparing the data presented in tables 2 and 5. However, resonances near an adiabatic boundary are significantly broadened and asymmetries are introduced, figure 15 and 16. This causes substate-averaged resonances to be far wider than isolated resonances, and also affects their shape.
The classical dynamics of this process is complicated and not understood. Because the experimentally observed resonances behave in a similar fashion, we expect a similar dynamical effect in any quantal calculation, but understanding this is harder than understanding the classical dynamics.
5) The shape of an isolated classical resonance and the value of F (j) s , with no substate averaging, can be affected by the field envelope, if the fields are switched on sufficiently slowly, see figures 18 and 19. These changes are caused by the phase line representing the initial state becoming tangled as the separatrix of the resonance island passes through it, see figure 20. We expect this behaviour to be seen in the quantum dynamics provided the principal quantum number is large enough, but how large is not known.
6) Classical resonances develop over a time scale that is much longer than that of the ionisation process operating away from resonance. This is because, at resonance, ionisation happens by transport around the resonance island and this is a relatively slow process. Further, across a resonances the ionisation time appears to reflect the island dynamics. For instance, at the edge of a resonance the ionisation time is longest because motion near the separatrix is very slow: this and other, not understood, features are seen in the lower panel of figure 17.
As with point 5 above, we should expect to see similar behaviour in a quantal calculation, provided the principal quantum number is sufficiently large, but how large is not known.

Acknowledgements
It is a pleasure to thank Professor Peter Koch for valuable discussions and encouragement during all stages of this work, and for helpful comments on drafts of this manuscript.

Appendix: action variables
The derivation of the required results is easiest if scaled variables are used. If a is any scale length and I 0 an action, suitable scaled variables arẽ Taking I 0 to be the initial value of I n and a the semi-major axis of the initial Kepler ellipse, a = I 2 n /µe 2 we havẽ F = I 4 0 µ 2 e 6 F,Ẽ = I 0 µe 4 E,α = α.
In the following these scaled units are used but for clarity the tilde is not shown. Most of these results are obtained using a Maple program to manipulate the series, which were computed to higher-order than quoted here.

Series for I 1
For I 1 set ξ 2 = y so equation 8 becomes I 1 = 1 2π y2 y1 dy y f 1 (y), f 1 (y) = −F y 3 − 2W y 2 + 2α 1 y − I 2 m , where W = −E > 0. In the parameter range of interest f 1 (y) has three real roots, two positive 0 ≤ y 1 < y 2 (with y 1 = 0 only when I m = 0) and one negative root, y 3 = −y − /F < 0. We may write f 1 (y) = (F y + y − )(y − y 1 )(y 2 − y), y 1 y 2 y − = I 2 m , (y 1 + y 2 )F = y − − 2W, so that y − = O(1) and, if F = 0, y − = 2W . Only the value of y − is needed because all quantities of interest can be expressed in terms of the combinations y 1 + y 2 and y 1 y 2 . Then the series for I 1 is where H k = 1 4π C dz z k−1 h(z) and a k = √ π 2 k! Γ(3/2 − k) , and where h(z) = − (z − y 1 )(y 2 − z) has a cut on the real axis between y 1 and y 2 such that on the real axis between these two roots on the upper branch h(x) < 0, and on the lower branch h > 0; for x real and x > y 2 , h(x) = i (x − y 1 )(x − y 2 ) and for x < y 1 , h(x) = −i (y 1 − x)(y 2 − x). The contour C encloses the branch cut between y 1 and y 2 but not the origin. If k = 0 the integrand has a pole at z = 0 and contributions from the circle z = Re iθ , as R → ∞, so For k ≥ 1 the only contribution is from the circle at infinity, C ∞ H k = i 4π C∞ dz z k 1 − y 1 z 1 − y 2 z = (−1) k 2 k+1 r=0 a r a k+1−r y r 1 y k+1−r 2 , k ≥ 1.
Thus y 1 + y 2 = α 1 W + F 4W 3 W I 2 m − 2α 2 1 − α 1 F 2 8W 5 3W I 2 m − 4α 2 1 + · · · . and 16W 5 W I 2 m − 4α 2 1 F 2 + · · · . Substituting these expressions for y 1 + y 2 and y 1 y 2 into the above expressions for H k gives Series for I 2 For I 2 we set η 2 = u, and equation 9 gives In the parameter range of interest f 2 has three real positive roots, 0 ≤ u 1 ≤ u 2 and u + /F , so we write f 2 (u) = (u + − F u)(u − u 1 )(u 2 − u) with u 1 u 2 u + = I 2 m and u 1 + u 2 = (2W − u + )/F . Thus where a k and H k are defined above. The perturbation expansion for u + is given, as before, by setting u = z/F to write the equation f 2 = 0 in the form z 3 − 2W z 2 + 2α 2 F z − F 2 I 2 m = 0, giving This then gives u 1 + u 2 = α 2 W − F 4W 3 W I 2 m − 2α 2 2 − α 2 F 2 8W 5 3W I 2 m − 4α 2 2 + · · · . and u 1 u 2 = I 2 m 2W + α 2 I 2 m 4W 3 F − I 2 m 16W 5 W I 2 m − 4α 2 1 F 2 + · · · and hence the expression for I 2 is, These series for I 1 and I 2 now need to be inverted to give α 1 , α 2 = 2 − α 1 and W as power series in F . The zero-order term is trivial so we substitute the series into the series for I 1 and I 2 and solve for the unknown coefficients to give the energy and the separation constant quoted in equations 10 and 11.
where σ k = I k (I k + I m )/I n , k = 1, 2. When F = 0 these expression lead to the formulae quoted in equation 14.

Evaluation of ∂S/∂F
The generating function S 1 (θ 1 , θ 2 , ξ, η) returns to its initial value when either ξ or η increases through a period, see Born (1960, page 82); here we use the notation of Goldstein (1980) to label generating functions. It follows that the Hamiltonian 15 is a periodic function of the angle variables with zero mean value. By differentiating the generating function S(I, ξ, η) with respect to F and using the angles (ψ, χ), defined in equations 47 we see that ∂S/∂F may be written in the form where W F = ∂W/∂F and α 1F = ∂α 1 /∂F . The two integrands can be expressed in terms of even Fourier series in ψ and χ respectively where the constant terms are missing because of the argument given at the beginning of this section and because of the general result ∂S/∂t = ∂S 1 /∂t, which gives ∂S/∂F = ∂S 1 /∂F . Hence the integral 51 leads to the Fourier series defined in equation 16, where the coefficients (A k , B k ) depend upon F s , F µ and the action variables. Because of relation 49 it follows that no term ofḞ ∂S/∂F is independent of the angle variables. The Fourier series representation of ∂S/∂F is obtained in the same manner as that for P k (ψ) and Q k (χ): to O(F ), we have A 1 = 1 2 I 4 n σ 1 (I 1 + 3I 2 + 2I m ) − 1 4 I 3 n F 14I 2 1 + 27I 1 I m + 26I 1 I 2 + 22I 2 2 + 35I 2 I m + 22I 2 m , A 2 = − 1 4 I 5 n σ 2 1 1 − F I 3 n (4I 1 + 6I 2 + 5I m ) , A 3 = − 1 12 F I 9 n σ 3 1 , B 1 = − 1 2 I 4 n σ 2 (3I 1 + I 2 + 2I m ) − 1 4 I 3 n F 22I 2 1 + 35I 1 I m + 26I 1 I 2 + 2I 2 2 + 15I 2 I m + 22I 2 m , B 2 = 1 4 I 5 n σ 2 2 1 + F I 3 n (3I 1 + 7I 2 + 5I m ) , B 3 = − 1 12 F I 9 n σ 3 2 .