1. Introduction

The tremendous penetrating capacity of gravitational waves is a reason why it is so difficult to detect them. Another reason is the very low intensities by which they are emitted. A great number of theoreticians and experimentalists work on gravitational radiation not only because of its important applications but also because of its fundamental role in science. In the domain of general relativity, the particular solution of the field equations leading to the wave solution is usually preferred in the domain of the post-Newtonian (or weak field) approximation. The reason is that the wave contributes to its own source, besides the effect of the relatively strong field of source itself; greatly complicate the detection of the wave close to its source. Hence it is better to study the simple properties of the waves by looking far into the wave zone where the fields are weak. The effects of incident plane gravitational waves on the orbital motion of gravitational bound two-body system were investigated by several authors with different approaches and various features of both orbits and the waves [1] - [9] . In this work we treat the effect of GW on the orbital motion by semi-analytical treatment and we determine the perturbation on the orbital elements using new approach differ than the previous works. We construct the Hamiltonian of GW and solve the canonical equations in terms of the Delaunay elements $\left(l,g,h,L,G,H\right)$ , using the numerical method of Runge-Kutta and algebraic manipulations programs MATHEMATICA V10. The resulting theory includes short-period as long-period for practical application on planet Jupiter.

2. Hamiltonian of GW

According to the weak field approximation of general relativity, there exist in empty space gravitational waves generated by different sources far in the universe, and propagated with velocity c. The field equations for empty space in the linear weak field approximation

${\square }^2{h}_{ik}=0$

Has the plane wave as the simplest solution. ${\square }^2$ is d’Alembertian operator is given by

${\square }^2=\frac{{\partial }^2}{\partial x_i^2}-\frac{1}{c^2}\frac{{\partial }^2}{\partial t^2}\left(i=1,2,3\right)$

With the Z axis chosen in the direction of propagation $h_{ik}$ , ${\square }^2$ becomes

${\square }^2=\frac{{\partial }^2}{\partial z^2}-\frac{1}{c^2}\frac{{\partial }^2}{\partial t^2}$

and then the partial derivatives of $h_{ik}$ with respect to x and y vanish. Thus a plane gravitational wave is determined by two quantities $h_{12}$ and $h_{11}=-h_{22}$ in other words, gravitational waves are transverse waves whose polarization is determined by a symmetric tensor of the second rank in the XY plane, the sum of whose diagonal terms $h_{11}+h_{22}$ is zero. So the plane wave can be written a the sum of two components

$h_{11}=h_{+}\text{cos}\left(n_{w}t+{\alpha }_1\right)$ (1)

$h_{12}=h_{\times }\text{cos}\left(n_{w}t+{\alpha }_2\right)$ (2)

where $n_w$ is the frequency of the wave, ${\alpha }_1$ and ${\alpha }_2$ are the phase difference, $h_{+}$ and $h_{\times }$ are the amplitude of the wave in the two orthogonal directions in the transverse plane. From the equation of geodesic deviation we obtain the acceleration components of GW as in [10] .

$\frac{{\text{d}}^2{n}^i}{\text{d}{t}^2}=\frac{1}{2}\frac{{\partial }^2{h}_{ik}}{\partial t^2}{n}^k$ (3)

From which we can obtained the potential like function U of GW in terms of the orbital elements

$U=\frac{1}{2}{h}_1\left(y^2-x^2\right)-h_{2}xy$ (4)

where $h_1=\frac{1}{2}\frac{{\partial }^2{h}_{11}}{\partial t^2}$ and $h_2=\frac{1}{2}\frac{{\partial }^2{h}_{12}}{\partial t^2}$ , and we have

$x=r\left(\cos \Omega \cos \left(f+\omega \right)-\cos i\sin \Omega \sin \left(f+\omega \right)\right)$ (5)

$y=r\left(\sin \Omega \cos \left(f+\omega \right)+\cos i\cos \Omega \sin \left(f+\omega \right)\right)$ (6)

where $i,\Omega ,\omega$ and $f$ are the inclination, longitude of node, argument of perigee and the true anomaly of an orbit respectively. Substituting Equations ((1), (2), (5) and (6)) into Equation (4) we obtain

