Numerical Modelling of Satellite Downlink Signals in a Finslerian-Perturbed Schwarzschild Spacetime

The work presented in this paper aims to contribute to the problem of testing Finsler gravity theories by means of experiments and observations in the solar system. Within a class of spherically symmetric static Finsler spacetimes we consider a satellite with an on-board atomic clock, orbiting in the Finslerian-perturbed gravitational field of the earth, whose time signal is transmitted to a ground station, where its receive time and frequency are measured with respect to another atomic clock. This configuration is realized by the Galileo 5 and 6 satellites that have gone astray and are now on non-circular orbits. Our method consists in the numerical integration of the satellite’s orbit, followed by an iterative procedure which provides the numerically integrated signals, i.e., null geodesics, from the satellite to the ground station. One of our main findings is that for orbits that are considerably more eccentric than the Galileo 5 and 6 satellite orbits, Finslerian effects can be separated from effects of perturbations of the Schwarzschild spacetime within the Lorentzian geometry. We also discuss the separation from effects of non-gravitational perturbations. This leads us to the conclusion that observations of this kind combined with appropriate numerical modelling can provide suitable tests of Finslerian modifications of general relativity.


Introduction
In this paper, we present a method of testing Finsler gravity, a generalization of Einstein's general theory of relativity, by means of experiments and observations in the solar system. We mainly focus on the case that the time signal of an atomic clock on board an Earth satellite is transmitted to a ground station, where it is compared with the time signal of another atomic clock. The Finslerian perturbations of the Schwarzschild spacetime affect the satellite orbit and the downlink signals and cause small shifts of the receive times and frequencies of the downlink signals, compared to the unperturbed Schwarzschild spacetime.
Let us first make some remarks concerning the role of such experiments in testing the Einsteinian theory or generalized theories of gravity. One of the "three classical tests of general relativity" consists in verifying the gravitational redshift, which can be described as follows. If a source of electromagnetic radiation is situated in a gravitational field, its spectrum measured outside or less deep in the field is shifted to longer wavelengths compared to its spectrum measured nearby the source. In other words, a standard clock in a gravitational field appears to run slower when observed from outside or less deep in the field. Assuming that atomic clocks are standard clocks, there is the possibility to test general relativity by comparing the time signals of two atomic clocks, one on board an Earth satellite and one on the ground.
Recently, the accuracy of such experiments has been considerably improved by using Galileo on-board atomic clocks. On 22 August 2014, the satellites Galileo 5 and 6 were unintentionally launched into eccentric orbits. From the point of view of gravitational physics, this mishap was a stroke of luck. Thanks to this, a little more than four years later, Delva et al. [1] and Hermann et al. [2] were able to report on confirmation of the prediction from general relativity with an accuracy in the order of 10 −5 at 1 σ level. This is a few times more accurate than the Gravity Probe A experiment in 1976, see Vessot et al. [3], which had previously provided the most accurate confirmation of the general relativistic redshift effect. The accuracy now achieved opens up the possibility of not only testing Einstein's theory, but also of providing quantitative bounds on alternative theories of gravity by means of such experiments.
Finsler gravity is a generalization of general relativity, where the components of the metric tensor can depend not only on the spacetime coordinates but also on the components of a tangent vector. More precisely, the Lagrangian of the Finsler geodesics must no longer be bilinear in the components of the tangent vector but only fulfil the weaker requirement of homogeneity of degree two. Since bilinearity is a special case of homogeneity of degree two, Lorentzian geometry is a special case of Finsler geometry. Consequently, general relativity is a special case of Finsler gravity.
With regard to the position of Finsler gravity among the alternative theories of gravity, it should first be mentioned that it breaks spatial isotropy even on the tangent spaces. Therefore, the principle of local Lorentz invariance (LLI) is violated and the propagation of light in vacuum may be anisotropic. Since LLI is part of the Einstein equivalence principle (EEP), the EEP is also violated in Finsler gravity. (As is well known, the validity of the EEP requires a pseudo-Riemannian spacetime geometry, see e.g., Will [4]). However, since in our class of Finsler spacetimes (see Section 2) there is a unique timelike geodesic for every timelike initial condition, the weak equivalence principle is valid.
Since length is not a Lorentz invariant quantity, the above-mentioned violation of LLI can also be expected in a still to be found theory of quantum gravity which asserts that there is a fundamental length scale given by the Planck length. This motivates to consider non-quantized theories of Finsler gravity as possible "interpolations" between general relativity and quantum gravity in some regime of length and energy scales while maintaining the weak equivalence principle. More about this and other motivations, so in connection with the Ehlers-Pirani-Schild [5] axiomatic approach to general relativity, can be found in the recent review on the present status and the perspectives of Finslerian spacetime theories by Lämmerzahl and Perlick [6] and in Section I of the article by Lämmerzahl, Perlick, and Hasse [7]. The latter publication also provides the class of spherically symmetric static Finsler spacetimes on which our studies are based. The Finsler metrics of this model are purely kinematical perturbations of the Schwarzschild metric, i.e., these were not derived as solutions to any field equation.
Gravitational perturbations by higher multipole moments of the central body (a planet or the sun) or by other bodies (e.g., the sun and the moon in the case of a satellite orbiting the earth, the planets in the case of a test particle in the gravitational field of the sun) as well as non-gravitational perturbations (e.g., solar radiation pressure) have not been taken into account. This is justified by the assumption that all these perturbations, as well as the Finslerian perturbations, are so small that in principle one may linearize the equations of motion with respect to the corresponding terms and treat all these perturbations separately, followed by addition of all their effects on measurable quantities. (In this context, we note that in the already mentioned analyses [1,2] of Galileo 5 and 6 data, to test the gravitational redshift predicted by Einstein's theory, from the here mentioned perturbations only Earth's oblateness, represented by the zonal coefficient J 2 , has been considered.) Following this strategy, we restrict our analysis to perturbations which can be modelled by three "perturbation functions" which respect the stationarity and the spherical symmetry of the spacetime. These perturbation functions will be introduced in the next section.
The paper is organized as follows. In Section 2, we will fix our notations and conventions and introduce the considered class of Finsler spacetimes. Section 3 gives a brief description of our model of a downlink from an Earth satellite to a ground station. The subsequent Section 4 deals with the mathematical representation of our model and its numerical treatment. Section 5 presents tests of our software against results published in [7]. In Section 6 we explain the iterative procedure which provides the null geodesics of the downlink. Readers who are less interested in the numerical aspects may skip Sections 4-6. In order to separate the effects of the orbit perturbations from the effects of the signal perturbations we consider in Section 7 "hybrid models", in which both kinds of perturbations can be "switched on" and "off" separately. The most important part of our paper is Section 8. It contains core findings obtained by means of our simulations. We discuss, exemplarily for the orbits of the Galileo 5 and 6 satellites, to which extent Finslerian effects can be separated from effects of perturbations of the Schwarzschild spacetime within the Lorentzian geometry. Based on these results, we consider, in Section 9, a model satellite with high orbit eccentricity. In the concluding Section 10, we discuss, among other things, the advantages of combining the two Galileo satellites with a future dedicated mission and provide an outlook how our numerical approach can be brought over to other situations of communication or even to other modifications of general relativity.

