New Applications of Elliptic Functions and Integrals in GPS Inter-Satellite Communications with Account of General Relativity Theory

1. Introduction

In the past $10-15$ years the problem about GPS (GPS means Global Positioning System) satellite-ground station communications has been replaced by the problem about autonomous navigation and inter-satellite communications (ISC) (links), which has been mentioned yet in 2005 in the monograph [1]. Autonomous navigation means that satellites should have the capability to transmit data between them via inter-satellite cross-link ranging [2] and thus, to ensure navigation control and data processing without commands from Earth stations in the course of six months. It is important that next generation space missions will attain sub-millimeter precision of measuring distances beyond ${10}^{-6}$ meters (1 $micrometer=1$$\mu m={10}^{-6}$$m)$ by means of ultrashort femtosecond pulse lasers. The theoretical description of such measurements is inevitably related with General Relativity Theory (GRT). This means that the relative ranging and relative velocity model should account also for the bending of the transmitting path of the signal, which is significant for such large distances between the satellites due to the action of the gravitational field. The bending of the light (laser) is related to the important and fundamental physical fact that due to the action of gravitational field, the signal travels a greater distance in comparison with the straight path distance in the absence of a gravitational field.

In the book [8] autonomous navigation is defined as referring to “processes in which the spacecraft without the support from the ground-based TT&C (Telemetry & Tracking Control) system for a long time, relying on its onboard devices, obtains all kinds of measurement data; determines the navigation parameters like orbit, time and attitude”. In this monograph DORIS and PRARE navigation systems are determined as non-autonomous, because “DORIS system can determine the spacecraft’s orbit with high accuracy, but both need to exchange information with the ground stations”. Thus, the main aim of autonomous navigation is to reduce the dependence of spacecrafts on the ground TT&C network and in this way to enhance the capability of the systems anti-jam and autonomous survivability. It is evident also that the determination of the satellite’s exact location is achieved through the reception, processing, and transmitting of ranging signals between the different satellites.

Autonomous GPS navigation system is necessary in view of the proposal of researchers from the University of Texas, Aerospace Corporation, National Bureau of Standards, International Business Machines Corporation (IBM) and Rockwell Automation Inc. to monitor nuclear explosions, based on GPS inter-satellite communication link.

The development of inter-satellite laser communication systems in space - orbit technology that enable super-high-speed data transfers at rates greater than 1 $Gbps$ is widely applied also in.cube-/nanosatellite platforms such as CubeLCT, AeroCube-7B/C, CLICK, LINCS-A/B, SOCRATES and LaserCube [9]. However, in order to establish laser communications, high performance of arc-second level pointing system is required, and this is a difficulty for the nanosatellite platforms. It can be supposed why there is such a difficulty - since the data-transmission rate of $1Gbps$ between the two nanosatellites is at an inter-satellite range of 1000 $km$, the effects of curving the trajectory of the laser signal may be considerable. In publications, related to small satellite optical links [10] it has been admitted that the s.c. "pointing error" arizes not only due to tracking sensors and mechanical vibrations, but also due to the base motion of the satellite. This fact, together with the ranging of signals and the pointing error illustrate the assumption that not only the motion of two satellites is important, but the effect of curving of the signal trajectory. The paper [10] also asserts that for a GPS constellation with 24 satellites, there are a total of $8-16$ links - they can be forward and backward links in the same orbit, but as well as lateral links between adjacent orbits. Moreover, the distance of GPS inter-satellite crosslink can reach 49 $465km$.

The establishment of inter-satellite links (ISL) and measurement communications is of primary importance for the relative ranging and relative velocity determination between the satellites from one constellation or from different constellations such as the GPS system, the Russian GLONASS, the European Galileo and the China Beidou second-generation system, all of them considered to be interoperable with each other.

The theory of inter-satellite communications (ISC) is developed in the series of papers by S. Turyshev, V. Toth, M. Sazhin [3,4] and S. Turyshev, N. Yu, V. Toth [5]. The theory in these three paper concerns the space missions GRAIL (Gravity Recovery and Interior Laboratory) [3], GRACE- FOLLOW-ON (GRACE-FO - Gravity Recovery and Climate Experiment - Follow On) mission [4] and the Atomic Clock Ensemble in Space (ACES) experiment [5,6,7] on the International Space Station (ISS). For the second model in [4], the relativistic change of the phase of the signal has been computed due to the changing (geodesic!) distance between the emitting satellite and the receiving satellite. The geodesic distance is the real distance, travelled by the light or laser distance and due to the relativistic effect of the gravitational field to curve the trajectory of the signal, it is greater than the Euclidean distance. Because of this effect, in the papers [3,4] and [5] the difference between the emitting time and the receiving time is determined to be a relativistic observable, dependent on the geodesic distance, travelled by light. In such a case this time difference contains also the relativistic logarithmic correction, typical for the Shapiro delay formulae. This formulae and its application in the current research shall be outlined further.

It should be stressed also that the above mentioned missions are realized by means of low-orbit satellites (about 450 $km$ above the Earth) and the distance between the satellites is about 220 $km$, while the satellites from the GPS, GLONASS and the Galileo constellations are on a much higher orbit (from 19140 $km$ for GLONASS and ranging to 26560 $km$ for GPS, even higher). Consequently, the distance between the satellites from one or from different constellations cannot be determined with the accuracy (1 micrometer), typical for the GRAIL and GRACE missions, accounting also of the relative motion between the satellites.

2. A Brief Review of the Shapiro Delay Formulae and the Aims of This Paper

2.1. Standard Shapiro Delay Formulae, Derived from the Null Cone Equation

The basic theoretical approach to find the propagation time of signals in the near-Earth space with account also of the gravitational field and in the framework of General Relativity Theory is based on the s.c. Shapiro delay formulae, proposed yet in 1964 [12] and frequently used also in VLBI (Very Large Baseline Interferometry) (see the review article by Sovers, Fanselow, Jacobs [13])

$$T_{AB}={\displaystyle \frac{R_{AB}}{c}}+{\displaystyle \frac{2G{M}_E}{c^3}}ln\left({\displaystyle \frac{r_A+r_B+R_{AB}}{r_A+r_B-R_{AB}}}\right),$$

(1)

where the coordinates of the emitting and of the receiving satellite are correspondingly $\mid x_A\left(t_A\right)\mid =r_A$ and $\mid x_B\left(t_B\right)\mid =r_B$, and $R_{AB}=$$\mid x_A\left(t_A\right)-x_B\left(t_B\right)\mid$ is the Euclidean distance between the signal - emitting satellite and the signal - receiving satellite. Note that the space points of emission $x_A\left(t_A\right)$ and arrival $x_B\left(t_B\right)$ of the signal are taken at the time $t_A$ of emission and of arrival $t_B$ correspondingly. In (1) $G_{\oplus }{M}_{\oplus }$ is the geocentric gravitational constant, $M_{\oplus }$ is the Earth mass and $G_{\oplus }$ is the gravitational constant, $t=$TCG denotes the s.c. Geocentric Coordinate Time (TCG). Further (and also in the above formulae) we shall denote by $T_{AB}=T_B-T_A$ the signal propagation time and thus we have replaced the times $t_A$ and $t_B$ by $T_A$ and $T_B$, so that to underline that these are times of emission and reception. The second term in formulae Eq. (1) is the Shapiro time delay term, accounting for the signal delay due to the curved space-time.

In order to calculate the propagation time $T_{AB}$ of the signal, it should be remembered that according to General Relativity (see the book [14]) the emitter of the first satellite emits a signal propagating on the gravitational null cone after setting up to zero the infinitesimal metric element $d{s}^2=0$

$$d{s}^2=0=g_{00}{c}^{2}d{T}^2+2g_{oj}cdTd{x}^j+g_{ij}d{x}^{i}d{x}^j.$$

(2)

For the null-cone equation (2), the solution of this quadratic algebraic equation with respect to the differential $dT$ can be given as [16],

$$dT=\pm {\displaystyle \frac{1}{c}}{\displaystyle \frac{1}{\sqrt[.]{-g_{00}}}}\sqrt[.]{\left(g_{ij}+{\displaystyle \frac{g_{0i}{g}_{0j}}{-g_{00}}}\right)d{x}^{i}d{x}^j}+$$

$$+{\displaystyle \frac{1}{c}}\left({\displaystyle \frac{g_{0j}}{-g_{00}}}\right)d{x}^j,$$

(3)

where the metric tensor components are determined for an Earth Reference System with an origin at the centre of the Earth.

Further, the propagation time should be found after 1. choosing the metric tensor components $g_{ij}$, $g_{0i}$ and $g_{00}$ in the given expression for the zero infinitesimal metric $d{s}^2=0$ (2) and 2. choosing the coordinates $x^i$ ($i,j=1,2,3$) in the metric, which additionally may depend on one, two or more parameters. Note that we have assumed that $x^0=1$. In the present investigation, the coordinates shall be chosen to be related with the satellites motion along the elliptic orbit. Usually, geostationary satellites are situated on nearly circular orbits. Satellites from the GPS constellation have a very small eccentricity of the orbit of the order of $e=0.01$, for the Russian constellation GLONASS the eccentricity is $e\approx 0.02$. Let us remind that if a and b are respectively the large and the small semi-axis of the ellipse, then by definition $e={\displaystyle \frac{\sqrt{a^2-b^2}}{a}}={\displaystyle \frac{c}{a}}$ where c is the distance between the foci of the ellipse. Thus the eccentricity measures to what extent the ellipse deviates from the circle. Smaller eccentricity will mean an ellipse, which is close to the circle, while high eccentricity will be related to strongly elongated ellipses. For elliptic orbits $0\le e<1$, for $e\ge 1$ the orbit will be a hyperbolical one and will not be a closed trajectory - these are trajectories of satellites and celestial bodies, leaving the Solar system. .

2.2. The Purpose of This Paper and the Relation to the Parametrization of the Coordinates in the Metric Tensor

The main goal of this paper will be the exact analytical calculation of the propagation time of a signal by using a solution of the null cone equation (2) and the parametrization of the space coordinates of the satellite for two different cases.

2.2.1. Case A - Plane Elliptical Coordinates of the Satellite

This is the most simple case, when the two-dimensional $(x,y)$ coordinates of the plane elliptical orbit of the satellite is parametrized in terms of the large semi-major axis a of the ellipse, the eccentricity eand the eccentric anomaly angle E, which is a dynamical parameter, related to the motion of the satellite on the orbit. The analytical expression for the parametrization of $(x,y)$ is

$$x=a(cosE-e),y=a\sqrt[.]{1-e^2}sinE$$

(4)

and it is well-known in all basic courses on celestial mechanics see especially [24,25].The time coordinate is chosen to be the celestial time $t_{cel}$, which is related to the eccentric anomaly through the Kepler equation

$$E-esinE=n(t_{cel}-t_p)=M.$$

(5)

In this equation e is the eccentricity of the orbit, $n=\sqrt[.]{{\displaystyle \frac{GM}{a^3}}}$ is the mean motion, M is the mean anomaly and $t_p$ is the time of perigee passage for the satellite. The mean anomaly M is an angular variable, which increases uniformly with time and changes by ${360}^0$ during one revolution. The mean motion n has the following geometrical meaning: this is the motion of the satellite along an elliptical orbit, projected onto an uniform motion along a circle with a radius equal to the large axis of the ellipse. Usually M is defined with respect to some reference time - this is the time $t_p$ of perigee passage, where the perigee is the point of minimal distance from the foci of the ellipse (the Earth is presumed to be at the foci). Since the eccentric anomaly angle E will play an important role in the further calculation, let us also remind how this notion is defined: if from a point on the ellipse a perpendicular is drawn towards the large axis of the ellipse, then this perpendicular intersects a circle with a centre O at a point P. Then the eccentric anomaly represents the angle between the joining line $OP$ and the semi-major axis (the line of perigee passage).

In [28] E is determined as an auxiliary angular variable such that $a-r=aecosE$, which has the following geometrical meaning with respect to the ellipse: if $r_p$ and $r_a$ are the corresponding radius-vectors at the perigee passage and at the apogee, then $E=0$ for $r=r_p$ and $E=\pi$ for $r=r_a$.

The first main purpose of this paper will be to solve the equation (3) for $dT$ as differential equation with respect to the eccentric anomaly angle E , for the given parametrization (4) of the coordinates $(x,y)$ and for the null-cone metric of the near-Earth in the standard form

$$d{s}^2=-c^2\left(1+{\displaystyle \frac{2V}{c^2}}\right){\left(dT\right)}^2+\left(1-{\displaystyle \frac{2V_{.}}{c^2}}\right)\left({\left(dx\right)}^2+{\left(dy\right)}^2+{\left(dz\right)}^2\right)=0,$$

(6)

where $V={\displaystyle \frac{G_{\oplus }{M}_{\oplus }}{r}}$ is the standard gravitational potential of the Earth. No account is taken of any harmonics due to the spherical form of the Earth since the GPS orbits are situated at a distance more than 20000 $km$ (considered from the centre of the Earth). In fact, the null-cone metric is the same, derived from the standard metric for the near-Earth space in basic monographs of gravity theory such as the one by Fock [29].

2.2.2. First Aim of the Paper - Expressing the Propagation Time in Terms of a Sum of Elliptic Integrals of Different Kinds (First, Second and Third) for the Case of a Signal Emission by a Satellite on a Plane Elliptical Orbit

For the first time in [18,20] and in the review paper [22], the solution for the propagation time T is found to be expressed in terms of elliptic integrals of zero-order, but of the first , second and of the third kind. Later it will become clear what is the difference between elliptic integrals of different orders and kinds. The imaginary unit appears in these elliptic integrals, so it is important to prove that the integrals are real-valued quantities, since the propagation time T is a real and not an imaginary quantity. The physical approximation which is used for finding the solution is $\beta ={\displaystyle \frac{2V}{c^2}}\ll 1$, which physically is related to the assumption about a weak gravity field. This assumption is consistent with the experimental data and the calculation of the numerical quantity $\beta =$${\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^{2}a}}=0.334.{10}^{-9}\ll 1$ by some authors [31] and also with the relativistic effects of the order of $c^{-2}$(further in the text this accuracy will be denoted as $O\left(c^{-2}\right)$). This is a consequence of the experimental fact that a GPS constellation has to keep time to an accuracy of about four nanoseconds per day (4 $nsec/d=4\times {10}^{-9}$$sec/d$) [31], which is compatible with the fractional time stability of $5\times {10}^{-14}$, maintained by atomic clocks. Naturally, the fractional stability is constantly improving in the past and present years with the perspective of becoming smaller and smaller.

Accuracy of the Relativistic Effects and the Possibility to Find the Solution of the Null Cone Equation

The last fact is mentioned due to the following interesting consistency between the experimentally measured accuracy of relativistic effects and the found solution of the null cone equation (6) for the propagation time T in terms of the eccentric anomaly angle E As it will be shown, this solution is possible to be found only under the assumption $\beta ={\displaystyle \frac{2V}{c^2}}\ll 1$. The necessity to implement such an assumption, physically expressing the fact about the weak gravitational field, is related also with the fulfillment of the Equivalence Principle only for weak gravitational fields and slow motions - an important topic, discussed in the monograph by Fock [29].

2.2.3. Second Aim of the Paper - Theoretical Justification of the Approach from the Viewpoint of General Relativity Theory and Differential Geometry

The second main purpose of the paper is to justify the theoretical approach and to provide the reasoning for choosing as parametrization variables in the null cone equation (6) the coordinates (4), parametrizing the satellites trajectory along the orbit in terms only of three variables - $a,e$ and E, among which only the eccentric anomaly angle E is the dynamical parameter. The need to justify the theoretical approach will refer not only for the simplified case of the Kepler coordinates (4), but also for the next (second) case of a space-oriented orbit, when the disposition of the orbit in space will be described by 6 Kepler parameters (they will be determined in the next section). From a conceptual point of view, this is a very non-trivial problem, since the presence of the Kepler coordinates (4) in the null cone equation (6) may raise the question whether it is correct to use coordinates, not related to the trajectory of the signal, but to the motion of the satellite. Following such a wrong logic, then it might turn out that the calculated time T is not a propagation time of a signal, but the time for motion of a satellite. Fortunately, this line of reasoning is not correct, because the essential fact is the null cone equation (6) and not the parametrization of the space coordinates $x,y$ or $x,y,z$. This will be confirmed by 1. A theorem in the monograph by Fock [29], which further will be formulated in details [22]. 2. A theorem from differential geometry [30] with a meaning, which is very close to the meaning in the monograph by Fock (cited also in the papers [22] and [23]).

In fact, here the problem can be stated in a different way, since the propagation time T will be found as a solution from the null cone equation (remember also that the Shapiro delay formulae (1) is also derived from the null cone equation), but it should be found from the null geodesics. So the question is - should the null geodesic equation be used (which, as known, is difficult to solve exactly) or one may use instead only the null cone equation? However, the null cone equation is compatible with the light-like geodesic equation in the sense that each solution of the null cone equation is a solution also of the light-like geodesic equation according to the theorem in the book by Fock [29], so it represents a geodesics. But then if it is a light-like geodesics, according to another theorem in [32], it represents the shortest path (length) from one point of the geodesics to another point on the same geodesics. This is consistent with the Fermat principle about the propagation of light, asserting that light travels between two points for a minimal time [33]. Since the light velocity is constant, light will travel along the shortest path between the two points . In the book by Mishtenko and Fomenko [30] the theorem is formulated in the sense that each extremal of the null cone functional

$$E\left(\gamma \left(s\right)\right)=\int \mid {\stackrel{.}{\gamma }}^2\left(s\right)\mid ds=\int g_{ij}\left(x\right){\stackrel{.}{x}}^i\left(s\right){\stackrel{.}{x}}^j\left(s\right)ds$$

(7)

is also an extremal of the length functional

$$L\left(\gamma \left(s\right)\right)=\int \sqrt[.]{g_{ij}\left(x\right){\stackrel{.}{x}}^i\left(s\right){\stackrel{.}{x}}^j\left(s\right)}ds,$$

(8)

but the statement in the inverse direction is not always fulfilled. This means that the extremals of the null-cone functional constitute of subset of the extremals of the length functional- see also the arguments in [22] and [23]. Length functional of the type (8) can be obtained also after the integration of the null-cone equation (6) for the case when the space coordinates $x,y,z$ are parametrized by several variables (see the next section) - in differential geometry well-studied are the s.c. first and second fundamental forms [30,32]. In subsequent papers this will also be studied in details, but this paper will focus only on finding the solutions for the propagation time in terms of elliptic integrals.

2.2.4. Third Aim of the Paper- Finding the Solution for the Propagation Time T for the Case of Space-Distributed Orbits in Terms of Elliptic Integrals of Higher Order (Second and Fourth)

Space-distributed (oriented) orbits are well-known from celestial mechanics (further we shall use the coordinate representation from the book [25]) and in this case the position of the orbit in space is characterized by 6 Kepler parameters $(f_{.},a_{.},e_{.},{\Omega }_{.},I_{.},{\omega }_{.})$, which are in fact the orbital coordinates. The analogue of the coordinate transformation (4) for this case of space-distributed orbits are

$$x={\displaystyle \frac{a(1-e^2)}{1+ecosf}}\left[cos\Omega cos(\omega +f)-sin\Omega sin(\omega +f)cosI\right],$$

(9)

$$y={\displaystyle \frac{a(1-e^2)}{1+ecosf}}\left[sin\Omega cos(\omega +f)+cos\Omega sin(\omega +f)cosI\right],$$

(10)

$$z={\displaystyle \frac{a(1-e^2)}{1+ecosf}}sin(\omega +f)sinI,$$

(11)

where $r={\displaystyle \frac{a(1-e^2)}{1+ecosf}}$ is the radius-vector in the orbital plane, the angle $\Omega$ of the longitude of the right ascension of the ascending node is the angle between the line of nodes and the direction to the vernal equinox, the argument of perigee (periapsis) $\omega$ is the angle within the orbital plane from the ascending node to the perigee in the direction of the satellite motion $(0\le \omega \le {360}^0)$, the angle I is the inclination of the orbit with respect to the equatorial plane.

Very often, instead of the true anomaly angle f, another variable is used - the argument of latitude $u=\omega +f$, being defined as the sum of the argument of perigee $\omega$ and the true anomaly f and geometrically representing the angle between the line of nodes and the position vector r. The argument of latitude will appear further in the calculation of the propagation time in terms of the celestial coordinates. A similar additive angular variable is the eccentric longitude $F=E+\omega +\Omega$ (equinoctial orbital characteristic), representing the sum of the eccentric anomaly angle E, the right ascension of the ascending node $\Omega$ and the argument of the perigee $\omega$.

From celestial mechanics it is also known that there is a relation between the eccentric anomaly angle E and the true anomaly angle f

$$tan{\displaystyle \frac{f}{2}}=\sqrt[.]{{\displaystyle \frac{1+e}{1-e}}}.tan{\displaystyle \frac{E}{2}}⟹f=2arctan\left[\sqrt[.]{{\displaystyle \frac{1+e}{1-e}}}.tan{\displaystyle \frac{E}{2}}\right].$$

(12)

Two Definitions of the True Anomaly Angle

If a satellite moves along an elliptical orbit, then the dynamical parameter, related to the motion of the satellite is the true anomaly angle f. The true anomaly fis defined as the geometric angle in the plane of the ellipse between periapsis (the closest approach to the central body) and the position of the orbiting satellite at any given time [26]. . The vector $\overrightarrow{r}$ has an initial point at the centre of the ellipse.

There is also another more "mathematical" definition, which makes use of the s.c. "Runge-Lenz" (or Laplace) vector $A_L$ [27], which lies in the orbital plane and thus is orthogonal to the angular-momentum vector $J=\mathbf{r}\times \stackrel{.}{\mathbf{r}}$ (called also areal velocity $\overline{h}=\mathbf{r}\times \stackrel{.}{\mathbf{r}}=const=J$)

$$\stackrel{.}{{\mathbf{A}}_L:=\mathbf{r}}\times \left(\mathbf{r}\times \stackrel{.}{\mathbf{r}}\right)-G_{\oplus }{M}_{\oplus }{\displaystyle \frac{\mathbf{r}}{r}}\stackrel{.}{=\mathbf{r}}\times \mathbf{J}-G_{\oplus }{M}_{\oplus }{\displaystyle \frac{\mathbf{r}}{r}}.$$

(13)

Thus, the true anomaly f is the angle between the Runge-Lenz vector ${\mathbf{A}}_L$ and the position vector $\mathbf{r}$. The vector $A_L$ appears as an additive constant after integrating the equation

$$\mathbf{h}\times \stackrel{..}{\mathbf{r}}=-G_{\oplus }{M}_{\oplus }{\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{\mathbf{r}}{r}}\right),$$

(14)

and $\Delta A={\displaystyle \frac{1}{2}}\mid \mathbf{r}\times \stackrel{.}{\mathbf{r}}\Delta t\mid ={\displaystyle \frac{1}{2}}\mid \overline{h}\mid \Delta t$ is the area, swept by radius-vector r during the time $\Delta t$. The square of the Runge-Lentz-Laplace vector can be calculated to be [27]

$${\mathbf{A}}_L^2=G_{\oplus }{M}_{\oplus }^2+2J^2({\displaystyle \frac{1}{2}}{\stackrel{.}{r}}^2-{\displaystyle \frac{G_{\oplus }{M}_{\oplus }}{r}})=G_{\oplus }{M}_{\oplus }^2+2J^2ϝ,$$

(15)

where $ϝ$ is the conserved energy per unit mass.

Elliptic Integrals of the Second and of the Fourth Order with an Integration Over the True Anomaly Angle f

Finding the solutions for the propagation time T of the signal in terms of elliptic integrals of the second and the fourth order (as shall be explained, the order of the elliptic integral is related to the order of the polynomial in the nominator) is the third main goal of this paper. Initially it was performed in the paper [20] and some results have been summarized (with some new small additions and interpretations) in the review paper [22]. The unusual mathematical properties of the space-distributed transformations (9) - (11) (they will be clarified in the next sections) make it possible to define these elliptic integrals in terms of the rapidly changing dynamical parameter during the motion of the satellite - the true anomaly angle f. In fact, this should be considered as a first step in expressing the propagation time in terms of all the Kepler parameters $(f_{.},a_{.},e_{.},{\Omega }_{.},I_{.},{\omega }_{.})$, characterizing the transformations (9) - (11). This is a much more complicated problem and for the moment it is not clear whether the integrals over the other five variables $a_{.},e_{.},{\Omega }_{.},I_{.},{\omega }_{.}$ will be elliptic ones or not. Naturally, during the propagation of the signal from one satellite to another, the true anomaly angle will be the most rapidly changing and the other five parameters will change very slowly.

The motivation for searching of such solutions is important due to several reasons.

Motivation for Finding the Solutions for the Propagation Time from the Point of View of the Satellite Laser Communications

A. One of the essential ingredients of the mentioned concept about the autonomous navigation [8] is the laser optical communications between the satellites [10], in which the relative ranging and the relative velocity between the satellites depend on their base motion. If we take as an example the nanosatellites, the s.c. FSO (Free Space Optical Communications) [9,11] permit a narrow signal beam divergence and a satellite system with a smaller antenna in comparison with the RF (Radio Frequency) satellite system. So if two satellites move on two different space-oriented orbits and they communicate between each other, the propagation time for the (laser) signal, send by the signal-emitting satellite and the second signal-receiving satellite should take into account the motion of both satellites and the curving of the signal trajectory, so that the signal is received at the exact proper space-point. So in fact, there will be no ideal "narrow beam divergence", if the effects of General Relativity Theory are taken into account. Unfortunately, it was not the goal in the cited references [8,9,10,11] to account for these relativistic effects of signal propagation.

The problem here is that the application of General Relativity Theory in (primarily) satellite-ground stations communications began in the pioneering papers of Neil Ashby [34,37], but the investigations concerned perturbations of orbital elements (including the ones in the transformations (9) - (11), see [34]), frequency shifts, gravitational redshifts, first-order Dopler effect. In all these references there was no need to consider propagation time for signals between moving satellites in space with respect to one another.

Yet in 2002 Neil Ashby in the paper [35] formulated the following interesting problem (it concerns however not the propagation times, but the atomic times): what will be the synchronization error between the atomic clocks, situated on a rapidly moving GPS-transmitter at about $x=20$ 000 $km$ above the surface of the Earth and a receiver in low orbit around the Earth? If the receiver is moving with a velocity $v=7.6$$km/sec$, the estimated synchronization error between the atomic clocks was estimated to be $1.7$$\mu sec$ (micro-seconds, 1 $\mu sec={10}^{-6}$sec) according to the simple formulae ${\displaystyle \frac{xv}{c^2}}$. Now the problem can be formulated in another way: will there be a synchronization error between moving with respect to one another satellites on two different space-oriented orbits?

Motivation for Finding the Solutions for the Propagation Time for Satellites on Space-Oriented Orbits from the Point of View of the Method of Two Intersecting Null Cones

B. In order to account for the General Relativity effects on the signal trajectory between moving satellites, in the papers [18,19] and [22] a new theoretical method was proposed, called "the method of intersecting null cones". The idea of this mathematical method is that at the signal-emitting satellite with coordinates $(x_1,y_1,z_1)$ and at the signal-receiving satellite with coordinates $(x_2,y_2,z_2)$ two null cones are situated, which intersect each other. The first null cone is

$$d{s}_1^2=0=-(c^2+2V_1){\left(d{T}_1\right)}^2+(1-{\displaystyle \frac{2V_1}{c^2}})\left({\left(d{x}_1\right)}^2+{\left(d{y}_1\right)}^2+{\left(d{z}_1\right)}^2\right)$$

(16)

and the second null cone is

$$d{s}_2^2=0=-(c^2+2V_2){\left(d{T}_2\right)}^2+(1-{\displaystyle \frac{2V_2}{c^2}})\left({\left(d{x}_2\right)}^2+{\left(d{y}_2\right)}^2+{\left(d{z}_2\right)}^2\right).$$

(17)

The changing positions (Euclidean distance) between the satellites mean that the two four-dimensional null cones have to be additionally intersected with the six-dimensional hyperplane in terms of the variables $d{x}_1$, $d{y}_1$, $d{z}_1$, $d{x}_2$, $d{y}_2$, $d{z}_2$, which is obtained after taking the differential of the formulae for the Euclidean distance

$$R_{AB}^2={(x_1-x_2)}^2+{(y_1-y_2)}^2+{(z_1-z_2)}^2.$$

(18)

and thus the hyperplane equation in terms of the differentials is

$$d{R}_{AB}^2=2(x_1-x_2)d(x_1-x_2)+2(y_1-y_2)d(y_1-y_2)+2(z_1-z_2)d(z_1-z_2).$$

(19)

Thus the condition of signal exchange and finding the propagation times $T_1$ and $T_2$ mathematically is equivalent to finding $d{T}_1$ and $d{T}_2$ as solutions for the intersecting variety of the two null cones (16) and (17) with the hyperplane (19). Essentially, this is a problem, related to algebraic geometry. In [18,19] and [22] this has been performed for the case of a plane elliptical orbit, when there is no z variable in the corresponding equations. Now it remains to substitute in the equations (16), (17) and (19) the expressions for $T_1$ and $T_2$, which will be found for the case of space-oriented orbits. However, since the implementation of the mathematical formalism, related to elliptic integrals in the formalism of two null intersecting cones is very complicated and requires the use of algebraic geometry, it this paper we shall restrict mainly to elliptic functions and integrals in the framework only of the one null cone approach. It is evident, however, that finding the propagation time in terms of elliptic functions is the first step towards implementing the obtained formulaes in the approach of two null intersecting cones.

The necessity to perform such calculations is determined also from some new physical notions, defined in [18,19] and [22] for the case of plane elliptical orbits, given by the transformations (4). The solution of the equation in differentials for the null cone enabled to define the space-time distance, which was proved to have the property to be less than zero, equal to zero or greater than zero. If the space-time distance is assumed to be comparable to the Euclidean distance, then the s.c. ""compatibility condition for inter-satellite communications" is obtained. If substituted in the expression for the space-time distance, geodesic distance is obtained, which is only positive, and this has been proved to be valid in the cited papers for all the cases.

The null-cone equations (16) and (17) may be considered as a generalization of the null cones

$$\mid r-r_i\mid =c(t-t_i)i=1,2,3,4,$$

(20)

proposed in [35,36] for the case of a flat space-time and a responder with coordinates $(t,r)$ (t is the time of reception) and four transmitters (with synchronized atomic clocks) with coordinates $(t_i,r_i)$. From the nonlinear system (20) of four equations the coordinates $(t,r)$ of the reception event have to be found. In our case of only two null cones, the two propagation times $T_1$ and $T_2$ have to be found, and also the variable space-time distance from the additional hyperplane equation.

However, in order to be able to consider the above theoretical scheme to be of a general character, it has to be proved also for the more complicated case of space-oriented orbits.

Some Unusual and General Mathematical Properties of the Space-Oriented Transformations (9) - (11)

The transformations (9) - (11) have some unusual mathematical properties, which make it possible to express the propagation time in terms of elliptic integrals in terms of the true anomaly angle. The transformation (4) for the case of the plane elliptic orbit in matrix notations can be written as

$$\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{c}-ae\\ 0\end{array}\right)+\left(\begin{array}{cc}a& 0\\ 0& a\sqrt[.]{1-e^2}\end{array}\right)\left(\begin{array}{c}cosE\\ sinE\end{array}\right),$$

(21)

which represents a simple mapping $f:(a,e,E)\to (x,y)$, where the semi-major axis a and the eccentricity are kept constant. So this is a simple function $f\left(E\right)$.

The situation for the space-oriented transformations (9) - (11) is quite different, since if all the six Kepler parameters $(f_{.},a_{.},e_{.},{\Omega }_{.},I_{.},{\omega }_{.})$ are taken into consideration, it represents a nontrivial mapping

$$F:M(f_{.},a_{.},e_{.},{\Omega }_{.},I_{.},{\omega }_{.})\to N(x,y,z)$$

(22)

from a manifold with $n_1=6$ dimensions to a manifold with $n_2=3<n_1$ dimensions. Such a mapping of a manifold of more dimensions than the dimensions $n_2$ of the image manifold is called a submersion [39]. Some aspects of such a more general theory are considered in the recent paper [23]. Also, if one would like to find the propagation time of the signal, sent from one satellite on a space-oriented orbit with Kepler parameters $f_1,a_1,e_{1.},{\Omega }_{1.},I_{1.}{,}_{.}{\omega }_1$ to another satellite on another orbit with parameters $f_2,a_2,e_{2.},{\Omega }_{2.},I_{2.}{,}_{.}{\omega }_{2.}$, then two manifolds $M_1(f_1,a_1,e_{1.},{\Omega }_{1.},I_{1.}{,}_{.}{\omega }_{1.})$ and $M_2(f_2,a_2,e_{2.},{\Omega }_{2.},I_{2.}{,}_{.}{\omega }_{2.})$ should be defined, given by two sets of space-oriented celestial transformations (9) - (11), written for the two indices " 1" and "2". A basic problem is whether a mapping can be defined with respect to all the variables from the two sets of the Kepler parameters [23]

$${\overline{F}}_1:{\overline{M}}_1(f_1,a_1,e_{1.},{\Omega }_{1.},I_1,{\omega }_1,f_2,a_2,e_{2.},{\Omega }_{2.},I_{2.},{\omega }_{2.})\to {\overline{N}}_1(x_1,y_1,z_1,x_2,y_2,z_2)$$

(23)

and where $(x_1,y_1,z_1)$ and $(x_2,y_2,z_2)$ are the corresponding Cartesian coordinates of the two satellites. Such a mapping again will be a submersion of a 12 dimensional manifold ${\overline{M}}_1$ into the manifold ${\overline{N}}_1$, defined on the two sets of Cartesian coordinates for the two satellites $(x_1,y_1,z_1,x_2,y_2,z_2)$ . Such a possibility exists due to the Whitney theorem [39], according to which each connected smooth and closed n-dimensional manifold can be smoothly embedded in an $R^{2n+1}-$dimensional Euclidean space or it can be immersed in an $R^{2n}-$dimensional Euclidean space. In the given case, the $6-$dimensional manifold $M_1(f_1,a_1,e_{1.},{\Omega }_{1.},I_1,{\omega }_1)$ and the $6-$ dimensional manifold

