#  A torsion-based solution to the hyperbolic regime of the $J_{2}$-problem
††thanks: Preprint submitted to _Nonlinear Dynamics_ July 23, 2022

Martin Lara GRUCACI, University of La Rioja, and Space Dynamics Group –
<EMAIL_ADDRESS>Alessandro Masat Department of Aerospace Science and
Technology, Politecnico di Milano, Milano, ItalyPhD Candidate,
<EMAIL_ADDRESS>Camilla Colombo Department of Aerospace Science
and Technology, Politecnico di Milano, Milano, ItalyAssociate Professor,
<EMAIL_ADDRESS>
(Manuscript of )

###### Abstract

A popular intermediary in the theory of artificial satellites is obtained
after the elimination of parallactic terms from the $J_{2}$-problem
Hamiltonian. The resulting quasi-Keplerian system is in turn converted into
the Kepler problem by a torsion. When this reduction process is applied to
unbounded orbits the solution is made of Keplerian hyperbolae. For this last
case, we show that the torsion-based solution provides an effective
alternative to the Keplerian approximation customarily used in flyby
computations. Also, we check that the extension of the torsion-based solution
to higher orders of the oblateness coefficient yields the expected convergence
of asymptotic solutions to the true orbit.

### Keywords

Hamiltonian mechanics - canonical perturbation theory - perturbed Keplerian
motion - torsion - unbounded orbits - flyby

##  1 Introduction

The hyperbolic regime of an artificial satellite of an oblate body was solved
in closed form by Hori [13], who showed in this way that the perturbation
approach may be useful without restriction to the typical case of quasi-
periodic motion.111As in the case of bounded motion [15], the invariant
manifold of equatorial orbits of the main problem is integrable for unbounded
motion too [24]; see [3] for a general discussion. However, at difference from
the classical elliptic regime, perturbed hyperbolic motion is subjected to
strict boundary conditions that must be imposed to the analytical approach.
More precisely, the effects of the oblateness perturbation vanish when the
artificial satellite approaches infinity in its hyperbolic-type motion. In
consequence, the transformation between mean and osculating variables that
encapsulates the perturbation solution must become the identity when evaluated
at infinity.

In Hori’s first-order approach the $J_{2}$-problem Hamiltonian is directly
reduced to the Kepler Hamiltonian in (mean) hyperbolic Delaunay variables.
With this choice, the solution of the mean-to-osculating transformation that
yields the complete Hamiltonian reduction depends unavoidably on non-periodic
terms of the true anomaly. In spite of that, the explicit appearance of the
true anomaly in Hori’s solution is not of concern since its value is bounded
by the asymptotes of a Keplerian hyperbola. Therefore, these non-periodic
terms do not play any secular role. After application of the boundary
conditions to Hori’s transformation these kinds of terms appear always like
coefficients involving the difference between the instantaneous true anomaly
and the constant true anomaly of the asymptote of the Keplerian hyperbola.
While non-periodic, this type of coefficient resembles the equation of the
center, given by the difference between the true and mean anomalies, which
comprises essential short-period terms of the elliptic motion. But the
resemblance goes beyond a formal likeness, and we find that these terms can be
handled advantageously using known techniques for efficiently dealing with the
equation of the center.

In particular, it is known that the equation of the center does not play any
role in the elimination of the parallax simplification, which, on the
contrary, depends only on trigonometric functions of the true anomaly, cf.
§6.2.1 of [19]. That this last feature is not exclusive of the case of bounded
perturbed Keplerian motion and also applies to unbounded motion becomes
evident from Deprit’s original formulation of the elimination of the parallax
in polar variables [6]. Furthermore, Deprit’s radial intermediary, which is
obtained after the elimination of the parallax up to first order effects of
the oblateness coefficient, provides an integrable solution of the
$J_{2}$-problem, also called the _main problem_ in the theory of artificial
satellites [1], that is valid for both bounded and unbounded orbits.
Nevertheless, the nature of the solution is different in the two distinct
regimes, and Deprit’s reduction of the main problem to a quasi-Keplerian
system must be amended in the second case in order to take the boundary
conditions that characterize hyperbolic-type motion into account.

We compute the additional terms that are needed in the elimination of the
parallax transformation for converting the hyperbolic regime of the main
problem into a quasi-Keplerian system, which, in turn, is converted into a
pure Keplerian system by the standard torsion transformation [6]. When the new
transformation is applied to convert the mean hyperbolic Delaunay variables
into osculating ones, the transformation equations are simpler than those
provided by Hori. Still, the new transformation shares analogous deficiencies
to those of Hori’s solution. Namely, it fails in the limit case of parabolic
orbits, and the accuracy of the perturbation solution deteriorates for values
of the eccentricity close to 1. In this regards, it is worth mentioning that
existing perturbation solutions of bounded motion suffer analogous limitations
too.

Moreover, we show that the reduction to a quasi-Keplerian hyberbola is easily
extended to higher orders, yet the appearance of non-periodic functions of the
true anomaly in the transformation equations is now unavoidable. While these
higher order refinements of the perturbation solution may not have relevant
effects in practice —for flyby design or gravity assist maneuvers— they
radically reduce the errors of the perturbation solution close to, and after
perigee passage, and serve us to illustrate how the perturbation approach
follows the expected convergence beyond its customary application to the case
of quasi-periodic motion. Needless to say that analogous improvements are
expected if Hori’s perturbation solution were eventually extended beyond the
first order of the oblateness coefficient explicitly computed by himself.

Due to the hyperbolic character of unbounded perturbed Keplerian motion, one
might expect high sensitivity of the main problem dynamics with respect to
initial conditions. However, for common values of the $J_{2}$ coefficient of
solar system bodies the perturbation solution succeeds in capturing these
dynamics, and can even predict the qualitative features of branches of the
classical Newtonian rosetta when the energy becomes positive.222In the case of
zero energy, these branches are sometimes called “fish” orbits [24] In this
extreme case, high orders of the perturbation solution are indeed required.

##  2 Deprit’s radial intermediary

The main problem of the artificial satellite admits Hamiltonian formulation.
When using polar canonical variables, it is written in the form

$\mathcal{M}=\frac{1}{2}\left(R^{2}+\frac{\Theta^{2}}{r^{2}}\right)-\frac{\mu}{r}-\frac{1}{4}J_{2}\frac{\mu}{r}\frac{\alpha^{2}}{r^{2}}\left(2-3\sin^{2}I+3\sin^{2}I\cos
2\theta\right),$ (1)

where the physical parameters $\mu$, $\alpha$, and $J_{2}$ characterize the
gravitational field of the attracting body and denote the gravitational
parameter, equatorial radius, and oblateness coefficient, respectively. The
polar canonical coordinates $(r,\theta)$ stand for the radius and the polar
angle, respectively, and their conjugate momenta $(R,\Theta)$ for the radial
velocity and the specific angular momentum, respectively;
$I=\arccos(N/\Theta)$ is the inclination of the instantaneous orbital plane,
and the canonical pair $(\nu,N)$ represents the right ascension of the
ascending node and the third component of the angular momentum vector per mass
unit, respectively. Because $\nu$ is an ignorable coordinate of Hamiltonian
(1), we immediately infer that $N$ is an integral of the main problem
dynamics.

Up to first order effects of $J_{2}$, Hamiltonian (1) can be reduced to the
integrable quasi-Keplerian Hamiltonian in new, prime variables

$\mathcal{M}=\frac{1}{2}\left(R^{\prime 2}+\frac{\Theta^{\prime 2}}{r^{\prime
2}}\right)-\frac{\mu}{r^{\prime}}-\frac{1}{2}\frac{\Theta^{\prime
2}}{r^{\prime 2}}\frac{1}{2}J_{2}\frac{\alpha^{2}}{(\Theta^{\prime
2}/\mu)^{2}}\left[2-3\sin^{2}I(\Theta^{\prime},N)\right],$ (2)

which is known as Deprit’s radial intermediary [6]. In the case of bounded
motion this intermediary has been proposed for different applications [2, 10,
11]. The reduction of the main problem to Deprit’s radial intermediary is
achieved by an infinitesimal contact transformation
$(r,\theta,\nu,R,\Theta,N,J_{2})\mapsto(r^{\prime},\theta^{\prime},\nu^{\prime},R^{\prime},\Theta^{\prime},N)$
of generating function

$\mathcal{U}_{1}=-\Theta\frac{1}{8}\frac{\alpha^{2}}{p^{2}}\left[\left(4\kappa+3\right)s^{2}\sin
2\theta+(4-6s^{2}-2s^{2}\cos 2\theta)\sigma\right]+\mathcal{C}_{0},$ (3)

where $s$ abbreviates the sine of the inclination, $p=\Theta^{2}/\mu$ is the
parameter of the conic, and the non-dimensional functions

$\sigma=\frac{pR}{\Theta},\qquad\kappa=\frac{p}{r}-1,$ (4)

represent the projections of the eccentricity vector in the orbital frame.
Namely, $\sigma=e\sin{f}$ and $\kappa=e\cos{f}$, with $e$ denoting orbital
eccentricity and $f$ is the true anomaly, from which we readily obtain
$\lim_{r\to\infty}\cos{f}=-1/e$.

Differences between Eq. (3) and the equivalent term in [6] are a consequence
of a distinct arrangement of the involved functions, on the one hand, and the
addition of the “integration constant” $\mathcal{C}_{0}$, on the other. The
latter must fulfill the condition of being free from the explicit appearance
of $\theta$ when it is formulated in the specific algebra of the functions

$C=\frac{\Theta}{p}\left(\kappa\cos\theta+\sigma\sin\theta\right),\qquad
S=\frac{\Theta}{p}\left(\kappa\sin\theta-\sigma\cos\theta\right),$ (5)

used in [6]. However, their use may obscure the real nature of this
integration constant, which is no other than being free from short-period
effects [22]. For this reason, but also for simplicity, carrying out the
elimination of the parallax in Delaunay variables, rather than in the set of
polar variables, is nowadays encouraged [21].

Delaunay variables are traditionally denoted $(\ell,g,h,L,G,H)$ and represent
the mean anomaly, the argument of the periapsis, the longitude of the
ascending node, the so-called Delaunay action, the specific angular momentum,
and the third component of the angular momentum vector per mass unit,
respectively. When the generating function $\mathcal{U}_{1}$ in Eq. (3) is
reformulated in Delaunay variables, one must distinguish between the cases of
bounded and unbounded motion. For the former the classical set of elliptic
Delaunay variables is used [4, 7, 8], whereas hyperbolic Delaunay variables
are used for the latter instead [13, 9]. While in both cases Eq. (3) takes the
form

$\displaystyle\mathcal{U}_{1}=$
$\displaystyle-G\frac{1}{8}\frac{\alpha^{2}}{p^{2}}\left\\{s^{2}[3e\sin(f+2g)+3\sin(2f+2g)+e\sin(3f+2g)]\right.$
$\displaystyle\left.-(6s^{2}-4)e\sin{f}\right\\}+\mathcal{C}_{0}(g,L,G,H),$
(6)

