Brouwer's satellite solution redux

Brouwer's solution to the artificial satellite problem is revisited to show that the complete Hamiltonian reduction is rather achieved in the plain Poincar\'e's style, through a single canonical transformation, than using a sequence of partial reductions based on von Zeipel's alternative for dealing with perturbed degenerate Hamiltonian systems.


Introduction
Brouwer's analytical solution to the artificial satellite problem [7] based on von Zeipel's partial reduction method for dealing with perturbed degenerate Hamiltonians [42] fiercely resists obsolescence sixty years after publication. Indeed, in spite of the spectacular increase of the computational power, widespread software packages for approximate ephemeris prediction still rely on Brouwer's seminal results [18,41]. Furthermore, the success of Brouwer's closed-form solution among practitioners as well as the reputation gained among theorists by Brouwer's stepwise normalization, make that these days some authors designate the method as "the von Zeipel-Brouwer theory" [15].
Since then, the merits of Brouwer's decomposition of the solution of perturbed Keplerian motion into secular, long-, and short-period effects, seem not to have been questioned. Moreover, after the invention of Hamiltonian simplification methods [10], it has been suggested that carrying out additional decompositions, thus increasing the number of canonical transformations, could be the proper way to success in the search for separable perturbation Hamiltonians of celestial mechanics problems [11]. Conversely, it has been recently pointed out that the use of Hamiltonian simplification procedures could be merely optional in the construction of higher-order analytical approximations to the satellite problem [30]. Then, it emerges the question of which is the real value of splitting a normalization procedure, either partial or complete, into several different stages, a topic that may well deserve additional study.
We walk a step in that direction and, in the classical style of Poincaré's "new method" [37], undertake the construction of Brouwer's second-order completely reduced Hamiltonian of the main problem in the artificial satellite theory (the J 2 problem) by means of a single canonical transformation. The difficulties stemming from the degeneracy of the Kepler Hamiltonian, who has two null frequencies, are easily overcome with the addition of suitable integration constants to the generating function of the transformation that yields the complete Hamiltonian reduction.
The use of arbitrary functions in the construction of perturbation solutions is not new at all [35]. In fact, it can be traced back to Poincaré's efforts in approaching degenerate perturbation problems [37,Chap. XI]. They also play a fundamental role in a reverse partial normalization process in which the angular momentum is normalized in first place [3,28,38,30]. On the other hand, in spite of average perturbation Hamiltonians do not exist in general [14], the use of arbitrary constants to guarantee the purely periodic nature of the generating function became customary in attempts to bring the mean elements dynamics as close as possible to the true average dynamics [33,39,27].
To the first order, the construction of Brouwer's closed-form solution by means a single transformation amounts to the sum of the two transformations computed by Brouwer for the short-and long-period elimination. This is due to the linearization provided by the first order truncation of a perturbation theory. In view of no differences arise between Brouwer's and the current approach when the periodic corrections are constrained to first order effects, we feel compelled to supplement Brouwer's analytical solution with second-order periodic corrections, yet limited to the J 2 contributions. We compare our results with corresponding corrections obtained in the traditional way, in which the normalization is split into the elimination of the parallax, the elimination of the perigee, and the Delaunay normalization [3,26]. At the second order the single transformation is no longer the addition of the different canonical transformations. As expected, the (single) periodic corrections are now much more involved than those corresponding to each of the partial reductions or simplifications, and are also more intricate than the composition of all of them. However, the length of the series defining the solution is only a part of the whole picture, and we found clear computational advantages in the evaluation of the single transformation. The improvements stem from the fact that many inclination polynomials pertaining to the periodic corrections admit factorization. Because common factors repeat many times throughout the corrections, the com-piler is able to perform a higher optimization of the code in the case of the single transformation than in the case in which the transformation is split in different stages.
In the process of computing the second-order corrections, we will recognize how artificial the controversy created about the integration of the equation of the center was. Indeed, difficulties confronted by researchers involved in the automatization of celestial mechanics computations were, in fact, derived from their own programming strategies [12,21]. On the contrary, the trouble had been easily sidestepped by celestial mechanicians relying on traditional hand computations [23,1]. We will show that standard integration by parts reduces the equation of the center issue to the well known integration of cosine functions in elliptic motion [22,40]. We hasten to say that the controversy was in no way futile since it provoked the appearance of Hamiltonian simplification methods and led to the development of sophisticated computational strategies [16]. 1 In order to fully determine the second-order term of the generating function of Brouwer's theory, the third-order term of the completely reduced Hamiltonian needs to be previously specified. The use of higher-order secular terms should improve further the long-term performance of Brouwer's solution. However, to be effective in the propagation of an initial state vector, the initialization constants of the analytical solution, and, in particular, the secular frequencies, must be computed within comparable accuracy to that of the secular terms. Rather than carrying out the long and tedious computations required in the determination of the thirdorder generating function, we take the clever shortcut proposed by Breakwell and Vagners [6]. That is, we limit the computation of third-order corrections to the case of the secular mean motion, which, besides, is directly obtained from the secular Hamiltonian. With this effortless procedure the addition of third-order secular terms clearly improves the performance of Brouwer's solution.