$N_1(f_2,a_2,e_{2.},{\Omega }_{2.},I_{2.},{\omega }_{2.})$ can be either embedded in an $R^{2\times 6+1}=R^{13}$ Euclidean manifold of 13 dimensions (both sets of the Kepler parameters plus one more parameter s) or immersed in an $R^{2\times 6}=R^{12}-$dimensional Euclidean space. In particular, the validity of the Whitney theorem for this general case, which will not be considered in this paper, justifies the implementation of the approach of two null intersecting cones, based on the equations (16), (17) and (19). So this approach will be correct even when the initial propagation time $T_1$ is calculated only in terms of the first set of coordinates $(f_1,a_1,e_{1.},{\Omega }_{1.},I_1,{\omega }_1)$ or $(x_1,y_1,z_1)$ and the final propagation time $T_2$ is calculated only in terms of the second set of coordinates $(f_2,a_2,e_{2.},{\Omega }_{2.},I_{2.},{\omega }_{2.})$ or $(x_2,y_2,z_2)$.

General Mathematical Properties of the Space-Oriented Transformations (9) - (11) with Only the True Anomaly Angle f as the Dynamical Variable

Unlike the simplified transformation (21), the space-oriented transformations (9)-(11) [41] [40]

$$\left(\begin{array}{c}x\\ y\\ z\end{array}\right)=R_z(-\Omega )R_x(-I)R_z(-\omega )\left(\begin{array}{c}rcosf\\ rsinf\\ 0\end{array}\right),$$

(24)

is more complicated. The meaning of the above formulae is that expressions (9) - (11) can be obtained after performing three successive rotations $R_z(-\omega )$, $R_x(-I)$ and $R_z(-\Omega )$ with respect to the orbital vector $(rcosf$, $rcosf$, ${0)}^T$ (the transponed vector to the vector-column in (24)), where $R_z(-\omega )$ is the matrix of rotation at an angle $(-\omega )$ in the counterclockwise direction around the z axis, $R_x(-I)$ is the matrix of rotation at an angle $(-I)$ around the $x-$axis (again counterclockwise), $R_z(-\Omega )$ is the matrix of rotation at an angle $(-\Omega )$ around the $z-$axis [42,43]. Again, it is reminded that the $x-$axis is defined to have the direction of the vernal equinox and along the line of intersection of the equatorial plane of the Earth and the orbital plane of the Earth (i.e. the ecliptics).

In [41] the corresponding matrices $R_z(-\Omega )$, $R_x(-I)$ and $R_z(-\omega )$ have been denoted as $P_3,P_2$ and $P_1$. They have the following explicit forms

$$P_1\left(\omega \right):\equiv \left(\begin{array}{ccc}cos\omega & -sin\omega & 0\\ sin\omega & cos\omega & 0\\ 0& 0& 1\end{array}\right),P_2\left(I\right):\equiv \left(\begin{array}{ccc}1& 0& 0\\ 0& cosI& -sinI\\ 0& sinI& cosI\end{array}\right),$$

(25)

$$P_3(\Omega ):\equiv \left(\begin{array}{ccc}cos\Omega & -sin\Omega & 0\\ sin\Omega & cos\Omega & 0\\ 0& 0& 1\end{array}\right).$$

(26)

Since the above matrices do not depend on the true anomaly angle f, from the equation in matrix notations (24)

$$\left(\begin{array}{c}x\\ y\\ z\end{array}\right)=P_3(\Omega )P_2\left(I\right)P_1\left(\omega \right)\left(\begin{array}{c}rcosf\\ rsinf\\ 0\end{array}\right)$$

(27)

the expression $\left({\left(dx\right)}^2+{\left(dy\right)}^2+{\left(dz\right)}^2\right)$ in the null cone equation (6) can be found for the case when the orbital coordinates $(x,y,z)$ contain only the true anomaly f as the dynamical variable. For that purpose, the derivative of the vector-column

$${\displaystyle \frac{\partial }{\partial f}}\left(\begin{array}{c}x\\ y\\ z\end{array}\right)=P_3(\Omega )P_2\left(I\right)P_1\left(\omega \right){\displaystyle \frac{\partial }{\partial f}}\left(\begin{array}{c}rcosf\\ rsinf\\ 0\end{array}\right)$$

(28)

can be calculated, but it will not be necessary. The interesting fact is that the right-hand side contains the true anomaly fonly in the vector-column and the angles $\Omega ,I$ and $\omega$ in the $3\times 3$ matrix $A:=P_3(\Omega )P_2\left(I\right)P_1\left(\omega \right)$, but the square of the radius-vector $x^2+y^2+z^2$ will not contain any of the angles $\Omega ,I,\omega$. In fact, it can easily be calculated from the expressions (9) - (11) for $x,y,z$ that

$$x^2+y^2+z^2=r^2,\mathrm{where}r={\displaystyle \frac{a(1-e^2)}{1+ecosf}}$$

(29)

and also

$$dr={\displaystyle \frac{a(1-e^2)esinf}{{\left(1+ecosf\right)}^2}}.$$

(30)

Thus, by its construction r is a radius- vector in the orbital plane, depending only on fas a dynamical variable. This made possible the calculation of the elliptic integral for the propagation time T in terms only of the variable f.

The Celestial Time and its Relation to the Eccentric Anomaly Angle E and to the True Anomaly Angle f

There is one more physical reason for the propagation time to be expressed only through the dynamical variable f. As shall be shown further, for the plane elliptical orbit finding the propagation time T as a function of the eccentric anomaly E is physically justifiable because the eccentric anomaly E is related to the celestial time $t_{cel}$ through the Kepler equation (5) $E-esinE=n(t_{cel}-t_p)=M$. Similarly for space-oriented orbits, such a celestial time can be introduced by means of the expression

$$t_{cel}={\displaystyle \frac{1}{n}}\int {\displaystyle \frac{{\left(1-e^2\right)}^{\frac{3}{2}}}{{\left(1+ecosf\right)}^2}}df.$$

(31)

Thus, in both cases there is a correspondence $E\iff t_{cel}$ and $f\Rightarrow t_{cel}$. For the first case the correspondence is "two way", because for given $t_{cel}$ or M, the eccentric anomaly E can also be found as a solution of the transcendental Kepler equation. For the second case, the correspondence is "one-way", because the integral (31) can be solved analytically, although in some monographs on celestial mechanics it was claimed that only a numerical solution exists. The analytical solution of the integral (31) is

$$t_{cel}={\displaystyle \frac{\sqrt[.]{1-e^2}}{n}}[-e{\displaystyle \frac{sinf}{(1+ecosf)}}+$$

$$+{\displaystyle \frac{2}{sin\delta }}arctan\left(cotan{\displaystyle \frac{\delta }{2}}tan{\displaystyle \frac{f}{2}}\right)],$$

(32)

where $\delta$ is the following numerical parameter, depending on the eccentricity of the orbit

$$\delta =arccose.$$

(33)

However, an approximate correspondence $t_{cel}⟹f$ in this case cannot be established since f cannot be expressed from (32) if $t_{cel}$ is known. In this paper the solutions for the propagation time will be expressed through complicated elliptic integrals, depending on the eccentric anomaly E or on the true anomaly f. However, the above formulaes show that there is such a possibility to express the propagation time as complicated functions of the celestial time $t_{cel}$ by means of the Kepler equation (5) or the formulaes (32) and (33). This shall not be performed in this paper. If v is the celestial velocity and $d{t}_{cel}$ -the infinitesimal celestial time, then $vd{t}_{cel}$ can be replaced by

$$vd{t}_{cel.}=d{t}_{cel}\sqrt{{\left({\displaystyle \frac{dx}{d{t}_{cel}}}\right)}^2+{\left({\displaystyle \frac{dy}{d{t}_{cel}}}\right)}^2+{\left({\displaystyle \frac{dz}{d{t}_{cel}}}\right)}^2}$$

$$=\sqrt{{\left({\displaystyle \frac{dx}{df}}\right)}^2+{\left({\displaystyle \frac{dy}{df}}\right)}^2+{\left({\displaystyle \frac{dz}{df}}\right)}^2}df=v_{f}df,$$

(34)

where

$$v_f=\sqrt{1+e^2+2ecosf}$$

(35)

is the velocity, associated with the true anomaly angle f. Note also that the above expression will be correct if it (32) has no critical points, i.e. ${\displaystyle \frac{\partial v_f}{\partial f}}\ne o$. This shall not be investigated in this paper.

2.2.5. Fourth Aim of the Paper-Finding New Analytical Methods for Expressing Elliptic Integrals of the Second and of the Fourth Order

Essentially, in view of the application of elliptic integrals in GPS inter-satellite communications and also in several other areas of physics (cosmology, nonlinear equations and etc.), in the paper [21] general types of elliptic integrals were considered

$$J_n^{\left(4\right)}\left(y\right)=\int {\displaystyle \frac{y^{n}dy}{\sqrt{a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4}}},$$

$$J_n^{\left(3\right)}\left(x\right)=\int {\displaystyle \frac{x^{n}dx}{\sqrt{a{x}^3+b{x}^2+cx+d}}}$$

(36)

and some more general types of transformations

$$x={\displaystyle \frac{a}{y^2}}+b$$

(37)

were applied so that a and b represent complicated functions of the modulus of the elliptic integrals. The transformation (37) has been applied only with respect to a special type of elliptic integrals - the s.c. elliptic integrals of zero order (i.e. $n=0$) and in the Legendre form - this is the left integral in (38). The advantage of such an approach is that there is no restriction for the modulus of the elliptic integral to be small. Also, the theory of elliptic functions uses the standard parametrization of an elliptic curve of the type $y^2=4x^3-g_{2}x-g_3$ in terms of the Weierstrass function and its derivative, meaning that the function $x=\rho \left(z\right)$ and $y={\rho }^{\prime }\left(z\right)$ satisfy the above equation.

In this paper it will be proved that the propagation time T will be expressed by means of elliptic integrals of the types (with some coefficients in front)

$$\int {\displaystyle \frac{dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}},\int {\displaystyle \frac{y^{2}dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}.$$

(38)

Since it is known that integrals in the Weierstrass form $\int {\displaystyle \frac{dx}{\sqrt{4x^3-g_{2}x-g_3}}}$ can be solved analytically, it is evident that the transformation (37) can be used to find analytical expressions for the first integral in (38), which is called an elliptic integral of the first kind and of zero order in the Legendre form.

The second integral in (38) can also be solved analytically due to the recurrent chain of elliptic integrals of different orders, which usually is written as [44]

$$y^{n}Z=a_0(n+2)J_{n+3}+2a_1(2n+3)J_{n+2}$$

$$+6a_2(n+1)+2a_3(2n+1)J_n+n{a}_4{J}_{n-1},$$

(39)

where Z denotes the square of the fourth degree polynomial in the integral $J_n^{\left(4\right)}\left(y\right)$ in (36)

$$Z:\equiv \sqrt{a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4}.$$

(40)

Another important approach in the paper [21] is that one and the same integral can be represented by means of the transformation (37) in two different forms $J_n^{\left(4\right)}\left(y\right)$ and $J_n^{\left(3\right)}\left(x\right)$. Comparing the coefficient in the under-integral expressions in the denominators, the analytical expressions for the Weierstrass invariants $g_2$ and $g_3$ can be found.

3. Propagation Time for the Case of a Signal-Emitting Satellite, Moving along a Plane Elliptical Orbit

3.1. Propagation Time Without Any Approximations

As mentioned in the introduction of the paper, the null-cone metric is taken in the standard form (6)

$$d{s}^2=-c^2\left(1+{\displaystyle \frac{2V}{c^2}}\right){\left(dT\right)}^2+\left(1-{\displaystyle \frac{2V_{.}}{c^2}}\right)\left({\left(dx\right)}^2+{\left(dy\right)}^2+{\left(dz\right)}^2\right)=0,$$

where $V={\displaystyle \frac{G_{\oplus }{M}_{\oplus }}{r}}$ is the standard gravitational potential of the Earth. No account is taken of any harmonics due to the spherical form of the Earth since the GPS orbits are situated at a distance more than 20000 $km$ (considered from the centre of the Earth). The Kepler parametrization of the space coordinates $x=x(a,e,E)$ and $y=y(a,e,E)$ for the plane elliptical satellite orbit is (4) $x=a(cosE-e)$ , $y=a\sqrt[.]{1-e^2}sinE$, where a is the semi-major axis of the orbit, e is the eccentricity and E is the eccentric anomaly angle, related to the motion of the satellite along the plane elliptical orbit.

If the right-hand side of Eq. (6) is divided and multiplied by ${\left(dt\right)}^2$, where $t=t_{cel}$ is the celestial time from the Kepler equation $E-esinE=n(t_{cel}-t_p)$ ($t_p$ is the time of perigee passage), one can obtain

$$(c^2+2V){\left(dT\right)}^2=(1-{\displaystyle \frac{2V}{c^2}}){\mathbf{v}}^2,$$

(41)

where ${\mathbf{v}}^2$ is the square of the satellite velocity along the orbit

$${\mathbf{v}}^2=v_x^2+v_y^2={\left({\displaystyle \frac{dx}{dt}}\right)}^2+{\left({\displaystyle \frac{dy}{dt}}\right)}^2.$$

(42)

Taking into account the Kepler parametrization Eq. (4), the above expression for the velocity for the case of plane motion can be rewritten as

$$\mathbf{v}=\sqrt{v_x^2+v_y^2}={\displaystyle \frac{na}{\left(1-ecosE\right)}}\sqrt[.]{1-e^2{cos}^{2}E}$$

(43)

and $n=\sqrt[.]{{\displaystyle \frac{G_{\oplus }{M}_{\oplus }}{a^3}}}$ is the mean motion. Expressing $dT$ from (41), finding the celestial time $t_{cel}$ in terms of E from the Kepler equation and performing the integration over the eccentric anomaly angle E, one can obtain the expression for the propagation time (for the case of no restrictions imposed)

$$T=\int {\displaystyle \frac{\mathbf{v}}{c}}\sqrt[.]{{\displaystyle \frac{(c^2-2V)}{(c^2+2V)}}}dE+C$$

$$={\displaystyle \frac{a}{c}}\int \sqrt[.]{{\displaystyle \frac{(1-e^2{cos}^{2}E)\left[a\left(1-ecosE\right)-\beta \right]}{\left[a\left(1-ecosE\right)+\beta \right]}}}dE+C$$

(44)

The constant $\beta$ is determined as

$$\beta ={\displaystyle \frac{2G{M}_{\oplus }}{a{c}^2}}.$$

(45)

The other constant C can be determined as a result of the integration of the null cone equation from some initial condition - for example, the propagation time T should be equal to zero, when $E=0$ (if $E_{init}=0$). Further, making the substitution

$$1-ecosE=y\Rightarrow dE={\displaystyle \frac{dy}{esinE}}$$

(46)

in the integral (44) and denoting $\overline{\beta }={\displaystyle \frac{\beta }{a}}$, one can represent (44) as

$$T={\displaystyle \frac{a}{c}}\int \sqrt[.]{{\displaystyle \frac{y(2-y)(y-\overline{\beta })}{(y+\overline{\beta })\left[e^2-{\left(1-y\right)}^2\right]}}}dy.$$

(47)

This integral is of the form

$$T={\displaystyle \frac{a}{c}}\int \sqrt[.]{{\displaystyle \frac{a_1{y}^3+a_2{y}^2+a_{3}y}{b_1{y}^3+b_2{y}^2+b_{3}y+b_4}}}dy,$$

(48)

where the numerical constants $a_1$, $a_2$, $a_3$, $b_1$, $b_2$, $b_3$, $b_4$ can easily be calculated. Their analytical expressions are not given here, because the above expression will not be used further.

3.1.1. Is the Integral (48) an Abelian One?

If the square root, representing an irrational function of the variable y, is denoted as

$$x=\sqrt[.]{{\displaystyle \frac{a_1{y}^3+a_2{y}^2+a_{3}y}{b_1{y}^3+b_2{y}^2+b_{3}y+b_4}}}=\sqrt[.]{P\left(y\right)},$$

(49)

then the integral

$$\int xdy=\int \sqrt[.]{P\left(y\right)}dy$$

(50)

is a partial case of a class of integrals (often denoted in mathematical literature as $\int R(x,y)dy$). In case $P\left(y\right)$ is an algebraic polynomial, integrals $\int \sqrt[.]{P\left(y\right)}dy$ are known in mathematics as abelian integrals (see the monograph by Prasolov and Solovyev [44]), related to the curve

$$F(x,y):=x^2-P\left(y\right)=0.$$

(51)

Since in the case of the expression (49) $P\left(y\right)$ is a complicated irrational function of an rational expression of y, the integral (48) is not an abelian one. Further the expression Eq. (44) will not be used, because a justifiable and reasonable physical approximation shall be applied, which will enable us to express the solution in terms of elliptic integrals. Nevertheless, in view of the constantly improving accuracy for measuring the propagation time, this might represent an interesting problem for future mathematical research. In the next sections, when defining the elliptic integrals for the more general case of arbitrary polynomials of third- and fourth- degree in the denominator, it will be explained in details why integrals of the form (48) cannot belong to the class of abelian integrals.

3.2. The Useful Approximation for a Small Gravitational Potential, Compared to the Square of the Velocity of Light

Further we shall be interested in the case

$$\beta ={\displaystyle \frac{2V}{c^2}}={\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^{2}a}}\ll 1,$$

(52)

which can be assumed to be fulfilled. For the parameters of the $GPS$ orbit - $a=26561$$\left[km\right]$ , in the review paper [31] the constant $\beta$ can exactly be calculated to be $0.334\times {10}^{-9}$, with the velocity of light taken to be $c=299792458$$\left[{\displaystyle \frac{m}{sec}}\right].$ The geocentric gravitational constant $G_{\oplus }{M}_{\oplus }$ (obtained from the analysis of laser distance measurements of artificial Earth satellites) can be taken to be equal to $G_{\oplus }{M}_{\oplus }=\left(3986004.405\pm 1\right)\times {10}^8\left[{\displaystyle \frac{m^3}{{sec}^2}}\right]$, but the value of $G_{\oplus }{M}_{\oplus }$ can vary also in another range from $G_{\oplus }{M}_{\oplus }=3986056.75236\times {10}^8\left[{\displaystyle \frac{m^3}{{sec}^2}}\right]$ to the value $G_{\oplus }{M}_{\oplus }=3987999.07898\times {10}^8\left[{\displaystyle \frac{m^3}{{sec}^2}}\right]$ due to the uncertainties in measuring the Newton gravitational constant $G_{\oplus }$. The mass of the Earth can be taken approximately to be $M_E\approx 5.97\times {10}^{24}$ $\left[kg\right]$. One of the latest values for $G_{\oplus }$ from deep space experiments was reported in the paper [62] to be $(6.674+0.0003).{10}^{-11}$ $\left[{\displaystyle \frac{m^3}{kg.{sec}^2}}\right]$. The latest value for G from Encyclopedia Britanica (updated 23 September 2024) is $6.67\pm 0.00015\times {10}^{-11}$$\left[{\displaystyle \frac{m^3}{kg.s^2}}\right]$.

3.3. Derivation of the Propagation Time Under the Approximation $\beta ={\displaystyle \frac{2V}{c^2}}\ll 1$ and Physical Justification of the Obtained Result

After decomposing the under-integral expression with the square root in Eq. (48) and leaving only the first-order term in ${\displaystyle \frac{2V}{c^2}}$, one can obtain

$$T=\int {\displaystyle \frac{\mathbf{v}}{c}}\sqrt[.]{{\displaystyle \frac{(1-\frac{2V}{c^2})}{(1+\frac{2V}{c^2})}}}dt\approx \int {\displaystyle \frac{\mathbf{v}}{c}}(1-{\displaystyle \frac{2V}{c^2}})dt=I_1+I_2=$$

(53)

$$={\displaystyle \frac{a}{c}}\int \sqrt[.]{1-e^2{cos}^{2}E}dE-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}\int \sqrt[.]{{\displaystyle \frac{1+ecosE}{1-ecosE}}}dE.$$

(54)

This expression is consistent from a physical point of view due to the following reasons:

1. The coefficient ${\displaystyle \frac{a}{c}}$ as a ratio of the large semi-major axis of the orbit and the velocity of light $c=299792458$$\left[{\displaystyle \frac{m}{sec}}\right]$ will have a dimension $[m/{\displaystyle \frac{m}{sec}}]=[sec]$, as it should be. The second coefficient ${\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}$ has a corresponding dimension $[{\displaystyle \frac{m^3}{{sec}^2}}:{\displaystyle \frac{m^3}{{sec}^3}}]=[sec]$, which clearly proves that formulae Eq. (54) has the proper dimensions.

2. The formulae in fact gives the propagation time T of the signal, emitted by the satellite at some initial position (given by the eccentric anomaly angle $E_{init}$), and the final point (given by $E_{fin}$) on the trajectory of the signal. However, there is no real transmission of signals, because the final point $E_{fin}$ should be given. Let us clarify this moment-if there is a transmission of signals (emission+reception) the point of reception should be determined, depending on the equations for the two null cones in the framework of the two intersecting null cones, which was briefly outlined in the Introduction. In this sense, the integrals (53) and (54) represent curvilinear integrals (depending on the determined in advance initial and final points of integration), derived after the application of the null-cone equation (in a general form written as $F=g_{\alpha \beta }d{x}^{\alpha }d{x}^{\beta }=0$). This propagation time in fact can be found from the expression (3) after expressing the differentials $dT$ and $d{x}_i$ by means of the differential $dE$ of the eccentric anomaly angle E. Thus from (3) it is obtained for the propagation time T after integration

$$T=\pm {\displaystyle \frac{1}{c}}\int {\displaystyle \frac{1}{\sqrt[.]{-g_{00}}}}\sqrt[.]{\left(g_{ij}+{\displaystyle \frac{g_{0i}{g}_{0j}}{-g_{00}}}\right){\displaystyle \frac{d{x}^i}{dE}}{\displaystyle \frac{d{x}^j}{dE}}}dE+$$

$$+{\displaystyle \frac{1}{c}}\int \left({\displaystyle \frac{g_{0j}}{-g_{00}}}\right){\displaystyle \frac{d{x}^j}{dE}}dE.$$

(55)

In the same manner, the space variables $x^i$ can be parametrized by some other variable. For example, in the next section by using the celestial transformations (9) - (11) for the case of space-oriented orbits, the true anomaly angle f will be used and the formulae for T will be the same as in (55) but with the eccentric anomaly E replaced by the true anomaly f, i.e. $E\to f$. This case will be investigated in the next sections. In order to obtain the propagation time T in the form of the expression (54), the additional assumption $\beta ={\displaystyle \frac{2V}{c^2}}\ll 1$ has to be taken into account, when transforming expression (55).

3. The expression for the propagation time T (54) establishes the relation between the eccentric anomaly angle E and the propagation time $T-$$E⟹T$. Since T will be expressed through elliptic integrals and their analytical solutions will be complicated, it might turn out to be impossible to establish the correspondence in the other direction $T\Rightarrow E$, as was in the case of the "dual" correspondence $t_{cel}\iff E$. The meaning of the correspondence $E⟹T$ is the following: if some initial and final eccentric anomaly angles $E_{init}$ and $E_{fin}$ for the position of the satellite are given, to them will correspond a propagation time of the signal, found as the definite integral

$$\Delta T=T_{fin}-T_{init}={\displaystyle \frac{a}{c}}{\int }_{E_{init}}^{E_{fin}}\sqrt[.]{1-e^2{cos}^{2}E}dE$$

$$-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}{\int }_{E_{init}}^{E_{fin}}\sqrt[.]{{\displaystyle \frac{1+ecosE}{1-ecosE}}}dE.$$

(56)

This expression establishes the correspondence ($E_{init},E_{fin}$) $⟹T$, due to the curved signal trajectory and the fact that there is no real transmission of the signal, this does not mean that the final point of propagation of the signal should lie necessarily on the orbit.

4. It should be stressed that T was found as a solution of the null cone equation (6) and because of that, it is the time for propagation of the signal and definitely not related to the real motion of the satellite. This statement is also confirmed by the fact the null cone equation $F=0$ is proved to be compatible with the geodesic equation [29] for light-like geodesics

$${\displaystyle \frac{d^2{x}_{\nu }}{d{\tau }^2}}+{\Gamma }_{\alpha \beta }^{\nu }{\displaystyle \frac{d{x}^{\alpha }}{d\tau }}{\displaystyle \frac{d{x}^{\beta }}{d\tau }}=0,$$

(57)

where $\tau$ is the proper time along the light ray and $\nu ,\alpha ,\beta =0,1,2,3$. In other words, the zero-length (light-like) geodesics (representing the trajectory of the signal) are determined also by the null cone equation $g_{\alpha \beta }d{x}^{\alpha }d{x}^{\beta }=0$, because the null cone equation is a first integral of the geodesic equation Eq. (57). This point will be clarified in more details in a separate chapter, since it is confirmed also by a theorem in differential geometry [30,32]. This was explained also in the introduction of this paper and also in the recent paper [23].

5. Let now $x,y,z$ be the space coordinates, reached by the signal after it has been emitted. From the null cone equation (6), after using the approximation $\beta ={\displaystyle \frac{2V}{c^2}}\ll 1$, it can be obtained

$$\left({\left(dx\right)}^2+{\left(dy\right)}^2+{\left(dz\right)}^2\right)\approx c^2\left[1+{\displaystyle \frac{4V}{c^2}}\right]{\left(dT\right)}^2\ll 3c^2{\left(dT\right)}^2.$$

(58)

Since $cdT$ is the infinitesimal distance, travelled by the light signal, the above inequality means that the distance travelled by light is much greater than the Euclidean distance. In other words, the light trajectory signal is a curved one, in accord with the meaning of the Shapiro delay formulae (1) about the delay of the signal under the action of the gravitational field.

6. The trigonometric function ${cos}^{2}E$ has the same values for $E=\alpha$ and $E=180-\alpha$, so the correspondence eccentric anomaly angle E and the propagation time T is not reliable for angles in the second quadrant. The same refers for values of E in the first and the fourth quadrant since $cosE=cos(360-E)$.

7. The propagation time T Eq. (54) should be a real-valued expression, since time cannot be complex. This fact is not evident from the beginning of the calculations, but will be proved in the next subsections. Remarkably, this important property will turn out to be valid also for the next case of a signal, emitted and percepted by satellites on a space-oriented orbits.

8. The expressions for T (53) and (54) under the approximation $\beta ={\displaystyle \frac{2V}{c^2}}\ll 1$ with account of the vector r in the orbital plane is

$${\displaystyle \frac{2V}{c^2}}={\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^{2}a(1-ecosE)}}\ll 1,$$

(59)

which can be assumed to be fulfilled, because $\beta =$${\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^{2}a}}\ll 1$. For the parameters of a $GPS$ orbit, the constant $\beta$ can exactly be calculated to be $0.33.{10}^{-9}$, which justifies the above strong inequality. In the strict mathematical sense and for the given case, the above inequality (59) takes place when $cosE\ll {\displaystyle \frac{1}{e}}-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^{2}ae}}$. This is always fulfilled for the case of the small eccentricity of the GPS orbits ($e=0.01323881349526$), because the first term will be a large number of the order ${\displaystyle \frac{1}{13}}.{10}^3$, while the second term is a very small number of the order $1.4.{10}^{-8}$. Another values for the parameter $\beta$ have been obtained in the literature. For example, in the review article by J. Pascual Sanchez [31], it was obtained for the approximate radius of the satellite orbit $r_s=26561$$\left[km\right]$

$${\displaystyle \frac{G_{\oplus }{M}_{\oplus }}{r_s{c}^2}}=0.167.{10}^{-9}.$$

(60)

If multiplied by 2 (in order to obtain the constant $\beta =$${\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^{2}a}}$), this will give $0.334.{10}^{-9}$, which is in full accord (up to the second digit) with the estimate in this paper $\beta =0.33.{10}^{-9}$.

4. Mathematical Structure of the Expression for the Propagation Time (Signal, Emitted by a Satellite on a Plane Elliptical Orbit), Related to Zero-Order Elliptic Integrals of the First, Second and of the Third Rank

The first term in expression (54) after performing the simple transformation ${\displaystyle \frac{\pi }{2}}-E=\overline{E}$ can be brought to the following form

$$T_1=\underset{0}{\overset{E}{\int }}\sqrt[.]{1-e^2{cos}^{2}E}dE=\underset{0}{\overset{E}{\int }}\sqrt[.]{1-e^2{sin}^2({\displaystyle \frac{\pi }{2}}-E)}dE$$

$$=-\underset{{\displaystyle \frac{\pi }{2}}}{\overset{{\displaystyle \frac{\pi }{2}}-\overline{E}}{\int }}\sqrt[.]{1-e^2{sin}^2\overline{E}}d\overline{E}.$$

(61)

This represents an elliptic integral of the second kind. Note that the s.c.modulus of the elliptic integral is equal to the eccentricity e of the orbit.

The second term in Eq. (54) after performing the substitution

$$\sqrt[.]{{\displaystyle \frac{1+ecosE}{1-ecosE}}}=\overline{y}$$

(62)

and introducing the notations

$${\tilde{k}}^2={\displaystyle \frac{1-e}{1+e}}=q,{\displaystyle \frac{\overline{y}}{\tilde{k}}}=y,{\left({\displaystyle \frac{\overline{y}}{\tilde{k}}}\right)}^2={\tilde{y}}_1,{\tilde{y}}_1=-\tilde{\tilde{y}}$$

(63)

can be written as a sum of two integrals

$$T_2=-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}\int \sqrt[.]{{\displaystyle \frac{1+ecosE}{1-ecosE}}}dE=I_2^{A)}+I_2^{\left(B\right)}.$$

(64)

It should be noted that $\tilde{k}$ in eq. (63) is simply a notation, representing the real expression $\tilde{k}=\sqrt{{\displaystyle \frac{1-e}{1+e}}}=\sqrt{q}$, since $0<{\displaystyle \frac{1-e}{1+e}}<1$ (the eccentricity of the orbit e is always less than one, $0<e<1$). We shall call this second part (64) of the expression for the propagation time the $O\left(c^{-3}\right)$ correction.

4.1. Mathematical Proof of the Real-Valuedness of the Separate Expressions for the Propagation Time for the Case of Plane Elliptic Orbits

4.1.1. The First Real-Valued Elliptic Integral $I_2^{\left(A\right)}$ of the First Kind in the Weierstrass Form in the $O\left(c^{-3}\right)$ Correction of the Propagation Time

Using this notation, the first term $I_2^{\left(A\right)}$ in eq. (64) can be represented as

$$I_2^{\left(A\right)}={\displaystyle \frac{4GM}{c^3}}{\displaystyle \frac{1}{\tilde{k}\sqrt[.]{e^2-1}}}\int {\displaystyle \frac{d{\tilde{y}}_1}{\sqrt[.]{{\tilde{y}}_1\left({\tilde{y}}_1-\frac{1}{{\tilde{k}}^4}\right)\left({\tilde{y}}_1-1\right)}}}.$$

(65)

In terms of the variable

$${\tilde{y}}_1=\left({\displaystyle \frac{1+e}{1-e}}\right).\left({\displaystyle \frac{1+ecosE}{1-ecosE}}\right),$$

(66)

the integral (65) represents an elliptic integral of zero order and of the first kind, written in the Weierstrass form (i.e. a polynomial of third order under the square root in the denominator). Due to the presence of the $\sqrt[.]{e^2-1}$ term in the denominator (again, $0<e\le 1$), it might seem that expression (65) is imaginary. However, this is not true, because it can be rewritten as

$$I_2^{\left(A\right)}={\displaystyle \frac{4GM}{c^3}}{\displaystyle \frac{1}{\tilde{k}i\sqrt[.]{1-e^2.}\left(-i\right)}}\int {\displaystyle \frac{d\tilde{\tilde{y}}}{\sqrt[.]{\tilde{\tilde{y}}\left(\tilde{\tilde{y}}+1\right)\left(\tilde{\tilde{y}}+\frac{1}{{\tilde{k}}^4}\right)}}},$$

(67)

thus representing an elliptic integral of zero order and of the first kind. The classification of elliptic integrals, depending on the degree of the polynomial in the denominator is given in the monographs [44,63] and [64]. Since $\tilde{k}$ is real-valued expression, the two imaginary units i in the denominator appear from $\sqrt[.]{e^2-1}=i\sqrt[.]{1-e^2}$ and also from the square root, after performing the variable transformation from ${\tilde{y}}_1$ to $\tilde{\tilde{y}}$ in ${\tilde{y}}_1\left({\tilde{y}}_1-{\displaystyle \frac{1}{{\tilde{k}}^4}}\right)\left({\tilde{y}}_1-1\right)$ in eq. (65). As a result, a factor $\sqrt{{(-1)}^3}=\sqrt{i^6}=i^3=-i$ will appear in the denominator of eq. (67). The whole expression eq. (67) turns out to be a real-valued one and with a negative sign, because $d{\tilde{y}}_1=-d\tilde{\tilde{y}}$ . Since this expression is a part of the expression for the propagation time, this is physically reasonable.

4.1.2. The Second Real-Valued Elliptic Integral $I_2^{\left(B\right)}$ of the Third Kind

The second term $I_2^{\left(B\right)}$ in Eq. (64) can also be written in the form of a real-valued expression, as it should be

$$I_2^{\left(B\right)}={\displaystyle \frac{4GM}{c^3{q}^2}}{\displaystyle \frac{1}{\sqrt[.]{1-e^2}}}\int {\displaystyle \frac{d{\tilde{y}}_1}{({\tilde{y}}_1-\frac{1}{q})\sqrt[.]{{\tilde{y}}_1({\tilde{y}}_1+1)({\tilde{y}}_1+\frac{1}{q^2})}}}.$$

(68)

This is an elliptic integral of the third kind. Consequently, the whole expression for the propagation time can be represented as a sum of elliptic integrals of the second, of the first and of the third kind [44].

5. Propagation Time for a Signal, Emitted and Intercepted by Satellites on a Space-Oriented Orbit

5.1. Three-Dimensional Orbit Parametrization and the General Formulae for the Orbit Parametrization

The formulaes for the three-dimensional parametrization of space-oriented orbits previously were given by the expressions (9) - (11). Again, after using the null-cone equation (53) $T=\int {\displaystyle \frac{\mathbf{v}}{c}}\sqrt[.]{{\displaystyle \frac{(1-\frac{2V}{c^2})}{(1+\frac{2V}{c^2})}}}$$dt\approx \int {\displaystyle \frac{\mathbf{v}}{c}}(1-{\displaystyle \frac{2V}{c^2}})dt$, remembering relation (34) $vd{t}_{cel}=v_{f}df$ and also expression (35) $v_f={\displaystyle \frac{na}{\sqrt{1-e^2}}}\sqrt{1+e^2+2ecosf}$ for the velocity $v_f$, associated with the true anomaly angle f, the following expression can be obtained for the propagation time $\overline{T}$ for the three-dimensional case

$$\overline{T}=\int {\displaystyle \frac{v}{c}}(1-{\displaystyle \frac{2V}{c^2}})dt={\tilde{T}}_1+{\tilde{T}}_2={\displaystyle \frac{1}{c}}\int vdt-{\displaystyle \frac{2}{c^3}}\int vVdt.$$