the meaning of the symbols used is slightly different. Thus, while
$p=G^{2}/\mu$ and $I=\arccos{H}/{G}$ have the same meaning in both regimes of
motion, the eccentricity is $e=(1-\eta^{2})^{1/2}<1$, with $\eta=G/L$, for
bounded motion, whereas $e=(1+\eta^{2})^{1/2}>1$ and $\eta=-G/L$ in the
hyperbolic case. The true anomaly $f$ is an implicit function of the mean
anomaly and the eccentricity, whose computation requires the preliminary
solution of Kepler’s equation. For unbounded motion, Kepler equation takes the
form

$\ell=e\sinh{u}-u,$ (7)

where the hyperbolic anomaly $u$ is related to the true anomaly by means of

$\sqrt{e+1}\tanh\left(\mbox{$\frac{1}{2}$}u\right)=\sqrt{e-1}\tan\left(\mbox{$\frac{1}{2}$}f\right).$
(8)

##  3 Polar variables mapping to the quasi-Keplerian system

In his perturbation approach Hori chose that the new Hamiltonian be just the
Keplerian part of the perturbation problem. That is, the transformed
Hamiltonian in the prime hyperbolic Delaunay variables reads
$\mathcal{M}^{\prime}=+\mu/(2a)$, where $a={L^{\prime}}^{2}/\mu$ is the semi-
transverse axis of the osculating hyperbola. Because of this choice, non-
trigonometric terms of the true anomaly show up in the generating function,
which Hori solved up to an integration constant, as we did in Eq. (6). At
difference from the elliptic case the latter is in no way arbitrary, and Hori
determined it by imposing the required boundary conditions at infinity to his
transformation of the hyperbolic Delaunay variables.333Note a typo in Hori’s
final transformation equation for the longitude of the node, where the
exponent of $G^{\prime}$ in the first coefficient must be 4 instead of 3, as
readily shown by checking dimensions.

Hori had not yet invented his perturbation method based on Lie series [14] and
relied on von Zeipel’s algorithm for computing his first order solution [25].
Rather, we resort to the method of Lie transforms [5] in following
derivations. Moreover, in contrast with Hori’s choice of the Keplerian, we are
satisfied with finding the transformation that yields Deprit’s radial
intermediary, in this way simplifying the computation of the mapping.

At the end, the quasi-Keplerian system (2) is converted into a pure Keplerian
system by means of a torsion transformation [6]. Since the torsion is carried
out in polar variables, we find advantageous to derive the transformation
equations of the polar variables, rather than the hyperbolic Delaunay ones,
directly from Eq. (6). That is, to the first order

$\xi=\xi^{\prime}+J_{2}\xi_{0,1}|_{\xi=\xi^{\prime}}$ (9)

where $\xi\equiv(r,\theta,\nu,R,\Theta,N)^{\tau}$, $\tau$ means transposition,
$\xi_{0,1}\equiv\\{\xi,\mathcal{U}_{1}\\}$, and the curly brackets denote the
Poisson bracket operator. Still, since the boundary conditions are naturally
formulated in terms of the true anomaly, we find convenience in writing the
transformation in terms of the hyperbolic Delaunay variables. To this end, we
need to avail ourselves with the partial derivatives of the polar variables
with respect to the hyperbolic Delaunay ones. They are listed in Table 1. Note
that, while mostly analogous to the case of elliptic Delaunay variables, sign
changes appear in the partial derivatives with respect to the hyperbolic
Delaunay action, cf. Table 5.1 of [19].

| $\partial/\partial\ell$ | $\partial/\partial{L}$ | $\partial/\partial{G}$
---|---|---|---
$\theta$ | $\displaystyle\frac{p^{2}}{r^{2}}\frac{1}{\eta^{3}}$ | $\displaystyle\eta\frac{\partial\theta}{\partial{G}}$ | $\displaystyle-\frac{1}{e^{2}}\frac{pR}{\Theta^{2}}\left(\frac{p}{r}+1\right)$
$r$ | $\displaystyle\frac{p^{2}R}{\eta^{3}\Theta}$ | $\displaystyle\eta\left(\frac{\partial{r}}{\partial{G}}-2\frac{r}{\Theta}\right)$ | $\displaystyle\frac{1}{e^{2}}\frac{p}{\Theta}\left(\frac{p}{r}-1\right)$
$R$ | $\displaystyle\frac{e^{2}}{\eta^{3}}\frac{\Theta^{2}}{r^{2}}\frac{\partial{r}}{\partial{G}}$ | $\displaystyle\eta\left(\frac{\partial{R}}{\partial{G}}+\frac{R}{\Theta}\right)$ | $\displaystyle-\frac{1}{e^{2}}\frac{p^{2}}{r^{2}}\frac{R}{\Theta}$
Table 1: Non-vanishing partial derivatives of the polar variables with respect
to the hyperbolic Delaunay variables ($\partial\Theta/\partial{G}=1$).

We readily find that the correction to the radial velocity automatically
fulfills the boundary condition $\\{R,\mathcal{U}_{1}\\}|_{\infty}\equiv 0$.
Next, we compute

$\\{\Theta,\mathcal{U}_{1}\\}|_{\infty}=G\frac{1}{4}\frac{\alpha^{2}}{p^{2}}\frac{s^{2}}{e^{2}}\left[2\eta^{3}\sin
2g-(3\eta^{2}+1)\cos 2g\right]-\frac{\partial{\mathcal{C}_{0}}}{\partial{g}},$

from which $\mathcal{C}_{0}$ is solved to make the right hand of the equation
to vanish. We obtain

$\mathcal{C}_{0}=-G\frac{1}{8}\frac{\alpha^{2}}{p^{2}}\frac{s^{2}}{e^{2}}\left[2\eta^{3}\cos
2g+(3e^{2}-2)\sin 2g\right]+\mathcal{C}_{1}(L,G,H),$

which is the same as the function $Z$ computed by Hori except for a sign due
to the different convention in the Hamiltonian formulation used by Hori. It
follows the computation of

$\\{\nu,\mathcal{U}_{1}\\}|_{\infty}=\frac{3}{2}c\frac{\alpha^{2}}{p^{2}}\eta+\frac{\partial{\mathcal{C}_{1}}}{\partial{H}},$

where $c=H/G$. The constant $\mathcal{C}_{1}$ is then solved after imposing
the vanishing of the right side. We obtain

$\mathcal{C}_{1}=-G\frac{3}{4}\frac{\alpha^{2}}{p^{2}}\eta(1-s^{2})+\mathcal{C}_{2}(L,G).$

Analogously,

$\\{\theta,\mathcal{U}_{1}\\}|_{\infty}=\frac{1}{2}\eta\frac{\alpha^{2}}{p^{2}}+\frac{\partial{\mathcal{C}_{2}}}{\partial{G}},$

from which, recalling that $\eta=-G/L$ and $p=G^{2}/\mu$, we obtain

$\mathcal{C}_{2}=G\frac{1}{4}\frac{\alpha^{2}}{p^{2}}\eta+\mathcal{C}_{3}(L).$

Finally, we check that the choice $\mathcal{C}_{3}=0$ fulfills the last
condition $\\{r,\mathcal{U}_{1}\\}|_{\infty}=0$. In summary,

$\mathcal{C}_{0}=G\frac{1}{4}\frac{\alpha^{2}}{p^{2}}\left\\{(3s^{2}-2)\eta-\frac{s^{2}}{e^{2}}\left[\eta^{3}\cos
2g+\frac{1}{2}(3e^{2}-2)\sin 2g\right]\right\\},$ (10)

thus making Eq. (6) fully determined.

It is worth noting that Eq. (10) can be rewritten in polar variables using the
alternative form

$\mathcal{C}_{0}=\Theta\frac{1}{4}\frac{\alpha^{2}}{p^{2}}\left\\{\left(3s^{2}-2\right)\eta-\frac{p^{2}s^{2}}{e^{4}\Theta^{2}}\left[\eta^{3}\left(C^{2}-S^{2}\right)+\left(3e^{2}-2\right)CS\right]\right\\},$

with $e=(p/\Theta)\sqrt{C^{2}+S^{2}}$, whereas $C\equiv{C}(r,\theta,R,\Theta)$
and $S\equiv{S}(r,\theta,R,\Theta)$ are given in Eq. (5). In this way it is
shown that $\mathcal{C}_{0}$ is free from the _explicit_ appearance of
$\theta$, like it must be when using Deprit’s algebra based on these functions
[6].

Once the generating function has been fully determined, we compute the
corrections

