{'fraction_non_alphanumeric': 0.04688081809842974, 'fraction_numerical': 0.024148798009245385, 'mean_word_length': 4.327044025157233, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 25, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Light propagation in the gravitational field of N arbitrarily moving bodies in the 1.5PN approximation for high-precision astrometry

High-precision astrometry on sub-micro-arcsecond level in angular resolution requires accurate determination of the trajectory of a light-signal from the celestial light source through the gravitational field of the Solar system toward the observer. In this investigation the light trajectory in the gravitational field of N moving bodies is determined in the 1.5 post-Newtonian approximation. In the approach presented two specific issues of particular importance are accounted for: (1) According to the recommendations of International Astronomical Union, the metric of the Solar system is expressed in terms of intrinsic mass-multipoles and intrinsic spin-multipoles of the massive bodies, allowing for arbitrary shape, inner structure and rotational motion of the massive bodies of the Solar system. (2) The Solar system bodies move along arbitrary worldlines which can later be specified by Solar system ephemeris. The presented analytical solution for light trajectory is a primary requirement for extremely high-precision astrometry on sub-micro-arcsecond level of accuracy and associated massive computations in astrometric data reduction. An estimation of the numerical magnitude for time delay and light deflection of the leading multipoles is given.


I. INTRODUCTION
A substantial advancement in astrometric measurements has been achieved by the astrometry mission Hipparcos (launch: 8 August 1989) of European Space Agency (ESA), which has reached an accuracy of a milliarcsecond (mas) in determining the angular positions of about 10 5 stars [1,2]. Meanwhile, the state-of-the-art angular observations have finally arrived at the level of a few micro-acrseconds (µas) [3,4]. Especially, the stunning progress in astrometry has proceeded with the ESA mission Gaia [5] (launch: 19 December 2013) which aims at an all-sky survey of more than 10 9 stars of our galaxy and targets angular accuracy of up to a few µas for bright stars in the final catalog scheduled for publication in 2022.
In view of these advancements it becomes obvious that future astrometry is going to force into the exciting areas of sub-µas or even nano-arcsecond (nas) level of accuracy. To step up efforts toward sub-µas-astrometry is of fundamental importance in astronomy and astrophysics. For example, an accuracy of about 10 nas in angular resolution would allow for direct measurement of trigonometric parallaxes of stars belonging to galaxies of the Local Group which spans a diameter of about 10 7 light-years, that means would enable to determine spatial distances of extra-galactic objects independently of dynamical models of the Universe. Moreover, also extremely high-precision tests of relativity, detection of dark-matter distributions within or outside of our galaxy, determination of stellar and galactic kinematics, and finally even the discovery of one-Earth-mass exoplanets in the habitable zone of nearby Sun-like stars would be possible by means of subµas-astrometry. Recently, there are several mission proposals in this respect. For instance, the mission NEAT [6,7] has been proposed to ESA which intends to reach a precision of about 50 nas in angular resolution for being able to detect Earth-like exoplanets surrounding stars in the stellar neighborhood of the Sun. Further space missions like ASTROD [8,9], LATOR [10,11], ODYSSEY [12], SAGAS [13], or TIPO [14] have been proposed to ESA which imply the determination of light trajectory through the Solar system on sub-µas or even at nas level of accuracy. Also earth-bound telescopes are under consideration which aim at angular resolutions of about 10 nas [15].
But, although in view of the recent impressive achievements, the step from µas-astrometry toward sub-µaslevel or even nas-level of accuracy in angular resolution will surely be a long-term goal in the astronomical science. This is because the envisaged advancement toward space-based nas-astrometry implies many subtle effects and new kind of challenges which have not been encountered before: What kind of optical technology would allow for nas-astrometry? Is it technologically possible to measure the velocity of spacecraft (observer) with sufficient accuracy allowing for a precise determination of aberrational effects? How accurate do we have to determine the ephemeris of the Solar system bodies and could such precise ephemeris be provided? Is it possible to model accurate enough the influence of interstellar medium on light propagation? How strong is the effect of gravitational waves on light propagation on nanoarcsecond level? How is it possible to account for the gravitational light deflection caused by massive bodies located outside the Solar system?
Each of these and many other problems have to be clarified before nas-astrometry becomes feasible. But certainly, the fundamental assignment in astrometry remains to trace a lightray observed in the Solar system back to the celestial light source. The importance of this fact has also been underlined recently by the ESA-Senior-Survey-Committee (SSC) in response of the selection of science themes for the L2 and L3 launch opportunities, where it has been stated that "SSC recommends that proper modeling tools, most notably the availability of a General Relativistic framework able to model photon trajectories to the accuracy required should be given the proper attention to prove feasibility" of high-precision astrometry [16]. According to this, the primary effort in any astrometrical framework concerns the precise description of the light trajectory, that is to say the determination of the spatial coordinates of a light-signal as function of coordinate time, x (t), in some global coordinate system. Accordingly, the principal purpose of this investigation is the determination of the trajectory of a light-signal propagating through the Solar system. In the following four subsections it will be enlightened how one has to proceed in order to arrive that goal: (A) the theory of light propagation, (B) the post-Newtonian expansion, (C) the state-of-the-art, and (D) the primary objective of this investigation.

A. Theory of light propagation
The determination of spatial coordinates of the lightray takes the most simple form in the flat Minkowskian space-time and assuming a Cartesian coordinate system which covers the entire space, implying the metric tensor η αβ = diag (−1, +1, +1, +1). Suppose the light-signal is emitted at some initial time t 0 by a light source located at some space-point x 0 , then the light trajectory is simply given by a straight line which is also called unperturbed light trajectory, where the unit-vector σ determines the direction of lightpropagation and the sublabel "N" denotes Newtonian approximation.
In general relativity the four-dimensional space-time in the presence of matter is curved, that means is described by a semi-Riemannian manifold with non-vanishing curvature tensor rather than a flat Minkowskian space-time, and a light trajectory is no longer a straight line but propagates along a so-called null geodesic, which generalizes the concept of a straight light trajectory. The four-coordinates x α (λ) of a light trajectory depend on some affine curve-parameter λ, and are determined by the geodesic equation [17,18], where (2a) represents the geodesic equation, while the isotropic condition (2b) is an additional constraint for a null geodesic, a term which refers to the fact that the invariant line element vanishes, ds 2 = dx α (λ) dx α (λ) = 0, at any point along the light trajectory. The Christoffel symbols in (2a) are related to the metric of curved spacetime as follows: where g αβ and g αβ are the contravariant and covariant components of the metric tensor, respectively, where the metric signature (−, +, +, +). The geodesic equation (2a) represents a second-order differential equation, hence an unique solution implies the need of two initial values for the lightray: The equations in (2a) -(2b) are valid in any reference system. But in practical astrometry one is necessarily enforced to specify the reference systems for concrete observational data. In line with the recommendations of the International Astronomical Union (IAU) [19,20], the Barycentric Celestial Reference System (BCRS) with coordinates (ct, x) is the standard global chart to be used in modern-day astrometry, where t is the BCRS coordinatetime and x are Cartesian-like spatial coordinates from the barycenter of the Solar system to some field-point. Consequently, it becomes much preferable to exploit the freedom in the choice of scalar curve-parameter λ and to rewrite the affinely parametrized geodesic equation (2a) and the isotropic condition (2b) in terms of BCRS coordinate-time [17,18,21]: The zeroth component in (6a) does not carry any new information because it vanishes identically. In order to determine the solution of (6a) it is advantageous to transform the initial conditions in (4) -(5) into initialboundary conditions [18]: with (7) being the position of the light source at the moment t 0 of emission of the light-signal and (8) being the unit-direction of the lightray at past-null infinity. Then, the exact solution of (6a) for the light trajectory from the light source through the Solar system toward the observer can be written as follows, where the term ∆x (t, t 0 ) denotes gravitational corrections to the unperturbed light trajectory (1).

B. Post-Newtonian expansion
The correction terms ∆x (t, t 0 ) in Eq. (9) are highly complicated expressions which cannot be determined exactly and one has to resort on approximation schemes. Such an approximation scheme is provided by the post-Newtonian expansion of the metric of Solar system, which represents an expansion in terms of inverse power of the speed of light, up to terms of the order O c −5 given by: where h (n) αβ = O (c −n ) with n = 2, 3, 4. The justification of such an expansion is based on the fact that the gravitational fields in the Solar system are weak, (G M A ) / c 2 P A 1, as well as the velocities of the Solar system bodies are slow, v A /c 1, where M A , P A , and v A means mass, radius, and velocity, respectively, of some massive body A. For these reasons the post-Newtonian expansion is also called weak-field slowmotion expansion. As outlined in [17,[21][22][23], such an expansion is valid inside the near-zone of the Solar system, |x| λ gr , where λ gr ∼ 10 17 meter is a characteristic wavelength of gravitational waves emitted by the Solar system. The near-zone of Solar system is so large that it still contains all Solar system bodies and even encompasses the nearest stars of the stellar neighborhood of the Sun.
Inserting the expansion (10) into (6a) yields the geodesic equations for lightrays up to terms of the order O c −5 . Accordingly, the expansion of the metric in (10) inherits a corresponding expansion of the lightray, that means the corrections to the unperturbed lightray can formally be written as follows: (11) where ∆x 1PN = O c −2 are 1PN corrections, ∆x 1.5PN = O c −3 are 1.5PN corrections, and ∆x 2PN = O c −4 are 2PN corrections to the unperturbed lightray. In view of the fact that the post-Newtonian expansion of the metric (10) is only valid within the near-zone of the Solar system, the post-Newtonian expansion of the lightray (11) allows for near-zone astrometry, in particular for reduction of astrometric observations of all Solar system objects. The unique interpretation of astrometrical data of far objects, like stars or quasars, is the subject of farzone astrometry and necessitates the determination of light trajectory outside the near-zone of the Solar system. That especially means, the light trajectory in the near-zone has to be aligned with the light trajectory in the far-zone by means of a so-called matching procedure as described in detail in [21,24,25] which, however, will not be a topic of this investigation.

C. State-of-the-art in the theory of light propagation
A brief survey about the present status in the theory of light propagation in the gravitational field of massive bodies has recently been presented [26]. Here we will summarize and update that survey. In particular, we will restrict our review on those investigations where the explicit time-dependence of the photon's spatial coordinate, x (t), has been determined, a prerequisite for interpreting real astrometrical observations.

Monopoles at rest
The case of light propagation in the Schwarzschild metric, i.e. in the gravitational field of one spherically symmetric massive body at rest, where x A = const is the constant position of the body, is the most simple case and has been determined long time ago in 1PN approximation, e.g. [18,24,[26][27][28].
The solution for the light trajectory is given by Eq. (J7). Besides its simplicity, the determination of the photon's spatial coordinate in the Schwarzschild-field is the initial point in the theory of light propagation in astrometry.

Monopoles in motion
In reality, the bodies A = 1, ..., N of the Solar system move along their time-like worldlines x A (t) and for todays extremely high-precision in astrometric measurements the gravitational field of some Solar system body can not any longer be treated as static and spherically symmetric. In a first approximation, the motion of one massive body A can be considered as translational motion with constant velocity v A : where x A = x A (t A ) and v A = v A (t A ) are the spatial position and velocity of the body A at some initial time-moment t A . The light-trajectory in the field of one massive body in translational motion has completely been solved in 1PN approximation in [29]. This solution has later been rederived by means of a suitable Lorentz transformation [34]. Following a suggestion in [30], in the investigation [24] it has been shown that the free parameter t A in Eq. (13) should be chosen as the time-moment of closest approach (given by Eqs. (B14) and (B15)) between the massive body A and the photon in order to minimize the residual effects caused by the approximation of the real motion by a translational motion of the massive body. With the aid of advanced integration methods, originally introduced in [31] and further developed in [32], a rigorous solution for the trajectory of a light-signal through the gravitational field of an arbitrarily moving body has thoroughly been solved in [33] in the first post-Minkowskian approximation. The first post-Minkowskian approximation takes into account all terms proportional to the gravitational constant and especially all terms to any power in v A /c, hence the body can even be in ultra-relativistic motion and, therefore, the post-Minkowskian approximation is often called weak-field approximation opposite to the post-Newtonian approximation which is called weak-field slow-motion approximation. Comparing the solution in [33] with [29,34], it has been demonstrated in [35] that the simpler solution for the light-trajectory in the field of a uniformly moving body is actually sufficient for high-precision astrometry on sub-µas-level provided the free parameter t A is chosen either as time-moment of closest approach or as retarded time-moment (given by Eq. (143)) between the photon and the position of the massive body. All these results agree with our investigation in [26] for the case of bodies in slow-motion.

Spin-dipoles at rest
The light trajectory in the gravitational field of one body at rest having spin-dipole S A = const has first been solved in [36] and later confirmed in [24]. The magnitude of light deflection due to the rotational motion of Solar system bodies has been determined in [28,36] and turns out to be significant for astrometry on sub-µas-level of accuracy.

Spin-dipoles in motion
In [37] an explicit solution for the light-trajectory in the field of N uniformly moving bodies with intrinsic spin has been obtained. A comprehensive solution in 1PM approximation for the light-trajectory in the field of N arbitrarily moving bodies with individual spin-structure has been derived in [38] using the already mentioned advanced integration methods originally developed in [31,32].

Mass-quadrupoles at rest
The solution for the light-trajectory in the field of mass-quadrupoles at rest in 1PN approximation was given in [36] and later in [24,26,28]. Especially, in [36] the magnitude of light deflection caused by the massquadrupole structure of Solar system bodies has been determined, where it was figured out that astrometry on µas-level of accuracy is able to detect this light deflection effect. In fact, the light deflection due to the quadrupolestructure of Jupiter is presently under investigation by the ESA astrometry-mission Gaia [5].

Mass-quadrupoles in motion
The light trajectory in the field of N arbitrarily slowlymoving bodies with quadrupole structure has been determined in [26]. Recently, the light-trajectory in the field of N uniformly moving bodies with mass-quadrupole structure has also been obtained in [37] by integrating the geodesic equations for the lightray. Another interesting approach has been found in [42], which is based on the Time Transfer Function (TTF) which avoids to solve the geodesic equations and hence circumvents some of its involved peculiarities.

Higher mass-multipoles and spin-multipoles at rest
A fruitful and systematic approach which allows to integrate analytically the geodesic equations in 1.5 approximation in the field of one body at rest having full timeindependent mass-multipoles M A L and spin-multipoles S A L to any order in the multi-index L has been introduced in [31].
The advanced integration method in [31] has been developed further in [32] for the case of time-dependent mass-multipoles M A L (t) and spin-multipoles S A L (t) in 1PM approximation. Using this advanced approach the analytical solution in 1PM approximation for the lighttrajectory in the field of one massive body at rest with the full set of time-dependent multipoles has been determined in [39,40]. One comment should be in order at this stage. Namely, it is of course possible to interpret the Solar system just as one global massive body A which consists of many individual small massive bodies. But then the solution in [39,40] has to be interpreted as still expressed in terms of global multipoles m L (t) and s L (t) which characterize the entire multipole structure of the Solar system as a whole. However, physically meaningful multipoles can only be defined in the local reference system of each individual massive body. This important issue will later be further considered in some more detail.
Another approach is based on the solution for the TTF and its spatial derivative. A corresponding multipole decomposition of the TTF has been applied in [41] in order to determine the coordinate travel time and the light deflection of a lightray in the gravitational field of one axisymmetric body at rest expressed in terms of massmultipoles M A L .

Higher mass-multipoles in uniform motion
In [42] the TTF approach in 1PM approximation has been applied for the case of light propagation in the field of one axisymmetric body in uniform motion. Especially, an expression for the TTF and its spatial derivative is obtained for this case, which allows to determine astrometric observables like the coordinate travel time of the lightray, the direction of an incident lightray, and the gauge-invariant angle between the direction of two incoming photons. A similar investigation has been done in [43], where the TTF approach has been used in order to determine the coordinate travel time of a lightray in the field of one slowly and uniformly moving extended body with full mass-multipole and spin-multipole structure. 9. 2PN light propagation in the field of monopoles Light propagation 2PN approximation is not on the scope of the presented investigation, but for reasons of completeness some results obtained in 2PN approximation will briefly be mentioned, not only because of its relevance for future high-precision astrometry on sub-µaslevel of accuracy but also for its importance in todays high-precision astrometry on µas-level.
An important progress has been made in [18,27], where an analytical solution of the light-trajectory in 2PN approximation has been determined with explicit time-dependence of the photon's spatial coordinates by solving the null geodesic equations. This solution has later been confirmed by several progressing and ongoing investigations [24,37,[44][45][46], and has also been determined in this investigation, see Eqs. (J8) -(J10). Furthermore, in [46] the time-derivative of the light trajectory in the field of two pointlike bodies at rest has been obtained, allowing to determine the light deflection in such a system. An important new result of this investigation is the fact that the 2PN two-body effect in the Solar system is less than 0.1 nas which considerably simplifies future analytical investigations for high-precision astrometry on sub-µas-level of accuracy.
In [47][48][49] the general formalism of how to determine the TTF and its derivatives has been extended up to the second post-Newtonian (2PN) and second post-Minkowskian (2PM) order, that means including all terms to order O G 2 . The formalism has finally been specified for the case of light propagation in the gravitational field of one spherically symmetric body at rest where the 2PM and 2PN approximations become identically. Especially, explicit expressions for the coordinate travel time of lightray, for the direction of the lightray, and for the angular separation between two incident lightrays have been obtained.
Finally, we also mention another approach which is based on the eikonal concept [50], where the light trajectory in 2PN approximation in the field of one spherically symmetric body at rest has also been derived. The results of this work completely agree with [45].

D. Primary objective of this investigation
According to the survey given above about the present situation in the theory of light propagation, thus far there is no analytical solution available for the light trajectory in the field of arbitrarily moving extended bodies in 1.5PN approximation which, however, is of decisive importance in future high-precision astrometric measurements on sub-µas-level of accuracy and its foreseen involved massive computations, see also [26]. In respect thereof, two important aspects must carefully be treated: (1) The metric perturbations in the exterior of the massive bodies can be decomposed in terms of global mass-multipoles m L and global spin-multipoles s L [51][52][53][54]: These global mass and spin multipoles describe the gravitational field of the Solar system as a whole. However, from the theory of relativistic reference systems it is clear that physically meaningful multipole moments of a massive body A have to be defined in the body's local reference system (cT A , X A ) tied to that body under consideration. Such multipoles are called intrinsic massmultipoles M A L and intrinsic spin-multipoles S A L . Then the question arises about how to express the global BCRS metric in terms of such intrinsic multipoles, that is to say how to determine the global metric perturbations: Such a framework has been elaborated by the approach of Damour-Soffel-Xu (DSX) [55][56][57][58] and within the Brumberg-Kopeikin (BK) formalisms [18,[59][60][61][62], both of which became a part of the IAU resolutions [19,20].
(2) The second issue concerns the motion of the massive Solar system bodies. While in first approximation these bodies orbit the barycenter of the Solar system along ellipse-shaped trajectories, in reality their orbital motion x A (t) is highly complicated due to the mutual interactions among these bodies. The worldlines of all massive bodies can be concretized by Solar system ephemeris [63] at any stage of the calculations. One might prefer to series expand these worldlines as follows: where are the position, velocity and acceleration of body A at some time-moment t A . However, such an approach is problematic mainly for two reasons: (i) all terms of the infinite series expansion (16) contribute on 1PN or 1.5PN level, because the expansion in (16) is not performed with respect to the inverse powers of the speed of light.
(ii) the time-moment t A remains an open parameter as long as no additional arguments are introduced, which would uniquely allow to identify that parameter with the time of closest approach or with the retarded time.
These both aspects (1) and (2) enforce to determine the light trajectory x (t) of a light-signal from the celestial light source toward the observer as function of intrinsic multipoles M A L and S A L as well as function of the arbitrary worldlines x A (t) of these massive Solar system bodies.
In a previous investigation [26] a solution for the light trajectory in 1PN approximation in the gravitational field of N massive bodies in arbitrary motion and expressed in terms of their intrinsic multipoles has been obtained: However, as outlined in more detail in [26], such 1PN solution is not sufficient for astrometry on sub-µas-level of accuracy. For instance, the rotational motion of the massive bodies cannot be taken into account in 1PN approximation. However, the impact of the spin-dipole structure of the massive bodies on light deflection amounts to be about 0.7 µas , 0.2 µas , and 0.04 µas for a grazing lightray at Sun, Jupiter, and Saturn, respectively [28,36]. Moreover, also higher spin-multipoles have a significant impact on sub-µas-level [26,64]. Furthermore, in 1PN approximation there are no terms proportional to v A c M ab where M ab is the mass-quadrupole term. Already a straightforward estimate reveals that such terms become relevant on sub-µas-level of accuracy [26], see also Table III. In order to scrutinize the impact of such terms one is necessarily enforced to determine the 1.5PN solution for the light trajectory.
In view of these facts, the primary goal of this investigation is to determine a solution for the light trajectory in 1.5PN approximation, which includes all terms up to the order O c −4 , where both of the important aspects (1) and (2) addressed above are fully taken into account: Especially, the massive bodies of the Solar system are allowed to move along arbitrary worldlines x A (t) and they are having arbitrary shape and inner structure and rotational motion, given in terms of time-dependent intrinsic mass-multipoles M A L (t) and spin-multipoles S A L (t), in accordance with the IAU recommendations [19,20] and the theory of relativistic reference systems [18,[55][56][57][58][59][60][61][62]. The given solution for the light trajectory is considered as a further step towards a consistent model of generalrelativistic theory of light propagation in the gravitational field of the Solar system, which finally aims at accuracies on sub-µas-level and even on nas-level.
The article is organized as follows: In section II the geodesic equation in 1.5PN approximation is considered. A compendium of the DSX framework is presented in section III. The transformation of geodesic equation in terms of new variables, which are more efficient than the standard parametrization, is given in section IV. The first and second integration of geodesic equation is determined in section V and section VI, respectively. The important case of light-propagation in the gravitational field of moving spin-dipoles is investigated in section VII. Finally, the expressions for the observables of time delay and light deflection are obtained in section VIII and IX. Especially, numerical values for the impact of the leading mass-multipoles and spin-multipoles on time delay and light deflection are given in Table II and Table III, respectively. A summary and outlook can be found in section X. The used notations and conventions and further details and several checks of the calculations are shifted into appendix.

II. GEODESIC EQUATION IN 1.5PN APPROXIMATION
The Solar system is composed of N arbitrarily shaped, rotating and deformable massive bodies which move under the influence of their mutual gravitational interaction among their common barycenter. It is clear, that the metric of such a highly complicated N -body system is not known in its exact form and can only be determined within an approximative scheme. In view of the weak gravitational fields and slow motions of the bodies, the metric tensor of the Solar system in the BCRS coordinate system x µ = (ct, x) can be expanded in terms of inverse powers in the light-velocity, called post-Newtonian expansion [17]: where η αβ is the metric tensor of flat Minkowski spacetime and the metric perturbations are of the order h   (19) into (6a) yields the geodesic equation in 1.5PN approximation, which in terms of global coordinate time reads where a dot means total time-derivative. Note that the constraint in (6b) results inẋ (t) ·ẋ (t) c 2 = 1 + O c −2 , hence will not change the form of geodesic equation in 1.5PN approximation in (20).
In (20) we have taken into account that in general h In order to find an unique solution of the geodesic equation in (20), so-called mixed initial-boundary conditions can be imposed, which have extensively been used in the literature, e.g. [18,24,27,31,32,39,45]: The first condition (22) defines the spatial coordinates of the photon at the moment t 0 of emission of light. The second condition (23) defines the unit-direction (σ · σ = 1) of the lightray at past null infinity, that means the unittangent vector along the light path in the infinite past hence at infinite spatial distance from the origin of the global coordinate system. In the flat space-time there is no gravitational field, h (n) αβ = 0, hence the geodesic equation (20) simplifies to the formẍ (t) = 0, having the solution which is nothing else than just the unperturbed light trajectory in Eq. (1). The exact light trajectory x (t) in (9) deviates from the Newtonian approximation in (24) by terms of the order O c −2 , that means x (t) = x N (t) + O c −2 . Accordingly, in (20) we may replaceẋ (t) by its Newtonian approximation,ẋ N = c σ, and (20) simplifies as follows: ij,k σ j σ k which agrees with Eq. (3) in [32]; recall h 0i,0 = O c −4 . Furthermore, in 1.5PN approximation the metric perturbations in (25) can be taken at the spatial coordinates of the unperturbed lightray. That means, in (25) one has first to perform the differentiations with respect to BCRS coordinates x µ = (ct, x) and afterwards to insert the unperturbed lightray: In this investigation we will determine the solution of the geodesic equation (25) in 1.5PN approximation, which can formally be written as follows (cf. Eq. (18)): The 1PN corrections ∆x 1PN (t, t 0 ) in Eq. (27) are terms of the order O c −2 and have already been determined in our recent analysis [26]. Here, the primary goal is the determination of the 1.5PN corrections ∆x 1.5PN (t, t 0 ) in Eq. (27) which are terms of the order O c −3 .

III. COMPENDIUM OF DSX FRAMEWORK
The DSX framework represents a well-established formalism in the general-relativistic celestial mechanics of a N -body system of arbitrarily shaped, rotating and deformable bodies, and has been introduced and thoroughly formulated in [55][56][57][58]. The original intension of DSX was the description of the dynamics of N massive bodies, that is the equations of motion in celestial mechanics for N extended bodies under the influence of their mutual gravitational interaction.
The basic assumption is to introduce N + 1 reference systems: one global chart (BCRS) with coordinates x µ = (ct, x) having its origin of the spatial axes at the barycenter of the Solar system, and N local charts with coordinates X α A = (cT A , X A ), one for each individual body A = 1, ..., N and having their origins at the barycenter of these massive bodies and co-moving with them. The local coordinate systems are tied to each individual massive body and are defined very similar to the Geocentric Celestial Reference System (GCRS) which is in use for the Earth and, therefore, they are called GCRSlike reference systems. A central result of the DSX approach is the form of the global metric g µν of BCRS and the form of the local metric G A αβ for each GCRS-like system, and the coordinate transformation among all these reference systems (ct, x) ↔ (cT A , X A ). Another central achievement in the DSX formalism is the decomposition of the global metric in terms of intrinsic mass-multipoles M A L and intrinsic spin-multipoles S A L . In this section we will present a compendium of the DSX theory, which has become a basic part of IAU resolution B1. 3 (2000) [19] and which are of upmost relevance for our own considerations aiming at applications of the DSX approach in the astrometrical science.

A. BCRS
The harmonic BCRS coordinates are denoted by x µ = ct, x i , where t = TCB is the BCRS coordinate time; about a practical synchronization of a set of clocks distributed somewhere in the Solar system we refer to [65]. The origin of the spatial axes of BCRS is located at the barycenter of the Solar system and cover the entire threedimensional space and can therefore be used to model light trajectories from distant celestial objects to the observer. The IAU Resolution B2 (2006) [20] recommends the spatial axes of BCRS to be oriented according to the spatial axes of the International Celestial Reference System (ICRS) [66]. Furthermore, according to IAU resolution B1. 3 (2000) [19] the Solar system is assumed to be isolated and the space-time is asymptotically flat, that means the BCRS metric g µν (t, x) at infinity reads: where r = |x|. The BCRS is completely characterized by the form of its metric tensor which, however, is not known in its exact form. According to the geodesic equation in 1.5PN approximation (25), for our intentions the metric is required to be known only up to terms of the order O c −4 , which are given by [19]: The gravitational potentials in (29) -(31) are given by the integrals where the integrals in (32) and (33) run over the entire Solar system, and t µν is the energy-momentum tensor of the Solar system in global BCRS coordinates; recall the components of energy-momentum tensor scale as follows: The global gravitational potentials in (32) -(33) admit an expansion in terms of global Blanchet-Damour (BD) mass-multipoles and spin-multipoles, m L and s L , [32,[51][52][53], which characterize the multipole structure of the Solar system as a whole, The global mass-multipoles and global spin-multipoles in (34) - (35) are Cartesian symmetric and trace-free (STF) tensors, and up to order O c −2 given by (cf. Eqs. (2.34a) and (2.34b) in [52]): where m 0 = const. is the mass of the entire Solar system, the mass-dipole m i = 0 because the origin of BCRS is located at the barycenter of the Solar system, and the spin-dipole s i = const describes the spin of the entire Solar system which safely can be assumed to be timeindependent. The spatial derivative operator in (34) - (35) is defined by and a dot means derivative with respect to global coordinate-time. The expansion in (34) -(35) has two specific features, which do not allow a straightforward application in our investigations: (1) The expansion in (34) -(35) is valid outside a sphere which encloses the N -body system [51][52][53][54]67]. It is quite obvious that for modeling of light propagation through the Solar system we need to have a metric which is valid in spatial domains between these N massive bodies.
(2) It has already been underlined in the introductory section that according to the theory of reference systems [18,19,[55][56][57][58][59][60][61][62] physically meaningful multipole moments of some massive body A have to be defined in the local reference system (cT A , X A ) tied to that body and comoving with it.
For these reasons, the global gravitational potentials in (34) -(35) must have to be expressed by intrinsic massmultipoles M A L and intrinsic spin-multipoles S A L , which are defined in the local reference system (cT A , X A ) of each individual massive body A. The prototype of all these GCRS-like coordinate systems is the GCRS especially designed for the Earth and which will be considered now.

B. GCRS
For the description of physical problems nearby the Earth the GCRS is the appropriate reference system. The harmonic GCRS coordinates are denoted by X α = cT, X i , where T = TCG is the GCRS coordinate time. According to IAU resolution B1. 3 (2000) [19], the origin of the spatial axes of GCRS is located at the center-ofmass of the Earth and co-moving with it. The spatial axes of GCRS are kinematically non-rotating with respect to the BCRS, that means the GCRS is a spacefixed reference system and is not a local inertial system. The GCRS is completely characterized by the form of its metric tensor, up to order O c −4 given by [19,55,56], The gravitational potentials in (39) -(41) can uniquely be separated into two components: a local component, W loc , W i loc which originates from the body A itself and an external component, W ext , W i ext , which is associated with inertial effects and tidal forces [19,55,56]: Explicit expressions for the external potentials are given in [55,56], while the local potentials are defined by the following integrals: where the integrations run over the entire volume of the Earth, and where T µν are the components of the energy-momentum tensor in GCRS coordinates; recall the components of energy-momentum tensor scale as follows: The local potentials (44) -(45) generated by the Earth can be expanded into a series of STF multipole moments, M L and S L . In the harmonic skeletonized gauge they are given by [19,[51][52][53][54][55]: where R = |X| is the spatial distance from the origin of GCRS to some field point outside the Earth, and and a dot in (47) denotes a derivative with respect to GCRS coordinate time T . The intrinsic STF multipoles in (46) and (47) can be approximated by their Newtonian expressions, that means up to terms of the order O c −2 they are given by: where the integration runs over the volume of the Earth, and T αβ is the energy-momentum tensor in the local system of the Earth. The intrinsic mass-monopole term M = const in (49) is the Newtonian mass of the Earth. Actually, the mass-dipole vanishes, M i = 0, because the origin of the GCRS is assumed to be located at the barycenter of the Earth, but in real measurements of celestial mechanics the center-of-mass of massive Solar system bodies can not be determined exactly, so it is meaningful to keep this term and to assume M i = const. in general. The spin-dipole S i (T ) of the Earth is not constant but timedependent due to inner forces of the Earth and due to the gravitational interaction of the Earth with other massive bodies.

C. Metric of Solar system in terms of intrinsic multipoles
In order to describe the light trajectory through the Solar system, one needs to introduce one global chart (BCRS) x µ = (ct, x) but expressed in terms of intrinsic multipoles, M A L and S A L , of each massive body A = 1, ..., N . For being able to define the multipole structure of each individual body in a physically meaningful manner, the DSX formalism [55,56] introduces N local GCRS-like reference systems X α A = (cT A , X A ), each one very similar to the GCRS in Eqs. (39) - (41). These N +1 coordinate systems are linked with each other via coordinate-transformations. The DSX theory [55][56][57][58] provides the theoretical framework for such an approach, and has originally been established for celestial mechanics and for deriving the equations of motion of a system of N massive bodies with full multipole structure. Consequently, one central result of DSX theory are the coordinate transformations among these reference systems, which are given by where x µ A (T A ) is the worldline of body A in BCRS coordinates. The inverse coordinate transformations could be found in [19], but is not of relevance here for our purposes. The tetrads e µ a along the worldline of this body are explicitly given by (cf. Eqs. (2.16) in [55]): where in (52) a dot means derivative with respect to the local coordinate time of body A. That means,ẋ A (T A ) is the three-velocity of body A in the global system and given in terms of the body's local coordinate time T A , which could easily be transformed into terms of global BCRS coordinate-time.
The contravariant components of the BCRS metric tensor g µν (t, x) in Eqs. (29) -(31) and the contravariant components of the metric tensor G αβ (T A , X A ) in Eqs. (39) - (41) in the local GCRS-like coordinate system of body A are related via the following transformation: Using (51) in virtue of (54), the global potentials w, w i in (32) - (33) can be expressed in terms of intrinsic STF multipoles M A L and S A L as follows [19,55,56]: In (55) and (57) the sum runs over all bodies of the Nbody system, R A = |X A | is the spatial distance from the origin of local coordinate system to some field point located outside the massive body, and The local mass-multipoles and spin-multipoles of some massive body A in Newtonian approximation, i.e. up to terms of the order O c −2 , are given by (cf. Eqs. (49) - (50): where the integration runs over the spatial volume of massive body A, and T αβ A is the energy-momentum tensor of body A in the local coordinate system of that body.
Finally, using coordinate transformations (51), the spatial derivatives in (56) and (58) must be transformed into the BCRS, where they read as follows: where the second term in (62) as well as in (63) generate terms which are beyond the order of O c −4 in the global metric, that means these terms will actually not contribute in the final results for the light trajectory. From (51) follows the relation [19,21,55,56]: where we recall that some massive body A moves along the arbitrary worldline x A (t), which can later be concretized by Solar system ephemeris [63] at any stage of the calculations. Because of the fact that the BCRS coordinate-time and the coordinate time T A of local system of body A are related as follows [19,21,55,56], we obtain for the time-dependence of the intrinsic multipoles the following relation: that means the neglected terms in (66) and (67) are beyond 1.5PN approximation for lightrays. Summarizing the conclusions in Eqs. (51) - (67), the metric perturbation in the near-zone of the Solar system and expressed in terms of local multipoles is given by: where the summmation in (68) and (71) runs over all massive bodies of the Solar system, while the metric perturbations caused by one individual body are given by (69) and (72). The dot in the first term of expression (72) means here derivative with respect to global BCRS coordinate time, and the spatial derivatives in (68) - (72) are derivatives in the global system and given by and is the distance between some field-point with spatial coordinate x and the spatial position x A (t) of massive body A in the global reference system at BCRS time t. The metric perturbations in (68) - (72) have to be implemented into the geodesic equation in (20) and, therefore, the field-point x in (74) will be identified with the photons position x (t) according to Eq. (21). In view of this fact we will use the same notation for the distance in (74) and for the absolute value of (B10). Before going further, we underline the absence of terms (72). Such terms are of the order O c −4 in the metric, that means they are beyond the 1.5PN approximation for lightrays.

IV. TRANSFORMATION OF GEODESIC EQUATION
As it has been discussed above, instead of (20) we actually may consider the simpler form of geodesic equation in (25), which is integrated along the unperturbed light trajectory (1). That means, according to Eq. (26), the field-point x in Eq. (74) can be approximated by the unperturbed photon-trajectory x N (t) in (1), so that we get the following expression for the vector pointing from the center of massive body A toward the spatial position of the photon at time t: where the unperturbed lightray is given by Eq. (1) or Eq. (24). It especially means, that all derivatives in geodesic equation (25) and in the metric perturbations in (68) -(72) act on the unperturbed lightray. In view of this important fact, it is highly effective to embark on a strategy, where all expressions in the geodesic equation (25) are expressed in terms of new parameters which fully characterize the unperturbed light trajectory from the very beginning of the integration procedure. This strategy especially implies, that we will transform the spatial derivatives in (73), the derivatives in the geodesic equation (25), the distance in (76) and the time-argument of the multipoles in terms of these new parameters.  78) and points from the origin of global system to the point of closest approach of the unperturbed lightray to that origin, and is time-independent. The impact vector dA (τ + t * ) is defined by Eq. (85) and points from the origin of local system of body A toward the point of closest approach of unperturbed lightray to that origin, and is time-dependent due to the motion of the body. Furthermore, x (τ + t * ) is the global spatial coordinate of the photon of the light trajectory, while xN (τ + t * ) is the unperturbed lightray. The worldline of massive body A in the global system is given by xA (τ + t * ), and rA (τ + t * ) points from the origin of local system toward the exact photon's position, while r N A (τ + t * ) points from the origin of local system toward the unperturbed lightray.
The problem and the need for introducing new variables is namely the following. The variables t and x are field variables of the gravitational field and, therefore, they are of course independent of each other. But since the integration of geodesic equation proceeds along the lightray (cf. Eq. (21)) these field variables have to be replaced by the photon trajectory, x (t), and then these variables become dependent on each other. A drastical simplification is achieved in view of Eq. (26) which states that the geodesic equation in (25) can be integrated along the unperturbed lightray. Therefore, we are looking for new time-variable and spatial-variable, which fully parametrize the unperturbed lightray and which are independent on each other. In this way the integration of geodesic equation becomes feasible. Just for that reason, the following independent variables τ and ξ have been introduced in [32,39,40]: where P i j is the operator of projection onto the plane perpendicular to vector σ, where the covariant and contravariant positions of spatial indices is insignificant: P i j = P ij = P ij . According to (78), the three-vector ξ is the impact vector of the unperturbed lightray, see also Eq. (B2). Especially, ξ is time-independent and directed from the origin of global coordinate system toward the point of closest approach of the unperturbed light trajectory and the absolute value is denoted by d = |ξ|. For a graphical elucidation see Fig. 1. Another important parameter is the time of closest approach of unperturbed lightray to the origin of the global coordinate system, defined by which differs from (B14) which is the time of closest approach of the lightray to the origin of the local coordinate system of some massive body A. Notice that dt = dτ for the total differentials, because t * is a constant for each particular lightray, and τ = t − t * and τ 0 = t 0 − t * .
With the aid of these new variables ξ and τ , the mixed initial-boundary conditions (22) and (23) take the form where a dot means derivative with respect to variable τ . In the new variables the interpretation of these initialboundary conditions remains the same: the first condition (81) defines the spatial coordinates of the photon at the moment of emission of light, while the second condition (82) defines the unit-direction (σ · σ = 1) at infinite past and infinite distance from the origin of global coordinate system, that means at the so-called past null infinity. The unperturbed lightray in (24) transforms as follows [26,32,33,39,40]: while its derivative with respect to variable τ readṡ x N (τ + t * ) = c σ. The vector pointing from the spatial position of the arbitrarily moving body toward the unperturbed lightray in these new variables transforms as follows: with the absolute value r N A (τ + t * ) = r N A (τ + t * ) , and the impact parameter in (B13) for arbitrarily moving bodies in these new variables reads: In virtue of Eqs. (77) and (78) two new variables, τ and ξ, have been introduced and in addition the auxiliary variable t * by Eq. (80). As next, the partial derivatives with respect to space and time in the geodesic equation (25) have to be expressed in terms of these new variables. In the pioneering investigations in [32,39,40] it has been shown by chain rule that these partial derivatives transform in the following way: Two remarks are in order to interpret these relations correctly.
First, we notice that the explicit time-dependence of the metric tensor, h Second, it should be realized, that in the left-hand side in (86) and (87) one has first to perform the differentiations and afterwards the field-point x has to be substituted by the unperturbed lightray x N (t) = x 0 + c σ (t − t 0 ). Opposite, in the right-hand side in (86) and (87) one has first to substitute t * + τ and x N (τ + t * ) = ξ + c τ σ and afterwards to perform the differentiations. By means of these relations (86) and (87), the geodesic equation in 1.5PN approximation in (25) transforms as follows: which agrees with Eq. (36) in [32] and Eq. (19) in [33]; note that P ab σ b = 0. The double-dot on the left-hand side in (88) means twice of the total differential with respect to the new variable τ . Subject to relation (70), the geodesic equation in (88) simplifies further: Let us note that the first two terms are of order O c −2 and agree with Eq. (95) in [26], while the last two terms are of order O c −3 . This fact implies that if one integrates the geodesic equation (89) then the first two terms in (89) give rise to terms of the order O c −2 as well as to terms of order O c −3 , while the last two terms generate only terms of the order O c −3 . The mathematical structure of (89) is considerably simpler than the original form in (25), but of more decisive importance in the integration procedure is the fact that the time-variable τ and the space-variable ξ are independent of each other.
As final step in the transformation, the metric perturbations in (68) - (72) have to be transformed in terms of these new variables ξ and τ . One obtains and where the sum in (90) and (92) runs over all massive bodies of the Solar system. The expressions in Eqs. (91) and (93) contain the STF spatial derivative operation ∂ L , which also has to be expressed in terms of these new variables. That issue is considered in detail in appendix C and yields the following expression for the STF partial derivative operation in Eqs. (91) and (93): (94) These expressions in (90) -(94) have to be inserted into the geodesic equation (89), which finally yields the geodesic equation for lightrays which propagate in the gravitational field of one arbitrarily moving body A in terms of these new variables τ and ξ: where the indices M and S stand for mass-multipole and spin-multipole component, respectively. That means, the linearity of geodesic equation in 1.5PN approximation allows simply to sum over all N arbitrarily moving bodies just straight away. The contributions due to the massmultipole structure of one body A is given bÿ and the contribution due to the spin-multipole structure of one body A reads where the derivative operator is given by (94). By Eqs. (95) -(97) the transformation of geodesic equation in 1.5PN approximation in terms of these new variables τ and ξ has been accomplished, which describes the propagation of a light-signal through the field of N massive bodies in arbitrary motion and having arbitrary shape and inner structure and which can also rotate arbitrarily. Before we proceed further, three comments should be in order: (i) First, let us note that the spatial derivative operator in (94) depends on time-variables τ and t * , but in such a way that it does not act on time-dependent multipoles or the velocity of the body, that means: because the construction of the derivative operator in (94) is such that the derivatives with respect to variable τ cancel exactly the derivatives with respect to t * in all those functions which depend on the combination τ + t * . But of course ∂ L r N A (τ + t * ) = 0. (ii) Second, let us also remark that in (96) the STF notation for the derivative operator has been kept. But we recall the following relation, which is a specific example of the more general relation Eq. (A1) in [68]: The relation in (101) has allowed to replace the STF derivative operator ∂ L by ∂ L in Eqs. (100) -(102) in [26]. Here, in view of relation (101) and (102) we may also replace the STF derivative operator ∂ L by ∂ L in all terms in (96), and correspondingly in the first integral in (110) and (111), as well as in the second in (118) and (119). On the other side, such replacement is not possible for the spin-multipole terms in (97), because of (iii) Third, it should also be mentioned that in the limit of one massive body at rest with the origin of the coordinate-system located at the center-of-mass and with time-independent multipoles then the geodesic equation (95) -(97) agrees with the geodesic equation given in [31], notice the comment in [69].

V. FIRST INTEGRATION OF GEODESIC EQUATION
The coordinate velocity of the photon is determined by the first integral of geodesic equation (89). In terms of the new variables we may separate the first integral of geodesic equation (89) into 1PN and 1.5PN terms as follows:ẋ That means, according to (95) we may consider the lightpropagation in the field of one arbitrarily moving body A and finally we have to build the sum over all massive bodies A = 1, ..., N in order to obtain the light trajectory in the entire Solar system. Furthermore, according to Eq. (95) we split these expressions into mass-multipole terms and spin-multipole contributions as follows: where we have taken into account that in (105) there are no spin-multipoles because they are of the order O c −3 , hence they do appear only in (106). We shall consider mass-multipole and spin-multipoles in the next both subsections separately.
A. First integration for mass-multipoles The first integral of geodesic equation (95) for the mass-multipole component of one massive body A reads: where the integrand up to the required order is given by Eq. (96). Let us underline that the integration of the first expression on the r.h.s. in (96) yields terms of the order O c −2 as well as terms of the order O c −3 . For that reason, the integral in (107) is written as sum of 1PN and 1.5PN terms. In particular, for the integration of geodesic equation the following rules are important (cf. Eqs. (4.9) and (4.10) in [39] or Eqs. (4.38) and (4.39) in [40]): where the function C (ξ) in (108) depends only on variable ξ, thence disappears in case of definite integrals. The rule in (108) and (109) [40].
The integration of the first expression on the r.h.s. in (96) is shown in more detail in appendix E, while in view of relation (108) the integrals of the second, third, and fourth expression in (96) are straightforward. The fifth term in (96) can be integrated by parts using relation (D9) and is shown in more detail in appendix F, while the integration of the sixth term goes very similar. Altogether, for the 1PN terms one obtains: In the first term on the r.h.s. in (110) we have used relation (D5). For the 1.5PN terms one obtains: where we recallṀ A = 0. In the second and sixth term on the r.h.s. in (111) [26]. In [26] it has been demonstrated that in case of bodies at rest (v A = 0) having timeindependent mass-multipoles and located at the origin of coordinate system (x A = 0) our result in (110) agrees with the time-derivative of Eqs. (33) and (36) in [31]. It should also be noticed that the derivative of (111) with respect to variable cτ yields the expression in (96).

B. First integration for spin-multipoles
The first integral of geodesic equation (95) for the spinmultipole component of one massive body A reads: where the integrand up to the required order is given by Eq. (97). The integration in (112) can be performed straightforward and one obtains: Let us remark that the second term in (113) is obtained by integration by parts, using (D9) and afterwards mak-ing use of relation (D5). Note that the derivative operator ∂ L in terms of the new variables ξ, τ , t * is given by (94). In the appendix H it is shown that in the limit of bodies at rest and stationary spin-multipoles our result in (113) agrees with Eqs. (32) and (37) in [31], up to an overall sign which has been clarified [69]. We also note that the derivative of (113) with respect to variable cτ yields the expression in (97). Let us remark that neglecting terms of the order in Eqs. (110) and (111)

VI. SECOND INTEGRATION OF GEODESIC EQUATION
The light trajectory of the photon is determined by the second integration of geodesic equation (89), and can be written as follows: where the sum runs over all massive bodies A = 1, ..., N of the Solar system. Like in the case of first integration, we split these expressions into mass-multipole terms and spin-multipole contributions as follows: where in (115) there are no spin-multipoles because they are terms of the order O c −3 and consequently they do appear only in (116). We will consider the massmultipole and the spin-multipole components separately.

A. Second integration for mass-multipoles
The mass-multipole terms in (115) and (116) read For the mass-multipole terms to order O c −3 one ob-tains: Notice, that the derivative operator ∂ L in (118) and (119) in terms of the new variables ξ, τ , t * is given by (94). One may demonstrate, that (118) and (119)  The 1PN solution in Eq. (118) coincides with Eq. (137) in [26]. Recall that in 1PN approximation the derivative operator in Eq. (94) can be replaced by the simplified expression in Eq. (C6) (cf. Eq. (101) in [26]), because derivatives with respect to variable t * generate terms of the order O c −3 ; see also comments below Eq. (C5). Furthermore, in [26] it has already been shown that in case of bodies at rest and located at the origin of coordinate system our result in (118) agrees with Eqs. (33) and (36) in [31].

B. Second integration for spin-multipoles
The spin-multipole terms in (116) read where the integrand in (120) is given by the expressions in Eq. (113). The second expression on the r.h.s. in Eq. (113) is rewritten by means of relation (D5) and then, by means of relations (D10) and (D11), we may integrate by parts. We obtain the following solution: where for the second expression we also have used relation (D6). The derivative operator ∂ L in terms of the new variables ξ, τ , t * is given by (94). One may easily check, that (121) is consistent with (113), in the sense that the derivative of (121) with respect to variable τ just yields the expression in (113) up to terms of the order O c −4 . For that proof simply apply the relations (D9) and (D12) and take into account (G3). Furthermore, in appendix H it is shown that in the limit of bodies at rest and timeindependent spin-multipoles our result in (121) agrees with Eqs. (33) and (38) in [31], up to an overall sign which has been clarified [69]. We underline again that neglecting terms of the order in Eqs. (118) and (119)

VII. LIGHT TRAJECTORY IN THE FIELD OF SPIN-DIPOLES
In our previous investigation [26] the light trajectory in the field of N arbitrarily moving mass-monopoles, massdipoles, and mass-quadrupoles has been considered as specific examples of the general solution, see Eqs. (139), (140), and (143) -(148) in [26], respectively. Here we will consider the light trajectory in the field of N arbitrarily moving spin-dipoles as specific example of the general solution. It may also serve as a further instructive example about how the presented approach runs.
A. Light trajectory in the field of N arbitrarily moving spin-dipoles The rotational motion of a real body like the Sun, Earth, or Jupiter, is a highly complicated physical subject, because these bodies are not rigid monopoles and the rotational motion can therefore not be described by a simple spin-dipole, but must be expressed by the full set of time-dependent spin-multipoles S A L (t) with l = 1, 2, 3, .... On the other side, the main impact among all spin-multipoles on light deflection is of course given by the first summand in (121) which is proportional to the intrinsic spin vector S A (t) of body A and which is called spin-dipole. It is also well-known that for submicro-arcsecond astrometry the light trajectory in the field of a spin-dipole is of specific importance, because the light deflection caused by the spin-dipole of a body at rest amounts to be 0.7 µas for grazing rays at the Sun, 0.2 µas for grazing rays at Jupiter, and 0.04 µas for grazing rays at Saturn [28,36]. Therefore, we will consider the light trajectory in the field of one arbitrarily moving body with time-dependent spin-dipole in more detail in this section.
According to Eq. (114) with (120) and (121), the light trajectory in the field of N arbitrarily moving spindipoles reads: where .
The derivative operator in terms of the variables ξ, τ, t * is given by (94), which for one index reads: By inserting (124) into (123), we encounter the following individual terms: and and and we recall jab σ a σ j = 0. Inserting (125) -(127) into (123) yields where the notation ijk a j b k = (a × b) i has been used. The complete expression for the light trajectory in 1.5PN approximation in the field of N arbitrarily moving and time-dependent intrinsic spin-dipoles is finally obtained by inserting (128) into (122).
As mentioned in the introductory section, in [38] the light trajectory in post-Minkowskian approximation in the field of N arbitrarily moving pointlike spin-dipoles has been determined. That means, the pointlike objects in [38] may even be in ultra-relativistic motion, while our 1.5PN solution in (128) is valid for extended bodies with spin-dipole but in slow-motion along arbitrary worldlines. In appendix I it is shown that our result in (128) agrees with the results in [38] for the light trajectory up to terms of the order O c −4 . One may also verify that in the limit of time-independent spin-dipoles, S A = const, and in the limit of uniform motion, v A = const, our result in (128) agrees with Eq. (26) in [37] in GR, noticing that constant terms cancel each other according to Eq. (122).

B. Light trajectory in the field of N bodies at rest with spin-dipole
In this section we will consider the case of light propagation in the field of N spin-dipoles at rest and compare with results in the literature. For time-independent spin-dipole S A = const. and for one body at rest located at x A = const. in the global reference system we have where d A is the time-independent impact vector defined by Eq. (B2). From (122) we obtain the light trajectory in the field of N bodies at rest with time-independent spin-dipoles: where from (128) we obtain the following expression for the correction-term: where a time-independent term 2 G c 3 has been omitted because this term will be cancelled in view of (129). The time-dependence of (129) and (130) is solely caused by the time-dependence of the unperturbed lightray in (83). In order to obtain the form of the expression in (130) we have also used d 2 The expression in (129) -(130) agrees with the solution in Eq. (56) in [36], where the trajectory of a photon as function of time has been determined in the field of N bodies at rest in post-Newtonian approximation for the lightrays. It is straightforward to show that the time-derivative ∂ c τ ∆x S A (τ ) coincides with Eq. (59) in [36].

VIII. TIME-DELAY
In the previous sections we have determined the light trajectory of a light-signal which propagates through the metric field of the Solar system, that means through the gravitational field of N arbitrarily moving massive bodies. However, the light trajectory is not an observable at all. In real astrometric measurements one of the most important observable quantity concerns the time delay of some light-signal propagating in the Solar system. The considerations here are similar to what has been discussed in [26] about observable effects, but with the extension to 1.5PN approximation. Especially, we will assume that the light source is located at x 0 = x (t 0 ) where t 0 is the moment of emission of the light-signal, and the observer is located at x 1 = x (t 1 ) where t 1 is the moment of reception of the light-signal by the observer. Furthermore, both the light source and the observer are assumed to be at rest with respect to the global reference system.
In the pioneering work [70], Shapiro has considered the general-relativistic effect of time delay of a light-signal which propagates through the gravitational field of a static and spherically symmetric massive body. Especially, Shapiro has drawn the attention to the fact about the measurability of that additional test of relativity by radar technology. In fact, the Shapiro time delay was discovered soon afterwards [71]. It might be useful to realize that the reason for the time delay is not only laying upon the fact that the light-trajectory is curved but also because the speed of a photon is decelerated in the gravitational field of a monopole at rest. While the classical Shapiro effect is originally related to a time delay of a light-signal in the monopole-field, it became a matter of common knowledge to call the time delay of a light-signal in any gravitational field just Shapiro effect.
For describing the Shapiro-effect, we introduce a vector pointing from the light source at the moment of emission toward the observer at the moment of reception, which in terms of the new variables reads where k is just the corresponding unit vector with R = |R| being the absolute value of R. Using very similar steps as in [32], we obtain from Eq. (114) the following expression for the time delay in the gravitational field of N arbitrarily moving massive bodies in 1.5PN approximation that means up to terms of the order O c −4 : where the sum runs over all massive bodies and the expressions for ∆x A 1PN and ∆x A 1.5PN are given by Eqs. (115) and (116) with (118), (119) and (121), respectively. The 1.5PN relation (133) generalizes the 1PN relation (154) in [26].
We will consider the time delay in (133) of a light-signal caused by N arbitrarily moving bodies in some more detail, but will restrict ourselves on the case of N moving bodies with monopole-structure (M) , quadrupolestructure (J 2 ), and spin-dipole-structure (S). Higher multipoles are so tiny that they are negligible in the time delay effect. These first terms in the general formula (133) Furthermore, as mentioned in the introductory section, there are several proposals to ESA for future spacebased missions, like ASTROD [8,9], LATOR [10,11], ODYSSEY [12], SAGAS [13], TIPO [14], which aim at time-transfer accuracies of two separated clocks within the Solar system of up to 10 ps. The question arises about the ability of such extremely-precise astrometry missions, especially designed for tests of relativity in the Solar system, to detect some 1.5PN terms in the Shapiro effect which will be discussed in this section.
In general, the light-signal will be assumed to be emitted at a space-time point with BCRS coordinates x 0 , t 0 and received by an observer at a space-time point with BCRS coordinates x 1 , t 1 . We also introduce the following notations: . Furthermore, we notice that σ = k + O c −2 according to Eq. (167) given below, that means we may replace the vector σ in favor of vector k whenever it is reasonable.

A. Moving mass-monopole
We will consider the time delay in (133) of a light-signal caused by an arbitrarily moving monopole.
In the limit of monopoles at rest only the term in Eq. (138) remains which then represents the well-known classical Shapiro effect [17,18,21,22] which is growing logarithmically with R, while in our result (139) the argument of the logarithm depends on the worldline of the arbitrary-moving body x A (t). One may verify that our result for the Shapiro delay for arbitrarily moving monopoles in Eq. (138) -(139), agrees in the limit of uniform motion with Eq. (20) in [72], with Eq. (45) in [42], and with Eq. (33) in [43] up to terms of the order O c −4 . In this respect we recall that the term in the last line in (139) can be written as follows: where the term in the last line is proportional to the acceleration of the massive body A and vanishes in case of uniform motion. The neglect of this term, as suggested in [33], is well-justified because a simple estimate reveals that such terms are extremely small and far out of detectability even for future astrometry missions. An estimate of the absolute value of the 1PN time delay formula in Eq. (138) for one body A and assuming an astrometric configuration with σ ·r 0 A −r 0 A and σ ·r 1 A r 1 A , is given by [22]: A very similar estimate of the absolute value of the 1.5PN correction in Eq. (139) for one body A and same configuration yields The second term in (142) is proportional to ∼ r 1 A /d 1 A , which for grazing rays becomes a large quantity. For instance, for Jupiter we would get r 1 A /d 1 A ∼ 10 4 which spoils the effect of the tiny factor v A /c ∼ 10 −5 which is typical for 1.5PN corrections. This large term is solely caused by the term in the last line in Eq. (139). Below, we will consider the expressions for light deflection where we will encounter this large term again, cf. text below Eq. (34) in [28]. As we will show in the next subsection, this large factor r 1 A /d 1 A is related to the retardation of gravitational action.

In terms of retarded time
Gravitational action travels with the finite speed of light and this effect cannot be ignored in high-precision astrometry, as it has been outlined long time ago [28,30,33,35,38]. In order to take account for that effect we follow the arguments of the investigations in [28,33,35,38,73], which have shown that the position of the massive body must not be taken at the time of observation, x A (t 1 ), but at the retarded time-moment, x A (t ret 1 ). In general, the retarded time is defined by an implicit relation, where t is the coordinate time. For the special case where t is the time of emission t 0 or the time of reception t 1 see Eq. (149). Actually, the retarded time is a function of the position of body under consideration and, therefore, an index A should also be attached at t ret but for simpler notation such label is omitted. According to Eqs. (47) - (48) in [74], the retarded position can be series-expanded and leads to the following relations for any instant of time: These relations allow one to rewrite identically the expressions in (137) -(139) into the following form up to terms of the order O c −4 : where (cf. Eq. (143)): The solution for the time delay in (146) -(148) agrees with Eq. (51) in [33]. Especially, we notice that the term in the last line of Eq. (139) has been absorbed in (147). Consequently, if one uses the expression for the time delay in terms of retarded time, Eqs. (146) -(148), then one obtains the following correct estimate for the time delay in 1.5PN approximation: For numerical values of the upper bound in Eq. (150) and Eq. (151) see Table II.

B. Moving spin-dipole
Now let us consider the time delay in (133) of a lightsignal caused by N arbitrarily moving spin-dipoles.

In terms of coordinate time
From (121) we obtain in the field of arbitrary-moving spin-dipoles (l = 1) the expression for ∆x A S 1.5PN , as given by Eq. (128). According to Eq. (135) we obtain for the   [63]. The value for J A 2 for the Sun is taken from [75], while J A n with n = 2, 4, 6 for Jupiter and Saturn are taken from [76], while J A n with n = 8, 10 for Jupiter and Saturn are taken from [77] and [78], respectively. The spin angular momenta SA are determined from the moment of inertia IA with the ratio IA MAP 2 A = 0.059, 0.254, 0.210 for Sun, Jupiter, Saturn, respectively from NASA planetary fact sheets. For the distance between light-source and body we assume r 0 A = 10 13 m so that the light-source is within the near-zone of the Solar system, while r 1 A is computed under assumption that the observer (spacecraft) is located at Lagrange point L2, i.e. 1.5 × 10 9 m from the Earth's orbit.
Shapiro-delay the following expression up to terms of the order O c −4 : where S 1 A = S A (t 1 ) and S 0 A = S A (t 0 ) are the spindipoles of body A at time observation-time t 1 and at emission-time t 0 respectively. It can be checked that in the limit of bodies at rest our result in (152) agrees with Eq. (72) in [36]. Furthermore, by very similar steps as used in appendix I one may verify an agreement of our solution in Eq. (152) with Eqs. (48) - (50) in [38] in case of slow motion; note that the global spin-tensor in [38] has to be reexpressed in terms of intrinsic spin-dipole, for instance by means of the relations Eqs. (B.8) and (C.10) in [74] and the retarded time has to be series-expanded in terms of global coordinate-time. An estimate of the upper bound of Eq. (152) yields which agrees with the estimate in Eq. (75) in [36] for grazing rays and spin-dipoles at rest.

In terms of retarded time
In view of relations (144) -(145) and up to terms of the order O c −4 one may perform the following replacements in Eq. (152): The upper bound is then given by For numerical values of the upper bound in Eq. (157) see Table II.   [8,9], LA-TOR [10,11], ODYSSEY [12], SAGAS [13], or TIPO [14], which aim at an accuracy in the determination of time delay for a light-signal better than ∆t ∼ 0.1 ns. Accordingly, 1.5PN effects in time delay will surely not be detectable even within the very next generation of high-precision space-based astrometry missions.
Actually, the expression ∆x J2 1PN has already been presented in its explicit form by Eq. (144) in [26]. In view of their involved structure, ∆x J2 1PN as well as ∆x J2

1.5PN
will not be given here. The estimate of (159) and (160) proceeds very similar to what has been done in detail in [79]. For an axisymmetric body one obtains after some amount of algebra: where J 2 is the actual coefficient of second zonal harmonics. The estimate in (161) agrees with the estimate for quadrupoles at rest, cf. Eq. (26) in [79]. Like in (142), we encounter in (162) we encounter a large term which is proportional to ∼ r 1 A /d 1 A .

In terms of retarded time
With the aid of relations (144) -(145) one rewrites ∆x J2 1PN (t) and ∆x J2 1.5PN (t) in terms of retarded time. Formally, one may also replace M A ab (t n ) by M A ab (t ret n ) , n = 0, 1, but the impact of such replacement on time delay is negligible. Then, after considerable amount of algebra, one obtains the correct estimates in 1.5PN correction, which are given by: The numerical magnitude of the 1PN correction in (163) is given in Table II, while the 1.5PN correction in (164) is by far much below the detectability of future astrometry missions and will not be given in Table II.
In view of the tininess of ∆ t J2 1PN it becomes obvious that higher multipole terms are negligible in the time delay and, therefore, will not be considered here.

IX. LIGHT-DEFLECTION
The light deflection is of fundamental importance in astrometric measurements. Like in the previous section, we assume the light source to be located at x 0 = x (t 0 ) where t 0 is the moment of emission of the light-signal, and the observer is located at x 1 = x (t 1 ) where t 1 is the moment of reception of the light-signal by the observer. Both the light source and the observer are assumed to be at rest with respect to the global reference system.
The light deflection is defined by the angle ϕ between unit vector k and the unit tangent vector n of the lightray at the observers position: ϕ = arcsin |k × n| [28]. In 1.5PN approximation the unit tangent vector at the observer is given by By inserting Eq. (104) into (165), we obtain where ∆ẋ A 1PN and ∆ẋ A 1.5PN are given by (105) and (106), respectively, with the expressions in Eqs. (110), (111) and (113). The 1.5PN relation (166) generalizes the 1PN relation (156) in [26].
The expression in (166) for the unit tangent vector along the light trajectory at observers position is valid in case of stars, which means in case of light sources which are at far distances from the observer. For astrometry within the Solar system we need to obtain an expression which is valid for light sources at finite distances from the observer. In order to obtain such an expression we use the following relation among the vectors k and σ, where ∆x A 1PN and ∆x A 1.5PN are given by Eqs. The relation follows from the definitions (131) and (132) and with the aid of the expression for the light trajectory in (114) and for the Shapiro effect in (133). The 1.5PN expression in (167) generalizes the 1PN relation (157) in [26]. We also notice that the first line in (167) agrees with Eq. (66) in [28]. By inserting (167) into (166) we finally arrive at the following expression for the unit tangent vector at the observers position: The 1.5PN relation in (168) generalizes the 1PN relation (158) in [26]. The formula (168) is valid for light sources at finite distance. In the limit of infinite spatial distances, R → ∞, the relation (168) changes into the expression in (166). In summary of this section, the expression for the time delay in (133) and for the unit tangent vector in (168) are valid for a light-signal which has been emitted by a source located at finite spatial distances, and which propagates through the Solar system, that means through the gravitational field of N arbitrarily moving bodies and having arbitrary shape and inner structure and which can be in arbitrary rotational motion.
If the light-source is located at infinity, i.e. R → ∞ and in a good approximation realized by stars or quasars, then the light deflection angle of a light-signal in the field of N arbitrarily moving bodies in 1.5PN approximation is determined by where n 1.5PN is given by Eq. (166). If the light-source is located at finite distance, i.e. R is finite and in a good approximation realized by Solar system objects, then the light deflection angle of a light-signal in the field of N arbitrarily moving bodies in 1.5PN approximation is defined by where n 1.5PN is given by Eq. (168). The relation (169) is of simpler structure than (170), but which equation can be utilized depends on how far the light-source is. For our preliminary considerations here it will be sufficient to consider light-source at infinity, that means to apply just relation (169). Like in case of Shapiro delay, we will consider the light deflection caused by N arbitrarily moving bodies in some more detail, by considering bodies with mass-multipole structure and spin-dipole-structure. An estimate is also given for spin-octupole. The terms which are of relevance for nas-accuracy read In what follows we will consider these terms in some detail and give some estimates of their magnitude.
A. Light deflection for moving mass-monopole

In terms of coordinate time
From (110) and (111) we obtain for the coordinate velocity of the photon in the field of N arbitrarily moving monopoles: In the limit of uniformly moving bodies our result in (172) -(173) agrees with Eq. (6.3) and (6.5) in [24] and with Eq. (20) in [37] up to terms of the order O v 2 A /c 2 . By inserting (172) -(173) into (165) and then into (169) we obtain the light deflection angle, which for one massive body A can be estimated as follows: Like in Eqs. (142) and (162), we encounter again the typical large term in (175) which is proportional to ∼ r 1 A /d 1 A and originates from the last two terms in (173). This large term is solely caused by the retardation of gravitational action. That means, the use of the timemoment of reception at the body's position, x A (t 1 ) in Eq. (172) causes a significant error in the determination of light deflection for moving bodies. This peculiarity has been recognized long time ago, for instance see text below Eq. (34) in [28]. Especially, this issue has thoroughly and comprehensively been solved for moving pointlike bodies in the investigations [28,33,35,38]. In the next subsection we will further elucidate this fact.

In terms of retarded time
From the physical point of view, it is obvious that instead of t 1 one has to use the retarded time-moment for the position of the massive body in (172). That means, with the aid of relations (144) -(145) one may show that Eqs. (172) -(173) can be rewritten as follows: .
The last two terms in (173) do not explicitly appear in (177), because they are absorbed in (176). Accordingly, instead of (175) we obtain the following correct estimates for the 1PN and 1.5PN corrections in (176) and (177), respectively: The given upper limit in (179) agrees with Eq. (42) and (46) in [73], and with the results in [28]. For numerical values of the upper bound in Eq. (178) and (179) see Table III.
B. Light deflection for moving spin-dipole

In terms of coordinate time
The coordinate velocity of a light-signal propagating in the field of arbitrarily moving spin-dipoles can either be obtained from (113) using (124), or simply by timedifferentiation of Eq. (128), and reads: One may verify that in the limit of bodies at rest our result agrees with Eq. (59) in [36]. An upper bound for the magnitude of the light deflection is given by in agreement with the estimate given by Eq. (65) in [36] for lightrays which propagate in the equatorial plane of a rotating body at rest.

In terms of retarded time
Following the same arguments as in the above considerations, we may replace all expression in (180) by their retarded expressions according to Eqs. (154) -(156) for n = 1. Then the estimate of light deflections yields:    Table I. The given numerical values are determined for grazing lightrays, that means the impact parameter equals the radius of the massive body: dA = PA. For the light deflection in the field of spin-octupole, ϕ SO 1.5PN , we take the results of Ref. [64] where the light deflection in the field of one rotating body at rest and having constant mass density has been determined. Blank entries indicate that the effect is smaller than 1 nas. In view of the fact that astrometry on sub-µas-level implies an accuracy for ϕ at least better than 0.1 µas, the 1.5PN effects in light deflection become detectable within the very next generation of high-precision space-based astrometry missions.
C. Light deflection for moving mass-quadrupole

In terms of coordinate time
The 1PN correction to the coordinate velocity of the lightray in the field N arbitrarily-moving bodies with quadrupole structure, ∆ẋ J2 1PN (t), has already been given Eq. (117) in [26] and can also be deduced from Eq. (110), while the 1.5PN correction ∆ẋ J2 1.5PN (t) from Eq. (111). In view of the complexity of these terms, we will not present these expressions in their explicit form. We just mention that the estimation of these terms proceeds sim-ilar to the procedure performed in [79]. After some considerable amount of algebra one obtains: The estimate of the 1PN quadrupole term in (183) is equal to the much simpler case of quadrupoles at rest, cf. Eq. (41) in [36] and Eq. (13) in [79]. The second term in (184) is proportional to ∼ r 1 A /d 1 A which for grazing rays becomes large. Like in Eqs. (142), (162), and (175), this term is caused by the finite speed of gravitational action.

In terms of retarded time
One may rewrite the expression for ∆ẋ J2 1PN (t) and ∆ẋ J2 1.5PN (t) in terms of retarded time by means of Eqs. (144) -(145), and formally one may also replace M ab (t 1 ) → M ab (t ret 1 ). Then the estimation of the 1PN and 1.5PN correction terms in the quadrupole light deflection becomes which agrees with Eqs. (44) and (46) in [73]. The numerical magnitude of these upper bounds (185) and (186) can be found in Table III.

D. Light deflection for higher mass-multipoles
The 1PN solution (110) and the 1.5PN solution (111) for moving bodies with full mass-multipole structure allow to determine the light deflection in the field of moving mass-multipoles to any order in l. However, the expressions for ∆ẋ Jn 1PN and ∆ẋ Jn 1.5PN (J n are the actual zonal harmonic coefficients of the massive body) become more and more involved the higher the order of the massmultipoles are and imply a considerable amount of algebra. The investigation of these terms will be postponed for awhile. In meanwhile let us consider an educated guess that the light deflection in the field of higher massmultipoles is determined by the following relation: which in case of n = 2 agrees with Eqs. (185) -(186). The suggestion in Eqs. (187) -(188) is based on the considerations above and triggered by the fact that in the limit of bodies at rest formula (187) agrees with the results in [41]. Numerical values for (187) are presented in Table III, while (188) yields values below 1 nas for n ≥ 3. A detailed proof of (187) and (188) and a comparison of formula (188) with [42] will be presented in a subsequent investigation.

X. SUMMARY AND OUTLOOK
During the last 25 years, astrometric measurements have made an impressive advancement from milliarcsecond level of accuracy by the ESA astrometry mission Hipparchos [1, 2] toward micro-arcsecond level of accuracy by the ESA astrometry mission Gaia [5]. Ever since, applied relativity has evolved into one of the basic components of modern astrometry, the branch of science which includes the whole machinery of advanced astrometric measurements, especially: (1) theory of reference systems, (2) precise description of light trajectory from the celestial light source toward the observer, (3) relativistic modeling of real observations, (4) determination of the metric of the Solar system in post-Newtonian approximation (weak-field slow-motion approximation) or post-Minkowskian approximation (weak-field approximation) and beyond, (5) multipole expansion of metric tensor of Solar system, (6) relativistic data reduction of astrometric measurements, and (7) determination of ephemeris of the Solar system bodies and of the observer accurate enough for a given accuracy.
But for all that stunning progress, the step from microarcsecond toward nano-arsecond astrometry will be a long-term ambition, which implies many challenges on theoretical as well as technological side. While a few of these issues have been mentioned in the introductory section, most of these challenges and especially their elaborated details cannot be foreseen at present. But for any actual ambitions about sub-micro-arcsecond astrometry two of these problems are of decisive importance: first to establish a set of accurate reference systems and reference frames for exact data reduction, and second to provide an accurate modeling of light trajectory from the celestial light source through the Solar system toward the observer. As it has been mentioned in the introductory section, especially these two highly important issues have also been emphasized by the ESA-Senior-Survey-Committee (SSC) in response of the selection of science themes for future space-based astrometry missions [16]. The presented investigation is mainly devoted to these two specific subjects. Especially, in order to arrive at a precise modeling of light-propagation through the Solar system, two difficult aspects have carefully to be treated: (1) First, in compliance with the requirements of the IAU recommendations [19,20], one has to introduce one global reference system (BCRS) and N local reference systems (GCRS-like), one for each massive body, which allow to describe the global metric of the Solar system in terms of intrinsic mass-multipoles and intrinsic spinmultipoles the massive bodies, that means for the metric perturbations h αβ M A L , S A L , as mentioned by Eq. (15). (2) Second, for sub-micro-arcsecond or even nanoarcsecond-astrometry one has to describe the light trajectory in the field of arbitrarily moving massive bodies, that means as a function of their worldlines x A (t), because a series expansion like in Eq. (16) is unsuitable for several reasons discussed in the introductory section. The worldlines can be concretized by Solar system ephemeris [63] at any stage of the calculations.
In a previous investigation [26] we have obtained a solution in 1PN approximation for the light trajectory through the Solar system in full agreement with these both requirements (1) and (2). As outlined in more detail in [26] and also mentioned in the introductory section, for high-precision astrometry on sub-µas-level or nas-level of accuracy the 1PN approximation is not sufficient at all. Instead, it is inevitable to determine the light trajectory through the Solar system in 1.5PN approximation and to reconcile the entire approach with the important requirements (1) and (2). Such an approach has been developed here in the presented investigation. Accordingly, the main results of our investigation are given by the first integration of geodesic equation in Eq. (104) and by the second integration of geodesic equation in Eq. (114): where the time-argument τ + t * have been omitted here for simpler notation. The terms in (189) for one body A are given by Eqs. • time-derivative of (104) yields (95).
• our results agree with [31] for bodies at rest and time-independent mass-multipoles in 1PN approximation.
• our results agree with [31] for bodies at rest and time-independent spin-multipoles in 1.5PN approximation.
• our results agree with [38] for arbitrarily moving bodies with spin-dipole in 1.5PN approximation.
• our results agree with [36] for bodies at rest with spin-dipole.
Further cross-checks in 1PN approximation have already been done in [26] for the case of light-propagation in the field of bodies with mass-monopole, mass-dipole, mass-quadrupole structure and bodies at rest with full mass-multipole structure.
The numerical magnitude about the impact of massmultipoles and spin-multipoles on light deflection, presented in Table III, reveal that the first mass-multipoles up to order l = 10 and the first spin-multipoles up to order l = 3 have to be taken into account for astrometry on nano-arcsecond level of accuracy. This fact is important in view of the complexity of the 1.5PN solution for the light trajectory, because it allows to simplify that solution considerably. However, more detailed investigations are very necessary in order to simplify the massive computations in astrometric data reduction as much as possible.
The approach presented has further to be developed into several directions before the conditions are complied for a complete modeling of light-propagation through the Solar system on sub-µas or nas-level of accuracy. In particular, the following issues may serve as minimal supplement to the list of aspects which have already been mentioned in the introductory section: A. The model for the light trajectory has to implement some terms in 2PN approximation, which can formally be written as follows: where (191) and (192) represents the coordinate velocity and the trajectory of the light-signal, respectively. The 2PN corrections have been determined for the case of monopoles at rest [18,27] and later recalculated in progressing investigations in [24,45,46] and also within this work, see appendix J. It is clear that for a comprehensive theory of light propagation aiming at sub-µas-level of accuracy it needs carefully to be scrutinized which 2PN corrections beyond the monopole part are of relevance for such extremely-precise astrometry. B. A fundamental prerequisite in order to gain further progress in the theory of light propagation in 2PN approximation, one necessarily needs to determine the space-space part of the BCRS as well as of the GCRS metric tensor including all terms of the order O c −4 . However, an extension of these global and local reference systems to the post-post-Newtonian order is a highly involved assignment of a task and is presently an active field of research [80][81][82][83] and far from being completed.
C. In the first instance, the post-Newtonian approach of the DSX formalism allows for astrometry in the nearzone of the Solar system. However, astrometric measurements of stars or extragalactic celestial objects are subject to far-zone astrometry, which requires a matching procedure of two asymptotic solutions: the near-zone solution and the far-zone solution for the light trajectory [21]. Such matching approach has been proposed in [24,25], which has to be further developed in such a way to be in line with the requirements of nas-astrometry.
D. The unique interpretation of observational data implies a hierarchy of several reference systems [18,21]: a) BCRS x 0 , x 1 , x 2 , x 3 for description of the light trajectory in the Solar system, b) GCRS-like X 0 A , X 1 A , X 2 A , X 3 A , one for each body A = 1, ..., N of the Solar system in order to define the intrinsic multipoles, c) CoMRS X 0 , X 1 , X 2 , X 3 which is co-moving with the observer, d) ToRS z 0 a , z 1 a , z 2 a , z 3 a , one for each ground-station a = 1, ..., n on Earth which are involved in data reduction, where CoMRS stands for co-moving reference system and ToRS denotes topocentric reference system. The light trajectory in our investigation is given in the BCRS, but that is of course not sufficient for a comprehensive astrometric model of light propagation. In particular, the presented solution has to be transformed into the reference system which is co-moving with a free-falling observer (CoMRS) [84]. This transformation takes account for aberrational effects. Especially, it has to be clarified whether or not the CoMRS in [84], which was primarily intended for the Gaia mission, is also sufficient for the requirements on nas-level of accuracy.
E. The basic assumption of post-Newtonian expansion is that all retardations of the gravitational actions are small. In the model presented the effect of retardation has been implemented in a more or less heuristic manner, in order to provide a proper estimation for the upper limit of time delay and light deflection. This procedure needs to be scrutinized in considerably more detail. Especially, it has to be clarified how the retardation of gravitational action has to be implemented based on clear theoretical foundation in the entire approach. The solution of this problem is related to the far-zone astrometry about how the presented solution in the near-zone can be matched with the solution for the lightray in the far-zone of the Solar system [21,24].
In summary, a precise determination of light trajectory up to a given accuracy is of fundamental importance in the theory of any astrometric measurements. Besides considerable effort which has still to be done in near future, we come to the conclusion that a complete modeling of light trajectory from celestial light sources through the Solar system toward the observer is accomplishable also for extremely high-precision astrometry on sub-µas and even on nano-arcsecond level of accuracy.

XI. ACKNOWLEDGMENT
This work was supported by the Deutsche Forschungsgemeinschaft (DFG).

Appendix A: Notations
Throughout the article the following notations are in use: • G is the Newtonian constant of gravitation.
• c is the vacuum speed of light in flat Minkowski space.
• The three-dimensional coordinate quantities ("three-vectors") referred to the spatial axes of the corresponding reference system are set in boldface: a.
• The contravariant components of "three-vectors" are a i = a 1 , a 2 , a 3 .
• Repeated indices imply the Einstein's summation irrespective of their positions (e.g. • The absolute value (Euclidean norm) of a "threevector" a is denoted as |a| or, simply, a and can be computed as a = |a| = (a 1 a 1 + a 2 a 2 + a 3 a 3 ) 1/2 .
• The scalar product of any two "three-vectors" a and b with respect to the Euclidean metric δ ij is denoted by a · b and can be computed as a · b = δ ij a i b j = a i b i .
• The vector product of any two "three-vectors" a and b is designated by a × b and can be computed is the fully antisymmetric Levi-Civita symbol.
• The global coordinate system is denoted by lowercase letters: (ct, x).
• The local coordinate system of a massive body A is denoted by upper-case letters: (cT A , X A ).
• The photon trajectory is denoted by x (t). In order to distinguish the photon's spatial coordinate x (t) from the spatial coordinate x of the global system, the time-dependence of photon's spatial coordinate will everywhere be shown explicitly throughout the article.
• L = i 1 i 2 ...i l is a Cartesian multi-index of a given tensor T , that means T L ≡ T i1i2 . . . i l , and each index i 1 , i 2 , ..., i l runs from 1 to 3 (i.e. over the Cartesian coordinate label).
• Two identical multi-indices imply summation, e.g.: • The symmetric tracefree (STF) part of a tensor T L is defined by Eq. (A2) in [26] and denoted by T L .

Appendix B: Notation of impact vectors:
Before we distinguish between the case of massive bodies at rest and massive bodies in motion, we consider the unperturbed lightray in flat Minkowskian space-time, which in Cartesian coordinates is given by the expression in (1), which describes a straight line and where the subscript N stands for Newtonian limit. By Eq. (78) we have introduced the following impact vector: The impact vector in (B2) points from the origin of the global system (BCRS) toward the point of closest approach of the unperturbed lightray to that origin. The impact vector in (B2) is time-independent, both in case of massive bodies at rest as well as in case of massive bodies in motion.

Massive bodies at rest:
Massive bodies at rest means their positions remain constant with respect to the global reference system: x A = const. We will make use of the following notation for the vector from the massive body at rest toward the photon propagating along the exact light trajectory: with the absolute value r A = |r A |. The vector from the massive body at rest toward the photon along the unperturbed light trajectory reads: with the absolute value r N A = r N A , and obviously r A = r N A + O c −2 . We also need the vector from the massive body at rest toward the photon at the moment of signalemission: with the absolute value r 0 A = r 0 A . Note that in case of massive bodies at rest there will be no time-argument in r A and r N A , irrespective of the fact that the distance between the photon and the body actually depends on time due to the propagation of the photon. In case of massive bodies at rest we introduce the following impactvector: The impact-vector in (B7) is time-independent,ḋ A = 0, and points from the origin of local coordinate system of massive body A toward the unperturbed lightray at the time of closest approach to that origin, defined by Notice that the term weak gravitational field implies

Massive bodies in motion:
In case of massive bodies in motion, their positions become time-dependent: x A (t). Then we will make use of the following notation for the vector from the massive body toward the photon propagating along the exact light trajectory: with the absolute value r A (t) = |r A (t)|. The vector from the massive body in motion toward the photon along the unperturbed light trajectory reads: with the absolute value r N A (t) = r N A (t) and obviously r A (t) = r N A (t) + O c −2 . We also will need the vector from the massive body toward the photon at the timemoment of emission of the light-signal, given by with the absolute value r N A (t 0 ) = r N A (t 0 ) . In case of massive bodies in motion we introduce the following impact vector: with the absolute value d A (t) = |d A (t)|. The impactvector in (B13) is time-dependent,ḋ A = 0, and points from the origin of local coordinate system of massive body A toward the unperturbed lightray at the time of closest approach to that origin. The time-dependence of the impact-vector in (B13) is solely caused by the motion of the massive body, that means a time-derivative of (B13) is proportional to the orbital velocity of this body, G M A c 2 for the time of closest approach of the lightray to the massive body, which are given by and which slightly differ from the expressions (B8) and (B9) by the time-argument of the spatial coordinates of the massive body.
Appendix C: Partial derivative operator The spatial derivative in terms of the new variables τ and ξ has been given by relation (86) which is valid for any smooth function F (t, x), that means According to the metric perturbations in (68) -(72) we have to consider the STF partial derivative operation in Eq. (73), which reads In order to express the spatial derivative operation in (C2) in terms of these new variables, we apply the binomial theorem: where the binomial coefficients are defined by In virtue of the binomial theorem in (C3), we obtain for the STF partial derivative operator in (C2) in terms of the new variables τ and ξ the following expression: (C5) The same expression for ∂ L has been used in [40] (cf. Eqs. (4.42) -(4.43) ibid.); note the symmetry p − q ↔ q of the expression in (C5). The derivatives with respect to variable c t * act only on M A L (τ + t * ) and x A (τ + t * ), hence the partial derivatives ∂ ∂c t * q produce terms of the order O (c −q ). For that reason it was possible to neglect all derivatives with respect to variable c t * in 1PN approximation which has been investigated in [26]. If one neglects such derivatives (i.e. take only the terms with q = 0 in Eq. (C5)), then we would obtain the simpler derivative operator: which coincides with the expression as given by Eq. (24) in [31] or by Eq. (101) in [26] where it was allowed to omit the STF operation because of relation (101).

Appendix D: Derivatives
In this appendix we will summarize some useful spatial-derivatives and time-derivatives. Throughout this appendix all time arguments are omitted in order to simplify the notations, that means

Spatial-derivatives
In this appendix some relevant relations for spatialderivatives are summarized. The vector r N A depends only on the variables ξ, τ and x A (τ + t * ). Since variable ξ is independent of τ and x A (τ + t * ), we consider partial derivatives with respect to variable ξ. We obtain the following relations: Notice that ∂ ∂ξ j σ · r N A = 0.

Time-derivatives
In this appendix some relations of time-derivatives are summarized which are of relevance for integrations by part. Since the position of massive body depends on time-variable, x A (τ + t * ), this vector is not independent of τ . Therefore, in order to perform integration by parts, we need to consider the derivatives with respect to variable τ . Taking into account r N A = d A + σ σ · r N A , we find the following relations: (D10) The first integration of the expression in the first line in (96) reads: where ∂ L is given by Eq. (C5) where cτ is formally replaced by the integration variable cτ . In order to solve that integral, we use relation (D9) and obtain where we also have used that r N A = d A +σ σ · r N A . Note that Eq. (E2) is an exact expression for the integral in Eq. (E1). Now the expression in the first line in (E2) will be integrated by parts. For the integral in the second line in (E2) we use relation (D9) again, while for the integral in the third line in (E2) we will use relation (D12) and obtain by means of integration by parts: By means of relation (C5) in [26] one may show that the lower integration limit τ → −∞ in the first line of Eq. (E2) vanishes. In order to determine the integral in the second line in (E3), we use relation (D11), and obtain finally: The integration of the fifth term in Eq. (96) reads as follows: Finally, by means of relation (D5), we just obtain the expression in the last line of Eq. (111). .

(G2)
In order to get the expression in (G2), we have performed an integration by parts which results in the expressions in the first and second line. Furthermore, for the expression in the third line we have used relation (D6). Consequently, (G2) represents an exact expression for the integral in (G1). Now we are going to proceed with the consideration of the remaining three integrals in the third, fourth, and fifth line in (G2). For the integral in the third and fifth line we will use relation (D13) and integrate by parts; note that we also need relation (D9) and the facts that For the integral in the fourth line we use relation (D11) and integrate by parts and afterwards we apply relation (D6). Altogether, we obtain: where in the first line we have used relation (D6).

Integral I4
The integral I 4 reads .

Integral I5
The integral I 5 reads .

Integral I7
The integral I 7 reads .

Integral I10
The integral I 10 , using relation (D5), reads Integration by parts using relation (D11) and (G4) yields, with the aid of relation (D6): I 10 (τ + t * , τ 0 + t * ) = I 10 (τ + t * ) − I 10 (τ 0 + t * ) , where where S γ ( α u β ) = S γ α u β + S γ β u α / 2 means the symmetrization with respect to the indices α and β. Thereby, u β = γ v (c , v) where γ −1 v = 1 − v 2 /c 2 is the Lorentz factor, and all time-dependent quantities depend on the retarded time-variable τ ret , that means for the global spin-tensor S αβ = S αβ (τ ret ), for the four-velocity u α = u α (τ ret ) and for the three-velocity v = v (τ ret ). Furthermore, r α = (−r , r) with r = |r| and r being the vector pointing from the spatial position of the body at retarded time, x A (t ret ), toward the spatial position of the photon at global coordinate-time, x (t). That means, in (I3) we may replace the new variables ξ, τ ret by the old variables x, t ret (see also text below Eq. (128)): where the retarded time in terms of the old variables is given by Eq. (143); see also Eq. (11) in [33] or Eqs. (12) in [38]. The solution in (I3) is valid for the light trajectory in post-Minkowskian approximation in the gravitational field of an arbitrarily moving point-like body carrying a spin-dipole, while our result in Eq. (122) and Eq. (128) is valid for the light trajectory in post-Newtonian approximation in the gravitational field of an arbitrarily moving extended body carrying a spin-dipole. In order to compare both results we have to expand all expressions in (I3) with respect to variable v/c 1 and neglect all terms of the order O c −4 , and afterwards we have to express the global spin-tensor S αβ in terms of the intrinsic spin-vector S. Especially, we find where in (I5) and (I6) we have attached an index A in order to indicate that actually the body A is meant here, while the spin-tensor in (I7) describes still the global spin, besides the fact that this spin-tensor originates from the intrinsic spin S A of that single body. Let us consider one specific example by performing a series-expansion of the first term in (I3), for which we obtain: where in the second line we have determined the symmetrization, while in the third line we have used (I5) -(I7) and P ij r j = d i A (t). Very similar steps for the other terms in (I3) yield the following expression: where we have used σ α σ γ S γα = 0 because the spintensor is anti-symmetric, S γα = −S αγ , and we have introduced r N γ (t) = −r N A (t) , r N A (t) . In order to compare (I9) with our result in (128) we have to express the global spin-tensor S αβ in (I9) in terms of the intrinsic spinvector S A , where the index refers to body A. Recalling relations (24) and (C.10) in [74] we have By inserting (I10) -(I11) into (I9) we arrive at Finally, using r N A (t) = d A (t) + σ σ · r N A (t) we obtain where the last two terms in (I13) comprise the last two terms in (I12). The expression in (I13) agrees with our result in (128); note that all derivatives according to (124) have been performed, hence the replacement τ + t * → t in (128) is possible.
where r N A = x N (t) − x A and r 0 A = x (t 0 ) − x A . The expression in (J7) corrects some typos in Eq. (B6) in [26].

Light trajectory in 2PN approximation for monopole at rest
The light trajectory in post-post-Newtonian approximation (2PN) reads: where r 1PN A = x 1PN (t) − x A . The vectorial functions B 1 and B 2 are given by (cf. Eqs. (3.2.41) and (3.2.42) in [18] or Eqs. (50) and (51) in [45]): It should be mentioned that in B 1 the coordinate of the photon in 1PN approximation, x 1PN , can be replaced by the exact coordinate x of the photon, and in B 2 the coordinate of the photon in Newtonian approximation, x N , can be replaced by the exact coordinate x, because such replacements are correct up to terms of the order O c −6 .
Sometimes it is useful to perform a series-expansion of the vectorial function B 1 in terms of the small parameter m A and to express the 2PN solution (J8) in terms of unperturbed lightray as follows: where we recall that r N A = d A +σ σ · r N A and r 0 A = d A + σ σ · r 0 A . The expression in (J11) corrects some typos in Eq. (37) in [26]. Let us notice that in the second line in (J11) it is not allowed to replace the coordinate of the photon in Newtonian approximation, x N , by the exact coordinate of the photon, x, because such a relacement would cause an error of the order O c −4 . This is the reason for the fact that the form in (J11) is usually not in use in favor of the expression in Eqs. (J8) -(J10).