Evolution to Mirror-Symmetry in Rotating Systems

: A natural example of evolution can be described by a time-dependent two degrees-of-freedom Hamiltonian. We choose the case where initially the Hamiltonian derives from a general cubic potential, the linearised system has frequencies 1 and ω > 0. The time-dependence produces slow evolution to discrete (mirror) symmetry in one of the degrees-of-freedom. This changes the dynamics drastically depending on the frequency ratio ω and the timescale of evolution. We analyse the cases ω = 1,2,3 where the ratio’s 1,2 turn out to be the most interesting. In an initial phase we ﬁnd 2 adiabatic invariants with changes near the end of evolution. A remarkable feature is the vanishing and emergence of normal modes, stability changes and strong changes of the velocity distribution in phase-space. The problem is inspired by the dynamics of axisymmetric, rotating galaxies that evolve slowly to mirror symmetry with respect to the galactic plane, the model formulation is quite general.


Introduction
There are many non-autonomous dynamical systems with evolutionary aspects but the necessary theory for understanding them is still restricted to linear systems with equations that have a limit as time increases (see [1] Chapter 3). The main reason for this is probably that understanding steady state, autonomous, nonlinear systems is already a formidable task. Physical examples where explicit time-dependence plays a part are the evolution to spherical structures in nature and on solar system scale planetary satellite systems evolving by long-term tidal forces to more symmetric orbits (for references see [2]). Explicit time-dependence also is ubiquitous in engineering systems but there it is usually concerned with short-time transient phenomena or time-periodic influences.
We will consider Hamiltonian systems where explicit time-dependence slowly vanishes producing the transition from a system that is asymmetric in some sense to a Hamiltonian system with symmetry. In [3] a 1-dimensional cartoon problem is considered where time-dependent asymmetric terms vanish with time; it is of the form: x + x = a(εt)x 2 . (1) A dot means differentiation with respect to time, 2 dots second order differentiation with respect to time, ε is a positive small pameter. The function a(εt) is monotonically decreasing from the value 1 to zero. Equation (1) models slow evolution to symmetric dynamics in the sense that the asymmetric potential vanishes with time; although the equation is simple it shows already a relation with dissipative systems and the construction of a global adiabatic invariant.
As before the asymmetric term vanishes slowly, the system evolves to a 2 -dimensional harmonic oscillator. Again a relation with dissipative systems can be established, 2 adiabatic invariants can be found. The system displays overall dynamics that keeps some information from its asymmetric past.
In the sequel we will study evolution in two dof Hamiltonian systems keeping nonlinearities. That makes the models both more realistic and more complicated.
In this paper we describe an evolutionary system (14) inspired by the evolution of a rotating galaxy. It is important to note that this system in its general formulation of a time-dependent Hamiltonian system extends to many other physical problems apart from those derived from this particular application.

Set-Up of the Paper
In Section 2 we describe the dynamics of rotating stellar systems and we formulate the equations for the axisymmetric case. These formulations have been known a long time in astrophysics but we summarise this for a more general audience. In the sequel we describe the dynamics of prominent resonances producing periodic solutions, adiabatic invariants and stability changes during evolution.
In Section 3 we formulate the time-dependent Hamiltonian where the terms that are not discrete (mirror) symmetric with respect to the galactic plane are slowly vanishing with time. This leads to the slow-fast system (14) where the fast motion is described by a family of discrete symmetric two dof freedom Hamiltonian systems with resonances.
Section 4 describes the reformulation of system (14) suitable for averaging-normalisation while using a modification of averaging developed in an appendix (Section 9). It is shown that the system for evolution to symmetry is equivalent with a dissipative system.
The prominent resonance 1:2 is studied in Section 5. The results are surprising as first order averaging-normalisation producing approximations valid on an interval of size O(1/ε) yields no epicyclic normal mode and an unstable vertical normal mode. Second order averaging-normalisation producing approximations valid on an interval of size O(1/ε 2 ) show that after a long time the epicyclic and vertical normal modes arise and are stable. These changes of presence and stability have drastic consequences for the evolution of the overall dynamics of the system.
In contrast to the preceding case, the 1:3 resonance is dynamically less interesting (Section 6). We have in this case basically the dynamics of the 2:6 resonance that is described by higher order resonance techniques from [4].
In Section 7 the dynamics of the 1:1 resonance turns out to be the most complicated case. Second order averaging-normalisation is necessary to describe the dynamics. Depending on the system parameters we find again normal modes and stability changes during evolution. Explicit examples include the evolution to the Hénon-Heiles potential.
After the concluding discussion we present in Section 9 a brief survey of the averaging theorems we have used in this paper. For the proofs we indicate the relevant literature.
Numerics. For the computation of the figures we have used MATCONT (free access) under MATLAB (licensed). The algorithm used is ode 78 with relative and absolute errors exp.(−15).