$\displaystyle r_{0,1}=$
$\displaystyle\;p\frac{1}{4}\frac{\alpha^{2}}{p^{2}}\Big{\\{}(3s^{2}-2)\left(1+\frac{e}{\eta}\sin{f}\right)+\frac{s^{2}}{2e^{3}}\left[(e^{2}-4)\eta\sin(f-2g)\right.$
$\displaystyle-3e^{2}\eta\sin(f+2g)+(3e^{2}-4)\cos(f-2g)+3e^{2}\cos(f+2g)$
$\displaystyle\left.+2e^{3}\cos(2f+2g)\right]\Big{\\}}$ (11)
$\displaystyle\theta_{0,1}=$
$\displaystyle\;\frac{1}{16}\frac{\alpha^{2}}{p^{2}}\Big{(}\frac{1}{\eta}\Big{\\{}6\left[2(5s^{2}-4)-(7s^{2}-6)e^{2}\right]+8e(3s^{2}-2)\cos{f}+2e^{2}$
$\displaystyle\times(3s^{2}-2)\cos
2f\Big{\\}}+\frac{\eta}{e^{3}}\Big{\\{}(e^{2}-4)es^{2}\cos(2f-2g)+4(e^{2}-4)s^{2}$
$\displaystyle\times\cos(f-2g)+2e\left[e^{2}(7s^{2}-4)-4(4s^{2}-1)\right]\cos
2g-12e^{2}s^{2}$
$\displaystyle\times\cos(f+2g)-3e^{3}s^{2}\cos(2f+2g)\Big{\\}}+\frac{1}{e^{3}}\Big{\\{}(4-3e^{2})es^{2}\sin(2f-2g)$
$\displaystyle-4(3e^{2}-4)s^{2}\sin(f-2g)+2e\left[3e^{2}(5s^{2}-2)-4(4s^{2}-1)\right]\sin
2g$
$\displaystyle-8e^{4}(6s^{2}-5)\sin{f}+4e^{2}\left[e^{2}(5s^{2}-3)-3s^{2}\right]\sin(f+2g)$
$\displaystyle+e^{3}(11s^{2}-12)\sin(2f+2g)+4e^{4}(s^{2}-1)\sin(3f+2g)\Big{\\}}\Big{)}$
(12) $\displaystyle\nu_{0,1}=$
$\displaystyle\;c\frac{1}{4}\frac{\alpha^{2}}{p^{2}}\Big{\\{}\frac{1}{e^{2}}\left[(3e^{2}-2)\sin
2g+2\eta^{3}\cos 2g\right]-6(\eta+e\sin{f})$
$\displaystyle+3e\sin(f+2g)+3\sin(2f+2g)+e\sin(3f+2g)\Big{\\}}$ (13)
$\displaystyle R_{0,1}=$
$\displaystyle\;\frac{G}{p}\frac{1}{32}\frac{\alpha^{2}}{p^{2}}\Big{\\{}\frac{e}{\eta}(3s^{2}-2)\left[2e^{2}\cos
3f+8e\cos 2f+(6e^{2}+8)\cos{f}+8e\right]$
$\displaystyle+\eta\frac{s^{2}}{e^{3}}\left[(e^{2}-4)e^{2}\cos(3f-2g)+4(e^{2}-4)e\cos(2f-2g)-(e^{4}+4e^{2}\right.$
$\displaystyle+16)\cos(f-2g)-8(e^{2}+2)e\cos 2g-(5e^{2}+16)e^{2}\cos(f+2g)-12$
$\displaystyle\left.\times
e^{3}\cos(2f+2g)-3e^{4}\cos(3f+2g)\right]-\frac{s^{2}}{e^{3}}\left[(3e^{2}-4)e^{2}\sin(3f-2g)\right.$
$\displaystyle+4(3e^{2}-4)e\sin(2f-2g)+(3e^{4}+4e^{2}-16)\sin(f-2g)+4(e^{4}+4)$
$\displaystyle\times e\sin
2g+(19e^{2}+16)e^{2}\sin(f+2g)+4(2e^{2}+7)e^{3}\sin(2f+2g)$
$\displaystyle\left.+19e^{4}\sin(3f+2g)+4e^{5}\sin(4f+2g)\right]\Big{\\}}$
(14) $\displaystyle\Theta_{0,1}=$
$\displaystyle\;\Theta\frac{1}{4}\frac{\alpha^{2}}{p^{2}}s^{2}\Big{\\{}\frac{1}{e^{2}}\left[(3e^{2}-2)\cos
2g-2\eta^{3}\sin 2g\right]$
$\displaystyle+3e\cos(f+2g)+3\cos(2f+2g)+e\cos(3f+2g)\Big{\\}}.$ (15)

The appearance of $\eta$ in denominators of $r_{0,1}$, $\theta_{0,1}$ and
$R_{0,1}$ will deteriorate the efficiency of these corrections for orbits
close to parabolic. Note that the dependence of Eqs. (11)–(15) on the true
anomaly is only through trigonometric functions, in clear contrast with Hori’s
corrections [13]. On the other hand, the explicit appearance of $f$ in Hori’s
solution is not criticizable as far as the possible values of the true anomaly
are constrained by the values taken by the asymptotes of the Keplerian
hyperbola.

The analytical solution is completed by applying a torsion to the quasi-
Keplerian Hamiltonian (2). Disregarding the prime notation without risk of
confusion, Eq. (2) is rearranged like

$\mathcal{M}=\frac{1}{2}\left(R^{2}+\frac{\tilde{\Theta}^{2}}{r^{2}}\right)-\frac{\mu}{r},$
(16)

where

$\tilde{\Theta}^{2}=\Theta^{2}\left[1-\frac{1}{2}J_{2}\frac{\alpha^{2}}{p^{2}}(3c^{2}-1)\right]$
(17)

is function only of $\Theta$ and the parameters of the simplified problem;
namely the integral $N$ and the physical parameters defining the gravitational
field. The objective of the torsion is to find a canonical change of variables
$(\theta,\nu,\Theta,N)\mapsto(\theta^{*},\nu^{*},\Theta^{*},N^{*})$ such that
in the new variables $\tilde{\Theta}=\Theta^{*}$. That is, in the asterisk
variables Eq. (16) turns into the Kepler problem Hamiltonian.

Working in polar variables there is no distinction between bounded or
unbounded motion, so the torsion transformation is exactly the same as in the
elliptic case originally described in §9 of [6]. Rather, we follow the
notation in [11] and denote $\Phi=\Phi(\Theta,N)\equiv\tilde{\Theta}/\Theta$,
$\epsilon=\epsilon(\Theta)\equiv-\frac{1}{2}J_{2}(\alpha/p)^{2}<0$, and define
the torsion by the sequence,

$\displaystyle\Theta^{*}=$ $\displaystyle\;\Theta\Phi,$ (18)
$\displaystyle\theta^{*}=$
$\displaystyle\;\theta\Phi\left(\Phi^{2}-2\epsilon\frac{\partial\Phi^{2}}{\partial\epsilon}-\frac{1}{2}c\frac{\partial\Phi^{2}}{\partial{c}}\right)^{-1},$
(19) $\displaystyle\nu^{*}=$
$\displaystyle\;\nu-\frac{1}{2}\frac{\theta^{*}}{\Phi}\frac{\partial\Phi^{2}}{\partial{c}},$
(20)

whereas not only $N$, but $r$ and $R$ remain unaltered. The inverse
transformation requires first to solve $\Theta=\Theta(\Theta^{*},N)$ from the
implicit function (18), which is then substituted in Eqs. (19) and (20) to be
trivially solved in $\theta$ and $\nu$. Alternatively to the numerical
solution by root-finding procedures, Eq. (18) can be inverted analytically to
be given in the form of a power series in $\epsilon$. If we limit to the first
order of $J_{2}$ of the perturbation solution computed hitherto, we obtain

$\Theta=\Theta^{*}\left[1-\frac{1}{2}(3c^{2}-1)\epsilon\right],$ (21)

in which, now, $c=N/\Theta^{*}$, and $p=\Theta^{*2}/\mu$ must be replaced in
the definition of $\epsilon$.

After written in asterisk variables, Eq. (16) becomes a Keplerian that is next
completely reduced by the exact transformation from polar to hyperbolic
Delaunay variables. Thus, $\mathcal{M}=+\mu^{2}/(2L^{2})$, from which

$\ell=\ell_{0}+nt,\quad g=g_{0},\quad h=h_{0},\quad L=L_{0},\quad
G=G_{0},\quad H=H_{0},$

with $n=-\mu^{2}/L^{3}$.

Note that the Keplerian orbit obtained with the torsion may have no physical
sense, showing complex values of the inclination. This fact is clearly
observed in the limit case obtained when making $\Theta=N$ in Eq. (21).
Namely,

$\frac{N}{\Theta^{*}}\equiv{c}=1-(3c^{2}-1)\epsilon,$

a quadratic equation in $c$ that exists only in the interval
$-\frac{1}{2}+\frac{1}{\sqrt{6}}\leq\epsilon<0$, in which $c>1$ always.
Irrespective of the meaning that $c$ may have in Eq. (21) in these particular
cases, the torsion transformation still remains valid.

##  4 Higher-order extensions

In the case of bounded orbits, a second order solution of the main problem
based on a quasi-Keplerian Hamiltonian would be unavoidably harmed by the
appearance of secular terms in the transformation equations. However, these
kinds of terms are not at all of worry in the case of unbounded motion, in
which their range of variation is constrained by the limit values of the
asymptotes of a hyperbola. Therefore, we can safely extend the computed first
order solution to higher orders following exactly the same approach. Namely,
first computing the mapping that yields a quasi-Keplerian Hamiltonian while
accepting the boundary conditions at infinity, and then applying a torsion to
obtain the final Keplerian.

Thus, from the known terms obtained at the second order of the elimination of
the parallax, we select those terms that are free from the true anomaly but
also from the argument of the periapsis and the eccentricity. The latter, we
recall, is a function of the radius and radial velocity in addition to the
specific angular momentum. We obtain

$\mathcal{M}_{0,2}=-\frac{\Theta^{2}}{r^{2}}\frac{\alpha^{4}}{p^{4}}\frac{1}{16}(21s^{4}-42s^{2}+20),$
(22)

The simple consideration of this additional term in the reduction of the
quasi-Keplerian Hamiltonian
$\mathcal{M}=\mathcal{M}_{0,0}+J_{2}\mathcal{M}_{0,1}+\frac{1}{2}J_{2}^{2}\mathcal{M}_{0,2}$
yields refinements in the computation of the final Keplerian mean motion, as
well as in the solution of the torsion transformation given by Eqs. (18)–(20),
where, now,

$\Phi=\left[1-\epsilon(3c^{2}-1)+\frac{1}{4}\epsilon^{2}(1-21c^{4})\right]^{1/2}.$
(23)

In addition, the analytical inversion of Eq. (18) provided by Eq. (21) must be
replaced by the second order of $J_{2}$ solution given by

$\Theta=\Theta^{*}\left[1-\frac{\epsilon}{2}(3c^{2}-1)-\frac{3\epsilon^{2}}{4}(2c^{2}-1)(5c^{2}-1)\right],$
(24)

where $c$, and $p$ are now functions of the asterisk variables.

Major improvements will be obtained if, besides, we complement the analytical
solution with second order corrections to the elimination of the parallax.
That is, Eq. (9) is replaced by

$\xi=\xi^{\prime}+J_{2}\xi_{0,1}|_{\xi=\xi^{\prime}}+\frac{1}{2}J_{2}^{2}\xi_{0,2}|_{\xi=\xi^{\prime}},$
(25)

where
$\xi_{0,2}\equiv\\{\xi_{0,1},\mathcal{U}_{1}\\}+\\{\xi,\mathcal{U}_{2}\\}$. To
do that, we need to compute the generating function term $\mathcal{U}_{2}$ of
the Lie transformation that yields $\mathcal{M}_{0,2}$. It is obtained in the
usual way up to an arbitrary constant $Q_{0}\equiv{Q}_{0}(g,h,L,G,H)$.

The new corrections need to match the boundary condition at infinity too.
Analogously to the first order, we readily find that $R_{0,2}|_{\infty}\equiv
0$, and hence the boundary conditions are automatically fulfilled by the
corrections to the radial velocity. On the contrary, the equation
$\Theta_{0,2}|_{\infty}=0$ involves the partial derivative
$\partial{Q}_{0}/\partial{g}$, which is solved by indefinite integration to
give

$\displaystyle Q_{0}=$
$\displaystyle\;G\frac{\alpha^{4}}{p^{4}}\frac{1}{256e^{2}}\Big{\\{}-2s^{2}\left[3e^{4}(17s^{2}-18)+8e^{2}(75s^{2}-68)+8(s^{2}+6)\right.$
$\displaystyle\left.+96\eta^{3}(5s^{2}-4)(\pi+\arctan\eta)\right]\sin
2g-3(3e^{4}+6e^{2}-16)s^{4}\sin 4g$
$\displaystyle+\frac{4s^{2}}{\eta^{2}}\left\\{6\eta^{2}\left[e^{4}(15s^{2}-14)+4(3e^{2}-2)(5s^{2}-4)\right](\pi+\arctan\eta)\right.$
$\displaystyle\left.+\eta\left[e^{4}(278-329s^{2})+e^{2}(298s^{2}-284)+4(s^{2}+6)\right]\right\\}\cos
2g$ $\displaystyle+6\eta(e^{2}+8)s^{4}\cos 4g\Big{\\}}+Q_{1}(L,G,H).$

