This notebook has been created with Wolfram’s Mathematica 12.3.1.0

The cells with Light Gray background color contain information about the calculus that is being performed.
The executable cells (that must be executed in order) have been assigned a Light Blue background color.

Notations:
As usual in Wolfram’s Mathematica, when two or more characters are written together, without spaces between them, they correspond to names of variables. (Sub-indexes will not be used in this notebook)
Instead of the multiplication sign “*” a [space] between the names of variables will often be written since Mathematica interprets such spaces as multiplication.
Examples:
abc    is a possible name of a variable.
ab c    means the product of the variable ab by the variable c.

Gravitation  with out a Tempo Rate Law (TRL).

Abstract:
Applying the paradigms of the EVE model of the aether, see
https://www.eterinica.net
https://www.eterinica.net/EVEG/Force_between_neutral_bodies.pdf
https://www.eterinica.net/EVE9/Eve9.pdf
this notebook calculates:
1) the gravitation force FG between two hydrogen atoms disregarding the influence of any Tempo Rate Law.

1) Gravitation force FG between two hydrogen atoms disregarding the influence of any Tempo Rate Law.

According to the model, the gravitation force that an electrically neutral body #1 exerts on an electrically neutral body #2 is just the net residual (non fully cancelled) electrodynamic force that the electrons and the protons of the body number #1 exert on the electrons and the protons of the body number #2
Supposing for example that the interacting bodies are simply two hydrogen atoms, then the net force F12 suffered by the atom #2 due to the presence of the atom #1 can be considered  the sum of the 4 following forces:
F12 = FPp + FPe + FEp + FEe 
where
FPp   is the force exerted by the proton of the atom #1 on the proton of the atom #2.
FPe   is the average force exerted by the proton of the atom #1 on the electron of the atom #2. (The “average” is over all directions of the velocity w of the electron #2).
FEp    is the average force exerted by the electron of the atom #1 on the proton of the atom #2. (The “average” is over all directions of the velocity w of the electron #1).
FEe   is the average force exerted by the electron of the atom #1 on the electron of the atom #2. (The “average” is over all directions of the velocities w of both electrons #1 and #2).

The force f12x that a moving electron #1 exerts on a distant moving electron #2 will be calculated first.
Consider the following Fig[1]

Fig[1]

$HistoryLength=0;

SetOptions[EvaluationNotebook[],DynamicUpdating->False]

Force f12x  that a “moving” atomic  electron #1 exerts on a “moving” atomic electron #2 
See Fig[1]
The electron #1 always remains approximately at the same location #1 although its velocity w1 can be found with equal probability pointing in any direction of space. Its speed w1 always remains the same.
The electron #2 always remains approximately at the same location #2 although its velocity w2 can be found with equal probability pointing in any direction of space. Its speed w2 always remains the same.

Those two electron “systems” are considered adequate to represent  two electrons of randomly oriented orbits of two distant atoms.

d is the distance at all epochs between the location #1 of the electron #1 and the location #2 of the electron #2

Calling X the direction joining the locations #1 and #2 then, due to symmetry reasons, only the X component of the force f12 that the electron system #1 exerts on the electron system #2 does not cancel. The other two Cartesian components, X and Y, cancel. Therefore only the X component f12x of the force will be calculated.

The force f12x (that will deduced in the calculus below) corresponds to the FEe term of the net gravitation force but the expression of such force f12x  can be directly applied to calculate the other 3 forces as follows:
FPp  =  f12x   taking  w1=0 and w2=0
FPe  = -f12x  taking  w1=0 and w2=w
FEp  = -f12x  taking  w1=w  and w2=0

Clear[redis,σS,σI,aS,bS,aI,bI,VM,N0,d,v,α1,α2,φ2,vR,vRx,vRy,vE,vEx,vEy,t,tE,w1,w2,w,c,h1,igrdfx,Fα1, Fα2,Fα1α2,τ,μ]