The Dynamics of the Milky Way
In this section we will summarise the formulation for rotating stellar systems and we will specialise this to the axisymmetric case. The modelling was used in [5] and can also be found in [6].

The Collisionless Boltzmann Equation
The reasonable looking approximation of few collisions in stellar systems raises questions on the statistical mechanics of these systems. Another basic assumption often made in model building for galaxies is most of the systems being in equilibrium. A more realistic assumption is quasi-static equilibrium that leaves room for evolution to this stage.
In this paper the choice of a time-dependent Hamiltonian will be the basis for studying such models of evolution and the quite remarkable consequences.
The collisionless Boltzmann or Liouville equation describes the distribution of n particles, stars in this case, in 6n dimensional phase-space. One has to add the Poisson equation for self-gravitation. A good description with examples can be found in [6] Chapter 4. The equation for the distribution of particles f (t, x, v) where x indicates the position, v the velocity, is the continuity or Liouville equation in 6n dimensional phasespace with clearly the assumption that no particle escapes, is destroyed or is created in the system. With Lagrangian derivative d/dt the Liouville equation is d f /dt = 0. If the collective gravitational potential ruling the dynamics of the system is Φ(t, x), the Liouville equation becomes explicitly: The characteristics of this first order partial differential equation are given by the Hamiltonian equations of motion. The solutions of the Hamiltonian system produce according to Monge the solutions of the Liouville equation as they represent the geometric sets where the solutions of the Liouville equation are constant. This procedure is very effective if we find not only solutions but integrals of motion that contain sets of solutions. Assuming for instance a time-independent potential and axi-symmetry, we have already 2 independent integrals of motion, energy E and angular momentum L with respect to the axis of rotation. Any differential function f of E and L will satisfy the Liouville equation. It was noted very early that such solutions will produce velocity distributions that are symmetric perpendicular to and in the direction of the rotation axis. This does not agree with observations in our own galaxy, in fact the two dispersions differ considerably; for velocity dispersions in the plane of the galaxy see [7], for dispersions in the halo [8]. Velocity dispersions in galaxies is an ongoing research topic, the reader interested in recent observational results might google the authors Rachel Bezanson, M. Franx and P.G. van Dokkum.
The velocity asymmetry triggered off a long search for "a third integral of the galaxy" that came to an important point in the 1960's when it was shown in Arnold-Moser theory that such 3rd integrals in general do not exist but that an infinite number of so-called KAMtori corresponding with approximate first integrals may survive. This near-integrability caused a lot of confusion in the literature.
One of the aims of our study is to investigate whether evolution from an asymmetric state to a near-integrable symmetric dynamical state might influence the velocity distributions as if the system "remembers" its asymmetric past.
From [6] Chapter 4 we consider the dynamics in cylindrical coordinates R, φ, z witḣ Solving system (4) and obtaining (approximate) first integrals enables us to describe families of distribution functions that satisfy the Liouville Equation (3).
In the sequel we assume that the systems to be considered have reached the axisymmetric state. Initially the systems are not mirror-symmetric with respect to a plane, they evolve slowly in time to such a state.