The new integration constant $Q_{1}$ is solved from the equation
$\nu_{0,2}|_{\infty}=0$, which only involves the partial derivative
$\partial{Q}_{1}/\partial{H}$. We obtain

$\displaystyle Q_{1}=$
$\displaystyle\;G\frac{\alpha^{4}}{p^{4}}\frac{3}{128\eta}(1-s^{2})\left[2e^{2}\eta(5s^{2}+13)(\pi+\arctan\eta)\right.$
$\displaystyle\left.-e^{2}(11s^{2}+75)-85s^{2}+107\right]+Q_{2}(L,G).$

Analogously, the integration constant $Q_{2}$is solved from the equation
$\theta_{0,2}|_{\infty}=0$, which only involves the partial derivative
$\partial{Q}_{2}/\partial{G}$. We obtain

$\displaystyle Q_{2}=$
$\displaystyle\;G\frac{\alpha^{4}}{p^{4}}\frac{1}{128\eta}\left[81e^{2}-30e^{2}\eta(\pi+\arctan\eta)-49\right]+Q_{3}(L).$

We take $Q_{3}\equiv 0$ for simplicity, and finally obtain

$\displaystyle\mathcal{U}_{2}=$
$\displaystyle\;G\frac{\alpha^{4}}{p^{4}}\frac{3}{64e^{2}}\left\\{\left[2e^{4}(15s^{2}-14)+8(3e^{2}-2)(5s^{2}-4)\right]s^{2}\cos
2g\right.$ $\displaystyle\left.-16\eta^{3}(5s^{2}-4)s^{2}\sin
2g-e^{4}(5s^{4}+8s^{2}-8)\right\\}\psi$
$\displaystyle+G\frac{\alpha^{4}}{p^{4}}\frac{1}{256e^{3}\eta}\sum_{k=0}^{2}s^{2k}\bigg{\\{}\sum_{j=j_{0}}^{5}\sum_{i=0}^{3}q_{k,i,j}e^{2i+1-(j\bmod
2)}\cos(jf+2kg)$
$\displaystyle+\eta\sum_{j=j_{0}}^{j_{1}}\sum_{i=0}^{2}p_{k,i,j}e^{2i+1-(j\bmod
2)}\sin(jf+2kg)\bigg{\\}},$ (26)

where $j_{0}=-2(k\bmod 2)-1$, $j_{1}=6-2(k\bmod 2)$, $\psi=\pi-f+\arctan\eta$,
and the inclination polynomials $p_{k,i,j}$ and $q_{k,i,j}$ are listed in
Tables 2 and 3. The second order corrections $\xi_{0,2}$ to the polar
variables are then obtained from the straightforward evaluation of the
corresponding Poisson brackets.

Note that the appearance of the function $\psi$, which also exists in Hori’s
first order solution, was only delayed to the second order of our approach.
This postponement resembles analogous effects in the case of bounded motion,
in which the elimination of the parallax simplification delays by one order of
the perturbation approach the need of integration by parts of terms involving
the equation of the center in the computation of the generating function [12,
18].

$i,j$ | $k=0$ | $k=1$ | $k=2$
---|---|---|---
0,-2 | $0$ | $64(3s^{2}-2)$ | $0$
0,-1 | $0$ | $160(3s^{2}-2)$ | $0$
0,0 | $0$ | $16(s^{2}+6)$ | $-48$
0,1 | $0$ | $0$ | $-48$
0,2 | $-24s^{2}(15s^{2}-4)$ | $0$ | $-120$
0,3 | $-112s^{4}$ | $0$ | $0$
0,4 | $-72s^{4}$ | $0$ | $0$
1,-3 | $0$ | $8(3s^{2}-2)$ | $0$
1,-2 | $0$ | $-80(3s^{2}-2)$ | $0$
1,-1 | $-12s^{2}(13s^{2}-4)$ | $-8(45s^{2}-26)$ | $-12$
1,0 | $2(255s^{4}-576s^{2}+272)$ | $8(149s^{2}-142)$ | $42$
1,1 | $-12s^{2}(13s^{2}-4)$ | $96(2s^{2}-1)$ | $-36$
1,2 | $12s^{2}(45s^{2}-16)$ | $-96(6s^{2}-5)$ | $60$
1,3 | $4s^{2}(7s^{2}+8)$ | $0$ | $-120$
1,4 | $90s^{4}$ | $0$ | $-54$
1,5 | $-12s^{4}$ | $0$ | $0$
2,-3 | $0$ | $-10(3s^{2}-2)$ | $0$
2,-2 | $0$ | $16(3s^{2}-2)$ | $0$
2,-1 | $8(75s^{4}-72s^{2}+20)$ | $-2(153s^{2}-134)$ | $15$
2,0 | $6(11s^{4}+64s^{2}-48)$ | $-4(329s^{2}-278)$ | $6$
2,1 | $8(75s^{4}-72s^{2}+20)$ | $-6(127s^{2}-90)$ | $81$
2,2 | $4(27s^{4}-72s^{2}+32)$ | $24(9s^{2}-10)$ | $60$
2,3 | $s^{2}(145s^{2}-64)$ | $-2(195s^{2}-146)$ | $105$
2,4 | $-18s^{4}$ | $-36(3s^{2}-2)$ | $54$
2,5 | $15s^{4}$ | $0$ | $-9$
3,-3 | $0$ | $2(3s^{2}-2)$ | $0$
3,-1 | $-6(5s^{4}+4s^{2}-4)$ | $24(7s^{2}-6)$ | $-3$
3,1 | $-6(5s^{4}+4s^{2}-4)$ | $12(25s^{2}-22)$ | $3$
3,3 | $-25s^{4}-16s^{2}+16$ | $8(15s^{2}-14)$ | $15$
3,5 | $-3s^{4}$ | $-6(3s^{2}-2)$ | $9$
Table 2: Inclination polynomials $q_{k,i,j}$ in Eq. (26). $i,j$ | $k=0$ | $k=1$ | $k=2$
---|---|---|---
0,-2 | $0$ | $-64(3s^{2}-2)$ | $0$
0,-1 | $0$ | $-160(3s^{2}-2)$ | $0$
0,0 | $0$ | $-16(s^{2}+6)$ | $48$
0,1 | $0$ | $0$ | $48$
0,2 | $-24s^{2}(15s^{2}-4)$ | $0$ | $120$
0,3 | $-112s^{4}$ | $0$ | $0$
0,4 | $-72s^{4}$ | $0$ | $0$
1,-3 | $0$ | $-8(3s^{2}-2)$ | $0$
1,-2 | $0$ | $48(3s^{2}-2)$ | $0$
1,-1 | $12s^{2}(13s^{2}-4)$ | $24(5s^{2}-2)$ | $12$
1,0 | $0$ | $-16(75s^{2}-68)$ | $-18$
1,1 | $-12s^{2}(13s^{2}-4)$ | $-96(2s^{2}-1)$ | $60$
1,2 | $72s^{2}(5s^{2}-2)$ | $-32(15s^{2}-13)$ | $0$
1,3 | $4s^{2}(8-7s^{2})$ | $0$ | $120$
1,4 | $54s^{4}$ | $0$ | $66$
1,5 | $-12s^{4}$ | $0$ | $0$
2,-3 | $0$ | $6(3s^{2}-2)$ | $0$
2,-1 | $-2(27s^{4}+180s^{2}-128)$ | $18(17s^{2}-14)$ | $-9$
2,0 | $0$ | $-6(17s^{2}-18)$ | $-9$
2,1 | $2(27s^{4}+180s^{2}-128)$ | $-2(825s^{2}-742)$ | $-45$
2,2 | $6(5s^{4}+8s^{2}-8)$ | $-24(s^{2}-2)$ | $-15$
2,3 | $3s^{2}(39s^{2}-16)$ | $2(55s^{2}-34)$ | $-15$
2,4 | $0$ | $6(13s^{2}-10)$ | $-3$
2,5 | $9s^{4}$ | $0$ | $21$
2,6 | $0$ | $0$ | $3$
Table 3: Inclination polynomials $p_{k,i,j}$ in Eq. (26).

##  5 Performance of the analytical solution

For reference, in the first two examples we choose the physical parameters of
the Earth ($\mu=398600.44$ $\mathrm{km^{3}/s^{2}}$, $\alpha=6378.1363$ km,
$J_{2}=0.001082634$). The altitude of Earth flybys for gravity assist
typically range from about 300 km, as in the case of NASA’s Galileo
mission,444solarsystem.nasa.gov/missions/galileo/in-depth Retrieved July 18,
2022 to the more common several thousands of km, as for instance for NASA’s
Stardust mission.555solarsystem.nasa.gov/stardust/news/ega Retrieved July 18,
2022 Accordingly, we fix a perigee height of 1000 km over the Earth’s surface
in the examples, which is low enough to clearly undergo the oblateness
disturbing effect.

The first test case is for a high eccentricity hyperbolic-type orbit, for
which the flyby happens in a relatively short time, and both the Kepler
approach and the first order of the perturbation solution are expected to
provide acceptable results. The second example deals with a quasi-parabolic
orbit, for which the spacecraft remains much more time in the close vicinity
of the Earth, and, therefore, it is expected to be a more challenging test for
the analytical solution.

In the first case we choose initial conditions corresponding to the orbital
elements $a=2459.38$ km, $e=4$, $I=23.5^{\circ}$, $\Omega=60^{\circ}$,
$\omega=90^{\circ}$, $M=-21400^{\circ}$, that correspond approximately to the
entrance of the satellite into the Earth’s sphere of influence (SOI). The
“true” reference solution is computed from a numerical integration of these
initial conditions in the main problem dynamics that guarantees the
preservation of both the energy integral and the third component of the
angular momentum within at least 14 digits along the propagation interval of
36 hours. The reference orbit is shown in Fig. 1.

Figure 1: Earth’s flyby with $e=4$ and 1000 km periapsis distance. Left: true
orbit. Right: magnification on the closest approach. Distance in km.

Then, the initial conditions are propagated analytically, first in the
Keplerian dynamics, next, in the “common” interpretation of Deprit’s radial
intermediary, in what follows abbreviated DRI, in which Eqs. (16)–(17) are
assumed to remain valid in osculating elements, and finally in the “natural”
interpretation of DRI [6]; that is, the approximation of the main problem
dynamics provided by the new perturbation solution. The Root Sum Square (RSS)
errors of each analytical solution with respect to the true orbit are shown in
Fig. 2. The Keplerian propagation starts with zero error, and, because the
disturbing effect of $J_{2}$ is very small at SOI entrance, errors grow
slowly, reaching the meter level at a distance of approximately 25 Earth’s
radius, just a few hours before the closest approach. Errors grow notably
during the perigee passage, and, due to the different mean Keplerian motion
and the one of the main problem dynamics, increase almost linearly after that.
At the end of the propagation, the RSS errors of the Keplerian approximation
reach several hundreds of km.