(69)

We denote the propagation time for the three-dimensional case (9) - (11) by $\overline{T}$, in order to distinguish it from the propagation time T for the case of a satellite, moving along a plane elliptical orbit.

5.2. Analytical Calculation of the First Elliptic Integral (First $O\left({\displaystyle \frac{1}{c}}\right)$ Correction) Without the Use of Elliptic Integrals

We shall present a calculation, which shows that some integrals can be calculated both analytically, and also by the use of elliptic integrals. Although the theory of elliptic integrals is considered to be well-developed, there are problems, not well studied. One such problem, as it will be demonstrated further is that after performing some variable transformations, the elliptic integrals will be possible to be solved analytically, while after some other transformations, this will not be possible.

Let us take the first $O\left({\displaystyle \frac{1}{c}}\right)$ propagation time correction

$${\tilde{T}}_1={\displaystyle \frac{1}{c}}\int v_{f}df={\displaystyle \frac{na}{c\sqrt{1-e^2}}}\int \sqrt{1+e^2+2ecosf}df={\displaystyle \frac{na}{c\sqrt{1-e^2}}}{\tilde{T}}_1^{\left(1\right)}$$

(70)

and let us perform the series of six subsequent variable transformations

$$\sqrt{1+e^2+2ecosf}=y,\overline{y}={(e+1)}^2-y^2,$$

(71)

$$\sqrt{1-{\displaystyle \frac{\overline{y}{q}^2}{\overline{y}+B_1}}}=z,z^2=z_1,1-z_1=m,\overline{m}=m-{\displaystyle \frac{q^2}{2}}.$$

(72)

The first transformation $\sqrt{1+e^2+2ecosf}=y$ in (71) will be used only in this paragraph, further another transformation $\tilde{y}(1+e)=\sqrt{1+e^2+2ecosf}$ will be applied.

After applying the sequence of transformations (71) and (72), the following sum of two integrals for ${\tilde{T}}_1$ will be obtained

$${\tilde{T}}_1=-{\displaystyle \frac{nai}{c\sqrt{(1-e^2)q}}}\left[\int {\displaystyle \frac{d\overline{m}}{\left(\overline{m}-\frac{q^2}{2}\right)\sqrt{\overline{m}-\frac{q^2}{2}}}}+\int {\displaystyle \frac{d\overline{m}}{\left(\overline{m}+\frac{q^2}{2}\right)\sqrt{\overline{m}-\frac{q^2}{2}}}}\right],$$

(73)

where again the notation $q={\displaystyle \frac{1-e}{1+e}}$ has been used and each of the integrals inside the square bracket will be denoted correspondingly as ${\tilde{T}}_1^{\left(1\right)}$ and ${\tilde{T}}_1^{\left(2\right)}$. Making use of the analytically calculated integral from the book [66]

$$\int {\displaystyle \frac{dx}{(x\pm p)\sqrt{a+2bx+c{x}^2}}}$$

$$=-{\displaystyle \frac{1}{\sqrt{a\mp 2bp+c{p}^2}}}ln\left\{{\displaystyle \frac{\sqrt{a+2bx+c{x}^2}+\sqrt{a\mp 2bp+c{p}^2}}{x\pm p}}+{\displaystyle \frac{b\mp cp}{\sqrt{a\mp 2bp+c{p}^2}}}\right\},$$

(74)

one can derive for the second integral ${\tilde{T}}_1^{\left(2\right)}$ in Eq. (73) the following expression

$${\tilde{T}}_1^{\left(2\right)}=-{\displaystyle \frac{na(1+e)\sqrt{2}}{c\sqrt{1-e^2}\sqrt{q}\sqrt{3e^2+2e+3}}}ln\left[\left({\displaystyle \frac{\sqrt{2}{m}_1(f_b;r_b)+\frac{1}{2}{m}_2(f_b;r_b)}{\sqrt{2}{m}_1(f_a;r_a)+\frac{1}{2}{m}_2(f_a;r_a)}}\right)\left({\displaystyle \frac{m_3(f_a;r_a)}{m3(f_b;r_b)}}\right)\right]$$

(75)

and $m_1(f;r)$, $m_2(f;r)$, $m_3(f;r)$ are expressions, written in terms either of the initial and final true anomaly angles $f_a$ and $f_b$ or, of the initial distance $r_a$ (at which the emission of the signal takes place) and the final point $r_b$ of reception, corresponding to the propagation time ${\tilde{T}}_1^{\left(2\right)}$. The first integral ${\tilde{T}}_1^{\left(1\right)}$ can be calculated analogously. It should be noted that both ${\tilde{T}}_1^{\left(1\right)}$ and ${\tilde{T}}_1^{\left(2\right)}$ are real-valued expressions and moreover, the logarithmic term is again present, as in the original Shapiro delay formulae (1).

Further this analytical method will not be used because a transformation will be proposed, enabling to find analytical expressions for zero-order elliptic integrals in the Weierstrass or in the Legendre form.

5.3. Analytical Calculation of the First $O\left({\displaystyle \frac{1}{c}}\right)$ Correction by Means of Elliptic Integrals

Let us calculate the integral (70) ${\tilde{T}}_1={\displaystyle \frac{na}{c\sqrt{1-e^2}}}\int \sqrt{1+e^2+2ecosf}df$ by making use of the substitution

$$y=\sqrt{{\displaystyle \frac{1}{q}}{\displaystyle \frac{(1+ecosE)}{(1-ecosE)}}},q={\displaystyle \frac{1-e}{1+e}}$$

(76)

and also of the well-known relation from celestial mechanics between the eccentric anomaly angle E and the true anomaly angle f [25]

$$tan{\displaystyle \frac{f}{2}}=\sqrt{{\displaystyle \frac{1-cosf}{1+cosf}}}=\sqrt{{\displaystyle \frac{1+e}{1-e}}}tan{\displaystyle \frac{E}{2}}.$$

(77)

Then we can write the integral ${\tilde{T}}_1$ in the form of an elliptic integral of the second order and of the first kind in the Legendre form

$${\tilde{T}}_1=-2i{\displaystyle \frac{na}{c}}{q}^{{\displaystyle \frac{3}{2}}}\int {\displaystyle \frac{y^{2}dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}=-2i{\displaystyle \frac{na}{c}}{q}^{{\displaystyle \frac{3}{2}}}{\tilde{J}}_2^{\left(4\right)}(y;q).$$

(78)

Thus the integral (70) is very interesting from a mathematical point of view, because it becomes evident that it is an elliptic one after performing the transformations (76) and (77). The integral ${\tilde{J}}_2^{\left(4\right)}(y;q)$ in (78) can also be represented as

$${\tilde{J}}_2^{\left(4\right)}(y;q)=\int {\displaystyle \frac{y^{2}dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}$$

(79)

$$=-{\displaystyle \frac{1}{q^2}}\int {\displaystyle \frac{(1-q^2{y}^2)dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}+{\displaystyle \frac{1}{q^2}}\int {\displaystyle \frac{dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}$$

(80)

$$=-{\displaystyle \frac{1}{q^2}}\int \sqrt{1-q^2{sin}^2\phi }d\phi +{\displaystyle \frac{1}{q^2}}\int {\displaystyle \frac{dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}.$$

(81)

The first integral in Eq. (81) is an elliptic integral of the second kind (denoted usually by $E\left(\phi \right)=\int \sqrt{1-q^2{sin}^2\phi }d\phi$, where $x=sin\phi$). So we obtain a relation between the second-order elliptic integral ${\tilde{J}}_2^{\left(4\right)}(y;q)$ of the first kind in the Legendre form, the zero-order elliptic integral of the first kind in the Legendre form ${\tilde{J}}_0^{\left(4\right)}(y;q)$ (the second integral in (81)) and the elliptic integral $E\left(\phi \right)$ (see the monograph [44])

$${\tilde{J}}_2^{\left(4\right)}(y;q)=\int {\displaystyle \frac{y^{2}dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}=-{\displaystyle \frac{1}{q^2}}E\left(\phi \right)+{\displaystyle \frac{1}{q^2}}{\tilde{J}}_0^{\left(4\right)}(y;q).$$

(82)

Since in the preceding section an analytical expression was obtained for ${\tilde{T}}_1$ without the use of any elliptic functions, the obtained result in Eq. (81) means that for the investigated case, elliptic integrals of second order can be expressed through elementary functions.

Using the relation between the eccentric anomaly angle E and the true anomaly angle f (derived from Eq. (77) ), it is interesting to express the first $\left(O\left({\displaystyle \frac{1}{c}}\right)\right)$ propagation time correction Eq. (78) also in terms of the variable E

$${\tilde{T}}_1={\displaystyle \frac{na}{c}}\sqrt{1-e^2}\int {\displaystyle \frac{1}{(1-ecosE)}}\sqrt{{\displaystyle \frac{1+ecosE}{1-ecosE}}}dE.$$

(83)

This integral resembles the second integral $-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}\int \sqrt[.]{{\displaystyle \frac{1+ecosE}{1-ecosE}}}dE$ in the $O\left({\displaystyle \frac{1}{c^3}}\right)$ time correction Eq. (54) for the case of plane elliptical orbit, but in the case the integral is with another coefficient and is modified with the term ${\displaystyle \frac{1}{(1-ecosE)}}$, multiplying the square root.

5.4. Second Analytical Calculation of the First Time $O\left({\displaystyle \frac{1}{c}}\right)$ Correction in Terms of Second-Order Elliptic Integrals and Proof of the Real-Valuedness of the $O\left({\displaystyle \frac{1}{c}}\right)$ Correction

This second calculation will not make any use of the eccentric anomaly angle variable E. The main purpose of this second analytical calculation will be prove the real-valuedness of the expression (78) ${\tilde{T}}_1=-2i{\displaystyle \frac{na}{c}}{q}^{{\displaystyle \frac{3}{2}}}\int {\displaystyle \frac{y^{2}dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}=-2i{\displaystyle \frac{na}{c}}{q}^{{\displaystyle \frac{3}{2}}}{J}_2^{\left(4\right)}(y,q)$, which means that the second-order elliptic integral in the Legendre form should be imaginary-valued.

Let us apply the variable transformation

$$\tilde{y}={\displaystyle \frac{\sqrt{1+2ecosf+e^2}}{1+e}}$$

(84)

in expression (70) for ${\tilde{T}}_1$, after which the integral ${\tilde{T}}_1$ acquires the form

$${\tilde{T}}_1={\displaystyle \frac{na}{c\sqrt{1-e^2}}}\int \sqrt{1+e^2+2ecosf}df$$

(85)

$$=i{\displaystyle \frac{2na(1+e)}{cq\sqrt{1-e^2}}}\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-\frac{{\tilde{y}}^2}{q^2}\right)}}}$$

$$={\displaystyle \frac{2na(1+e)}{cq\sqrt{1-e^2}}}\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(\frac{{\tilde{y}}^2}{q^2}-1\right)}}}$$

(86)

$$={\displaystyle \frac{2na(1+e)}{icq\sqrt{1-e^2}}}\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-\frac{{\tilde{y}}^2}{q^2}\right)}}}={\displaystyle \frac{2na(1+e)}{icq\sqrt{1-e^2}}}{\tilde{J}}_2^{\left(4\right)}(\tilde{y};{\displaystyle \frac{1}{q}})$$

(87)

Note that the coefficient in front of the integral is modified in comparison with the one in Eq. (78). More importantly, the resulting expression is again a real-valued one due to the property of the elliptic integral

$${\tilde{J}}_2^{\left(4\right)}(\tilde{y},{\displaystyle \frac{1}{q}})=-i\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(\frac{{\tilde{y}}^2}{q^2}-1\right)}}},$$

(88)

which is easily established after comparing expressions (86) and (87). So it can be concluded that ${\tilde{T}}_1$in fact is expressed in two different variables - the first one is (78)

$${\tilde{T}}_1=-2i{\displaystyle \frac{na}{c}}{q}^{{\displaystyle \frac{3}{2}}}\int {\displaystyle \frac{y^{2}dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}=-2i{\displaystyle \frac{na}{c}}{q}^{{\displaystyle \frac{3}{2}}}{\tilde{J}}_2^{\left(4\right)}(y;q),$$

where the variable y is expressed as (76) $y=\sqrt{{\displaystyle \frac{1}{q}}{\displaystyle \frac{(1+ecosE)}{(1-ecosE)}}}$ , $q={\displaystyle \frac{1-e}{1+e}}$, and the second one is (86) in terms of the variable (84) $\tilde{y}={\displaystyle \frac{\sqrt{1+2ecosf+e^2}}{1+e}}$.

Formulae (88) also is valid, because due to the inequality $cosf\le 1$ and the choice of the variable $\tilde{y}$ in Eq. (84), it can easily be proved that

$$\tilde{y}\le 1,⊳1-{\tilde{y}}^2\ge 0,1-{\displaystyle \frac{{\tilde{y}}^2}{q^2}}\ge -{\displaystyle \frac{(1-q^2)}{q^2}}=-{\displaystyle \frac{4e}{{(1-e)}^2}}.$$

(89)

Consequently, since $1-{\displaystyle \frac{{\tilde{y}}^2}{q^2}}$ can take negative values (note that $q^2={\left({\displaystyle \frac{1-e}{1+e}}\right)}^2<1$) and thus $0<{\displaystyle \frac{{\tilde{y}}^2}{q^2}}-1\le {\displaystyle \frac{(1-q^2)}{q^2}}$, the representation (88) gives a real-valued integral $\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(\frac{{\tilde{y}}^2}{q^2}-1\right)}}}$ and thus an imaginary integral ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},{\displaystyle \frac{1}{q}})$ is obtained.

Let us present another more elegant and simple proof (see also the review paper [22]) that the inequalities (89) are fulfilled. Inverting the sign of the last inequality in (89), it means that the inequality

$${\displaystyle \frac{{\tilde{y}}^2}{q^2}}-1\le {\displaystyle \frac{4e}{{(1-e)}^2}}$$

(90)

should be proved. Keeping in mind the definition for the variable $\tilde{y}$ (84) $\tilde{y}={\displaystyle \frac{\sqrt{1+2ecosf+e^2}}{1+e}}$and also for $q^2$, the inequality (90) is transformed to

$${\displaystyle \frac{(1+2ecosf+e^2)}{{(1+e)}^2}}-{\displaystyle \frac{{(1-e)}^2}{{(1+e)}^2}}$$

$$\le {\displaystyle \frac{4e}{{(1-e)}^2}}{q}^2={\displaystyle \frac{4e}{{(1+e)}^2}}.$$

(91)

From here it follows

$$1+2ecosf+e^2\le 4e+{(1-e)}^2\Rightarrow 2ecosf\le 2e,$$

(92)

which is fulfilled for all eccentricities e, because of the simple trigonometric inequality $cosf\le 1$ for the function cos and also because $0<e<1$.

The other case, which has to be considered is when $1-{\displaystyle \frac{{\tilde{y}}^2}{q^2}}$ takes only positive values

$$1-{\displaystyle \frac{{\tilde{y}}^2}{q^2}}\ge 0\ge -{\displaystyle \frac{(1-q^2)}{q^2}}$$

(93)

and consequently ${\displaystyle \frac{{\tilde{y}}^2}{q^2}}-1\le 0$. However, in view of the future calculations, we will be interested in the case when the transformation $q\to {\displaystyle \frac{1}{q}}$ is applied in (88), so that the integral ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$ is obtained. Then from (93) after the transformation $q\to {\displaystyle \frac{1}{q}}$ it follows that ${\tilde{y}}^2{q}^2-1\le 0$ or $1-{\tilde{y}}^2{q}^2\ge 0$. Thus one will have the integral transformation

$$\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(\frac{{\tilde{y}}^2}{q^2}-1\right)}}}\stackrel{q\to {\displaystyle \frac{1}{q}}}{⟹}\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}}},$$

(94)

where the integral in the right is real-valued. Multiplying both integrals by ${\displaystyle \frac{1}{i}}$ and remembering the definition (88) for ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},{\displaystyle \frac{1}{q}})$, one obtains

$${\tilde{J}}_2^{\left(4\right)}(\tilde{y},{\displaystyle \frac{1}{q}})\stackrel{q\to {\displaystyle \frac{1}{q}}}{⟹}{\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)={\displaystyle \frac{1}{i}}\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}}},$$

(95)

which means that ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$ is imaginary valued. After presenting the calculation for the second $O\left(c^{-3}\right)$ propagation time correction in the next sections, it will become clear that due to the important transformation (95) it will be possible to prove the real-valuedness of terms with second-order elliptic integrals in the final expression for the propagation time for the case of space-oriented orbits.

Since in this final expression there will be also a term with a fourth-order elliptic integral, first some basic definitions and facts will be presented about higher-order elliptic integrals.

6. Definitions of Elliptic Integrals of Higher Order and Comparison with Some Statements from Standard Textbooks

6.1. Higher-Order Elliptic Integrals - Basic Definitions

This definition is necessary to be given because further, when calculating the second part of the propagation time, we shall encounter elliptic integrals of the fourth-order.

According to the general definition, elliptic integrals are of the type

$$\int R(y,\sqrt{G\left(y\right)})dy\mathrm{or}\int R(x,\sqrt{P\left(x\right)})dx,$$

(96)

where $R(\tilde{y},y)$ and $R(\tilde{x},x)$ are rational functions of the variables $\tilde{y},y$ or $\tilde{x},x$ and $\tilde{y}=\sqrt{G\left(y\right)}$, $\tilde{x}=\sqrt{P\left(x\right)}$ are square roots of the arbitrary polynomials $G\left(y\right)$ and $P\left(x\right)$ of the fourth or of the third degree respectively. In accord with the notations in the monograph [44], the fourth-order and the third-order polynomials will be of the form

$${\tilde{y}}^2=G\left(y\right)=a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4,$$

$${\tilde{x}}^2=P\left(x\right)=a{x}^3+b{x}^2+cx+d.$$

(97)

So in the general case, elliptic integrals of the $n-$th order and of the first kind are defined as

$$J_n^{\left(4\right)}\left(y\right)=\int {\displaystyle \frac{y^{n}dy}{\sqrt{a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4}}},$$

$$J_n^{\left(3\right)}\left(x\right)=\int {\displaystyle \frac{x^{n}dx}{\sqrt{a{x}^3+b{x}^2+cx+d}}}.$$

(98)

Elliptic integrals of the $n-$th order and of the third kind are defined as

$$H_n^{\left(4\right)}\left(y\right)=\int {\displaystyle \frac{dy}{{(y-c_1)}^n\sqrt{a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4}}},$$

$$H_n^{\left(3\right)}\left(x\right)=\int {\displaystyle \frac{dx}{{(x-c_2)}^n\sqrt{a{x}^3+b{x}^2+cx+d}}}.$$

(99)

In the above definitions, the upper indices "$\left(3\right)$" or "$\left(4\right)$" denote the order of the algebraic polynomial under the square root in the denominator, while the lower indice "n" denotes the order of the elliptic integral, related to the "$y^n$" or "$x^n$" terms in the nominator of (98) and in the denominator of (99) (for n-th order elliptic integrals of the third kind).

6.2. Comparison with Some Definitions from Standard Textbooks

In the textbook [65] a statement is expressed (although without any proof) that "higher-order elliptic integrals in the Legendre form can be expressed by means of the three standard integrals

$$\int {\displaystyle \frac{dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}},\int {\displaystyle \frac{y^{2}dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}},$$

$$\int {\displaystyle \frac{dy}{(1+h{y}^2)\sqrt{(1-y^2)(1-q^2{y}^2)}}},$$

(100)

where the parameter h can be also a complex number". As asserted in [65], these integrals "cannot be expressed through elementary functions". This assertion is not precise, because in the preceding sections an example was found, when the second integral in (100) in fact can be expressed through an elementary function - for this case, by the logarithmic term in Eq. (75). A similar statement to the one in [65] is given also in the monograph [64]: "Integrals of the type (96) $\int R(y,\sqrt{G\left(y\right)})dy$ and $\int R(x,\sqrt{P\left(x\right)})dx$ cannot be expressed by elementary functions in their final form, even if an extended understanding of this notion is considered."

In the contemporary monograph [44] a statement is expressed in a correct and precise manner, namely: "Every elliptic integral can be represented in the form of a linear combination of a rational function of the variables y and $\tilde{y}$, integrals of a rational function of the variable y and also the integrals

$$\int {\displaystyle \frac{dy}{\sqrt{G\left(y\right)}}},\int {\displaystyle \frac{ydy}{\sqrt{G\left(y\right)}}},$$

$$\int {\displaystyle \frac{y^{2}dy}{\sqrt{G\left(y\right)}}},\int {\displaystyle \frac{dy}{(y-c)\sqrt{G\left(y\right)}}}.\text{"}$$

(101)

In this statement it is not mentioned that elliptic integrals cannot be expressed by elementary functions. One new analytical algorithm for expressing elliptic integrals in the Legendre form will be proposed in the next sections. Since the monographs [65] and [64] are older monographs, this example clearly shows that some basic facts about elliptic functions and integrals cannot be considered to be established forever and the theory may change in time. One basic and contemporary monograph on elliptic functions is [67] (see especially ch. 8, dedicated to elliptic integrals). A short summary of the various types of elliptic integrals and functions is given in Appendix A of the book [68], as well as applications of elliptic functions in some problems of biophysics.

6.3. Elliptic Integrals Versus Abelian Integrals

The contemporary definition of the elliptic integrals (98) for $J_n^{\left(4\right)}\left(y\right)$ and $J_n^{\left(3\right)}\left(x\right)$ and (99) for $H_n^{\left(4\right)}\left(y\right)$ and $H_n^{\left(3\right)}\left(x\right)$ in the monograph [44] is fully consistent with the definition in older monographs, such as [45]. The elliptic integrals (98) and (99) cannot be considered to be a generalization of the s.c. abelian integrals, defined in the monograph [64] as

$$\int R(x,\sqrt[m]{{\displaystyle \frac{\alpha x+\beta }{\gamma x+\delta }}})dx,\int R(x,\sqrt{a{x}^2+bx+c})dx.$$

(102)

The reason, according to the definition in [64] is that the integrals (102) are considered to be related to algebraic curves of the type

$$\left(\gamma x+\delta )y^m-(\alpha x+\beta \right)=0,$$

(103)

$$y^2-(a{x}^2+bx+c)=0,$$

(104)

obtained after setting up

$$y=\sqrt[m]{{\displaystyle \frac{\alpha x+\beta }{\gamma x+\delta }}},y=\sqrt{a{x}^2+bx+c}.$$

(105)

It has been pointed out in [64] that abelian integrals of the type

$$\int {\displaystyle \frac{P\left(x\right)dx}{\sqrt[.]{a{x}^2+bx+c}}},\int {\displaystyle \frac{Adx}{{(x-\alpha )}^k\sqrt[.]{a{x}^2+bx+c}}},$$

(106)

$$\int {\displaystyle \frac{\left(Mx+N\right)dx}{{(x^2+px+q)}^m\sqrt[.]{a{x}^2+bx+c}}}$$

(107)

($P\left(x\right)$-polynomial, $A,M,N$-numerical constants ) can be treated analytically, the second integral for example by means of the known Euler substitutions. Evidently, this is possible since the degree of the polynomial under the square root in the denominator of these integrals is two.

Now let us define the elliptic integrals of the first-, second- and the third- kind and of the zero - order (see ch. 8 of the monograph by Armitage and Eberlein [67] and also [44]), which are partial cases of the already defined integrals (98) and (99)

$$\int {\displaystyle \frac{dx}{\sqrt[.]{\left(1-x^2\right)\left(1-q^2{x}^2\right)}}}=\underset{0}{\overset{\varphi }{\int }}{\displaystyle \frac{d\varphi }{\sqrt[.]{1-q^2{sin}^2\varphi }}}=F(q,\varphi ),$$

(108)

$$\int {\displaystyle \frac{\sqrt[.]{1-q^2{x}^2}}{\sqrt[.]{1-x^2}}}dx=\int \sqrt[.]{1-q^2{sin}^2\varphi }d\varphi =E(q,\varphi ),$$

(109)

$$\int {\displaystyle \frac{dx}{(1+n{x}^2)\sqrt[.]{\left(1-x^2\right)\left(1-q^2{x}^2\right)}}}=\int {\displaystyle \frac{d\varphi }{(1+n{sin}^2\varphi )\sqrt[.]{1-q^2{sin}^2\varphi }}}=\Pi (q,\varphi ).$$

(110)

In (108) $F(q,\varphi )$ is the commonly accepted notation for an elliptic integral of the first kind, in (109) $E(q,\varphi )$ is the notation for an elliptic integral of the second kind and in (110) $\Pi (q,\varphi )$ is the notation for the elliptic integral of the third kind. In the above equalities the transformation $x=sin\varphi$ has been used.

In the same way as with the algebraic curves (103) and (104) and setting up

$$y={\displaystyle \frac{\sqrt[.]{1-q^2{x}^2}}{\sqrt[.]{1-x^2}}},$$

(111)

$$y=(1+n{x}^2)\sqrt[.]{\left(1-x^2\right)\left(1-q^2{x}^2\right)},$$

(112)

it is easily seen that these algebraic curves are different from the curves (103) and (104), consequently the elliptic integrals (108), (109), (110) are not abelian integrals. Similarly, it can be proved that the found integral (48)

$$T={\displaystyle \frac{a}{c}}\int \sqrt[.]{{\displaystyle \frac{a_1{y}^3+a_2{y}^2+a_{3}y}{b_1{y}^3+b_2{y}^2+b_{3}y+b_4}}}dy$$

is not also an abelian integral. It remains also to find the corresponding integral for the case of space-oriented orbits and without the approximation $\beta ={\displaystyle \frac{2V}{c^2}}\ll 1$.

6.4. An Earlier Example for an Elliptic Integral, Expressed Also Analytically

However, from the formulaes in the preceding section it cannot be concluded that only abelian integrals can be treated analytically (in the sense of being expressed by elementary functions) and the elliptic integrals cannot be expressed by analytical functions. A typical example was given earlier, related with the calculation of the integral (70)

$${\tilde{T}}_1={\displaystyle \frac{na}{c\sqrt{1-e^2}}}\int \sqrt{1+e^2+2ecosf}df,$$

which was proved also to be the elliptic integral (86) ${\tilde{T}}_1={\displaystyle \frac{2na(1+e)}{cq\sqrt{1-e^2}}}\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(\frac{{\tilde{y}}^2}{q^2}-1\right)}}}$ of second order in the Legendre form in terms of the variable (84) $\tilde{y}={\displaystyle \frac{\sqrt{1+2ecosf+e^2}}{1+e}}$.

Also, in terms of the variables (76)

$$y=\sqrt[.]{{\displaystyle \frac{(1+ecosE)}{q(1-ecosE)}}},(q={\displaystyle \frac{1-e}{1+e}})$$

and (77)

$$tan{\displaystyle \frac{f}{2}}=\sqrt[.]{{\displaystyle \frac{1-cosf}{1+cosf}}}=\sqrt[.]{{\displaystyle \frac{1+e}{1-e}}}tan{\displaystyle \frac{E}{2}},$$

the integral ${\tilde{T}}_1$ was expressed as an elliptic integral of the second order and in the Legendre form (78)

$${\tilde{T}}_1=-2i{\displaystyle \frac{na}{c}}{q}^{{\displaystyle \frac{3}{2}}}\int {\displaystyle \frac{y^{2}dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}=-2i{\displaystyle \frac{na}{c}}{q}^{{\displaystyle \frac{3}{2}}}{\tilde{J}}_2^{\left(4\right)}(y;q).$$

Earlier by means of the transformations (71) and (72) it was proved that the same integral can be expressed analytically by means of the integrals (73)

$${\tilde{T}}_1=-{\displaystyle \frac{nai}{c\sqrt{(1-e^2)q}}}\left[\int {\displaystyle \frac{d\overline{m}}{\left(\overline{m}-\frac{q^2}{2}\right)\sqrt{\overline{m}-\frac{q^2}{2}}}}+\int {\displaystyle \frac{d\overline{m}}{\left(\overline{m}+\frac{q^2}{2}\right)\sqrt{\overline{m}-\frac{q^2}{2}}}}\right],$$

which can be integrated. Thus it can be concluded that when expressed in two different variables $\tilde{y}$ and $\overline{m}$, one and the same integral (70) can be considered to be an elliptic one, but also it can be written in another form and even calculated exactly. From the transformations (71) and (72) the relation between the variables f and $\overline{m}$ can be found to be

$$cosf=1+{\displaystyle \frac{1}{2e}}{\displaystyle \frac{(2\overline{m}+q^2)B_1}{(2\overline{m}-q^2)}}.$$

(113a)

If account is taken also of (84) for $\tilde{y}$, the relation between the variables $\tilde{y}$ and $\overline{m}$ is calculated to be a quadratic algebraic curve with respect to the variable $\tilde{y}$

$${\tilde{y}}^2=1+{\displaystyle \frac{1}{{(1+e)}^2}}.{\displaystyle \frac{(2\overline{m}+q^2)B_1}{(2\overline{m}-q^2)}}.$$

(114)

Then if the integral (70) in terms of the variable $\overline{m}$ can be integrated, it can be asked if there is any motivation for the assertion in the monographs [64,65,67] that elliptic integrals cannot be integrated in elementary functions? Further it shall be proved that analytical formulaes can be found for the determination of the s.c. Weierstrass invariants $g_2$ and $g_3$, depending on complicated polynomials of the modulus q of the elliptic integral in the Legendre form. But then the s.c."four-dimensional parametrization" of an elliptic curve will allow the parametrization of an fourth-order polynomial in terms of expressions, depending ot the Weierstrass function $\rho \left(z\right)$ and its derivative ${\displaystyle \frac{d\rho \left(z\right)}{dz}}$, which of course are not elementary functions. This will be proved further, but evidently depending on the variables used in the calculations, a certain elliptic integral can be expressed analytically in elementary functions, but also can be expressed in an analytical form, containing non-elementary functions.

7. Fourth-Order Elliptic Integrals and the Second $O\left({\displaystyle \frac{1}{c^3}}\right)$ Correction for the Propagation Time for the Case of a Space-Oriented Orbit

Taking into account the general equality Eq. (69) $\overline{T}={\tilde{T}}_1+{\tilde{T}}_2={\displaystyle \frac{1}{c}}\int vdt-{\displaystyle \frac{2}{c^3}}\int vVdt$ , the second $O\left({\displaystyle \frac{1}{c^3}}\right)$ correction can be written as

$${\tilde{T}}_2=-\int {\displaystyle \frac{v}{c}}.{\displaystyle \frac{2V}{c^2}}dt=-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}\int {\displaystyle \frac{v_f}{r}}df,$$

(115)

where the radius-vector r in the orbital plane, the velocity $v_f$, associated with the true anomaly angle f and the potential V of the gravitational field around the Earth are defined in the usual way

$$v_f={\displaystyle \frac{na}{\sqrt{1-e^2}}}\sqrt{1+2ecosf+e^2},r={\displaystyle \frac{a(1-e^2)}{1+ecosf}},$$

$$V={\displaystyle \frac{GM}{r}}={\displaystyle \frac{GM}{a(1-e^2)}}(1+ecosf).$$

(116)

The final integral for the second correction is of the form

$${\tilde{T}}_2=-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}.{\displaystyle \frac{na}{a\sqrt{1-e^2}}}\int (1+ecosf)\sqrt{1+2ecosf+e^2}df=T_2^{\left(1\right)}+T_2^{\left(2\right)},$$

(117)

where

$$T_2^{\left(1\right)}=-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}.{\displaystyle \frac{n}{{(1-e^2)}^{\frac{3}{2}}}}\int \sqrt{1+2ecosf+e^2}df$$

$$=-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}.{\displaystyle \frac{n}{{(1-e^2)}^{\frac{3}{2}}}}{\tilde{T}}_1$$

(118)

and ${\tilde{T}}_1$ is the previously calculated integral ${\tilde{T}}_1=\int \sqrt{1+2ecosf+e^2}df$ . The second term $T_2^{\left(2\right)}$ in the second correction $T_2$ is

$$T_2^{\left(2\right)}=-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}.{\displaystyle \frac{ne}{{(1-e^2)}^{\frac{3}{2}}}}\int cosf\sqrt{1+2ecosf+e^2}df$$

$$=-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}{\displaystyle \frac{ne}{{(1-e^2)}^{\frac{3}{2}}}}{\tilde{T}}_2^{\left(2\right)},$$

(119)

where ${\tilde{T}}_2^{\left(2\right)}$ is the notation for the more complicated integral

$${\tilde{T}}_2^{\left(2\right)}=\int cosf\sqrt{1+2ecosf+e^2}df.$$

(120)

7.1. The Second Part of the Second $O\left({\displaystyle \frac{1}{c^3}}\right)$ Correction, Expressed in Terms of Elliptic Integrals of the Second and Fourth Order

It can be proved that

$$T_2^{\left(2\right)}=-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}.{\displaystyle \frac{ne}{{(1-e^2)}^{\frac{3}{2}}}}{\tilde{T}}_2^{\left(2\right)}$$

(121)

$$=-in{q}^{{\displaystyle \frac{5}{2}}}{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}(1+e^2){\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$$

$$+in{q}^{{\displaystyle \frac{3}{2}}}{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}{\displaystyle \frac{(1+e^2)}{(1-e^2)}}{\tilde{J}}_4^{\left(4\right)}(\tilde{y},q),$$

(122)

where

$${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)=\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}}}$$

$$={\displaystyle \frac{1}{i}}\int {\displaystyle \frac{q^3{\widehat{y}}^{2}d\widehat{y}}{\sqrt{\left({\widehat{y}}^2-1\right)\left(1-q^2{\widehat{y}}^2\right)}}}={\displaystyle \frac{q^3}{i}}{\tilde{J}}_2^{\left(4\right)}(\widehat{y},q)$$

(123)

is an elliptic integral of the second order in the Legendre form and $\widehat{y}$ is the notation for the variable $\widehat{y}={\displaystyle \frac{\tilde{y}}{q}}$. The last formulae can also be compared to (88). Analogously, the formulae for the fourth-order elliptic integral ${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)$ can be written as

$${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)={\displaystyle \frac{q^5}{i}}\int {\displaystyle \frac{{\widehat{y}}^{4}d\widehat{y}}{\sqrt{\left({\widehat{y}}^2-1\right)\left(1-q^2{\widehat{y}}^2\right)}}}.$$

(124)

7.2. Relation Between the Fourth-Order and the Second-Order Elliptic Integrals in the Expression for $T_2^{\left(2\right)}$

It can very easily be proved that

$$\int \sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}d\tilde{y}={\displaystyle \frac{2}{3}}{\tilde{J}}_0^{\left(4\right)}$$