Brouwer's complete reduction at once
Constraining the dynamical model of the artificial satellite problem to the J 2 perturbation (main problem), Brouwer's gravitational Hamiltonian takes the form [7] where the Earth's gravitational field is materialized by the physical constants µ, the gravitational parameter, R ⊕ , the equatorial radius, and C 2,0 = −J 2 , the non-dimensional oblateness coefficient. 2 The symbols a, r, f and ω, stand for semimajor axis, radius from the Earth's center of mass, true anomaly, and argument of the perigee, respectively, whereas s ≡ sin I abbreviates the sine of the inclination I. Since we are dealing with Hamiltonian mechanics, these symbols must be understood as functions of some set of canonical variables. In particular we assume, with Brouwer, that the Hamiltonian is written in terms of the Delaunay coordinates , the mean anomaly, g = ω, and h, the longitude of the ascending node, and their conjugate momenta L = √ µa, G = L(1 − e 2 ) 1/2 , with e denoting the eccentricity, and H = G cos I, standing for the Delaunay action, the total angular momentum, and the projection of the angular momentum vector along the Earth's rotation axis, respectively. That H is an integral of Eq. (1) becomes evident from the cyclic character of h. Besides, the Hamiltonian itself is constant because the time does not appear explicitly on it.
The small value of the Earth's J 2 coefficient identifies Eq. (1) like a case of perturbed Keplerian motion, which, therefore, can be reduced to a separable Hamiltonian by perturbation methods. This is achieved by finding a canonical transformation T : ( , g, h, L, G, H, ) → ( , g , h , L , G , H ), from osculating to mean variables, depending on the small parameter ∼ J 2 , such that the transformed Hamiltonian in mean (prime) variables becomes a function of only the momenta, namely H • T = K(−, −, −, L , G , H ; ). The transformation T , we learned from Poincaré [37], is derived from a determining function that is solved in the form of a Taylor series up to some truncation order of the small parameter .
Brouwer, for his part, after introducing the method of solution, seems to refuse approching the direct computation of the transformation T since the beginning, by simply declaring that "[. . . ] it is more convenient to choose a determining function in such a manner that the mean anomaly is not present in the transformed Hamiltonian while the argument of the perigee is permitted to appear." Next, after invoking von Zeipel, he proceeded stepwise by partial reduction, first computing a canonical transformation that only removes the short-period terms from the Hamiltonian, and then carrying out a second canonical transformation that removes the long-period terms. In this way Brouwer outstandingly achieves the complete Hamiltonian reduction in closed form.
Conversely, we ignore the presumed convenience of Brouwer's procedure and approach the perturbation problem searching for a single determining function in the original style of Poincaré, yet we better rely on the equivalent but more functional method of Lie transforms [19,9,13]. Thus, we write Eq. (1) in the usual form of a perturbation Hamiltonian H = m≥0 ( m /m!)H m,0 , with At each step m, terms H 0,m in Eq. (2) are known, coming either from the original Hamiltonian or stemming from intermediate computations at previous orders.
Terms H 0,m are selected in such a way that they cancel those terms of H 0,m pertaining to the kernel of the Lie operator. Finally, the integration "constants" C m -arbitrary functions of the Delaunay variables fulfilling the condition ∂C m /∂ = 0will be chosen like such trigonometric functions of g that they prevent the appearance of purely long-period terms at the next order of the perturbation approach, in this way making feasible the complete normalization at once. The method is standard these days, and the required details can be found in textbooks as, for instance, [34,4].
In preparation of the solution, the equivalence where p = aη 2 is the orbit parameter and η = (1 − e 2 ) 1/2 is the so-called eccentricity function, is applied to the instances j > 2 in H 0,1 = H 1,0 , which is then written in the convenient form where B 0 = 1 − 3 2 s 2 , B 1 = 3 4 s 2 , and we abbreviate j ≡ j mod 2. On account of j ≥ k in Eq. (4), we immediately verify that H 0,1 is not affected of purely long-period terms. Then, the complete reduction is achieved at the first order by choosing the new Hamiltonian term H 0,1 like the average of H 1,0 over the mean anomaly.
The average is obtained in closed form with the help of the Keplerian differential relation between the true and mean anomalies ηa 2 d = r 2 df . It is equivalent to removing all the terms with j > 0 from Eq. (4) after multiplied by the factor r 2 /(a 2 η). We trivially obtain the usual result which is, of course, the same expression obtained by Brouwer. Then, Eq. (2) is rearranged in the form where φ = f − denotes the equation of the center, and the integrand in Eq. (5) only embraces periodic functions of f . We obtain where the first term of the right hand member is the same as Brouwer's first order determining function of the short-period elimination, and C 1 is an integration constant that is left undetermined by the time being.
On account of H 2,0 ≡ 0, the known terms at the second order of the Lie transforms approach are H 0,2 = {H 1,0 ; W 1 } + {H 0,1 ; W 1 }, from which the terms pertaining to the kernel of the Lie operator must be cancelled by the adequate selection of H 0,2 . The usual choice is H 0,2 = H 0,2 , yet additional terms could be left in the new Hamiltonian in particular cases [10,29]. However, this process would leave purely long-period terms in the new Hamiltonian in addition to the secular terms, both certainly pertaining to the kernel of the Lie operator. Since this is against the total normalization criterion, purely long-period terms should vanish identically in H 0,2 , a requirement that is achieved with the proper selection of C 1 , whose partial derivatives with respect to g, G, and L, appear formally in H 0,2 .
We attack the computation of the second-order of the perturbation theory by parts. To that effect, we make Straightforward evaluation of the Poisson brackets, followed by the use of Eq. (3) and standard trigonometric reduction, yields where the needed coefficients B i,j,k (s) are listed in Table 1. Table 1: Inclination polynomials B i,j,k in Eq. (7).
We intentionally split H 0,2 into three different blocks. Namely, all the terms on the first row of Eq. (7) are free from the equation of the center and factored by a 2 /r 2 , hence being of trivial integration in the true anomaly. Terms of the second row are free from both φ and r; the integration of terms of this type reduces to the well-known case of the integration of cosine functions in elliptic motion [7,22,40]. Terms on the third row are of the form (p/r) 2 φ sin(mf + α), with m integer, and are readily integrated by parts. That is, on account of in this way leading to the previous case of integration of cosine functions. Particularization for definite integration follows from the fundamental theorem of calculus.
It is worth noting that if, on the contrary, we replace r by the conic equation in the third row of Eq. (7), in order to arrange this row like the Fourier series with q 0 = 3e 2 +2, q 1 = 3 4 (e 2 +4), q 2 = 3 2 , and q 3 = 1 4 , then the equation of the center shows as an isolated function of the mean anomaly when the summation index takes the value j = 0. This arrangement brings no problem in the computation of definite integrals, which can be carried out using the general rules for computing φ sin mf and φ cos mf provided in [32]. On the contrary, while indefinite integration is still possible, it requires the sophisticated use of special functions, which could make notably difficult to progress in the perturbation approach [36].
On the other hand, the evaluation of the Poisson brackets involving the integration constant C 1 yields where j = (|j − 2i| − 2)j , (j + 1) ≡ (j + 1) mod 2, and the inclination polynomials b i,j,k (s) and the eccentricity polynomials q i,j (e) are provided in Table 2. Table 2: Non-null inclination b i,j,k and eccentricity polynomials q i,j in Eq. (9).
In the same way as we did in the first order, we choose H 0,2 = H 0,2 = H 0,2 + H * 0,2 to guarantee that it cancels all the terms of H 0,2 perteining to the kernel of the Lie derivative. Firstly, we compute H 0,2 as follows. To average the first row of Eq. (7) over the mean anomaly, it is first multiplied by the factor r 2 /(a 2 η) to carry out the integration in the true anomaly, and then those terms that are free from f , which are those with j = 0, are selected. The term free from f in the second row averages to itself while the remaining terms in this row are averaged using the rule 1 2π cf. [22]. Finally, the terms on the third row of Eq. (7) are averaged by parts with the help of Eqs. (8) and (10). We finally obtain, which is precisely Brouwer's second-order Hamiltonian after the elimination of short-period terms. The average of Eq. (9) is readily obtained with analogous procedures, to obtain Visual inspection of Eqs. (11) and (12), immediately shows that if we complete the computation of the first-order term of the generating function in Eq. (6) choosing then Eq. (12) turns into the opposite of the term in the last row of Eq. (11), the only one that depends on g, thus mutually canceling out. Hence which is the same as the second-order term of Brouwer's Hamiltonian after the elimination of long-period terms. In this way we have achieved Brouwer's total Hamiltonian reduction of the main problem at once, with a single canonical transformation.
It is not a big surprise that C 1 is the same integration constant obtained in Alfriend and Coffee's elimination of the perigee [3,28] or in the author's reverse normalization of the angular momentum [30], since the motion in the orbital plane is decoupled from the rotation of that plane in each case.
The computation of first-order periodic corrections is now straightforward from the simple evaluation of Poisson brackets, namely ξ − ξ = J 2 ∆ξ, where ∆ξ ≡ {ξ, W 1 } and ξ denotes either a canonical variable or some wanted function of the canonical variables [9]. For instance, for the first-order periodic corrections to the semi-major axis we obtain , and the coefficients B i are the same as those in Eq. (4). Recall that Eq. (15) must be evaluated in mean (prime) variables in the direct transformation from mean elements to osculating ones, and in original (unprimed) variables in the inverse transformation from osculating to mean variables.

