#### 1. 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 *B* – mode 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–17], TAMA [18–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 Flanagan and Hughes [1]. The space-based LISA mission^{1} [24–26], now evolved into eLISA^{2} [27], aimed to detect gravitational waves in the frequency range 10^{−5}−1 Hz, while accurate timing measurements of pulsars may detect signals in the range 10^{−9}−10^{−7} Hz [28–30]. Section 6.2 of Flanagan and Hughes [1] yields an overview of typical LISA sources: they include, e.g., equal mass binaries, in which the member black holes are of roughly equal mass (~10^{5}−10^{8} *M*_{⊙}), and extreme mass ratio binaries made by white dwarfs, neutron stars and 10−100*M*_{⊙} black holes captured by 10^{5}−10^{7} *M*_{⊙} black holes, located at quite large distances. As far as the very low-frequency band is concerned (10^{−9}−10^{−7} Hz), refer Section 6.3 of Flanagan and Hughes [1]. It may be due to several unresolved coalescing massive black holes; binaries which are either too massive to emit in the LISA band, or else are in-spiralling and will finally merge in several centuries or millennia.

In this paper, we will analytically work out the long-term orbital variations of all the six standard Keplerian orbital elements of a solar system planet due to the action of an externally incident monochromatic plane gravitational wave in the green–black^{3} part of the spectrum, i.e. with frequency ${\nu}_{g}\ll {10}^{-7}-{10}^{-10}$ Hz. Such kinds of gravitational waves are important since they carry information about how galaxies and black holes co-evolved over the history of the universe [31, 32], the early universe and related exotic physical processes like, e.g., inflation and cosmic strings, and possible physics beyond the standard model of particles and fields [33–37]. In particular, a primordial background of ultra-low frequency stochastic gravitational waves with a characteristic spectral shape should be produced due to the coupling of the gravitational field with the exponential expansion driven by the inflation [34, 35, 38–41]. Concerning the possible existence of a background of gravitational waves dating back to the origin of the universe, see, e.g., Weber [42]; Wheeler [43]; Zel'dovich and Novikov [44]; Carr [45], and the discussion in Mashhoon et al. [46].

The effects of incident gravitational waves on the orbital motion of gravitationally bound systems were inspected by several authors with a variety of approaches and approximations pertaining various features of both the waves and the orbits themselves [46, 47–61]. The idea of using the solar system to try to detect a stochastic background of gravitational waves of wavelengths much larger than about 1 au was first suggested by Bertotti [47].