$$-{\displaystyle \frac{1}{3}}(1+q^2){\tilde{J}}_2^{\left(4\right)}(\tilde{y},q))+{\displaystyle \frac{1}{3}}\left[\tilde{y}\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}\right]{\mid }_{\tilde{y}={\tilde{y}}_0}^{\tilde{y}={\tilde{y}}_1}.$$

(125)

We have already proved that ${\tilde{J}}_2^{\left(4\right)}$ can be expressed in elementary functions, but yet we do not know whether ${\tilde{J}}_0^{\left(4\right)}$ or the second-rank integral

$$\int \sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}d\tilde{y}$$

(126)

can also be expressed through elementary functions. This will be proved further.

Calculating the derivatives

$${\displaystyle \frac{d}{d\tilde{y}}}\left(\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}\right)\mathrm{and}{\displaystyle \frac{d}{d\tilde{y}}}\left(\tilde{y}\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}\right),$$

(127)

integrating the resulting equations from ${\tilde{y}}_0$ to ${\tilde{y}}_1$, and combining all the three equations, the following two equations can be derived for ${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)$ and ${\tilde{J}}_3^{\left(4\right)}(\tilde{y},q)$:

$${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)={\displaystyle \frac{1}{3q^2}}\left[\tilde{y}\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}\right]{\mid }_{\tilde{y}={\tilde{y}}_0}^{\tilde{y}={\tilde{y}}_1}$$

$$+{\displaystyle \frac{2(1+q^2)}{3q^2}}{\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)-{\displaystyle \frac{1}{3q^2}}{\tilde{J}}_0^{\left(4\right)}(\tilde{y},q),$$

(128)

$${\tilde{J}}_3^{\left(4\right)}(\tilde{y},q)={\displaystyle \frac{1}{2q^2}}\left[\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}\right]{\mid }_{\tilde{y}={\tilde{y}}_0}^{\tilde{y}={\tilde{y}}_1}+{\displaystyle \frac{(1+q^2)}{2q^2}}{\tilde{J}}_1^{\left(4\right)}.$$

(129)

In expressions (125) - (129) the symbol "${\mid }_{\tilde{y}={\tilde{y}}_0}^{\tilde{y}={\tilde{y}}_1}$" means that the corresponding expression is taken at the value $\tilde{y}={\tilde{y}}_1$ of the upper boundary and from it the value of the expression at $\tilde{y}={\tilde{y}}_0$(the value at the lower boundary) is substracted.

We should also add to this system of equations the earlier derived equation for ${\tilde{J}}_2^{\left(4\right)}$

$${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)=\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{(1-{\tilde{y}}^2)(1-q^2{\tilde{y}}^2)}}}$$

$$=-{\displaystyle \frac{1}{q^2}}E\left(\phi \right)+{\displaystyle \frac{1}{q^2}}{\tilde{J}}_0^{\left(4\right)}({\tilde{y}}^2;q).$$

(130)

The second equation (129) for ${\tilde{J}}_3^{\left(4\right)}$ clearly suggests that ${\tilde{J}}_3^{\left(4\right)}$ is expressed in elementary functions because the first-order elliptic integral ${\tilde{J}}_1^{\left(4\right)}$ $=\int {\displaystyle \frac{\tilde{y}d\tilde{y}}{\sqrt{(1-{\tilde{y}}^2)(1-q^2{\tilde{y}}^2)}}}$ can be analytically calculated after performing the Euler substitution

$$\sqrt{{\tilde{y}}^4-{\displaystyle \frac{(1+q^2)}{q^2}}{\tilde{y}}^2+{\displaystyle \frac{1}{q^2}}}=u+{\tilde{y}}^2.$$

(131)

The result is

$${\tilde{J}}_1^{\left(4\right)}=-{\displaystyle \frac{1}{2q^3}}ln\mid \sqrt{(1-{\tilde{y}}^2)(1-q^2{\tilde{y}}^2)}-{\tilde{y}}^2+{\displaystyle \frac{(1+q^2)}{2q^2}}{\mid }_{\tilde{y}={\tilde{y}}_0}^{\tilde{y}={\tilde{y}}_1}.$$

(132)

Note that the variable $\tilde{y}$ is not related to any elliptic function! Concerning the fourth-order elliptic integral ${\tilde{J}}_4^{\left(4\right)}$, taking into account the already proved fact that ${\tilde{J}}_2^{\left(4\right)}$ can be expressed in elementary functions, it would follow that ${\tilde{J}}_4^{\left(4\right)}$ can also be expressed by elementary functions, but it is necessary to prove that the zero-order elliptic integral ${\tilde{J}}_0^{\left(4\right)}$ can be expressed in elementary functions. In another publication, it will be proved that the integral can be expressed by an analytical formulae, but also it can be numerically calculated, since the modulus of the elliptic integral is equal to the eccentricity $e=0.01323881349526$ of the GPS orbit, which is exactly known [72]. The general mathematical theory for the recurrent relations between $n-$th order elliptic integrals of any kind and the lower-order elliptic integrals has been developed in the monographs [44,63] and [64].

7.3. The General Expression for the Propagation Time for the Case of Space-Oriented Orbits

Now we shall collect all the separate terms, which constitute the expression for the propagation time T. These terms are (86), earlier proved to be a real-valued integral

$${\tilde{T}}_1={\displaystyle \frac{2na(1+e)}{cq\sqrt{1-e^2}}}\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(\frac{{\tilde{y}}^2}{q^2}-1\right)}}},$$

the expression (118) for $T_2^{\left(1\right)}$

$$T_2^{\left(1\right)}=-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}.{\displaystyle \frac{n}{{(1-e^2)}^{\frac{3}{2}}}}{\tilde{T}}_1,$$

expression (122) for $T_2^{\left(2\right)}$

$$T_2^{\left(2\right)}=-in{q}^{{\displaystyle \frac{5}{2}}}{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}(1+e^2){\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$$

$$+in{q}^{{\displaystyle \frac{3}{2}}}{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}{\displaystyle \frac{(1+e^2)}{(1-e^2)}}{\tilde{J}}_4^{\left(4\right)}(\tilde{y},q),$$

expression (123) for ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$

$${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)=\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}}}$$

$$={\displaystyle \frac{q^3}{i}}\int {\displaystyle \frac{{\widehat{y}}^{2}d\widehat{y}}{\sqrt{\left({\widehat{y}}^2-1\right)\left(1-q^2{\widehat{y}}^2\right)}}}={\displaystyle \frac{q^3}{i}}{\tilde{J}}_2^{\left(4\right)}(\widehat{y},q),$$

expression (124) for $J_4^{\left(4\right)}(\tilde{y},q)$

$${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)={\displaystyle \frac{q^5}{i}}\int {\displaystyle \frac{{\widehat{y}}^{4}d\widehat{y}}{\sqrt{\left({\widehat{y}}^2-1\right)\left(1-q^2{\widehat{y}}^2\right)}}}.$$

Note that this expression for ${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)$ resembles the previous expression ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$, since both of them have one and the same denominators $\sqrt{\left({\widehat{y}}^2-1\right)\left(1-q^2{\widehat{y}}^2\right)}$, which determine whether both formulaes are imaginary- or real- valued. Earlier by formulaes (85) - (89) and (93), (94), (95) it was proved that ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$ is an imaginary-valued expression. Consequently, it can be concluded that expression (124) for ${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)$ is also imaginary-valued.

This conclusion can be made also from the chain of elliptic integrals, establishing the relation between elliptic integrals of zero-, second- and forth- order (128)

$${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)={\displaystyle \frac{1}{3q^2}}\left[\tilde{y}\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}\right]{\mid }_{\tilde{y}={\tilde{y}}_0}^{\tilde{y}={\tilde{y}}_1}$$

$$+{\displaystyle \frac{2(1+q^2)}{3q^2}}{\tilde{J}}_2^{\left(4\right)}-{\displaystyle \frac{1}{3q^2}}{\tilde{J}}_0^{\left(4\right)}.$$

If the integrals ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$ and ${\tilde{J}}_0^{\left(4\right)}(\tilde{y},q)$ are imaginary-valued, then since this is related to the imaginary -valuedness of the expression $\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}$ (the first term in (128)), then the left-hand side of (128) will also be imaginary valued. This precludes the proof that ${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)$ is also imaginary-valued. Note also that if the variable $\tilde{y}$ in ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$ and ${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)$ is replaced by y, the integrals ${\tilde{J}}_2^{\left(4\right)}(y,q)$ and ${\tilde{J}}_4^{\left(4\right)}(y,q)$ will also be imaginary valued.

7.3.1. The Entire Expression for the Propagation Time of the Signal for the Case of Space-Oriented Orbits

Making use of all the above mentioned formulaes, the total propagation time (both the $O\left({\displaystyle \frac{1}{c}}\right)$ and $O\left({\displaystyle \frac{1}{c^3}}\right)$ corrections) can be written as

$$\overline{T}=\left[-2i{\displaystyle \frac{na}{c}}{q}^{{\displaystyle \frac{3}{2}}}+4i{\displaystyle \frac{G_{\oplus }{M}_{\oplus }{n}^{2}a}{c^4{(1-e^2)}^{\frac{3}{2}}}}{q}^{{\displaystyle \frac{3}{2}}}\right]{\tilde{J}}_2^{\left(4\right)}(y,q)$$

$$-i{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }na{q}^{\frac{5}{2}}(1+e^2)}{c^3}}{\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)+i{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }n{q}^{\frac{3}{2}}(1+e^2)}{c^3(1-e^2)}}{\tilde{J}}_4^{\left(4\right)}(\tilde{y},q),$$

(133)

where all the integrals ${\tilde{J}}_2^{\left(4\right)}(y,q)$, ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$ and ${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)$ (76), written in terms of the variable $y=\sqrt{{\displaystyle \frac{1}{q}}{\displaystyle \frac{(1+ecosE)}{(1-ecosE)}}}$ or (84) $\tilde{y}={\displaystyle \frac{\sqrt{1+2ecosf+e^2}}{1+e}}$ are proved to be imaginary-valued. This precludes the proof that the propagation time $\overline{T}$ is a real-valued expression, because all the imaginary-valued integrals are multiplied by imaginary coefficients.

7.4. Real-Valuedness and Complex-Valuedness of Elliptic Integrals of Zero - Order in the Legendre Form - Basic Knowledge about the Christoffel-Schwartz Integral from Complex Analyses

It has been proved that the properties of higher-order elliptic integrals (especially of the second- and of the fourth- order) are very important in order to find analytically the propagation time for a signal, emitted by satellites and influenced by the gravitational field of the Earth. At the same time, as already mentioned, higher-order integrals are investigated only in the monographs [44,63] and [64] and are neglected in many other monographs. For the moment, what is known for such integrals is that they are obtained from the recurrent chain of elliptic integrals (39) and (40), in which the zero-order elliptic integral ${\tilde{J}}_0^{\left(4\right)}(\tilde{y},q)$ participates. So whether the fourth-order integral ${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)$ (previously written in formulaes (125) and (128)) and the second-order integral ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$ (written in formulae (130)) are real-valued or imaginary-valued will depend also on the properties of ${\tilde{J}}_0^{\left(4\right)}(\tilde{y},q)$, because it can be real- or imaginary- valued.

Zero-order elliptic integrals in the Legendre form such as $J_0^{\left(4\right)}(y,q)=\int {\displaystyle \frac{dy}{\sqrt{\left(1-y^2\right)\left(1-q^2{y}^2\right)}}}$ are analyzed by means of the well-known Christoffel - Schwartz integral (see especially [63,69,70,71], but also many other monographs on analytical functions)

$$w=f\left(z\right)=C\underset{0}{\overset{z}{\int }}{(s-a_1)}^{{\alpha }_1-1}{(s-a_2)}^{{\alpha }_2-1}.......{(s-a_n)}^{{\alpha }_n-1}ds+C_1,$$

(134)

which represents a conformal mapping of the upper complex half-plane onto the inner part of an $n-$polygon. In the last formulae $a_1$, $a_2$,$......a_n$ are points on the real axis and ${\alpha }_1$, ${\alpha }_2....{\alpha }_n$ denote the inner angles, represented by real numbers. Each one of these numbers is not greater than 2 and for them the following equality

$${\alpha }_1+{\alpha }_2+....+{\alpha }_n=n-2$$

(135)

is fulfilled. For the four-dimensional case of the zero-order integral $J_0^{\left(4\right)}(y,q)$, formulae (134) represents a mapping of the upper complex half-plane onto the rectangle (i.e. the $4-$polygon, which for the case turns out to be the rectangle). So the integral $J(y,q)$ is a partial case of the integral (134) for the following values of the parameters ${\alpha }_n$, $a_n$ and the constants $C_1$ and C:

$${\alpha }_1={\alpha }_2={\alpha }_3={\alpha }_4={\displaystyle \frac{1}{2}},a_1=-{\displaystyle \frac{1}{q}},a_2=-1,$$

(136)

$$a_3=1,a_4={\displaystyle \frac{1}{k}},C_1=0,C={\displaystyle \frac{1}{q}}.$$

(137)

Therefore, the integral $J_0^{\left(4\right)}(y,q)$ can be represented in the form of the following Christoffel-Schwartz integral

$$J_0^{\left(4\right)}(y,q)=f\left(z\right)={\displaystyle \frac{1}{q}}\underset{0}{\overset{z}{\int }}{\left(y+{\displaystyle \frac{1}{q}}\right)}^{{\displaystyle \frac{1}{2}}-1}{\left(y+1\right)}^{{\displaystyle \frac{1}{2}}-1}{\left(y-1\right)}^{{\displaystyle \frac{1}{2}}-1}{\left(y-{\displaystyle \frac{1}{q}}\right)}^{{\displaystyle \frac{1}{2}}-1}dy.$$

(138)

At the points $a_1$, $a_2$, $a_3$, $a_4$ the under-integral function in the integral $J_0^{\left(4\right)}(y,q)$ (similarly - also in the integral (134)) is a branching multi-valued function, which maps the points on the real axis onto the points on the segments of the rectangle, situated on the complex plane. Each expression ${(s-a_k)}^{{\alpha }_k-1}$ in (134) should be interpreted as the branch of the multi-valued function which on the real axis (when $z=s>a_k$) takes real positive values. Indeed, if s belongs to the segment $[0,1]$, then the function $f\left(z\right)$ in (138) takes values from zero to the value of the elliptic integral of the first kind:

$$K=\underset{0}{\overset{1}{\int }}{\displaystyle \frac{ds}{\sqrt[.]{\left(1-s^2\right)\left(1-q^2{s}^2\right)}}}.$$

(139)

In fact, this is an equivalent formulation of the period of the elliptic integral $K=F({\displaystyle \frac{\pi }{2}},e)=\underset{0}{\overset{{\displaystyle \frac{\pi }{2}}}{\int }}{\displaystyle \frac{d\overline{E}}{\sqrt[.]{1-e^2{sin}^2\overline{E}}}}$. Since the function $\left(1-s^2\right)\left(1-q^2{s}^2\right)$ is a two-valued function at the branching points $s=\pm 1$ and $s=\pm {\displaystyle \frac{1}{q}}$, one has to choose the under-integral function so that to ensure the continuity of this function, when $1<s<{\displaystyle \frac{1}{q}}$. This means that for this case it should be written as [69]

$${\displaystyle \frac{1}{\pm i\sqrt[.]{\left(s^2-1\right)\left(1-q^2{s}^2\right)}}}.$$

(140)

Consequently, for $s\subset [1,{\displaystyle \frac{1}{q}}]$ the function $w=f\left(z\right)$ (138) maps the interval $[1,{\displaystyle \frac{1}{q}}]$ into the interval $[K,K+i{\tilde{K}}^{,}]$ on the complex plane, where

$${\tilde{K}}^{,}=\underset{1}{\overset{{\displaystyle \frac{1}{q}}}{\int }}{\displaystyle \frac{ds}{\sqrt[.]{\left(s^2-1\right)\left(1-q^2{s}^2\right)}}}$$

(141)

and K is the integral (139). Then the integral in (138) can be written as:

$$\underset{0}{\overset{y}{\int }}{\displaystyle \frac{ds}{\sqrt[.]{\left(1-s^2\right)\left(1-q^2{s}^2\right)}}}=K+i{\tilde{K}}^{,}=$$

(142)

$$=\underset{0}{\overset{1}{\int }}{\displaystyle \frac{ds}{\sqrt[.]{\left(1-s^2\right)\left(1-q^2{s}^2\right)}}}+i\underset{1}{\overset{y}{\int }}{\displaystyle \frac{ds}{\sqrt[.]{\left(s^2-1\right)\left(1-q^2{s}^2\right)}}}.$$

(143)

In the same way,similarly the analysis can be performed for the other intervals on the real axis. Thus the following theorem about the Christoffel-Schwartz integral can be proved.

Theorem 1.

The intervals $[0,1]$, $[1,{\displaystyle \frac{1}{q}}]$, $[{\displaystyle \frac{1}{q}},+\infty )$ on the real axis are mapped by the Christoffell - Schwarz integral (138) onto the rectangle with endpoints correspondingly $(0,0)$, $(K,0)$, $(K,K+i{\tilde{K}}^{,})$, $(0,i{\tilde{K}}^{,})$ on the complex plane, where ${\tilde{K}}^{,}$ is represented by the integral (141) and K is the integral (139).

The first significance of this theorem is that since the mentioned points are the points of a rectangle, it gives the opportunity to predict when the integral (138) is real-valued or imaginary.

Secondly, it may be noted that the Christoffel-Schwarz theorem is not formulated for elliptic integrals of higher than zero order. So taking into account the results of the theorem for the zero-order integral ${\tilde{J}}_0^{\left(4\right)}(\tilde{y},q)$ and combining the higher-order integrals (125), (128) for the integral ${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)$ and for the integral (130) ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$ respectively from the recurrent system of equations, it may be established in another way whether these integrals are real-valued or complex-valued and how the interval $[0,1]$, $[1,{\displaystyle \frac{1}{q}}]$, $[{\displaystyle \frac{1}{q}},+\infty )$ is mapped by higher-order integrals. This is a problem, not solved for the moment in the theory of elliptic functions. Another problem, not solved is: how to find the period of an elliptic integral of higher order?

8. Numerical Calculation of the Propagation Time of a Signal, Emitted by a Satellite on a Plane Elliptical Orbit

8.1. Numerical Calculation of the First Six Iterations of the Eccentric Anomaly Angle E

We shall perform some numerical calculations of the propagation time, based on the derived formulae (54). The numerical data about the parameters of the $GPS$ orbit are taken from the PhD dissertation [72] and are known with great precision. From all the parameters, listed below only the first three will be used in the calculation, performed by the online program web2.0 scientific calculator [74]

semi-major axis a$26560.25169632944$$\left[km\right]$ ,

eccentricity e$0.01323881349526$ ,

mean anomaly M$-0.3134513508155$$\left[rad\right]$ ,

inclination I$0.9614884100802$$\left[rad\right]$ ,

longitude of the ascending node $\Omega$$-0.4495096737336$$\left[rad\right]$ ,

argument of perigee $\omega$$-3.001488651204$$\left[rad\right]$ .

In the following numerical calculations we shall use only the first three parameters $a,e$ and M, but in other publications the approach will be extended and it will be required to find the propagation time dependent on the full set of 6 Kepler parameters $(a,e,M,I,\Omega ,\omega )$. Instead of M, either the eccentric anomaly angle E will be used or the true anomaly angle f. In fact, all the Kepler parameters might be used for finding the propagation time, and this is a problem to be solved in future.

For finding the propagation time for the case of a signal, emitted by a satellite, moving along a plane elliptical orbit with an eccentricity e, most important is to find the eccentric anomaly angle E from the Kepler equation $M=E-esinE$. The calculation will be based on formulaes $(2.54)$ and $(2.55)$ in Ch. 2 of the monograph [41].

The iterative sequence of the formulaes are

$$E_{(i+1)}=M+esin{E}_{\left(i\right)},i=0,1,2,.....,$$

(144)

where for the first approximation it is assumed $E_0=M=-0.3134513508155$$\left[rad\right]$ The first three iterative solutions are given according to the following formulaes:

$$E_{\left(1\right)}=M+esinM,$$

(145)

$$E_{\left(2\right)}=M+esin{E}_{\left(1\right)}=M+esin(M+esinM),$$

(146)

$$E_{\left(3\right)}=M+esin{E}_{\left(2\right)}=M+esin[M+esin(M+esinM)].$$

(147)

Thus from the formulaes for $E_{\left(3\right)}$, $E_{\left(5\right)}$ and $E_{\left(6\right)}$ the following numerical values can be found

$$E_{\left(3\right)}=M+esin{E}_{\left(2\right)}=-0.31758547588467897473\left[rad\right],$$

(148)

$$E_{\left(4\right)}=M+esin{E}_{\left(3\right)}=-0.31758548401096719083\left[rad\right],$$

(149)

$$E_{\left(5\right)}=M+esin{E}_{\left(4\right)}=-0.31758548411317083102\left[rad\right],$$

(150)

$$E_{\left(6\right)}=M+esin{E}_{\left(5\right)}=-0.31758548411445499592\left[rad\right].$$

(151)

The approximation $E_{\left(5\right)}$ up to the ninth digit is identical with $E_{\left(4\right)}$ and $E_{\left(6\right)}$ up to the eleventh digit is identical with $E_{\left(5\right)}$. The approximation $E_{\left(6\right)}$ also up to the seventh digit is identical with $E_{\left(3\right)}$, given by (148). However, if the approximate value for $E_{\left(3\right)}$ is compared with the value for $E_{\left(6\right)}$ , the coincidence is only up to the fifth digit.

It is therefore necessary to check what is the impact of the approximations in the eccentric anomaly E on the approximations of the time of propagation of the signal.

8.2. Numerical Calculation of the First $O\left({\displaystyle \frac{1}{c}}\right)$ Correction in the Propagation Time

The corresponding elliptic integral (61) of the second kind for the eccentric anomaly approximation $E_{\left(3\right)}$ is

$$T_1^{\left(E_3\right)}=\underset{0}{\overset{E_{\left(3\right)}}{\int }}\sqrt[.]{1-e^2{cos}^{2}E}dE=-0.317557268125933936045.$$

(152)

If we compare this value with the value (148) $E_{\left(3\right)}=-0.31758547588467897473$$\left[rad\right]$ , then it can be noted that the integration changes the value of $E_{\left(3\right)}$ after the fourth digit after the decimal dot.

The same integral with $E_{\left(6\right)}$ as an upper integration limit is

$$T_1^{\left(E_6\right)}=\underset{0}{\overset{E_{\left(6\right)}}{\int }}\sqrt[.]{1-e^2{cos}^{2}E}dE=-0.317558568963886638536.$$

(153)

This value is different from the value (152) after the fifth digit, so it is a better approximation. However, if compared with $E_{\left(6\right)}=$$-0.317585484114454995929$$\left[rad\right]$, the integration in (153) changes the value of $E_{\left(6\right)}$ after the fourth digit after the decimal dot.

The calculation of the first $O\left({\displaystyle \frac{1}{c}}\right)$ time correction${\displaystyle \frac{a}{c}}{T}_1$ for the values of the $GPS$ orbit and for the eccentric anomaly $E_{\left(3\right)}$ gives the following numerical value

$${\displaystyle \frac{a}{c}}{T}_1^{\left(E_3\right)}=-0.0281341332790419\left[sec\right].$$

(154)

Correspondingly, the first $O\left({\displaystyle \frac{1}{c}}\right)$ time correction ${\displaystyle \frac{a}{c}}{T}_1^{\left(E_6\right)}$ for the eccentric anomaly $E_{\left(6\right)}$ is

$${\displaystyle \frac{a}{c}}{T}_1^{\left(E_6\right)}=-0.0281342485273829\left[sec\right].$$

(155)

So the two $O\left({\displaystyle \frac{1}{c}}\right)$ time corrections (154) and (155) are identical up to the sixth digit after the decimal dot. Consequently, the important conclusion is that the eccentric anomaly at the sixth approximation level can ensure microsecond stability of the $O\left({\displaystyle \frac{1}{c}}\right)$ time correction. This means also that the third approximation $E_{\left(3\right)}$ can be reliably used, at least at the microsecond level.

One more argument in support of the microsecond approximation is that for a satellite, moving along circular orbit (when $e=0$) the above expression acquires the form:

$${\displaystyle \frac{a}{c}}{t}_{1\left(circular\right)}^{\left(E_3\right)}={\displaystyle \frac{a}{c}}{E}_{\left(3\right)}=-0.028136517830159266[sec]$$

$$=-28136.517830159266[\mu sec].$$

(156)

Evidently, the coordinate time for the circular orbit differs from the coordinate time for the $GPS$ elliptic orbit only after the fifth digit after the decimal dot. This is however a clear proof that in order to account for changes in the propagation time in the micro- ($\mu$) and the nano-range, the ellipticity of the orbit should be accounted in the calculation of the propagation time according to

8.3. Numerical Calculation of the Shapiro Delay Term (the $O\left({\displaystyle \frac{1}{c^3}}\right)$ Time Correction)

Now let us calculate numerically the whole under-integral expression in formulae (54) for the propagation time, taking into account the $25-$th digit after the decimal dot

$${\displaystyle \frac{a}{c}}\sqrt{1-e^2{cos}^{2}E}-{\displaystyle \frac{2GM}{c^3}}\sqrt[.]{{\displaystyle \frac{1+ecos{E}_{\left(3\right)}}{1-ecos{E}_{\left(3\right)}}}}$$

$$=-0.0885884561709019072629880-0.0000000000299618094618168$$

(157)

$$=-0.0885884562008637167248048.$$

(158)

The second Shapiro delay term contains ten zeroes after the decimal dot, so the overall result of the calculation of both terms in (157) is changed by the Shapiro term at and after the nanosecond level (1 $nsec={10}^{-9}$sec). Here it should be clarified that the calculation in (157) and in (158) is of an approximate value because one should take the numerical values for T after performing the integration of the elliptic integral, which as already mentioned changes the result after the fourth number after the decimal dot. For the moment, we do not have a numerical estimation of the second integral term in (54), consequently we take the whole under-integral expression.

It should be noted also that in the previous calculations of the iterations $E_{\left(3\right)}$ and $E_{\left(6\right)}$ of the eccentric anomaly angle it was taken with a negative sign, due to the accepted convention. After obtaining the result of the calculation of $T_1^{\left(E_3\right)}$ and $T_1^{\left(E_6\right)}$ from the integrals (152) and (153), they were with a negative sign. This means that after integrating of the first term in (157), there will be a negative sign before the first number in (157) and the negative sign of the second term in (158) remains. This is so because after calculating ${\displaystyle \frac{2GM}{c^3}}\sqrt[.]{{\displaystyle \frac{1+ecos{E}_{\left(3\right)}}{1-ecos{E}_{\left(3\right)}}}}$, the obtained number is positive. Thus, the two numbers in (157) and (158) should be added, giving a total number $-0.0885884562008637167248048$, greater than the first number in (157) after the ninth digit after the decimal dot.

Although approximate, this result, based on the simple application of the celestial mechanics approach in the calculation of the propagation time is interesting, if compared with the calculations, based on the formulae (1) $T={\displaystyle \frac{R_{AB}}{c}}+2{\displaystyle \frac{G_{\oplus }{M}_E}{c^3}}ln\left({\displaystyle \frac{r_A+r_B+R_{AB}}{r_A+r_B-R_{AB}}}\right)$. From the literature it can be seen that relativistic effects on light propagation from Satellite Laser Ranging (SLR) data are measured with an accuracy of 1 $\mu m$ (1 micrometer), which corresponds to $0.01$$psec$ (see also paper [17]). The relativistic observables, defined for the $GRAIL$ mission in [3] are also accurate to 1 $\mu m$. No doubt that the numerical calculation of the integral over the Shapiro delay term in (157) will give a more realistic result.

8.4. Propagation Time of the Signal Compared to the Celestial Time - Consistency of the Numerical Calculation

Let us compare the calculated propagation time of the signal with the celestial time $t_{cel}$, which can be calculated from $M=n(t-t_P)$ (M is the mean anomaly, given in the numerical data), where $t_P$ is the initial time of perigee passage. The idea is to calculate the celestial time of motion for the given value M of the mean anomaly and thus to compare the distance, travelled by the satellite with the distance, travelled by the light (laser or radio) signal. Since this calculation also will be approximate, for the propagation time we have taken only the value $T_1^{\left(E_6\right)}$ for the propagation time from the integral (153) and (155)

$$T_1^{\left(E_6\right)}={\displaystyle \frac{a}{c}}\underset{0}{\overset{E_{\left(6\right)}}{\int }}\sqrt[.]{1-e^2{cos}^{2}E}dE=-0.0281342485273829[sec],$$

which does not account for the second term, responsible for the Shapiro delay. Let us take the ratio of the propagation time to the celestial time

$${\displaystyle \frac{\mathrm{propagation}\mathrm{time}\mathrm{for}\mathrm{electromagnetic}\mathrm{signals}}{\mathrm{celestial}\mathrm{time}\mathrm{from}\mathrm{the}\mathrm{Kepler}\mathrm{equation}}}={\displaystyle \frac{0.0281342485273829\left[sec\right]}{37.508256148\left[sec\right]}}$$

(159)

$$=0.0007500814864964887.$$

(160)

In other words, the propagation time for the signal is about ${10}^4$ times smaller than the celestial time, which will mean that for this celestial time $37.5082561\left[sec\right]$ the satellite moves $145.125\left[km\right]$ (the velocity of the satellite is taken to be $v_S=3.874$$[km/sec]$) and the light signal will move at a distance

$$26558.1510169176263505735093268096\left[km\right].$$

(161)

This value is obtained after multiplying the light velocity $c=299792.458$$\left[{\displaystyle \frac{km}{sec}}\right]$ with the calculated propagation time, taking into account the sixth iteration $T_1^{\left(E_6\right)}$ according to (155).

It is interesting to compare this value for the distance, travelled by the signal with the value, obtained after taking only the first term in (157) for the propagation time, where the term ${\displaystyle \frac{a}{c}}\sqrt{1-e^2{cos}^{2}E}$ has not been integrated. This term does not account for the second term for the Shapiro delay, which contained ten zeroes after the decimal dot and was of the order $0.299$$nsec$ (nanoseconds). The result will be

$$26558.151025899950855259224944504\left[km\right].$$

(162)

This number is greater (after the fourth digit after the decimal dot) than the preceding number (161), so the integration of the expression (155)

$${{\displaystyle \frac{a}{c}}}_0^{E_{\left(6\right)}}\sqrt[.]{1-e^2{cos}^{2}E}dE$$

in fact diminishes the propagation time (only the first $O\left({\displaystyle \frac{1}{c}}\right)$ correction) in comparison with the first term ${\displaystyle \frac{a}{c}}\sqrt{1-e^2{cos}^{2}E}$ in the expression (157). The difference between the numbers (162) and (161) is

$$0.0000089823245046857156176944\left[km\right]\approx 0.898\left[cm\right]\approx 8.98\left[mm\right].$$

(163)

In both cases of the numerical calculation of the numbers (161) and (162), the numerical result confirms that the distance, travelled by the light signal (numbers (161) and (162)) is much greater than the distance of $145.125\left[km\right]$, travelled by the satellite. Respectively, for the calculated propagation time $0.0281342485273829$$[sec]$ (only the first $O\left({\displaystyle \frac{1}{c}}\right)$ correction $T_1^{\left(6\right)}$ (155)) the satellite will move at a distance

$${\displaystyle \frac{145.125\left[km\right]}{37.5082561\left[sec\right]}}\times 0.0281342485273829\left[sec\right]$$

(164)

$$=0.1088555758671073849925\left[km\right].$$

(165)

This distance is considerably smaller than the distance, travelled by the signal. In fact, it is about $108.8$$meters$, which is $0.000004098763343794755826290507\sim 0.409\times {10}^{-5}$ times smaller than the distance, given by the number (161)

$$26558.151016917626350573509326809\left[km\right],$$

travelled by the light or radio signal for the calculated propagation time of $0.0281342485273829$$[sec]$.

8.5. Comparison of the Result in Formulae (163) With A Result in Other Papers. Approximations in the Calculations of the Shapiro Delay Term

The result of $8.98$$\left[mm\right]$ in (163) was obtained by substracting the propagation distance, obtained from the first term in the formula (54), first by taking into account the integrated term and afterwards - by calculating the numerical value only for the under-integral expression (see also (157)).

However, in the known review article by Ashby [34] it was proved that at the level of several millimeters, spatial curvature effects have to be considered, which is evident from the simple equality

$$\int dr\left[1+{\displaystyle \frac{G{M}_E}{c^{2}r}}\right]\approx r_2-r_1+{\displaystyle \frac{G{M}_E}{c^2}}ln{\displaystyle \frac{r_2}{r_1}}.$$

(166)

In this paper the second term was estimated to be $4.43\times ln(4.2)\approx 6.3$$mm$. It might be concluded that yet in the first term in (163) spatial curvature effects are important and have to be included. Since in the Shapiro formulae the first term is the Euclidean distance, divided by the light velocity and from a formal point of view it might be believed that the first terms ${\displaystyle \frac{a}{c}}\sqrt{1-e^2{cos}^{2}E}$ in (157) and ${\displaystyle \frac{a}{c}}\underset{0}{\overset{E_{\left(6\right)}}{\int }}\sqrt[.]{1-e^2{cos}^{2}E}dE$ in (155) should correspond to the first term ${\displaystyle \frac{R}{c}}$ in the Shapiro formulae. Therefore, since it might be suspected that curvature effects may appear yet in these formulaes, it is not quite correct to consider any correspondence ${\displaystyle \frac{a}{c}}\sqrt{1-e^2{cos}^{2}E}\Rightarrow {\displaystyle \frac{R}{c}}$ or ${\displaystyle \frac{a}{c}}\underset{0}{\overset{E_{\left(6\right)}}{\int }}\sqrt[.]{1-e^2{cos}^{2}E}dE\Rightarrow {\displaystyle \frac{R}{c}}$. In view of this, in a recent paper [23] a new modification of the Shapiro formulae has been proposed, where the first term ${\displaystyle \frac{R_{.}}{c}}$ is replaced by the term ${\displaystyle \frac{1}{c}}\underset{path}{\int }\sqrt{{\left(\stackrel{.}{x}\right)}^2+{\left(\stackrel{.}{y}\right)}^2+{\left(\stackrel{.}{z}\right)}^2}ds$. In this term s is some parameter, which parametrizes some curve, for example the space curve, determining the space coordinates on the satellite orbit.

It should be reminded also that the integration in formulae (166) is of approximate character. The reason is that some approximations are made when calculating the second Shapiro delay term in (1)

