Exact Solution for Relativistic Trajectories Using Modal Transseries

: In this article, we design a novel method for ﬁnding the exact solution of the geodesic equation in Schwarzschild spacetime, which represents the trajectories of the particles. This is a fundamental problem in astrophysics and astrodynamics if we want to incorporate relativistic effects in high precision calculations. Here, we show that exact analytical expressions can be given, in terms of modal transseries for the spiral orbits as they approach the limit cycles given by the two circular orbits that appear for each angular momentum value. The solution is expressed in terms of transseries generated by transmonomials of the form e − n θ , n = 1, 2, . . ., where θ is the angle measured in the orbital plane. Examples are presented that verify the effect of the solutions.


Introduction
General Relativity is one of the most fundamental and successful theories of physics [1]. It predicts both the large-scale expansion of the Universe [2], the behavior of small gyroscopes in orbit around the Earth [3], the so-called de Sitter and Lense-Thirring precessions, and the propagation of gravitational waves from distant black hole sources, among many other phenomena [4]. However, the field and motion equations of the theory are nonlinear, and many of these areas require the implementation of numerical methods to yield results capable of being tested by experiments and observations [5]. This is a handicap not only of General Relativity, but also of many physical theories whose equations are also nonlinear: Navier-Stokes equations in hydrodynamics [6], Yang-Mills equations in quantum field theory [7], or the Korteweg-de Vries equation for shallow waves [8], only to name a few. In mathematical finance, we have, for example, the Black-Scholes equation [9], or the Fitzhugh-Nagumo model for neurons in biology [10,11]. Therefore, it needs no further discussion to realize that nonlinear differential equations are an indispensable tool in our knowledge of the world. Numerical methods are not always reliable, especially when we are interested in the asymptotic behavior of the solutions, we are looking for very accurate predictions, or we are simply interested in the classification of the different behavior predicted by the equations or the existence and uniqueness of solutions under certain conditions. These are situations in which an analytical solution certainly has an advantage.
The interest in analytical and semi-analytical methods has increased in the last few decades with the development of the δ-expansion and Adomian's decomposition methods, as well as the homotopy method, which includes the other as particular cases [12]. The homotopy method is based on a continuous mapping of the solution to a function Φ(t; q) depending on an embedding parameter, q, ranging from zero to one in such a way that for q = 1, we get the solution of the nonlinear problem. The nonlinear equation is then mapped into a family of equations involving a linear operator, i.e., a set of deformation equations that correspond to the original equation for q = 1. However, the solution series involve derivatives of the auxiliary function Φ(t; q) that may become very complicated, and resorting to Padé's approximants is usual for many problems to accelerate convergence.
In the domain of linear differential equations, a special role is played by Fourier series and Laurent series for which a general theory slowly emerged in the works of Taylor, Fourier, Cauchy, and Laurent in the XVIIIth and XIXth Centuries, and it allowed the development of an analytical approach for many equations appearing in the sciences, which now carry the name of the mathematicians that investigated them: the Bessel, Legendre, Laguerre, Hermite, or Mathieu equations [13,14].
In recent times, a generalization of the theory of Laurent series has been formalized by several authors following a work that dates back to Levi-Civita and Hahn more than a century ago and that was extended later by the model theorists Dahn and Göring [15] in connection with Tarski's exponential function problem and by Écalle in his proof of the Dulac conjecture on plane analytic vector fields [16]. More recent reviews were given by Edgar [17] and in books by van der Hoeven [18] and Aschenbrenner et al. [19]. These authors approached the topic from the point of view of model theory, and they proved model completeness among many other results.
In this paper, we will not need all this machinery, but only a subfield of the field of transseries known as modal transseries. This can be defined as a trascendental extension of the field of reals, is a continuous and differentiable function [20]. This structure has proven very useful in the solution of nonlinear differential equations with polynomial nonlinearities (appearing as powers of the function and/or its derivatives). These techniques are in the spirit of traditional analysis of linear differential equations, and they provide solutions in closed analytical form with a series of coefficients that can be obtained by recurrence. For example, by using the modal transseries method, we solved the SIRmodel of epidemiology [21,22], the Michaelis-Menten equation for enzyme kinetics [23], the Lorenz system in the laminar regime [24], and Einstein's field equations for planar gravitational waves [4].
These efforts, unsystematic to the present day, show a path towards a general technique for finding analytical solutions of nonlinear differential and integral equations with polynomial nonlinearities. The adequate context for them is the real analysis of transseries, and in Appendix A, we give some formal definitions that could help in the formalization of the modal series method.
In particular, in this paper, we will discuss the orbital equation for a particle in Schwarzschild spacetime, and it will be shown that e −θ is a resolvent generator (see Appendix A) for it, providing a method for an analytical discussion of spiral orbits that, up to now, have only been studied numerically [25,26].
There are techniques to tackle this classic problem. For example, Friedman and Steiner [27] used perturbation methods to analyze the precession of Mercury's perihelion at a high order of correction. Exact solutions are also available using Jacobi's elliptical functions [28]. However, we consider that this problem is ideal for testing the modal series method, which has also been fruitful in other areas of research since it is computationally very simple to program. Spiral orbits are a special solution of the orbital equation in general relativity not found in Newtonian mechanics. A curious anecdote is the letter written by Isaac Newton and addressed to Robert Hooke in which a drawing of a spiral orbit followed by a projectile towards the center of the Earth appears. Later on, Hooke bitterly critiqued Newton on this subject spurred by his unbounded claiming of his priority in the discovery of Universal Gravitation.
Newton replied that the spiral has only been a "negligent stroke whit his pen". Surprisingly, spiral orbits recuperated in the General Theory of Gravitation, which superseded Newton's Theory of Gravity as action at a distance.
In contrast with the amplitudes of Newtonian orbits in the two-body problem, solutions for the orbits in General Relativity are often treated numerically or implicitly.
Einstein's equations can also be solved on a computer using sophisticated numerical methods. Given sufficient computer power, such solutions can be more accurate than post-Newtonian solutions. However, such calculations are demanding because the equations must generally be solved in a four-dimensional space. Nevertheless, beginning in the late 1990s, it became possible to solve difficult problems such as the merger of two black holes, which is a very difficult version of the Kepler problem in general relativity.
This paper is organized as follows. In Section 2, we present the equation of study. Section 3 is devoted to the proposed methodology and its numerical analysis. In Section 4, scattering and bound precessing orbits are presented. In Section 5, the simulation results and applications are shown, and the last section is devoted to the discussions. Finally, some formal definitions are given in Appendix A.