Figure 2: RSS position errors of the flyby in Fig. 1 when using the Keplerian
dynamics, DRI-common, and the 1st order of the perturbation solution. Note the
logarithmic scale.

When using the “common” version of DRI, we get closer to the $J_{2}$ dynamics
although of a slightly different orbit due to the lack of transformation
between original and prime variables. By this reason, the errors in the
arrival branch of the predicted orbit grow faster than the Keplerian
approximation until the accumulation of nonlinear effects due to the $J_{2}$
perturbation is significant enough, which happens at about 1.5 Earth’s radius
for this orbit. Then, the better representation of the main problem dynamics
by DRI results in clearly smaller errors than those of the Keplerian
approximation during the flyby, and also in the departure branch of the orbit,
as illustrated with the dashed gray lines in both plots of Fig. 2. The
computational burden of DRI is notably alleviated when the root finding
procedure needed in the solution of the implicit function Eq. (18) is replaced
by the direct evaluation of Eq. (21), yet in this case the analytical
propagation starts with an initial error of half a meter —which is clearly
irrelevant at the precision of the analytical propagation.

Things notably improve when using the natural version of DRI, labelled “1st
order” in Fig. 2. Now we unavoidably deal with the mandatory transformation of
the initial osculating state to the mean elements that initialize the
constants of the perturbation solution —using the corrections in Eqs.
(11)–(15) with opposite signs— and the consequent recovery of osculating
elements using Eqs. (11)–(15) in order to compute the RSS error.666Recall
that, while direct and inverse transformations are just opposite for first
order corrections, they are fed with different sets of elements, either mean
or osculating, and hence their composition is affected by the truncation error
of the perturbation theory Errors provided by the 1st order solution grow fast
at perigee passage, but their following almost linear increase is,
approximately, $J_{2}$ times smaller than in the Keplerian and DRI cases due
to the much better approximation provided by the perturbation solution, and
barely reaches 100 m at the end of the propagation.

For the second case we choose the quasi-parabolic orbit $a=1.47563\times
10^{6}$ km, $e=1.005$, $I=23.5^{\circ}$, $\Omega=60^{\circ}$,
$\omega=90^{\circ}$, $M=-1^{\circ}$, which will fly about 1000 km over the
Earth’s surface too, but now during a much longer time than in the previous
example. The reference orbit, which is shown in Fig. 3, was numerically
integrated in the main problem dynamics up to 24 hours. As before, it is
compared with equivalent solutions provided by the Keplerian approximation, on
the one hand, and the torsion-based solution of the main problem, both in the
common and natural versions, on the other.

Figure 3: Quasi-parabolic flyby 1000 km above Earth’s surface. Distance in km.

Now the spacecraft undergoes the disturbing effects of the oblateness
perturbation for a much longer time than in the previous flyby, and hence the
improvements of DRI over the Keplerian approximation are much more evident in
this example, as shown in Fig. 4 —where the approximate inversion in Eq. (21)
has been used and hence the apparent initial error of the torsion-based
solution. Also, analogous improvements to the previous example are also
obtained by the 1st order perturbation solution in this case, for which the
errors remain $\mathcal{O}(J_{2})$ with respect to corresponding ones of the
Kepplerian approximation. However, as shown in Fig. 4, the errors of the first
order solution now clearly peak along the close visit to Earth reaching a RSS
error of about 700 m at perigee (about 10% of the minimum distance to the
Earth’s surface), which immediately falls by one order of magnitude, finally
reaching about 200 m at the end of the propagation. The latter figure means a
$\mathcal{O}(J_{2})$ improvement when compared to about the same number of km,
instead of m, provided by the Keplerian approach. The peak, we recall, is due
to the appearance of $\eta$ in denominators of some of the first order
corrections, which in this more challenging case takes the value $\eta\approx
0.1$. On the other hand, it deserves mentioning that the fact that the errors
of the 1st order solution decrease form the initial 10 meter level to just a
few meters when approaching Earth is just an apparent paradox that is easily
explained by the approximate truncation of the inversion of the torsion to the
first order of $J_{2}$, and the decreasing value of $\theta$ in Eq. (19).
Eventually, errors in the dynamical modeling provided by the natural version
of DRI overcome those due to the approximation in Eq. (21), and grow again
irrespective of the $\theta$ value.

Figure 4: RSS position errors of the flyby in Fig. 3 when using the Keplerian
dynamics and the 1st order of the perturbation solution. Note the logarithmic
scale.

Like in the case of bounded motion, a general improvement of the 1st order
solution is obtained when including second-order terms in the mean elements
propagation, a strategy that negligibly increases the computational burden [1,
2]. In our case, this means including the term in Eq. (22) in the quasi-
Keplerian Hamiltonian, as well as Eqs. (23) and (24) in the solution of the
torsion transformation. As shown in Fig. 5, while these additional terms
generally reduce the errors of the analytical propagation (labelled “1st order
+”), they do not show improvements during the close approach to the Earth.
Certainly, computing the second order corrections derived from Eq. (26) makes
a full second order solution that, as expected, generally improves the
accuracy to $\mathcal{O}(J_{2}^{2})$, as is also shown in Fig. 5 (black-dashed
curve labelled “2nd order”). The fact that the expected improvement does not
happen at perigee to the whole extent, is due to the appearance in
denominators of the second order corrections of higher powers of $\eta$ than
in the first order approach, which now reach the fourth power in the case of
$\theta_{0,2}$ and $R_{0,2}$. Still, the errors of the 2nd order solution at
perigee passage remain at least 20 times smaller than those of the first order
solution. Finally, the RSS error curve labelled “2nd order +” includes third-
order secular terms obtained from the elimination of the parallax in the
perturbation solution. Like in the first order case, these additional terms
slightly improve the computation of the “mean” mean motion —the mean motion of
the quasi-Keplerian Hamiltonian— yet their effects are only noted in the
departure branch of the orbit.

Figure 5: RSS position errors of the flyby in Fig. 3 when using higher orders
of the perturbation solution. Note the logarithmic scale.

Because perturbation theories of this kind are rather uncommon, the final
example is intended to show that the convergence of the successive orders of
the perturbation solution, which is commonly observed in the case of bounded
perturbed Keplerian motion (see [20, 17], for instance), also happens in the
current case of perturbed hyperbolic Keplerian motion. To do that, we test our
solution with the challenging case of a fish-type orbit [24]. However, since
parabolic orbits are necessarily excluded from the current perturbation
solution due to the $\eta$ divisor that affects the elimination of the
parallax transformation, we rather base our tests on an equivalent hyperbolic
branch of the Newtonian rosetta.

Thus, for the gravitational parameters of Jupiter ($\mu=1.268\times
10^{8}\,\mathrm{km^{3}/s^{2}}$, $\alpha=71492$ km, $J_{2}=0.01475$) we choose
the initial conditions of an equatorial orbit with $a=724920$ km, $e=1.1$,
$g=270^{\circ}$, and $\ell=-120^{\circ}$. The fact that this orbit would
impact Jupiter does not worry us in our aim of validating the perturbation
theory. Results are summarized in Fig. 6, where the true orbit, computed
numerically, is represented with a black curve, approaching to Jupiter from
the top-left corner, surrounding the attraction center by more than
$180^{\circ}$ to first cross the approaching branch of the orbit and then
leave the scene by the left side of the picture. The Keplerian approximation
(gray curve of Fig. 6) only mimics the true trajectory in the approaching
part, and then takes a clearly divergent path. DRI (gray dashed curve) is able
to turn the departure branch towards the true solution, a result that is
clearly improved by the first order theory (gray dotted curve), yet not enough
enhanced to produce the crossing between the approaching and departing
branches of the orbit. This intersection is obtained with the second order
solution (black dotted curve of Fig. 6), yet it happens out of the frame of
Fig. 6, clearly far away from the true event. The third-order approximation
(black dot-dashed curve of Fig. 6) notably improves the results, and the
crossing now occurs inside the displayed area in Fig. 6. Finally, the black-
dashed curve, which was computed with a fifth-order truncation of the
elimination of the parallax complemented with additional secular terms coming
from a sixth-order truncation of the quasi-Keplerian Hamiltonian, reasonably
approximates the true orbit. It is expected that higher order truncations
continue to converge slowly to the true solution, yet we did not compute them
due to the important growth in the size of the perturbation series and the
corresponding exponential increase of the computational burden. We are
satisfied with these results, which show the correctness of Hori’s approach
for dealing with unbounded perturbed Keplerian motion.

Figure 6: True hyperbolic rosetta branch (black curve) and successive
approximations of the perturbation solution truncated to the zeroth
(Keplerian, gray), first (gray, dotted), second (black, dotted), third (black,
dot-dashed) and fifth order (black, dashed) of Jupiter’s main problem. The
gray-dashed curve was generated with the “common” version of DRI. Distances
are km.

##  6 Remarks on the boundary condition: the third body perturbation case

Similarly to Hori’s approach, the transformation that yields Deprit’s radial
intermediary depends on the choice made for the integration constant. The main
problem straightforwardly suggests to impose null perturbing effects at
infinity, as at great distances from the primary the non-uniform gravity field
of the planet flown by is negligible.

Yet, in practical examples for the Solar or a planetary system, the great
distances from the flyby body are characterized by stronger effects of the
main attractor of the whole system. We therefore proceed at assessing whether
perturbation approaches of this kind, with Keplerian boundary condition, are
suitable to study the perturbed hyperbolic Kepler problem with third body
effects only.

Identifying with the subscript $B$ the perturbing body, the full Hamiltonian
of the third body perturbation case is

$\mathcal{M}=\mathcal{M}_{0}+\mathcal{M}_{B}=-\frac{\mu}{2a}-\mu_{B}\bigg{(}\frac{1}{|\mathbf{r}-\mathbf{r}_{B}|}-\frac{\mathbf{r}\cdot\mathbf{r}_{B}}{r_{B}^{3}}\bigg{)}$
(27)

Following Lidov and Kozai [23, 16], we replace
$\mathbf{r}\cdot\mathbf{r}_{B}=rr_{B}\cos\psi$ and model the perturbing
potential $\mathcal{M}_{B}$ as a series expansion on the Legendre polynomials
of $\cos\psi$:

$\mathcal{M}_{B}=-\frac{\mu_{B}}{r_{B}}\sum_{j=2}^{+\infty}\frac{\varepsilon^{j}}{j!}j!\,\frac{r^{j}}{r_{B}^{j}}P_{j}(\cos\psi)$
(28)