$$T={\displaystyle \frac{R_{AB}}{c}}+2{\displaystyle \frac{G_{\oplus }{M}_E}{c^3}}ln\left({\displaystyle \frac{r_A+r_B+R_{AB}}{r_A+r_B-R_{AB}}}\right).$$

For example, in the thesis [73] the approximation

$$cdt=\left(1+{\displaystyle \frac{G{M}_E}{c^{2}r}}\right)\Vert dr\Vert$$

(167)

has been used, where the infinitesimal vector $dr$ is not related to the signal trajectory. The vector $\overrightarrow{r}$ in (167) is related to the vector $\overrightarrow{R}=\overrightarrow{r}-{\overrightarrow{r}}_A$, where ${\overrightarrow{r}}_A$ is the vector, associated with the position of emission and $\overrightarrow{r}$ - with the position of a space point in a Geocentric Reference Frame. In other words, the integration is not along $\Vert dR\Vert$, because the approximation $\Vert dR\Vert \approx \Vert dr\Vert$ has been assumed and thus the integration is along $dr$. A similar approximation has been used in the PhD thesis [72].

That is why, in the papers [19,20,21,22] and [23] it has been asserted that the derivation of the Shapiro formulae should start from the null cone equation.

8.6. Comparison of the Propagation Time with the Atomic Time of the Atomic Clocks of the Satellites in the Near-Earth Space - Consistency of the Results

Now let us point out one another numerical fact, which is interesting, because it confirms the consistency of the approach of calculating the propagation time on the base of using the parametrization of the satellite orbit. In the book [38] it was pointed out that at an altitude of 20184 $km$, due to the difference in the gravitational potential, the atomic time of the satellite clock runs faster by 45 $\mu sec/d$ (microseconds per day, 1 $\mu sec={10}^{-6}$sec). Consequently, for one second the atomic time will run faster by $0.5208333.{10}^{-9}$$[sec]$. In the preceding sections it was proved that the dominant part in the propagation time $dT$ (i.e. the part without the Shapiro delay term, which is very small) is $dT=0.0281341332790419$$\left[sec\right]$, given by the number (155). The corresponding to this propagation time interval $dT$ atomic time is $d\tau =0.0146531934.{10}^{-9}$$[sec]$ and can easily be found from 45 $\mu sec/d$ (keeping in mind how many seconds are contained in 24 hours, i.e. 1 day). The ratio of the two time intervals is equal to

$${\displaystyle \frac{d\tau }{dT}}={\displaystyle \frac{0.0146531934.{10}^{-9}}{0.028134248527382}}=0.520831163687866092.{10}^{-9}.$$

(168)

The very small atomic time interval compared to the propagation interval means that the atomic time can serve as a standard for measuring the propagation time, because it will be able to detect changes even at the nanosecond level.

It should be noted that this value is of the order of ${10}^{-9}$, but it is a little greater than the value $\beta ={\displaystyle \frac{2V}{c^2}}=0.334\times {10}^{-9}$.

Similarly, if a clock on an orbit is compared to a clock on the Earths surface, due to the net effect of time dilation and gravitational redshift the satellite clock will appear to run fast by 38 $\mu sec/d$ [38], which also is considered to be an enormous difference for an atomic clock with precision of a few nanoseconds. So for this value, a similar number will be obtained for the ratio ${\displaystyle \frac{d\tau }{dT}}$, again confirming that this nanosecond $({10}^{-9}$$sec)$ precision of the atomic clock (no matter whether it is on the Earths surface or on a satellite in a near-Earth orbit) is sufficient for measuring the propagation time T of the signal.

9. Justification of the Theoretical Approach in This Paper from the Viewpoint of General Relativity Theory and Differential Geometry

As mentioned in the Introduction, the basic approach of calculation of the propagation time in this paper is based on finding the solution for the propagation time from the null cone equation and parametrizing the space coordinates for both cases, considered of plane elliptical orbits or space-oriented orbits by the orbital coordinates. So in fact, the light-like geodesics are not used in this approach. In the previous chapter, it was demonstrated that the solution for the propagation time does not give any time, related to the motion of the satellite, but the calculated propagation time is of the order of microseconds (1 $\mu sec={10}^{-6}$sec) or nanoseconds (1 $nsec={10}^{-9}$sec), which is the "starting" sensitivity of the atomic clocks. This is the first fact in favour of the correctness of using the approach of the null cone equation. The second fact is related to the proper dimension of seconds of the expression for the propagation time for the case of a signal, emitted by a satellite on a plane elliptical orbit. The third fact is very important, in view of the application of elliptic integrals of the first-, second and the third kind (for the case of a satellite on a plane elliptical orbit) or elliptic integrals of the first-, second- and the fourth- order (for the case of a satellite on a space-distributed orbit).

Now we shall present the fourth argument, related to a theorem, proved in the book on General Relativity Theory by Fock [29], which turns out to be related with known theorems from basic textbooks on differential geometry [32] and [30]. Let us first begin with the theorem in [29]

Theorem 2.

Let $F={\displaystyle \frac{1}{2}}{g}_{\alpha \beta }{\displaystyle \frac{d{x}^{\alpha }}{ds}}$${\displaystyle \frac{d{x}^{\beta }}{ds}}$ , $F=0$ is the null cone equation and the extremum of the integral $s=\underset{s_1}{\overset{s_2}{\int }}Lds$ is searched $(L=\sqrt{2F})$, then the Lagrange equation of motion will give

$${\displaystyle \frac{dF}{ds}}=0⟹F=const⟹{\displaystyle \frac{d}{ds}}{\displaystyle \frac{\partial F}{\partial {\stackrel{.}{x}}_{\alpha }}}-{\displaystyle \frac{\partial F}{\partial x_{\alpha }}}=0.$$

(169)

From the last equation the geodesic equations $(\alpha ,\beta ,\nu =0,1,2,3)$

$${\displaystyle \frac{d^2{x}_{\nu }}{d{s}^2}}+{\Gamma }_{\alpha \beta }^{\nu }{\displaystyle \frac{d{x}_{\alpha }}{ds}}{\displaystyle \frac{d{x}_{\beta }}{ds}}=0$$

(170)

can be obtained. Since the constant in (169) can be set up to zero, this would mean that the geodesic equation (170) is compatible with the null-cone equation $F={\displaystyle \frac{1}{2}}{g}_{\alpha \beta }{\displaystyle \frac{d{x}^{\alpha }}{ds}}$${\displaystyle \frac{d{x}^{\beta }}{ds}}=0$.

It is important to acquire the proper understanding of this statement - it means that every solution of the null equation (including in the proper parametrization of the differentials, related to the trajectory of motion of the satellite) is also a solution of the light-like geodesic equation (170). In the "backward direction", however the statement is not true, meaning that there might be another solutions of the geodesic equation, which are not solutions of the null-cone equation.

The meaning of the geodesic as the shortest path, joining two given points (but only in the local sense, i.e a sufficiently short part of the geodesics) is confirmed also by the following theorem in the monograph by Prasolov [32]. Let us present the exact formulation of this theorem:

Theorem 3.

Any point p on a surface S has a neighborhood U such that the length of a geodesic γ going from the point p to some point q and lying entirely in the neighborhood U does not exceed the length of any curve α on S, joining the points p and q. Moreover, if the lengths of the curves γ and α are equal, then γ and α coincide as non-parametrizable curves.

In fact, what does it mean that "the length of the geodesics from point p to point q does not exceed the length of any other curve" (but only locally, i.e. in the vicinity of the given neighborhood)? This means that this geodesic line is minimal. So this confirms the theorem in another, third well-known book on differential geometry by Mishtenko and Fomenko [30]:

Theorem 4.

The geodesic line $\gamma :\left[a,b\right]\to M^n$ is minimal if it is not longer than any smooth path, joining the endpoints $\gamma \left(a\right)$ and $\gamma \left(b\right)$.

So if we have for example the functional $E\left(\gamma \left(s\right)\right)=\int \mid {\stackrel{.}{\gamma }}^2\left(s\right)\mid ds=\int g_{ij}\left(x\right){\stackrel{.}{x}}^i\left(s\right){\stackrel{.}{x}}^j\left(s\right)ds$ for the null cone equation, and we are searching for the extremals, these extremals will be derived from the geodesics - the multitude of the solutions for the extremals from the geodesics turns out to be larger in comparison with the one for the other local curves in the vicinity, in particular for the functional $E\left(\gamma \right(s\left)\right)$ for the null-cone equation.

On the base of combining the results in the previously mentioned monographs by Fock [29] and [32], one easily establishes the validity of the following theorem in [30]:

Theorem 5.

The extremals of the length functional $L\left(\gamma \left(s\right)\right)=\int \sqrt[.]{g_{ij}\left(x\right){\stackrel{.}{x}}^i\left(s\right){\stackrel{.}{x}}^j\left(s\right)}ds$ are smooth functions, derived from the geodesics (by means of smooth change of the parameter along the curve). Then each extremal of the functional $E\left(\gamma \left(s\right)\right)=\int \mid {\stackrel{.}{\gamma }}^2\left(s\right)\mid ds=\int g_{ij}\left(x\right){\stackrel{.}{x}}^i\left(s\right){\stackrel{.}{x}}^j\left(s\right)ds$ (called also action functional of the trajectory) is also an extremal of the functional $L\left(\gamma \right(s\left)\right)=$$\int \sqrt[.]{g_{ij}\left(x\right){\stackrel{.}{x}}^i\left(s\right){\stackrel{.}{x}}^j\left(s\right)}ds$. However, each extremal of the functional $L\left(\gamma \right(s\left)\right)$ is not necessarily an extremal of the functional $E\left(\gamma \right(s\left)\right)$.

There is a simple inequality $L^2\left(\gamma \left(s\right)\right)\le$$E\left(\gamma \right(s\left)\right)$ [30] in confirmation of the above theorem, which can be immediately proved by means of the Schwartz inequality

$$\left(\underset{0}{\overset{1}{\int }}f\left(s\right)g\left(s\right)ds\right)\le \left(\underset{0}{\overset{1}{\int }}{f}^2\left(s\right)ds\right)\left(\underset{0}{\overset{1}{\int }}{g}^2\left(s\right)ds\right)$$

(171a)

for the following functions [30]

$$f\left(s\right)\equiv 1g\left(s\right)\equiv \mid \stackrel{.}{\gamma }\left(s\right)\mid .$$

(172)

It should be noted that the parameter s in the book [29] is defined simply as a "space parameter along a curved line". In the above formulaes, s is also a space parameter. But as noted in [30], when $g\left(s\right)=const$, i.e. $\gamma \left(s\right)=const$ and the parameter s is proportional to the length of the arc of the curved line, an equality can be derived $L^2\left(\gamma \left(s\right)\right)=E\left(\gamma \left(s\right)\right)$. Thus, every solution of the null cone equation $g_{ij}\left(x\right){\stackrel{.}{x}}^i\left(s\right){\stackrel{.}{x}}^j\left(s\right)=0$ (i.e. a zero-extremal of $E\left(\gamma \right(s\left)\right)$) will also be a zero extremal of $L^2\left(\gamma \left(s\right)\right)$. As clarified in the introduction, if this zero-extremal corresponds to the propagation of the signal along the geodesics of the minimal length, the propagation time will satisfy the Fermats principle about the least time of propagation of light between two space points. Consequently, even if there are other extremals for the length functional, of interest are only those, related to the minimal geodesics, obtained from the action functional $E\left(\gamma \left(s\right)\right)=\int g_{ij}\left(x\right){\stackrel{.}{x}}^i\left(s\right){\stackrel{.}{x}}^j\left(s\right)ds$, which ensure the fulfillment of the Fermat principle.

10. New Approach for Analytical Calculation of Elliptic Integrals. Analytic Relation Between Elliptic Integrals in the Weierstrass and in the Legendre Form

In this section some new analytical algorithms for calculation of elliptic integrals will be presented, which have been published in the paper [21]. The main advantage of the theoretical method is that 1. It will be applicable for various values of the modulus q of the elliptic integral, no matter whether q is a small or a big number. 2. The integrals, encountered in the calculation of the propagation time in the preceding sections were in the Legendre form, which contains the specific fourth-order polynomial $(1-y^2)(1-q^2{y}^2)$. However, in previous sections the definition for elliptic integrals with arbitrary polynomials of the third and of the fourth degree were given. Such integrals are very often encountered in cosmology, nonlinear evolution equations, black hole physics and etc. For example, in the definition for the luminocity in cosmology the more general polynomial of the fourth degree under the square root in the denominator will be used. Also, there is an interesting interplay between the methods of the standard three-dimensional uniformization and of the four-dimensional uniformization, described in the monograph by Whittaker and Watson [60]. This method will be exposed in details further and will allow not only the analytical computation of the Weierstrass invariants $g_2$ and $g_3$ for four-dimensional elliptic integrals of the type $J_n^{\left(4\right)}\left(y\right)$ , but also the comparison of elliptic integrals in the Weierstrass form $J_n^{\left(3\right)}\left(x\right)$, both of them defined by formulaes (96), (97), (98), (99) in section "Higher-order Elliptic Integrals-Basic Definitions".

Some physical applications of elliptic integrals were presented in the paper [21].

10.1. Transforming an Elliptic Integral in the Legendre Form into an Elliptic Integral in the Weierstrass Form

The purpose of this section will be to propose a transformation, which will transform the elliptic integral