Equation of the Orbit
In this paper, we are interested in studying the exact solutions of the nonlinear differential equation for the orbit of a test particle in Schwarzschild spacetime. For the sake of clarity, the ordinary differential equation is presented in the following form, where u = 1 r and r is the inverse of the radial coordinate towards the center of the large body.
The derivative is with respect to the angle θ, β = M h 2 (h being the angular momentum, M the black hole mass in terms of distance), and α = 3M. Both M and h are normalized to have units of length; see [29]. Multiplying both sides of the equation by Schwarzschild's radius r S = 2 G M/c 2 and redefining u = r S /r, we obtain an equation for an adimensional variable with adimensional parameters.

Proposed Methodology
Following a technique proposed in [21,22], we consider a series solution of the form: where w and the coefficients u n will be determined after. Therefore, the second order derivative in (2) is then given by: Now, on the other hand, the expression u 2 (θ) can be rearranged as a Cauchy product, i.e., the discrete convolution of the two is: By direct substitution of (2), (3), and (4) into Equation (1), we obtain the following expression: In the expression (5), we impose the condition that: and that Therefore, from (6), u 0 is given by: We know that u 0 is real if 4αβ > 1. Then, 12M 2 > h 2 . Thus, there are two possible circular orbits, for each angular momentum h > 2 √ 3 m. Now, for h = 2 √ 3 m, i.e., 4αβ = 1, we find the smallest possible circular orbit; that is, On the other hand, from (7) and using the technique proposed in [23], we have that: Thus, for n = 1: If we choose u 1 = 0, then it follows that w = √ 2αu 0 − 1. For the minus sign in Equation (2), we have that: Thus, for n = 2, 3, 4, ..., from (9), it follows that:

Analysis of Convergence
In this subsection, we analyze the convergence of series given by (2), where the coefficients u n are calculated by (12).
We put U 1 = |u 1 | and U k = |u k | for k ≥ 2. Therefore, we obtain the following expression to (12): Thus, one gets the recurrence relation: For k = 2 in (13), Therefore, Now, for k = 3 in (13), Thus, We can continue, and finally, it follows that: Let C = α U 1 w 2 , and we put Equation (16) as: the iterating (17), we have that: Thus, using the method of mathematical induction, one gets that: i.e., If we choose U 1 = |u 1 | ≤ 3w 2 α , then U k is a monotone decreasing sequence, and the series ∞ ∑ k=1 (−1) k U k converges absolutely. The above analysis allows us to state the following theorem: converges absolutely.