At first sight, the calculations presented here may be regarded just as an academic exercise with respect to empirical celestial mechanics, although such an allegation may sound somewhat bizarre in view of the large amount of more or less analogous studies existing in literature concerning all sort of putative modified models of gravity based on much less solid theoretical and/or empirical support with respect to a GTR prediction like gravitational waves. Actually, it is not so for a variety of reasons. Indeed, recent developments in the field of planetary orbital determination, pursued by independent teams of astronomers [62–63], provided us with an increasing number of empirically estimated corrections $\text{\Delta}\dot{\text{\Psi}}$ to the standard Newtonian long-term precessions of some Keplerian orbital elements Ψ for almost all the major bodies of the solar system. The times when the focus was solely on the perihelion of Mercury are definitively waned. Such corrections $\text{\Delta}\dot{\text{\Psi}}$ are, in general, determined in a purely phenomenological way, and may account for any unmodeled/mismodeled dynamical feature not included in the usual dynamical force models fit to the observations. The availability of increasingly extended data records of higher quality, and the forthcoming adoption of more accurate observational techniques [64, 65] will allow to reach unprecedented accuracies in knowing $\text{\Delta}\dot{\text{\Psi}}$. This would allow to effectively put more and more stringent constraints on passing ultra-low frequency waves. This aspect is particularly important after the recent BICEP2 detection of B-mode polarization of the CMB at large angular scales [7]. Indeed, it is of crucial importance to look for independent tests of the gravitational waves causing it. Moreover, the availability of exact analytical expressions for the wave-induced long-term orbital effects may allow to set suitable linear combinations of corrections $\text{\Delta}\dot{\text{\Psi}}$ in order to separate them from unwanted, confusing orbital changes caused by other subtle standard physical effects like, e.g., solar and planetary oblateness, tides, minor bodies, etc. In this respect, it is also remarkable the fact that we have at our disposal empirically determined corrections $\text{\Delta}\dot{\text{\Psi}}$ for more than one planet. Moreover, as far as the search for gravitational waves in the solar system is concerned, our results may be extended also to spacecraft-based missions in the solar system like ASTROD-GW [66, 67], constituted of a number of probes in wide heliocentric orbits, and LISA [68, 69] envisaging the use of three spacecrafts orbiting the Sun at 1 au distance in a quasi-equilateral triangle formation 20° behind the Earth. It has to be noticed that, in this cases, quite accurate devices to measure relative positions among the spacecraft would be used. It is also the case of recalling that the possibilities offered by Doppler tracking of interplanetary drag-free spacecraft to detect cosmological gravitational waves with wavelength of the order of, or larger than 1 au were studied in the past [70–72]; for a recent review, see Armstrong [73] and references therein. Our calculations are valid, in principle, also for other natural systems like extrasolar planets^{4} [74] many of which have orbital frequencies of the order of ${10}^{-4}$ Hz, orbiting their parent stars at distances as small as ≃ 0.01 au. In such cases, our findings are 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 $\widehat{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.

#### 2. 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 $\{x,y,z,t\}$, 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}_{\hspace{0.17em}0j0}^{i},\hspace{0.17em}i,j=1,2,3$ of the local Riemann curvature tensor **R** as [76].

**K**which are, dimensionally, $[{\mathcal{K}}_{ij}]={\text{T}}^{-2}$.

In the linearized weak-field and slow-motion approximations, and for the case of a propagating gravitational wave, one has:

^{6}

*h*

_{00}and

*h*

_{0i},

*i*= 1, 2, 3 directly in the transverse traceless

^{7}(TT) gauge [76] finding that they vanish. At this point, a clarification about a subtle issue pertaining the form of eq. (1) is in order. If, on the one hand, the tidal matrix elements ${\mathcal{K}}_{ij}$ are gauge-invariant, on the other hand, the matrix product yielding eq. (1) is gauge-dependent. In GTR, a gauge transformation brings with it a coordinate transformation. In general, the Fermi coordinates are different from the TT coordinates. Suffice it to say that, in the TT gauge, it is

*h*

_{00}= 0, thus eliminating the Newtonian potential and making impossible any description of orbital mechanics because of lacking of the central body. Nonetheless, as a general remark by Baskaran and Grishchuk [79], even if one starts from the general form of the wave's field, which includes also the non-TT components, one would end up with equations of motion containing only the TT components. It is so because the equation for the geodesic deviation involves the curvature tensor (and its derivatives) in which the non-TT components automatically cancel out. Moreover, the transformation between the Fermi and the TT coordinates [79] is proportional to the ratio of the binary's spatial extension to the wavelength and to its time derivative, which in our case are negligible. This would allow to neglect also the “magnetic”-type components entering the equations of motion in addition to the “electric” ones yielding just eq. (1).

The tidal matrix **K** is not only symmetric but also traceless. It has five independent components, so that the acceleration $A$ of eq. (1) felt by a test particle becomes, quite generally:

For a wave not traveling in the reference $\{x,y\}$ plane, i.e. for ${\widehat{k}}_{z}\ne 0$, eq. (5) yields

When the wave propagates in the reference $\{x,y\}$ plane, i.e., for ${\widehat{k}}_{z}=0$, the second equation in eq. (6) becomes:

*x*axis, i.e. for ${\widehat{k}}_{y}=0$.

Finally, let us notice that when $\widehat{k}=\{\pm 1,0,0\}$, eq. (5) tells us that

so that it can be posed:In order to make contact with realistic situations occurring in typical solar system data analyses, we remark that the unit vector $\widehat{k}$ can be written, in general, as:

*β*and longitude

*λ*, which are to be considered as unknown: in this case, the reference $\{x,y\}$ 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

*β*and

*λ*through the angle Θ [59]. It is defined from $\widehat{k}\cdot \widehat{N}=\mathrm{cos}\hspace{0.17em}\Theta $ and $\left|\widehat{k}\times \widehat{N}\right|=\mathrm{sin}\hspace{0.17em}\Theta $, where $\widehat{N}$ is the unit vector directed along the test particle's out-of-plane direction coinciding with the direction of the orbital angular momentum.

Inserting eq. (6) with eq. (13) in eq. (4) yields^{8}

*z*axis, i.e. for $\beta =\pm 90$°, eq. (5) and eq. (6) yield:

For $\beta =0$, i.e. for ${\widehat{k}}_{z}=0$, eq. (4) reduces to^{9}

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) ($\widehat{k}=\{0,0,\pm 1\}$), widely treated in literature.