A Class of Spherically Symmetric Static Finsler Spacetimes
We work in the same class of Finsler spacetimes as in the paper of Lämmerzahl, Perlick, and Hasse [7]. This class is defined in Schwarzschild coordinates t, r, ϑ, and ϕ by the Finsler Lagrangian see (14) in [7] or (37) in [8], which defines a perturbed Schwarzschild spacetime. Here h tt and h rr are the Schwarzschild metric coefficients, given by where c is the speed of light, G is the gravitational constant, M is the mass of the gravitating body, and r S = 2GM/c 2 is its Schwarzschild radius. The metric convention is (-+ + +), i.e., the unperturbed metric has signature +2. φ 0 , φ 1 , and φ 2 are "perturbation functions" which depend differentiably only on the radial coordinate r. For stationary observers, φ 0 and φ 1 perturb the measurement of proper time and radial length, respectively. Whereas φ 0 and φ 1 describe perturbations within the Lorentzian geometry, φ 2 introduces a spatial anisotropy which is a genuine Finsler feature.
Let us remark that it results from sections III and IV of [7] that the choice of the Finsler Lagrangian (1) is in some sense natural. To be more precise, it results from adding the leadingorder spherical-harmonics term in a general Finslerian perturbation of the spatial part of the metric to the Lagrangian of the Schwarzschild metric, followed by linearisation with respect to the perturbing terms. Thereby, "general" means only restricted to compatibility with the spherical symmetry.
Throughout this paper, as in [7], we assume that the perturbation functions are so small that we may linearize all equations with respect to the φ n (r) and their derivatives φ n (r), n = 0, 1, 2. The outputs of our numerical codes are also based on this linearisation procedure.
In [7] only problems are considered that can be reduced to the equatorial plane ϑ = π/2 by spherical symmetry. But in our case, when it comes to modelling the downlink, all three spatial dimensions must be taken into account.
The first order equations of motion in the plane ϑ = π/2, with the constants of motion as parameters, are derived in Section IV of reference [7]. From that one can easily derive the equations of motion for given initial values in the four-dimensional spacetime by means of a suitable spatial rotation of the coordinate system.
For later reference we note how φ 0 , φ 1 , and φ 2 affect the propagation of material bodies and light rays as seen from a stationary observer at infinity (at least in linear approximation). This can be seen from the equations of motions in the Appendixes A and B. In general, φ 0 , φ 1 , and φ 2 affect both the radial and the tangential velocity components. However, there are some special cases: • φ 0 , φ 1 , and φ 2 never cause tangential acceleration components if the velocity is pure tangential or pure radial.