For instance, $P_{0}(x)=1$, $P_{1}(x)=x$, $P_{2}(x)=\frac{1}{2}(3x^{2}-1)$,
$P_{3}(x)=\frac{1}{2}(5x^{3}-3x)$, and so on. The small parameter
$\varepsilon$ is just a token to identify the strength of the perturbation
effects with the ratio $r/r_{B}$, that is, retaining few orders only provide
an accurate approximation if $r\ll{r}_{B}$. More generally, the series
converges if $r<r_{B}$. Note the absence of the first two orders, consequence
of $P_{0}$ leading to the constant expression $\mu_{B}/r_{B}$ and the first
order term canceling with the tidal part of the perturbing potential.

Following the same principle of Hori’s approach for the main problem, we seek
a condition where the third body perturbation vanishes by itself. We
approximate

$\mathcal{M}_{B}\approx-\frac{\mu_{B}}{r_{B}}2\,\frac{r^{2}}{r_{B}^{2}}\frac{1}{2}\big{(}3\cos^{2}\psi-1\big{)}$
(29)

and observe that $\cos^{2}\psi=1/3$ makes the second term vanish. We could be
tempted to use this particular point to set the Keplerian boundary condition,
however, this unique setting corresponds to the second order term only being
null, rather than the full original potential. In addition, trajectories such
that $\cos^{2}\psi=1/3$ is never verified may exist. Considering the higher
order terms, we do not find a value of $\cos\psi$ that makes all the terms in
the expansion to vanish.

These considerations suggest that the boundary condition derived from a
Keplerian hyperbola may not be suitable for the third body perturbation case.
This aspect remarks that this method remains valid for the main problem only,
or more in general for vanishing perturbations at infinity. We strengthen our
claim with a numerical proof of the just mentioned observation, whose
analytical equations of motion in prime variables we obtained with the help of
a symbolic processor, following the same steps of Hori’s approach but applying
the Lie transform method. We do not include them in this article for the sake
of brevity. We study the case of a a planetary flyby of Venus modeling the
third body perturbation of the Sun, within Venus’ SOI. We consider the SOI
entrance state, with position
$\mathbf{r}^{T}=[-104551.25,-597237.74,-110314.51]$, km, and velocity
$\mathbf{v}^{T}=[3.25,17.76,3.67]$, km/s, at the time $t_{0}=8119.84$ MJD2000,
in the Ecliptic J2000 reference frame centered on Venus. We simulate the flyby
motion using an ephemeris model that includes the Solar System planets, the
Moon and general relativity effects up to the "critical" point where
$\cos\psi=-1/\sqrt{3}$. Then, we compare in Fig. 7 the consequent forward
simulations, starting from this critical point and up to the SOI exit, of the
Lie transform-based approach, the numerical simulation of the second order
term only of the third body potential, the Keplerian integration, and the full
third body effect (not expanded).

Figure 7: Evolution of the error for the third body models, relative to the
ephemeris solution.

In clear contrast with the improvements obtained for the main problem case, we
observe from Fig. 7 the divergent behavior of the analytical perturbation
solution. Setting the boundary condition as the Keplerian solution where the
perturbing potential vanishes, in analogy with Hori’s approach, the so
obtained analytical solution provides even worse predictions than the
Keplerian only case, without converging toward the numerically simulated
trajectory for the same physical model. Other types of boundary conditions
should be explored in dedicated works, for instance the magnitude difference
between the perturbation and the Keplerian term in specific regions of the
hyperbolic trajectory.

##  7 Conclusions

Most main problem intermediaries in the theory of artificial satellites were
originally computed in polar variables, and hence are not constrained to the
case of bounded orbits. Therefore, the improvements they provide in their
_common_ realization over the Keplerian approximation encompass also the case
of unbounded motion, and are clearly leveraged in their _naturalized_
versions. However, in the latter case, the integration “constant” on which the
perturbation theory depends upon can no longer remain arbitrary. On the
contrary, it must be determined in such a way that the perturbation solution
fulfills the boundary conditions at infinity derived from the Keplerian
hyperbola.

In particular, we have shown that Deprit’s radial intermediary provides an
efficient alternative to the Keplerian approximation commonly used in flyby
design. Moreover, because the appearance of secular terms is not at all of
concern for hyperbolic-type motion, higher order refinements of Deprit’s
quasi-Keplerian approach are readily computed. Because of that, the analytical
computation of unbounded perturbed Keplerian motion can take full advantage of
the convergence to the true solution provided by the computation of
consecutive higher orders of the perturbation approach, in which case the
accuracy of the analytical solution is notably increased.

## Acknowledgements

The research has received funding from the European Research Council under the
European Union’s Horizon 2020 research and innovation program (grant agreement
679086 - COMPASS). ML also acknowledges partial support from the Spanish State
Research Agency and the European Regional Development Fund (Projects
PID2020-112576GB-C22 and PID2021-123219OB-I00, AEI/ERDF, EU).

## References

  * [1] D. Brouwer.  Solution of the problem of artificial satellite theory without drag.  The Astronomical Journal, 64:378–397, 1959.
  * [2] S. Coffey and K. T. Alfriend.  An analytical orbit prediction program generator.  Journal of Guidance, Control and Dynamics, 7(5):575–581, 1984.
  * [3] Z. Dang, J. Luo, P. Shi, and H. Zhang.  General Characteristics of the Motion on $J_{2}$-Perturbed Equatorial Orbits.  Journal of Guidance Control Dynamics, 42(10):2319–2324, 2019.
  * [4] C. E. Delaunay.  La Théorie du Mouvement de la Lune, Premier volume, volume 28 of Mémoires de l’Academie des Sciences de l’Institut Impérial de France. Mallet-Bachellier, Paris, 1860.
  * [5] A. Deprit.  Canonical transformations depending on a small parameter.  Celestial Mechanics, 1(1):12–30, 1969.
  * [6] A. Deprit.  The elimination of the parallax in satellite theory.  Celestial Mechanics, 24(2):111–153, 1981.
  * [7] A. Deprit and A. Rom.  The main problem of artificial satellite theory for small and moderate eccentricities.  Celestial Mechanics, 2(2):166–206, 1970.
  * [8] S. Ferrer and M. Lara.  Families of Canonical Transformations by Hamilton-Jacobi-Poincaré Equation. Application to Rotational and Orbital Motion.  Journal of Geometric Mechanics, 2(3):223–241, 2010.
  * [9] L. Floria.  A Simple Derivation of the Hyperbolic Delaunay Variables.  The Astronomical Journal, 110:940–942, 1995.
  * [10] P. Gurfil and M. Lara.  Satellite onboard orbit propagation using Deprit’s radial intermediary.  Celestial Mechanics and Dynamical Astronomy, 120(2):217–232, 2014.
  * [11] D. Hautesserres and M. Lara.  Intermediary LEO propagation including higher order zonal harmonics.  Celestial Mechanics and Dynamical Astronomy, 127:505–526, 2017.
  * [12] L. M. Healy.  The main problem in satellite theory revisited.  Celestial Mechanics and Dynamical Astronomy, 76(2):79–120, 2000.
  * [13] G.-I. Hori.  The motion of a hyperbolic artificial satellite around the oblate earth.  The Astronomical Journal, 66:258–263, 1961.
  * [14] G.-I. Hori.  Theory of general perturbation with unspecified canonical variables.  Publications of the Astronomical Society of Japan, 18(4):287–296, 1966.
  * [15] D. J. Jezewski.  An analytic solution for the $J_{2}$ perturbed equatorial orbit.  Celestial Mechanics, 30(4):363–371, 1983.
  * [16] Y. Kozai.  Secular perturbations of asteroids with high inclination and eccentricity.  The Astronomical Journal, 67:591–598, 1962.
  * [17] M. Lara.  Solution to the main problem of the artificial satellite by reverse normalization.  Nonlinear Dynamics, 101(2):1501–1524, 2020.
  * [18] M. Lara.  Brouwer’s satellite solution redux.  Celestial Mechanics and Dynamical Astronomy, 133:47, 2021.
  * [19] M. Lara.  Hamiltonian Perturbation Solutions for Spacecraft Orbit Prediction. The method of Lie Transforms, volume 54 of De Gruyter Studies in Mathematical Physics.  De Gruyter, Berlin/Boston, 1 edition, 2021.
  * [20] M. Lara, J. F. Palacián, and R.P. Russell.  Mission design through averaging of perturbed Keplerian systems: The paradigm of an Enceladus orbiter.  Celestial Mechanics and Dynamical Astronomy, 108(1):1–22, 2010\.
  * [21] M. Lara, J. F. San-Juan, and L. M. López-Ochoa.  Delaunay variables approach to the elimination of the perigee in Artificial Satellite Theory.  Celestial Mechanics and Dynamical Astronomy, 120(1):39–56, 2014.
  * [22] M. Lara, J. F. San-Juan, and L. M. López-Ochoa.  Proper averaging via parallax elimination.  In S. B. Broschart, J. D. Turner, K. C. Howell, and F. R. Hoots, editors, Astrodynamics 2013, volume 150 of Advances in the Astronautical Sciences, pages 315–331, P.O. Box 28130, San Diego, California 92198, USA, January 2014. American Astronautical Society, Univelt, Inc.
  * [23] M. L. Lidov.  Evolution of artificial planetary satellites under the action of gravitational perturbations due to external bodies.  Iskusstviennye Sputniki Zemli, 8:5–45, 1961.
  * [24] V. Martinusi and P. Gurfil.  Analytical solutions for $J_{2}$-perturbed unbounded equatorial orbits.  Celestial Mechanics and Dynamical Astronomy, 115(1):35–57, 2013.
  * [25] H. von Zeipel.  Research on the Motion of Minor Planets.  NASA TT F-9445, 1965.  (NASA Translation of: Recherches sur le mouvement des petites planètes, Arkiv för matematik, astronomi och fysik, vol. 11, 1916, vol. 12, 1917, vol. 13, 1918).


# Legal Search in Case Law
and Statute Law

Julien Rossi Corresponding Author<EMAIL_ADDRESS>Evangelos Kanoulas
Amsterdam Business School, University of Amsterdam  Institute of Informatics,
University of Amsterdam

###### Abstract

In this work we describe a method to identify document pairwise relevance in
the context of a typical legal document collection: limited resources, long
queries and long documents. We review the usage of generalized language
models, including supervised and unsupervised learning. We observe how our
method, while using text summaries, overperforms existing baselines based on
full text, and motivate potential improvement directions for future work.

###### keywords:

Legal Search Case Law Statute Law Language Models Transfer Learning

,

##  1 Introduction

The increasing amount of legal data available requires the ability to search
and identify relevant information in this data, it calls for the ability to
automate or assist in specialized retrieval tasks, as a service to either the
public or practitioners. Legal Search is the specialized Information Retrieval
(IR) task dealing with legal information relevant to a situation described
factually and legally. The legal documents collections differ from generic IR
collections by the level of professional knowledge needed to produce a labeled
corpus [1]. This results in limited and specialized collections, each covering
a narrow field of interest: lease contracts, financial products, court
decisions, etc…Existing literature shows how Legal Search is an essential tool
for practitioners and citizens in need of critical information [2, 3, 4],
whereas it introduces difficult challenges that our work addresses.