Scattering and Bound Precessing Orbits
It is clear that the series in Equation (2) corresponds to a spiral orbit. Obviously, a series of exponentials cannot describe other kinds of orbits to be found around a star or black hole [26]. However, we can look for another resolvent generator in consonance with Definition A3. Sine powers can also play this role, and consequently, we propose another family of solutions in the form: We will show that this series encompasses two different families of orbits: (i) the unbounded scattering orbits and (ii) the bounded precessing orbits.

Scattering Orbits
In this kind of orbit, the particle comes from infinity and goes to infinity. As u is the reciprocal of the distance to the mass center, we should choose u 0 = 0. For θ = 0 and θ = π/w, we have then that u(0) = u(π/w) = 0, and the orbit is unbounded. Another physical condition comes from the fact that the incoming and outgoing radial velocities must be finite (indeed, they should be strictly smaller than the speed of light), and so: because lim θ→0 r(θ) = ∞. This implies that u 1 = 0. By substitution of Equation (21) into the equation of the orbit in (1), we also find that: and the recurrence relation: for n ≥ 2. In this case, the parameter w can be chosen freely, and it determines the number of loops the particle performs around the black hole before escaping again to infinity. An example of an orbit of this kind will be given in Section 5.

Bounded Precessing Orbits
By using the same series as given in Equation (21), but taking the following conditions on the coefficients: u 2 = u 3 = 0, u 1 = 0, we obtain a bounded orbit. This is the consequence of a nonzero u 0 value obtained as the solution of: The physical solution corresponding to a precessing orbit at a large distance from the black hole will be obtained by taking the minus sign in Equation (8): The first recurrence for the coefficients with an odd index is: and this gives us the value of the frequency w: The coefficients u n , n = 4, 5, . . . are then found by using the general recurrence relation in Equation (23).
A special case corresponds to α = 0, i.e., to Newtonian mechanics. In this particular case, we obtain u 0 = β, ω = 1, u n = 0 for n ≥ 2, and u 1 is a free parameter. In celestial mechanics, we usually define u 1 = β, yielding: This is the equation for a Keplerian elliptic orbit, being the orbital eccentricity.

Raabe-Duhamel's Convergence Test
In this case, the study of convergence is more complicated than for the exponential series as discussed in Section 3.1. A general expression for the series coefficients cannot be obtained, but we can check numerically that lim n→∞ u n+2 /u n = 1; consequently, we cannot prove convergence by Cauchy's test.
Some insight is provided by the particular series given by Equation (23) with α = 0: Therefore, we have two independent series for the odd and even subscripts, and each of them can be proven to be convergent for any real value of ω. Firstly, if √ 1/ω is an integer, the number of nonzero terms defined by Equation (29) is finite. For other values of ω, we can write: and we conclude that the series is convergent by Raabe-Duhamel's test [27]. From Equation (30), we also have: u n /u n+2 = 1 + 3/n + O(n −2 ). By assuming that this expression can be generalized for the general case with α = 0, we can prove convergence for small values of α. (23) for u 1 = 0, u 0 = 0 converges absolutely if u n /u n+2 = 1 + f (α)/n + O(n −2 ) for n 1, f (α) being a positive real value that depends on α and α |u 0 u 2 | e f (α)/2 < 4 ω.

Theorem 2. The series defined by the recurrence relation in Equation
Proof. The idea is to find a bound for the Cauchy product: which for n even can be written as: C n = 2u 0 u n + 2u 2 u n−2 + . . . + u n/2 u n/2 , and a similar expression for n odd with u (n−1)/2 = 0, u (n+1)/2 for the last term. We will assume n to be even (the proof runs in parallel for n odd). We propose the hypothesis that u n is a decreasing series, and consequently: By the triangular inequality, we constrain the absolute value of C n as follows: where K is a constant and [ ] denotes the integer part. Now, from Equation (23), we have: where we used the bound in Equation (34) and the assumption that all u n , n = 0, 1, . . . have the same sign. Consequently, Raabe-Duhamel's test implies convergence for any value of α sufficiently small to verify: This gives us a very stringent condition for convergence, and in practice, it is found numerically that values of α larger than those provided by Equation (34) also provide convergent solutions.