Our Model of a Downlink from an Earth Satellite to a Ground Station
Let us assume that the geometry outside the earth is determined by the Finsler Lagrangian (1) where, in this case, M is the earth's mass. We consider an Earth satellite whose orbit is modelled by a spatially bounded timelike geodesic. An on-board clock is modelled by the satellite's Finsler proper time (with arbitrary zero point). Furthermore, we consider a ground station on the terrestrial surface. The rotation of the earth around its own axis is taken into account by assuming that the worldline of the ground station is a spatially circular timelike curve with constant value of the angular coordinate ϑ, given by the station's latitude π/2 − ϑ.
The effect of the earth's rotation on the spacetime geometry is not taken into account. The model also contains a clock at the ground station which gives the station's Finsler proper time.
Finally, the downlink from the satellite to the ground station is described by null geodesics connecting the world lines of the satellite and the ground station.

Equations and Numerical Treatment of Our Model
The equations of the causal geodesics, proper times, and frequency shifts that are the basis of the simulations of the satellite orbits and the signal propagation presented in Sections 8 and 9 are given in Appendix B. The iterative procedure for determining the null geodesics from the satellites to the ground stations will be outlined in Section 6. The calculations were performed by one of us (I.A.) with the software Mathematica R , version 12.0.0.0, making use of its inbuilt procedures for numerical (differential) equation solving and integration. The numerical effort was reduced by taking advantage of the fact that every causal geodesic r(t) is restricted to the plane spanned by the vectors r(t 0 ) and dr/dt(t 0 ) at arbitrary time t 0 . Therefore, all satellite orbits and signals (null geodesics) were calculated by choosing two orthonormal basis vectors k and l in the plane spanned by the respective initial position and velocity vectors and then solving the two-dimensional analogue of (73). The working precision (maximum number of digits for internal computations) was 105 for the satellite orbits and 65 for all other calculations.

Comparison with the Results of Lämmerzahl, Perlick, and Hasse
In the article by Lämmerzahl, Perlick, and Hasse [7] the effects of φ 0 , φ 1 , and φ 2 on three known phenomena (time delay of light rays, light deflection by the sun, perihelion precession of the planet Mercury) are presented under the assumption that with regard to these effects the perturbation functions can be approximated as (see (39) in [7]) These results were obtained by numerical treatment of integrals rather than by numerical treatment of the geodesic equations. It is therefore an obvious approach, among other tests, to calculate these effects also with the software described in the preceding section and to compare the results.
In general the calculations described in Section 4 confirm the results presented in [7], but we found two discrepancies regarding light deflection and perihelion precession. Therefore, we additionally calculated these two effects on basis of the equations in Appendix A. These calculations were performed again by one of us (I.A.) with Mathematica and independently by another of us (M.P.) with the software MATLAB R , version 5.3.0.10183 (R11). The MATLAB calculation made use of its available solver functions for ordinary differential equations and were performed with the relative tolerance of the solver output set to 10 −12 .

Time Delay of Light Rays
Consider a radio signal that is sent from Earth (r = 1 au), then passes the sun at r = 0.0074 au and finally reaches a spacecraft at r = 8.43 au. According to Equation (82) in [7] the difference between the perturbed one-way time delay δt and the unperturbed one-way time delay δt 0 is The calculations described in Section 4 exactly confirm this result.

Light Deflection
Consider a light ray from a distant star (r → ∞) that grazes the sun at r = 0.0046 au and finally reaches Earth (r = 1 au). According to Equation (73) in [7] the difference between the perturbed deflection angle ∆ϕ and the unperturbed deflection angle ∆ϕ 0 is Regarding the coefficients of φ 01 and φ 11 our calculations confirm Equation (5), with maximum deviation ± 1 in the last reported digit. But for the coefficient of φ 21 the results agree less well. For this coefficient, both Mathematica programs give the value 3.6 × 10 −11 (for e.g., 10 −8 ≤ φ 21 ≤ 10 −2 ). An exemplary calculation with MATLAB gives the value 2.9 × 10 −11 . We cannot explain this discrepancy at the moment.

Perihelion Precession
Equation (107) in [7] for the perihelion motion of Mercury is In this equation ω and ω 0 denote the precession rate ("angular velocity") of the perihelion in the perturbed and unperturbed case, respectively. Whereas our calculations confirm the coefficients of φ 11 and φ 21 , with maximum deviation ± 1 in the last reported digit, for the coefficient of φ 01 we get the value of −0.5, which differs by more than seven powers of ten from the value in Equation (6). We cannot uncover the cause of this discrepancy, but can only guess that it has to do with the orbital period. Namely, since the perihelion motion is calculated with fixed apsidal distances, φ 0 , φ 1 , and φ 2 cause a change ∆T s of the sidereal period, whereby the effect of φ 0 clearly dominates (on circular orbits the effects of φ 1 and φ 2 even vanish, cf. Equation (30) in [7]). ∆T s in turn makes the largest contribution to the change ∆T a of the anomalistic period. Ultimately, ∆T a makes the dominating contribution to ∆ω. This contribution can be calculated analytically in the framework of the PPN formalism, which also results in a value of approximately −0.5 for the coefficient of φ 01 .