Rotating Axi-Symmetric Models
We will consider models that are inspired by axisymmetric rotating galaxies. One can think of disk galaxies or rotating flattened elliptical galaxies. The axisymmetry expressed by produces the angular momentum integral J (in [6] called L z ) enabling us to reduce spatial 3-dimensional motion to 2 degrees of freedom (dof). We have with Equations (4) and (5) leading to the angular momentum integral The equations of motion (4) can be written as In the equatorial plane that is perpendicular to the axis of rotation we have circular orbits at R = R 0 where the rotation speed matches constant angular momentum: We will expand the Hamilitonian around the circular orbits. In most models, the potential Φ is assumed to be mirror (discrete) symmetric with respect to the equatorial plane. In cylindrical coordinates R, z with z in the direction of the rotation axis and z = 0 corresponding with the equatorial plane we will put in a final stage of evolution Φ = Φ(R, z 2 ). The evolution towards this symmetric state is caused by mechanisms unknown to us, maybe contraction combined with rotation or dynamical friction plays a part. We propose to avoid the speculative description of complicated mechanisms by introducing a function of time slowly destroying the asymmetric potential.
Around the circular orbits in the galactic plane we find epyciclic orbits in the Rdirection nonlinearly coupled to bounded vertical motion in the z-direction. An early study of such orbits in a steady state galaxy model is [5], for a systematic evaluation of the theory see [6]. A detailed analysis of the orbits can be found in [9]. A simplifying transformation from [5] is to use the reduced potential leading to the equations of motion:

Shortcut to a Model with Evolution
We want to describe the dynamical consequences of evolution to an axi-symmetric state by expanding around the circular orbits at position (R, z) = (R 0 , 0) putting R = R 0 + x. To study evolution to symmetry we consider as a model the time-dependent two dof Hamiltonian: The epicyclic frequency in the galactic plane has been scaled to 1, the vertical frequency is ω, so x corresponds with deviations of R 0 . The function α(δt) is continuous and monotonically decreasing from α(0) = 1 to zero; we will usually choose α(δt) = e −δt with 0 < δ 1. If a 1 = 1, a 2 = −1, a 3 = a 4 = 0 we have the famous Hénon-Heiles problem [10]. A detailed modelling of the evolution of a large contracting cloud to an axisymmetric galaxy is not available. The model presented here ignores all special types of dynamics as for instance the formation of new stars, dissipation by gas cloud motions or the spectacular emergence of high-velocity stars by supernovae in binary stars.
The evolution to a time-independent Hamiltonian is also clear when considering the Lagrangian derivative: The energy given by (10) decreases slowly. Near the circular orbits, the time-dependent perturbation factor decreases monotonically, the change of the energy will be small and depends on the signs of a 3 , a 4 and the behaviour of z(t). On the other hand the quantitative dynamical consequences of the asymmetric potential in the transitional stage can be considerable and depends on the size of a 3 , a 4 and δ.
The quadratic part is called H 2 , the cubic part H 3 . The coefficients a 1 , . . . , a 4 are free parameters with a 4 = 0, ω represents the frequency ratio of the 2-dof oscillators with prominent resonances ω = 1, 2, 3. The possible values of ω depend on the galactic potential constructed. An example describing an axisymmetric rotating oblate galaxy can be found in [6] eq. (3-50), leading to frequency ratios given by eq. (3-56).
A time-independent example where in the centre of the galaxy a very massive nucleus is found produces with mass M large the potential: where, at least near the centre, Φ 1 is small with respect to Φ 0 . Assuming rotation and expanding in a neighbourhood of the centre and near the circular orbits using (4) and (7) we find with The dots represent the neglected Φ 1 and nonlinear terms; this system in 1:1 resonance will dominate the dynamics near the centre of a galaxy with a very massive nucleus.
The terms with coefficients a 1 , a 2 in Hamiltonian (10) are discrete symmetric in z, de terms with a 3 , a 4 are not symmetric in z but the asymmetry vanishes as t → ∞.
The dynamical system induced by (10) is not reversible as time-independent Hamiltonians are, the main question of interest is then whether after a long time the system induced by the Hamiltonian shows traces of the original asymmetry. The answer will be positive and the traces significant.
If δ = 0 (the time-independent case) the origin of phase-space in (x, z) coordinates is Lyapunov-stable, the energy manifold is bounded in an O(ε) neighbourhood of the origin.
Consider now the time-dependent dynamics. To make the local analysis more transparent we rescale the coordinates x = εx etc. Dividing by ε 2 and leaving out the bars we obtain the equations of motion: Because of the localisation near the circular orbits represented by the origin of phasespace we assume that ε is a small positive parameter. The parameter δ will also be small, the ratio of ε and δ influences the dynamics and analysis. For instance if we want to study the influence of dynamical friction in a galaxy we can choose δ = ε n , this implies that for n = 1 and so ε = δ we will consider a very small neighbourhood of the origin. On choosing ε = √ δ or n = 2, the neighbourhood will be larger. Again a larger neigbourhood is obtained for n = 3, these different cases complicate the analysis and will produce different local dynamics.