The switch cross section of a unit charge particle (in this case the electron #1)  that is responsible for the aetherino redistribution created by the particle (in this case the electron #1), is according to the model:

σS[vE_]:=aS Exp[-bS (vE^2)]

where
vE is the speed of the colliding aetherino relative to the redistributing particle (here the electron #1)
aS  is a constant with the dimension of surface
bS  is a constant with the dimension of (speed^(-2))

The canonical distribution of the aether, that is assumed to be the distribution of aetherino speeds of the aether local to the redistributing particles, is according to the model:

ρ[v_]:=((4 N0)/((Sqrt[π])(VM^3))) (v^2) Exp[-((((v/VM)))^2)]

where
v    is the speed of an aetherino of the local aether relative to the aether as a whole.
 (Note: It can be seen that as far as the speed of the electron relative to the aether as a whole is not much bigger than a few times  c,  then the redistribution created by the electron does not vary noticeably if it is supposed that such source electron is at rest relative to the aether).
 N0    is the average number of aetherinos in unit volume of the canonical aether. 
 VM   is the speed for which the canonical aether has a maximum number of aetherinos. It is supposed that VM >> c
 
 Note: writing the canonical distribution ρ[v] as above, it can be seen that N0 is always the total number of aetherinos in unit volume whatever the value of VM.

The redistribution of aetherinos created by a unit charge particle (like the electron #1) was shown (in the paper  https://www.eterinica.net/EVEG/Force_between_neutral_bodies.pdf ) to be:

redis[v_]:=σS[v]((ρ[v])/(8π)) v

Remember that in the model the redistribution of aetherinos created by a particle is the excess or deficit of aetherinos of a given speed in relation to those that would emerge the volume assignable to the particle if the particle was not there. 
Remember also that the model considers that there are two types (p-type and n-type) of aetherinos and assumes that a positive charge switches the n-type aetherinos that collide with it into p-type aetherinos while a negative charge switches the p-type aetherinos that collide with it into n-type aetherinos.

redis[vE]

An aetherino that, in the reference frame of description (in which the interacting atoms are at rest) travels at a speed v from the atom #1 to the atom #2 has a speed vE relative to the electron #1 of velocity w1 given by: 
(see Fig[1])

vE=((((((v-w1 Cos[α1]))^2)+(((w1 Sin[α1]))^2)))^(1/2))

Note: such speed vE (that depends on the pertinent velocity w1 of the electron) has an important influence on the redistribution emerging from the electron as can be acknowledged considering the switch cross section σS[vE] assumed above.
Note: although, whatever the velocity w1 of the electron, it can be supposed that it is approximately bathed by a canonical aether so that its redistribution of aetherino speeds relative to the electron itself  does not depend on w1, the redistribution seen by the target electron #2 does depend significantly on the velocity w1 of the electron #1 causing such redistribution (as explained in the paper of the model that analyzes de force between neutral bodies)

redis[vE]

The impulse cross section  σI  of a unit charge particle, is according to the model: 

(Note: the n-type aetherinos give impulse (instead of switching their type) to the negative charges with which they collide. Similarly, the p-type  aetherinos give impulse (instead of switching their type) to the positive charges with which they collide).

σI[vR_]:=aI Exp[-bI (vR^2)]

where vR is the speed of the colliding aetherino relative to the target particle (in this case the electron #2).
Considering, see Fig[1], that the target electron #2 has the Cartesian coordinates:
{w2x=w2*Cos[α2], w2y=w2*Sin[α2]Cos[φ2], w2z=w2*Sin[α2]Sin[φ2]} 
and considering a generic aetherino whose velocity in the reference frame of description is {v,0,0} 
then the velocity vR of the aetherino relative to the target electron has the Cartesian components
vR = {v-w2*Cos[α2], 0-w2*Sin[α2]Cos[φ2], 0-w2*Sin[α2]Sin[φ2]} 
And its speed is:
vR = ((((((v-w2 Cos[α2]))^2)+(((0-w2 Sin[α2]Cos[φ2]))^2)+(((0-w2 Sin[α2]Sin[φ2]))^2)))^(1/2))

vR=Simplify[((((((v-w2 Cos[α2]))^2)+(((0-w2 Sin[α2]Cos[φ2]))^2)+(((0-w2 Sin[α2]Sin[φ2]))^2)))^(1/2))]

σI[vR]

In the paper https://www.eterinica.net/EVEG/Force_between_neutral_bodies.pdf  it was deduced that integrating for all aetherino speeds v and averaging for all directions of the velocities w1 of the electron #1 and w2 of the electron #2 the force exerted by the electron #1 on the electron #2 is given by the expression [G-22b]

f12x = (h1/(4 (d^2)))(∫_0^π)Sin[α1](∫_0^π)Sin[α2](∫_0^∞)(((v-w2 Cos[α2])σI[vR]  vR  redis[vE]  (v^2) (v-w1 Cos[α1]))/(vE^4))dvdα2dα1

where h1 is the constant from the relation i1 = h1 vR  that by hypothesis gives the impulse given by an aetherino to a particle when it collides with it at a relative velocity vR.
Calling  igrdfx  the following part of the integrand of the above integral:

igrdfx=(h1/(4(d^2)))(((v-w2 Cos[α2])*σI[vR]*vR redis[vE] (v^2)(v-w1 Cos[α1]) )/(vE^4))

The α2 dependent part of the triple integral 
 f12x = (h1/(4 (d^2)))(∫_0^π)Sin[α1](∫_0^π)Sin[α2](∫_0^∞)(((v-w2 Cos[α2])σI[vR]  vR  redis[vE]  (v^2) (v-w1 Cos[α1]))/(vE^4))dvdα2dα1 =
  =  (∫_0^π)Sin[α1](∫_0^π)Sin[α2](∫_0^∞) igrdfx dvdα2dα1
 is integrated by Mathematica 12.3 giving a closed form solution (in less than 2 minutes in a standard personal computer of the year 2024). Calling Fα2 this integration of the variable α2
 Fα2 = (∫_0^π)Sin[α2]  igrdfx  dα2

Fα2=Integrate[Sin[α2]*igrdfx,{α2,0,π}]

that since the variables of the force that is being evaluated satisfy the conditions of the above solution, its result can be accepted, copied and simplified according to

Fα2=Simplify[(aI aS h1 N0 (v^2) (v-w1 Cos[α1]) (-(1/(16 (bI^(5/2)) (v^2) w2))(E^(-(((1+bS (VM^2)) ((v^2)+(w1^2)))/(VM^2))-bI (((v-w2))^2)+((2 v (1+bS (VM^2)) w1 Cos[α1])/(VM^2)))) (-2 (Sqrt[bI]) (3+4 bI v (v-w2)) (Sqrt[(((v-w2))^2)])+(E^(bI (((v-w2))^2))) (Sqrt[π]) (3+2 bI ((v^2)-(w2^2))) Erf[(Sqrt[bI]) (Sqrt[(((v-w2))^2)])])+(1/(16 (bI^(5/2)) (v^2) w2))(E^(-(((1+bS (VM^2)) ((v^2)+(w1^2)))/(VM^2))-bI (((v+w2))^2)+((2 v (1+bS (VM^2)) w1 Cos[α1])/(VM^2)))) (-2 (Sqrt[bI]) (Sqrt[(((v+w2))^2)]) (3+4 bI v (v+w2))+(E^(bI (((v+w2))^2))) (Sqrt[π]) (3+2 bI ((v^2)-(w2^2))) Erf[(Sqrt[bI]) (Sqrt[(((v+w2))^2)])])))/(8 (d^2) (π^(3/2)) (VM^3) (Sqrt[(v^2)+(w1^2)-2 v w1 Cos[α1]])),{w1 >0,w2>0,v>0}]

Furthermore  Mathematica 12.3 is able to integrate (also in lees than 2 minutes)  the α1 variable of the  force
f12x  = (∫_0^π)Sin[α1](∫_0^π)Sin[α2](∫_0^∞) igrdfx dvdα2dα1
from the knowledge of the integral Fα2 = (∫_0^π)Sin[α2]  igrdfx  dα2   solved above
 I.e. Mathematica 12.3 is able to solve (give an analytical expression) to
   (∫_0^π)Sin[α1]  Fα2 dα1    =    (whose solution will be called    Fα2α1)

Fα2α1=Integrate[Sin[α1]Fα2,{α1,0,π}]

that  again, since the variables w1, w2, v of the force that is being evaluated satisfy the conditions of the above solution, its result can be accepted, copied and simplified according to

Fα2α1=Simplify[(1/(128 (bI^(5/2)) (d^2) (π^(3/2)) (VM^3) w2))aI aS h1 N0 (-(((E^(-2 bI ((v^2)+(w2^2)))) VM (-((2 (E^(-(((1+bS (VM^2)) (((v-w1))^2))/(VM^2)))) VM (Sqrt[(((v-w1))^2)]))/(1+bS (VM^2)))+(((Sqrt[π]) (2 (v^2) (1+bS (VM^2))-2 (w1^2)+(VM^2) (1-2 bS (w1^2))) Erf[(((Sqrt[1+bS (VM^2)]) (Sqrt[(((v-w1))^2)]))/VM)])/(((1+bS (VM^2)))^(3/2)))))/(8 (v^2) w1))+(((E^(-2 bI ((v^2)+(w2^2)))) VM (-((2 (E^(-(((1+bS (VM^2)) (((v+w1))^2))/(VM^2)))) VM (Sqrt[(((v+w1))^2)]))/(1+bS (VM^2)))+(((Sqrt[π]) (2 (v^2) (1+bS (VM^2))-2 (w1^2)+(VM^2) (1-2 bS (w1^2))) Erf[(((Sqrt[1+bS (VM^2)]) (Sqrt[(((v+w1))^2)]))/VM)])/(((1+bS (VM^2)))^(3/2)))))/(8 (v^2) w1))) (2 (Sqrt[bI]) (E^(bI (((v+w2))^2))) (3+4 bI v (v-w2)) Abs[v-w2]+(E^(2 bI ((v^2)+(w2^2)))) (Sqrt[π]) (3+2 bI ((v^2)-(w2^2))) Erf[(Sqrt[bI]) (v+w2)]-(E^(bI (((v-w2))^2))) (2 (Sqrt[bI]) (v+w2) (3+4 bI v (v+w2))+(E^(bI (((v+w2))^2))) (Sqrt[π]) (3+2 bI ((v^2)-(w2^2))) Erf[(Sqrt[bI]) Abs[v-w2]])),{w1 >0,w2>0,v>0}]

But Mathematica is not able to find a closed form solution for the integral of the remaining variable (i.e. is unable to integrate for the aetherino speeds v) of the triple integral that gives the force f12x  exerted by the electron #1 on the electron #2, and therefore the evaluation of such force must rely on a final numerical integration (using Mathematica’s NIntegrate) that needs to previously assign values to all the constants. 
The numerical integration that will be done, introducing in its integrand the result Fα2α1 of the previous integration  of  the angle variables α2 and α1 will then be:
f12x = (∫_0^∞) Fα2α1  dv
but having observed that Mathematica has trouble solving that numerical integration with the limits vMin=0, vMax=∞ 
and having observed that in the present scenario (of aetherino speeds, particle speeds, ...) very good approximations of the results of such integral can be obtained narrowing the integration range to, for example, the  limits vMin=c/3000, vMax= c*100
Hence the numerical integral that will be applied in the following evaluations of the gravitation forces will be
 f12x = (∫_(c / 3000)^(300 c))  Fα2α1  dv
implemented by Mathematica as follows

The product of constants (aI aS h1 N0) that enter the expression of Fα1α2 will be first factored out of the expression since in theory they do not affect the result of the  integral f12x (see below) if they are moved out of the integrand (thus facilitating Mathematica to solve the numerical integral with the required precision). 
Therefore, calling:
fα2α1 = Fα2α1/(aI aS h1 N0)

fα2α1=Simplify[((Fα2α1)/(aI aS h1 N0))];

the expression that will be applied to obtain very approximate numerical values of 
f12x = (∫_0^∞) Fα2α1  dv
will be:
f12x = (aI aS h1 N0) NIntegrate[fα2α1, {v,c/3000,100 c}, WorkingPrecision->50, AccuracyGoal->50]

f12x:=(aI aS h1 N0)NIntegrate[fα2α1,{v,c/3000,100 c},WorkingPrecision->50,AccuracyGoal->50]

Note: Although the control, in Mathematica, of the precision of the numerical results seems confusing and is not mastered by  the author of this notebook, it is considered that many of the numerical evaluations made below with NIntegrate require to be calculated with a  WorkingPrecision and an AccuracyPrecision at least  48 digits.
Fur such high working precision to be applied by Mathematica in its numerical integrations, the constants needed for the calculus (h1, aS, aI, bS, bI, N0, c, VM, w1, w2, d, ...) must be input as Rational (and not Real) numbers, since Rational numbers are considered to have an infinite precision.

Example of the gravitation force between two hydrogen atoms.
(without including the effects of a Tempo Rate Law)
As explained above, the force F12 suffered by the atom #2 due to the presence of the atom #1 can be considered  the sum of the 4 following forces:
F12 = FPp + FPe + FEp + FEe 
where the expression f12x obtained above can be directly applied to calculate the 4 forces as follows:
FPp  =  f12x   taking  w1=0 and w2=0
FPe  = -f12x  taking  w1=0 and w2=w
FEp  = -f12x  taking  w1=w  and w2=0
FEe  =  f12x   taking  w1=w  and w2=w

The following values assigned to the constants of this example are somewhat arbitrary:
Let:

d=1;c=3*(10^8);w=(10^(-9))c;VM=(10^20)c;

aI=1;aS=1;h1=(10^(-40));N0=(10^30);

(* The real number 1.255 when written as 1255/1000 is considered a Rational number of infinite precision *)
bI=(1255/1000)(1/(c^2));
bS=(1255/1000)(1/(c^2));

(* to avoid the appearance of complex infinities due to division by 0, the proton speeds are assigned a value of (10^(-5)) instead of 0 *)

w1=(10^(-5));w2=(10^(-5));FPp=+f12x;

w1=(10^(-5));w2=w;FPe=-f12x;

w1=w;w2=(10^(-5));FEp=-f12x;

w1=w;w2=w;FEe=+f12x;

(* The gravitation force between the two atoms is then *)
F12=FPp+FPe+FEp+FEe

to be continued