The Iterative Procedure for Determining the Downlink Signals
In our model the satellite transmits a signal to the ground station every 2000 s with respect to its proper time τ . Therefore, after the calculation of the satellite orbit r S (t), the next step is the calculation of the coordinate transmit times t n corresponding to the proper transmit times τ n by integration of (50). Then the null geodesics from the satellite positions r S (t n ) to the rotating ground station with position r B (t) are calculated by means of (73). To this aim the main task is the iterative determination of the unknown emission directionsd(t n ). We will explain the procedure for a signal emitted at t 0 with unknown emission directiond(t 0 ). The first approximation d 1 of the emission direction is calculated according to the formula By making use of (48) for the coordinate light speed c at the emission position in the emission direction, the signal r L (t) emitted at the position r S (t 0 ) in the direction d 1 is then calculated, together with its approach distance ∆ := |r L (t) − r B (t)| min from the ground station. Next the gradient of the approach distance with respect to the emission direction is calculated, i.e., and the emission direction is updated according to for a series of steps of increasing width s, until the approach distance stops decreasing at a step width s max . This defines the improved emission direction This strategy is repeated several times and after the calculation of typically about 40 null geodesics the approach distance is typically only about 5 µm. The time t a of closest approach is used as coordinate receive time t r of the signal at the ground station. By integration of (50) the coordinate receive times are converted to the proper receive times τ r at the ground station. Finally, for every signal the receive frequency shift is calculated according to (80).

Hybrid Models
Since the perturbations influence both the orbital motion of the satellite and the signal propagation, the interesting question arises which of the two influences predominates in typical cases. This is not at all clear from the outset, as the following considerations show.
The orbital speeds are about five powers of ten smaller than the speed of light, so the satellites "feel" the φ 1 -and the φ 2 -perturbations (affecting only the spatial geometry) much less than the signals. One would therefore expect that the observable effects of the φ 1 -and the φ 2 -perturbations are quite predominantly caused by perturbations of the signal propagation, while a φ 0 -perturbation should affect the satellite orbit and the signals approximately equally. However, this expectation does not take into account that the perturbations act on the signals only for a short time and that the null geodesics describing the signals and the timelike geodesics describing the satellite orbits must fulfil significantly different initial and boundary conditions. While for the satellite orbit the initial position and the initial velocity are given, for the signals the initial position and the condition to reach the ground station are given. Thus, in a certain sense the geodesic of the signal is "more closely tied" than that of the satellite and therefore less "bent" by the perturbations. In addition, the effects of the perturbations on the signals are never accumulative. The situation is different with orbital perturbations because even if a disturbed orbit is compared with an undisturbed orbit with the same sidereal period, there is an accumulative effect due to the generally different rates of perihelion precession.
Further consideration shows that the so far open question of the predominant effect is relevant with regard to the separation of Finslerian effects from effects of other perturbations. In our class of models, Finslerian (φ 2 -) perturbations do not affect satellites and signals in radial motion and their effects are greatest in a certain angular range between radial and tangential motion. Therefore, if the orbit perturbations dominate, Finslerian effects would be greatest in certain sections of the orbit between the apsides and would increase with the eccentricity of the orbit. Conversely, if the signal perturbations dominate, Finslerian effects would be greatest when the satellite, seen from the ground station, is in a certain middle range of altitude above the horizon.
In order to answer this question we additionally studied "hybrid models", in which only the satellite orbit is and is not, respectively, affected by the perturbations. These models are not intended as a description of nature. Rather, they help to understand the results obtained with the actual model and to identify those orbits, orbit segments, and constellations (relative positions) of the satellite and the ground station that are most promising for the detection of possible genuine Finslerian effects.

Simulation Parameters
We will now study the effects of the first order perturbations φ n (ρ) = φ n1 r S /ρ (n = 0, 1, 2) on the time signals that are transmitted from the satellites Galileo 5 and 6 and received at a ground station on Earth. To that aim we simulated the motion of two model satellites with the following orbital parameters, that are similar to the parameters of the satellites Galileo 5 and 6: inclination i = 50 • , perigee distance ρ p = 23,500 km, and perigee speed v p = 4.436 km/s and thus apogee distance ρ a = 32,462 km, semi-major axis a = 27,981 km, and eccentricity e = 0.1601. As the arguments of perigee ω of the satellites Galileo 5 and 6 (GSAT0201 and GSAT0202) change with a rate of 0.034 • /d [9] the simulations were performed both for ω = 90 • and ω = 180 • .
The signals from the satellites were received at two ground stations that followed the rotation of the earth. For the description of the (spherical) earth we used the parameter values radius ρ E = 6371 km, angular speed ω E = dϕ/dt = 7.292115 × 10 −5 rad/s, standard gravitational parameter GM E = 398,600.44 km 3 /s 2 , and Schwarzschild radius r S = 2GM E /c 2 . The first ground station was on the equator and received the signals from the satellite with argument of perigee ω = 90 • . A second ground station was at a latitude of 45 • north and received the signals from the satellite with argument of perigee ω = 180 • . The detailed simulation parameters are listed in Table 1.
As we will find out, each perturbation causes characteristic patterns in the received signals, independent of the argument of perigee of the satellite and the latitude of the ground station.
To check the reliability of the results, the simulations were performed with two different values of each parameter φ n1 . We will see, that all results linearly depend on the parameters φ n1 in the time range under study, to a very good approximation. This confirms both the assumed smallness of the perturbations and the reliability of the results. Table 1: Simulation parameters: argument of perigee ω and orbital state vectors r and v of satellite S and ground station B at initial time t 0 .x,ŷ, andẑ denote the unit vectors along the coordinate axes.