Second-order periodic corrections
The second-order term of the generating function is now computed making m = 2 in Eq. (2). Namely The needed integrals in the computation of V 2 are either trivial, solved with the help of Eq. (8) for those terms involving the equation of the center, or using the differential relation between the mean and true anomalies for those other that are free from φ but depend on trigonometric functions of f . Straightforward computations yield  where j min = 2(i + 1) − 1, j max = 4 + i + 1 2 (i − 1) , and the inclination polynomials β i,j,k are listed in Table 3.
As before, the integration constant C 2 will be determined by imposing to the known terms of the next order where H 1,1 = H 0,2 + {H 0,1 , W 1 }, the condition of being free from pure longperiodic terms. Again, the known terms are split into terms directly computable and those depending on the arbitrary function C 2 . That is, It follows the customary computation of H 0,3 so that it cancels the terms of H 0,3 pertaining to the kernel of the Lie operator; namely After straightforward evaluation of the Poisson brackets, we obtain where i = i mod 2 and the inclination polynomials β i,k are given in Table 4. Analogously, In this process, we only found integrals of the same type as we did at the second order, and hence there were no special difficulties in solving them, yet in this case we needed to deal with notably longer series than in previous orders.  Once more, the simple inspection of Eqs. (17) and (18) shows that if we now choose then Eq. (18) cancels the terms of Eq. (17) depending on the argument of the perigee out, to yield which is completely reduced as desired.
Beyond the first order, direct and inverse transformations are no longer opposite. At the second order the direct transformation is given by ξ = ξ + J 2 ∆ξ + where the inclination coefficients A i,j,k are provided in Table 21.