First Order Averaging-Normalisation
To characterise the dynamics induced by Hamiltonian (10) we will use averagingnormalisation, see [11] and Section 9. We transform to slowly varying polar coordinates r, ψ by: (15) leading to the slowly varying system: Near the normal modes r 1 = 0 and r 2 = 0 we have to use a different coordinate transformation.
We can put τ = δt and treat τ as a new variable. It will also be useful to introduce the actions E 1 , E 2 by: As we shall see in subsequent sections, for each choice of ω ≥ 1 the average of the terms in system (16) with coefficients a 1 , a 2 vanish, so we can use the near-identity transformation (45) of Section 9. The implication is that to first order in ε and on time intervals of size 1/ε system (14) is described by the intermediate normal form equations: We recognize the presence of the z,ż normal mode solution as x =ẋ = 0 satisfies system (18); this can be found for slowly moving variables by using a coordinate system different from polar coordinates.
As the nonlinear terms are homogeneous in the coordinates we can remove the time-dependent term by a transformation involving α(δt). For instance if α(δt) = e −δt we put x = α(δt)y 1 , z = e δt y 2 (note that such a transformation exists for any positive sufficiently differentiable function of time that decreases monotonically to zero). System (18) transforms to: So the time-dependence that is removing the asymmetry transforms to a dissipative system with friction coefficient 2δ. Another interesting aspect of system (19) is that the transformation decouples the parameters ε and δ. This is useful if δ is not a small parameter as we can solve the system for ε = 0 and start perturbing the resulting dissipative system.
We will average the right-hand sides of system (16) over t keeping r 1 , r 2 , ψ 1 , ψ 2 , τ fixed. We haveτ = δ, so to match the size of the other equations of system (16) we choose δ = ε n with a suitable choice of n ≥ 1; for simplicity we restrict n to natural numbers. As stated above, by first-order averaging the terms with coefficients a 1 , a 2 vanish, we can use system (18). The subsequent averaging results depend strongly on the choice of ω. To make the calculations more explicit we put in the sequel On choosing polynomial decrease of α(δt) we would have slower decrease with as a consequence that we have to retain more small perturbation terms.

The 1:2 Resonance
The prominent case for 2 dof conservative systems is the 1:2 resonance (ω = 2). We analyse the time-dependent system for different choices of δ.

First Order Averaging
Averaging system (18), see Section 9, we find: with combination angle χ = 2ψ 1 − ψ 2 and for χ the equation: The vertical z normal mode (x =ẋ = 0) exists and is unstable, see [11]; the epicyclic x normal mode does not exist on the interval of time O(1/ε). We will see what happens on longer intervals of time in the next subsection.
Remarkably, system (47) admits families of solutions with constant amplitude on intervals O(1/ε) if: The solutions with χ = 0 are called in-phase, the solutions with χ = π out-phase. If δ = ε, the corresponding phases are slowly decreasing with rate exp (−εt). If we choose δ = ε n , n ≥ 2, the amplitudes r 1,2 and phases ψ 1,2 will be constant with error O(ε) on intervals of time of size 1/ε, changes will arise on longer intervals of time. System (47) admits 2 time-independent integrals of motion: with constants E 0 , I 3 ; In the original coordinates we have: The solutions and integrals (adiabatic invariants) of system (47) have the error estimate O(ε) on time intervals of size 1/ε. On this long interval of time and longer ones we expect the terms O(ε 2 ) to play a part as the solutions of system (47) with coefficients a 3 , a 4 will vanish and other terms of system (14) will become important.
The result is surprising as we would expect τ-dependent terms for the amplitudes at second order; such terms arise only for the angles. Also combination angles for ψ 1 , ψ 2 are not present at this level of averaging-normalisation. The 2nd order system (25) was computed with time-independence in [4], eq. (4.2); leaving out the time-dependent terms the results agree. The implication is that near the 1:2 resonance the dynamics after mirror symmetry has set in behaves like a system in higher order resonance.
We give a more precise formulation of the interesting implications: 1.
During an interval of time of order 1/ε the integrals (adiabatic invariants) (24) are active, the system is dominated by the asymmetric a 4 term. On this interval of time it will govern the orbital dynamics and accordingly the corresponding distribution function in phase-space. This holds for n = 1, 2, . . .