Sidereal Period
In the following sections we will compare the received signals at the ground stations in the unperturbed spacetime and in the perturbed spacetimes. For the underlying satellite orbits we require the same value of the "area coordinate" r at perigee (same value of ρ p ) and the same sidereal period, measured with respect to the proper time of a stationary clock far away from the earth. The proper time of such a clock can be sufficiently approximated by the coordinate time t. With respect to the coordinate time the sidereal period in the unperturbed spacetime is T 0 = 46,580.92981316189761 s. Keeping the orbital state vectors of the satellites at t 0 fixed according to Table 1, the perturbation functions φ n result in shifts ∆T = T − T 0 of the sidereal period that are listed in Table 2.
Therefore, to fulfil the constraint of same sidereal periods in the unperturbed spacetime and in the perturbed spacetimes, the simulations were performed with a slightly reduced satellite speed at t 0 in the perturbed spacetimes according to Table 3. Table 3: Residual shifts ∆T of the sidereal periods of the model satellites, divided by T 0 , with reduced satellite speeds v p + ∆v p at t 0 .

Satellite Orbits
In this section we report about the effects of the perturbation functions on the satellite orbits.
To start with, Figure 1 shows the coordinate distance between the satellite and the ground station in the unperturbed spacetime as function of the coordinate time.
In comparison, Figures 2-4 show the shifts of the satellite positions due to the perturbation functions φ n , divided by the magnitudes φ n1 of the perturbations. In combination with the values of the parameters φ n1 listed in Table 3 it can be seen from these figures, that in the depicted time span the maximum shifts of the satellite positions caused by the perturbations are 20 cm (φ 01 = 10 −8 ) and 7 cm (φ 11 = φ 21 = 1), respectively.     The root of the very similar effects of the two perturbations can be read from (73) and (54). According to these equations, the differences between the accelerations caused by the perturbations φ 11 = 1 and φ 21 = 1 depend on the factors f = (r · v) 2 /(r · r v · v) and f 2 . But along the satellite orbits defined in Table 1, the angle between r and v is always between 80.8 • and 99.2 • . Therefore, the peak values of f and f 2 are only 2,6% and 0,066%, respectively. Only for orbits of higher eccentricity these factors become significant.
As our model is based on a point-like Earth, the signals reach the ground station both when the satellite is above the horizon and when it is below the horizon. The latter is not possible in reality. In the figures below, the regions in which the satellite is under the horizon are therefore shaded, since comparison with observational data is not possible there.