Simulation Results and Applications
In this section, we will discuss some examples for orbits calculated with the modal series obtained in this paper. We will show that modal series have an advantage over other numerical methods for calculating high precision trajectories in the two-body problem for General Relativity. In Figure 1, we show two examples of spiral orbits calculated by summing 1000 terms in Equation (2). Both cases correspond to α = 1/5 and β = 1, so the radius for the inner unstable circular orbit is given by: We studied two spiral orbits by choosing different values of u 1 and taking u 0 = 1/r in . For u 1 = −5, the inner orbit in Figure 1 is found. This orbit starts for θ = 0 at a distance r 1.03686 from the center of the massive body. Similarly, for u 1 = −10, the outer trajectory plotted in the polar graph starts at r 2.98215 for θ = 0 (this point is outside the plot of Figure 1). The convergence of the series in Equation (2) is very fast for both cases, and the coefficient u 1000 is of order 10 −124 for the second case and even smaller for the first case. This means that a very accurate calculation of the spiral orbit is possible with a relatively small number of coefficients and that precision will be limited by the accuracy of the physical parameters of the system of interest and finite arithmetic instead of the truncation of the series.
Scattering orbits were studied in Section 4.1. They correspond to a particle with sufficient energy to be unbound from the central body, so it may come from a large distance, perform a flyby with one or several loops around the main body, and escape again to infinity. In classical Newtonian physics, only hyperbolic and parabolic orbits are found, but General Relativity shows a richer casuistry.
As an example, we chose the parameters β = 1/100, α = 1/5, w = 2/5, u 0 = u 1 = 0, and u 2 = β/(2w 2 ) = 1/32 as fixed by Equation (22). The resulting one-loop orbit is shown in Figure 2. The reason that the orbit exhibits a single loop is a consequence of the chosen value of w: Notice that for θ = 0, we have u = 0 (see Equation (20)), and this means that the particle comes from infinity. The same value is found again for θ = π/w = 5π/2 = 2π + π/2. Therefore, the particle escapes at a right angle from the initial direction after performing a loop around the black hole. Convergence is not so fast for scattering orbits as in the case of spiral ones because the coefficients' ratio tends to unity for large n. We have that u 1000 −4.4386 × 10 −7 , and double precision is not achieved with this number of terms. Anyway, uncertainty in the masses of the bodies or their distances is indeed a most important limiting factor for increasing the accuracy of astronomical observations. Finally, we considered a bounded precessing orbit for α = 0.05, β = 0.1, u 1 = 0.1, u 0 = 0.100505, and w 0.994962. The values of w and u 0 are fixed by Equations (11) and (25), respectively. Remember that u 2 = u 3 = 0 in this case, and the rest of the coefficients are obtained by the recurrence relation in Equation (23). The result is a precessing elliptic orbit with an advance of the perihelion per revolution given by: corresponding to 1.82284 degrees in the case of Figure 3. The anomalous advance of the perihelion of Mercury was already known by astronomers since the XIXth Century [28], and this was the first prediction of General Relativity that convinced Einstein himself of its validity.
The result obtained from approximate perturbation methods gives the well-known formula: where the gravitational constant, G, and the speed of light, c, have been explicitly restored and we have taken into account that β = M/h 2 = 1/(a(1 − e 2 )) is the reciprocal of the semi-latus rectum [29]. The accepted values for the parameters in Equation (37) for Mercury's orbit [30] are G = 6.67408 ( yielding the value of 42.98118333 arcsec per century. This corresponds to 0.00858 milliarcseconds per century more than the prediction of the approximate formula usually given in textbooks [28]. We should notice that, by using elliptic functions, Saca obtained an improved formula for the perihelion advance, which leads to similar results to ours [31], but his approach is more cumbersome.