2.
If δ = ε, on time intervals of order 1/ε 2 (and longer) the time-independent system involving the coefficients a 1 , a 2 dominates the dynamics. In [4] it is shown that for this system, depending on a 1 , a 2 , 2 resonance manifolds can exist on the energy manifold (see also Section 9). Introducing the combination angle: we find according to [4] that on intervals of time asymptotically larger than 1/ε we have: Resonance manifolds exist if the right-hand side of Equation (27) has a zero, the combination angle χ 2 is locally not timelike. In this case the resonance manifold with 4ψ 1 − 2ψ 2 = 0 has stable 2:4 resonant periodic orbits surrounded by tori, for 4ψ 1 − 2ψ 2 = π the 2:4 resonant periodic orbits also exist but are unstable. The resonance manifolds have size O(ε), the dynamics takes place on intervals of time of order 1/ε 3 ; for details on higher order resonance see [4]. Outside the resonance manifolds the dynamics is characterised by the quadratic integrals E 1 , E 2 for each of the 2 dof. For the interaction of the 2 modes on these long time intervals coefficient a 2 is essential, but note that if a 1 = 0, the resonance manifolds do not exist as the right-hand side of Equation (27) has no zero.

3.
An interesting consequence arises when considering orbits that change stability during the process of evolution to mirror symmetry. The instability of the z,ż normal mode persists at second order but leads to stability in the final stage of mirror symmetry after a time dependent on 1/ε n . The implication is that this will generally produce a drastic change of actions. Another qualitative consequence is that the epicyclic x,ẋ normal mode does not exist in the first stage (time order 1/ε) but the normal mode exists in the final stage of mirror symmetry; again we can expect action changes. Similar considerations hold for the unstable out-phase periodic solutions derived before. For illustrations see Figures 1 and 2 and the discussion below. In Figure 1 the corresponding action E 1 (t) is not shown; it exchanges energy met the z-mode.

Consequences for the Distribution Function at the 1:2 Resonance
Suppose we start with a collection of particles (stars) characterised by a distribution function satisfying the Boltzmann equation that is nearly collisionless as it has small dynamical friction added. The system is already in axi-symmetric state but the evolution to mirror symmetry to the galactic plane is still going on. On an interval of time O(1/ε), the first stage of evolution, we have, apart from angular momentum J, two active integrals of motion: E 0 , I 3 and a prominent resonance manifold for each value of the energy described by Equation (23). The family of solutions with constant amplitude (constant to O(ε 2 )) on the energy manifold will be stable for combination angle χ = 0. The distribution function will be a function of the 3 integrals. The velocity distribution v 1 , v 2 and their dispersion will depend on the presence of these integrals.
On intervals of time asymptotically larger than O(1/ε), the primary resonance manifolds described by Equation (23) vanish, they are replaced by the smaller resonance manifolds (of size ε) located by the zeros of Equation (27) and χ 2 = 0, π. The distribution function will now evolve to a function of the actions E 1 , E 2 but with an initial velocity distribution obtained by the resonance on time intervals O(1/ε). See Figure 1 (left) where n = 2 so that the symmetric state develops on intervals of time O(1/ε 2 ). On an interval of time order 1/ε 2 there is first a typical 1:2 resonance exchange between the 2 dof after which the dynamics settles at different actions. If n = 3 we find, as expected, that it takes much longer to settle in a permanent state, see Figure 1. With the choice of ε = 0.1 the effect of evolution to symmetry is much less.
In Figure 2 we have n = 3 so that on intervals of size O(1/ε 2 ) we have still resonant interaction, the symmetric state develops on intervals of size O(1/ε 3 ). The exchanges and time evolution are first much more frequent as it takes longer for the asymmetric terms to vanish. Then the system settles at a much lower oscillation amplitude, in a sense experiencing more of its asymmetric past.
In Figure 3 we show the behaviour (left) near the z normal mode for different values of the energy in the case n = 2. Interesting is that starting at E 2 (0) = 1.624 the resulting values of E 2 are lower than for a few orbits starting at still lower energy. This will be caused by passage of more complicated structures like resonance manifolds in phase-space. Figure 3. Left the behaviour of the action E 2 (t) of system (14) in 1:2 resonance for the case n = 3 and 4 initial conditions, successively E 2 (0) = 2, 1.624, 1.283, 0.982, starting near the normal z-mode, and, E 1 (0) small. For n = 3 it takes longer for the time-dependent interaction to vanish. Velocity v 2 (t) in the case E 2 (0) = 2, v 2 (0) = 0 is given in Figure 2 right.
When the system is close to mirror symmetry, we will find for each value of the energy in a resonance manifold a family of tori around the stable periodic solutions, so for varying energy values this will be a 2-parameter family. Outside these resonance manifolds the orbits will move quasi-linearly, only the phases are position dependent.
The choice of δ, the parameter determining the timescale of destroying the asymmetry of the potential field, together with the energy level will determine the resulting positions and velocities.

