We study the perturbedfromsynchronous librational state of a double asteroid modeled by the Full Two Rigid Body Problem F2RBP with primary emphasis on deriving analytical formulas which describe the systems evolution after deflection by a kinetic impactor. To this end both linear and nonlinear canonical theories are developed. We make the simplifying approximations to be relaxed in a forthcoming paper of planar binary orbit and axisymmetric shape of the primary body. To study the effect of a DARTlike hit on the secondary body the momentum transfer enhancement parameter inlineformula idIEq1mmlmathmmlmimmlmimmlmathinlineformula is introduced and retained as a symbolic variable throughout all formulas derived either by linear or nonlinear theory. Our approach can be of use in the context of the analysis of the postimpact data from kinetic impactor missions by providing a precise modeling of the impactors effect on the seconadrys librational state as a function of inlineformula idIEq2mmlmathmmlmimmlmimmlmathinlineformula.
['Binary asteroid', 'Kinetic impactor', 'Dynamical evolution', 'Canonical perturbation theory', 'Astrophysics - Earth and Planetary Astrophysics']
Analytical theory of the spinorbit state of a binary asteroid deflected by a kinetic impactor
{'Analytical theory of the spin-orbit state of a binary asteroid deflected by a kinetic impactor': "Michalis Gaitanas 1* , Christos Efthymiopoulos 2 , Ioannis Gkolias 1 , George Voyatzis 1 , Kleomenis Tsiganis 1  \n1  \n2 Trieste 63, Padova, 35121, Italy.  \nDepartment of Physics, Aristotle University of Thessaloniki, Thessaloniki, 54124, Greece. Department of Mathematics, Universit'a degli Studi di Padova, Via  \n*Corresponding author(s). E-mail(s): mgaitana@physics.auth.gr; Contributing authors: cefthym@math.unipd.it; igkoli@physics.auth.gr; voyatzis@auth.gr; tsiganis@auth.gr;", 'Abstract': "We study the perturbed-from-synchronous librational state of a double asteroid, modeled by the Full Two Rigid Body Problem (F2RBP), with primary emphasis on deriving analytical formulas which describe the system's evolution after deflection by a kinetic impactor. To this end, both a linear and nonlinear (canonical) theory are developed. We make the simplifying approximations (to be relaxed in a forthcoming paper) of planar binary orbit and axisymmetric shape of the primary body. To study the effect of a DART-like hit on the secondary body, the momentum transfer enhancement parameter β is introduced and retained as a symbolic variable throughout all formulas derived, either by linear or nonlinear theory. Our approach can be of use in the context of the analysis of the post impact data from kinetic impactor missions, by providing a precise modeling of the impactor's effect on the seconadry's librational state as a function of β .  \nKeywords: binary asteroid, kinetic impactor, dynamical evolution, canonical perturbation theory", '1 Introduction': "The historic DART (Double Asteroid Redirection Test) experiment was completed in September 2022, when the DART spacecraft impacted asteroid Dimorphos, the small secondary of (65803) Didymos, with the purpose of changing its 12-hour orbital period. Before impact, the system was assumed to be in a 'relaxed state', i.e. a singlysynchronous configuration of the secondary on a circular orbit, around the swiftly rotating primary (see Richardson et al. (2022)). As expected, the DART impact induced a significant orbit period change, measured to -33 ± 1 minutes (Thomas et al., 2023). Simulations of the post-impact dynamics suggest that the momentum enhancement parameter β 1 was even higher than expected, of order β = 3 . 6 (Cheng et al., 2023). The imparted eccentricity on the new orbit was estimated to be of order 0.03 - later, Meyer et al. (2023) showed that the eccentricity may be slowly decreasing, a possible indication of shape deformation (Raducan et al., 2024) and/or chaotic libration (Agrusa et al., 2021). It is interesting to note that the observed pre-impact shape, as derived in Daly et al. (2023) was much more oblate than originally thought, a factor that may play a critical role in the post-impact evolution of such systems.  \nIn this paper, motivated by the DART experiment, we delve into the dynamics of a double asteroid system, perturbed from the exact single-synchronous state due to collision with a DART-like kinetic impactor. Our study is in the framework of the Full Two Rigid Body Problem (F2RBP) (Maciejewski, 1995). We model analytically the post-impact librational state by two types of perturbation theories, hereafter referred to as i) linear, and ii) nonlinear canonical. In the former, the normal modes of a linearized system of equations of motion around a suitably defined equilibrium (representing the new synchronous state after the impact) are computed. In the latter, we use the Lie series method to compute series terms of order higher than linear in the small parameters of the problem.  \nTo set clear the main steps in the formulation of our theories, we presently make two simplifying assumptions, to be relaxed in a forthcoming paper. The first is planar motion, restricting the two rigid bodies to move in the x -y plane and to rotate only around their z -axis. This simplification reduces the problem's number of degrees of freedom, making the theory more manageable in terms of algebraic manipulations. The methods used, however, are straightforward to generalise in the full 3D problem, while the results found with the planar approach are already of practical utility in the applications. Secondly, we approximate the primary body as axisymmetric. As discussed in Section 2, such an approximation is justified whenever the primary rotates rapidly with respect to the orbital motion (and rotation of the synchronous secondary), as for example in the case of the single synchronous equilibrium state assumed to have been holding in the case of 65803 Didymos before DART's impact. Under the above assumptions, we then show how to derive analytical functions of time, which describe the binary's evolution, retaining as a free parameter in all formulas the momentum transfer enhancement β , which plays a crucial role in the timeseries of the observable quantities to be recovered after the mission.  \nOver the years, noteworthy advancements have been made in both the analytical and numerical study of the dynamics of binary asteroids. Notable contributions include (Kinoshita, 1972) who studied a binary system composed of a spherical and a triaxial body. Bou'e, G. and Laskar, J. (2009) studied the secular problem, using potential expansion up to order four and performed averaging over the fast angles. Scheeres (2009) explored the stability of the F2RBP under the planar assumption of two triaxial bodies. McMahon, J. W. and Scheeres, D. J. (2013) incorporated the axisymmetry of the primary body in the planar case and studied the bounds of librational motion of the secondary body. All these analytical results are obtained by the hypothesis of regular orbits. In the case of chaotic orbits, one has to rely on numerical methods. The numerical integrator GUBAS (Hou et al., 2017; Davis, A. B. and Scheeres, D. J. , 2020; Agrusa et al., 2021) has been widely employed over the years for the purpose of numerically investigating the F2RBP. In one of these studies, Meyer et al. (2021) showed that potentially observable, post-impact orbit period variations, whose characteristics depend on the physical parameters of the system, are induced to the mutual orbit.  \nWhile most previous analytical studies of the librational regime of binary asteroids are based on some form of linear theory around the synchronous state, in the present paper we discuss a higher precision, nonlinear analytical theory, capable of addressing even complex waveforms of the times series of the various observable quantities corresponding to the librational state after a DART-like hit. A key element of our approach is the employment of a 'book-keeping' method, which, at each order, collects in one (book-keeping) symbol all the various small quantities of the problem under study, grouped at similar orders of 'smallness' (see Section 4). Such small quantities are the eccentricity of the post-impact orbit, the amplitude of libration of the secondary, the effective J 2 term of the primary, the asphericity of the secondary, etc. We then show how, through the use of a unique book-keeping parameter, one can arrive easily at a suitable Birkhoff normalization process, by which a 'normal form' i.e., a solvable Hamiltonian is constructed.  \nIn comparison, we examine also classical linear theory based on the representation of all motions as a superposition of two linear normal modes. Of these modes, one represents a torque-induced libration combined with a nearly circular orbit, while the other an eccentricity-induced libration combined with the corresponding eccentric orbit. Our nonlinear theory then describes the coupling between these two fundamental modes. Overall, examining both nonlinear and linear approaches gives us the chance to compare the accuracy of analytical theories with respect to numerical simulations, as obtained, e.g. by GUBAS.  \nThe structure of the paper is as follows: Section 2 defines the problem in its general form and presents the two reductions mentioned earlier: planar motion and axisymmetric form of the primary body. We start with the definition of the initial conditions associated with the single synchronous equilibrium state, which serve as the point of departure for our study. Finally, we give some basic formulas, as well as a basic description of the effect on the system by the kinetic impactor, illustrated with numerical examples. Section 3 describes our linear perturbation theory. The problem's normal modes are derived in their generic form. Then we discuss how to reconfigure  \nthe final solution as a function of the β parameter. Section 4 introduces our canonical (normal form) approach. We expand the Hamiltonian into series in the above mentioned book-keeping parameter up to a desired order. Through the application of Birkhoff normalization by Lie series, we then construct a Hamiltonian normal form by which the analytical solution to the problem is obtained. Again, the numerical results are compared with those of the canonical theory. Section 5 gives a summary of our conclusions.", '2.1 Hamiltonian': "Let B 1 and B 2 be two rigid bodies of arbitrary mass distributions interacting gravitationally, embedded in a global inertial frame Oxyz as in Figure 1.  \nFig. 1 Definition of the F2RBP in the generic 3-dimensional form.  \n<!-- image -->  \nThe generic form of the mutual potential function is  \nV = -∫ B 1 ∫ B 2 GdM 1 dM 2 d , (1)  \nwhere G is the gravitational constant. The integration takes place over both mass distributions, and d is the distance between the finite mass elements dM 1 of B 1 and dM 2 of B 2 . Equation (1) can be expanded in the moments of inertia and up to the second order it reads (Scheeres, 2009)  \nV ( ⃗r, A 1 , A 2 ) = -GM 1 M 2 r -G 2 r 3 [ M 2 tr( I 1 ) + M 1 tr( I 2 ) ] + 3 G 2 r 5 ⃗r · ( M 2 A 1 I 1 A T 1 + M 1 A 2 I 2 A T 2 ) · ⃗r , (2)  \nwhere M 1 , M 2 are the primary's and secondary's masses, ⃗r the relative position vector between the two centers of mass, A 1 , A 2 the rotation matrices of B 1 , B 2 and I 1 , I 2 the (constant) inertia matrices expressed in each body frame, centered around the corresponding center of mass. The system admits the energy and angular momentum integrals given by  \nE = T + V = 1 2 M 1 M 2 M 1 + M 2 ⃗υ · ⃗υ + 1 2 ⃗ ω 1 · I 1 · ⃗ ω 1 + 1 2 ⃗ ω 2 · I 2 · ⃗ ω 2 + V ( ⃗r, A 1 , A 2 ) , (3)  \n⃗ L = M 1 M 2 M 1 + M 2 ⃗r × ⃗υ + A 1 · ( I 1 · ⃗ ω 1 ) + A 2 · ( I 2 · ⃗ ω 2 ) , (4)  \nwhere ⃗υ is the relative velocity vector between the two centers of mass and ⃗ ω 1 , ⃗ω 2 are the body frame angular velocity vectors of B 1 , B 2 .  \nWe now introduce two approximations:  \n1) Reduction to planar motion : assuming the mutual orbit to be restricted to the x -y plane and the rotation of both bodies to be restricted around their axes normal to the plane, the mutual potential function becomes  \nV ( r, ϕ 1 , ϕ 2 ) = -GM 1 M 2 r + GM 2 ( I 1 x + I 1 y -2 I 1 z ) + GM 1 ( I 2 x + I 2 y -2 I 2 z ) 4 r 3 + 3 GM 2 ( I 1 x -I 1 y ) cos (2 ϕ 1 ) 4 r 3 + 3 GM 1 ( I 2 x -I 2 y ) cos (2 ϕ 2 ) 4 r 3 . (5)  \nThe corresponding Hamiltonian is  \nH = T + V = p 2 r 2 m + p 2 ϕ 1 2 I 1 z + p 2 ϕ 2 2 I 2 z + ( p θ -p ϕ 1 -p ϕ 2 ) 2 2 mr 2 + V ( r, ϕ 1 , ϕ 2 ) . (6)  \nFig. 2 Schematic representation of the planar motion constraint. The two bodies are restricted to move in the x -y plane and to rotate only around their z -axis.  \n<!-- image -->  \nHere, m = ( M 1 M 2 ) / ( M 1 + M 2 ) is the reduced mass of the system and ϕ 1 , ϕ 2 are the relative yaw angles with respect to the orbital displacement angle θ . We refer below also to the angle θ 1 z = θ + ϕ 1 (see Figure 2). Finally, p r , p θ , p ϕ 1 , p ϕ 2 are the canonical momenta, conjugate to the coordinates r, θ, ϕ 1 , ϕ 2 . We have p r = m ˙ r , p ϕ 1 = I 1 z ( ˙ θ + ˙ ϕ 1 ), p ϕ 2 = I 2 z ( ˙ θ + ˙ ϕ 2 ), while p θ = mr 2 ˙ θ + p ϕ 1 + p ϕ 2 which is the total angular momentum of the system.  \n2) Approximation of the primary body as axisymmetric : such an approximation is justified whenever the primary is a fast rotator, so that all its non-axisymmetric multipole moments lead to fast oscillations which can be eliminated from the equations of motion. Formally, this implies adopting as value for both I 1 x , I 1 y the mean of the real values, I s = ( I 1 x + I 1 y ) / 2. The latter assumption eliminates the Hamiltonian's dependence on the angle ϕ 1 (or θ 1 z ), thus giving rise to four equilibrium points, corresponding to the values ϕ 2 = 0 , π, ± π/ 2. The final Hamiltonian then takes the form  \nH = p 2 r 2 m + p 2 ϕ 1 2 I 1 z + p 2 ϕ 2 2 I 2 z + ( p θ -p ϕ 1 -p ϕ 2 ) 2 2 mr 2 + V ( r, ϕ 2 ) , (7)  \nwhere  \nV ( r, ϕ 2 ) = -GM 1 M 2 r + GM 1 ( I 2 x + I 2 y -2 I 2 z ) + 2 GM 2 ( I s -I 1 z ) 4 r 3 + 3 GM 1 ( I 2 x -I 2 y ) cos (2 ϕ 2 ) 4 r 3 . (8)", '2.2 Pre-impact synchronous equilibrium state': 'Given the above approximations, and assuming a (fast) rotation of the primary with angular velocity ν 1 = ˙ θ 1 z (0), the following set of initial conditions define the so-called single synchronous equilibrium state of the binary at a separation (distance between the centers of mass) equal to r eq :  \nFig. 3 Representation of the pre-impact planar equilibrium configuration.  \n<!-- image -->  \nr (0) = r eq ϕ 2 (0) = ϕ 2 eq = 0 (9)  \n˙ r (0) = ˙ r eq = 0 ˙ θ (0) = ˙ θ eq = √ G ( M 1 + M 2 ) r 3 eq [ 1 + 3 2 r 2 eq ( I 1 z -I s M 1 + I 2 y + I 2 z -2 I 2 x M 2 )] ˙ ϕ 1 (0) = ˙ ϕ 1 eq = ν 1 -˙ θ eq ˙ ϕ 2 (0) = ˙ ϕ 2 eq = 0 (10)  \np r (0) = p r eq = 0 p θ (0) = p θ eq = p ϕ 1 eq + p ϕ 2 eq + mr 2 eq ˙ θ eq p ϕ 1 (0) = p ϕ 1 eq = ν 1 I 1 z p ϕ 2 (0) = p ϕ 2 eq = ˙ θ eq I 2 z . (11)  \nWe note that the entire set of the above initial conditions are determined by two free parameters, namely r eq and ν 1 (see Figure 3).', '2.3 Post-impact evolution': "Assume, now, that the secondary undergoes a DART-like head-on collision, with the impactor's velocity vector lying always in the x -y plane, being opposite to the velocity vector of the secondary, and pointing towards the secondary's center (Figure 3). If M D , υ D are the mass and velocity magnitude of the impactor, then the velocity change ∆ υ induced to the secondary can be expressed in terms of the momentum enhancement factor β through the equation (Rivkin et al., 2021)  \n∆ υ = -βM D υ D M 2 . (12)  \nBecause the recoiled ejecta are expelled in the opposite direction of the impactor's velocity, the momentum change is enhanced, yielding β > 1. If no ejecta are expelled, the collision is considered perfectly inelastic, resulting in β = 1.  \n̸  \nDepending on the parameters of the impactor and of the binary, we distinguish two possible evolution cases after the impact: regular or chaotic. In a chaotic case, little can be done in terms of analytical theories as indicated in some cases found by Agrusa et al. (2021) and hence, one must resort to numerical integration. If the collision is not 'violent' enough to provoke chaotic evolution, the post impact state remains regular and can be described approximately by analytical formulas, which we focus on the present paper. Physically, after the impact, the trajectory turns from circular to eccentric, with the apocenter initially coinciding with the point of impact, while the secondary also develops a libration ϕ 2 ( t ) = 0.  \nFigures 4 to 7 show four examples of the numerical evolution of the functions r ( t ) and ϕ 2 ( t ) post impact (blue curves), obtained by assuming two different sets of physical parameters for the binary system, as well as two possible values of β for each set. Typically, we find that r ( t ) just follows a low amplitude epicyclic oscillation due to the eccentricity induced to the orbit. On the other hand, the angle ϕ 2 ( t ) starts exhibiting a librational behavior. As shown in the next section, this is dominated by two harmonics, corresponding to two fundamental frequencies, and several linear combinations thereof. In general, the amplitude of the libration increases with β and it can reach several degrees, leading to a substantially nonlinear behavior of the system. In the case of regular motion, nonlinearity manifests itself mostly in the fact that the fundamental frequencies of libration, as derived by a numerical analysis, differ from the corresponding linear (normal mode) frequencies. This fact causes a slow dephasing of the analytical solutions from the numerical ones, as eye-evident, for  \nexample, in Figure 7. The orange-dashed curves in each panel of Figures 4 to 7 show the corresponding analytical prediction for each orbit, obtained by the different theories discussed in detail in Sections 3 and 4 below. While deferring a detailed discussion to the corresponding Sections, we observe, already at this stage, that the reduction of the error of analytical theories at higher order is due, mostly, to the improvement in the predictions by which a theory recovers the numerical values of the fundamental frequencies as the order of the theory increases, preserving the parameter β at the same time. We now turn our attention to the analytical methods used to achieve such improvement.  \nFig. 4 Post impact evolution of the functions r ( t ) (left column) and ϕ 2 ( t ) (right column) for the case β = 1 and for the first set of physical parameters ( G = 0 . 0864989, M 1 = 5 . 15045, M 2 = 0 . 0392647, I s = 0 . 263844, I 1 z = 0 . 337921, I 2 x = 8 . 20357 · 10 -5 , I 2 y = 8 . 88678 · 10 -5 , I 2 z = 1 . 18976 · 10 -4 , r eq = 1 . 18, ν 1 = 2 π/ 2 . 26, M D = 5 . 79434 · 10 -9 , υ D = 22121 . 6). The units of measurement are [ hr ] , [ km ] , [ kg ∗ ], where 1[ kg ∗ ] = 10 11 [ kg ]. The blue curves correspond to the results derived via numerical integration, while the dashed orange curves correspond to the results obtained via each of the examined analytical theory i.e., as indicated on top of each panel. The green curves represent the initial (pre-impact) equilibrium state of the binary.  \n<!-- image -->  \nFig. 5 Same as in Figure 4, but for the case β = 3.  \n<!-- image -->  \nFig. 6 Same as in Figure 4, but for the case β = 1 and for the second set of physical parameters ( G = 0 . 0864989, M 1 = 5 . 14073, M 2 = 0 . 048272, I s = 0 . 254128, I 1 z = 0 . 326872, I 2 x = 9 . 28891 · 10 -5 , I 2 y = 1 . 27131 · 10 -4 , I 2 z = 1 . 29876 · 10 -4 ), r eq = 1 . 18, ν 1 = 2 π/ 2 . 26, M D = 5 . 79434 · 10 -9 , υ D = 22121 . 6). The units remain the same.  \n<!-- image -->  \nFig. 7 Same as in Figure 6, but for the case β = 3.  \n<!-- image -->", '3.1 New equilibrium': "As soon as the impact occurs, the binary is deflected away from its pre-impact synchronous equilibrium state. Formally, the momentum transfer by the impactor affects the value of the integral of motion p θ in the Hamiltonian (7), a head-on collision leading to a decrease of p θ . In terms of the parameter β , the post-impact state then starts from the initial conditions  \n˙ r imp (0) = 0 , p r imp (0) = p r eq = 0 ˙ θ imp (0) = ˙ θ eq -M D υ D β r eq M 2 , p θ imp (0) = p θ eq -mr eq M D υ D β M 2 ˙ ϕ 1 imp (0) = ν 1 -˙ θ imp (0) , p ϕ 1 imp (0) = p ϕ 1 eq ˙ ϕ 2 imp (0) = M D υ D β r eq M 2 , p ϕ 2 imp (0) = p ϕ 2 eq . (13)  \nSince the ratio M D /M 2 in (13) is a small quantity, the resulting new state can be studied as a small deviation from a new equilibrium state, corresponding to the new value of the integral p θ = p θ imp . We specify the (constant in time) distance r eq,new of the two bodies at the new equilibrium by solving (for r eq,new ) the equation  \n˙ θ eq,new = ( I 2 z + mr 2 eq ) ˙ θ eq I 2 z + mr 2 eq,new -mr eq M D υ D β M 2 ( I 2 z + mr 2 eq,new ) = √ G ( M 1 + M 2 ) r 3 eq,new [ 1 + 3 2 r 2 eq,new ( I 1 z -I s M 1 + I 2 y + I 2 z -2 I 2 x M 2 )] . (14)  \nNote that, in order to keep β as a symbol in all the expressions related to the new equilibrium state, one has to perform a Taylor expansion of equation (14) in the small quantity δr = r eq -r eq,new around the value δr = 0, which corresponds to no change in the kinetic state. We find  \nr eq,new = r eq + C r eq β + O ( δr 2 ) , (15)  \nwhere  \nC r eq = 2 mr 9 / 2 eq √ 6 C I +4 r 2 eq M D υ D M 2 √ G ( M 1 + M 2 ) ( 6 I 2 z r 2 eq -2 mr 4 eq +3 C I (5 I 2 z + mr 2 eq ) ) C I = I 1 z -I s M 1 + -2 I 2 x + I 2 y + I 2 z M 2 . (16)  \nNote that even if β > 1, the quantity C r eq β is O ( M D /M 2 ), hence small. Substituting r eq,new , as given by equation (15), whenever necessary in the formulas of the linear theory leads to a method referred to as 'Linear Taylor' (see Figures 4 to 7). Using a root-finding algorithm to numerically determine the exact solution of equation (14) is also possible, but only marginally alters the final results. Figure 8 depicts the absolute error between the new equilibrium point obtained via root-finding and the one via equation (15). For the expected values of β , the maximum error is within the order of a few decimetres.  \nFig. 8 Absolute error between the new equilibrium point, computed via root-finding and the one computed via the 'Linear Taylor' approximation.  \n<!-- image -->", '3.2 Linear analytical solution for the post-impact state': "Since the angles θ and ϕ 1 are ignorable in the Hamiltonian (7), the quantities p θ and p ϕ 1 are constants of motion, whose values are determined only by the initial conditions, i.e., by equations (9) - (11), with r eq replaced by r eq,new , found in turn through equation (15) or numerically. We then consider the 4x4 reduced system of Hamilton's equations in the variables ( r, ϕ 2 , p r , p ϕ 2 )  \n˙ r = ∂H ∂p r ˙ p r = -∂H ∂r ˙ ϕ 2 = ∂H ∂p ϕ 2 ˙ p ϕ 2 = -∂H ∂ϕ 2 , (17)  \nwhere H = H ( r, ϕ 2 , p r , p ϕ 2 ; p θ = p θ imp , p ϕ 1 = p ϕ 1 eq ). Setting r = r eq,new + δr , ϕ 2 = 0 + δϕ 2 , p r = 0+ δp r , p ϕ 2 = p ϕ 2 eq,new + δp ϕ 2 , with p ϕ 2 eq,new = ˙ θ eq,new ( r eq,new ) I 2 z , and linearizing the 4x4 equations of motion around the new equilibrium point, leads to the linear system  \n˙ -→ δX = J · -→ δX , (18)  \nwhere  \n-→ δX = [ δr, δϕ 2 , δp r , δp ϕ 2 ] T and J =     0 0 j 13 0 j 21 0 0 j 24 j 31 0 0 -j 21 0 j 42 0 0     . (19)  \nThe Jacobian's non-zero elements j nm (evaluated at the new equilibrium point) depend only on the adopted values of the physical parameters of the binary and on the value of  \nβ (through their dependence on r eq,new in equation (15). The characteristic polynomial is  \nP ( λ ) = λ 4 + ξ 2 λ 2 + ξ 0 , (20)  \nwhere the analytical expressions for j nm , ξ 0 and ξ 2 are given in the Appendices A and B respectively. The roots of the characteristic polynomial are λ 1 , 2 = ∓ iω 1 , λ 3 , 4 = ∓ iω 2 , where the two fundamental frequencies  \nω 1 , 2 = √ ξ 2 ± √ ξ 2 2 -4 ξ 0 √ 2 , (21)  \ncorrespond, respectively, to the normal modes associated with the eccentricity-induced and with the torque-induced librational motion of the secondary. General expressions for the solution of the system (18) are easy to compute for an arbitrary initial excitation ( δr (0) , δϕ 2 (0) , δp r (0) , δp ϕ 2 (0) ) , by superposing the solutions corresponding to the two normal modes. However, the resulting expressions are cumbersome and will not be presented here. We focus, instead, on the simple case of initial conditions corresponding to a DART-like hit, for which δr (0) = r eq -r eq,new , δϕ 2 (0) = 0, δp r (0) = 0, δp ϕ 2 (0) = p ϕ 2 eq -p ϕ 2 eq,new (Appendix D). We then obtain the expressions  \nδp  \nδr ( t ) = A 1 r cos ( ω 1 t ) + A 2 r cos ( ω 2 t ) δϕ 2 ( t ) = ( j 13 ( j 2 21 + j 24 j 31 ) + j 24 ω 2 1 ) A 1 r j 13 j 21 ω 1 sin ( ω 1 t ) + ( j 13 ( j 2 21 + j 24 j 31 ) + j 24 ω 2 2 ) A 2 r j 13 j 21 ω 2 sin ( ω 2 t ) (22) δp r ( t ) = -ω 1 A 1 r j 13 sin ( ω 1 t ) -ω 2 A 2 r j 13 sin ( ω 2 t ) ϕ 2 ( t ) = ( j 13 j 31 + ω 2 1 ) A 1 r j 13 j 21 cos ( ω 1 t ) + ( j 13 j 31 + ω 2 2 ) A 2 r j 13 j 21 cos ( ω 2 t ) ,  \nwhere the coefficients A 1 r , A 2 r are functions of j nm , ω 1 , ω 2 , δr (0) and δp ϕ 2 (0), and are given in the Appendix C. Due to the explicit dependence of r eq,new on β (see equation (15)), it turns out that the linear amplitudes, frequencies and hence, the whole linear solution can be written using explicit formulas as functions of β . We then obtain an analytical (called 'Linear-Taylor') approximation ( r lin ( t ; β ) , ϕ 2 lin ( t ; β ) , p r lin ( t ; β ) , p ϕ 2 lin ( t ; β ) ) of the solution of Hamilton's equations for a DART-like hit given by  \nr lin ( t ; β ) = r eq,new ( β ) + δr ( t ; β ) ϕ 2 lin ( t ; β ) = δϕ 2 ( t ; β ) p r lin ( t ; β ) = δp r ( t ; β ) p ϕ 2 lin ( t ; β ) = p ϕ 2 eq,new ( r eq,new ( β )) + δp ϕ 2 ( t ; β ) . (23)  \nRegarding the relative orbital displacement θ , and the primary's relative rotation ϕ 1 , explicit functions of time and β can be derived for both, assuming linear approximation to the right hand side of ˙ θ (detailed derivation is presented in Appendix E). We find  \nθ ( t ; β ) = θ (0) + A 0 θ t + A 1 θ sin ( ω 1 t ) + A 2 θ sin ( ω 2 t ) + A 3 θ sin (2 ω 1 t ) + A 4 θ sin (2 ω 2 t ) + A 5 θ sin [ ( ω 1 -ω 2 ) t ] + A 6 θ sin [ ( ω 1 + ω 2 ) t ] , (24)  \nwhere the coefficients A jθ (given in Appendix E) depend on the physical parameters and on β . The analytical formula of θ ( t ; β ) can be used to directly evaluate the orbital period and hence, its variation over time T ( t ). The latter is in good agreement with the numerical computations presented in (Meyer et al., 2021). The angle ϕ 1 is given by  \nϕ 1 ( t ; β ) = ϕ 1 (0) + p ϕ 1 I 1 z t -θ ( t ; β ) . (25)  \nThe left panels in Figures 4 to 7 compare the solution (23) with the numerical one for two different values of β and two parameter sets, as indicated in each Figure's caption. One can deduce that for β = 1, the solutions stay close to each other for several orbital periods (even months), but for β = 3, the solutions start deviating after a few periods. In fact, a gradual dephasing of the analytical from the numerical curves is present in all four cases, but eye-evident only in Figures 6 and 7. Consequently, while the essential features of the motion and its basic parameters (amplitudes and frequencies) are determined relatively accurately by our linear model, they cannot be used always to express the long-term behavior of the perturbed system. For this purpose, a higher-order analytical theory is needed. In the following Section we present a more accurate model of the perturbed motion, using canonical perturbation theory.", '4.1 Hamiltonian preparation': "In the previous Section, we discussed linear approximations to the post-deflection state computed via the equations of motion linearized around the new state ( r eq,new , p θ imp ). Nonlinear corrections to the above thereby can be obtained by considering terms of order higher than linear in the equations of motion, or higher than quadratic in the expansion of the Hamiltonian in powers of the quantities ( δr, δϕ 2 , δp r , δp ϕ 2 ). In the present Section, we examine how to produce such nonlinear theories using the Lie series method of canonical perturbation theory. However, for algorithmic reasons, and also due to its necessity in the normal form theory developed for the full 3D triaxial model in a forthcoming paper, it turns convenient to first slightly modify the point of departure of such theories, i.e. the initial point around which the series expansions are to be made. To this end, consider a decomposition of the Hamiltonian (7) as  \nH = H 0 + H 1 , (26)  \nwhere  \nH 0 = p 2 r 2 m + p 2 ϕ 1 2 I 1 z + p 2 ϕ 2 2 I 2 z + ( p θ -p ϕ 1 -p ϕ 2 ) 2 2 mr 2 -GM 1 M 2 r (27)  \nH 1 = GM 1 ( I 2 x + I 2 y -2 I 2 z ) + 2 GM 2 ( I s -I 1 z ) 4 r 3 + 3 GM 1 ( I 2 x -I 2 y ) cos (2 ϕ 2 ) 4 r 3 . (28)  \nThe Hamiltonian H 0 describes the motion of two perfect spheres of masses M 1 , M 2 and moments of inertia I 1 z , I 2 z . The corresponding synchronous (Keplerian) equilibrium at separation r = r ∗ of the binary is given by:  \nr ∗ = r (0) ϕ ∗ 2 = 0 p ∗ r = 0 p ∗ θ = p ∗ ϕ 1 + p ∗ ϕ 2 + mr ∗ 2 ν ∗ θ p ∗ ϕ 1 = ν 1 I 1 z p ∗ ϕ 2 = ν ∗ θ I 2 z , (29)  \nwith ν 1 = ˙ θ 1 z (0) and ν ∗ θ = √ G ( M 1 + M 2 ) /r ∗ 3 . Consider now an expansion of the Hamiltonian (27), up to a chosen truncation order N + 2, N ≥ 1 in the variables ( δr, δϕ 2 , δp r , δp ϕ 2 ) defined by the (canonical) transformation  \nr = r ∗ + ϵ · δr ϕ 2 = ϵ · δϕ 2 p r = ϵ · δp r p θ = p ∗ θ + ϵ 2 · δp θ p ϕ 1 = p ∗ ϕ 1 + ϵ 2 · δp ϕ 1 p ϕ 2 = p ∗ ϕ 2 + ϵ · δp ϕ 2 , (30)  \nwhere ϵ = 1 is a 'book-keeping' symbol. Following the 'book-keeping' method (see Efthymiopoulos (2012)), in all series expansions defined below we keep ϵ as a symbol, whose various powers collect terms representing a similar order of smallness of the problem. Since ( r ∗ , ϕ ∗ 2 = 0 , p ∗ r = 0 , p ∗ ϕ 2 ; p ∗ θ , p ∗ ϕ 1 ) is an equilibrium of the 4x4 reduced set of Hamilton's equations, introducing (30) into the Hamiltonian H 0 and expanding in powers of ϵ , leads to:  \nH 0 = GI 2 z ( M 1 + M 2 ) -GM 1 M 2 r ∗ 2 + ν 2 1 I 1 z r ∗ 3 2 r ∗ 3 + ϵ 2 ( Π 2 ( δr, δϕ 2 , δp r , δp ϕ 2 ) + ν ∗ θ δp θ +( ν 1 -ν ∗ θ ) δp ϕ 1 + O ( ϵ ) ) , (31)  \nwhere Π 2 is a polynomial of second degree in its arguments. We find  \nΠ 2 = GM 1 M 2 2 r ∗ 3 δr 2 + 1 2 m δp 2 r + ( 1 2 I 2 z + 1 2 mr ∗ 2 ) δp 2 ϕ 2 + 2 ν ∗ θ r ∗ δrδp ϕ 2 . (32)  \nSince ϵ = 1, we can lower at no cost all the book-keeping powers in the above expression by two, thus re-expressing H 0 as  \nH 0 = GI 2 z ( M 1 + M 2 ) -GM 1 M 2 r ∗ 2 + ν 2 1 I 1 z r ∗ 3 2 r ∗ 3 + GM 1 M 2 2 r ∗ 3 δr 2 + 1 2 m δp 2 r + ( 1 2 I 2 z + 1 2 mr ∗ 2 ) δp 2 ϕ 2 + 2 ν ∗ θ r ∗ δrδp ϕ 2 + ν ∗ θ δp θ +( ν 1 -ν ∗ θ ) δp ϕ 1 + ∞ ∑ j =1 ϵ j H 0 ,j +2 ( δr, δϕ 2 , δp r , δp ϕ 2 ; δp θ , δp ϕ 1 ) . (33)  \nSince δp θ , δp ϕ 1 are integrals of motion, the constant terms ν ∗ θ δp θ +( ν 1 -ν ∗ θ ) δp ϕ 1 can be omitted in the study of the reduced 4x4 dynamics in the variables ( δr, δϕ 2 , δp r , δp ϕ 2 ). Keeping them emphasises, however, the fundamental difference between the degrees of freedom ( θ, δp θ ) and ( ϕ 1 , δp ϕ 1 ), on one hand, and ( δr, δp r ), ( δϕ 2 , δp ϕ 2 ) on the other. Namely, the system is a rotator in the former, with frequencies ν θ , ν 1 -ν ∗ θ , while it is an oscillator in the latter, with frequencies to be recovered by a diagonalization + normalization procedure as described in the sequel. Also, keeping the term ( ν 1 -ν ∗ θ ) δp ϕ 1 becomes essential once we relax the condition of axisymmetry of the primary, since then, the Hamiltonian becomes dependent on the angle ϕ 1 , and the above term properly enters the kernel of the normalization of the Hamiltonian (the term δp θ still remains integral of motion in the full triaxial 3D case).  \nReturning, to the Hamiltonian term H 1 , substituting (30) and expanding in powers of ϵ , leads to:  \nH 1 = -G ( M 2 ( I 1 z -I s ) + M 1 ( -2 I 2 x + I 2 y + I 2 z ) ) 2 r ∗ 3 + ϵ 3 G ( M 2 ( I 1 z -I s ) + M 1 ( -2 I 2 x + I 2 y + I 2 z ) ) 2 r ∗ 4 δr + ϵ 2 [ -3 G ( M 2 ( I 1 z -I s ) + M 1 ( -2 I 2 x + I 2 y + I 2 z ) ) r ∗ 5 δr 2 + 3 GM 1 ( I 2 y -I 2 x ) 2 r ∗ 3 δϕ 2 2 ] + ∞ ∑ j =3 ϵ j H 1 ,j ( δr, δϕ 2 , δp r , δp ϕ 2 ; δp θ , δp ϕ 1 ) . (34)  \nSince ϵ = 1, we lower again by two the book-keeping powers in the terms of H 1 ,j , but this time only for j ≥ 2. This leads to the book-kept Hamiltonian H 1  \nH 1 = -G ( M 2 ( I 1 z -I s ) + M 1 ( -2 I 2 x + I 2 y + I 2 z ) ) 2 r ∗ 3 -3 G ( M 2 ( I 1 z -I s ) + M 1 ( -2 I 2 x + I 2 y + I 2 z ) ) r ∗ 5 δr 2 + 3 GM 1 ( I 2 y -I 2 x ) 2 r ∗ 3 δϕ 2 2 + ϵ 3 G ( M 2 ( I 1 z -I s ) + M 1 ( -2 I 2 x + I 2 y + I 2 z ) ) 2 r ∗ 4 δr + ∞ ∑ j =1 ϵ j H 1 ,j +2 ( δr, δϕ 2 , δp r , δp ϕ 2 ; δp θ , δp ϕ 1 ) . (35)  \nThe final book-kept Hamiltonian has the form  \nH = H 0 + H 1 = Z 0 + ∞ ∑ j =1 ϵ j h j , (36)  \nwhere Z 0 and h j are obtained by adding the corresponding terms of H 0 and H 1 of the same book-keeping powers of ϵ . Explicit expressions of the terms Z 0 , h 1 and h 2 are given in the Appendix F.  \nThe next step is to change our set of variables, from generalized positions and momenta ( δr, δϕ 2 , δp r , δp ϕ 2 ) to Birkhoff complex canonical variables ( Q 1 , Q 2 , P 1 , P 2 ), derived by diagonalizing the quadratic part of the Hamiltonian (36), i.e. the term Z 0 . The relationship between ( δr, δϕ 2 , δp r , δp ϕ 2 ) and ( Q 1 , Q 2 , P 1 , P 2 ) is a linear canonical transformation  \n    δr δϕ 2 δp r δp ϕ 2     =     m 11 m 12 m 13 m 14 m 21 m 22 m 23 m 24 m 31 m 32 m 33 m 34 m 41 m 42 m 43 m 44         Q 1 Q 2 P 1 P 2     , (37)  \nwhere the elements m ij depend on the physical parameters of the system, as well as on the constant r ∗ . Detailed derivation of m ij is given in the Appendix G. Substitution of (37) into (36) leads to  \nH ( Q 1 , Q 2 , P 1 , P 2 ) = Z 0 + N ∑ j =1 ϵ j h j ( Q 1 , Q 2 , P 1 , P 2 ) , (38)  \nwhere  \nZ 0 = C Z 0 + ν ∗ θ δp θ +( ν 1 -ν ∗ θ ) δp ϕ 1 -iω 1 k Q 1 P 1 -iω 2 k Q 2 P 2 , (39)  \nwith  \nC Z 0 = 1 2 ( ν 2 1 I 1 z + G ( M 1 (2 I 2 x -I 2 y -M 2 r ∗ 2 ) + M 2 ( -I 1 z + I 2 z + I s ) ) r ∗ 3 ) , (40)  \nand ω 1 k , ω 2 k being the kernel frequencies (see expressions in the Appendix G).", '4.2 Hamiltonian normal form and solution': "In order to extract analytical solutions from the book-kept Hamiltonian (38), we apply a Birkhoff canonical normalization process using the method of Lie series (Hori, 1966; Deprit, 1969) up to N normalization steps. Detailed description of the process, in accordance with Efthymiopoulos (2012) is presented in Appendix H. After performing N steps of the normalization algorithm upon the Hamiltonian (38), we obtain the normal form Hamiltonian  \nZ ( N ) = N 1 ≤ N ∑ j 1 =0 N 2 ≤ N ∑ j 2 =0 N 3 ≤ N ∑ j 3 =0 N 4 ≤ N ∑ j 4 =0 C ( N ) j 1 j 2 j 3 j 4 ( Q ( N ) 1 ) j 1 ( P ( N ) 1 ) j 2 ( Q ( N ) 2 ) j 3 ( P ( N ) 2 ) j 4 . (41)  \nIn equation (41), the constant coefficients C ( N ) j 1 j 2 j 3 j 4 encapsulate the integrals of motion δp θ , δp ϕ 1 , i.e. C ( N ) j 1 j 2 j 3 j 4 = C ( N ) j 1 j 2 j 3 j 4 ( δp θ , δp ϕ 1 ), as well as the physical parameters. The variables Q ( N ) j and P ( N ) j appear only in the form of the products Q ( N ) j P ( N ) j , j = 1 , 2. The real quantities I j = iQ ( N ) j P ( N ) j are then integrals of motion under the flow of Z ( N ) , since  \n{ I j , Z ( N ) } = { iQ ( N ) j P ( N ) j , Z ( N ) } = 0 , j = 1 , 2 . (42)  \nDue to the canonical nature of transformations, the preservation of the form of Hamilton's equations is achieved, thereby yielding  \n˙ Q ( N ) 1 = ∂Z ( N ) ∂P ( N ) 1 = iω 1 c Q ( N ) 1 ˙ Q ( N ) 2 = ∂Z ( N ) ∂P ( N ) 2 = iω 2 c Q ( N ) 2 ˙ P ( N ) 1 = -∂Z ( N ) ∂Q ( N ) 1 = -iω 1 c P ( N ) 1 ˙ P ( N ) 2 = -∂Z ( N ) ∂Q ( N ) 2 = -iω 2 c P ( N ) 2 ˙ θ ( N ) = ∂Z ( N ) ∂δp θ = ω θc ˙ ϕ ( N ) 1 = ∂Z ( N ) ∂δp ϕ 1 = ω ϕ 1 c , (43)  \nwhere the quantities ω 1 c , ω 2 c , ω θc and ω ϕ 1 c contain only products of the form Q ( N ) j P ( N ) j , i.e. ω 1 c = ω 1 c ( I 1 , I 2 ), ω 2 c = ω 2 c ( I 1 , I 2 ), ω θc = ω θc ( I 1 , I 2 ) and ω ϕ 1 c = ω ϕ 1 c ( I 1 , I 2 ). It follows that equations (43) can be directly integrated, giving the solutions  \nQ ( N ) 1 ( t ) = Q ( N ) 1 (0) e iω 1 c t Q ( N ) 2 ( t ) = Q ( N ) 2 (0) e iω 2 c t P ( N ) 1 ( t ) = P ( N ) 1 (0) e -iω 1 c t P ( N ) 2 ( t ) = P ( N ) 2 (0) e -iω 2 c t θ ( N ) ( t ) = θ ( N ) (0) + ω θc t ϕ ( N ) 1 ( t ) = ϕ ( N ) 1 (0) + ω ϕ 1 c t , (44)  \nwhere ( Q ( N ) 1 (0) , Q ( N ) 2 (0) , P ( N ) 1 (0) , P ( N ) 2 (0) ) are the initial conditions of any trajectory transformed to the Birkhoff variables, and subjected to N transformations with the Lie series. In Hamilton's equations, we include the variables ( θ ( N ) , ϕ ( N ) 1 ) , which, on the one hand are subjected to the Lie series transformations, but on the other hand do not affect the Hamiltonian normalization process since they are ignorable coordinates. Mapping the solution (44) into the original variables ( δr ( t ) , δϕ 2 ( t ) , δp r ( t ) , δp ϕ 2 ( t ) ) and ( θ ( t ) , ϕ 1 ( t ) ) , yields the final solution  \nω  \nr can ( t ; β ) = r ∗ + P r, 0 ( β ) + r,tot ∑ j =1 C r,j β k r,j cos [ ( n j ω 1 c + m j ω 2 c ) t ] ϕ 2 can ( t ; β ) = ω ϕ 2 ,tot ∑ j =1 C ϕ 2 ,j β k ϕ 2 ,j sin [ ( n j ω 1 c + m j ω 2 c ) t ] p r can ( t ; β ) = ω pr,tot ∑ j =1 C p r ,j β k pr,j sin [ ( n j ω 1 c + m j ω 2 c ) t ] p ϕ 2 can ( t ; β ) = p ∗ ϕ 2 + P p ϕ 2 , 0 ( β ) + ω p ϕ 2 ,tot ∑ j =1 C p ϕ 2 ,j β k pϕ 2 ,j cos [ ( n j ω 1 c + m j ω 2 c ) t ] θ can ( t ; β ) = θ (0) + ω θc ( β ) t + ω θ,tot ∑ j =1 C θ,j β k θ,j sin [ ( n j ω 1 c + m j ω 2 c ) t ] ϕ 1 can ( t ; β ) = ϕ 1 (0) + p ϕ 1 I 1 z t -θ can ( t ; β ) . (45)  \nIn equation (45), the frequencies ω 1 c , ω 2 c , ω θc and ω ϕ 1 c are polynomial functions of β  \nω 1 c ( β ) = Ω 10 +Ω 11 β + ... +Ω 1 N β N + ... ω 2 c ( β ) = Ω 20 +Ω 21 β + ... +Ω 2 N β N + ... ω θc ( β ) = Ω θ 0 +Ω θ 1 β + ... +Ω θN β N + ... ω ϕ 1 c ( β ) = Ω ϕ 1 0 +Ω ϕ 1 1 β + .. +Ω ϕ 1 N β N + ... , (46)  \nwith Ω 1 j , Ω 2 j , Ω θj , Ω ϕ 1 j depending on the physical parameters. Also, the polynomials P r, 0 ( β ) and P p ϕ 2 , 0 ( β ) have the form  \nP r, 0 ( β ) = D r, 0 + D r, 1 β + ... + D r,N β N + ... P p ϕ 2 , 0 ( β ) = D p ϕ 2 , 0 + D p ϕ 2 , 1 β + ... + D p ϕ 2 ,N β N + ... (47)  \nwhere D r,j , D p ϕ 2 ,j depend on the physical parameters as well. The numbers n j , m j are integers that constitute combinations of ω 1 c and ω 2 c . The counters ω r,tot , ω ϕ 2 ,tot , ω p r ,tot , ω p ϕ 2 ,tot and ω θ,tot are the total number of combinations of ( n j , m j ), correspondingly for the functions r can ( t ; β ), ϕ 2 can ( t ; β ), p r can ( t ; β ), p ϕ 2 can ( t ; β ) and θ can ( t ; β ). The amplitude coefficients C r,j , C ϕ 2 ,j , C p r ,j , C p ϕ 2 ,j and C θ,j depend on the physical parameters. Comparing equations (23) and (45), we can see that the canonical perturbation theory adds corrections to the functions r ( t ; β ), ϕ 2 ( t ; β ), p r ( t ; β ) and p ϕ 2 ( t ; β ) through higher order polynomials of β . Those corrections affect the amplitudes ( C r,j β k r,j , C ϕ 2 ,j β k ϕ 2 ,j , C p r ,j β k pr,j , C p ϕ 2 ,j β k pϕ 2 ,j ), the frequencies ( n j ω 1 c ( β ) + m j ω 2 c ( β )) and the equilibrium point ( r ∗ + P r, 0 ( β ), p ∗ ϕ 2 + P p ϕ 2 , 0 ( β )) for which our initial guess is the Keplerian one. What's more, the canonical theory corrects the orbital angle θ (by which the orbital period can be measured), by changing the Keplerian orbital rate from ν θ to ω θc ( β ) and by adding appropriate trigonometric variations.", '4.3 Numerical comparisons and precision tests': "So far we have examined two types of theories by which we can obtain a semi-analytic representation of the solutions δr ( t ) , δϕ 2 ( t ) , δp r ( t ) , δp ϕ 2 ( t ) , θ ( t ) , ϕ 1 ( t ) for perturbed librational motions around the isochronous equilibrium state: i) the 'Linear Taylor' theory of Section 3, and ii) the nonlinear theories of the present section. These theories can be diversified in two aspects. First, the linear theory of Section 3 differs from the normal form theories of the present section in that the former is computed around the new, post-impact, equilibrium, of radius r eq,new (see Section 3), while the latter introduce corrections to a theoretical equilibrium of the Keplerian problem at r = r ∗ (see Subsection 4.1). This implies that the linear theory of the previous section can be better in accuracy than the low order normal form theories (e.g. N = 1 , 2). However, normal form theories of higher truncation orders eventually take over as regards the accuracy. Below we examine the accuracy of the series solutions of orders N = 1 , 2 , 4 , 6. A visual comparison, in particular of Figures 5, 6 and 7, shows that the nonlinear theories can fit the waveforms of the curves of the numerical solutions, better than the linear theory, in most cases already at truncation N = 2. A more detailed comparison is shown in Figures 9 and 10 (errors in r ( t ) on the left and errors in ϕ 2 ( t ) on the right). As regards the 'Linear-Taylor' and the N = 2 theories, in general the latter performs better than the former, but with noticeable exceptions (e.g. in the curves r ( t ), ϕ 2 ( t ) in Figure 9 for β = 1). We attribute this to the better definition of the equilibrium state for the linear theory, which (depending on the magnitude of the perturbation and the physical parameters) can only be superseded by orders higher than second in perturbation theory.  \nIt is clear that the systematic growth of errors in Figures 9 and 10 reflects mostly the error in the determination of the fundamental frequencies ω 1 ,num and ω 2 ,num through any of the examined theories. The latter can be observed in Figure 11 which depicts the frequencies themselves as functions of the β parameter. Apart from the black curve (numerical), all others refer to the frequency evaluation by each of the  \nFig. 9 Absolute differences between the numerical and analytical solutions of the functions r ( t ) and ϕ 2 ( t ) for a total duration of 100 days. For the first set of physical parameters, the cases for β = 1 , 3 are shown.  \n<!-- image -->  \ntheories in consideration. The values that correspond to the black curves (numerical), are obtained by performing Laskar's NAFF (Laskar, J. et al., 1992) on the timeseries of Figures 4 to 7. We now see more clearly the improvement of the precision with the order of the theory. Note, in particular, that for the examined range of β , we typically have small relative errors, implying that the corresponding analytical formulas can be used to fit β from observed time series of length, even as high as of several months or years. Linear theory, instead, can be used for fitting β from observed time series of a limited span of few days.  \nAn estimate of the convergence properties, as well as the error of the normal form approximation, can be obtained without comparison with the numerical simulations, by comparing several quantities, depending on the order of the normal form truncation, as shown in Figure 12. This figure shows the relative differences in the prediction  \nFig. 10 Same plots as in Figure 9, but for the second set of physical parameters.  \n<!-- image -->  \nof the analytical frequencies ω 1 , ω 2 , as derived from the 4 th and 6 th order normal form theories, for the two parameter sets, plotted against β and the secondary asteroid's asphericity (1 -b 2 /a 2 ). The figure serves to distinguish the regions where the theory correctly predicts (up to a predefined accuracy) the frequencies, from those where it does not. Although the check performed in Figure 12 can be repeated for a variety of physical parameters which enter into the problem, such as the masses and the moments, our focus on β and the secondary's asphericity stems from their critical influence on the dynamics of the DART experiment. The yellow regions in the colormap indicate areas of high frequency errors. These regions correspond to resonances between the natural frequencies of the system, which in turn may lead to chaotic evolution (Agrusa et al. (2021)). The polynomial expressions (46) diverge at these parameter regions, causing the predicted frequencies to become unbounded. As a result, the relative differences increase significantly, indicating the overall failure of  \nthe approximation based on the perturbation theory. Figure 13 depicts the cumulative distributions of the percentage of the area (grid points) of the parameter space of Figure 12, where the relative error of the method is smaller than the value indicated in the abscissa. By selecting an upper bound for the error on the x -axis, one may determine the proportion of colormap grid points, represented on the y -axis, that fall within this error threshold. As shown in these plots, for both our chosen sets of parameters, the normal form method produces a relative error of 10 -5 or smaller for about 20% of the grid points, while this percentage grows to about 80% accepting a relative error smaller than 10 -2 .", '5 Conclusions': 'In the present paper we developed linear and nonlinear (canonical perturbation) theories able to describe by explicit formulas in time the evolution of the librational state of a binary deflected via a DART-like hit. We also gave examples of application of such theories for parameter values pertinent to the case of the 65803 Didymos system. The main steps and conclusions in our study are given below.  \nWe assumed a second-moment multipole expansion of the mutual potential and introduced two basic assumptions: planar motion and axisymmetry of the primary. The first of these assumptions can be removed, at the cost of larger complexity in the computations, as will be discussed in a forthcoming work. The second assumption can also be removed by adding the degree of freedom ϕ 1 to the system. However, we argued that this is not really necessary when the primary asteroid is a fast rotator, since, then, at formal level the assumption of axisymmetry is similar to first order averaging theory with respect to the fast angle ϕ 1 . On the other hand, the formulas introduced here are straightforward to generalize in the genuine case of a triaxial primary.  \nIn all examined theories, the corresponding explicit formulas representing the time series of any post-impact observable quantity of the trajectories (i.e. depending on r ( t ) and ϕ 2 ( t )) were given in parametric form in the momentum transfer parameter β . This is crucial for the above-reported methods to be of utility in data fitting obtained by observations of the system after the impact. Of course, the obvious advantage of such an analytical approach is that an analytical formula can be introduced at nearly zero computational cost in the least-square fitting algorithm, and can lead even to direct formulas for the optimal estimator ˆ β of β , while, on the contrary, any fitting method based on numerical trajectories requires running the whole trajectory anew per trial value of β . We may also imagine several additional applications of the analytical formulas, i.e., in the modelling of Earth-observed light curves of the binary asteroid, the detection of non-uniform motions, variations in the orbital eccentricity, etc.  \nAs would be expected, nonlinear normal form theories give better levels of precision as the order of the theory increases. On the other hand, truncation beyond a certain order, quickly becomes non-practical, as the number of series terms quickly proliferates. We found that our theories give small errors in the determination of the fundamental frequencies of the librational and epicyclic motion, which are the main factors affecting the systematic growth of errors in the curves r ( t ) and ϕ 2 ( t ), caused by systematic drifts in the phase of the oscillation. On the other hand, even linear theory  \nFig. 11 Comparison of the frequencies derived by NAFF ( ω 1 ,num , ω 2 ,num ) with those of analytical theories ( ω 1 , ω 2 , ω 1 c , ω 2 c ) for different values of β .  \n<!-- image -->', 'Parameter set 1': "<!-- image -->  \n<!-- image -->  \nParameter set 2  \n<!-- image -->  \nFig. 12 Relative errors in the analytical frequencies, as derived via the 4 th and the 6 th order theory.  \n<!-- image -->  \ncan give small errors in the frequencies determination, provided that the determination of the initial synchronous state around which the theory is developed is precise. Finally, the errors grow in general with β , as larger β implies in general larger postimpact libration amplitudes for the system, hence, a larger importance of the nonlinear corrections to the equations of motion 2 .  \nFig. 13 Cumulative distributions of the relative errors of the analytical frequencies that correspond to Figure 12.  \n<!-- image -->  \nAcknowledgements. M.G., I.G., G.V. and K.T acknowledge funding support from the European Union's Horizon 2020 research and innovation program under grant agreement No. 581870377 (project NEO-MAPP).  \nConflicts of interest. M.G., I.G. and G.V. have no conflicts of interest to declare that are relevant to the content of this article. C.E. and K.T. are Editorial Board Members of Celestial Mechanics and Dynamical Astronomy.", 'Appendix A Nonzero elements of the Jacobian matrix J': 'j 13 = 1 m (A1)  \nj 21 = √ 2 G ( 3 M 2 ( I 1 z -I s ) + 3 M 1 ( -2 I 2 x + I 2 y + I 2 z ) + 2 M 1 M 2 r 2 eq,new ) √ mr 7 eq,new (A2)  \nj 24 = 1 I 2 z + 1 mr 2 eq,new (A3)  \nj 31 = G ( 3 M 2 ( I 1 z -I s ) + 3 M 1 ( -2 I 2 x + I 2 y + I 2 z ) -2 M 1 M 2 r 2 eq,new ) 2 r 5 eq,new (A4)  \nj 42 = 3 GM 1 ( I 2 x -I 2 y ) r 3 eq,new (A5)', 'Appendix B Coefficients of the characteristic polynomial': 'ξ 2 = 3 G ( M 2 ( I s -I 1 z ) + M 1 ( I 2 y -I 2 z ) ) 2 mr 5 eq,new + G ( M 1 ( -3 I 2 x +3 I 2 y + I 2 z ) + M 2 I 2 z ) I 2 z r 3 eq,new (B6)  \nξ 0 = -3 G 2 ( I 2 x -I 2 y )( M 1 + M 2 ) 2 I 2 z M 1 M 2 2 r 10 eq,new × [ -3 M 1 M 2 r 2 eq,new ( M 2 ( I 1 z +2 I 2 z -I s ) + M 1 ( -2 I 2 x + I 2 y +3 I 2 z ) ) -15 I 2 z ( M 1 + M 2 )( M 2 ( I 1 z -I s ) + M 1 ( -2 I 2 x + I 2 y + I 2 z ) ) +2 M 2 1 M 2 2 r 4 eq,new ] (B7)', 'Appendix C Amplitude coefficients of the linear solution': 'A 1 r = δp ϕ 2 (0) j 13 j 21 -ω 2 2 δr (0) -j 13 j 31 δr (0) ω 2 1 -ω 2 2 (C8)  \nA 2 r = -δp ϕ 2 (0) j 13 j 21 + ω 2 1 δr (0) + j 13 j 31 δr (0) ω 2 1 -ω 2 2 (C9)', 'Appendix D DART-like disturbance at t = 0': 'δr (0) = r eq -r eq,new (D10)  \nδp ϕ 2 (0) = [ ˙ θ eq ( r eq ) -˙ θ eq,new ( r eq,new )] I 2 z (D11)  \nδp r (0) = δϕ 2 (0) = 0 (D12)', 'Appendix E Evaluation of the functions θ ( t ) and ϕ 1 ( t )': 'At any moment, the orbital displacement rate ˙ θ is given by  \n˙ θ = p θ -p ϕ 1 -p ϕ 2 ( t ) mr ( t ) 2 ⇒ θ ( t ) = ∫ p θ -p ϕ 1 -p ϕ 2 ( t ) mr ( t ) 2 dt (E13)  \nwhere p θ and p ϕ 1 are constants in time. Assuming a DART-like excitation, the linear solutions are those of equations (23). Substituting the latter in (E13), yields  \nθ ( t ) = ∫ p θ -p ϕ 1 -( p ϕ 2 eq,new + A 1 p ϕ 2 cos ( ω 1 t ) + A 2 p ϕ 2 cos ( ω 2 t ) ) m ( r eq,new + A 1 r cos ( ω 1 t ) + A 2 r cos ( ω 2 t ) ) 2 dt (E14)  \nwhere the coefficients A 1 p ϕ 2 , A 2 p ϕ 2 are assumed to be the amplitudes of δp ϕ 2 ( t ) of equations (22), i.e.  \nA 1 p ϕ 2 = ( j 13 j 31 + ω 2 1 ) A 1 r j 13 j 21 , A 2 p ϕ 2 = ( j 13 j 31 + ω 2 2 ) A 2 r j 13 j 21 . (E15)  \nSince the direct evaluation of the integral of equation (E14) proves challenging, we rewrite it as  \nθ ( t ) = ∫ p θ -p ϕ 1 -( p ϕ 2 eq,new + A 1 p ϕ 2 cos ( ω 1 t ) + A 2 p ϕ 2 cos ( ω 2 t ) ) m ( r eq,new + δr ( t ) ) 2 dt . (E16)  \nSubstituting δr ← ϵ · δr , expanding up to first order terms in powers of ϵ and back substituting the function δr ( t ), renders the integral analytically computable. By separating the polynomial and trigonometric terms, we obtain  \nθ ( t ; β ) = θ (0) + ( A 1 r A 1 p ϕ 2 + A 2 r A 2 p ϕ 2 + mr 3 eq,new ˙ θ eq,new mr 3 eq,new ) t + -2 A 1 r mr eq,new ˙ θ eq,new -A 1 p ϕ 2 mr 2 eq,new ω 1 sin ( ω 1 t ) + -2 A 2 r mr eq,new ˙ θ eq,new -A 2 p ϕ 2 mr 2 eq,new ω 2 sin ( ω 2 t ) + A 1 r A 1 p ϕ 2 2 mr 3 eq,new ω 1 sin (2 ω 1 t ) + A 2 r A 2 p ϕ 2 2 mr 3 eq,new ω 2 sin (2 ω 2 t ) + A 1 r A 2 p ϕ 2 + A 2 r A 1 p ϕ 2 mr 3 eq,new ( ω 1 -ω 2 ) sin [ ( ω 1 -ω 2 ) t ] + A 1 r A 2 p ϕ 2 + A 2 r A 1 p ϕ 2 mr 3 eq,new ( ω 1 + ω 2 ) sin [ ( ω 1 + ω 2 ) t ] . (E17)  \nThe coefficient of each of the above terms corresponds to A jθ of equation (24). On the other hand ˙ ϕ 1 , obeys the equation  \n˙ ϕ 1 = p ϕ 1 I 1 z -˙ θ . (E18)  \nDirect integration yields  \nϕ 1 ( t ; β ) = ϕ 1 (0) + p ϕ 1 I 1 z t -θ ( t ; β ) . (E19)', 'Appendix F Full form of the terms Z 0 , h 1 and h 2': 'Z 0 = 1 2 ( ν 2 1 I 1 z + G ( M 1 (2 I 2 x -I 2 y -M 2 r ∗ 2 ) + M 2 ( -I 1 z + I 2 z + I s ) ) r ∗ 3 ) + ν ∗ θ δp θ +( ν 1 -ν ∗ θ ) δp ϕ 1 + 1 2 m δp 2 r + ( 1 2 I 2 z + 1 2 mr ∗ 2 ) δp 2 ϕ 2 + 2 ν ∗ θ r ∗ 2 δrδp ϕ 2 + G ( 6 M 2 ( I s -I 1 z ) + 6 M 1 (2 I 2 x -I 2 y -I 2 z ) + M 1 M 2 r ∗ 2 ) 2 r ∗ 5 δr 2 -3 GM 1 ( I 2 x -I 2 y ) 2 r ∗ 3 δϕ 2 2 (F20)  \nh 1 = δp ϕ 1 -δp θ mr ∗ 2 δp ϕ 2 -3 ν ∗ θ r ∗ 3 δp ϕ 2 δr 2 + G ( 5 M 2 ( I 1 z -I s ) + 5 M 1 ( -2 I 2 x + I 2 y + I 2 z ) -M 1 M 2 r ∗ 2 ) r ∗ 6 δr 3 + ( 3 G ( M 2 ( I 1 z -I s ) + M 1 ( -2 I 2 x + I 2 y + I 2 z ) ) 2 r ∗ 4 + 2 ν ∗ θ ( δp ϕ 1 -δp θ ) r ∗ 2 ) δr -1 mr ∗ 3 δrδp 2 ϕ 2 + 9 GM 1 ( I 2 x -I 2 y ) 2 r ∗ 4 δrδϕ 2 2 (F21)  \nh 2 = I 1 z ( δp θ -δp ϕ 1 ) 2 + mr ∗ 2 δp 2 ϕ 1 2 I 1 z mr ∗ 2 + 2( δp θ -δp ϕ 1 ) mr ∗ 3 δrδp ϕ 2 + 4 ν ∗ θ r ∗ 4 δr 3 δp ϕ 2 + 3 G ( 5 M 2 ( I s -I 1 z ) + 5 M 1 (2 I 2 x -I 2 y -I 2 z ) + M 1 M 2 r ∗ 2 ) 2 r ∗ 7 δr 4 + GM 1 ( I 2 x -I 2 y ) 2 r ∗ 3 δϕ 4 2 + 3 ν ∗ θ ( δp θ -δp ϕ 1 ) r ∗ 3 δr 2 + 3 2 mr ∗ 4 δr 2 δp 2 ϕ 2 -9 GM 1 ( I 2 x -I 2 y ) r ∗ 5 δr 2 δϕ 2 2 (F22)', 'Appendix G Conversion to Birkhoff complex canonical variables ( Q 1 , Q 2 , P 1 , P 2 )': "We write down the Hamilton's equations that correspond to the term Z 0 of equation (36)  \n˙ δr = ∂Z 0 ∂δp r ˙ δp r = -∂Z 0 ∂δr ˙ δϕ 2 = ∂Z 0 ∂δp ϕ 2 ˙ δp ϕ 2 = -∂Z 0 ∂δϕ 2 (G23)  \nand then we evaluate the corresponding Jacobian matrix  \nK =     0 0 k 13 0 k 21 0 0 k 24 k 31 0 0 -k 21 0 k 42 0 0     (G24)  \nwhere  \nk 13 = 1 m (G25)  \nk 21 = 2 ν ∗ θ r ∗ 2 (G26)  \nk 24 = 1 I 2 z + 1 mr ∗ 2 (G27)  \nk 31 = G ( 6 M 2 ( I 1 z -I s ) + 6 M 1 ( -2 I 2 x + I 2 y + I 2 z ) -M 1 M 2 r ∗ 2 ) r ∗ 5 (G28)  \nk 42 = 3 GM 1 ( I 2 x -I 2 y ) r ∗ 3 . (G29)  \nK resembles the Jacobian J of equation (19), but in fact K is evaluated for the twospheres equilibrium point of equation (29). The eigenvalues of K are obtained from the characteristic polynomial  \nΠ( λ ) = λ 4 + ζ 2 λ 2 + ζ 0 (G30)  \nζ 2 = GM 1 r ∗ 2 I 2 z mr ∗ 5 × [ ( -3 I 2 x m +3 I 2 y m + I 2 z M 2 ) -3 GI 2 z ( 2 M 2 ( I 1 z -I s ) + M 1 ( -3 I 2 x + I 2 y +2 I 2 z ) ) ] (G31)  \nζ 0 = 3 G 2 M 1 ( I 2 x -I 2 y ) I 2 z m 2 r ∗ 10 × [ 4 M 1 M 2 I 2 z r ∗ 2 +( I 2 z + mr ∗ 2 ) ( 6 M 2 ( I 1 z -I s ) +6 M 1 ( -2 I 2 x + I 2 y + I 2 z ) -M 1 M 2 r ∗ 2 ) ] (G32)  \nThe roots are λ 1 , 2 = ∓ iω 1 k , λ 3 , 4 = ∓ iω 2 k and ω 1 k , ω 2 k are the kernel frequencies.  \nω 1 k, 2 k = √ ζ 2 ± √ ζ 2 2 -4 ζ 0 √ 2 . (G33)  \nThe corresponding eigenvectors are  \nV =     ⃗ v 1 ⃗ v 2 ⃗ v 3 ⃗ v 4     =     f ( λ 1 ) g ( λ 1 ) h ( λ 1 ) 1 f ( λ 2 ) g ( λ 2 ) h ( λ 2 ) 1 f ( λ 3 ) g ( λ 3 ) h ( λ 3 ) 1 f ( λ 4 ) g ( λ 4 ) h ( λ 4 ) 1     (G34)  \nwhere  \nf ( λ ) = k 13 k 21 k 13 k 31 -λ 2 g ( λ ) = -k 13 ( k 2 21 + k 24 k 31 ) + k 24 λ 2 -k 13 k 31 λ + λ 3 h ( λ ) = k 21 λ k 13 k 31 -λ 2 . (G35)  \nThe transformation between the variables ( δr, δϕ 2 , δp r , δp ϕ 2 ) and ( Q 1 , Q 2 , P 1 , P 2 ) is canonical if the symplectic condition is satisfied, i.e.  \nM T SM = S (G36)  \nwhere M is the 4 × 4 matrix  \nM =     P 1 f ( λ 1 ) P 2 f ( λ 3 ) -P 1 f ( λ 2 ) -P 2 f ( λ 4 ) P 1 g ( λ 1 ) P 2 g ( λ 3 ) -P 1 g ( λ 2 ) -P 2 g ( λ 4 ) P 1 h ( λ 1 ) P 2 h ( λ 3 ) -P 1 h ( λ 2 ) -P 2 h ( λ 4 ) P 1 P 2 -P 1 -P 2     (G37)  \nand S is the 4 × 4 symplectic matrix  \nS =     0 0 1 0 0 0 0 1 -1 0 0 0 0 -1 0 0     . (G38)  \nBy substituting (G37) into (G36) and recalling that λ 2 = -λ 1 , λ 4 = -λ 3 , f ( -λ ) = f ( λ ), g ( -λ ) = -g ( λ ) and h ( -λ ) = -h ( λ ), we obtain the two equations  \n2 ( g ( λ 1 ) -f ( λ 1 ) h ( λ 1 ) ) P 2 1 +1 = 0 , 2 ( g ( λ 3 ) -f ( λ 3 ) h ( λ 3 ) ) P 2 2 +1 = 0 . (G39)  \nWe choose one of the corresponding set of roots  \nP 1 = -i u ( λ 1 ) , P 2 = -i u ( λ 3 ) , where u ( λ ) = √ 2 g ( λ ) -2 f ( λ ) h ( λ ) (G40)  \nSubstituting (G40) back at (G37), yields  \nM ≡     m 11 m 12 m 13 m 14 m 21 m 22 m 23 m 24 m 31 m 32 m 33 m 34 m 41 m 42 m 43 m 44     =             -if ( λ 1 ) u ( λ 1 ) -if ( λ 3 ) u ( λ 3 ) if ( λ 1 ) u ( λ 1 ) if ( λ 3 ) u ( λ 3 ) -ig ( λ 1 ) u ( λ 1 ) -ig ( λ 3 ) u ( λ 3 ) -ig ( λ 1 ) u ( λ 1 ) -ig ( λ 3 ) u ( λ 3 ) -ih ( λ 1 ) u ( λ 1 ) -ih ( λ 3 ) u ( λ 3 ) -ih ( λ 1 ) u ( λ 1 ) -ih ( λ 3 ) u ( λ 3 ) -i u ( λ 1 ) -i u ( λ 3 ) i u ( λ 1 ) i u ( λ 3 )             (G41)  \nand the inverse matrix M -1 is  \nM -1 =              iu ( λ 1 ) q 1 ih ( λ 3 ) u ( λ 1 ) q 2 ig ( λ 3 ) u ( λ 1 ) q 3 -if ( λ 3 ) u ( λ 1 ) q 4 -iu ( λ 3 ) q 1 ih ( λ 1 ) u ( λ 3 ) q 3 ig ( λ 1 ) u ( λ 3 ) q 2 if ( λ 1 ) u ( λ 3 ) q 4 -iu ( λ 1 ) q 1 ih ( λ 3 ) u ( λ 1 ) q 2 ig ( λ 3 ) u ( λ 1 ) q 3 if ( λ 3 ) u ( λ 1 ) q 4 iu ( λ 3 ) q 1 ih ( λ 1 ) u ( λ 3 ) q 3 ig ( λ 1 ) u ( λ 3 ) q 2 -if ( λ 1 ) u ( λ 3 ) q 4              (G42)  \nwhere  \nq 1 = 2 ( f ( λ 1 ) -f ( λ 3 ) ) q 2 = 2 ( -g ( λ 3 ) h ( λ 1 ) + g ( λ 1 ) h ( λ 3 ) ) q 3 = 2 ( g ( λ 3 ) h ( λ 1 ) -g ( λ 1 ) h ( λ 3 ) ) q 4 = 2 ( f ( λ 1 ) -f ( λ 3 ) )  \nUltimately, M and M -1 set the transformations (forward and inverse) between ( δr, δϕ 2 , δp r , δp ϕ 2 ) and ( Q 1 , Q 2 , P 1 , P 2 ) as functions of the physical parameters and r ∗ .", 'Appendix H Evaluation of the normal form Hamiltonian': 'The Birkhoff normalization process defines a sequence of N consecutive near-toidentity canonical transformations to the variables F = ( Q 1 , Q 2 , P 1 , P 2 ), namely  \nF ≡ F (0) → F (1) → ... → F ( N ) , (H43)  \nsuch that, in the variables F ( N ) , the Hamiltonian, resumes the form  \nThe quantity  \nZ ( N ) = Z 0 + ϵZ 1 + ... + ϵ N Z N , (H45)  \nis in Birkhoff normal form, implying that Z ( N ) is integrable. The quantity  \nR ( N ) = ∞ ∑ j = N +1 ϵ j H ( N ) j , (H46)  \nis called the remainder function. Various theoretical arguments (see e.g. Efthymiopoulos (2012)) imply that the remainder yields, in general, a small correction which can be ignored in the equations of motion. The functional form of the transformations (H43) is obtained by the method of composition of Lie series (see Appendix I). Up to the N th order we get  \nF ≡ F (0) = exp( L χ N ) exp ( L χ N -1 ) · · · exp( L χ 1 ) F ( N ) (H47)  \nwhich expresses the variables F (0) in terms of F ( N ) . The operator exp ( L χ ), acting on an arbitrary function f , is defined as  \nexp( L χ ) f = f + L χ f + L 2 χ f 2 + L 3 χ f 6 + · · · , (H48)  \nwhere L χ denotes the Poisson bracket operator, L χ f = { f, χ } . The functions χ 1 , ..., χ N are called Lie generating functions and at each normalization step, they are evaluated through the solution of the homological equation  \n{ Z 0 , χ n +1 } + ϵ n +1 h ( n ) n +1 = 0 , (H49)  \nwhere Z 0 is the kernel of the normal form, and h ( n ) n +1 contains all the terms of bookkeeping order n + 1 in the Hamiltonian H ( n ) obtained after n normalization steps, for which the Poisson bracket with Z 0 is non-zero. After performing N steps of the normalization algorithm, we obtain the normal form Hamiltonian (41). Note that during the normalization steps, all terms of book-keeping order greater than N are truncated ( O ( ϵ q>N ) = 0).  \nIn order to map the solution (44) into the original variables ( δr ( t ) , δϕ 2 ( t ) , δp r ( t ) , δp ϕ 2 ( t ) ) and ( θ ( t ) , ϕ 1 ( t ) ) , we implement the following two steps: 1) Assuming an initial disturbance ( δr (0) = 0, δϕ 2 (0) = 0, δp r (0) = 0, δp θ = p θ imp -p ∗ θ , δp ϕ 1 = 0, δp ϕ 2 (0) = p ϕ 2 imp -p ∗ ϕ 2 ) , we use the inverse transformation (37) to compute the initial conditions F (0) (0), and then the inverse of the Lie series (H47) (see Appendix I)  \nF ( N ) = exp( -L χ 1 ) · · · exp( -L χ N -1 ) exp ( -L χ N ) F (0) (H50)  \nin order to compute the initial conditions F (4) (0). By doing so, we can express ω 1 c , ω 2 c , ω θc , ω ϕ 1 c as functions of β , as noted in equation (46). 2) We transform the  \nH ( N ) = ( Z 0 + ϵZ 1 + ... + ϵ N Z N ) + ∞ ∑ j = N +1 ϵ j H ( N ) j . (H44)  \nwhole solution in time back to the original variables. This is accomplished by utilizing the Lie series (H47) to convert equation (44) to F (0) ( t ) and ( θ ( t ) , ϕ 1 ( t ) ) . Finally the transformation (37) yields the desired functions ( δr ( t ) , δϕ 2 ( t ) , δp r ( t ) , δp ϕ 2 ( t ) ) .', 'Appendix I Composition of Lie series up to fourth order': 'e L χ 4 e L χ 3 e L χ 2 e L χ 1 f = ( 1 + L χ 1 + L 2 χ 1 2 + L 3 χ 1 6 + L 4 χ 1 24 + L χ 2 + L χ 2 L χ 1 + L χ 2 L 2 χ 1 2 + L 2 χ 2 2 + L χ 3 + L χ 3 L χ 1 + L χ 4 ) f = f + { f, χ 1 } + 1 2 {{ f, χ 1 } , χ 1 } + 1 6 {{{ f, χ 1 } , χ 1 } , χ 1 } + 1 24 {{{{ f, χ 1 } , χ 1 } , χ 1 } , χ 1 } + { f, χ 2 } + {{ f, χ 1 } , χ 2 } + 1 2 {{{ f, χ 1 } , χ 1 } , χ 2 } + 1 2 {{ f, χ 2 } , χ 2 } + { f, χ 3 } + {{ f, χ 1 } , χ 3 } + { f, χ 4 } (I51)  \ne -L χ 1 e -L χ 2 e -L χ 3 e -L χ 4 f = ( 1 -L χ 1 + L 2 χ 1 2 -L 3 χ 1 6 + L 4 χ 1 24 -L χ 2 + L χ 1 L χ 2 -L 2 χ 1 2 L χ 2 + L 2 χ 2 2 -L χ 3 + L χ 1 L χ 3 -L χ 4 ) f = f -{ f, χ 1 } + 1 2 {{ f, χ 1 } , χ 1 } -1 6 {{{ f, χ 1 } , χ 1 } , χ 1 } + 1 24 {{{{ f, χ 1 } , χ 1 } , χ 1 } , χ 1 } - { f, χ 2 } + {{ f, χ 2 } , χ 1 } -1 2 {{{ f, χ 2 } , χ 1 } , χ 1 } + 1 2 {{ f, χ 2 } , χ 2 } - { f, χ 3 } + {{ f, χ 3 } , χ 1 } - { f, χ 4 } (I52)', 'References': "Agrusa, H.F., Gkolias, I., Tsiganis, K., Richardson, D.C., Meyer, A.J., Scheeres, D.J., ' Cuk, M., Jacobson, S.A., Michel, P., Karatekin, O., Cheng, A.F., Hirabayashi, M., Zhang, Y., Fahnestock, E.G., Davis, A.B.: The excited spin state of Dimorphos resulting from the DART impact. Icarus 370 , 114624 (2021) https://doi.org/10. 1016/j.icarus.2021.114624  \nBou'e, G., Laskar, J.: Spin axis evolution of two interacting bodies. Icarus 201 (2), 750-767 (2009) https://doi.org/10.1016/j.icarus.2009.02.001  \nCheng, A.F., Agrusa, H.F., Barbee, B.W., Meyer, A.J., Farnham, T.L., Raducan, S.D., Richardson, D.C., Dotto, E., Zinzi, A., Della Corte, V., Statler, T.S., Chesley, S., Naidu, S.P., Hirabayashi, M., Li, J.-Y., Eggl, S., Barnouin, O.S., Chabot, N.L., Chocron, S., Collins, G.S., Daly, R.T., Davison, T.M., DeCoster, M.E., Ernst, C.M., Ferrari, F., Graninger, D.M., Jacobson, S.A., Jutzi, M., Kumamoto, K.M., Luther, R., Lyzhoft, J.R., Michel, P., Murdoch, N., Nakano, R., Palmer, E., Rivkin, A.S., Scheeres, D.J., Stickle, A.M., Sunshine, J.M., Trigo-Rodriguez, J.M., Vincent, J.-B., Walker, J.D., Wunnemann, K., Zhang, Y., Amoroso, M., Bertini, I., Brucato, J.R., Capannolo, A., Cremonese, G., Dall'Ora, M., Deshapriya, P.J.D., Gai, I., Hasselmann, P.H., Ieva, S., Impresario, G., Ivanovski, S.L., Lavagna, M., Lucchetti, A., Epifani, E.M., Modenini, D., Pajola, M., Palumbo, P., Perna, D., Pirrotta, S., Poggiali, G., Rossi, A., Tortora, P., Zannoni, M., Zanotti, G.: Momentum transfer from the DART mission kinetic impact on asteroid Dimorphos. Nature 616 (7957), 457-460 (2023) https://doi.org/10.1038/s41586-023-05878-z  \nDaly, R.T., Ernst, C.M., Barnouin, O.S., Chabot, N.L., Rivkin, A.S., Cheng, A.F., Adams, E.Y., Agrusa, H.F., Abel, E.D., Alford, A.L., Asphaug, E.I., Atchison, J.A., Badger, A.R., Baki, P., Ballouz, R.-L., Bekker, D.L., Bellerose, J., Bhaskaran, S., Buratti, B.J., Cambioni, S., Chen, M.H., Chesley, S.R., Chiu, G., Collins, G.S., Cox, M.W., DeCoster, M.E., Ericksen, P.S., Espiritu, R.C., Faber, A.S., Farnham, T.L., Ferrari, F., Fletcher, Z.J., Gaskell, R.W., Graninger, D.M., Haque, M.A., Harrington-Duff, P.A., Hefter, S., Herreros, I., Hirabayashi, M., Huang, P.M., Hsieh, S.-Y.W., Jacobson, S.A., Jenkins, S.N., Jensenius, M.A., John, J.W., Jutzi, M., Kohout, T., Krueger, T.O., Laipert, F.E., Lopez, N.R., Luther, R., Lucchetti, A., Mages, D.M., Marchi, S., Martin, A.C., McQuaide, M.E., Michel, P., Moskovitz, N.A., Murphy, I.W., Murdoch, N., Naidu, S.P., Nair, H., Nolan, M.C., Ormo, J., Pajola, M., Palmer, E.E., Peachey, J.M., Pravec, P., Raducan, S.D., Ramesh, K.T., Ramirez, J.R., Reynolds, E.L., Richman, J.E., Robin, C.Q., Rodriguez, L.M., Roufberg, L.M., Rush, B.P., Sawyer, C.A., Scheeres, D.J., Scheirich, P., Schwartz, S.R., Shannon, M.P., Shapiro, B.N., Shearer, C.E., Smith, E.J., Steele, R.J., Steckloff, J.K., Stickle, A.M., Sunshine, J.M., Superfin, E.A., Tarzi, Z.B., Thomas, C.A., Thomas, J.R., Trigo-Rodr'ıguez, J.M., Tropf, B.T., Vaughan, A.T., Velez, D., Waller, C.D., Wilson, D.S., Wortman, K.A., Zhang, Y.: Successful kinetic impact into an asteroid for planetary defence. Nature 616 (7957), 443-447 (2023) https://doi.org/10.1038/s41586-023-05810-5  \nDeprit, A.: Canonical transformations depending on a small parameter. Celestial Mechanics 1 (1), 12-30 (1969)  \nDavis, A. B., Scheeres, D. J. : Doubly synchronous binary asteroid mass parameter observability. Icarus 341 , 113439 (2020) https://doi.org/10.1016/j.icarus.2019. 113439  \nEfthymiopoulos, C.: Canonical perturbation theory; stability and diffusion in Hamiltonian systems: applications in dynamical astronomy. Workshop Series of the Asociacion Argentina de Astronomia 3 , 3-146 (2012)  \nRivkin, A.S., Chabot, N.L., Stickle, A.M., Thomas, C.A., Richardson, D.C., Barnouin, O., Fahnestock, E.G., Ernst, C.M., Cheng, A.F., Chesley, S., Naidu, S., Statler, T.S., Barbee, B., Agrusa, H., Moskovitz, N., Terik Daly, R., Pravec, P., Scheirich, P.,  \nDotto, E., Della Corte, V., Michel, P., Kuppers, M., Atchison, J., Hirabayashi, M.: The Double Asteroid Redirection Test (DART): Planetary Defense Investigations and Requirements. Planetary Science Journal 2 (5), 173 (2021) https://doi.org/10. 3847/PSJ/ac063e  \nRaducan, S.D., Jutzi, M., Cheng, A.F., Zhang, Y., Barnouin, O., Collins, G.S., Daly, R.T., Davison, T.M., Ernst, C.M., Farnham, T.L., Ferrari, F., Hirabayashi, M., Kumamoto, K.M., Michel, P., Murdoch, N., Nakano, R., Pajola, M., Rossi, A., Agrusa, H.F., Barbee, B.W., Syal, M.B., Chabot, N.L., Dotto, E., Fahnestock, E.G., Hasselmann, P.H., Herreros, I., Ivanovski, S., Li, J.-Y., Lucchetti, A., Luther, R., Ormo, J., Owen, M., Pravec, P., Rivkin, A.S., Robin, C.Q., S'anchez, P., Tusberti, F., Wunnemann, K., Zinzi, A., Epifani, E.M., Manzoni, C., May, B.H.: Physical properties of asteroid Dimorphos as derived from the DART impact. Nature Astronomy (2024) https://doi.org/10.1038/s41550-024-02200-3  \nScheeres, D.J.: Stability of the planar full 2-body problem. Celestial Mechanics and Dynamical Astronomy 104 (1-2), 103-128 (2009) https://doi.org/10.1007/ s10569-009-9184-7  \nThomas, C.A., Naidu, S.P., Scheirich, P., Moskovitz, N.A., Pravec, P., Chesley, S.R., Rivkin, A.S., Osip, D.J., Lister, T.A., Benner, L.A.M., Brozovi'c, M., Contreras, C., Morrell, N., Ro ˙ Zek, A., Kuˇsnir'ak, P., Hornoch, K., Mages, D., Taylor, P.A., Seymour, A.D., Snodgrass, C., Jørgensen, U.G., Dominik, M., Skiff, B., Polakis, T., Knight, M.M., Farnham, T.L., Giorgini, J.D., Rush, B., Bellerose, J., Salas, P., Armentrout, W.P., Watts, G., Busch, M.W., Chatelain, J., Gomez, E., Greenstreet, S., Phillips, L., Bonavita, M., Burgdorf, M.J., Khalouei, E., Longa-Pe˜na, P., Rabus, M., Sajadian, S., Chabot, N.L., Cheng, A.F., Ryan, W.H., Ryan, E.V., Holt, C.E., Agrusa, H.F.: Orbital period change of Dimorphos due to the DART kinetic impact. Nature 616 (7957), 448-451 (2023) https://doi.org/10.1038/s41586-023-05805-2"}