In this paper, we focus on limitations of the applicability of generalized
language models [5, 6, 7, 8] in the context of the handling of legal
information tasks. We observe that specific features of typical legal
documents adversely affects those models, such as long texts, mixing of formal
abstract concepts with casual verbatims in layman language, or dealing with an
unclear formulation of relevance.

The development of Generalized Language Models such as BERT [7] has allowed
the application of complex neural models to tasks even with limited data
available, significantly improving over lexical methods, while past neural
models were having difficulties to deal with smaller datasets, due to
vanishing gradients, and were difficult to compute and implement on hardware
[9].

The core of these generalized language models is the separation between
learning the language (pre-training) and adapting to a task (fine-tuning). The
adaptation to legal tasks is complicated by the difference between the generic
language and the legalese language, a difference in word usage and meaning,
information structuring, or sentence complexity [10].

Our contribution is a method to improve Legal Information Retrieval, where the
ranking problem is reformulated as a pairwise relevance score problem, modeled
as a fine-tuning task of a generalized language model, and where we adapt long
documents to limited input sequence length by summarizing. Furthermore we show
how additional pre-training based on a narrow selection of texts can improve
performance.

We will focus on the following questions:

  * •

RQ1 Can we use summary encoding as a dense representation of long documents ?

  * •

RQ2 Can we use pre-training and fine-tuning of neural language model with
limited data to learn a specific legal language ?

In section 2, we will introduce our methods to answer the research questions,
then in section 3 the tasks we will use for our work, then the research for
each task will be described with results and analysis in sections 4 and 5. We
will present our global analysis and conclusion in section 6, along with
motivation for future research.

##  2 Methods

Our work reformulates the ranking problem into a pairwise relevance
classification problem, where we train a classifier to separate relevant pairs
(made of a query case and a noticed case) from irrelevant pairs of documents.
For each query, we let the model infer the probability for the positive class
for each pair made of the query and a candidate document, and then rank the
candidate documents according to that score.

This method follows after the well-established probability ranking principle
in IR [11]. For the tasks we consider, we argue that the relevance of a
document with regards to a given query is independent from other documents of
the collection, confirming the applicability of this framework to our tasks.

Within this framework, we introduce input pre-processing as well as model
training steps that will help us address the research questions. We refer to
[7] for details on the topology of BERT as well as for the setup of the pre-
training task or the fine-tuning task.

###  2.1 Pairwise Embeddings

We make use of the capacity of BERT models to encode single dense embeddings
for a pair of texts. These embeddings encode the interaction between the texts
in the pair, so that it is a suitable input for a downstream pair classifier.
In a BERT system, this corresponds to the embeddings of [CLS] token, inserted
at the beginning of each input sequence.

###  2.2 Fine-Tuning of Pre-Trained Language Model

The fine-tuning is a supervised learning task for a complete model comprising
a pre-trained BERT and a fully-connected layer used for classification. The
learning is driven by the true classification labels of the examples, and all
the weights of the combined model are updated during this training. The usual
loss for classification is the Cross-Entropy. During this phase, the combined
model learns how to perform the classification based on the pairwise
embeddings.

###  2.3 In-Domain Additional Pre-Training of Language Model

This task is an unsupervised learning task, where an already pre-trained BERT
model receives further pre-training with additional texts. We make the remark
that BERT is a NLP system, with a high capacity of understanding at the
language level. For example [12, 13] demonstrate how the layers in BERT work
as a hierarchical system with a capacity to discover syntactic features of the
text, while the semantic features are assumed to emerge naturally from co-
location and word contexts. In that regard, legal texts create a variety of
domain specific languages, considering different domains such as contracts,
court decisions, legislation, etc.

Legal languages differ from casual languages with unusual vocabulary (rare
words, latin words, etc.), semantics (words where the legal meaning differs
from casual usage), and syntactic features. We suggest that a pre-trained
model would learn specific knowledge from a limited additional pre-training on
legal domain-specific texts. We base our suggestion on the fact that within a
specific legal domain, the text features will be fairly uniform. Following
[14], we expect pre-training on the available material being more beneficial
to the capacities of the system than only fine-tuning for the task. We expect
to answer RQ2 with this method.

###  2.4 Summarization of Long Documents

This method will allow us to address our question RQ1, on long documents. As
we base our models on BERT, we are constrained by a maximum input length of
512 WordPiece [15] tokens. WordPiece will subdivide each word into multiple
tokens, we consider that this will nearly double the number of tokens observed
under a standard Punkt [16] tokenizer as implemented in NLTK [17].

For this reason, we introduce an extractive summarization of the texts in the
corpus using TextRank [18]. TextRank is an extractive summarization model,
based on a graph ranking model operating at sentence level. As it extracts
full sentences, it preserves the structure of a natural language and makes the
summary fit for input to a Neural Language Model. We choose to limit the size
of the summary to 180 words, so a concatenated pair of texts will not be
longer than 512 WordPiece tokens. In the case this limit of 512 is not
respected, then the sequence will drop the last tokens from the text of
candidate case. We use the implementation from gensim [19].

##  3 Legal Retrieval Tasks

We will illustrate our work with 2 different tasks and their associated
document collection, taken from COLIEE [20]. The 2 tasks illustrate different
aspects of our research:

  * •

Case Law Retrieval: Find past cases similar to a query case. In this task, the
query and the documents are both long texts, and the amount of available data
is limited. We will address RQ1 and RQ2 with this task.

  * •

Statute Law Retrieval: Find law articles relevant to a situation. In this
task, the query and documents are short snippets of text, and only a small
amount of labeled data is available. We will address RQ2 with this task.

###  3.1 Case Law Retrieval

In this task, a given court case is considered as a query to retrieve
supporting cases (also named ’noticed cases’) for the query case. A noticed
case supports the decision taken in the query case, although the final
decision itself is irrelevant in our retrieval. What matters is the proximity
of the themes that are tackled by the query case and the noticed cases, either
legal themes or narrative themes.

The dataset is drawn from an existing collection of Federal Court of Canada
case law, provided by vLex Canada111http://ca.vlex.com. Each query case is
given a collection of 200 potential supporting cases (also named ’candidate
cases’), these collections are provided labeled for the training dataset, and
unlabelled for the unknown test dataset. The training dataset contains 285
query cases. All cases are relative to Citizenship and Immigration
proceedings.

We consider the corpus as the complete collection of query cases, and
candidate cases from the training dataset, Figures 2 and 2 illustrate the
distribution of document length across the corpus, measured in number of
words. With a median length of 2500 words, this dataset illustrates well the
challenge of dealing with long documents and long queries, as formulated in
our RQ1. The corpus is a limited resource with only 285 queries, which fits
well with our RQ2.

Figure 1: Distribution of the number of words per document

Figure 2: Cumulative distribution for the number of words per document

###  3.2 Statute Law Retrieval

In this task, a brief query has to be disputed with regards to the existing
legislation, in this case limited to the Japanese Civil Code. The retrieval
task identifies which articles of the Civil Code are relevant to the
evaluation of the legal validity of the query.

The data for this task is compiled from Japanese Bar Exams, and provided both
in Japanese and English. We focus on the English dataset from the 2018
edition. The training dataset contains 651 queries, while the test dataset is
made of 69 queries. This amount of data will be a good fit for our RQ2. The
queries are short snippets with an average of 40 words, with a maximum of
under 120 words, while the average law article has 60 words, with 98% of
articles having less than 200 words. We will not face the issue of long texts.
We consider it a limited resource with few examples of relevant matches, as
94% of queries have 1 or 2 relevant articles, among around 1000 candidate
articles.

###  3.3 Evaluation Metrics

In the COLIEE setting, the system decides for each query the number of
candidates it returns, and it is evaluated based on the F1 or F2 score of that
sub-list. In this paper, we will consider a more traditional approach with
respect to Information Retrieval, grounded in a realistic use case for the
Case Law Retrieval Task: the end user is a staff from a judge’s office, and
the queries are submitted to a computer software and returned as a ranked
list. We use the same metrics for the Statute Law Retrieval Task.

We choose to report Precision at R (P@R), where R is the total number of
relevant documents, Recall and Precision at k (R@k, P@k) and Mean Average
Precision (MAP). For each metric, we will consider the macro average222$m_{i}$
is the value of the metric for sample $i\in[1,N]$, then
$macro(m)=\frac{1}{N}\sum_{i=1}^{N}m_{i}$. Metrics are computed with
trec_eval333https://github.com/usnistgov/trec_eval, and Precision-Recall curve
is plotted with plot-trec_eval444https://github.com/hscells/plot-trec_eval

##  4 Case Law Retrieval

###  4.1 Experimental Setup

Following the ”Methods” section 2, we pre-process the collection by
summarizing all texts. The final dataset is a collection of triplets made of
query text, candidate text and relevance judgment. We will fine-tune a pre-
trained BERT model for a binary classification task, using the relevance
labels from the labeled training dataset. We will also perform additional in-
domain pre-training, and then fine-tune this model for the same classification
task. In the results section 4.3, the models appear under the names FineTuned
and PreTrained.

For the in-domain pre-training, we use the entire corpus of court decisions as
the pre-training corpus. This corpus has 18000 documents, and a total of 45
million tokens. The training will start with a pre-trained BERT model (bert-
base-uncased), using the Megatron-
LM555https://github.com/cybertronai/Megatron-LM library. We run the pre-
training on 1 GPU nVidia Tesla P40, with 24GB onboard RAM during 24 hours, or
70000 iteration steps. We export this model for fine-tuning with fast-
bert666https://github.com/kaushaltrivedi/fast-bert, during 10 hours for 4
epochs.

The labeled dataset is split 75%/25% between training data and evaluation
data. The dataset is split according to the cases, so that the cases in the
evaluation dataset are not in the training dataset. This split reflects
properly the capacity of the system to generalize to unseen data. In this
setting, there will also be candidate cases that are only in the evaluation
dataset, and therefore unseen data.

For both models, after the fine-tuning step is finalized, the trained model
computes the pairwise relevance score for each pair of query case and
candidate case in the test dataset, the score for the positive class is used
to rank the candidate cases of each query case.

Our implementation is based on libraries pytorch-
transformers777https://github.com/huggingface/pytorch-transformers, using
PyTorch888https://pytorch.org/ framework.

###  4.2 Baselines

We consider 2 lexical features baselines based on BM25. As we setup our BERT-
based models to learn from the summaries of the full documents, we propose to
have a baseline that represents as well this limitation in the availability of
information:

  * •

One baseline will use BM25 score as the pairwise relevance score, when the
corpus contains all full texts documents

  * •

Another baseline will use BM25 score as the pairwise relevance score, with a
corpus made of the summaries of all documents