Initialization of the secular constants and performance tests
Soon after Brouwer's solution appeared in print, different reports pointed out an apparent contradiction between the accuracy expected from the series truncation order and the comparatively large in-track errors obtained in a variety of tests against numerical integrations [5]. The issue, however, did not happen when fitting Brouwer's solution to observational data. Hence, the apparent discrepancy was easily identified with an inconsistency in the use of Brouwer's theory. Indeed, to get the expected accuracy provided by the secular terms, the initialization of the constants of Brouwer's solution should be done with analogous accuracy. However, Brouwer only provided the periodic corrections up to the first order of J 2 , and hence the direct initialization of the secular mean motion for given initial conditions yields analogous accuracy. The trouble is, of course, solved if the inverse periodic corrections are computed up to the same order as the secular terms.
On the other hand, since the trouble arises from an inaccurate computation of the secular mean motion, the theory can be patched by supplementing Brouwer's first order corrections only with the inverse second-order correction to the semimajor axis, either using Eq. (21) or in the different much shorter formulation given by [31]. Alternatively, the errors in the initialization procedure can be palliated by fitting the secular frequencies to data obtained from a preliminary numerical integration over several revolutions, or by a calibration of the secular mean motion n = µ 2 /L 3 from the energy equation [6].
The latter approach is particularly appealing because it totally avoids the need of carrying out additional computations to those already carried out by Brouwer. Thus, for given initial conditions ( 0 , g 0 , h 0 , L 0 , G 0 , H 0 ), the initial Hamiltonian in osculating elements evaluates to H( 0 , g 0 , −, L 0 , G 0 , H 0 ) = E 0 . On the other hand, after the complete Hamiltonian reduction However, the constants L and G are computed from the osculating initial conditions through the inverse periodic corrections only up to O(J k−1

2
). While this fact does not compromise the accuracy of Eq. (22) in what respects to the terms H 0,m (m ≥ 1) because they are multiplied by corresponding factors J m 2 , it certainly does in the case of the Keplerian term. What Breakwell and Vagners [6] propose is then to replace L by the calibrated value obtained by solving the Keplerian term from Eq. (22). If now L is replaced in Eq. (22) by the calibrated value L then the energy equation will remain certainly accurate to O(J k+1 2 ). Therefore, the initialization of the secular frequencies is notably improved using the values Obviously, the use of Breakwell and Vagners' calibration procedure is not constrained to the second-order of Brouwer's theory, and also applies to any truncation order. In our particular case, in which we had already computed the second-order direct and inverse periodic corrections of Brouwer's theory, the calibration of the (mean) Delaunay action allowed us to improve the accuracy of Brouwer's secular terms to the third order of J 2 without need of computing the long series that comprise the third-order term of the generating function.
We checked that the new extended third-order Brouwer's solution, in which the second-order corrections consist of a single transformation, enjoys the same accuracy as a third-order solution computed in the traditional way of splitting the Hamiltonian reduction into the sequence provided by the elimination of the parallax, followed by the elimination of the perigee, and ending with a Delaunay normalization [3]. At the third order, the initialization of the constants of the later was also calibrated in Breakwell and Vagners' style. In both cases, the required direct and inverse transformations were computed in polar variables to avoid singularities in the case of circular orbits, but also for efficiency reasons [20,2,24,25].
An example of the accuracy obtained with the new single-transformation approach is presented in Fig. 1  For efficiency in the evaluation of perturbation solutions, arrangement of the series that comprise the periodic corrections for optimal evaluation is an important consideration [8,17]. In this task, we limited our efforts to minor arrangements of the code, like the factorization of the inclination polynomials involved in the different summations and the following use of Horner's algorithm, and left the code optimization job to the compiler. Because we did the same for both analytical solutions (Brouwer's with single periodic corrections, and the traditional parallax-perigee-Delaunay solution), even if optimal evaluation is not achieved, the comparisons are not expected to be biased towards a particular theory.
After repeated evaluation of the periodic corrections for a variety of initial conditions, we found that the evaluation of the periodic corrections of the traditional analytical solution spends roughly twice the time needed by the singletransformation approach in the evaluation of the periodic corrections. This result was a priori unexpected because the series comprising the corrections of the new approach, which only involves a single transformation, are clearly longer than the composition of those involved in the three transformations needed in the traditional approach. The improved evaluation efficiency is then attributed to the fact that the compiler is able to carry out a better optimization of the code in the case of the single-transformation approach. This fact may be understood when taking into account that, for instance, the coefficient (5s 2 − 4) appears about 300 times in our arrangement of the periodic corrections of the single transformation, but only 73 times in the classical parallax-perigee-Delaunay transformation arranged with the same factorization criterion, where, in particular, this factor only appears in the corrections related to the elimination of the perigee. Thus, cancelling this common factor by the compiler is roughly four times more efficient in the first case. Undeniably, making a smarter organization of the code before sending it to the compiler might modify these figures. However, the balance is so radical on the side of the single transformation that these presumed improvements due to an additional preprocessing of the code are not expected to be relevant enough to revert the figures.