The 1:3 Resonance
The 1:3 resonance is for Hamiltonian (10) dynamically very different and less interesting. Most of the analysis can be deduced from [9], we will summarise the results. Putting ω = 3 in (10) we find after 2nd order averaging: The result is quite remarkable as this shows that for the 1:3 resonance the slowly vanishing asymmetry of the potential plays no part for the amplitudes to order 3 in ε. For the phases we find no asymmetric terms either but nontrivial variations depending on the symmetric potential: 12 a 2 1 r 2 1 + ( 1 2 a 1 a 2 − 1 35 a 2 2 )r 2 2 , ψ 2 = −ε 2 ( 1 6 a 1 a 2 + 1 105 a 2 2 )r 2 1 + 23 140 a 2 2 r 2 2 . (29) if δ = O(ε) or O(ε 2 ) the asymmetry plays no significant part. The theory of higher order resonance, see [11] and for extension and examples [4], shows that the combination angle χ 3 = 6ψ 1 − 2ψ 2 plays a crucial role. The equation for χ 3 becomes: Zero solutions of the right-hand side of Equation (30) produce, together with the condition χ 3 = 0, π resonance manifolds of size O(ε 2 ) with characteristic timescale O(1/ε 4 ) for the orbits and motion on the corresponding tori; these estimates follow from the analysis in [4]. With high precision the distribution function will depend on J and the 2 actions.
It is clear that if a 1 = 0, system (29) has no nontrivial zeros and so the resonance manifolds of this type are not present.
The solutions of the averaged system to second order with appropriate initial values produce an O(ε 2 ) approximation on an interval of time order 1/ε of system (16) with ω = 1 but an O(ε) approximation on an interval of time O(1/ε 2 ). The second estimate is a kind of "trade-off" of error estimates valid under special conditions formulated in [12]; see also [11]. It is remarkable that system (32) shows full resonance involving exchange of energy between the 2 dof even when the asymmetric potential terms have become negligible.
We expect the dynamics of the symmetric case to describe the orbits of system (31) on intervals of time larger than 1/ε n . During evolution the in-phase and out-phase solutions are present with constant amplitude and approaching periodicity. The dynamics will after long time be governed by the integrals (35) and (36) producing a complicated velocity distribution and varying actions.
On the long (starting) interval of order 1/ε 2 system (32) has also 2 integrals. As noted before integral (35) holds for all time t ≥ 0.
As the 1:1 resonance occurs in applications and has been studied in many papers we will consider examples in more detail.

Avoiding Degenerate Cases
Special values of parameters may produce degenerations that produce dynamics that is not typical. We illustrate this by looking for solutions of system (31) of the form x = λz with λ a suitable parameter. This assumption implies a degeneration as it means that there is no exchange of energy between the 2 modes. Substituting this relation into system (31) we find the conditions: If conditions (37) are satisfied we find for z the equation: Using 2nd order averaging as in example 12.20 of [13] with n ≥ 2 we find: This type of special solutions shows only variation of the phases. Although interesting in itself, we have to avoid the parameter values of conditions (37) when studying more general behaviour. Other parameter values that may lead to degenerations and should be avoided are the bifurcation values of existence and stability of normal modes and other periodic solutions listed in the preceding subsection.

