The Maxwell equations trigger the question on the potential existence of magnetic charges and currents to obtain the implied better symmetry of the partial differential equations representing the electric and magnetic field. Works in that direction were presented by, for example, Heaviside and Dirac. Ref.  calculated the transition radiation from the hypothetical Dirac magnetic monopoles and showed that any radiation from monopoles should be several orders of magnitude larger than from particles with unit electric charge. Ref.  generalized Maxwell equations (consistent with special relativity) by assuming the existence of “positive electricity” for which another set of partial differential equations were set up to be solved simultaneously with the standard set of equations (“negative electricity”) of Maxwell. That theory was used to offer a possible explanation of the earth’s magnetic and gravitational fields and the maintenance of the earth’s charge. Ref.  derived the potentials and fields of dipoles consisting of standard electric charge and currents and of hypothetical magnetic monopoles and their currents.
In Appendix of this paper it is proposed that the magnetic charges and their currents are of pure electrical origin and the potentials are accordingly given. The new gravitational vector potential was used (in addition to the standard one) as the base for determining the total power radiated during in-spiral of a binary system. Tentative calculations show that the impact on the orbital decay is a small increase (~ 0.5% relative deviation in the binaries’ separation at 0.5 s into the transient) using initial conditions and masses approximately matching the GW150914 event. During that study it was noticed that by changing the standard gravitational vector potential by a factor of about 2, the orbital decay significantly increased. The work presented here was triggered by that finding.
Researchers as Maxwell and Heaviside applied the electromagnetic equations to gravitation in the 19 century. This is currently known as gravitoelectromagnetism (GEM). Ref.  derived, from the GEM vector potential, the equation for the total power radiated by binary systems moving at low speeds in circular orbits. That equation is in striking similarity with the linear theory of relativity (LGR) result, the only difference being the value of a proportional constant N (2/5 vs. 32/5 in LGR). Ref.  extended the work of  to relativistic speeds (using numerical methods) and applied GEM and a relativistic 2nd law of Newton (ENET), to derive a relativistic Kepler’s 3rd law (K3L), the results were equivalent to N ~ 8/5.
The finding about is interpreted as the possibility that the gravitational wave speed in vacuum (cg) be smaller than the speed of light in vacuum (c) and therefore the gravitational permeability of free space, in analogy with electromagnetism, would be . The mentioned factor of 2 then would imply that . Note that the gravitational permittivity is kept as .
The detection of the almost simultaneous GW170817/GRB170817A  events about 3 years ago put extremely high constraint on the deviation of cg from c but, as noted in  the constraint is valid under some assumptions with respect to the emission time of the sources and to the intergalactic medium and the reported bounds can be circumvented if the delay between the events can be greater than 10 sec. or if the electromagnetic wave (EW) can precede the GW (see more discussion in Section 4).
Even if is confirmed without doubts, the present work has teaching importance as a first introduction to GW without using the complicated physics and mathematics inherent in general relativity (GR) while at the same time recognizing the historical and remarkable scientific relevance of Einstein GR. Additionally the computationally efficiency of the numerical solution of many-body-problems can potentially be improved (with respect to numerical GR) by solving equations similar to the force balance equation presented in Subsection 2.1.
The main objective of this work is to assess the impact of the gravitational wave and force propagation speed on the orbital frequency of circular orbits under a GEM framework.
The rest of this work is structured as follows. In Section 2.1, K3L for circular orbits is extended considering GEM and ENET. In Section 2.2 the equations for the calculation of the power radiated during in-spirals of binaries are derived. In Section 2.3 the gravitational orbit decay calculation is described. In Section 3 a purely geometric method to calculate cg is proposed. Section 4 presents a comparison of the results obtained here with the ones from LGR for 3 simulated GW events. In this section cg is also calculated along with its deterministic constraints for 2 GW events. Section 5 presents a summary and some concluding remarks.
2. Gravitational Decay of Circular Orbits
2.1. Extension of Kepler’s 3rd Law
The balance between the Lorentz force extended to gravitation and the inertial force can be written as 
: Gravito-electromagnetic (GEM) force acting on a point mass due to another point mass ;
: Gravito-electric (GE) field acting on body a;
: Velocity of body a in the center of mass reference frame;
: Gravito-magnetic (GM) field acting on body a;
: Newtonian acceleration, in polar coordinates:
: Radial and azimuthal coordinates respectively.
: Radial (pointing from the origin of the laboratory reference frame to a) and azimuthal (counter clock wise from x to y) unit vectors respectively which can be written as
The number of dots on top of a variable represents the order of the time derivative.
Note that the Newtonian acceleration is corrected by where is the ubiquitous Lorentz factor. This correction yields the correct values of the intrinsic (two-body problem) perihelion precession of the planets of the solar system . Note also that the radiation reaction is ignored.
The Lienard-Wiechert potentials extended to gravitation (ELW) is written as
: Coefficient to consider the impact of the speed (delay) of the GEM force on the orbit decay.
: Coefficient to consider the impact of the GM permeability on the orbit decay .
: Gravitational wave speed.
: Distance from the position of the point mass to the observation point (separation).
To obtain ELW potentials the following conversions are made: and .
Following similar approach as the ones described in   the following is obtained for the fields of a point mass moving in an arbitrary motion:
where, : Separation vector, velocity and acceleration of the source of the fields at the retarded time.
The acceleration is to be calculated as : Newtonian acceleration.
Assuming circular orbits with piece-wise constant angular speed (ω) the fields acting on mass a are written as:
: Radial unit vector (pointing from the origin of the laboratory reference frame to b) and azimuthal (counter clock wise, from x to y) unit vector respectively which are determined as
Equation (1) for circular orbit with a constant angular speed can be written as:
Substituting Equation (6) and (7) into Equation (8), considering small ( ), a CM reference frame ( ) and ignoring in the left hand side the dependency on , leads to the following extended Kepler’s 3rd law for circular orbits:
where, , , .
The CM relations and were assumed to hold when deriving Equation (9).
Note that making in Equation (9), Kepler’s 3rd law for circular orbit is recovered.
2.2. Power Radiated during In-Spirals of Binaries
For a binary system the superposition of the ELW potentials at the observation point R are written as
where, are the distances from mass a and b respectively, to the observation point. The positions and velocities are evaluated at their respective retarded time: and .
The positions and velocities (in the CM) of the components of the binary system moving with piece-wise constant angular speed are:
Expanding in Taylor series the velocities around tR, using the binomial expansion in the law of cosines relating with , considering and assuming that the following is obtained:
Rearranging Equation (12) as , and expressing the Cartesian unit vectors in spherical coordinates with the choice of as in Ref. :
The assumption was made to obtain Equations (15), (16). The inertial accelerations were corrected with the relativistic factors , .
The GE radiation field is determined as which depends only on the components of the vector potential that are transverse to the observation direction:
Note that , were considered constants for the calculation of the derivatives in Equation (17).
The power radiated per unit area in the direction of a monochromatic wave propagation is determined as where the GM radiation field is included.
Note that the average value, which is optional in EM and required in GR  is not used in .
The total power radiated is calculated as
, , ,
Note that, does not depend on R.
2.3. Gravitational Orbit Decay
The total classical mechanical energy of the binary system can be written as , therefore
From Equation (9),
Taking the time derivative in both sides:
, , ,
Substituting into Equation (19):
Plug it into Equation (18):
3. Determination of the Speed of Gravitational Waves from Arrival-Time Delay Measurements
As pointed out before, under some assumptions relative to the emission time of the sources and to the intergalactic medium the coincident events GW170817/GRB170817A provided a very tight constrain on , the bounds of which can be circumvented if the GW-EW delay is allowed to be greater than 10 s or if the EW can precede the GW . Despite these extreme narrow bounds on that reference proposed an alternative (not limited by the same systematics and without expecting to compete with the constraint on those bounds) to measure cg based on the detection of continuous monochromatic signals emitted from, for example, rapidly rotating neutron stars.
In this section a method to measure cg is proposed that is also not limited by the same systematics without having the expectation to compete with simultaneous GW-EW event-based measurements. Assuming that accurate and model-independent GW arrival time delays among the detectors can be obtained, cg could be determined as follows:
The arrival time delays between the 3 detectors are related to the GW speed through the following equations:
: Arrival time delay between indicated detector sites.
: Distance travelled by the wave front (arriving first at the first indicated site).
: Angle between the wave propagation direction and the line joining Virgo and Livingston sites.
: Angle between the wave propagation direction and the line joining Livingston and Hanford sites.
D: Straight distance (through earth) between indicated detector sites.
Note that the wave direction is perpendicular to a vector originating at the second indicated site and ending at the wave vector. Note also that the projection of the wave vector on the detectors’ plane does not change those angles and that the angles are the same along their corresponding joining lines as long as the wave source is far enough from the detectors.
The angles are related through
: Angle between the line joining Virgo and Livingston and the line joining Livingston and Hanford.
From the law of cosines (generalized Pythagoras theorem)
Solving for in Equation (22) and substituting it (along with Equation (23)) into Equation (21) the following is obtained:
is then obtained from Equation (22) or (21).
Note that if the angles could be measured (or accurately inferred from measurements), detectors’ signals, not complying with the equality of in Equation (21) and in Equation (22): , could be rejected.
4. Computational Results and Analysis
The values of and will be determined in such a way that the period decay, , of the Hulse-Taylor (H&T) binary pulsar (PSR 1913+16) assuming a circular orbit, is reproduced (in another work the elliptical orbit case will be addressed). The value to be reproduced is the one from GTR  :
assuming . The GTR result (for elliptical orbit) is in very good agreement with experiment. The hypothetical H&T result for circular orbit is used as a target because the H&T value for elliptical orbit is very robust: direct measurements of orbital parameters are made and it has been systematically reproduced using experimental results from more than 40 years (measurements can still be made).
Relatively recent detected gravitational waves (e.g. GW150914 and GW170817) are not used as a target considering that they are not as robust as the H&T experimental results: the GW measured data are not direct measurements of orbital parameters, and they are not reproducible (the experiment cannot be repeated on demand).
Additionally in Ref.  information can be found about doubts that a group of physicist of the Niels Bohr institute in Copenhagen had, on the gravitational wave signals extracted from the LIGO detectors. Among other things it was pointed out that the extraction of the GW signals should use a model independent algorithm. Ref.  developed such an algorithm and was able to reproduce the results for the GW150914 but they found no statistical significant indication of any common signal for the GW151226 event. Also Ref.  reported that the detector noise correlations between the LIGO sites, at the time of the event, were maximized for the same time lag (delay) as the one found in GW150914.
The GW170817/GRB170817A coincidental events (which apparently implies that the speed of the gravitational wave is the same as the speed of light: ), have non trivial measurement/interpretation uncertainties due to the uniqueness of the events. For example Ref.  raised doubts on the correlation between the GRB and the GW. They proposed an alternative mechanism to explain the very low statistical probability (~10−10) of a coincidental origin of the events for the model of LIGO-Virgo which reports a probability of 5 × 10−8 that both events happened by chance.
Table 1, shows 2 sets of special values of and that approximately reproduce the GR value along with the set representing the standard GEM which yields significant smaller value of . The 2nd set ( , ) implies while the gravitational force news is transmitted at the speed c. The last set ( , ) represents a case for which both the GEM force news and the gravitational wave propagate at the same speed as it happens in classical EM. The shown in Table 1 are average values (the time profiles are very oscillatory). The simulation input are : , , , . An integration step of was used. The simulation time was 365 days.
The distance to the GRB170817A event that results in a time coincidence with the GW170817 event at the GW detectors for , was calculated (assuming no common distance to the source and assuming equal time of emission) For arrival time coincidence of the signals (moving in vacuum):
So for and . Therefore, in this example, an isotropic (or a disk/beam reaching the earth) GRB source located at about 12.6 Mpc behind the GW source will be detected in time-coincidence with a GW moving at about 0.76c. Note that time coincidences could also happen for astronomical events resulting in emission of EM waves at any distance if no
Table 1. Period decay vs. the speed (as a fraction of c) of the GW and the force propagation.
emission time coincidence is assumed. Note also that if a large difference of really exists no real coincidence will be detected on earth for simultaneous time emissions and common spatial origin of the events since as pointed out in  that would imply a large time offset.
Figures 1-4 show in-spiral simulations of three GW events. The integration step for all transients was .
Figure 1 shows the orbital frequency (Hz) calculations using initial conditions and masses corresponding to the GW150914 event and taking ( ). The 4 cases shown are:
K3L (Kepler’s 3rd law): Make in Equation (9) and in q (Equation (20)).
EM (Electro-Magnetism): Make and .
EM-ENET: Use Equation (9) and Equation (20) without approximation.
LGR (Linear General Relativity): . . Use K3L to determine , where, , , are the initial transient time, initial separation and the Schwarzschild radius respectively.
Figure 2 shows also the results for the GW150914 but with ( ). It is remarkable that by just reducing by 24%, the difference (w.r.t. LGR) of the average slopes of the curves is greatly reduced (by a factor greater than 3 for the K3L case). Note that the slope of EM-ENET is smaller than that of K3L (the opposite happens in Figure 1) but note also that the speed in the Newtonian acceleration correction was not changed (it remained as the speed of light in vacuum).
Figure 1. Orbital frequency (Hz) vs. Time (s) LGR vs. GEM (cg = c) for GW150914.
Figure 2. Orbital frequency (Hz) vs. Time (s) LGR vs. GEM (cg = 0.76c) for GW150914.
Figure 3. Orbital frequency (Hz) vs. Time (s) LGR vs. GEM (cg = 0.76c) for GW170814.
Figure 4. Orbital frequency (Hz) vs. Time (s) LGR vs. GEM (cg = 0.76c) for GW170817.
Figure 3 shows the orbital frequency for using initial conditions and masses corresponding to GW170814 : , , and . Figure 3 shows similar behavior as Figure 2 except for the small ripples which are visible on the plot depending on the axis scale, plotting frequency, etc. Note that LGR does not have those ripples since S needs to be averaged over several wave periods in GR.
Figure 4 shows the orbital frequency for the GW170817 event using The input used for this transient was: , , and . Note the close agreement for this relatively slow transient in all cases until about 5 sec. and until about 11 sec. for K3L and EM-ENET. This is consistent with the H&T simulation (extremely slow transient) where practically the same results were obtained for all cases.
For and ( ) it is expected to yield similar results since calculations for the GW150914 event indicated that. The following values of N were calculated (using Equation (4)) that approximately reproduced the results for the GW150914 event: for K3L, EM and EM-ENET respectively taking . Note that for LGR , K3L is used and .
Table 2 Shows the locations of the 3 GW detectors  in the Geodetic reference frame (GRF) along with the calculated distances (through earth) between the sites. To calculate the distances, the GRF data was converted to Cartesian reference frame according to :
where f: Latitude, λ: Longitude, h: Ellipsoidal height, a: Ellipsoid semi-major axis: 6378137.0 m, b: Ellipsoid semi-minor axis: 6356752.0 m. The axis values correspond to the World Geodetic System 1984. In this application h does not make significant impact on the results. The distance between the site i and j was calculated as
Equation (24) yields .
Table 3 shows the arrival time delay of the GW1708 (14 and 17) signals    along with calculated arrival angles and the GW speed (Equations (22)-(25)). Note that the obtained angles are consistent with the arrival time delays (the smaller the angles the larger the delays). Note that cg is about 0.99c and 1.08c for GW1708 (14 and 17) respectively.
Table 2. Calculated distances between GW detectors’ sites. The detectors’ coordinates were taken from .
Table 3. Gravitational speed calculation from the arrival time delay between detectors’ sites.
The same results are obtained if the Cartesian coordinates of Ref.  are used instead of converting from geodetic system to Cartesian coordinates (except for GW170814 where cg deviates in −0.1 km/s). If the approximation (light travel time) is made and as noted in , and 316418.3 km/s is obtained for GW1708 (14 and 17) respectively.
The uncertainty of the GW speed can be determined making further calculations based on the uncertainties of the arrival time delays. Uncertainty of ±1 msec. in the arrival time delays were estimated in  for GW170814 event. Calculations (using Equation (22)-(25)) considering those uncertainties result in a speed constraint of 0.89c < cg < 1.11c. For GW170817 assuming 10% uncertainty  in the arrival time delays, the GW speed constraint was 0.98c < cg < 1.20c. These results put significantly more constraint than the bounds of 0.55c < cg < 1.42c reported in . Perhaps combining both methods could further constrain cg.
Note that the calculated values of cg are model dependent because the determination of the time delays assumes in some part of the experimental procedure to detect and localize the GW  . Note that the signal templates used to detect/extract GW signals in LIGO have to be modified for the detection/extraction of CGW to account for .
It is noted that Ref.  reported a constraint of the speed of gravity of 0.85c < cg < 1.22c based on the Jovian deflection experiment (measurements of the propagation of quasar’s radio signal past Jupiter), in that reference a distinction is made on the meaning of c in different parts of the GR equations.
Note that if cg is given, the expected could be inferred for the GW150914 event (VIRGO was not operational during that time) using Equations (22)-(25). Perhaps the angular sky location could be improved as well. For example assuming and using  the following is obtained: , , .
5. Summary and Concluding Remarks
Speeds of GW and force propagation were obtained that account for the expected period decay of a hypothetical Hulse-Taylor binary pulsar moving in circular orbits. The so-obtained speeds were used to simulate the in-spiral process of 3 GW events. The results are in significantly better agreement (with the results of LGR) than the results obtained using the speed of light for both the GW and the gravitational force. A method to calculate the speed of GW was proposed and applied to 2 GW events. Uncertainties of the arrival time delays were used to deterministically constrain the GW speed. It is remarked the need for measuring and/or calculating model-independent delays to determine and constrain the GW speed.
It is noted that the accuracy level of the results obtained here (concerning the 3 GW events) is not determined necessarily by the agreement with the LGR results but by its agreement with, for example, model-independent measured profiles of the detectors’ strain (h) and/or with orbital parameters as, for example, the orbital frequency, calculated using experimental raw data from the LIGO-VIRGO detectors.
It is envisioned to extent the GEM model developed here to consider elliptical orbits.
It could be worthy to assess the impact of the implied kinetic energy of ENET on the orbital decay and also, if needed, the impact of the magnetic charges and currents using the potentials obtained in Appendix.
Ref.  modified the Maxwell equations (to include magnetic charge and current) as follows:
, , ,
The conservation equations were given by , . The electric and magnetic field were written as , .
Note that the variables with subscript “m” represent the new magnetic sources.
Assuming Lienard-Wiechert potentials for point magnetic charges also:
, , ,
Assuming that all magnetic charges are of electrical origins, the new potentials are expressed as , . Note the quadratic dependence of on the speed.
The vector potentials for gravitation are then written as , .
is calculated in the main text ( ). is similarly determined as , where,
The total radiation field (Coulomb gauge) is then written as
The Poynting vector is calculated as .
Similarly for .
The total power radiated is then determined as .
 Swann, W.F.G. (1927) A Generalization of Electrodynamics, Consistent with Restricted Relativity and Affording a Possible Explanation of the Earth’s Magnetic and Gravitational Fields and the Maintenance of the Earth’s Charge. Philosophical Magazine Series 7, 3, 1088-1136.
 Banales, M.I. (2019) Fundamental Physics in the Era of Gravitational Wave Astronomy: The Direct Measurement of Gravitational Wave Polarizations and Other Topics. PhD Thesis, California Institute of Technology, Pasadena.
 Taylor, J.H. and Weisberg, J.M. (1982) A New Test of General Relativity: Gravitational Radiation and the Binary Pulsar PSR 1913 + 16. The Astronomical Journal, 253, 908-920.