Receive Times
We assume that the satellite transmits signals at a constant rate with respect to its proper time τ t . The zero points of the time axes are fixed by the assumption that one signal is transmitted at coordinate time t 0 , corresponding to τ t = 0, and received at the ground station at its proper time τ r = 0. This applies to both the unperturbed and the perturbed scenarios. With the time axes fixed that way, the shift of the receive times τ r caused by the perturbations is calculated according to ∆τ r := τ r (φ n = 0) − τ r (φ n = 0).
The results are shown in Figures 5-7. We see that φ 0 , φ 1 , and φ 2 affect the receive times qualitatively different. φ 0 causes a simple and almost periodic pattern with ∆τ r < 0 for the most part, whereas φ 1 and φ 2 cause rich and irregular patterns with varying amplitudes.
Accordingly, φ 0 -perturbations can easily be distinguished from φ 1 -and φ 2 -perturbations by analysing the shifts of the receive times. But it seems that φ 1 -and φ 2 -perturbations cannot be easily separated this way, at least not in the case of the satellites Galileo 5 and 6.   To get more information and to find out how to possibly achieve complete separation, we calculated the receive time shifts in two hybrid models with the same orbit parameters as above, as introduced in Section 7. In the first model the perturbations affect only the satellite orbit. In the second model the perturbations affect everything but the satellite orbit, i.e., the signal propagation and the proper times of the satellite and the ground station. Figures 8-10 show the results.   In the case of a pure φ 0 -perturbation the receive time shift caused by the perturbation of the orbit ("orbit contribution") completely dominates the receive time shift caused by the perturbation of the signals and the proper times ("signal contribution"), as can be seen from Figure 8. In contrast, in the case of a pure φ 1 -perturbation both contributions are of about the same size, see Figure 9. In addition, there is a characteristic pattern: sharp positive peaks in the signal contribution tend to coincide with negative extrema in the orbit contribution. Thus, Figure 9 reveals the hidden compensating effects in Figure 6. Finally, according to Figure 10, in the case of a pure φ 2 -perturbation the orbit contribution dominates already during the first revolutions of the satellite, although less than with a φ 0 -perturbation. However, the receive time shift (red dashes in Figure 7) is well approximated by the orbit contribution (black dashes in Figure 10). It is also notable that the signal contributions of φ 1 and φ 2 tend to have the same sign in the shaded regions and to have opposite signs in the middles of the unshaded regions.
The following considerations make plausible that φ 0 , φ 1 , and φ 2 give qualitatively different results.
Under the constraint of same sidereal period in the unperturbed and in the perturbed models (cf. Section 8.2), a φ 0 -perturbation results in a considerable change of the satellite orbit's semi-major axis (this follows from a simple consideration similar to that in the last paragraph of Section 5). That change is in turn associated with a considerable shift of the receive times and receive frequencies. In contrast, φ 1 and φ 2 affect the semi-major axis and hence the orbital motion of the satellite only slightly (cf. Figures 2-4 for the effect of φ 0 , φ 1 , and φ 2 on the satellite orbit).
The relatively sharp and strong peaks in the receive time shift that occur only with φ 1 and only in the shaded regions (one peak every time the satellite is near the nadir) can be understood as the results of fictitious experiments of light deflection at a point mass, and it is already known from Section 5.2 that the effect is much stronger with φ 1 than with φ 2 .
As mentioned at the end of Section 2, φ 2 does not affect pure radial velocities. This makes plausible that the signal contributions of φ 1 and φ 2 (red lines in Figures 9 and 10) differ significantly also in the middles of the unshaded regions, i.e., where the satellite is near its maximum height above the horizon seen from the ground station and the signal velocity has the largest radial component.
In summary, Figures 9 and 10 demonstrate that in the case of the satellites Galileo 5 and 6 the signal contributions of φ 1 and φ 2 differ more than the orbit contributions. But in the solely accessible unshaded regions the orbit contributions dominate and hamper the reliable separation. However, in Section 9 we will demonstrate that for highly eccentric orbits the separation is possible and in this case the different effect of φ 1 and φ 2 on almost radial signals is crucial. Therefore, one has to conclude that the possibility of the separation of the receive time shifts caused by φ 1 and φ 2 depends on the orbit of the satellite.

Receive Frequencies
Next we consider the receive frequency shift at the ground station caused by φ 0 , φ 1 , and φ 2 . For comparison, Figure 11 shows the relative shift between transmit frequency f t and receive frequency f r in the Schwarzschild spacetime. In the Schwarzschild spacetime the frequency shift is dominated by the Doppler effect resulting from the relative motion of the satellite and the ground station. To analyse the effects of the perturbation functions on the frequency shift we use the equation (13) Since ∆f r /f t as defined by (13) is essentially given by −d∆τ r /dτ t , analogous patterns as with ∆τ r are to be expected. Indeed, the graphical representation of our results in Figures 12-14 confirms this expectation. Δf r / f t / ϕ 11 Figure 13: Relative receive frequency shift ∆f r /f t divided by φ 11 . Black dashes φ 11 = 0.1, red dashes φ 11 = 1, left (right) curves ω = 90 • (180 • ).
As with ∆τ r , we see that φ 0 , φ 1 , and φ 2 affect the received signals qualitatively different. Again φ 0 causes a simple, almost periodic pattern with almost constant amplitude, whereas φ 1 and φ 2 cause patterns with varying amplitudes, and again only φ 1 causes sharp peaks in ∆f r /f t which are located in the shaded regions of Figure 13. However, unlike with ∆τ r , these peaks come in pairs oriented downward and upward, as is expected from the explanation following (13). Figures 15-17 show the corresponding results of our hybrid model calculations. As with ∆τ r , we see that only φ 1 perturbs the signal propagation to an extend that contributes significantly to the whole effect (and again only in the shaded regions).
Finally we observe that φ 2 causes a clearly increasing amplitude of ∆f r /f t (see Figure 14) and that φ 1 should also cause this effect. Namely, the black dashes in Figure 16 indicate that the "orbit contribution" of ∆f r /f t also increases in time and will eventually dominate the non-increasing (or at least much less increasing) "signal contribution". Simulations of ∆τ r and ∆f r /f t over time spans of up to 10 7 seconds confirm this expectation (and show no new features).
As expected, for the frequency shifts the situation regarding the desired separation of the effects of φ 1 and φ 2 is very similar to that for the receive time shifts.
In summary we have to state that at least in the examples considered here, based on the Galileo 5 and 6 satellites, the effects of a genuine Finslerian (φ 2 -) perturbation and of a φ 1 -perturbation can hardly be separated. This assessment takes into account that noise and additional gravitational or non-gravitational perturbations may cover or distort the specific structures in the signals.