An Example of Evolution
In Figures 4-8 the choice of a 1 , a 2 gives the system stable epicyclic and stable vertical normal modes if we would have the time-independent case. However, our choice of a 3 , a 4 keeps the normal modes but destabilises them; we have a 1 = −1.5, a 2 = 1, a 3 = 0, a 4 = 4. When the time-dependent perturbation has become negligible after roughly 200 timesteps if n = 2 and 3000 timesteps if n = 3, the orbits have moved into general position. The time-independent case will approximate the dynamics after a long initial period of time, the in-and out-phase periodic solutions exist if a 3 = a 4 = 0, the out-phase solutions are stable.
With the chosen parameter values we find from system (32) for the angles: leading to: For the actions: we have from system (32): with c(t) = 1 6 a 1 a 2 − a 2 2 + exp (−ε n t)( 1 6 a 3 a 4 − a 2 4 ). In our example we have c(t) = − 5 4 − 16 exp (−ε n t).   Parameter values a 1 = −1.5, a 2 = 1, a 3 = 0, a 4 = 4, ε = 0.1. Right we display the corresponding z behaviour by plotting E 2 (t). It takes around 3000 timesteps to settle in the stationary state which is drastically different from the initial state.  Parameter values a 1 = −1.5, a 2 = 1, a 3 = 0, a 4 = 4, ε = 0.1. Right we display the corresponding z behaviour by plotting E 2 (t). It takes around 3000 timesteps to settle in the stationary state which is drastically different from the initial state.
Solutions with constant amplitude are obtained when putting dχ/dt = 0; they are valid with error O(ε) on intervals of time O(1/ε n−1 ). On these intervals time-dependent terms are important. In our example, the in-phase and out-phase periodic solutions with constant amplitude exist if also a 4 = 0 for each small value of the energy where averagingnormalisation is valid. They represent, like the normal modes, continuous families of periodic solutions parametrised by the energy. Adding the time-dependent terms destroy the periodicity but these terms become negligible after a long time.
In Figures 4 and 6 we start near the epicyclic normal mode, in Figures 5 and 8 near the vertical normal mode; in both cases the solutions move into general position on tori around periodic solutions. The asymmetric past of the system has changed the overall dynamics drastically, but this depends strongly on the initial phase where the potential is still asymmetric.
If n ≥ 2 we will have that c(t) changes O(ε) on this interval, in our example c(t) = −16 5 4 + O(ε). The adiabatic invariant (33) will be active on an interval of size 1/ε (n−1) . As we observe in Figure 6, (enlarged in Figure 9) the initial phase takes around 500 timesteps if n = 3, ε = 0.1. In the initial phase we have for the dynamics: (43) Figure 9. Enlargement of E 1 (t) of Figure 8 (the 1:1 resonance for the case n = 3 near the epicyclic normal mode). The initial phase where the asymmetric potential is active takes around 500 timesteps.