$\begin{array}{c}U=ϵr^2\{\frac{s^2}{8}[\cos \left(2\Omega +n_{w}t+{\alpha }_1\right)-\cos \left(2\Omega -n_{w}t-{\alpha }_1\right)\\ -A_g\left[\sin \left(2\Omega +n_{w}t+{\alpha }_2\right)+\sin \left(2\Omega -n_{w}t-{\alpha }_2\right)\right]]\\ +\left(\frac{c^2}{16}+\frac{c}{8}-\frac{1}{16}\right)[\cos \left(2f+2\omega +2\Omega +n_{w}t+{\alpha }_1\right)+\cos \left(2f+2\omega +2\Omega -n_{w}t-{\alpha }_1\right)\\ -A_g[\sin \left(2f+2\omega +2\Omega +n_{w}t+{\alpha }_2\right)+\sin \left(2f+2\omega +2\Omega -n_{w}t-{\alpha }_2\right)\\ +\sin \left(2f+2\omega -2\Omega +n_{w}t+{\alpha }_2\right)+\sin \left(2f+2\omega -2\Omega -n_{w}t-{\alpha }_2\right)]]\\ +\left(\frac{c^2}{16}-\frac{c}{8}+\frac{1}{16}\right)\left[\cos \left(2f+2\omega -2\Omega +n_{w}t+{\alpha }_1\right)+\cos \left(2f+2\omega -2\Omega -n_{w}t-{\alpha }_1\right)\right]\}\end{array}$ (7)

where

$c=\cos i,s=\sin i,ϵ=\frac{1}{2}{n}_w^2{h}_{\times }\text{and}{A}_g=-\frac{h_{\times }}{h_{+}}$ (8)

The Hamiltonian will be in the form

$H=H_0+ϵH_1$

$H_0$ is the unperturbed Hamiltonian given by

$H_0=-\frac{{\mu }^2}{2L^2}$ (9)

$\mu =GM$ the gravitational constant and $L=\sqrt{\mu a}$ . The perturbed Hamiltonian $H_1$ due to the gravitational wave is given by

$\begin{array}{c}{H}_1=\frac{L^4}{{\mu }^2}\frac{r^2}{a^2}\{\frac{s^2}{8}[\cos \left(2\Omega +n_{w}t+{\alpha }_1\right)-\cos \left(2\Omega -n_{w}t-{\alpha }_1\right)\\ -A_g\left[\sin \left(2\Omega +n_{w}t+{\alpha }_2\right)+\sin \left(2\Omega -n_{w}t-{\alpha }_2\right)\right]]\\ +\left(\frac{c^2}{16}+\frac{c}{8}-\frac{1}{16}\right)[\cos \left(2f+2\omega +2\Omega +n_{w}t+{\alpha }_1\right)+\cos \left(2f+2\omega +2\Omega -n_{w}t-{\alpha }_1\right)\end{array}$

$\begin{array}{c}-A_g[\sin \left(2f+2\omega +2\Omega +n_{w}t+{\alpha }_2\right)+\sin \left(2f+2\omega +2\Omega -n_{w}t-{\alpha }_2\right)\\ +\sin \left(2f+2\omega -2\Omega +n_{w}t+{\alpha }_2\right)+\sin \left(2f+2\omega -2\Omega -n_{w}t-{\alpha }_2\right)]]\\ +\left(\frac{c^2}{16}-\frac{c}{8}+\frac{1}{16}\right)\left[\cos \left(2f+2\omega -2\Omega +n_{w}t+{\alpha }_1\right)+\cos \left(2f+2\omega -2\Omega -n_{w}t-{\alpha }_1\right)\right]\}\end{array}$ (10)

3. Equations of Motion

The Delaunay canonical variables $\left(l,g,h,L,G,H\right)$ are defined in terms of the Keplerian elements $\left(a,e,i,\Omega ,\omega ,M\right)$ by