$${\tilde{J}}_0^{\left(4\right)}\left(y\right)=\int {\displaystyle \frac{dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}$$

(173)

in the Legendre form into an elliptic integral in the Weierstrass form

$${\tilde{J}}_0^{\left(3\right)}\left(x\right)=-\sqrt{a}\int {\displaystyle \frac{dx}{\sqrt{4x^3-g_{2}x-g_3}}}.$$

(174)

For the purpose, the transformation

$$x={\displaystyle \frac{a}{y^2}}+b$$

(175)

shall be applied. It is possible to calculate the parameters a and b, so that the integral (174) is satisfied. These parameters represent complicated functions of the modulus parameter q of the integral

$$a={\displaystyle \frac{9g_3}{g_2}}K\left(q\right),$$

(176)

$$b=-{\displaystyle \frac{a}{3}}(1+q^2)=-{\displaystyle \frac{3g_3}{g_2}}K\left(q\right)$$

(177)

and $K\left(q\right)$ is a rational function, depending on the modulus parameter q of the elliptic integral in the Legendre form

$$K\left(q\right)\equiv {\displaystyle \frac{q^4-q^2+1}{2q^4-5q^2+2}}.$$

(178)

However, $g_2$ and $g_3$ remain undetermined. This will turn out to be an advantage of the applied formalism, since further another representation of the integral in the Weierstrass form shall be found with different Weierstrass invariants ${\overline{g}}_2$ and ${\overline{g}}_3$, which will depend explicitly on the modulus parameter q. Then by comparing both representations in the Weierstrass form and proving a new theorem, $g_2$ and $g_3$ will be expressed through ${\overline{g}}_2$ and ${\overline{g}}_3$ in a complicated way.

Further the following important property will be very important.

Corollary 1.

The parameter b satisfies the cubic equation, i.e.

$$4b^3-g_{2}b-g_3=0.$$

(179)

The proof is straightforward, if the expressions for $g_2$ and $g_3$ are expressed from (176) and (177) and are substituted into the above cubic equation.

10.2. Equivalence Between a Formulae from Integral Calculus and the Recurrent System of Equations for the Elliptic Integral in Terms of the $y-$Variable

Due to the property (179), the elliptic integral (174) can be represented in the equivalent way

$${\tilde{J}}_0^{\left(3\right)}\left(\overline{x}\right)=-{\displaystyle \frac{\sqrt{a}}{2}}\int {\displaystyle \frac{d\overline{x}}{{\overline{x}}^{\frac{1}{2}}\sqrt{{\overline{x}}^2+3b\overline{x}+(b^2-\frac{g_2}{4})}}},$$

(180)

where $\overline{x}\equiv x-b$. Such integrals in the monograph by Timofeev [66] are analytically computed according to the formulae

$$\int {\displaystyle \frac{d\overline{x}}{{\overline{x}}^m\sqrt{{\left(\tilde{a}+2\tilde{b}\overline{x}+\tilde{c}{\overline{x}}^2\right)}^{2n+1}}}}=-{\displaystyle \frac{\sqrt{{\left(\tilde{a}+2\tilde{b}\overline{x}+\tilde{c}{\overline{x}}^2\right)}^{2n-1}}}{(m-1)\tilde{a}{\overline{x}}^{m-1}}}$$

$$-{\displaystyle \frac{(2m+2n-3)\tilde{b}}{(m-1)\tilde{a}}}\int {\displaystyle \frac{d\overline{x}}{{\overline{x}}^{m-1}\sqrt{{\left(\tilde{a}+2\tilde{b}\overline{x}+\tilde{c}{\overline{x}}^2\right)}^{2n+1}}}}$$

$$-{\displaystyle \frac{(m+2n-2)\tilde{c}}{(m-1)\tilde{a}}}\int {\displaystyle \frac{d\overline{x}}{{\overline{x}}^{m-2}\sqrt{{\left(\tilde{a}+2\tilde{b}\overline{x}+\tilde{c}{\overline{x}}^2\right)}^{2n+1}}}}.$$

(181)

Since this is a formulae from integral calculus, it is interesting to prove the following theorem, which most surprisingly establishes the equivalence between the above relation from integral calculus with some of the recurrent system of equations for the case of a fourth-degree polynomial in the denominator of the elliptic integral [44]. The recurrent system for the case of a cubic polynomial under the square root has been investigated in the monographs of Fichtenholz [64] and Smirnov [63]. Below the formulation of the newly proved theorem is given.

Theorem 6.

The integral (181) is equivalent to the second equation

$${\overline{z}}^{.}\tilde{Z}=3\tilde{c}{J}_4^{\left(4\right)}\left(\overline{z}\right)+4\tilde{b}{J}_2^{\left(4\right)}\left(\overline{z}\right)+\tilde{a}{J}_0^{\left(4\right)}\left(\overline{z}\right).$$

(182)

from the recurrent system of equations

$${\overline{z}}^m\tilde{Z}=(m+2)\tilde{c}{J}_{m+3}^{\left(4\right)}\left(\overline{z}\right)+(m+1)2\tilde{b}{J}_{m+1}^{\left(4\right)}\left(\overline{z}\right)+m\tilde{a}{J}_{m-1}^{\left(4\right)}\left(\overline{z}\right),$$

(183)

where the fourth-degree polynomial $\tilde{Z}$ is defined as

$${\tilde{Z}}^2=a_0{\overline{z}}^4+4a_1{\overline{z}}^3+6a_2{\overline{z}}^2+4a_3\overline{z}+a_4.$$

(184)

The coefficients in the above polynomial have the values

$$a_0=\tilde{c},a_1=0,6a_2=2\tilde{b},a_3=0,a_4=\tilde{a}.$$

(185)

The corresponding elliptic integrals of the first kind and of the third kind and of the $n-$th order are

$$J_n^{\left(4\right)}\left(\overline{z}\right)=\int {\displaystyle \frac{{\overline{z}}^n}{\sqrt{\tilde{a}+2\tilde{b}{\overline{z}}^2+\tilde{c}{\overline{z}}^{4.}}}}d\overline{z}=\int {\displaystyle \frac{{\overline{z}}^{n}d\overline{z}}{\sqrt[.]{{\overline{z}}^4+3b{\overline{z}}^2+(3b^2-\frac{g_2}{4})}}},$$

(186)

$$H_n^{\left(4\right)}\left(\overline{z}\right)=\int {\displaystyle \frac{d\overline{z}}{{\overline{z}}^n\sqrt{\tilde{a}+2\tilde{b}{\overline{z}}^2+\tilde{c}{\overline{z}}^{4.}}}}=\int {\displaystyle \frac{d\overline{z}}{{\overline{z}}^n\sqrt[.]{{\overline{z}}^4+3b{\overline{z}}^2+(3b^2-\frac{g_2}{4})}}}.$$

(187)

Due to the established equivalence, it can be concluded that the integral (181) cannot be used for the calculation of the zero-order elliptic integral in the Legendre form (173), because the recurrent relation (182) contains also the higher-order elliptic integrals $J_2^{\left(4\right)}\left(\overline{z}\right)$ and $J_4^{\left(4\right)}\left(\overline{z}\right)$.

10.3. A New Theorem, Concerning the Equivalence between Some of the Equations in the Recurrent System

Theorem 7.

From (183), the corresponding equation for $m=-3$

$${\displaystyle \frac{\tilde{y}}{y^3}}=-q^2{J}_0\left(y\right)+2(1+q^2)H_2^{\left(4\right)}\left(y\right)-3H_4^{\left(4\right)}\left(y\right)$$

(188)

can be transformed into the equation for $m=1$

$$y{\tilde{y}}^{\left(q\right)}=J_0\left(y\right)-2(1+q^2)J_2\left(y\right)+3q^2{J}_4\left(y\right)$$

(189)

by means of the transformation $y\Rightarrow {\displaystyle \frac{1}{qy}}$ and the simple relations (for the concrete case for $m=2$)

$$H_m\left(y\right)=-q^m{J}_m\left({\displaystyle \frac{1}{qy}}\right),H_m\left({\displaystyle \frac{1}{qy}}\right)=-q^m{J}_m\left(y\right)$$

(190)

In such a way, only one independent equation can be obtained from the $4D$ recurrent system of equations

$$2(1+q^2)J_2^{\left(q\right)}\left(y\right)=3q^2{J}_4^{\left(q\right)}\left(y\right)+J_0^{\left(q\right)}\left(y\right)-y{\tilde{y}}^{\left(q\right)}.$$

(191)

10.4. Relations Between the Recurrent System of Equations in Terms of the y and x Variables

10.4.1. The Recurrent System of Equations in Terms of the x Variable

Now we shall write the recurrent system of equations for the elliptic integrals

$$I_m^{\left(3\right)}\left(x\right)\equiv \int {\displaystyle \frac{x^{m}dx}{\sqrt{4x^3-g_{2}x-g_3}}},$$

$$G_m^{\left(3\right)}\left(x\right)\equiv \int {\displaystyle \frac{dx}{x^m\sqrt{4x^3-g_{2}x-g_3}}},x={\displaystyle \frac{{\tilde{a}}_2}{y^2}}+{\tilde{b}}_2,$$

(192)

as this has been performed in the monographs [64] and [63]. Generally, the polynomial of a fourth-order under the square root is denoted as

$$Y^2=a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4,$$

(193)

but for the concrete case of the integrals (192), the coefficients will be given by

$$a_0=0,a_1=1,a_2=0,a_3=-{\displaystyle \frac{g_2}{4}},a_4=-g_3.$$

(194)

For this choice of the coefficients due to $a_0=0$ the polynomial becomes cubic, so we change the variables in (193) from $(Y,y)$ to $(X,x)$, i.e. $X^2=4x^3-g_{2}x-g_3$, where the Weierstrass invariants (coefficients) $g_2$ and $g_3$ will be defined in the next sections.

The recurrent system of equations is easily found to be (omitting the upper indice "$\left(3\right)$" in $I_m^{\left(3\right)}$)

$$x^{n}X=2(2n+3)I_{n+2}-{\displaystyle \frac{g_2}{2}}(2n+1)I_n-n{g}_3{I}_{n-1}.$$

(195)

It should be noted also that the variable transformation $x={\displaystyle \frac{{\tilde{a}}_2}{y^2}}+{\tilde{b}}_2$ in (192) contains another parameters ${\tilde{a}}_2$ and ${\tilde{b}}_2$, which are different from the parameters $a_2$ and $b_2$ in the transformation (175). Now we shall transform the three-degree polynomials under the square root in (192) into the known four-degree polynomials for the elliptic integral in the Legendre form. The transformation $x={\displaystyle \frac{{\tilde{a}}_2}{y^2}}+{\tilde{b}}_2$ (as well as the previously applied transformation (175) $x={\displaystyle \frac{a}{y^2}}+b$) are generalizations of the transformation $x=e_3+{\displaystyle \frac{(e_1-e_3)}{y^2}}$ in the paper [59], where $e_1,e_2,e_3$ are constant roots of a cubic polynomial. The transformation transforms the integral $\int {\displaystyle \frac{dx}{\sqrt{4(x-e_1)(x-e_2)(x-e_3)}}}$ into the integral $-{\displaystyle \frac{1}{\sqrt{e_1-e_3}}}\int {\displaystyle \frac{dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}$, where $q^2={\displaystyle \frac{(e_2-e_3)}{(e_1-e_3)}}$. Under another transformation $y\Rightarrow {\displaystyle \frac{1}{\tilde{k}y}}$, the $x-$variable will transform as

$$\tilde{x}\equiv {\tilde{a}}_2{\tilde{k}}^2{y}^2+{\tilde{b}}_2,$$

(196)

where a new variable $\tilde{x}$ has been introduced. The two recurrent system of equations (183) (written in terms of the $y-$variable) and the recurrent system (195) for the $x-$variable are not independent because

$$I_0^{\left(3\right)}\left(x\right)=-{\displaystyle \frac{J_0^{\left(q\right)}\left(y\right)}{{\left(\sqrt{a}\right)}^q}}={\displaystyle \frac{2{\tilde{a}}_2}{{\left(\sqrt{a}\right)}^q}}{J}_0^{\left(\tilde{k}\right)}\left(y\right),$$

(197)

where $\tilde{k}$ is determined so that $4x^3-g_2\left(q\right)x-g_3\left(q\right)$ is transformed into the expression $(1-y^2)(1-$${\tilde{k}}^2{y}^2)$. Consequently, the last expression is satisfied if $\tilde{k}$ is defined as

$$\tilde{k}=-1-{\displaystyle \frac{6\alpha }{\sqrt[3]{2}}},\alpha ={\tilde{b}}_2\mathrm{is}\mathrm{a}\mathrm{root}\mathrm{of}4{\tilde{b}}_2^3-g_2^{\left(q\right)}{\tilde{b}}_2-g_3^{\left(q\right)}=0,$$

$${\tilde{a}}_2={\displaystyle \frac{1}{\sqrt[3]{4}}}={\displaystyle \frac{{\tilde{k}}^4-{\tilde{k}}^2+1}{3g_2^{\left(q\right)}}}.$$

(198)

Two central problems remain to be solved:

1. How to combine the derived equations in terms of the x variable with the found relation (191)?

2. How to estimate the parameter $g_2$ in Eq. (198)? In fact, further in the next sections it will be proved that $g_2=g_2^{\left(q\right)}\equiv g_2\left(q\right)$ and $g_3=g_3^{\left(q\right)}\equiv g_3\left(q\right)$ will depend in a complicated way on the modulus parameter q, but this dependence will be found from another Weierstrass representation of the integral $J_0^{\left(4\right)}\left(y\right)$ (173). This new representation will not contain the conformal factor $\sqrt{a}$ in Eq. (174), but instead will be characterized by different Weierstrass invariants ${\overline{g}}_2^{\left(q\right)}$ and ${\overline{g}}_3^{\left(q\right)}$, which will be explicitly determined by the method of the s.c. "fourth-degree polynomial uniformization". This algebraic method will be considered in details in the following sections.

10.4.2. Combining the Two Recurrent System of Equations

From the two recurrent system of equations in terms of the $y-$variable and the $x-$variable, the following equation can be obtained

$${\tilde{X}}^{\left(\tilde{k}\right)}=N_1({\tilde{a}}_2,{\tilde{b}}_2,\tilde{k})J_2^{\left(\tilde{k}\right)}\left(y\right)+N_2({\tilde{a}}_2,{\tilde{b}}_2,\tilde{k})J_0^{\left(\tilde{k}\right)}\left(y\right)+N_3({\tilde{a}}_2,{\tilde{b}}_2,\tilde{k})=$$

(199)

$$=N_1^{\left(\tilde{k}\right)}{J}_2^{\left(\tilde{k}\right)}\left(y\right)+N_2^{\left(\tilde{k}\right)}\left[-{\displaystyle \frac{1}{2{\tilde{a}}_2}}{J}_0^{\left(q\right)}\left(y\right)\right]+N_3^{\left(\tilde{k}\right)},$$

(200)

where

$$N_1({\tilde{a}}_2,{\tilde{b}}_2,\tilde{k})\equiv N_1^{\left(\tilde{k}\right)}\equiv 8{\tilde{a}}_2^3{\tilde{k}}^2(1+{\tilde{k}}^2)+24{\tilde{a}}_2^2{\tilde{b}}_2{\tilde{k}}^2,$$

(201)

$$N_2({\tilde{a}}_2,{\tilde{b}}_2,\tilde{k},)\equiv N_2^{\left(\tilde{k}\right)}\equiv 12{\tilde{a}}_2{\tilde{b}}_2^2-{\displaystyle \frac{g_2{\tilde{a}}_2}{{\left(\sqrt{a}\right)}^{\left(q\right)}}}-4{\tilde{a}}_2^3{\tilde{k}}^2,$$

(202)

$$N_3({\tilde{a}}_2,{\tilde{b}}_2,\tilde{k})\equiv N_3^{\left(\tilde{k}\right)}\equiv 4{\tilde{a}}_2^3{\tilde{k}}^2{y}^2{\tilde{y}}^{\left(\tilde{k}\right)}$$

(203)

and ${\tilde{X}}^{\left(\tilde{k}\right)}\equiv 4{\tilde{x}}^3-g_2^{\left(\tilde{k}\right)}\tilde{x}-g_3^{\left(\tilde{k}\right)}$. The notations $g_2^{\left(\tilde{k}\right)}$ and $g_3^{\left(\tilde{k}\right)}$ mean that these Weierstrass invariants are expressed as complicated functions of the parameter $\tilde{k}$.

Now it can be noted that if it is possible to make the replacement $\tilde{k}⟺q$, then another equation will be obtained

$${\tilde{X}}^{\left(q\right)}=N_1^{\left(q\right)}{J}_2^{\left(q\right)}\left(y\right)+N_2^{\left(q\right)}{J}_0^{\left(q\right)}\left(y\right)+N_3^{\left(q\right)}$$

(204)

and in this way, the equations Eq. (200) and Eq. (204) will give the following solution for the zero-order elliptic integral in the Legendre form

$$J_0^{\left(4\right)}\left(y\right)\equiv J_0^{\left(q\right)}\left(y\right)={\displaystyle \frac{2{\tilde{X}}^{\left(q\right)}}{\sqrt[3]{2}{N}_2^{\left(q\right)}}}-{\displaystyle \frac{N_3^{\left(q\right)}}{N_2^{\left(q\right)}}}$$

$$-{\displaystyle \frac{\sqrt[3]{2}{\tilde{X}}^{\left(\tilde{k}\right)}\left(\sqrt[3]{2}-2N_2^{\left(q\right)}\right)}{N_2^{\left(q\right)}\left(\sqrt[3]{2}-2\right)}}-{\displaystyle \frac{N_1^{\left(q\right)}\left(2{\tilde{X}}^{\left(q\right)}-\sqrt[3]{2}{N}_3^{\left(q\right)}\right)}{N_2^{\left(q\right)}\sqrt[3]{2}}}.$$

(205)

As a consequence of the replacement $\tilde{k}⟺q$, one also has the following equality, which follows from the defining relations (198)

$${\displaystyle \frac{{\tilde{k}}^4-{\tilde{k}}^2+1}{3g_2\left(\tilde{k}\right)}}={\displaystyle \frac{q^4-q^2+1}{3g_2\left(q\right)}}.$$

(206)

If this is a quite a general condition (not imposing any restrictions), then the replacement $\tilde{k}⟺q$ will be justified. Evidently this will depend on the determination of $g_2^{\left(q\right)}\equiv g_2\left(q\right)$, which will be given in the next sections.

10.5. Another Integral in the Weierstrass Form After Applying the Four-dimensional Uniformization

Further by "four-dimensional uniformization" it shall be meant that the polynomial of the fourth degree under the square root in the denominator of the elliptic integral shall be "uniformized" by introducing complicated functions, depending on the Weierstrass function and its derivative, which represent functions of the complex variable z. The term "uniformization" shall be clarified in the next section.

10.5.1. The Standard Method for $3D$ Uniformization in Elliptic Functions Theory - Reminder of the Basic Facts

Let in accord with the standard notions in elliptic functions theory [44] and [46], the s.c. "fundamental parallelogram" is defined on the two-periodic lattice (see also the monograph [47])

$$\Lambda \equiv \left\{m{\omega }_1+n{\omega }_2;0\le m\le 1,0\le n\le 1,Im\left({\displaystyle \frac{{\omega }_1}{{\omega }_2}}\right)>0\right\}.$$

(207)

The two periods ${\omega }_1$ and ${\omega }_2$ constitute the basis on the complex plane and thus each complex number can be represented as $z=$$m{\omega }_1+n{\omega }_2$. Then the Weierstrass function

$$\rho \left(z\right)\equiv {\displaystyle \frac{1}{z^2}}+\sum \left({\displaystyle \frac{1}{{(z-\omega )}^2}}-{\displaystyle \frac{1}{{\omega }^2}}\right)$$

(208)

and its first derivative

$${\rho }^{\prime }\left(z\right)\equiv -{\displaystyle \frac{2}{z^3}}-2\sum {\displaystyle \frac{1}{{(z-\omega )}^3}}$$

(209)

satisfy the parametrizable form of the cubic algebraic equation

$$y^2=4x^3-g_{2}x-g_3\mathrm{i}.\mathrm{e}.x\equiv \rho \left(z\right),y\equiv {\rho }^{\prime }\left(z\right),$$

(210)

where the functions $\rho \left(z\right)$ and ${\rho }^{\prime }\left(z\right)$ are called also uniformization functions for the above cubic algebraic curve $y^2=4x^3-g_{2}x-g_3$. The summation in the defining formulaes (208) for the Weierstrass function $\rho \left(z\right)$ is over all points on the two-dimensional lattice on the complex plane, i.e. on all possible numbers m and n. The coefficient functions $g_2$ and $g_3$ are the s.c. "invariants of the Weierstrass function" $\rho \left(z\right)$

$$g_2=60G_4=60\sum {\displaystyle \frac{\prime 1}{{\omega }^4}},g_3=140G_6=140\sum {\displaystyle \frac{\prime 1}{{\omega }^6}},$$

(211)

defined as infinite convergent sums over the period $\omega$ of the two-periodic lattice. The prime "′" above the two sums means that the period $\omega =0$ is excluded from the summation, so that the expressions would not tend to infinity.

10.6. Application of the Weierstrass Integral and of the Weierstrass Elliptic Curve in the Parametrizable Form

This section has the purpose to clarify the importance of finding the Weierstrass invariants $g_2$ and $g_3$ for the calculation of an elliptic integral in the Weierstrass form $\int {\displaystyle \frac{dx}{\sqrt{4x^3-g_{2}x-g_3}}}$. In the next sections on the base of the methods for four-dimensional uniformization, a primary goal will be to find the Weierstrass invariants $g_2$ and $g_3$ not only for integrals in the Weierstrass form, but also for integrals in the Legendre form and also for more general integrals (98) (for $n=0)$

$$J_0^{\left(4\right)}\left(y\right)=\int {\displaystyle \frac{dy}{\sqrt{a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4}}},$$

with an arbitrary polynomial of the fourth-degree under the square in the denominator. Here it shall be clarified that finding the Weierstrass invariants will be possible after transforming the initially given integral in the Legendre form into an integral in the Weierstrass form, and in such a way the Weierstrass invariants $g_2$ and $g_3$ will depend in a complicated manner on the modulus q of the integral in the Legendre form. Thus, as it will be demonstrated by a theorem in the monograph by [48], for known $g_2$ and $g_3$ the solution (and the roots) of the cubic equation $y^2=4x^3-g_{2}x-g_3$ will be possible to be determined. Then, the periods ${\omega }_1$ and ${\omega }_2$ of the two-dimensional lattice of the fundamental parallelogram (207) can be calculated.

Let us first note that finding a solution $z-z_0$ of the integral $\int {\displaystyle \frac{dx}{\sqrt{4x^3-g_{2}x-g_3}}}$ in the Weierstrass form is equivalent to knowing the roots of the cubic polynomial $y^2=4x^3-g_{2}x-g_3$. This "equivalency" can be written as [71]

$$z-z_0=\underset{x_0}{\overset{x}{\int }}{\displaystyle \frac{dx}{\sqrt{4x^3-g_{2}x-g_3}}}\iff y^2=4x^3-g_{2}x-g_3.$$

(212)

In fact, the equality $x=\rho \left(z\right)$ is a solution also of the elliptic integral in (212) and this is the known problem of "inversion" for elliptic integrals. In (210) $\rho \left(z\right)$ is the Weierstrass function [71]

$$\rho \left(z\right)\equiv {\displaystyle \frac{1}{z^2}}+\sum \left({\displaystyle \frac{1}{{(z-\omega )}^2}}-{\displaystyle \frac{1}{{\omega }^2}}\right)={\displaystyle \frac{1}{z^2}}+{\displaystyle \frac{g_2}{20}}{z}^2+{\displaystyle \frac{g_3}{28}}{z}^4+.....$$

(213)

From the second representation for the Weierstrass function $\rho \left(z\right)$ it becomes clear why in the literature it is denoted by $\rho (z;g_2,g_3)$.The following theorem, proved in Ch. 6 of the monograph by Knapp [48] allows us to understand how the result of integration for the elliptic integral (173) ${\tilde{J}}_0^{\left(q\right)}\left(y\right)=\int {\displaystyle \frac{dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}$ can be represented in an analytical form by means of the Weierstrass integral, after calculating the Weierstrass invariants $g_2$ and $g_3$ by means of the method of four-dimensional uniformization. This method will be presented in the next sections.

10.7. A Theorem about the Unique Correspondence Between the Periods ${\omega }_1$ and ${\omega }_2$ on the Two-Dimensional Complex Lattice and the Weierstrass Invariants $(g_2,g_3)$ for Elliptic Integrals in the Weierstrass Form

Theorem 8.

[48] There exists an unique correspondence $(g_2\left(\Lambda \right),g_3\left(\Lambda \right))\iff \Lambda$ between the lattices Λ on the complex plane $\mathbb{C}$ and the pair $(g_2,g_3)$ of the complex numbers, the Weierstrass invariants such that the discriminant $\Delta \equiv g_2^3-27g_3^2$ of the cubic polynomial $4x^3-g_{2}x-g_3=0$ is different from zero. If $e_1,e_2,e_3$ are the roots of the above polynomial, i. e.

$$4x^3-g_{2}x-g_3=4(z-e_1)(z-e_2)(z-e_3),$$

(214)

then the periods $({\omega }_1,{\omega }_2)$ on the complex lattice Λ (207) can be found from the following integrals

$${\omega }_1=\underset{{\Gamma }_1}{\int }{\displaystyle \frac{dz}{2\sqrt{(z-e_1)(z-e_2)(z-e_3)}}},$$

(215)

$${\omega }_2=\underset{\Gamma 2}{\int }{\displaystyle \frac{dz}{2\sqrt{(z-e_1)(z-e_2)(z-e_3)}}}.$$

(216)

In the first integral the unique branch ${\Gamma }_1$ of the square root is chosen with cuts from $e_1$ to $e_2$ and from $e_3$ to ∞ and in the second integral the branch ${\Gamma }_2$ is with cuts from $e_2$ to $e_3$ and from $e_1$ to ∞.

Moreover, for every lattice $\Lambda$ on the complex plane, defined according to (207), a biholomorphic mapping is determined as $\phi :\mathbb{C}/\Lambda \mapsto E\left(\mathbb{C}\right)$, where $E\left(\mathbb{C}\right)$ denotes the elliptic curve $y^2=4x^3-g_2\left(\Lambda \right)x-g_3\left(\Lambda \right)$. The mapping $\phi$ is defined by means of the Weierstrass elliptic function $\rho \left(z\right)$ and its derivative ${\rho }^{\prime }\left(z\right)$ and the inverse mapping - by the corresponding elliptic integral.

Since $\rho \left(z\right)$ and ${\rho }^{\prime }\left(z\right)$ are uniformization functions, they will satisfy the cubic polynomial [44]

$${\left({\rho }^{\prime }\left(z\right)\right)}^2=4\left(\rho \left(z\right)-e_1\right)\left(\rho \left(z\right)-e_2\right)\left(\rho \left(z\right)-e_3\right),$$

(217)

from where the following relations are easily derived [44], following from the algebraic properties of the roots of a cubic equation

$$e_1+e_2+e_3=0,e_1{e}_2+e_2{e}_3+e_1{e}_3=-{\displaystyle \frac{g_2}{4}},e_1{e}_2{e}_3={\displaystyle \frac{g_3}{4}}.$$

(218)

Now it may easily be calculated that

$$g_2^3-27g_3^2=16{(e_1-e_2)}^2{(e_2-e_3)}^2{(e_1-e_3)}^2.$$

(219)

If the discriminant $\Delta \equiv g_2^3-27g_3^2$ (a well known notion from higher algebra) is different from zero, then the cubic equation $4x^3-g_2\left(\Lambda \right)x-g_3\left(\Lambda \right)=0$ will not have coinciding roots, i.e.

$$e_1\ne e_2\ne e_3$$

(220)

This property is very important also for the theory of elliptic curves (see also the monograph [49]) since if the property $g_2^3-27g_3^2\ne 0$ is fulfilled then a lattice $\Lambda$ on the complex plane exists, so that the Weierstrass invariants $g_2$ and $g_3$ can be defined as (211)

$$g_2=60G_4=60\sum {\displaystyle \frac{\prime 1}{{\omega }^4}},g_3=140G_6=140\sum {\displaystyle \frac{\prime 1}{{\omega }^6}}.$$

At this point it becomes clear why the correspondence $(g_2\left(\Lambda \right),g_3\left(\Lambda \right))\iff \Lambda$ in the cited monograph [48] is in both directions - previously we proved that if the lattice $\Lambda$ is given, then the Weierstrass invariants $g_2$ and $g_3$ can be defined, i.e. $\Lambda \Rightarrow (g_2\left(\Lambda \right),g_3\left(\Lambda \right))$.

10.7.1. The Consistency Problem for a Non-Zero Discriminant for the Case, When the Invariants in the Weierstrass Integrals are Expressed Through the Modulus Parameter q of the Integral in the Legendre Form

Now the statement is proved in the opposite directions - if the Weierstrass invariants are given, then a complex lattice can be defined $(g_2\left(\Lambda \right),g_3\left(\Lambda \right))\Rightarrow \Lambda$. Note that this statement is not trivial and it can be reformulated in the following way: if $g_2$ and $g_3$ are given (real or complex) numbers or even functions of some parameter (such examples will be shown further), then a two-dimensional lattice $\Lambda$ on the complex plane with periods ${\omega }_1$ and ${\omega }_2$ can be defined, so that the sums in the defining equalities for $g_2$ and $g_3$ in (211) will be convergent. This is also a problem, which needs further investigation.

Further it shall be understood why this is important in view of the fact that an elliptic integral in the Legendre form (with a modulus parameter q) shall be transformed to an elliptic integral in the Weierstrass form $\int {\displaystyle \frac{dx}{\sqrt{4x^3-g_{2}x-g_3}}}$. The Weierstrass integral is characterized by invariants $g_2$ and $g_3$, which shall be expressed by complicated polynomials, depending on the modulus parameter q, i.e. $g_2=g_2\left(q\right)$ and $g_3=g_3\left(q\right)$. So a natural question arizes: will the integral in the Weierstrass form with the calculated values for the Weierstrass invariants $g_2$ and $g_3$ be correctly defined? The answer will be affirmative, because from the above correspondence it will follow that the lattice $\Lambda$ can be defined. However, if $g_2$ and $g_3$ are defined as $g_2=g_2\left(q\right)$ and $g_3=g_3\left(q\right)$, then the condition $\Delta \equiv g_2^3-27g_3^2\ne 0$ for the discriminant should also be fulfilled, which in fact will mean that the higher-order $g_2^3-27g_3^2$ polynomial of q should have no zeroes! Further in the section for four-dimensional uniformization, we shall obtain the parametrization for $g_2$ and $g_3$ (254) in the form of the following higher-order polynomials of the modulus parameter q of the elliptic integral (173) in the Legendre form

$${\overline{g}}_2^{(1,2)}\equiv q^2+{\displaystyle \frac{1}{12}}{(1+q^2)}^2>0$$

and (255) (the notation in both equalities will be ${\overline{g}}_2$ and ${\overline{g}}_3$)

$${\overline{g}}_3^{(1,2)}\equiv {\displaystyle \frac{q^6}{4}}-{\displaystyle \frac{q^4}{6}}-{\displaystyle \frac{5q^2}{12}}-{\displaystyle \frac{1}{6^3}}{\left[5q^2-1\right]}^3.$$

Then will the higher-order polynomial ${\overline{g}}_2^{3(1,2)}-27{\overline{g}}_3^{2(1,2)}$ have the property to be without any roots? The same question appears with respect to the found second pair of Weierstrass invariants (268) (the same as (254)

$${\overline{g}}_2=q^2+{\displaystyle \frac{1}{12}}{(1+q^2)}^2>0$$

and also (269)

$${\overline{g}}_3=-{\displaystyle \frac{1}{6}}{q}^2(1+q^2)+{\displaystyle \frac{{(1+q^2)}^3}{6^3}}.$$

In view of the theorems, which have been proved it follow that in both cases the higher-order polynomial ${\overline{g}}_2^{3(1,2)}-27{\overline{g}}_3^{2(1,2)}$ should have no roots! Since both the polynomials ${\overline{g}}_2^{3(1,2)}-27{\overline{g}}_3^{2(1,2)}$ depend on higher degrees of q, which is defined for $0\le q<1$ it is known from higher algebra that there are many theorems [53,56], investigating the case when an arbitrary polynomial of the $n-$th order will have (or will not have) roots in the circle $q<1$. One such theorem, formulated in the monograph by Obreshkoff [53], is the theorem by Schur [55], published yet in 1918. The concrete application of this theorem requires cumbersome analytic calculations, because it should be applied to a chain of algebraic equations of diminishing degrees. For example, for the case of the polynomials (268) and (269), the polynomial ${\overline{g}}_2^{3(1,2)}-27{\overline{g}}_3^{2(1,2)}$ will be of the $12-$th degree! For polynomials of the fourth degree, the Schur theorem has been applied several times in the papers [18] and [19] and the results have shown excellent consistency with the physical content, implied in these polynomials. Moreover, in [18] it has been demonstrated that the Schur theorem can be used for the both cases of proving that a certain polynomial has roots in the circle $q<1$ or does not have such roots. In the next sections, we shall give only the formulation of this theorem.

On the other hand, if in one of the cases for the polynomial ${\overline{g}}_2^{3(1,2)}-27{\overline{g}}_3^{2(1,2)}$ it is proved that there are roots, then this will mean that either the corresponding parametrization should be rejected, or that the elliptic integral in the Legendre form is not defined for the values of q, satisfying the above polynomial ${\overline{g}}_2^{3(1,2)}-27{\overline{g}}_3^{2(1,2)}=0$. In standard monographs on elliptic functions and integrals, this interesting problem is not investigated.

Now it remains to prove that the three roots $e_1,e_2,e_3$ are different and they will correspond to three different values for.the Weierstrass function $\rho \left(z_1\right)$, $\rho \left(z_2\right)$, $\rho (z_1+z_2)$, corresponding to the values

$$z_1={\displaystyle \frac{{\omega }_1}{2}},z_2={\displaystyle \frac{{\omega }_2}{2}},z_1+z_2={\displaystyle \frac{{\omega }_1}{2}}+{\displaystyle \frac{{\omega }_2}{2}}.$$

(221)

The proof is standard but since it contains some nontrivial facts, related to functions on the complex plane, it shall be given below.

A detailed study whether the periods $({\omega }_1,{\omega }_2)$ are real and imaginary depending on the roots of the polynomial $4x^3-g_{2}x-g_3=0$ and on the properties of the integrals (215) and (216) is given in Ch. $6.5$ of the known monograph [69].

10.8. Definition for an Elliptic Function. A Proof that the Three Roots of the Cubic Equation $4x^3-g_{2}x-g_3=0$ are Different

10.8.1. Defining Properties of the Weierstrass Function as an Elliptic Function

This section does not contain any new material, but the properties of the Weierstrass function are important in view of the definition of a new problem for uniformization of higher order (higher than four) algebraic equations by means of complicated functions, depending on the Weierstrass function and its derivative.

Let us give first the definition for an elliptic function $f\left(z\right)$. By definition, $f\left(z\right)$ is an elliptic function if it is meromorphic and doubly periodic. The last means that the equality

$$f(z+n{\omega }_1+m{\omega }_2)=f\left(z\right)$$

(222)

is fulfilled, where n and m are integer numbers.

The Weierstrass function (208), defined by the equality

$$\rho \left(z\right)\equiv {\displaystyle \frac{1}{z^2}}+\sum \left({\displaystyle \frac{1}{{(z-\omega )}^2}}-{\displaystyle \frac{1}{{\omega }^2}}\right)$$

is an example for an elliptic function. This definition is related to the definition of the notion of elliptic integral in the sense that $x=\rho \left(z\right)$ is a solution of the integral in the Weierstrass form (212) $z-z_0=\underset{x_0}{\overset{x}{\int }}{\displaystyle \frac{dx}{\sqrt{4x^3-g_{2}x-g_3}}}$. However, the two definitions for an elliptic integral and for an elliptic function are different.

First, the notion of a meromorphic function should be clarified: a function is meromorphic, if it does not have singular points other than poles [44]. Singular points in complex analyses and for the function $f\left(z\right)$ of a complex variable can be of three different kinds [71]:

1. Removable singular points (one of them for example can be denoted as a), for which ${lim}_{z\to a}f\left(z\right)$ exists and is finite. This means also that $\mid f\left(z\right)\mid \le M$ (M is a constant), consequently the Laurent expansion around the point a is

$$f\left(z\right)=c_0+c_1(z-a)+.......+c_n{(z-a)}^n\Rightarrow \underset{z\to a}{lim}f\left(z\right)=c_0$$

(223)

and does not contain a main part, i.e. these are terms with poles in the Laurent decomposition. The constant Mbecause of the Cauchy inequalities can be chosen as $\mid c_n\mid \le M{\rho }^{-n}$ and since $\rho$ can be chosen to be very small, this would mean that all coefficients $c_n$ with negative indices will be zero.

2. Essentially singular points - these are points for which ${lim}_{z\to a}f\left(z\right)$ does not exist. For this case the Sohotsky theorem is valid, according to which if a is an essentially singular point, then for every complex number A there exists a sequence of points $z_k\to a$, for which ${lim}_{z\to a}f\left(z\right)=A$. The difference of this case from the previous one is that in the previous case the limit point was only one, while here there is an infinite amount of limit points. Note also that the infinite point $A\to \infty$ can also be a limit point for the sequence, since if it was only finite, then one can take the function $g\left(z\right)={\displaystyle \frac{1}{f\left(z\right)}}$, for which ${lim}_{z\to a}g\left(z\right)=0$, consequently $z_k\to a$ will be a removable singular point (but with many limit points for the sequence). Also, a point is essentially singular if and only if the main part in the Laurent decomposition contains infinitely many terms.

Consequently, the two cases cases can be united by writing the Laurent decomposition in the form

$$f\left(z\right)=f_1\left(z\right)+f_2\left(z\right)=\sum _{n=-\infty }^{\infty }{c}_n{(z-a)}^n,$$

(224)

where

$$c_n={\displaystyle \frac{1}{2\pi i}}\underset{\gamma }{\int }{\displaystyle \frac{f\left(\zeta \right)d\zeta }{{(\zeta -a)}^{n+1}}}(n=0,\pm 1,\pm 2.......\pm n,......)$$

(225)

and the integration is performed over the circle $\gamma :\mid z-a\mid =\rho$ . For n with plus signs and for $n\to \infty$ there will be infinitely many positive terms in the decomposition (223) - this is the first case of removable singular points. For n with a negative sign there will be infinitely many terms in the main part of the decomposition - the second case of essentially singular points.

3. The third case is the case with poles, for which ${lim}_{z\to a}f\left(z\right)=\infty$, which means that the point a is a pole if and only if the main part (these are the pole terms) in the Laurent decomposition contains a finite number of pole terms [71]:

$$f\left(z\right)={\displaystyle \frac{c_{-n}}{{(z-a)}^n}}+........+{\displaystyle \frac{c_{-1}}{z-a}}+\sum _{k=0}^{\infty }{c}_k{(z-a)}^k.$$

(226)

Proof that the Weierstrass Function is Meromorphic

Since the case with poles is characteristic for the Weierstrass function, we have to prove that the pole terms are finite. Since the fundamental parallelogram covers the entire complex plane and for the pole terms we have ${lim}_{z\to a}f\left(z\right)=\infty$, we have to establish that this property is valid for arbitrary poles, even for those periods ${\omega }_1$ and ${\omega }_2$, tending to infinity. Let us take the term inside the brackets in the defining equality (208) for the Weierstrass function and calculate the limit when $\omega \to \infty$ [44]

$${\displaystyle \frac{1}{{(z-\omega )}^2}}-{\displaystyle \frac{1}{{\omega }^2}}={\displaystyle \frac{2z\omega -z^2}{{\omega }^2{(z{\omega }^{-1}-1)}^2}}={\displaystyle \frac{\omega }{{\omega }^2}}.{\displaystyle \frac{2z-z^2{\omega }^{-1}}{{\omega }^2{(z{\omega }^{-1}-1)}^2}}.$$

(227)

When the number $\mid \omega \mid$ tends to infinity, the term

$${\displaystyle \frac{2z-z^2{\omega }^{-1}}{{\omega }^2{(z{\omega }^{-1}-1)}^2}}\approx 2z.$$

(228)

This means that when $z\to \infty$, the term ${\displaystyle \frac{1}{{(z-\omega )}^2}}-{\displaystyle \frac{1}{{\omega }^2}}$ is infinite, i.e. this is a pole term. Now it remains to prove that although every term ${\displaystyle \frac{1}{{(z-\omega )}^2}}-{\displaystyle \frac{1}{{\omega }^2}}$ is infinite, the sum of all these terms is finite, corresponding to the finite number of pole terms in the decomposition (226). From (227) and (228) and when $\mid \omega \mid \to \infty$ it can be obtained that [44]

$$\mid {\displaystyle \frac{1}{{(z-\omega )}^2}}-{\displaystyle \frac{1}{{\omega }^2}}\mid <{\displaystyle \frac{C}{\mid \omega {\mid }^3}}.$$

(229)

It can be checked that the sum $\sum {\displaystyle \frac{\prime 1}{\mid \omega {\mid }^3}}$ is convergent (see also [47]). It can be represented in the form

$$\sum {\displaystyle \frac{\prime 1}{\mid \omega {\mid }^3}}=\sum _{n=1}^{\infty }\sum _{max(\mid p_1\mid ,\mid p_2\mid )=n}\mid p_1{\omega }_1+p_2{\omega }_2{\mid }^{-3}\le \sum _{n=1}^{\infty }8n{\left(nh\right)}^{-3},$$

(230)

where h is the smaller height of the fundamental parallelogram. Thus, since the sum $\sum _{n=1}^{\infty }{\displaystyle \frac{1}{n^2}}$ is convergent, the sum $\sum \left({\displaystyle \frac{1}{{(z-\omega )}^2}}-{\displaystyle \frac{1}{{\omega }^2}}\right)$ of the infinite number of pole terms in the defining equality (208) for $\rho \left(z\right)$ is finite.

A Proof that the Weierstrass Function is Doubly-Periodic  For the proof of the double-periodicity of $\rho \left(z\right)$, one has to take its derivative (209), written as ${\rho }^{\prime }\left(z\right)\equiv -2\sum {\displaystyle \frac{1}{{(z-\omega )}^3}}$, where the summation includes also the zero-points of the lattice $\Lambda$. Since the derivative is also a doubly-periodic function, the difference

$${\rho }^{\prime }(z+{\omega }_i)-{\rho }^{\prime }\left(z\right)=0⟹\rho (z+{\omega }_i)-\rho \left(z\right)=C(C=const)$$

(231)

will be zero - a translation of z with a period ${\omega }_i$ of the lattice will not change ${\rho }^{\prime }\left(z\right)$ because

$${\displaystyle \frac{1}{z^3}}-{\displaystyle \frac{1}{{(z-{\omega }_i)}^3}}={\displaystyle \frac{-3z^2{\omega }_i+3z{\omega }_i^2+{\omega }_i^3}{{(z-{\omega }_i)}^3}}$$

(232)

and in the sum $\sum \left({\displaystyle \frac{1}{z^3}}-{\displaystyle \frac{1}{{(z-{\omega }_i)}^3}}\right)$ all the terms in the nominator $-3z^2{\omega }_i+3z{\omega }_i^2+{\omega }_i^3$ will cancel. The first term $-3z^2{\omega }_i$ and the third term ${\omega }_i^3$ will cancel because ${\omega }_i$ can be both positive and negative, the second term $3z{\omega }_i^2$ will become zero, because $z=m{\omega }_1+n{\omega }_2$ also can be both positive and negative. (${\omega }_1,{\omega }_2$ can be positive and negative, but the numbers m and n according to the definition (207) are $0\le m\le 1$ , $0\le n\le 1$). Setting up $z\equiv -{\displaystyle \frac{{\omega }_i}{2}}$ and taking into account that $\rho \left(z\right)$ is an even function, it is easily obtained [49]

$$\rho \left({\displaystyle \frac{{\omega }_i}{2}}\right)-\rho (-{\displaystyle \frac{{\omega }_i}{2}})=C=0.$$

(233)

The function $\rho \left(z\right)$ has at the points of the lattice double poles, no other singular points. According to a theorem from complex analysis, the sum of the poles and of the zeroes inside a closed contour (in the case this is the fundamental parallelogram) should be comparable to zero. This means that inside the fundamental parallelogram there are two zeroes u and v with a sum, comparable to zero by the module of the lattice, i.e. $u+v=0$$\left(mod\Lambda \right)$[44]. For every constant C the poles of the function $\rho \left(z\right)-C$ coincide with the poles of $\rho \left(z\right)$. Consequently, there are only two points $\rho \left(u\right)=$$\rho \left(v\right)=C$, for which $u+v=0$$\left(mod\Lambda \right)$. At the points, for which $u\equiv -u$$\left(mod\Lambda \right)$ and these two points coincide, the function $\rho \left(u\right)-C$ takes twice the value $z=u$ and the derivative ${\rho }^{\prime }\left(z\right)$ is zero. Then inside the fundamental parallelogram there are four points, for which the equality $u\equiv -u$$\left(mod\Lambda \right)$ is fulfilled - the pole $z_0=0$ of the function $\rho \left(z\right)$, the other two points $z_1={\displaystyle \frac{{\omega }_1}{2}}$ and $z_2={\displaystyle \frac{{\omega }_2}{2}}$ are zeroes of the derivative ${\rho }^{\prime }\left(z\right)$. Because of the equality (217) [44]

$${\left({\rho }^{\prime }\left(z\right)\right)}^2=4\left(\rho \left(z\right)-e_1\right)\left(\rho \left(z\right)-e_2\right)\left(\rho \left(z\right)-e_3\right),$$

this means that the roots $e_1,e_2,e_3$ are

$$e_1=\rho \left({\displaystyle \frac{{\omega }_1}{2}}\right)=\rho \left(z_1\right),e_2=\rho \left({\displaystyle \frac{{\omega }_2}{2}}\right)=\rho \left(z_2\right),e_3=\rho \left({\displaystyle \frac{{\omega }_1+{\omega }_2}{2}}\right)=\rho \left(z_3\right).$$

(234)

A Differential Equation for the Weierstrass Function as a Consequence of the Group-Theoretical Law for Summing up Points on the Cubic Curve

Since the group law for summing up points on the cubic curve and the resulting differential equation will be applied further for finding another formulaes for the Weierstrass invariants $g_2$ and $g_3$ in (270), we shall review briefly the derivation of this equation, given in [44].

The third root $e_3$ in (234) follows from the property (218) $e_1+e_2+e_3=0$$\left(mod\Lambda \right)$ of the cubic equation andthe fact that the Weierstrass function is an even function, which means that $e_3=\rho (-{\displaystyle \frac{{\omega }_1+{\omega }_2}{2}})=\rho \left({\displaystyle \frac{{\omega }_1+{\omega }_2}{2}}\right)$. The property (218) is not only an algebraic one, but also a consequence from the following well-known theorem in complex analysis, according to which the sum of the residues for the singular points (situated inside the fundamental parallelogram) of an elliptic function is equal to zero. This can be proved by taking the residues for the function $g\left(z\right)$

$$\sum _{\Pi }resg\left(z\right)={\displaystyle \frac{1}{2\pi i}}\underset{\partial \Pi }{\int }g\left(z\right)dz,$$

(235)

where $\Pi$ denotes the fundamental parallelogram, $\partial \Pi$ is the boundary contour and the function $g\left(z\right)$ is taken in the form $g\left(z\right)={\displaystyle \frac{f^{\prime }\left(z\right)}{f\left(z\right)}}$. Also, if $a_i$ are zeroes and poles of an elliptic function and $r_i$ are their orders (positive for the zeroes and negative for the poles), then it can easily be proved that

$$\sum r_i=0,\sum r_i{a}_i=0.$$

(236)

The three roots $e_1,e_2,e_3$ mean that if a straight line $y=fx+h$ ($f,h$ are coefficients) intersects the cubic curve $y^2=4x^3-g_{2}x-g_3$, then from the cubic equation

$${(fx+h)}^2=4x^3-g_{2}x-g_3,$$

(237)

also from the Wiets formulae for the roots and from the equalities (234) for the three values of the Weierstrass function at the points $z=z_1$, $z=z_2$ and $z=z_3$, it follows [44]

$$\rho \left(z_1\right)+\rho \left(z_2\right)+\rho (z_1+z_2)={\displaystyle \frac{f^2}{4}}.$$

(238)

Since $x=\rho \left(z\right)$ and $y={\rho }^{\prime }\left(z\right)$ are uniformization functions for the cubic curve and also for the straight line $y=fx+h$, let us write it in terms of the Weierstrass function

$${\rho }^{\prime }\left(z_1\right)=f\rho \left(z_1\right)+h,{\rho }^{\prime }\left(z_2\right)=f\rho \left(z_2\right)+h.$$

(239)

Expressing from here the coefficient f as

$$f={\displaystyle \frac{{\rho }^{\prime }\left(z_1\right)-{\rho }^{\prime }\left(z_2\right)}{\rho \left(z_1\right)-\rho \left(z_2\right)}}$$

(240)

and substituting into the cubic equation (237), the following differential equation is obtained, which is a consequence from the group-theoretic law for summing up points on a cubic curve

$$\rho (z_1+z_2)=\rho \left(z_1\right)+\rho \left(z_2\right)+{\displaystyle \frac{1}{4}}{\left({\displaystyle \frac{{\rho }^{\prime }\left(z_1\right)-{\rho }^{\prime }\left(z_2\right)}{\rho \left(z_1\right)-\rho \left(z_2\right)}}\right)}^2.$$

(241)

Closely related with the group theoretical law on elliptic curves is the problem about rational points on elliptic curves. Both topics are extensively discussed in the monographs [44,50] and [51].

10.9. Uniformization of Fourth-Degree Algebraic Equations and Application to the Theory of Elliptic Integrals. Another Representation of an Integral in the Legendre Form as an Integral in the Weierstrass Form

The derivation of the second representation of the elliptic integral in the Weierstrass form is based on the following theorem, proved in the monograph by Whittaker and Watson [60].

Theorem 9.

Let

$$J_0^{\left(4\right)}\left(y\right)=\int {\displaystyle \frac{dy}{\sqrt{a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4}}}$$

(242)

is a zero-order elliptic integral with the fourth-degree polynomial under the square root in the denominator

$$Y\equiv f\left(y\right)\equiv a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4,f\left(y_0\right)=0.$$

(243)

Then after the variable changes

$${\displaystyle \frac{1}{t-y_0}}={\displaystyle \frac{\sigma -\frac{1}{2}{A}_2}{A_3}},{\displaystyle \frac{1}{y-y_0}}={\displaystyle \frac{s-\frac{1}{2}{A}_2}{A_3}}$$

(244)

the integral can be rewritten as

$$I_0^{\left(3\right)}\left(x\right)={\int }_s^{\infty }{\displaystyle \frac{d\sigma }{\sqrt{4{\sigma }^3-{\overline{g}}_2\sigma -{\overline{g}}_3}}},$$

(245)

where the Weierstrass invariants ${\overline{g}}_2$ and ${\overline{g}}_3$ can be exactly calculated as

$${\overline{g}}_2=3A_2^2-4A_1{A}_3;{\overline{g}}_3=2A_1{A}_2{A}_3-A_2^3-A_0{A}_3^2$$

(246)

and $A_0,A_1,A_2,A_3$ are the introduced notations for the polynomial expressions

$$A_0\equiv a_0,A_1\equiv a_0{y}_0+a_1,$$

(247)

$$A_2\equiv a_0{y}_0^2+2a_1{y}_0+a_2,A_3\equiv a_0{y}_0^3+3a_1{y}_0^2+3a_2{y}_0+a_3.$$

(248)

It can easily be verified that if the expressions (247) and (248) for $A_0,A_1,A_2,A_3$ are used, then the invariant ${\overline{g}}_2$ in (246) can be calculated as

$${\overline{g}}_2=a_0{a}_4-4a_1{a}_3+3a_2^2.$$

(249)

Later in the other formulation of the theorem in the next paragraph this formulae will be used, but it will be taken from the theorem in the monograph [60]. Note that under the changes in the coefficients

$$4a_1\to a_1,6a_2\to a_2,4a_3\to a_3,a_4\to a_4$$

(250)

in the fourth-degree polynomial (243) the invariant ${\overline{g}}_2$ changes as

$${\overline{g}}_2=a_0{a}_4-{\displaystyle \frac{a_1{a}_3}{4}}+{\displaystyle \frac{a_2^2}{12}}.$$

(251)

Technically, the calculation of the invariant ${\overline{g}}_3$ is more difficult.

10.9.1. Four-Dimensional Uniformization of the (Zero-Order) Elliptic Integral in the Legendre Form and its Representation as an Integral in the Weierstrass Form

First we shall clarify that the application of the above formulaes (244) - (248) to the elliptic integral (242) with a polynomial of the fourth degree in the square root and transforming it to the elliptic integral $I_0^{\left(3\right)}\left(x\right)$ (245) will be called a "four-dimensional uniformization". If the fourth-degree polynomial (243)

$$Y\equiv f\left(y\right)\equiv a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4,f\left(y_0\right)=0.$$

is represented in a parametrizable form by functions y and Y, expressed in a complicated way by functions, depending on the Weierstrass function $\rho \left(z\right)$ and its derivative ${\rho }^{\prime }\left(z\right)$, then this will be called "an uniformization of the fourth-degree polynomial (243) with the uniformization functions $y=y(y_0,z,\rho \left(z\right),{\rho }^{\prime }\left(z\right))$ and $Y=Y(y_0,z,\rho \left(z\right),{\rho }^{\prime }\left(z\right))$. Further these uniformization functions will be shown.

In the concrete case for the polynomial in the Legendre form

$$f\left(y\right)\equiv (1-y^2)(1-q^2{y}^2)=1-(1+q^2)y^2+y^4,$$

(252)

one can compute ${\overline{g}}_2^{(1,2)}$ and ${\overline{g}}_3^{(1,2)}$ for the roots of the polynomial $f\left(y\right)$ (243) $y_{1,2}=\pm 1$ and $y_{3,4}=\pm {\displaystyle \frac{1}{q}}$ and also for the coefficients $a_0,a_1,a_2,a_3,a_4$, taking into account the identity of the polynomials (243) and (252)

$$a_0\equiv a_1\equiv a_3\equiv 0,a_2\equiv -{\displaystyle \frac{1}{6}}(1+q^2),a_4\equiv 1.$$

(253)

The two pairs of coefficient functions ${\overline{g}}_2^{(1,2)}$ and ${\overline{g}}_3^{(1,2)}$, calculated according to the formulaes (246) - (248) and depending on functions of the modulus parameter q are

$${\overline{g}}_2^{(1,2)}\equiv q^2+{\displaystyle \frac{1}{12}}{(1+q^2)}^2>0,$$

(254)

$${\overline{g}}_3^{(1,2)}\equiv {\displaystyle \frac{q^6}{4}}-{\displaystyle \frac{q^4}{6}}-{\displaystyle \frac{5q^2}{12}}-{\displaystyle \frac{1}{6^3}}{\left[5q^2-1\right]}^3$$

(255)

$$=\left({\displaystyle \frac{1}{4}}-{\displaystyle \frac{5^3}{6^3}}\right)q^6+\left(-{\displaystyle \frac{1}{6}}+{\displaystyle \frac{75}{6^3}}\right)q^4$$

$$-\left({\displaystyle \frac{1}{12}}+{\displaystyle \frac{15}{6^3}}\right)q^2+{\displaystyle \frac{1}{6^3}}.$$

(256)

Similar expressions for ${\overline{g}}_2^{(3,4)}$ and ${\overline{g}}_3^{(3,4)}$ can also be written, but they will not be given here. Here the important fact, related to the theory of elliptic functions is that if ${\overline{g}}_2^{(1,2)}$, ${\overline{g}}_3^{(1,2)}$ or ${\overline{g}}_2^{(3,4)}$, ${\overline{g}}_3^{(3,4)}$ are known, then from the theorem in the monograph [48] it will follow that the Weierstrass function $\rho (z;{\overline{g}}_2,{\overline{g}}_3)$ is uniquely determined.

10.9.2. Another Formulation of the Theorem for the Uniformization of Algebraic Equations of the Fourth-Order and Application to the Theory of Elliptic Integrals

Theorem 10.

If the Weierstrass elliptic function is defined as $\rho (z;{\overline{g}}_2,{\overline{g}}_3)$, then the two variables

$$y=y_0+{\displaystyle \frac{1}{4}}{\displaystyle \frac{f^{´}\left(y_0\right)}{\left[\rho (z;{\overline{g}}_2,{\overline{g}}_3)-\frac{1}{24}{f}^{´}\left(y_0\right)\right]}},$$

(257)

$$Y=-{\displaystyle \frac{1}{4}}{\displaystyle \frac{f^{´}\left(y_0\right){\rho }^{´}(z;{\overline{g}}_2,{\overline{g}}_3)}{{\left[\rho (z;{\overline{g}}_2,{\overline{g}}_3)-\frac{1}{24}{f}^{´}\left(y_0\right)\right]}^2}}$$

(258)

define an unique parametrization of the $4D$ algebraic curve (243)

$$Y\equiv f\left(y\right)\equiv a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4,f\left(y_0\right)=0.$$

Consequently, the method for "four-dimensional uniformization" can be considered to be the analogue of the known "three-dimensional uniformization", given by (210).

In the above formulaes the point $y_0$ is a root of the quartic polynomial (243), but this requirement can be relaxed and $y_0=a$ can be an arbitrary constant. The polynomial $f\left(y\right)$ is supposed not to have multiple roots. Then the uniformization function y can be represented in the form, given in the monograph [60] (the formulaes are shown as a part of two problems, after the formulation of the theorem for the two uniformization functions (257) and (258)

$$y=y_0+{\displaystyle \frac{F_1(y_0,z)}{F_2(y_0,z)}},$$

(259)

where $F_1(y_0,z)$ and $F_2(y_0,z)$ are the expressions

$$F_1(y_0,z)\equiv {\left[f\left(y_0\right)\right]}^{{\displaystyle \frac{1}{2}}}{\rho }^{\prime }(z;{\overline{g}}_2,{\overline{g}}_3)+{\displaystyle \frac{1}{2}}{f}^{´}\left(y_0\right)\{\rho (z;{\overline{g}}_2,{\overline{g}}_3)$$

$$-{\displaystyle \frac{1}{24}}{f}^{´}\left(y_0\right)\}+{\displaystyle \frac{1}{24}}f\left(y_0\right){\displaystyle \frac{d^{3}f\left(y_{.}\right)}{d{y}^3}}{\mid }_{y=y_0},$$

(260)

$$F_2(y_0,z)\equiv 2{\left[\rho (z;{\overline{g}}_2,{\overline{g}}_3)-{\displaystyle \frac{1}{24}}{f}^{´}\left(y_0\right)\right]}^2-{\displaystyle \frac{1}{48}}f\left(y_0\right)f^{IV}\left(y_0\right).$$

(261)

If again $y_0$ is supposed to be a root of the algebraic equation $f\left(y\right)=0$ and the Weierstrass function $\rho (z;{\overline{g}}_2,{\overline{g}}_3)$ is known, then every other value of the polynomial $f\left(y\right)$ can be found by means of the following concise formulae ([60], also formulated as a separate problem)

$${\left[f\left(y\right)\right]}^{{\displaystyle \frac{1}{2}}}=-{\displaystyle \frac{f^{´}\left(y_0\right){\rho }^{´}(z;{\overline{g}}_2,{\overline{g}}_3)}{4\left[\rho (z;{\overline{g}}_2,{\overline{g}}_3)-\frac{1}{24}{f}^{´}\left(y_0\right)\right]}}.$$

(262)

The formulaes (260), (261) and (262) will not be used further in this paper, but probably the interested reader can find some application of these expressions.

10.9.3. Second Formulation of the Theorem for Four-Dimensional Uniformization

This theorem is also exposed in the monograph by E.T. Whittaker and G. N.Watson [60]. We shall review the corresponding theorem by making a slight modification in the fourth-degree polynomial (243)- the numbers $4,6$ and 4 in the coefficients are removed.

Theorem 11.

If $f\left(s\right)$ is a fourth-order polynomial

$$f\left(y\right)=a_0{y}^4+a_1{y}^3+a_2{y}^2+a_{3}y+a_4$$

(263)

and the variable transformation $z=\underset{y_0}{\overset{y}{\int }}{\displaystyle \frac{ds}{\sqrt{f\left(s\right)}}}$ is performed in the $4D$ elliptic integral

$$J_0^{\left(4\right)}\left(y\right)=\int {\displaystyle \frac{dy}{\sqrt{a_0{y}^4+a_1{y}^3+a_2{y}^2+a_{3}y+a_4}}},$$

(264)

then the Weierstrass invariants${\overline{g}}_2$and${\overline{g}}_3$in the integral (245)$I_0^{\left(3\right)}\left(s\right)={\int }_s^{\infty }{\displaystyle \frac{d\sigma }{\sqrt{4{\sigma }^3-{\overline{g}}_2\sigma -{\overline{g}}_3}}}$are uniquely calculated as

$${\overline{g}}_2=a_0{a}_4-{\displaystyle \frac{1}{4}}{a}_1{a}_3+{\displaystyle \frac{1}{12}}{a}_2^2,$$

(265)

$${\overline{g}}_3={\displaystyle \frac{1}{6}}{a}_0{a}_2{a}_4+{\displaystyle \frac{1}{48}}{a}_1{a}_2{a}_3-{\displaystyle \frac{a_2^3}{6^3}}-{\displaystyle \frac{1}{16}}(a_0{a}_3+a_1^2{a}_4).$$

(266)

One may note that in this formulation of the theorem, the Weierstrass invariants ${\overline{g}}_2$ and ${\overline{g}}_3$ do not depend on any root $y_0$ of the polynomial (263).

10.9.4. Four-Dimensional Uniformization of the Elliptic Integral in the Legendre Form and its Representation in the Weierstrass Form - Another Calculation by Using the Second Formulation of the Theorem

Let us take as an example, which will be used further, the polynomial

$$f\left(y\right)\equiv (1-y^2)(1-q^2{y}^2)=1-(1+q^2)y^2+y^4.$$

(267)

The two Weierstrass invariants are calculated to be

$${\overline{g}}_2=q^2+{\displaystyle \frac{1}{12}}{(1+q^2)}^2>0,$$

(268)

$${\overline{g}}_3=-{\displaystyle \frac{1}{6}}{q}^2(1+q^2)+{\displaystyle \frac{{(1+q^2)}^3}{6^3}}.$$

(269)

Clearly, the result for ${\overline{g}}_2$ in inequality (268) coincides with the result for ${\overline{g}}_2^{(1,2)}\equiv q^2+{\displaystyle \frac{1}{12}}{(1+q^2)}^2>0$ in (254), but formulae (269) for ${\overline{g}}_3$ is quite different from formulaes (255) and (256) for ${\overline{g}}_3^{(1,2)}$. The same refers also to the formulae for ${\overline{g}}_3^{(3,4)}$, not presented here.

Important consequences: 1. For this case there is only one pair of Weierstrass invariants ${\overline{g}}_2$ and ${\overline{g}}_3$ instead of the two pairs ${\overline{g}}_2^{(1,2)},{\overline{g}}_3^{(1,2)}$ and ${\overline{g}}_2^{(3,4)},{\overline{g}}_3^{(3,4)}$ for the previous case, considered in the preceding Section 10.9.

2. The expression for ${\overline{g}}_2$ (268) coincides with the expression for ${\overline{g}}_2^{(1,2)}$ in (254), but expression ${\overline{g}}_3$ in (269) is different from expression (255) for ${\overline{g}}_2^{(1,2)}$. It is not clear whether this will be so for the case of an elliptic integral with an arbitrary polynomial in the denominator.

3. The two invariants ${\overline{g}}_2$ and ${\overline{g}}_3$ are symmetric with respect to the interchanges $a_0\rightleftarrows a_4$ and also $a_1\rightleftarrows a_3$. The polynomial $f\left(y\right)$ is also invariant with respect to these changes.

10.10. Comparison of Elliptic Integrals in Different Variables and with Different Weierstrass Invariants

In this section the method for calculating the Weierstrass invariants $g_2$ and $g_3$ will be presented, which is based on the combination of two approaches. In the first approach (see section "Transforming an Elliptic Integral in the Legendre Form into an Elliptic Integral in the Weierstrass Form") the initial integral ${\tilde{J}}_0^{\left(q\right)}\left(y\right)$ was brought to the Weierstrass form (174) (multiplied by the conformal factor $\sqrt{a}=\sqrt{{\displaystyle \frac{9g_3}{g_2}}.{\displaystyle \frac{\left(q^4-q^2+1\right)}{\left(2q^4-5q^2+2\right)}}}$) by means of the transformation (175) $x={\displaystyle \frac{a}{y^2}}+b$, as a result of which the integral acquires the form (174). We shall call this representation "equation A".

The second approach is based on the method of "four-dimensional uniformization", as a result of which the integral is brought to the second representation (245) $I_0^{\left(3\right)}\left(x\right)={\int }_s^{\infty }{\displaystyle \frac{d\sigma }{\sqrt{4{\sigma }^3-{\overline{g}}_2\sigma -{\overline{g}}_3}}}$ - "equation B". As a result of the two representations, the following two equalities can be written

$${\tilde{J}}_0^{\left(q\right)}\left(y\right)=\int {\displaystyle \frac{dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}=-\sqrt{a}\int {\displaystyle \frac{dx}{\sqrt{4x^3-g_{2}x-g_3}}}$$

(270)

$$=\int {\displaystyle \frac{d\sigma }{\sqrt{4{\sigma }^3-{\overline{g}}_2\sigma -{\overline{g}}_3}}}.$$

(271)

Now we shall use the following theorem, proved in the monograph by Hancock [61] in order to transform the integral (270).

Theorem 12.

The elliptic integral in the Weierstrass representation (174) (equation (A) in (270) in terms of the "x"-variable) can be written in terms of the variable σ as

$$\int {\displaystyle \frac{dx}{\sqrt{4x^3-g_{2}x-g_3}}}=-\int {\displaystyle \frac{d\sigma }{\sqrt{4{\sigma }^3-g_2\sigma -g_3}}}$$

(272)

after performing the variable transformation

$$x=x_0+{\displaystyle \frac{(x_1-x_0)(x_2-x_0)}{\sigma -x_0}},$$

(273)

where $x_0,x_1$ and $x_2$ are the roots of the cubic polynomial $4x^3-g_{2}x-g_3=0$.

In the original formulation of the theorem in [61] (page 29) the transformation $t=e_1+{\displaystyle \frac{(e_2-e_1)(e_3-e_1)}{s-e_1}}$ has been applied to the integral $\int {\displaystyle \frac{dt}{\sqrt{(t-e_1)}(t-e_2)(t-e_3)}}$. Note however the minus sign in front of the right term in Eq. (272), which is important and can easily be confirmed by a direct calculation. Applying again the above theorem with respect to equation (B) in (271) and also Eq. (245), the equality (272) can be written as

$$\int {\displaystyle \frac{d\sigma }{\sqrt{4{\sigma }^3-g_2\sigma -g_3}}}=-{\displaystyle \frac{1}{\sqrt[3]{a}}}\int {\displaystyle \frac{d\overline{\sigma }}{\sqrt{4{\overline{\sigma }}^3-{\widehat{g}}_2\overline{\sigma }-{\widehat{g}}_3}}},$$

(274)

where the variable transformation $\sigma ={\displaystyle \frac{1}{\sqrt[3]{a}}}\overline{\sigma }$ has been performed and the notations

$${\widehat{g}}_2\equiv {\overline{g}}_2{a}^{{\displaystyle \frac{2}{3}}},{\widehat{g}}_3\equiv {\overline{g}}_{3}a$$

(275)

are introduced. From equation (A) (272) and equation (B) (274) (combined with (245), we obtain

$$\int {\displaystyle \frac{d\sigma }{\sqrt{4{\sigma }^3-g_2\sigma -g_3}}}={\displaystyle \frac{1}{\sqrt[3]{a}}}\int {\displaystyle \frac{d\sigma }{\sqrt{4{\sigma }^3-{\widehat{g}}_2\sigma -{\widehat{g}}_3}}}.$$

(276)

Note the important peculiarity of eq. (276): the Weierstrass invariants in the left-hand side of (276) are $g_2$, $g_3$, while in the right-hand side they are ${\widehat{g}}_2\equiv {\overline{g}}_2{a}^{{\displaystyle \frac{2}{3}}},{\widehat{g}}_3\equiv {\overline{g}}_{3}a$.

10.10.1. Algebraic Roots of the Cubic Equations and Comparison of the Weierstrass Invariants - the Case of Positive Discriminant

From the Cardano formulae, we have for the roots of the cubic polynomial $4{\sigma }^3-g_2\sigma -g_3=0$ ($i=1,2,3$) the following formulae (see the monographs [53] and [54])

$${\alpha }_1=A_1+B_1=\sqrt[3]{{\displaystyle \frac{g_3}{8}}+{\displaystyle \frac{1}{4}}\sqrt{{\displaystyle \frac{g_3^2}{4}}-{\displaystyle \frac{1}{27.4}}{g}_2^3}}$$

$$+\sqrt[3]{{\displaystyle \frac{g_3}{8}}-{\displaystyle \frac{1}{4}}\sqrt{{\displaystyle \frac{g_3^2}{4}}-{\displaystyle \frac{1}{27.4}}{g}_2^3}}.$$

(277)

This is only the first root of the cubic equation $4{\sigma }^3-g_2\sigma -g_3=0$. Further we shall investigate the case of a positive discriminant $\Delta ={\displaystyle \frac{q^2}{4}}+{\displaystyle \frac{p^3}{27}}>0$ for the cubic equation $x^3+px+q=0$[53]. For the given cubic equation $4{\sigma }^3-g_2\sigma -g_3=0$ with the first root (277), the discriminant is

$$\Delta ={\displaystyle \frac{g_3^2}{4}}-{\displaystyle \frac{1}{27.4}}{g}_2^3={\displaystyle \frac{g_3^2}{4}}-{\displaystyle \frac{1}{108}}{g}_2^3>0\Rightarrow g_3^2>{\displaystyle \frac{1}{27}}{g}_2^3.$$

(278)

Let us assume that the roots of this polynomial are ${\alpha }_1,{\beta }_1,{\gamma }_1$ and the roots of the polynomial $4{\sigma }^3-{\widehat{g}}_2\sigma -\widehat{g}=0$ are ${\alpha }_2,{\beta }_2,{\gamma }_2$. Then the other two roots ${\beta }_1,{\gamma }_1$ can be found if the formulae for the s.c. $n-$th root of the number 1 is taken into account This simple formulae is

$$\sqrt[n]{1}=\sqrt[n]{ϵ}=cos{\displaystyle \frac{2\pi }{n}}+isin{\displaystyle \frac{2\pi }{n}}.$$

(279)

For $n=3$ and by means of the formulae

$$\sqrt[3]{1}=\sqrt[3]{ϵ}=cos{\displaystyle \frac{2\pi }{3}}+isin{\displaystyle \frac{2\pi }{3}}=-cos60+isin60=-{\displaystyle \frac{1}{2}}+{\displaystyle \frac{\sqrt[.]{3}}{2}},$$

(280)

the roots ${\beta }_1,{\gamma }_1$ can be expressed as

$${\beta }_1=A_1ϵ+B_1{ϵ}^2=A_1(-{\displaystyle \frac{1}{2}}+i{\displaystyle \frac{\sqrt[.]{3}}{2}})$$

$$+B_1(-{\displaystyle \frac{1}{2}}-i{\displaystyle \frac{\sqrt[.]{3}}{2}})=-{\displaystyle \frac{(A_1+B_1)}{2}}+i{\displaystyle \frac{\sqrt[.]{3}}{2}}(A_1-B_1),$$

(281)

$${\gamma }_1=A_1{ϵ}^2+B_1{ϵ}^{.}=A_1(-{\displaystyle \frac{1}{2}}-i{\displaystyle \frac{\sqrt[.]{3}}{2}})$$

$$+B_1(-{\displaystyle \frac{1}{2}}+i{\displaystyle \frac{\sqrt[.]{3}}{2}})=-{\displaystyle \frac{(A_1+B_1)}{2}}-i{\displaystyle \frac{\sqrt[.]{3}}{2}}(A_1-B_1).$$

(282)

From the denominators of the two integrals in (276) we should have the equality of two cubic polynomials under a square root, the first expression for the polynomial obtained from the square of the second cubic polynomial, multiplied by $\sqrt[3]{a}$

$$\sqrt[.]{4{\sigma }^3-g_2\sigma -g_3}=\sqrt[3]{a}\sqrt[.]{4{\sigma }^3-{\widehat{g}}_2\sigma -{\widehat{g}}_3}.$$

(283)

This expression can be written in the form

$$(\sigma -{\alpha }_1)(\sigma -{\beta }_1)(\sigma -{\gamma }_1)={\left(\sqrt[3]{a}\right)}^2(\sigma -{\alpha }_2)(\sigma -{\beta }_2)(\sigma -{\gamma }_2)$$

$$=(a^{{\displaystyle \frac{2}{9}}}\sigma -a^{{\displaystyle \frac{2}{9}}}{\alpha }_2)(a^{{\displaystyle \frac{2}{9}}}\sigma -a^{{\displaystyle \frac{2}{9}}}{\beta }_2)(a^{{\displaystyle \frac{2}{9}}}\sigma -a^{{\displaystyle \frac{2}{9}}}{\gamma }_2).$$

(284)

The first root $a^{{\displaystyle \frac{2}{9}}}{\alpha }_2$ of the cubic equation in the right-hand side of (284), analogously to the Cardano formulae (277) for the previous equation, can be expressed as

$$a^{{\displaystyle \frac{2}{9}}}{\alpha }_2=A_2+B_2=\sqrt[3]{{\displaystyle \frac{{\widehat{g}}_3}{8}}+{\displaystyle \frac{1}{4}}\sqrt{{\displaystyle \frac{{\widehat{g}}_3^2}{4}}-{\displaystyle \frac{1}{27.4}}{\widehat{g}}_2^3}}$$

$$+\sqrt[3]{{\displaystyle \frac{{\widehat{g}}_3}{8}}-{\displaystyle \frac{1}{4}}\sqrt{{\displaystyle \frac{{\widehat{g}}_3^2}{4}}-{\displaystyle \frac{1}{27.4}}{\widehat{g}}_2^3}}.$$

(285)

From here, after multiplication of the right-hand side of (285) with $a^{-{\displaystyle \frac{2}{9}}}$, only the root ${\alpha }_2$ can be expressed as

$${\alpha }_2=a^{-{\displaystyle \frac{2}{9}}}(A_2+B_2)=\sqrt[3]{a^{-{\displaystyle \frac{2}{27}}}{\displaystyle \frac{{\widehat{g}}_3}{8}}+{\displaystyle \frac{1}{4}}\sqrt{a^{-{\displaystyle \frac{1}{27}}}{\displaystyle \frac{{\widehat{g}}_3^2}{4}}-{\displaystyle \frac{1}{27.4}}{a}^{-{\displaystyle \frac{1}{27}}}{\widehat{g}}_2^3}}$$

$$+\sqrt[3]{a^{-{\displaystyle \frac{2}{27}}}{\displaystyle \frac{{\widehat{g}}_3}{8}}-{\displaystyle \frac{1}{4}}\sqrt{a^{-{\displaystyle \frac{1}{27}}}{\displaystyle \frac{{\widehat{g}}_3^2}{4}}-{\displaystyle \frac{1}{27.4}}{a}^{-{\displaystyle \frac{1}{27}}}{\widehat{g}}_2^3}}.$$

(286)

Now it is important to note that ${\alpha }_2$ is not another root of the cubic equation (284)

$$0=(\sigma -{\alpha }_1)(\sigma -{\beta }_1)(\sigma -{\gamma }_1)={\left(\sqrt[3]{a}\right)}^2(\sigma -{\alpha }_2)(\sigma -{\beta }_2)(\sigma -{\gamma }_2).$$

If one assumes that ${\alpha }_2$ is another root (different from the roots ${\alpha }_1,{\beta }_1,{\gamma }_1$), then this would mean that the polynomial in the left-hand side of (284) should be of the fourth degree, which evidently will be a contradiction. Consequently, the root ${\alpha }_2$ should coincide with one of the roots ${\alpha }_1,{\beta }_1,{\gamma }_1$. Let us first assume that ${\alpha }_2={\alpha }_1$. From the comparison of the corresponding expressions under the cubic and square roots of (277) and (286) it follows that the equality is fulfilled when the following relations between the Weierstrass invariants $g_2,g_3$ and ${\widehat{g}}_2,{\widehat{g}}_3$ are fulfilled

$$g_3=a^{-{\displaystyle \frac{1}{27}}}{\widehat{g}}_3,g_2^3=a^{-{\displaystyle \frac{1}{27}}}{\widehat{g}}_2^3.$$

(287)

The last relation can be rewritten, expressing the dependance of ${\widehat{g}}_2,{\widehat{g}}_3$ on $g_2,g_3$

$${\widehat{g}}_3=a^{{\displaystyle \frac{1}{27}}}{g}_3,{\widehat{g}}_2=a^{{\displaystyle \frac{1}{81}}}{g}_2.$$

(288)

It can be noted that these relations are different from the ones, obtained in a previous paper [21]

$$g_2^3=a^2{\widehat{g}}_2^3,g_3=-a{\widehat{g}}_3.$$

(289)

In fact, the correct relations are (287) and (288) (and not (289), which are wrong), because they are obtained from the relation for the roots (284). This relation takes into account the multiplication both of the roots and the variables in (284) with the conformal factor $a^{{\displaystyle \frac{2}{9}}}$ and not with $({\sigma }_0,{\sigma }_1,{\sigma }_2)=\sqrt[3]{a}({\sigma }_0,{\sigma }_1,{\sigma }_2)$, as this is applied in the paper [21].

10.10.2. Finding the Relation Between the Other Two Roots $({\beta }_2,{\gamma }_2)$ and $({\beta }_1,{\gamma }_1)$ for the Case When ${\alpha }_2={\alpha }_1-$ the Case of Positive Discriminant

The other two roots (${\beta }_2,{\gamma }_2)$ of the second equation (283) $4{\sigma }^3-{\widehat{g}}_2\sigma -{\widehat{g}}_3=0$ can be written similarly to the expressions (281) and (282) for ${\beta }_1$ and ${\gamma }_1$, but with indices $"2"$

$${\beta }_2=A_2ϵ+B_2{ϵ}^2=-{\displaystyle \frac{(A_2+B_2)}{2}}+i{\displaystyle \frac{\sqrt[.]{3}}{2}}(A_2-B_2),$$

(290)

$${\gamma }_2=A_2{ϵ}^2+B_2{ϵ}^{.}.=-{\displaystyle \frac{(A_2+B_2)}{2}}-i{\displaystyle \frac{\sqrt[.]{3}}{2}}(A_2-B_2).$$

(291)

Taking into account that ${\alpha }_1={\alpha }_2$ and also the equality (286), it can be obtained

$$A_2+B_2=a^{{\displaystyle \frac{2}{9}}}(A_1+B_1)\Rightarrow A_2-B_2=a^{{\displaystyle \frac{2}{9}}}(A_1-B_1).$$

(292)

If substituted in the formulaes (290) and (291) for ${\beta }_2$ and ${\gamma }_2$, the following relations are found between the roots (${\beta }_2,{\gamma }_2)$ and (${\beta }_1,{\gamma }_1)$

$${\beta }_2=a^{{\displaystyle \frac{2}{9}}}{\beta }_1,{\gamma }_2=a^{{\displaystyle \frac{2}{9}}}{\gamma }_1,$$

(293)

where $({\beta }_1,{\gamma }_1)$ are the complex roots, given by the formulaes (281) and (282).

10.10.3. Expressing the Ratio of the Weierstrass Invariants $g_2$ and $g_3$ as a Rational Function of Higher-Power Polynomials, Depending on the Modulus Parameter q of the Elliptic Integral in the Legendre Form- the Case of Positive Discriminant

Finding the relation between the Weierstrass invariants $g_2,g_3$ and ${\widehat{g}}_2,{\widehat{g}}_3$ in (287) and (288) is rather insufficient and does not give any valuable information because of the following two reasons:

A. The variable a is not just a parameter, but represents a complicated function, depending on the ratio ${\displaystyle \frac{g_3}{g_2}}$ and on the rational function $K\left(q\right)$ and is determined in (176) as $a={\displaystyle \frac{9g_3}{g_2}}K\left(q\right)$, where (178) $K\left(q\right)\equiv {\displaystyle \frac{q^4-q^2+1}{2q^4-5q^2+2}}$ .

B. The invariants ${\widehat{g}}_2,{\widehat{g}}_3$ do not contain any useful information - they appeared after the variable change (275)

$${\widehat{g}}_2\equiv {\overline{g}}_2{a}^{{\displaystyle \frac{2}{3}}},{\widehat{g}}_3\equiv {\overline{g}}_{3}a,$$

when the equality of the two integrals (276)

$$\int {\displaystyle \frac{d\sigma }{\sqrt{4{\sigma }^3-g_2\sigma -g_3}}}={\displaystyle \frac{1}{\sqrt[3]{a}}}\int {\displaystyle \frac{d\sigma }{\sqrt{4{\sigma }^3-{\widehat{g}}_2\sigma -{\widehat{g}}_3}}}.$$

was obtained. Although we found the relation between the roots of the polynomials in both sides of (276), it is more important to find the dependance of the invariants $g_2,g_3$ on the invariants ${\overline{g}}_2,{\overline{g}}_3$, expressed through polynomials of the modulus q and calculated in the process of the four-dimensional uniformization, presented in the previous sections.

That is why, now we shall combine the transformations (287) and (275). We obtain

$$g_3=a^{-{\displaystyle \frac{1}{27}}}{\widehat{g}}_3=a^{-{\displaystyle \frac{1}{27}}}{\overline{g}}_{3}a=a^{{\displaystyle \frac{26}{27}}}{\overline{g}}_3,$$

(294)

$$g_2={\left(a^{-{\displaystyle \frac{1}{27}}}\right)}^{{\displaystyle \frac{1}{3}}}{\widehat{g}}_2=a^{-{\displaystyle \frac{1}{81}}}{\overline{g}}_2{a}^{{\displaystyle \frac{2}{3}}}={\overline{g}}_2{a}^{{\displaystyle \frac{53}{81}}}.$$

(295)

Taking into account (176) for $a={\displaystyle \frac{9g_3}{g_2}}K\left(q\right)$ with (178) for $K\left(q\right)\equiv {\displaystyle \frac{q^4-q^2+1}{2q^4-5q^2+2}}$ and substituting in the above equalities, it follows that

$$g_3={\overline{g}}_3{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{26}{27}}}.{\left({\displaystyle \frac{g_3}{g_2}}\right)}^{{\displaystyle \frac{26}{27}}},$$

(296)

$$g_2={\overline{g}}_2{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{53}{81}}}.{\left({\displaystyle \frac{g_3}{g_2}}\right)}^{{\displaystyle \frac{53}{81}}}.$$

(297)

Dividing the above two equalities, it can be obtained

$${\displaystyle \frac{g_3}{g_2}}={\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}{\left(9K\left(q\right)\right)}^{\left[{\displaystyle \frac{26.3}{27.3}}-{\displaystyle \frac{53}{81}}\right]}{\left({\displaystyle \frac{g_3}{g_2}}\right)}^{\left[{\displaystyle \frac{26}{27}}-{\displaystyle \frac{53}{81}}\right]}.$$

(298)

After some transformations the ratio ${\displaystyle \frac{g_3}{g_2}}$ can be expressed as a function of powers of the ratio ${\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}$ (which is a rational function of polynomials of the modulus parameter q- see (173) and (255)) and powers of $\left(9K\right(q\left)\right)$, which is also a rational function of higher-order polynomials in the nominator and in the denominator according to (178)

$${\displaystyle \frac{g_3}{g_2}}={\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{56}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{56}}.}.$$

(299)

This is an important formulae, enabling to express $g_3$ as a linear function of $g_2$ and a complicated function $N\left(q\right)$

$$g_3=g_{2}N\left(q\right),N\left(q\right):={\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{56}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{56}}.}.$$

(300)

In such a way, substituting the above expression in formulae (277) for the (real-valued) first root ${\alpha }_1=A_1+B_1$ of the cubic polynomial $4{\sigma }^3-g_2\sigma -g_3=0$, one can obtain

$${\alpha }_1=A_1+B_1=\sqrt[3]{g_2}[\sqrt[3]{F\left(q\right)+{\displaystyle \frac{1}{4}}\sqrt[.]{\tilde{\Delta }}}$$

$$+\sqrt[3]{F\left(q\right)-{\displaystyle \frac{1}{4}}\sqrt[.]{\tilde{\Delta }}}]=\sqrt[3]{g_2}.\left[{\tilde{A}}_1+{\tilde{B}}_1\right],$$

(301)

where $F\left(q\right)$ and $\tilde{\Delta }$ (obtained from the discriminant formulae (278) $\Delta ={\displaystyle \frac{g_3^2}{4}}-{\displaystyle \frac{1}{27.4}}{g}_2^3>0$) denote the following functions

$$F\left(q\right):={\displaystyle \frac{1}{8}}{\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{56}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{56}}.}={\displaystyle \frac{1}{8}}N\left(q\right),$$

(302)

$$\tilde{\Delta }:={\displaystyle \frac{1}{4}}{\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{28}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{28}}.}-{\displaystyle \frac{1}{4.27}}{g}_2.$$

(303)

In (301) the following new notations were introduced

$${\tilde{A}}_1:=\sqrt[3]{F\left(q\right)+{\displaystyle \frac{1}{4}}\sqrt[.]{\tilde{\Delta }}},{\tilde{B}}_1:=\sqrt[3]{F\left(q\right)-{\displaystyle \frac{1}{4}}\sqrt[.]{\tilde{\Delta }}}.$$

(304)

Since the "modified" discriminant $\tilde{\Delta }$ (303) is linearly dependant on the Weierstrass invariant $g_2$, it is evident that $g_2$ can be found if the ratio ${\displaystyle \frac{{\tilde{A}}_1}{{\tilde{B}}_1}}$ is known. This will be used in the next sections.

10.10.4. A Simple Proof for the Positivity of the Weierstrass invariants $g_2$ for the One Pair of Weierstrass Invariants ${\overline{g}}_2$ (268) and ${\overline{g}}_3$ (269)

From (303) the condition for a non-zero discriminant can be written as

$$g_2\ne 27{\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{28}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{28}}.}.$$