Further in the Results section, we will also introduce the Perfect Ranker,
this ideal system ranks first all noticed cases, and then all unnoticed cases.
It will provide us with an upper bound for the performances at each rank.

###  4.3 Results and Analysis

System | R@10 | P@10 | R@1 | P@1 | P@R | MAP
---|---|---|---|---|---|---
BM25 Summaries | 0.14 | 0.07 | 0.02 | 0.11 | 0.11 | 0.14
BM25 | 0.76 | 0.36 | 0.32 | 0.70 | 0.68 | 0.73
PERFECT | 0.97 | 0.50 | 0.41 | 1.00 | 1.00 | 1.00
Fine-Tuned | 0.75 | 0.34 | 0.31 | 0.80 | 0.63 | 0.70
Pre-Trained | 0.81 | 0.39 | 0.34 | 0.90 | 0.73 | 0.79
Table 1: Table of Results for All Systems Figure 3: Precision-Recall curves

We first observe that all of our models outperform the baseline of BM25 with
summarized texts (model named ’BM25Summ’), and perform at a level comparable
to BM25 on full texts. We interpret this part as a validation of our approach,
and of the underlying assumption that summarization was a way to capture the
information necessary to establish relevance, therefore bringing a positive
answer to our research question RQ1.

The Precision-Recall curves in Figure 3, show the skills of the introduced
systems. It is noticeable that the Fine-Tuned system shows a consistent better
performance at all levels of Recall.

With regards to statistical significance, we computed the p-values of an OLS
regression model with the target metric as the response variable, and a binary
dummy (False for the Baseline, True for the tested model) as the independent
variable. We made use of the python library Statsmodels [21]. The p-value for
the coefficient of the dummy indicates whether the observed difference in mean
is significant or not. We use BM25 as the baseline, and observe a positive
effect for both PreTrained and FineTuned models, while this effect is
significant only for the PreTrained model.

We attribute this improvement to the knowledge gained during the additional
pre-training. The amount of data available for our additional pre-training is
orders of magnitude smaller that the material used for the initial pre-
training, if we compare millions of words (our data) with tens of billions of
words (Wikipedia, BookCorpus, etc.), so we consider this result as a positive
answer to our research question RQ2.

##  5 Statute Law Retrieval

###  5.1 Experimental Setup

Following the ”Methods” section 2, we pre-process the collection by
summarizing all texts. The final dataset is a collection of triplets made of
query text, law article text and relevance judgment. We will fine-tune a pre-
trained BERT model for a binary classification task, using the relevance
labels from the labeled training dataset. We will also perform additional in-
domain pre-training, and then fine-tune this model for the same classification
task. In the results section 5.3, the models appear under the names FineTuned
and PreTrained.

For the in-domain pre-training, we used a collection of english translations
of court decisions from the Japanese Supreme Court, scraped from their
website999http://www.courts.go.jp/app/hanrei_en/search. This is a rather small
dataset of 1500 texts, for a total of 230000 tokens. The training will start
with a pre-trained BERT model (bert-base-uncased). We run the pre-training on
1 GPU nVidia Tesla P40, with 24GB onboard RAM during 2 hours, or 70000
iteration steps. We export this model for fine-tuning of 2 hours for 4 epochs.

For both models, after the fine-tuning step is finalized, the trained model
computes the pairwise relevance score for each pair of query text and law
article text in the test dataset, the score for the positive class is used to
rank the law articles for each query.

###  5.2 Baselines

We use the SOTA from COLIEE 2018 as the baseline, namely UB3 [22]. This system
is based on Terrier101010http://terrier.org/ and uses
TagCrowd111111https://tagcrowd.com/ for query reduction by extracting
keywords.

###  5.3 Results and Analysis

We report less metrics than for the previous task. As 98% of queries have 1 or
2 relevant articles, we do not report Precision for ranks higher than 1.
Results from different systems are compiled in Table 2.

Name | R@5 | R@10 | R@30 | MAP | P@1
---|---|---|---|---|---
UB3 | 0.7978 | 0.8539 | 0.9551 | 0.7988 | unknown
FineTuned | 0.9010 | 0.9203 | 0.9686 | 0.8246 | 0.7971
PreTrained | 0.8913 | 0.9130 | 0.9444 | 0.8321 | 0.8261
Table 2: Table of Results for All Systems

Our work overperforms the baseline significantly. When focusing on only the
models we introduce, we observe no significant differences in the evaluation
metrics, with the exception of P@1 where the model with additional pre-
training significantly improves over the ”off the shelf” model.

In this context, the additional pre-training yield a better performance for
some metrics, but not others. The observed significant improvement on the P@1
metric is of high interest in the setting of a Question Answering task. We
will consider this as a conditional positive answer to our research question
RQ2, as the improvement might realize only on some specific metrics.

##  6 Conclusion

We have demonstrated how models based on Generalized Language Models could
perform in the context of Legal Information Retrieval, in the presence of
limited data. We introduce a document summarization step in order to
accomodate the sequence length limitations of BERT. In this setting, our
system significantly improves over a BM25 system operating on full text.

We introduced an additional step of pre-training for existing models, which
provided significant improvement in the task with the largest amount of
training material. We leave for future work the transfer of this method to
other specific legal domains.

We consider for future work the possibilities of other new training tasks,
considering Masked Language Modeling and Next Sentence Prediction as the way
to establish a comprehension at semantic level, while other tasks would
contribute to learn the deeper knowledge needed to achieve higher performance
on the retrieval tasks.

## Use of COLIEE Data

We present these research findings based on the COLIEE Dataset for the Legal
Case Retrieval Task, in accordance with the ”MEMORANDUM ON PERMISSION TO USE
ICAIL 2019 PARTICIPANT DATA COLLECTION”. While we use some competition models
as baseline, we do not claim that this paper is an entry to the official
COLIEE competition. In the ”Evaluation Metrics” section 3.3, we elaborate on
how the metrics used for competition ranking differ from the metrics we
present. We refer to the proceedings of ICAIL 2019 [23] for further reading.

## References

  * [1] P. Arora, M. Hossari, A. Maldonado, C. Conran, G.J. Jones, A. Paulus, J. Klostermann and C. Dirschl, Challenges in the development of effective systems for Professional Legal Search (2018).
  * [2] I. Nejadgholi, R. Bougueng and S. Witherspoon, A Semi-Supervised Training Method for Semantic Search of Legal Facts in Canadian Immigration Cases., in: JURIX, 2017, pp. 125–134.
  * [3] M. Maclean, J. Eekelaar and B. Bastard, Delivering family justice in the 21st century, Bloomsbury Publishing, 2015\.
  * [4] J.-P. Boyd, Alienated children in family law disputes in British Columbia, Canadian Research Institute for Law and the Family, University of Calgary, 2015\.
  * [5] A. Radford, J. Wu, R. Child, D. Luan, D. Amodei and I. Sutskever, Language Models are Unsupervised Multitask Learners (2019). https://d4mucfpksywv.cloudfront.net/better-language-models/language_models_are_unsupervised_multitask_learners.pdf.
  * [6] J. Howard and S. Ruder, Fine-tuned Language Models for Text Classification, CoRR abs/1801.06146 (2018). http://arxiv.org/abs/1801.06146.
  * [7] J. Devlin, M. Chang, K. Lee and K. Toutanova, BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding, CoRR abs/1810.04805 (2018). http://arxiv.org/abs/1810.04805.
  * [8] L. Weng, Generalized Language Models, 2019\. https://lilianweng.github.io/lil-log/2019/01/31/generalized-language-models.html.
  * [9] E. Culurciello, The fall of RNN / LSTM, 2018\. https://towardsdatascience.com/the-fall-of-rnn-lstm-2d1594c74ce0.
  * [10] P. Tiersma, The nature of legal language (2008), 7–25.
  * [11] S.E. Robertson, The probability ranking principle in IR, Journal of documentation 33(4) (1977), 294–304.
  * [12] Y. Lin, Y.C. Tan and R. Frank, Open Sesame: Getting Inside BERT’s Linguistic Knowledge (2019). http://arxiv.org/abs/1906.01698.
  * [13] K. Clark, U. Khandelwal, O. Levy and C.D. Manning, What Does BERT Look At? An Analysis of BERT’s Attention (2019). http://arxiv.org/abs/1906.04341.
  * [14] S. Ruder, M.E. Peters, S. Swayamdipta and T. Wolf, Transfer Learning in Natural Language Processing, in: Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Tutorials, 2019, pp. 15–18.
  * [15] Y. Wu, M. Schuster, Z. Chen, Q.V. Le, M. Norouzi, W. Macherey, M. Krikun, Y. Cao, Q. Gao, K. Macherey, J. Klingner, A. Shah, M. Johnson, X. Liu, L. Kaiser, S. Gouws, Y. Kato, T. Kudo, H. Kazawa, K. Stevens, G. Kurian, N. Patil, W. Wang, C. Young, J. Smith, J. Riesa, A. Rudnick, O. Vinyals, G. Corrado, M. Hughes and J. Dean, Google’s Neural Machine Translation System: Bridging the Gap between Human and Machine Translation, CoRR abs/1609.08144 (2016). http://arxiv.org/abs/1609.08144.
  * [16] T. Kiss and J. Strunk, Unsupervised Multilingual Sentence Boundary Detection, Comput. Linguist. 32(4) (2006), 485–525. doi:10.1162/coli.2006.32.4.485.
  * [17] S. Bird, E. Klein and E. Loper, Natural Language Processing with Python, O’Reilly Media, 2009\.
  * [18] F. Barrios, F. López, L. Argerich and R. Wachenchauzer, Variations of the Similarity Function of TextRank for Automated Summarization, CoRR abs/1602.03606 (2016). http://arxiv.org/abs/1602.03606.
  * [19] R. Řehůřek and P. Sojka, Software Framework for Topic Modelling with Large Corpora, in: Proceedings of the LREC 2010 Workshop on New Challenges for NLP Frameworks, ELRA, Valletta, Malta, 2010, pp. 45–50, http://is.muni.cz/publication/884893/en.
  * [20] C.O. Group, COLIEE Website, 2019\. https://sites.ualberta.ca/~rabelo/COLIEE2019/.
  * [21] S. Seabold and J. Perktold, Statsmodels: Econometric and statistical modeling with python, in: 9th Python in Science Conference, 2010\.
  * [22] M. Yoshioka, Y. Kano, N. Kiyota and K. Satoh, ‘Overview of Japanese Statute Law Retrieval and Entailment Task at COLIEE-2018,’, in: Twelfth International Workshop on Juris-informatics (JURISIN 2018), 2018\.
  * [23] Proceedings of the Seventeenth International Conference on Artificial Intelligence and Law, ICAIL 2019, Montreal, QC, Canada, June 17-21, 2019, ACM, 2019, ICAIL. ISBN 978-1-4503-6754-7. doi:10.1145/3322640.