High-Eccentricity Orbit Simulation
This section deals with a model satellite in a prograde equatorial Earth orbit with perigee distance ρ p = 23,500 km and perigee speed v p = 5.370 km/s and thus apogee distance ρ a = 133,227 km, semi-major axis a = 78,364 km, and high eccentricity e = 0.7001. The ground station is on the equator and we assume ϕ station = ϕ satellite + π/2 at t = 0.
This choice of parameter values is justified as follows. On the one hand, the perigee distance should not be significantly smaller than in the case of the Galileo 5 and 6 satellites since the satellite should not enter the inner Van Allen radiation belt. On the other hand, our analytical calculations (these are too extensive to be presented here) show that at an angle α (describing the direction of the satellite's velocity and defined by Equation (23)) of about 0.86 the absolute value of the α-dependent factor of the φ 2 -part ofφ is maximal, see Equation (29) (whereas φ 1 does not enter the expression forφ). It turns out that for α to be about 0.86 somewhere on the trajectory, the eccentricity must be at least about 0.7. This explains the relatively great apogee distance. (When comparing the following figures with the respective figures in the previous section please note that now the semi-major axis is larger than the radius of the geostationary orbit. Seen from the ground station the high eccentricity satellite is thus moving from east to west whereas the Galileo 5 and 6 satellites are moving from west to east).
As in Section 8 the initial values of the perigee speed are so adjusted that the sidereal period is the same with and without perturbation (absolute value of relative shift <10 −17 ).
Since with regard to φ 0 there are no substantially different patterns than in Section 8, we will discuss here only the disentanglement of the effects of φ 1 and φ 2 . Figures 18 and 19 show that the effects of φ 1 and φ 2 on the receive times are significantly different, even in the unshaded regions (please note the different scaling of the ordinate axes). Especially, unlike with the Galileo satellites, there are periods of time in which φ 1 and φ 2 shift the receive times in opposite directions. Remarkably, in Figures 20 and 21 the "orbit" and the "signal" contributions to the shifts calculated with the hybrid models are of the same order of magnitude (whereas with the Galileo satellites the orbit contributions clearly dominate, except in the shaded regions in the case of a φ 1 -perturbation). It is also noteworthy that the mentioned shifts in opposite directions present in some periods of time are solely caused by the signal contributions because again, as with the Galileo satellites, the orbit contributions of φ 1 and φ 2 to the shifts run in parallel (but unlike with these satellites the orbit contribution of φ 2 amounts to only about half the orbit contribution of φ 1 ).
As mentioned in Section 8.5, ∆f r /f t as defined by (13) is essentially given by −d∆τ r /dτ t . Accordingly, we expect φ 1 and φ 2 to cause distinguishable structures also in ∆f r /f t .

Conclusions
While purely analytical methods quickly reach their limits in the treatment of the problems considered in this paper, as can be seen from Section IV of the recent article of Hasse and Perlick [8], it is demonstrated here that valuable information about Finslerian effects can be obtained by means of numerical methods, even with today's standard computing systems. Assuming first order perturbations φ n (ρ) = φ n1 r S /ρ (n = 0, 1, 2) in the considered class of spherical symmetric Finsler spacetimes, our numerical simulations in combination with our hybrid models allow us to draw conclusions about the optimal types of orbits for testing models of Finsler gravity. As can be seen in Section 8, the most difficult task is to separate the effects of φ 1 and φ 2 . In the case of the Galileo 5 and 6 satellites the "orbit contributions" to the receive time and frequency shifts caused by φ 1 and φ 2 are practically indistinguishable and dominate the "signal contributions" when the satellites are above the horizon, even within the first revolution of the satellite. The "signal contribution" is significant only when the satellites are below the horizon (i.e., in situations of fictitious experiments of light deflection at a point mass) and only with φ 1 .
However, for highly eccentric orbits genuine Finslerian effects can be separated from effects of perturbations of the Schwarzschild spacetime within the Lorentzian geometry. This is one of the main results and important for the planning of possible future dedicated missions. Although in the illustrating example of Section 9 the "orbit contributions" to the receive time and frequency shifts caused by φ 1 and φ 2 run in parallel too (but with different heights, unlike with the Galileo satellites), the superimposed "signal contributions" differ significantly even when the satellite is above the horizon and make both effects separable.
Unfortunately, the orbit eccentricity of the Galileo 5 and 6 satellites is too small to allow this separation (see our conclusions at the end of Sections 8.4 and 8.5). Nevertheless, these two already available and appropriately equipped satellites in principle allow the determination of the sum of φ 11 and φ 21 as a valuable first step. Although, as stated above, along highly eccentric orbits the contributions of φ 1 and φ 2 to the receive time and frequency shifts differ to an extent that φ 11 and φ 21 can be separately determined, it helps with the analysis of such an experiment if the sum of φ 11 and φ 21 is additionally known from another experiment (the Galileo satellites).
We like to conclude with a few remarks that go beyond the work presented here. In principle, non-gravitational perturbations (see our remarks in Section 1) may be present that are unknown or insufficiently understood and therefore cannot be modelled. Nevertheless, one can assume that such perturbations cause no "Finsler-like" effects in the receive times and frequencies of the signals. The reason is that the various (non-Finslerian) perturbations typically show either a nearly diurnal or semidiurnal variation (if they are not of terrestrial origin, e.g., if they are caused by the sun or the moon) or a variation with a leading period that corresponds to the satellite's orbital period (if they are of terrestrial origin). In the latter case, it is not expected that minima of the non-Finslerian perturbations occur just in the apsides, especially in the perigee, unlike with Finslerian (φ 2 -) perturbations, regarded as additional "forces", that are very weak in the apsides and strongest somewhere between the apsides.
In view of tests of Finsler gravity, an obvious idea is to extend our approach, in particular the method for determining the signals between satellite and ground station and the concept of hybrid models, to the analysis of other situations of communication. So one could first think of a two-way time transfer procedure, i.e., signals that travel both ways between two clocks that are being compared.
Finally, the methods presented here may also be adapted to similar tests of other alternatives to general relativity. In principle, this should be possible at least for theories with well-defined second-order equations of motion for spinless test bodies and light signals in which one can consider a parametrized class of models that includes the Schwarzschild spacetime or other suitable reference solutions of Einstein's field equations. Examples are the BransDicke theory, Einstein-Cartan models of gravity, and Weyl manifolds as spacetime models. In this section we briefly describe the derivation of the second order equations of motion in Schwarzschild coordinates in the plane ϑ = π/2, with the coordinate time as the independent variable. The first steps are the elimination of the affine parameter by means of the constants of motion L, E = −∂L/∂ṫ, and L = ∂L/∂φ (cf. Equations (16)-(19) in [7]) and the linearisation of the resulting equations with respect to the perturbation functions. As we are only interested in causal geodesics, we will only consider geodesics withṫ = 0. The elimination of the affine parameter is based on the constant of motion and the resulting identity Restricted to the plane ϑ = π/2 the Lagrangian (1) reads and the affine parameter can be eliminated by multiplying the right side of (16) with the square of the left side of (15). This leads to Restricted to the plane ϑ = π/2 the constant of motion L = ∂L/∂φ reads and the affine parameter can be eliminated by multiplying the right side of (19) with the left side of (15). This leads to Linearising (20) with respect to the perturbation functions (especially, replacing (1 + φ 0 ) −1 with its linear approximation 1 − φ 0 ) then results in Since there are no derivatives with respect to the affine parameter in (18) and (21), from now on we will designate derivatives with respect to t by an overdot. Derivatives with respect to r will be designated by a prime. Furthermore, with this we define the "speed" u := h rrṙ 2 + r 2φ2 , the angle α uniquely determined by cos α = h rrṙ u and sin α = rφ u (23) and the coefficients p ik := cos i α sin k α (i, k ∈ N 0 ).
By differentiation of (25) and (26) with respect to t we get rid of the constants of motion and obtain the linear (in the algebraic sense) system of equations B Equations for the Three-Dimensional Simulations

B.1 Lagrangian in Isotropic Coordinates
For the numerical study of observable effects of Finslerian perturbations on earth-bound satellites in the general geometry, without a common orbital plane of the satellite and the ground station, we used isotropic coordinates t, x, y, and z. These coordinates treat all directions in space and all orbital planes on an equal footing and are routinely used in orbit calculation programs. Schwarzschild coordinates and isotropic coordinates are related by the equations x = ρ sin ϑ cos ϕ, y = ρ sin ϑ sin ϕ, z = ρ cos ϑ, r = ρ 1 + r S 4ρ With the above results we are now prepared to transform the Lagrangian (1), restricted to the plane ϑ = π/2, to isotropic coordinates: By use of the definitions f := (r ·ṙ) 2 (r · r)(ṙ ·ṙ) , the last line of (41) can be written in the more compact form By use of Lagrange's identity (a × b) · (c × d) = (a · c)(b · d) − (b · c)(a · d) the last line of (41) can also be written as 2L = h ttṫ 2 (1 + φ 0 ) + h ρρṙ ·ṙ 1 + f φ 1 + f (r ×ṙ) · (r ×ṙ) (r · r)(ṙ ·ṙ) φ 2 .
Equations (44) and (45) have been derived for geodesics in the plane ϑ = π/2, but due to their spherical symmetry they are invariant under spatial rotations of the isotropic coordinate system and therefore valid for arbitrary geodesics.
To determine the equation for the coordinate light speed c(r, d) at the position r in the direction d, we first note for the Lagrangian (44) and geodesics withṫ = 0 (e.g., causal geodesics) the identity L(t, r,ṫ,ṙ) =ṫ 2 L(t, r, 1, dr dt ).
we get