(305)

Now we shall prove that this expression is different from zero and also is positive. The nominator of the expression for (178) $K\left(q\right)\equiv {\displaystyle \frac{q^4-q^2+1}{2q^4-5q^2+2}}$ and also the expression for ${\overline{g}}_3$ (269)

$${\overline{g}}_3=-{\displaystyle \frac{1}{6}}{q}^2(1+q^2)+{\displaystyle \frac{{(1+q^2)}^3}{6^3}}.$$

are different from zero. For the nominator of $K\left(q\right)$ the positivity is easily established, because it can be represented in the form

$$q^4-q^2+1=q^4-2q^2+1={(q^2-1)}^2+q^2>0,$$

(306)

which evidently is greater than zero.

Consequently, whether $g_2$ will be positive depends on the ratio ${\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}$. The invariant ${\overline{g}}_2$ (268)

$${\overline{g}}_2=q^2+{\displaystyle \frac{1}{12}}{(1+q^2)}^2>0$$

evidently is positive.

Let us check the sign of expression (269) for ${\overline{g}}_3$, which can be rewritten as

$${\overline{g}}_3={\displaystyle \frac{(q^2+1)}{6^3}}(q^4-2{\displaystyle \frac{3.17}{2}}{q}^2+1)$$

(307)

$$={\displaystyle \frac{(q^2+1)}{6^3}}\left[{\left(q^2-{\displaystyle \frac{3.17}{2}}\right)}^2+1-{\left({\displaystyle \frac{3.17}{2}}\right)}^2\right].$$

(308)

This expression will be positive when the inequalities

$$q^2>{\displaystyle \frac{3.17}{2}}+{\displaystyle \frac{\sqrt[.]{6.{17}^2-4}}{2}}>0,$$

(309)

and also

$$q^2>{\displaystyle \frac{3.17}{2}}-{\displaystyle \frac{\sqrt[.]{6.{17}^2-4}}{2}}>0$$

(310)

are fulfilled. The first case (309) should be disregarded, because it means that $q>6.804$, which is impossible because q should be in the range $0<q<1$.

The second expression (310) is positive when $q>2.168$, which is also impossible. Consequently, the expression in the square brackets in (308) cannot be positive when $q>6.804$ and $q>2.168$.

However, when $q<6.804$ and $q<2.168$ (always fulfilled since $0<q<1$), the expression is positive, because

$$\left(q^2-{\displaystyle \frac{3.17}{2}}-{\displaystyle \frac{\sqrt[.]{6.{17}^2-4}}{2}}\right)\left(q^2-{\displaystyle \frac{3.17}{2}}+{\displaystyle \frac{\sqrt[.]{6.{17}^2-4}}{2}}\right)>0$$

(311)

both expressions in the round brackets are negative and (311) is with a plus sign.

Therefore, we proved that for all values of q the invariant $g_2$ (305) is positive and from the assumption for the positivity of the discriminant (278)

$$\Delta ={\displaystyle \frac{g_3^2}{4}}-{\displaystyle \frac{1}{108}}{g}_2^3>0\Rightarrow g_3^2>{\displaystyle \frac{1}{27}}{g}_2^3$$

it follows also that $\mid g_3\mid >{\displaystyle \frac{1}{3\sqrt{3}}}\mid g_2^{{\displaystyle \frac{3}{2}}}\mid$.

10.10.5. The Possibility for a positive and Real Invariant $g_2$, Complex Invariant $g_3$ and a Negative Discriminant- A Theorem in [49]

Formally, it may seem that if $g_3$ is negative but the condition for the moduli is fulfilled, then the discriminant $\Delta$ will be positive. In fact, here we may have the case, when $g_2$ is positive, but $g_3$ can be even complex. This problem is not investigated in the mathematical literature, but it is important that the possibility for $g_3$ to be complex does not contradict a theorem in problem 12, paragraph 6 of chapter 1 in the monograph by Koblitz [49]:

Theorem 13.

