Orbital effects of a monochromatic plane gravitational wave with ultra-low frequency incident on a gravitationally bound two-body system

We analytically compute the long-term orbital variations of a test particle orbiting a central body acted upon by an incident monochromatic plane gravitational wave. We assume that the characteristic size of the perturbed two-body system is much smaller than the wavelength of the wave. Moreover, we also suppose that the wave's frequency is much smaller than the particle's orbital one. We make neither a priori assumptions about the direction of the wavevector nor on the orbital geometry of the planet. We find that, while the semi-major axis is left unaffected, the eccentricity, the inclination, the longitude of the ascending node, the longitude of pericenter and the mean anomaly undergo non-vanishing long-term changes. They are not secular trends because of the slow modulation introduced by the tidal matrix coefficients and by the orbital elements themselves. They could be useful to indepenedently constrain the ultra-low frequency waves which may have been indirectly detected in the BICEP2 experiment. Our calculation holds, in general, for any gravitationally bound two-body system whose characteristic frequency is much larger than the frequency of the external wave. It is also valid for a generic perturbation of tidal type with constant coefficients over timescales of the order of the orbital period of the perturbed particle.


INTRODUCTION
Gravitational waves [1,2] are a key theoretical prediction of the general theory of relativity (GTR). Indeed, since GTR relies upon the Lorentz invariance, which carries with it the concept of a limiting speed for physical interactions, the existence of gravitational waves is a natural consequence of it. A direct measurement of them is still lacking; see, e.g., Cerdonio [3]; Giazotto [4]; Fairhurst et al. [5] for recent reviews of the status of gravitational wave detection. To date, only indirect evidences of their existence have been inferred from the orbital decay rate of the binary pulsar PSR B1913+16 [6] and, more recently, from the detection [7] of Bmode polarization at degree angular scales in the Cosmic Microwave Background (CMB) by the ground-based experiment BICEP2 [8] at the South Pole. The consequences of detecting gravitational waves for physics, astrophysics, and cosmology would be certainly remarkable; see, e.g., Sathyaprakash and Schutz [9]. The role of a direct detection of the gravitational waves for GTR and extended theories of gravity was pointed out by Corda [10,11]. The gravitational wave spectrum covers an interval of about 18 orders of magnitude in wavelengths, encompassing a very broad range of physics and astrophysics [12]. The frequencies in the range 10 1 10 4 Hz are the targets of several ground-based detectors like, LIGO [13,14], VIRGO [15][16][17], TAMA [18][19][20], GEO [21,22], etc. Typical sources of such high-frequency waves are coalescing binary systems hosting stellar-sized compact objects like neutron stars and/or black holes; see, e.g., Section 6.1 of technically valid for waves with higher frequencies with respect to the solar system ones; indeed, they might be as high as about 5 10 6 Hz corresponding to the green-light blue part of the spectrum in Figure 1 of Prince et al. [12]. Another scenario to which our analysis can be applied is the stellar system orbiting the supermassive black hole (SBH) hosted by the galactic center (GC) in Sgr A* [75], where the orbital periods of the stars discovered so far are larger than 16 years corresponding to frequencies smaller than 2 Â 10 9 Hz. On the other hand, our results are not necessarily limited to the very-low frequency waves case, being valid for any tidal force with constant (over particle's characteristic timescales) matrix coefficients as well, independently of its physical origin. The plan of the paper is as follows. In Section 2, we will discuss the analytical form of the acceleration experienced by the test particle due to an incident monochromatic plane gravitational wave traveling along a generic direction b k of an arbitrary local Fermi frame. We will also consider the simplified cases of a gravitational wave moving along the z axis, as it is a choice widely adopted in literature, and in the reference {x, y} plane. In Section 3, we will analytically workout the long-period changes occurring in the particle's orbital motion in the case of wave's frequencies much smaller than the orbital ones. We will make neither any simplifying assumptions about the inclination and the eccentricity of the orbit nor on the reference frame adopted. In Section 4, we briefly review some of the approaches followed in literature. Section 5 is devoted to the conclusions.

THE ACCELERATION IMPARTED ON AN ORBITING PLANET BY A PASSING MONOCHROMATIC PLANE GRAVITATIONAL WAVE
The action of an incoming gravitational wave on a planet of the solar system is of tidal type with respect to a suitably constructed local inertial frame, represented by a Fermi coordinate system fx; y; z; tg, whose origin is located at the Sun's position. In general, a tidal acceleration A experienced by a slowly moving test particle due to an external curved spacetime metric can be written in terms of the "electric" components R i 0j0 ; i; j ¼ 1; 2; 3 of the local Riemann curvature tensor R as [75]. x where we introduced the coefficients K ij ¼ : À R i 0j0 ; i; j ¼ 1; 2; 3 of the tidal matrix K which are, dimensionally, ½K ij ¼ T À2 . In the linearized weak field and slow motion approximations, and for the case of a propagating gravitational wave, one has: In the case of a propagating plane wave, of the form 6 The fact that both 7 h 00 and h 0i , i = 1, 2, 3 vanish in the transverse traceless (TT) gauge [76], which is tacitly assumed here 8 , implies that the tidal matrix K is not only symmetric but also traceless. Thus, it has five independent components, so that the acceleration A of eq. (1) felt by a test particle becomes, quite generally: 5 However, for waves with such relatively high frequencies other facilities like LISA would be available: if they will finally become operative at the expected level of accuracy, they would likely surpass the possibilities offered by the extrasolar planets. 6 In eq. (3) v ln ; l; n ¼ 0; 1; 2; 3 represent the gravitational radiation field, while I ¼ : ffiffiffiffiffiffi ffi À1 p . 7 More general expressions for the TT metric tensor, not limited to the plane-wave and radiative approximation and with h 00 6 ¼ 0; h 0i 6 ¼ 0; i ¼ 1; 2; 3, can be found in Kopeikin et al. [77,78]; see also Poisson [2]. It can be seen that the parts not directly related to the wave are essentially identical to those used in the standard harmonic gauge used in ordinary planetary data reduction. 8 R i 0j0 is gauge-invariant in the linearized theory; thus, one can use any convenient coordinate system to evaluate it [1].
In fact, the condition that the force exerted on the particle by the wave is orthogonal to its direction of propagation yields three further constraints so that the independent components of K ij ; i; j ¼ 1; 2; 3 are just two. In the case of eq. (3), they are harmonic functions proportional to n 2 g . For a wave not traveling in the reference fx; yg plane, i.e. for b k z 6 ¼ 0, eq. (5) yields : Thus, we pose for the two independent polarizations of the gravitational wave. When the wave propagates in the reference fx; yg plane, i.e., for b k z ¼ 0, the second equation in eq. (6) becomes: while from eq. (5) it turns out: Thus, in this case, we define: Notice that both eqs. (8) and (9) do not hold for a wave traveling along the reference x axis, i.e. for b k y ¼ 0. Finally, let us notice that when b k ¼ fAE1; 0; 0g, eq. (5) tells us that so that it can be posed: SOR-ASTRO L. Iorio: Orbital effects of a monochromatic plane gravitational wave on a gravitationally bound two-body system In order to make contact with realistic situations occurring in typical solar system data analyses, we remark that the unit vector b k can be written, in general, as: in terms of the ecliptic latitude β and longitude λ, which are to be considered as unknown: in this case, the reference fx; yg plane would typically coincide with the mean ecliptic at the epoch J2000.0. Explicit expressions of the tidal matrix coefficients for a generic wave's incidence can be found in Chicone et al. [59]. In addition to the amplitudes of the two independent wave's polarizations and of their mutual constant phase difference, they contain b and k through the angle Θ [59].
Finally, for a wave traveling along the reference x axis eq. (4), with eq. (11) and eq. (12), becomes In Section 3, we will analytically workout the effects of eq. (4) on the trajectory of a test particle orbiting a central body with gravitational parameter GM, where G is the Newtonian gravitational constant and M is its mass, located at the origin of the Fermi frame traversed by a monochromatic plane gravitational wave. We will also consider the particular case of eq. (16) ( b k ¼ f0; 0; AE1g), widely treated in literature.

THE LONG-TERM VARIATIONS OF THE KEPLERIAN ORBITAL ELEMENTS
The typical planetary orbital frequencies n b in the solar system vary from n b ¼ 1:3 Â 10 À7 Hz (Mercury) to n b ¼ 1:2 Â 10 À10 Hz (Pluto), corresponding to timescales P b ranging from 7 Â 10 6 s (Mercury) to 8 Â 10 9 s (Pluto). They are much larger than the time required by the light to travel across the spatial extensions of the Sun's planetary orbits, ranging from 2 Â 10 2 s (Mercury) to 2 Â 10 4 s (Pluto). Thus, if we consider the action of a monochromatic plane gravitational wave of frequency n g and wave vector k, with k ¼ n g =c, on a planetary orbit during a time interval Dt comparable to an orbital period P b , its phase U ¼ : n g t À kÁr ¼ n g ½t À ðr=cÞ cos a can be reasonably approximated by U ' n g t, independently of the orientation α of k with respect to the planet's position r. As previously stated, in the rest of this section we will assume n g =n b 51 as well.

Monochromatic plane gravitational wave propagating along a generic direction
A straightforward first-order perturbative calculation performed with the standard Gauss equations for the variations of the osculating Keplerian orbital elements [79] yields the long term, i.e. averaged over one orbital period, variations of all the osculating Keplerian orbital elements of the test particle due to eq. (4). We briefly recall the usual computational procedure. A standard Keplerian ellipse is assumed as reference, unperturbed orbit; any small additional acceleration like, e.g., eq. (4) is considered as a perturbation; the Gauss equations are valid for any kind of perturbation, irrespectively of its physical origin. The disturbing acceleration is inserted into the right-hand-sides of the Gauss equations, which are evaluated onto the unperturbed orbit. Then, an average over one full orbital revolution of the test particle is performed to obtain the long-term orbital variations. We notice that the analytical expression of eq. (4) was obtained in the TT gauge, which is a particular case of the harmonic gauge: from the point of view of the pertubative calculation using the Newtonian two-body problem as zeroth order term, it does not pose problems since the latter one is gauge-invariant [1]. It clearly appears also in the explicit expression for, say, h 00 in the TT gauge when a central mass is present in the local frame traversed by the external wave [2,77,78]. In principle, one could also take different reference orbits already including 1PN Schwarzschild-like terms [80,81]; in such a case additional mixed, Schwarzschild-radiative terms would arise, but they would be negligible because of higher order in powers of c À1 . As a final remark to adequately put the problem in context, we recall that we are in a linearized, weak-field and slow-motion scenario. We are not dealing here with the final merging stage of a two-body system made by highly relativistic, self-gravitating compact objects whose orbits are shrinking because of the emission of gravitational waves from the system itself. The long-term orbital variations due to eq. (4) are, thus In order to deal with manageable expressions of general validity, we did not display in eq. (19) the explicit expressions of the tidal matrix coefficients in terms of any specific representation of b k. We stress that eq. (19) is valid for a generic reference frame. Moreover, it is exact both in the eccentricity e and in the inclination I to the reference fxyg plane in the sense that no a priori simplifying assumptions about the eccentricity and the inclination of the perturbed test particle's orbit were assumed in the calculation. It can be noticed that the semi-major axis a is not affected by the passage of a very slowly varying gravitational wave; the variation of the eccentricity is proportional to e itself, so that a circular orbit does not change its shape. It is important to remark that, since in the real world exactly circular orbits do not exist, the previous result do not mean that the shape of the orbit is left unaffected by a passing ultra-low plane gravitational wave. Indeed, as previously stated, e characterizes just the orbit's appearance, and its change in time implies that, e.g., the aphelion and perihelion distances r max ¼ að1 þ eÞ; r min ¼ að1 À eÞ vary even if a, which is responsible of the orbit's size, does not. In other words, the ellipse would remain 11 inscribed within the same imaginary circle, but it would get more or less elongated as time passes. Moreover, eq. (19) does not contain genuine secular effects because of the presence of 12 I; X; x which, actually, experience slow changes 13 in time as far as the planets of the solar system are concerned. A further modulation is due to the coefficients K ij of the tidal matrix containing the (ultra-low) harmonic variation of the impinging gravitational wave. Since all such frequencies are much smaller than the typical orbital ones for the Sun's planet, the terms containing them were kept fixed in the integration yielding eq. (19). Finally, we notice that the temporal signatures of all the long-term variations of eq. (19) are peculiar to the action of just 14 an incident ultra-low plane gravitational wave. This would be of great help in recognizing its passage if and when orbital variations matching the specific patterns of eq. (19) will be, actually, detected over multidecadal analyses.

Monochromatic plane gravitational wave propagating along the z axis
By using the same procedure it is possible to work out the long-term variations of the Keplerian orbital elements for a direction of incidence of the wave coinciding with the reference z axis. By using eq. (16) one gets ð19Þ À 10e 2 sin I cos 2! 2K xy cos 2 þ K yy À K xx À Á sin 2 Â Ã þ þ 5e 2 sin 2! 4 cos 2I K yz cos À K xz sin 11 Actually, also its orientation within its plane itself would change because of the variation of -. 12 Here ω denotes the argument of pericenter, while -¼ : X þ x is the argument of latitude. 13 They are mainly due to the mutual N-body classical perturbations: the rates of change can be found at http://ssd.jpl.nasa.gov/ txt/p_elem_t1.txt. 14 This is true if the standard Newtonian-Einsteinian laws of gravitation are involved. Also in this case, no a priori assumptions on e and I were made, so that the rates of eq. (20) are exact in this respect. Notice that for I = 0, i.e. for incidence of the gravitational wave normal to the orbital plane, the inclination is left unaffected. Moreover, also in this case, the elements h 1 and h 2 of the tidal matrix, which are (slowly) time-varying harmonic functions whose amplitudes are proportional to n 2 g , were assumed constant in the integrations over the planet's orbital period. Notice that, for a plane wave traveling along the z axis, i.e. for eq. (15), eq. (19) reduces just to eq. (20). It is straightforward to obtain the formulas valid also for the other specific directions of propagation examined in Section 2 by suitably specializing eq. (19) to such cases (cfr. eq. (8), eq. (9) and eq. (11)).

A comment on the approximation used for the harmonic wave functions
In obtaining eqs. (19) and (20) we kept the tidal matrix coefficients, which include the time-dependent harmonic functions of the wave, constant in the integrations over one orbital revolution of the test particle. This implies that we considered waves having frequencies n g much smaller than the orbital ones n. As we will see in Section 4, it is in contrast with the other approaches followed so far in literature, in which no similar assumptions on n g were made. Apart from the fact that, from a computational point of view, our approach allows for exact calculations in e and I, our choice is justified also from a phenomenological point of view. Indeed, applying our results to the solar system implies that we could, in principle, put constraints over a part of the ultralow frequency spectrum of gravitational waves for which neither ground-based nor space-based dedicated experimental devices are available 15 . Moreover, in view of continuous tracking of planets by means of ranging either directly to their surfaces or to ever more numerous orbiting spacecraft it will be possible, in principle, to obtain more and more accurate bounds because of increasing data records over the years. On the other hand, it is, after all, of little practical utility to perform cumbersome calculations of the wave-induced orbital effects involving frequencies as large as, or even larger than, the planet's ones since much more accurate dedicated experiments already exist covering such windows of the spectrum of gravitational waves. Indeed, the accuracy reachable in future planned interplanetary laser ranging facilities [83][84][85], of the order of about 1 À 10 cm, is not comparable with that of the latest gravitational wave detectors; it can be acceptable in order to put upper bounds when no other, more accurate devices exist, as for the ultra-low frequency waves.

CONFRONTATION WITH OTHER APPROACHES IN LITERATURE
In this Section we briefly outline the approaches followed by some other researchers, with particular care to the orbital effects caused by the traveling wave.   15 A possible exception may be the proposed Inflation Probe by NASA [82], dedicated to map the polarization of the Cosmic Microwave Background at 10 À16 Hz (see Figure 1 of Prince et al. [12]). Notice also that precision timing of millisecond pulsars may be used for gravitational waves in the range 10 À7 À 10 À9 Hz [28][29][30].
Rudenko [48] uses polar coordinates in the orbital plane and considers an orbit disturbed by a monochromatic plane gravitational wave traveling along a generic direction; the additional force resulting on the test particle is obtained within the linear approximation of general relativity (within the framework of the "theory of gravitation in plane space" by Zel'dovich and Novikov [86]). Then, Rudenko [48] assumes a normally incident wave and solves the equations of motion at zero order in eccentricity. Finally, he studies the variation Dr of the orbital radius for various values of the wave's frequency.
Turner [52] works in the TT gauge by deriving the geodesic equations of motion in cartesian coordinates of both the binary system's constituents. Then, he takes their differences obtaining the equations for their relative motion. It turns out that they are formally different from those obtained by other authors like, e.g., Mashhoon [49] from the geodesic deviation equation, but Turner [52] shows that they are, actually, equal up to an unphysical coordinate transformation. In the scenario by Turner [52], the monochromatic plane wave travels along the negative z direction. He uses the Gaussian perturbative scheme to workout the temporal changes over one orbital revolution of the semi-major axis and the eccentricity of a circular orbit by setting X ¼ -¼ I ¼ 0. Then, he discards the limitation I = 0, but not the other ones. Finally, Turner [52] considers non-circular orbits, and computes the variation of a to the lowest order in e for X ¼ I ¼ -¼ 0. In all of such cases, he takes different values for n g =n b .
Mashhoon et al. [46], relying upon the methods developed in Mashhoon [49], look at the effects that various kinds of gravitational radiation, among which monochromatic plane waves are considered as well, have on a Newtonian binary system. Let us recall that Mashhoon [49], in considering slow particle motions and weak waves having wavelengths larger than the system's orbital size, adopts the linearized Jacobi equation for the geodesic deviation equation. Mashhoon [49] writes down the equations of motion in cylindrical coordinates and solve them for various values of n g with respect to n b . As far as the monochromatic plane wave case is concerned, also Mashhoon et al. [46] assume that its wavelength is much larger than the size of the binary system; no simplifying assumptions on the wave's frequency are made. Then, Mashhoon et al. [46] work out the resulting relative change Dr=r occurring in the test particle's distance from the primary in the low-eccentricity approximation. Moreover, Mashhoon et al. [46] consider also the case of a pair of planets moving in the same plane along circular orbits, and calculate the wave-induced relative change DR=R in their mutual distance. Nelson and Chau [56] look in the TT gauge at a monochromatic plane gravitational wave traveling along the z axis, and having wavelength much larger than the size of the system considered. They write down the components of the wave-induced acceleration experienced by the test particle in cylindrical coordinates for generic shape and inclination of the orbit. Then, they specialize them to the case of normally incident wave, i.e. for an orbit with zero inclination. At this point, Nelson and Chau [56] solve the resulting perturbed equations of motion by evaluating the radial change Dr over one orbital revolution. In their computation, they resort to some approximations in e, and do not consider the harmonic wave functions constant over typical orbital timescales. The resulting changes for different values of n g =n b are inspected. Then, Nelson and Chau [56] examine the case of a generic inclination of the test particle's orbit by working out the shift along the wave's propagation, i.e. along the z axis. Finally, Nelson and Chau [56] perform some numerical analyses of the radial shifts for different values of the eccentricity and of the phases of the wave. Also Ivashchenko [57] considers the case of a monochromatic plane gravitational wave, in the TT gauge, traveling along the z axis. He assumes that its wavelength is much larger than the characteristic size of the perturbed system, so that he can neglect the term k Á r in the wave's phase. On the other hand, Ivashchenko [57] does not make any a priori assumptions about the magnitude of the wave's frequency n g with respect to the orbital one n b . After having obtained the exact radial, transverse and normal components of the perturbing acceleration from its expression in cartesian coordinates, Ivashchenko [57] makes use of the Gauss equations for the variation of all the Keplerian orbital elements, apart from the mean anomaly. At this point, after having set ω = 0, he makes an expansion of the right-hand sides of the Gauss equations, evaluated onto the unperturbed Keplerian ellipse, to first order in e by using the mean anomaly M as independent variable. Thus, Ivashchenko [57] obtains approximate expressions of the form dn dt ¼ FðtÞ; n ¼ a; e; I; X; x; ð21Þ which he straightforwardly integrates. In doing that, he does not consider h 1 and h 2 as constant by keeping, instead, their harmonic time dependence. Then, Ivashchenko [57] discusses various particular cases for different orbital geometries and, especially, for different values of n g =n b . Also Kochkin and Sbytov [58] adopt a wave in TT gauge traveling along the reference z axis. Moreover, they assume for it a circular polarization, and do not make any assumptions about the wave's frequency and its wavelength. Kochkin and Sbytov [58] use cylindrical coordinates to write down the spacetime line element for a linearized superposition of the Schwarzschild and wave fields, from which the particle's equations of motion are obtained. Concerning the orbital geometry, they put the perturbed orbit, not assumed a priori circular, in the reference fx; yg plane, so that k Á r ¼ 0. Then, Kochkin and Sbytov [58] do not recur to any known perturbative schemes, like, e.g., the Gauss and the Lagrange equations [79], in dealing with the resulting equations of motion. After some changes of variables, they integrate the resulting modified equations over one orbital revolution by inferring the cumulative changes of the semi-major axis, the eccentricity and the pericenter after one orbital period. Finally, Kochkin and Sbytov [58] consider how the results they obtained vary for different values of n g =n b ; they consider the circular case as well. Chicone et al. [59] and Chicone et al. [60] look at the action of a monochromatic plane gravitational wave traveling along the reference z axis and normally impinging on an unperturbed two-body system. The wave's wavelength is assumed to be much larger than the semi-major axis of the orbit, and the characteristic velocities within the system are much smaller than c. No assumptions are made a priori on the wave's frequency n g . The resulting relative acceleration is obtained from the geodesic deviation equation in cartesian coordinates. Then, Chicone et al. [59] and Chicone et al. [60] write down the Hamiltonian of the perturbed system in polar coordinates, and adopt the Delaunay orbital elements L; G; '; g for a non-circular orbit. Fourier series expansions in terms of the mean anomaly are performed. Different values for n g =n b and wave's polarizations are, then, considered. Chicone and Mashhoon [87], after having generalized the Jacobi equation by taking into account also the relative velocity of the geodesics, consider a monochromatic plane gravitational wave traveling along the reference x axis of the local Fermi quasi-inertial frame. Chicone and Mashhoon [87] write down the resulting equations of motion of the test particle in cartesian coordinates; they are rather involved because of the presence of the non-linear, velocity-dependent terms. Then, Chicone and Mashhoon [87] study the motion of the test particle along the x axis itself. The general wave-binary system scenario adopted by Ismaiel and Saad [61] is analogous to that by, e.g., Ivashchenko [57]. However, after having written the components of the tidal-type wave acceleration in cartesian coordinates, Ismaiel and Saad [61] obtain a potential function R from them in order to use the perturbative scheme based on the Lagrange planetary equations [79] for all the Keplerian orbital elements, apart from M. Then, they evaluate R onto the unperturbed Keplerian ellipse, without fixing ω, and expand it to some, unspecified, order in e by using the mean anomaly a independent variable. In such a way, they obtain that the right-hand sides of the Lagrange equations depend only on t, which allows Ismaiel and Saad [61] to integrate them over one orbital period. Notice that they do not keep the wave's harmonic function constant in the integration. Finally, they compute the long-term changes of 16 a; e; I; X; x of Venus and Pluto for three, unspecified, different sources of gravitational waves.

SUMMARY AND CONCLUSIONS
We analytically worked out the long-term variations of all the six osculating Keplerian orbital elements of a test particle orbiting a central body due to the perturbing tidal acceleration induced by the passage through the system of a monochromatic plane gravitational wave. We assumed that its frequency is much smaller than the orbital frequency of the test particle, so that the time-dependent harmonic functions of the wave were kept constant in the integrations performed over one orbital period. We considered a generic direction of the wavevector. All the osculating Keplerian orbital elements, apart from the semi-major axis, undergo long-term variations. In particular, the eccentricity changes, so that the overall appearance of the orbit is altered as well: for example, the perihelion and aphelion distances vary accordingly. Such long-term variations are linear combinations of the non-vanishing elements of the tidal matrix, containing the slowly time-varying harmonic terms due to the wave, with coefficients which are complicated functions of the eccentricity, the inclination, the node and the pericenter of the test particle. We did not make any a priori simplifying assumptions concerning the inclination and the eccentricity of the test particle's orbit; thus, we obtained exact results as far as such a point is concerned. Moreover, their validity is not restricted to any specific reference frame. The variation of the eccentricity is proportional to the eccentricity itself, so that circular orbits do not change their shape. In the case of incidence normal to the orbital plane, also the inclination is left unaffected. From a practical point of view, in the most general case one has six unknowns: the two angles of the wavevector, the frequency of the wave, its two independent amplitudes, and a phase lag. In principle, it is possible to constrain all of them by determining the secular variations of some Keplerian orbital elements for, e.g., several planets of the solar system. Actually, this seems to be the current trend in astronomical research since extensive planetary data reductions, in which the corrections to the standard precessions of the perihelia and the nodes of some major bodies of the solar system are estimated, are currently ongoing by independent teams of astronomers. Future implementation of planned or proposed accurate interplanetary ranging facilities may further enhance such a program. It should also be considered that the resulting bounds would represent independent tests of the ultra-low frequency primordial gravitational waves which may have been indirectly detected in recent data analyses of BICEP2. A systematic analysis of the currently available planetary data for various possible orientations of the wavevector may be the subject of further, dedicated investigations. Our results are quite general since they are valid not only for the solar system's planets, but also for any generic gravitationally bound two-body systems whose characteristic orbital frequencies are quite larger than the incident wave's frequency. Typical 16 It is unclear if the figures in Table III of Ismaiel and Saad [61] refer to the rates of change of the Keplerian orbital elements, i.e. if they are m s 1 and deg s 1 , or to their shifts over one orbital period, i.e. if they are m and deg. Also the units adopted in Table III of Ismaiel and Saad [61] are unspecified.
examples are extrasolar planetary systems and the stars orbiting the Galactic black hole. Moreover, our findings can be extended to any generic perturbing acceleration of tidal type with constant coefficients over the characteristic timescales of the perturbed system. Exotic modified models of gravity and/or various types of external cosmological curved backgrounds may do so.