Evolution to the Hénon-Heiles Potential
The Hénon-Heiles potential is slightly degenerate as we have to normalise to terms of degree 6 (H 6 ) to obtain structurally stable results, see [14]. There are many details available on this system. If a 3 = a 4 = 0 and ε = 1 the system is chaotic but large-scale chaos sets in at E 0 ≥ 1/12 = 0.083, the energy manifold is bounded if 0 ≤ E 0 ≤ 1/6. In Figure 10 (left) we present for ε = 1 action E 1 (t) close to the chaotic stage. The first 600 timesteps shows irregular motion, then the dynamics settles in motion on tori. The normal forms obtained before will not be valid for these parameter values. Figure 10. Left E 1 (t) of the Hénon-Heiles system for the case n = 2, a 3 = 0, a 4 = 4 with ε = 1 close to the chaotic stage, x(0) = 0.4, z(0) = 0.1, E 1 (0) = 0.083, E 2 (0) = 0.005. Right the Hénon-Heiles system for the case n = 2, a 3 = 1, a 4 = 4, ε = 0.1 and for E 1 (0) = 0.005(x(0) = 0.1) 3 cases: ; we find small oscillations of E 1 (t) around a long-time oscillation (increase of E 2 (0) produces higher maxima and shorter return times.

1.
The resulting dynamics of a rotating axisymmetric galaxy to a state of mirror symmetry with respect to the galactic plane depends on at least three aspects: the local resonance ratio of epicyclic and vertical frequencies, the timescale of evolution to symmetric equilibrium characterised by ε n t, n = 1, 2, 3, . . . and the initial phase-space distribution.

2.
The dynamics of a potential producing locally the 1:2 resonance shows surprisingly that the main effects are described by an O(ε) approximation.

5.
An interesting conclusion is that in the final stage of evolution to mirror symmetry with respect to the galactic plane of a rotating axisymmetric galaxy their position in phase-space depends strongly on its asymmetric past. The epicyclic and vertical velocity dispersions will be different and are tied in with the resonance-dependent adiabatic invariants.
One of the implications is that a statistical mechanics description of systems of particles in near-equilibrium and held together by mainly gravitational forces can not be understood without considering earlier evolutionary states.

A Synopsis of Averaging Theorems
Averaging methods were used as a formal method in the 18th century. In the first half of the 20th century proofs and extensions of the methods were developed in France and in the Sowjet-Union. For a brief historical survey see appendix A of [11]. Sometimes the methods are referred to as averaging-normalisation as there is a definite relation with normal form methods, see [11] chs. 9-13.
Consider initial value problems of the form with z(t 0 ) given and ε a small parameter, ε ≥ 0. To analyse this initial value problem as a perturbation problem makes sense if we have knowledge about the system for ε = 0. Using this knowledge and applying variation of constants one can often obtain a system of slowly varying variables of the form where a T-periodic vector field plays a part. Consider explicitly the T-periodic vector fields f 1 , f 2 and the slowly varying ODE: With x, f 1 , f 2 ∈ R n and f 1 , f 2 twice continously differentiable in a bounded domain D in R n , continuously differentiable in t. Suppose in addition that: where x is kept constant during integration. We will call the terms represented by f 1 non-resonant to first order. We use the near-identity transformation: x(t) = y(t) + εu(t, y(t)), u(t, y(t)) = t 0 f 1 (s, y(t))ds.
I is the n × n unit matrix, the matrix I + ε∂u/∂y has a bounded inverse, so that: where the O(ε 2 ) terms can be computed explicitly. The procedure removes the non-resonant terms from the right-hand side of Equation (44); this is useful if we are able to perform analysis on the resonant part with explicit slow time as in Section 3. In our models a second small parameter δ plays a part. This is handled in various ways. We can assume a relation between ε and δ. In our models we have that δ occurs in the combination δt so we can introduce the independent slowly varying variable τ byτ = δ.
In formulating next a number of theorems we omit detailed assumptions regarding smoothness and dimensions. For a proof and more detailed conditions see [11] Chapter 4.2.
In many problems a second order approximation is necessary for qualitative and quantitative reasons. Introduce vector field f 1 (t, x) = D f (t, x)u 1 (t, x) − Du 1 (t, x) f 0 (x).
For a proof and more detailed conditions see [11] Chapter 4.5.
A useful extension was given in [12], see also [11] Chapter 2.9. It is: Theorem 3. Consider the formulation of theorem 2 and in addition assume f 0 (x) = 0. Then we have that |x(t) − v(t)| = O(ε) on time intervals of order 1/ε 2 .
The theory of higher order resonance for 2 dof Hamiltonian systems was started in [15]. It was extended in [4], see for a summary [11] Chapter 10.6.4. The theory is technically too involved to give a summary here. We mention two typical aspects. On studying 2 dof systems with basic k : l resonance, we have a higher order resonance if k + l ≥ 5; a lower order resonance is also described by the theory of higher order resonance if there are certain degenerations in the averaged system (this happens in our paper for the 1:2 and 1:3 resonance). In general the dynamics in phase-space is characterised by large domains without much interaction between the two degrees of freedom whereas interaction takes place in resonance manifolds with a fractal dimension in ε.
Funding: This research received no external funding.