$\begin{array}{l}l=M=n\left(t–\tau \right);L=\sqrt{\mu a}\\ g=\omega ;G=\sqrt{\mu a}\left(1-e^2\right)\\ h=\Omega ;H=G\cos i\end{array}$ (11)

The considered systems of differential equations in terms of canonical elements under the effect of GW are

$\frac{\text{d}{L}_i}{\text{d}t}=-\frac{\partial H}{\partial l_i};\frac{\text{d}{l}_i}{\text{d}t}=\frac{\partial H}{\partial L_i}\left(i=1,2,3,4\right)$ (12)

Then

$l^{·}=\frac{\partial H}{\partial L};g^{·}=\frac{\partial H_1}{\partial G};h^{·}=\frac{\partial H_1}{\partial H}$ (13)

$L^{·}=-\frac{\partial H_1}{\partial M};G^{·}=-\frac{\partial H_1}{\partial g};H^{·}=-\frac{\partial H_1}{\partial h}$ (14)

Using Equations ((11), thus (13) and (14)) yield

$a^{·}=2\sqrt{\frac{a}{\mu }}{L}^{·}$ (15)

$e^{·}=\frac{\left(1-e^2\right)}{n{a}^{2}e}{L}^{·}-\frac{\sqrt{\left(1-e^2\right)}}{n{a}^{2}e}{G}^{·}$ (16)

$i^{·}=\frac{\cot i}{n{a}^2\sqrt{\left(1-e^2\right)}}{G}^{·}-\frac{\csc i}{n{a}^2\sqrt{\left(1-e^2\right)}}{H}^{·}$ (17)

${\Omega }^{·}=\frac{\partial H_1}{\partial H}=\frac{\partial H_1}{\partial i}\frac{\partial i}{\partial H}$ (18)

${\omega }^{·}=\frac{\partial H_1}{\partial G}=\frac{\partial H_1}{\partial e}\frac{\partial e}{\partial G}+\frac{\partial H_1}{\partial i}\frac{\partial i}{\partial G}$ (19)

$M^{·}=\frac{\partial H_0}{\partial L}+\frac{\partial H_1}{\partial L}=\frac{\partial H_0}{\partial L}+\frac{\partial H_1}{\partial a}\frac{\partial a}{\partial L}+\frac{\partial H_1}{\partial e}\frac{\partial e}{\partial L}$ (20)

We solve this equations numerically using Runge-Kutta fourth order method, the mathematical program written by language of MATHEMATIC V10.

4. Solution and Results

We expand the trigonometric functions of the true anomaly $\frac{r^2}{a^2}$ , $\frac{r^2}{a^2}\cos 2f$ and $\frac{r^2}{a^2}\sin 2f$ of Hamiltonian (10) in series of the trigonometric functions of

mean anomaly $M$ up to order three in the eccentricity [11] . Changing the independent variable from the time $t$ to the mean anomaly $M$ , substituting the Hamiltonian into the equations of motion and solving the equations numerically to obtain the first order perturbation of GW on Jupiter. Table 1 represents the orbital elements of Jupiter in the epoch 2000, considering Bursts of GW with

frequency $56\times {10}^{-6}$ Hz and amplitude $15\times {10}^{-22}$ . ${\alpha }_1=0$ and ${\alpha }_2=\frac{\text{π}}{2}$ . The

variations of Jupiter’s elements due to the considering effect of gravitational waves are presented in Table 2, during the interval between zero and 4π as in Figure 1.

5. Conclusion

The resulting solution for our model from Equations (15) to (20) includes short-period as well as long-period terms for practical application it appears more expedient to calculate the perturbations on a revolution by revolution. The averaging process over the fast variable or short-period comes after the differentiated with respect to slow variable or long-period. The short-period perturbation on semi-major a, the eccentricity e, and the mean anomaly $M$ increase with each revolution during the time interval $\left(t-t_0\right)$ from zero to 4π. The long-period perturbation on the inclination i, Ω, and ω takes long time to change with revolution in the same interval as seen in Figure 1. This paper provides a good approach for long-period orbital evolution studies for planets orbiting under the effect of the gravitational waves.