Conclusions
Experience gained through the use of Hamiltonian simplification methods prompted us to question Brouwer's splitting normalization strategy. The convenience of dividing a normalization process into different stages has been taken for granted since the initial efforts in fully automatizing the computation of perturbation theories. Needless to say that we agree in which this way of proceeding may ease the construction of the perturbation solution. However, what is not so obvious is that the evaluation of the solution constructed this way must necessarily yield the less computational burden. On the contrary, results in this paper seem to point in the direction that the claimed benefits of partial normalization as well as Hamiltonian simplification procedures can be counterbalanced by other type of considerations, at least for the lower orders of normalization that suffice in many practical cases. Prospective application of the strategy proposed here to other instances of perturbed Keplerian motion, or to the computation of higher orders of the main problem of artificial satellite theory, should contribute to make clear the issue.
Brouwer's closed-form approach and full automatization of the computation of perturbation theories seem two legitimate aims in this epoch of computational plenty. However, as demonstrated by the equation of the center controversy, rather than running perturbation algorithms in batch processes, one should not disregard the power of modern hand computations carried out with the help of existing software tools. Indeed, as far as mathematical simplification remains in the category of arts, inspection of intermediate expressions turns into a convenient practice that may eventually lead the user to straightforward simplifications that make feasible or just simpler the next step of a partially automated procedure. Like chess players, celestial mechanicians are rarely able to anticipate more than a few moves in the outcome of a perturbation approach. On the contrary, they need to wait for the opponent's reaction in order to implement a winning strategy, which, in addition, is most times settled on an empirical basis. It was, in particular, the case of the current research, in which the help provided by the computer algebra system converted into a simple task the critical inspection of the seminal solution obtained by Brouwer.