The roots $e_1,e_2,e_3$ of the polynomial $4x^3-g_{2}x-g_3=0$ are all real one if and only if the invariants $g_2$ and $g_3$ are all real ones and the discriminant $\Delta =g_2^3-27g_3^2$$>0$ is positive.

In this paper, we are not investigating the case when all the roots of the polynomial $4x^3-g_{2}x-g_3=0$ are real, we deal with the case when only the first root ${\alpha }_1$ is a real one and the other two roots ${\beta }_1$ and ${\gamma }_1$ are imaginary. Consequently, it is not obligatory that the invariants $g_2$ and $g_3$ should be real ones. We obtained that the invariant $g_2$ is positive, but let us remember that for the invariant $g_3$ we found the formulae (300)

$$g_3=g_{2}N\left(q\right),N\left(q\right):={\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{56}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{56}}.}.$$

We already proved that the ratio ${\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}$ is positive and real, also we proved that the nominator of the function (178) $K\left(q\right)\equiv {\displaystyle \frac{q^4-q^2+1}{2q^4-5q^2+2}}$ is positive. The denominator can be decomposed as

$$2q^4-5q^2+2=2\left[\left(q^2-2\right)\left(q^2-{\displaystyle \frac{1}{2}}\right)\right].$$

(312)

Now if one assumes that this expression is positive,then this is possible only when the following inequalities are simultaneously fulfilled

$$q<\sqrt{2}\approx 1.41,q<{\displaystyle \frac{1}{\sqrt{2}}}\approx 0.707106.$$

(313)

The other two inequalities with the positive sign should be discarded because $0<q<1$. But if the inequalities (313) are fulfilled and $K\left(q\right)$, $N\left(q\right)$ and $g_3$ are positive, then from the real-valuedness of $g_2$ and $g_2$ and also the positivity and real-valuedness of the discriminant $\Delta =g_2^3-27g_3^2$$>0$ from the theorem in [49] it should follow that the roots $e_1,e_2,e_3$ of the polynomial $4x^3-g_{2}x-g_3=0$ should all be real. But this will be in contradiction to the case, investigated in this paper of one real root and two imaginary roots.

Consequently, an important conclusion can be made that there can be two cases:

1-st case. Positive invariant $g_2$, imaginary $g_3$ (or both imaginary invariants $g_2$ and $g_3$) and positive discriminant $\Delta =g_2^3-27g_3^2$$>0$- then not all the roots should be real ones (the presently investigated case).

2-nd case. Real roots ${\alpha }_1,{\beta }_1,{\gamma }_1$ (the case not investigated in this paper).Then the invariants $g_2$ and $g_3$ should be real and the discriminant $\Delta$ should be negative. From the point of view of practical applications for the inter-satellite communications, it is evident that the case of a negative discriminant and modulus parameter $q<{\displaystyle \frac{1}{\sqrt{2}}}\approx 0.707106$ (remember that for a GPS orbit $e\approx 0.01$ and $q={\displaystyle \frac{1-e}{1+e}}$ ), such a case seems to be realistic and should be investigated.

Note that in the strict mathematical sense, the theorem in the monograph by Koblitz [49] should be reformulated and proved.

The case of a negative discriminant for a general cubic equation $x^3+px+\tilde{q}=0$ is investigated in details in the monograph by Obreshkoff [53] and the roots have been expressed by the concise formulae

$$x_{0,1,2}=2\sqrt{-{\displaystyle \frac{p}{3}}}cos{\displaystyle \frac{\phi +2k\pi }{3}},k=0,1,2$$

(314)

and the angular variable $\phi$ is given by the formulae

$$tg\phi =-{\displaystyle \frac{2\sqrt{-\Delta }}{\tilde{q}}}.$$

(315)

We have not investigated the case of the other pair of invariants (254) ${\overline{g}}_2^{(1,2)}$ (identical with the invariant (268)) and (256) ${\overline{g}}_3^{(1,2)}$, when the second invariant is given by a polynomial of the sixth degree and depending on its roots, it can be either positive or negative.

10.11. Finding the Weierstrass Function $\rho \left(z\right)$ and the Weierstrass invariant $g_2$ From the Known Roots of the Cubic Equation and the Group-Theoretical Law for Summing up of Points on the Cubic Curve

The main goal of every investigation, related to the cubic algebraic equation in its parametrizable form is to find the expression for the Weierstrass function. $\rho \left(z\right)$ by known Weierstrass invariants $g_2$ and $g_3$. Since we started from a Weierstrass integral, derived from an elliptic integral in the Legendre form, previously we found how the invariant $g_3$ can be expressed through the other invariant $g_2$ (300)

$$g_3=g_{2}N\left(q\right),N\left(q\right):={\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{56}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{56}}.}.$$

Now it remains to find the other Weierstrass invariant $g_2$ and the Weierstrass function itself. This shall be performed by using the found roots ${\alpha }_1,{\beta }_1,{\gamma }_1$ for the case of a positive discriminant (278) $\Delta ={\displaystyle \frac{g_3^2}{4}}-{\displaystyle \frac{1}{27.4}}{g}_2^3>0$, the relation (234) between these roots with the values of the Weierstrass function at the points on the complex plane ${\omega }_1,{\omega }_2,{\omega }_3={\omega }_1+{\omega }_2$ (for convenience we omit the factor ${\displaystyle \frac{1}{2}}$ in the argument of the Weierstrass function $\rho \left(z\right)$)

$$e_1=\rho \left({\omega }_1\right)={\alpha }_1,e_2=\rho \left({\omega }_2\right)={\beta }_1,e_3=\rho ({\omega }_1+{\omega }_2)={\gamma }_1.$$

(316)

and also the group-theoretical law for summation of the points ${\omega }_1,{\omega }_2,{\omega }_3={\omega }_1+{\omega }_2$ on the elliptic curve, expressed analytically by the formulae (241)

$$\rho (z_1+z_2)=\rho \left(z_1\right)+\rho \left(z_2\right)+{\displaystyle \frac{1}{4}}{\left({\displaystyle \frac{{\rho }^{\prime }\left(z_1\right)-{\rho }^{\prime }\left(z_2\right)}{\rho \left(z_1\right)-\rho \left(z_2\right)}}\right)}^2.$$

Taking into account (316), let us take the ratio ${\displaystyle \frac{\rho \left(z_1\right)}{\rho \left(z_2\right)}}$ of the Weierstrass function at the two different values $z_1={\omega }_1$ and $z_2={\omega }_2$ and also the ratio ${\displaystyle \frac{\rho \left(z_1\right)}{\rho \left(z_3\right)}}$ at the two values $z_1={\omega }_1$ and $z_3={\omega }_3={\omega }_1+{\omega }_2$

$${\displaystyle \frac{\rho \left(z_1\right)}{\rho \left(z_2\right)}}={\displaystyle \frac{{\alpha }_1}{{\beta }_1}},{\displaystyle \frac{\rho \left(z_1\right)}{\rho \left(z_3\right)}}={\displaystyle \frac{{\alpha }_1}{{\gamma }_1}},$$

(317)

where the root ${\alpha }_1=A_1+B_1$ is given by formulae (277) and the complex roots ${\beta }_1$ and ${\gamma }_1$ are given by formulaes (281) and (282) respectively. Since ${\beta }_1$ is a complex root, the Weierstrass function $\rho \left(z_2\right)={\beta }_1$ is also a complex function and can be represented as $\rho \left(z_2\right)={\rho }_1\left(z_2\right)+i_2\rho \left(z_2\right)$, where ${\rho }_1\left(z_2\right)$ is the real part of $\rho \left(z_2\right)$ and ${\rho }_2\left(z_2\right)$ is the imaginary part.

10.11.1. The First Equation

The first equation in (317) can be written as

$${\displaystyle \frac{\rho \left(z_1\right)}{{\rho }_1\left(z_2\right)+i{\rho }_2\left(z_2\right)}}={\displaystyle \frac{A_1+B_1}{\frac{-(A_1+B_1)}{2}+i\frac{\sqrt[.]{3}}{2}(A_1-B_1)}}.$$

(318)

Denoting

$$m:={\displaystyle \frac{{\tilde{A}}_1-{\tilde{B}}_1}{{\tilde{A}}_1+{\tilde{B}}_1}},$$

(319)

the above equation can be represented as

$$-{\displaystyle \frac{{\rho }_{.}\left(z_1\right)}{2}}-{\rho }_1\left(z_2\right)+i{\displaystyle \frac{\sqrt[.]{3}}{2}}{\rho }_{.}\left(z_1\right)m-i{\rho }_2\left(z_2\right)=0.$$

(320)

Setting up equal to zero the real and imaginary part, one can obtain

$${\rho }_{1.}\left(z_2\right)=-{\displaystyle \frac{{\rho }_{.}\left(z_1\right)}{2}},{\rho }_{2.}\left(z_2\right)={\displaystyle \frac{\sqrt[.]{3}}{2}}{\rho }_{.}\left(z_1\right)m.$$

(321)

10.11.2. The Second Equation

Let us now make use of the second equation in (317), which can be written as

$${\displaystyle \frac{\rho \left(z_1\right)}{{\rho }_{.}\left(z_3\right)}}={\displaystyle \frac{A_1+B_1}{-\frac{(A_1+B_1)}{2}-i\frac{\sqrt[.]{3}}{2}(A_1-B_1)}}.$$

(322)

Using the group-theoretical law for summing up of points, represented by formulae (241) for expressing ${\rho }_{.}\left(z_3\right)={\rho }_{.}(z_1+z_2)$, the above equation can be written as

$${\displaystyle \frac{\rho \left(z_1\right)}{\rho \left(z_1\right)+\rho \left(z_2\right)+\frac{1}{4}{\left(\frac{{\rho }^{\prime }\left(z_1\right)-{\rho }^{\prime }\left(z_2\right)}{\rho \left(z_1\right)-\rho \left(z_2\right)}\right)}^2}}=$$

$$={\displaystyle \frac{1}{-\frac{1}{2}-i\frac{\sqrt[.]{3}}{2}m}}.$$

(323)

Using the two equalities in (321), the equation (323) can be represented as $Q_1+i{Q}_2=0$, where the real part $Q_1$ and the complex part $Q_2$ are

$$Q_1:=-3{\rho }^3\left(z_1\right)(5m^2+3)-{\displaystyle \frac{9}{4}}{\left({\rho }^{\prime }\left(z_1\right)\right)}^2(1-{\displaystyle \frac{1}{3}}{m}^2)$$

$$+{\displaystyle \frac{3}{4}}\rho \left(z_1\right){\displaystyle \frac{\partial m}{\partial q}}.{\displaystyle \frac{\partial q}{\partial z}}+{\displaystyle \frac{3}{2}}\rho \left(z_1\right){\rho }^{\prime }\left(z_1\right)m{\displaystyle \frac{\partial m}{\partial q}}.{\displaystyle \frac{\partial q}{\partial z}}=0,$$

(324)

$$Q_2:={\rho }^3\left(z_1\right)m(m^2-1)+{\displaystyle \frac{1}{2}}{\rho }^{\prime }\left(z_1\right)\left[{\rho }^{\prime }\left(z_1\right)m+\rho \left(z_1\right){\displaystyle \frac{\partial m}{\partial q}}.{\displaystyle \frac{\partial q}{\partial z}}\right]=0.$$

(325)

Since the equation $Q_1+i{Q}_2=0$ is fulfilled when the real and imaginary part are equal to zero, we have to set up $Q_1=0$ and $Q_2=0$. Also, we have assumed that the modulus parameter q of the elliptic integral in the Legendre form might depend on the complex coordinate z, on which by definition depends the Weierstrass function $\rho \left(z\right)$. However, in view of the complexity of the equations, for now we shall consider the simplified case ${\displaystyle \frac{\partial q}{\partial z}}=0$.

10.11.3. First Way for Representation of $m^2$ and of ${\left({\rho }^{\prime }\left(z_1\right)\right)}^2$

Further, from the second equation (325) for $Q_2$, one may express the ratio (319) $m^2:={\left({\displaystyle \frac{{\tilde{A}}_1-{\tilde{B}}_1}{{\tilde{A}}_1+{\tilde{B}}_1}}\right)}^2$ as

$$m^2={\displaystyle \frac{{\rho }^3\left(z_1\right)-\frac{1}{2}{\left({\rho }^{\prime }\left(z_1\right)\right)}^2}{{\rho }^3\left(z_1\right)}}.$$

(326)

Substituting $m^2$ from (326) in the first equation (324) for $Q_1$ and denoting ${\left({\rho }^{\prime }\left(z_1\right)\right)}^2:=X$, one obtains the following quadratic equation with respect to the variable X

$$9X^2-18.8{\rho }^3\left(z_1\right)X+72{\rho }^6\left(z_1\right)=0.$$

(327)

It is known that the roots of the simple quadratic equation $x^2+px+q=0$ are found from the formulae [53]

$$x=-{\displaystyle \frac{p}{2}}\pm \sqrt[.]{{\displaystyle \frac{p^2}{4}}-q}.$$

(328)

Consequently, the roots of the quadratic equation (327) are

$${\left({\rho }^{\prime }\left(z_1\right)\right)}^2:=X=18.4{\rho }^3\left(z_1\right)\pm \sqrt[.]{{\displaystyle \frac{{\left(18.8\right)}^2{\rho }^6\left(z_1\right)}{4}}-72}.$$

(329)

This relation can be substituted in formulae (326) for $m^2$ and for the moment investigating only the case of a positive sign in front of the square root, $m^2$ can be written as

$$m^2=-35-3\sqrt[.]{{12}^2-{\displaystyle \frac{2}{{\rho }^6\left(z_1\right)}}}.$$

(330)

Respectively, the solution (329) for ${\left({\rho }^{\prime }\left(z_1\right)\right)}^2$ can be represented as

$${\left({\rho }^{\prime }\left(z_1\right)\right)}^2=72{\rho }^3\left(z_1\right)+6{\rho }^3\left(z_1\right)\sqrt[.]{{12}^2-{\displaystyle \frac{2}{{\rho }^6\left(z_1\right)}}}.$$

(331)

The last two relations are the first way for representing $m^2$ and ${\left({\rho }^{\prime }\left(z_1\right)\right)}^2$. Now it shall be proved that there is a second way for representing these expressions.

10.11.4. Second Way for Representation of $m^2$ and of ${\left({\rho }^{\prime }\left(z_1\right)\right)}^2$

The second way is based on the standard formulae for the parametrizable form of a cubic curve

$${\left({\rho }^{\prime }\left(z_1\right)\right)}^2=4{\rho }^3\left(z_1\right)-g_2\rho \left(z_1\right)-g_3.$$

(332)

Substituting into (326) for $m^2$, we obtain the second representation

$$m^2=-2+{\displaystyle \frac{g_2}{2}}.{\displaystyle \frac{1}{{\rho }^2\left(z_1\right)}}+{\displaystyle \frac{g_3}{2}}{\displaystyle \frac{1}{{\rho }^3\left(z_1\right)}}.$$

(333)

Next we shall use also expression (300) for $g_3$ through $g_2$

$$g_3=g_{2}N\left(q\right),N\left(q\right):={\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{56}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{56}}.}.$$

Taking this equality into account and also that according to (302) $N\left(q\right)=8F\left(q\right)$, the second representation (333) for $m^2$ can be written also in the form

$$m^2={\displaystyle \frac{-{\rho }^3\left(z_1\right)+\left[\frac{\rho \left(z_1\right)}{2}+4F\left(q\right)\right]g_2}{{\rho }^3\left(z_1\right)}}.$$

(334)

Setting up equal the two representations (330) and (333) for $m^2$, one can obtain

$$-34{\rho }^3\left(z_1\right)-3{\rho }^3\left(z_1\right)\sqrt[.]{{12}^2-{\displaystyle \frac{2}{{\rho }^6\left(z_1\right)}}}$$

$$={\displaystyle \frac{g_2}{2}}(\rho \left(z_1\right)+8F\left(q\right)).$$

(335)

10.11.5. Another Interpretation of the Equation (335) for the Equality of the Two Representations of $m^2$. A Theorem from Higher Algebra for Investigating Whether There are Roots of an Algebraic Polynomial

After some transformations, the above equation can be represented in the form of a sixth-order algebraic equation with respect to $\rho \left(z_1\right)$

$$\left[9.{12}^2\right]{\rho }^6\left(z_1\right)-34g_2{\rho }^4\left(z_1\right)-68.4F\left(q\right){\rho }^3\left(z_1\right)$$

$$-g_2^2\left[{\rho }^2\left(z_1\right)-4F\left(q\right)\rho \left(z_1\right)+16F^2\left(q\right)\right]=0.$$

(336)

This equation of the $6-$th degree is interesting, because it does not always have roots. It will, however, have roots for certain values of the invariant $g_2$. Thus, the values of $\rho \left(z_1\right)$ are interrelated to the values of $g_2$. It is interesting to investigate the following problem: according to (335), for a given function $\rho \left(z_1\right)$ there should be only one value of $g_2$. However, the algebraic equation (336), if solved as a quadratic equation with respect to $g_2$ and for a given function $\rho \left(z_1\right)$, should have two roots with respect to $g_2$, which will contradict the unique expression of $g_2$ from (335).

Now let us give the formulation of the theorem of Schur from higher algebra, which gives the possibility to establish whether a polynomial has roots. The formulation of the theorem, originally published by Schur [55] yet in 1918, is taken from the monograph [53], where many other theorems, concerning roots of polynomials can be found. The theorem of Schur has been applied in the paper [18] for establishing whether two different algebraic equations of the fourth degree have roots.

The first of these equations was related to the introduced in [18] physical notion of "space-time interval", which for a particular case was proved to be positive, negative or even equal to zero. By means of the Schur theorem, it was proved that the corresponding algebraic equation has roots.

The second equation was related to the notion of "geodesic distance", which is related to the propagation of a light or radio signal in a gravitational field and therefore, should be positive. Respectively, again by applying the Schur theorem it was proved that the corresponding algebraic equation does not have roots. The confirmation of the physical definitions has been commented in the paper [19] and also the review paper [22]. This consistency between the mathematical results and the corresponding physical definitions give some grounds to believe that the Schurs theorem is a reliable mathematical tool for investigation also of equation (336).

Theorem 14.

([55]) The necessary and sufficient conditions for the polynomial of $n-$th degree

$$f\left(y\right)=a_0{y}^n+a_1{y}^{n-1}+....+a_{n-2}{y}^2+a_{n-1}y+a_n$$

(337)

to have roots only in the circle $\mid y\mid <1$ are the following ones:

1. The fulfillment of the inequality

$$\mid a_0\mid >\mid a_n\mid .$$

(338)

2. The roots of the polynomial of $(n-1)$ degree

$$f_1\left(y\right)={\displaystyle \frac{1}{y}}\left[a_{0}f\left(y\right)-a_n{f}^{*}\left(y\right)\right]$$

(339)

should be contained in the circle $\mid y\mid <1$, where $f^{*}\left(y\right)$ is the s.c. "inverse polynomial", defined as

$$f^{*}\left(y\right)=y^{n}f\left({\displaystyle \frac{1}{y}}\right)=a_n{y}^n+a_{n-1}{y}^{n-1}+....+a_2{y}^2+a_{1}y+a_0.$$

(340)

In the case of fulfillment of the inverse inequality

$$\mid a_0\mid <\mid a_n\mid ,$$

(341)

the $(n-1)$ degree polynomial $f_1\left(y\right)$ (again with the requirement the roots to remain within the circle $\mid y\mid <1$) is given by the expression

$$f_1\left(y\right)=a_{n}f\left(y\right)-a_0{f}^{*}\left(y\right).$$

(342)

In the paper [18] the Schur theorem has been generalized for a "chain" of algebraic equations, when it is applied to a chain of algebraic equations with diminishing degrees. Many other theorems for investigating when a polynomial has roots or does not have can be found also in the contemporary monograph [57].

10.12. Finding the Weierstrass Function $\rho \left(z_1\right)$ and the Invariant $g_2$

Equation (335) further will play an important role, since it allows expressing the invariant $g_2$ as a function of the Weierstrass function $\rho \left(z_1\right)$. More importantly, it shows that for a given function $\rho \left(z_1\right)$, there is only one value of $g_2$.

Yet, eq. (335) is not sufficient for finding both the Weierstrass invariant $g_2$ and the Weierstrass function $\rho \left(z_1\right)$. We need one more equation and for the purpose we will use the defining equality for m (319) $m:={\displaystyle \frac{{\tilde{A}}_1-{\tilde{B}}_1}{{\tilde{A}}_1+{\tilde{B}}_1}}$, where ${\tilde{A}}_1$ and ${\tilde{B}}_1$ were previously defined by the formulaes (304)

$${\tilde{A}}_1:=\sqrt[3]{F\left(q\right)+{\displaystyle \frac{1}{4}}\sqrt[.]{\tilde{\Delta }}},{\tilde{B}}_1:=\sqrt[3]{F\left(q\right)-{\displaystyle \frac{1}{4}}\sqrt[.]{\tilde{\Delta }}}.$$

Thus it can easily be obtained that ${\displaystyle \frac{{\tilde{A}}_1}{{\tilde{B}}_1}}={\displaystyle \frac{1+m}{1-m}}$ and also, with account of (304) that

$${\displaystyle \frac{F\left(q\right)+\frac{1}{4}\sqrt[.]{\tilde{\Delta }}}{F\left(q\right)-\frac{1}{4}\sqrt[.]{\tilde{\Delta }}}}={\displaystyle \frac{{(1+m)}^3}{{(1-m)}^3}},$$

(343)

where $\tilde{\Delta }$ has been defined by the equality (303)

$$\tilde{\Delta }:={\displaystyle \frac{1}{4}}{\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{28}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{28}}}-{\displaystyle \frac{1}{4.27}}{g}_2$$

Denoting the first term in the above equality as ${\tilde{\Delta }}_1:=$ ${\displaystyle \frac{1}{4}}{\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{28}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{28}}}$ and remembering the defining equality (302)

$$F\left(q\right):={\displaystyle \frac{1}{8}}{\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{56}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{56}}}={\displaystyle \frac{1}{8}}N\left(q\right),$$

one can obtain the simple equality

$${\displaystyle \frac{{\tilde{\Delta }}_1}{F\left(q\right)}}=2{\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{56}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{56}}.}=16F\left(q\right).$$

(344)

Consequently, the discriminant $\tilde{\Delta }$ (303) can be represented as

$$\tilde{\Delta }:={\tilde{\Delta }}_1-{\displaystyle \frac{1}{4.27}}{g}_2=16F^2\left(q\right)-{\displaystyle \frac{1}{4.27}}{g}_2.$$

(345)

Making use of (343), one can find an expression for $\sqrt[.]{\tilde{\Delta }}$

$$\sqrt[.]{\tilde{\Delta }}=4F\left(q\right){\displaystyle \frac{\left[{(1+m)}^3-{(1-m)}^3\right]}{\left[{(1+m)}^3+{(1-m)}^3\right]}}.$$

(346)

Transforming this equality and using (345), one obtains the final formulae for $\tilde{\Delta }$, from where $g_2$ can be expressed and $m^2$ can be expressed by two different ways

$$\tilde{\Delta }=16F^2\left(q\right)-{\displaystyle \frac{1}{4.27}}{g}_2={\displaystyle \frac{16F^2\left(q\right)m^2{(3+m^2)}^2}{{(1+3m^2)}^2}}.$$

(347)

In this formulae we substitute $m^2$, expressed in (333) through the Weierstrass invariant $g_2$ and higher powers of the Weierstrass function $\rho \left(z\right)$. After some long calculations, the following complicated expression is obtained

$${\displaystyle \frac{1}{4.27}}{g}_2=16F\left(q\right)+{\displaystyle \frac{16F^2\left(q\right){\rho }^6\left(z_1\right)}{{\rho }^9\left(z_1\right)}}{\displaystyle \frac{P_1(\rho \left(z_1\right),q)}{P_2(\rho \left(z_1\right),q)}},$$

(348)

where$P_1(\rho \left(z_1\right),q)$ and $P_2(\rho \left(z_1\right),q)$ are the expressions (note that everywhere the dot "." means multiplication, for example $35.{32}^2=35\times {32}^2$)

$$P_1(\rho \left(z_1\right),q):=-35.{32}^2{\rho }^9\left(z_1\right)+32.306{\rho }^6\left(z_1\right)\sqrt{{12}^2{\rho }^6\left(z_1\right)-2}$$

$$+\left(99.9\right){\rho }^3\left(z_1\right)\left({12}^2{\rho }^6\left(z_1\right)-2\right)+27{({12}^2{\rho }^6\left(z_1\right)-2)}^{{\displaystyle \frac{3}{2}}},$$

(349)

$$P_2(\rho \left(z_1\right),q):={104}^2{\rho }^6\left(z_1\right)+\left(104.18\right){\rho }^3\left(z_1\right)\sqrt{12{\rho }^6\left(z_1\right)-2}$$

$$+81(12{\rho }^6\left(z_1\right)-2).$$

(350)

This equation and also eq. (335)

$$-34{\rho }^3\left(z_1\right)-3{\rho }^3\left(z_1\right)\sqrt[.]{{12}^2-{\displaystyle \frac{2}{{\rho }^6\left(z_1\right)}}}$$

$$={\displaystyle \frac{g_2}{2}}(\rho \left(z_1\right)+8F\left(q\right)).$$

can be considered as a system of two equations with respect to $g_2$ and $\rho \left(z_1\right)$. If from (335) the Weierstrass invariant $g_2$ and is expressed and then substituted into (348), then a complicated equation is obtained only with respect to $\rho \left(z_1\right)$. It is not clear how to find its roots (and also how many are they), but if $\rho \left(z_1\right)$ is found, then $g_2$ is uniquely expressed from (335).

There is one more possibility, namely- to use the second representation for $m^2$ by formulae (334), which is to be substituted in equation (347) for $\tilde{\Delta }$. We shall not present here this method, which is more complicated and will lead to the derivation of a cubic equation with respect to $g_2$, which after solving has to be substituted in the earlier derived algebraic equation (336) of sixth order with respect to $\rho \left(z_1\right)$. However, again a higher-orded algebraic equation has to be solved, so this method is not effective.

11. Discussion

In this paper the propagation time of a signal was analytically calculated by making us of the null cone equation from the metric element for the near-Earth space. Two different cases were considered, corresponding to two different parametrizations of the satellite orbit: the first parametrization is in terms of the eccentric anomaly angle E, which is the dynamical parameter for the case of a plane elliptic orbit and the second parametrization in terms of all the 6 Kepler parameters $(f_{.},a_{.},e_{.},{\Omega }_{.},I_{.},{\omega }_{.})$ for the case of space-oriented orbit with only the true anomaly angle f considered to be the dynamical parameter, associated with the motion of the satellite. The propagation time in the second case is expressed in more complicated elliptic integrals of the second- and of the fourth- order, but in both cases, the sum of the $O\left(c^{-1}\right)$ and $O\left(c^{-3}\right)$ propagation time corrections are real-valued expressions. This proves the correctness of the theoretical approach of parametrizing the space coordinates $x,y$ or $x,y,z$ in the null cone equation with the orbital coordinates, known from celestial mechanics. Numerical calculations for the propagation time for the first case of plane elliptical orbit confirm that during the propagation time the signal travels much greater distance than the distance travelled by the satellite. It turned out also that theorems form General Relativity Theory [29] and from differential geometry [30,32] confirm the conclusion that various parametrizations of the space coordinates, related to the orbital motion of the satellite can be used.

It is important to clarify that the calculated propagation time gives the correspondence between the initial and final values of the dynamical parameter (the eccentric anomaly angle E or the true anomaly angle f) and the initial $T_{init}$ and final $T_{fin}$ moments of propagation of the signal. However, in such a theoretical set-up $T_{init}$ and $T_{fin}$ are not related to any real processes of transmission of signals between satellites on one and the same orbit or on different space-oriented orbits. In such a modified setting, the initial position of the signal-emitting satellite and the final position of the signal-receiving satellite have to be given and they will correspond to the initial $T_{init}$ and final $T_{fin}$ moments of propagation of the signal. They have to be calculated in view of the fact that the trajectory of the signal is curved due to the action of the gravitational field. Due to this, the propagation time will not be the (presumed to be known at each moment of time) Euclidean distance between the satellites, divided by the velocity of light. The nontrivial moment in such a problem is that both satellites are moving (at the moments of emission and receiving of the signal and also during the time of propagation of the signal) so that the moments of emission and reception of the signal (they are denoted by $T_1$ and $T_2$) should incorporate these important peculiarities. In other words, the curvature of the signal trajectory will have to be such, that at a certain space point of the second satellite the signal will be intercepted, provided also that the signal is emitted at the proper position by the first (signal emitting) satellite.

In the papers [18,19] and the review paper [22] it was suggested that $T_1$ and $T_2$ should be related to one another in the framework of a more general model, called "the two intersecting null cones model". In the Introduction of this paper a short summary is given of this formalism, based on algebraic geometry and elliptic functions and integrals. However, this paper concentrates only on the methods for calculation of the propagation time (related to elliptic functions and integrals) and not on the two null cones formalism. This is the first step for creating the theoretical formalism in a more general setting with the purpose to establish whether the introduced in the papers [18] and [19] physical notions of "space-time distance", "condition for inter-satellite communications" and "geodesic distance" will be valid also for the more general case of space-oriented orbits.

The last section, named "New Approach for Analytical Calculation of Elliptic Integrals. Analytic Relation Between Elliptic Integrals in the Weierstrass and in the Legendre Form" is dedicated to finding new solutions of elliptic integrals in the Legendre form by means of transforming them to elliptic integrals in the Weierstrass form, which after inversion have as solution the Weierstrass elliptic function. The necessity to study (in the pure mathematical aspect) such integrals is dictated not only by their important application in GPS inter-satellite communications with account also of the General Relativity effects, but also by the wide application of elliptic integrals in various problems of mechanics, cosmology, black hole physic,integrable systems and etc. A short review of these different applications has been given in the paper [21].

In this paper (as well in [21]) the representation of the integral in two different ways is exposed, one of them is based on the transformation $x={\displaystyle \frac{a}{y^2}}+b$ ($a,b$ are proved to be functions of the modulus parameter q of the integral in the Legendre form), and the second method is related to the s.c. "four-dimensional uniformization", exposed in the monograph by Whittaker, Watson. [60]. After comparison of the Weierstrass invariants in the two representations, the main result is the representation (299) for the ratio of the two invariants

$${\displaystyle \frac{g_3}{g_2}}={\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{56}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{56}}},$$

where $K\left(q\right)$ is the function (178) $K\left(q\right)\equiv {\displaystyle \frac{q^4-q^2+1}{2q^4-5q^2+2}}$, ${\overline{g}}_2$ and ${\overline{g}}_3$ are the Weierstrass invariants, obtained after the four-dimensional uniformization [60]. This result is not trivial and previously not investigated in the monographs on elliptic functions. Further it was proved also that the Weierstrass invariant $g_2$ and the Weierstrass function $\rho \left(z\right)$ can be found from a complicated system of two nonlinear equations.

Another interesting result was obtained in this paper, suggested by a theorem in the monograph by Koblitz [49]. If $g_2$ is a real number (or a real and positive function) and there are one real root and two imaginary ones for the cubic polynomial $4x^3-g_{2}x-g_3$ (the case investigated in this paper), then the other invariant $g_3$ should be imaginary and the discriminant $\Delta =g_2^3-27g_3^2$ should be negative (a case not investigated in this paper). It should be kept in mind also that the case "negative discriminant and imaginary $g_3$" does not follow directly from the formulation of the theorem, but rather than from the "non-fulfillment" of the theorem, which is an unusual moment. No doubt, more strict mathematical treatment is necessary.

The case of a negative discriminant would correspond to three real roots of the cubic polynomial (see [53] and [54]) and this case needs a separate investigation. However, the theorem in [49] has a very specific and intricate formulation: from it, for example follows that for a positive discriminant one may have also a case with three real roots of the polynomial (contrary to the standard case in the algebraic books [53] and [54]), but then the invariants $g_2$ and $g_3$ should also be real! So the interesting fact is that the Weierstrass invariants $g_2$ and $g_3$ have a purely "algebraic" meaning whether the roots are only real or one real and two complex. Consequently, the theorem in [49] should be reformulated for all the possible cases and investigated again.

Transformations of the elliptic integral with a general fourth-order polynomial under the root in the denominator into an integral in the Legendre form have been studied also in the monograph [52] (also in [67]), but in [52] transformation to an integral in the Weierstrass form has not been performed. However, in [67] the inverse transformation from an integral in the Weierstrass form with the polynomial

$$Y^2=(x-\alpha )(x-\beta )(x-\gamma )$$

(351)

to the integral in the Legendre form [67]

$$\sqrt{\alpha -\gamma }{\displaystyle \frac{dx}{Y}}={\displaystyle \frac{2dy}{\sqrt{(1-y^2)(1-{\tilde{q}}^2{y}^2)}}},$$

(352)

where the transformation

$$x-\alpha ={\displaystyle \frac{(\alpha -\beta )y^2}{(1-y^2)}}$$

(353)

has been applied and the modulus of the elliptic integral $\tilde{q}$ in the Legendre form was determined by the formulae

$${\tilde{q}}^2={\displaystyle \frac{\beta -\gamma }{\alpha -\gamma }}.$$

(354)

An interesting problem for future research is the following: the ratio (299) of the invariants $g_2$ and $g_3$ can be used for calculating the roots $\alpha ,\beta ,\gamma$. The resulting expressions will depend on the invariant $g_2$ and can be substituted in the formulae (354) for ${\tilde{q}}^2$. Then it is interesting to see whether the invariant $g_2$ can be determined in such a way that the calculated ${\tilde{q}}^2$ modulus according to (354) will coincide with the initially given q in the integral ${\displaystyle \frac{dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}$.

There are also another problems, related to the implementation of the Schur theorem from higher algebra. Besides the ratio (299), another ratio

$${\displaystyle \frac{g_3}{g_2}}=\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right){\left(9K\left(q\right)\right)}^{-{\displaystyle \frac{1}{3}}.}$$

(355)

can be obtained, if assumed that not the roots ${\alpha }_1$ and ${\alpha }_2$ coincide ( ${\alpha }_1={\alpha }_2$), but the roots ${\beta }_1$ and ${\beta }_2$ coincide, i.e. ${\beta }_1={\beta }_2$ This case also is not considered in this paper, but the investigation should be performed in the same manner.

Also, an important problem for further investigation is the algebraic treatment of the higher-order algebraic equation for the discriminant (278) $\Delta ={\displaystyle \frac{g_3^2}{4}}-{\displaystyle \frac{1}{108}}{g}_2^3$ of the cubic polynomial in the Weierstrass elliptic integral, which in this investigation was assumed to be positive.