Conclusions
Nonlinear differential equations have been known and used for a long time. Perhaps the canonical example is that of the three-body problem in celestial mechanics necessary to explain the motion of the Moon under the mutual influence of the Earth and the Sun. With this idea in mind, Euler developed his famous numerical method widely used since [32]. The difficulties in solving these equations have promoted the preponderance of numerical methods throughout the history of the subject. Particular analytical solutions have been found for some problems, but they are restricted to special cases or convergence is very slow. In the case of the three-body and N-body problems, the series solutions by Sundman [33] and Wang [34] provide almost a complete analytical treatment of these problems, but their practical application is very limited because we will need to add millions of terms to obtain predictions even for short periods of time. These series are given in powers of t 2/3 , and they are an example of transseries.
However, solutions in terms of transseries can be useful in many problems of interest for physics or engineering, and moreover, they can be studied systematically. In this paper, we provided such solutions for the orbits in a Schwarzschild spacetime. These solutions are found in terms of powers of exponential or sine functions showing that different functions can be used as a basis to solve the same problem depending on the initial conditions. The existence of a resolvent generator as defined in (A2) gives the chance to find a modal transseries solution, but we have shown that several series can be found for the same equation and can describe different classes of trajectories. Another advantage of transseries over traditional numerical or approximate methods is the calculation of an exact expression for the advance of the perihelion in bounded trajectories as given in Equation (38). This expression is usually approximated by Equation (37), but this yields the wrong predictions for very massive objects such as neutron star binaries or black holes.
Apart from the problems discussed at the Introduction, we also successfully applied the modal transseries method to another problem in General Relativity: gravitational waves [4], for which explicit solutions corresponding to plane waves are given. This proves the fecundity of the transseries approach to nonlinear differential equations.
We hope that these works stimulate the research in modal transseries solutions for nonlinear differential equations of interest in physics, chemistry, or engineering by providing an alternative to traditional numerical methods in finite differences, which can be used to back the numerical computations or to provide accurate results and a deeper analytical understanding of the solutions and behavior of these systems.
if and only if L 1 > L 2 as x → ∞ or L 1 = L 2 and b 1 > b 2 . As an example, we have [17].
The transseries T[G] forms a differential field under the usual operations of summation and derivation term by term and multiplication. For example: In this paper, we consider a particular subfield of this rich structure that we call the modal transseries subfield: Definition A1 (Modal transseries:). Given a transmonomial f[x], we define a modal transseries as the sum: The set of all modal transseries of the form in Equation (A2) is denoted by T f , and we say that it is generated by f[x] or that f[x] is the generator of the modal transseries. Notice that for f[x] → 0 as x → ∞, the terms in Equation (A2) are in descending order.
As stated above, we have also that: Lemma A1. The set of modal transseries generated by f[x], T f , is a differential subfield of T.
Proof. T f is isomorphic to a trascendental extension of R, and its elements can be differentiated term by term.
Modal transseries have proven useful in the solution of many nonlinear differential equations. In the particular case studied in this paper, the recurrence relations found for the exponential set of modal functions are given in Equations (5) and (23). We will also discuss the results for other problems in order to clarify the way in which the method is applied in a general case. In [4], for the case of plane gravitational waves, the analytical solution of the differential equation corresponding to a type of Sturm-Liouville problem:p (u) + V(u)p(u) = 0, is given by where is a constant, p(u) is a metric function, and V(u) is the fourth Weyl scalar Ψ4; a n holds the following recurrence relation: a n+1 = 2 (n 2 − n + 2) (n + 1)(2n + 1) a n−1 − 3 + ( − 2)n + 2( − 1)n 2 (n + 1)(2n + 3)a n , with Now, in [21], an epidemiological model of the SIRStype, which is a system of three first-order differential equations, is transformed into a first-order integral differential equation as a function of the variable I(t) that represents infectious individuals: where α, β, ν, R 0 are constants given and a type of solution is proposed in the form of transseries A k e −kωt , ω > 0 and with the following recurrence formula: with a = α − β = ν + µ − β, b = βR 0 , c = β, d = βν, ς = a + (2c + d δ )A 0 and ω satisfies ω 2 − (δ + ς)ω + δς + A 0 .
The set of modal functions necessary in a particular problem has been found by trial and error. There is no method, for the moment, to choose the adequate set for a general nonlinear differential equation. Now, on the other hand, there are several ways to combine both real and complex exponentials; for example in [24], the Lorenz system was solved by transforming it into a single integro-differential equation for x(t) given by: Thus, for the point relaxation towards the origin, it proposes a modal expansion series for x(t) in: whereas that for relaxation towards a non-trivial equilibrium point and the case when ω is a complex number, the solution for x(t) is given by: A j,k e −jωt−kω * t , with A j,k = A * k,j , (ω + ω * ) ∈ R. In [23,24], in those cases, the generator and any combination of its powers satisfy that under the application of the nonlinear operator, it transforms into modal transseries of the same kind. Therefore, to generalize these concepts, we define: Definition A2 (Nonlinear polynomial differential operator:). A nonlinear polynomial differential operator O[y [x]] acting on single one-variable functions y[x] is any finite multivariate polynomial in y[x] and its derivatives. Nonlinear polynomial differential equations in one variable are given by: An example is the trajectory equation for a particle in the spacetime around a Schwarzschild black hole we solve in this paper. A minimum condition for a modal transseries to be useful in the solution of nonlinear differential equations is the existence of a resolvent generator for the corresponding nonlinear operator defined as follows: Definition A3 (Resolvent generator). A transmonomial f[x] is called a resolvent generator for the nonlinear polynomial differential operator O if for any modal transseries, t ∈ T f , we have that O[t] ∈ T f , i.e., if O is an invariant differential operator for the elements of T f .