#### 3. THE LONG-TERM VARIATIONS OF THE KEPLERIAN ORBITAL ELEMENTS

The typical planetary orbital frequencies *n*_{b} in the solar system vary from ${n}_{\text{b}}=1.3\times {10}^{-7}$ Hz (Mercury) to ${n}_{\text{b}}=1.2\times {10}^{-10}$ Hz (Pluto), corresponding to timescales ${P}_{\text{b}}$ ranging from $7\times {10}^{6}$ s (Mercury) to $8\times {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\times {10}^{2}$ s (Mercury) to $2\times {10}^{4}$ s (Pluto). Thus, if we consider the action of a monochromatic plane gravitational wave of frequency ${\nu}_{g}$ and wave vector **k**, with $k={\nu}_{g}/c$, on a planetary orbit during a time interval $\Delta t$ comparable to an orbital period *P*_{b}, its phase $\Phi \doteq {\nu}_{g}t-k\cdot r={\nu}_{g}[t-(r/c)\mathrm{cos}\hspace{0.17em}\alpha ]$ can be reasonably approximated by $\Phi \simeq {\nu}_{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 ${\nu}_{g}/{n}_{\text{b}}\ll 1$ as well.

##### 3.1. 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 [80] 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 [81, 82]; 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 $\widehat{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 $\left\{xy\right\}$ 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}_{\mathrm{max}}=a(1+e),{r}_{\mathrm{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^{10} 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^{11} $I,\Omega ,\omega $ which, actually, experience slow changes^{12} in time as far as the planets of the solar system are concerned. A further modulation is due to the coefficients ${\mathcal{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^{13} 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.

##### 3.2. 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

*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 ${\nu}_{g}^{2}$, 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)).

##### 3.3. 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 ${\nu}_{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 ${\nu}_{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 ultra-low frequency spectrum of gravitational waves for which neither ground-based nor space-based dedicated experimental devices are available.^{14} 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 [84–86], 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.

#### 4. 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.

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 [87]). Then, Rudenko [48] assumes a normally incident wave and solves the equations of motion at zero order in eccentricity. Finally, he studies the variation $\Delta r$ 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 work out the temporal changes over one orbital revolution of the semi-major axis and the eccentricity of a circular orbit by setting $\Omega =\varpi =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 $\Omega =I=\varpi =0$. In all of such cases, he takes different values for ${\nu}_{g}/{n}_{\text{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 ${\nu}_{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 $\Delta r/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 $\Delta \mathcal{R}/\mathcal{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 $\Delta r$ 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 ${\nu}_{g}/{n}_{\text{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\cdot 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 ${\nu}_{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 *ℳ* as independent variable. Thus, Ivashchenko [57] obtains approximate expressions of the form

*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 ${\nu}_{g}/{n}_{\text{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 $\{x,y\}$ plane, so that $k\hspace{0.17em}\cdot \hspace{0.17em}r=0$. Then, Kochkin and Sbytov [58] do not recur to any known perturbative schemes, like, e.g., the Gauss and the Lagrange equations [80], 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 ${\nu}_{g}/{n}_{\text{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 ${\nu}_{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,\ell ,g$ for a non-circular orbit. Fourier series expansions in terms of the mean anomaly are performed. Different values for ${\nu}_{g}/{n}_{\text{b}}$ and wave's polarizations are, then, considered.

Chicone and Mashhoon [88], 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 [88] 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 [88] 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 *ℛ* from them in order to use the perturbative scheme based on the Lagrange planetary equations [80] for all the Keplerian orbital elements, apart from *ℳ*. Then, they evaluate *ℛ* 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^{15} $a,e,I,\Omega ,\omega $ of Venus and Pluto for three, unspecified, different sources of gravitational waves.

#### 5. 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 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.