May 2026
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.

Rotation curves of galaxies predicted by
	Gravitation with 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.
2) the gravitation force FGtrl between two hydrogen atoms accounting for the influence of a Tempo Rate Law
3) the quotients  FGtrl[d]/FG[d]  for a set of  “astronomical” distances d between the atoms.
4) the orbiting speeds V of the outer stars around the center of a galaxy as a function of their distance d to the center of mass of the galaxy according to the gravitation force FG that disregards the influence of any Tempo Rate Law.
5) the orbiting speeds VTRL of the outer stars around the center of a galaxy as a function of their distance d to the center of mass of the galaxy according to the gravitation force FGtrl that accounts for the existence of a 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 FG suffered by the atom #2 due to the presence of the atom #1 can be considered  the sum of the 4 following forces:
FG = 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,vE,t,tE,w1,w2,w,c,h1,igrdfx,Fα1, Fα2,Fα1α2,Fα2α1,τ,μ]

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)^(100 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 of at least  50 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 FG suffered by the atom #2 due to the presence of the atom #1 can be considered  the sum of the 4 following forces:
FG = 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;N0=(10^30);h1=(10^(-40));

(* 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, in the evaluations of the forces, the average speeds of the protons have here been assigned a value of  (Subscript[w, P]) = (10^(-15)) c instead of 0. But it can be seen that assigning to such proton speeds any other value of the same order or smaller no significant changes are predicted in the intensities of the gravitation forces.

wp=(10^(-15))c;

w1=wp;w2=wp;FPp=+f12x;

w1=wp;w2=w;FPe=-f12x;

w1=w;w2=wp;FEp=-f12x;

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

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

And it can be seen from the above expressions that the gravitation force FG predicted by the model, when the Tempo Rate Law is ignored, decays with the distance as the inverse square law 1/(d^2) .
Therefore when multiplying such gravitation force FG by (d^2), the product  FG*(d^2)  does not depend on the distance and is, of course, equal to the value of FG when the distance between the atoms is assigned, as in the above example, a value d=1

This value  FG = -2.42697391064367385818102`13.052095511366787*^-92 , just obtained, corresponding to the gravitation force between two atoms separated by a distance d=1 , and having assumed the values of the constants
c=3*(10^8);w=(10^(-9))c;  VM=(10^20)c;aI=1;aS=1;N0=(10^30);h1=(10^(-40)); bI=bS=1.255/(c^2)
will be needed below and so is now saved with the name FGd1:

FGd1=-2.426973910644×(10^(-92));

2) Gravitation force FGTRL  between two hydrogen atoms accounting for the influence of a Tempo Rate Law.

As explained in the paper
https://www.eterinica.net/EVE9/Eve9.pdf
the gravitation force predicted by the model when accounting for the influence of a Tempo Rate Law  f[τ]  can be calculated applying the same expressions of  Fα2α1  and of   f12x = (∫_0^∞)  Fα2α1  dv  
deduced above but replacing the constants as follows:
Considering that the gravitation-type disturbance of the aether, reaching today (τ=0) the target atom #2, was “created” by the source atom #1 at the approximate epoch τ=(Subscript[τ, E]) = -d/c then:
*    VM (speed of the aetherinos for which the canonical distribution of the aether reaches its maximum) must be replaced by   VM*f[(Subscript[τ, E])]
*    w1  (average speed of the electrons #1 of the source atoms) must be replaced by  w1*f[(Subscript[τ, E])]
*    bS  (the constant appearing in the  switch cross section of the source particle) must be replaced by bS/((f[τE])^2) or more precisely the value bS=1.255/(c^2) that the model considers reasonable must be replaced by bS = 1.255/(((c f[(Subscript[τ, E])]))^2)
    
And the Tempo Rate Law (TRL) that the model considers plausible since it predicts ad hoc the correct value for Hubble’s constant is:
f[τ] = (1/(1-μ τ))      
and therefore for τ=(Subscript[τ, E]) = -d/c   the TRL takes the form 
f[(Subscript[τ, E])]  = (1/(1 +  μ  d/c))
where  μ  is equal to Hubble’s constant, i.e.    μ = (Subscript[H, 0]) = 2.5 ×(10^(-18))  (sec^(-1))

3) Quotients  FGtrl[d]/FG[d]  for a set of  “astronomical” distances d between the atoms.

With those replacements of the constants VM, w1 and bS imposed by the Tempo Rate Law (TRL), it will now be obtained, for a sample of 100 astronomical distances, a Table of quotients FGtrl[d]/FG[d] where 
  FG[d]   is the gravitation force, when omitting the influence of the TRL, between 2 atoms separated a distance d
  FGtrl[d] is the gravitation force, accounting for the influence of the TRL, between 2 atoms separated a distance d
  
But an equivalent way of finding the quotient FGtrl[d]/FG[d] is to multiply both forces by (d^2) before calculating their quotient, i.e. evaluating first   (d^2) FGtrl[d]   and  (d^2) FG[d] and from them calculating their quotient ((d^2) FGtrl[d]) / (((d^2) FG[d]) that is of course equal to  FGtrl[d] / FG[d].
But it was seen above that the only dependence of FG[d] on the distance is implemented by a factor 1/ (d^2) and therefore the products  (d^2) FG[d] do not depend on the distance and are all equal to FG[1] that is the quantity calculated above and saved as
FGd1=-2.426973910644×(10^(-92))

Therefore the sought Table of quotients FGtrl[d]/FG[d] = ((d^2) FGtrl[d]) / (((d^2) FG[d]) = ((d^2) FGtrl[d]) / FGd1
will have the same elements as the easier table to create finding the quantities  ((d^2) FGtrl[d]) / FGd1  for a sample of  30 astronomical distances.

The numerical integral that needs to be applied to find the elements   ((d )^2) FGtrl[d]   is 
f12xTRL:=(aI aS h1 N0)NIntegrate[(d^2)fα2α1,{v,c/3000,100 c},WorkingPrecision->50,AccuracyGoal->50]
(Note: the product  (d^2)fα2α1 will below be called   fα2α1dd)
(notice that it is just the product of  (d^2) by the numerical integral called above f12x but now with some replaced variables) 
and as before 
(d^2) FGtrl[d]  = (d^2) FPp + (d^2) FPe +(d^2) FEp + (d^2) FEe 
where 
(d^2) FPp = + f12xTRL/.{w1->0,w2->0}  that will be called FPpd
(d^2) FPe = -f12xTRL/.{w1->0,w2->w}   that will be called FPed
(d^2) FEp =-f12xTRL/.{w1->w*f,w2->0}  that will be called FEpd
(d^2) FEp = + f12xTRL/.{w1->w*f,w2->w}  that will be called FEed
 where the f in the  w1->w*f is here the TRL :   f =  f[(Subscript[τ, E])]  = (1/(1 +  μ  d/c))

Clear[σS,σI,aS,bS,aI,bI,VM,N0,d,v,α1,α2,w1,w2,w,c,h1,Fα2α1,fα2α1,fα2α1dd,τ,μ,f,FG,FGtrl,FPp,FPe,FEp,FEe]

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}];

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

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

f:=(1/(1+μ d/c))  (* f is the TRL *)

μ=(5/2)*(10^(-18));

c=3*(10^8);w=(10^(-9))c;aI=1;aS=1;N0=(10^30);h1=(10^(-40));
 (* the same values used above to evaluate FGd1 in the non TRL scenario *)

VM:=(10^20)c*f;

bI=(1255/1000)(1/(c^2));
bS=(1255/1000)(1/(((c*f))^2));

The sought table of  elements ((d^2) FGtrl[d]) / FGd1 will later be deduced from the 4 following tables corresponding to the 4 terms of the sum    (d^2) FGtrl[d]  = (d^2) FPp + (d^2) FPe +(d^2) FEp + (d^2) FEe

wp=(10^(-15))c; (* very small average speeds of the protons *)

w1=wp;w2=wp;
FPpd=Table[{d,f12xTRL},{d,(10^22),(10^25),(10^23)}];

w1=wp;w2=w;
FPed=Table[{d,-f12xTRL},{d,(10^22),(10^25),(10^23)}];

w1:=w*f;w2=wp;
FEpd=Table[{d,-f12xTRL},{d,(10^22),(10^25),(10^23)}];

w1:=w*f;w2=w;
FEed=Table[{d,f12xTRL},{d,(10^22),(10^25),(10^23)}];

ln=Length[FPpd] (* all the 4 tables have the same length *)

FGd1=-2.426973910644×(10^(-92)); (* deduced above in the non TRL scenario *)

and the sought table with elements   { d,  FGtrl[d]/FG[d] }

FGtrlFG=Table[{Part[Part[FPpd,i],1],(Part[Part[FPpd,i],2]+Part[Part[FPed,i],2]+Part[Part[FEpd,i],2]+Part[Part[FEed,i],2])/FGd1},{i,1,ln,1}];

Prates=ListPlot[FGtrlFG,PlotRange->{{0,1.05*(10^25)},{0.95,1.65}},AxesLabel->{ d,(Subscript[FG, TRL])/FG}]

It is evident from the graphic that the quotient FGtrl[d]/FG[d] increases linearly with the distance d
A simple algebraical function that gives a good fit of such table of quotients FGtrl[d]/FG[d] can easily be found, and is 
g[d] = 1+ k d
where the constant k depends on the values assumed for the constants h1, N0, VM, w1, w2, aI, aS, bI, bS... affecting the forces.
With the values assigned to the constants in the above example, the function g[d] fits the table of quotients if 
 k = 5.6 (10^(-26)) (m^(-1))   as can be seen in the 2 following plots:

k=5.6*(10^(-26));
PQd=Plot[1+k d,{d,0,1.05*(10^25)},AxesLabel->{ d,Q[d]},PlotStyle->Hue[0.82]]

Show[Prates,PQd]

It will then be asserted that:
“to obtain the gravitation force exerted by a distant body #1 of mass m1 on another body #2 of mass m2, the mainstream (Newton) gravitation force   F = G ((m1 m2)/(d^2)) must be multiplied by a (non dimensional) factor g[d] = 1 + k d
where, for the values assigned to the constants in the above examples, the constant k of g[d] must be approximately equal to:  
 k = 5.6 (10^(-26)) (m^(-1))

4)  Rotation curves of galaxies according to the gravitation force FG that disregards the influence of any Tempo Rate Law.

By rotation curve of a galaxy is understood the relation that gives the orbiting speeds v  of the outer stars around the center of a galaxy as a function of their distance d to the center of mass of the galaxy.
In  Newton’s dynamics, that ignores the influence of a Tempo Rate Law, the orbiting speed v of a star in a circular orbit at a distance d from the nucleus of a galaxy of mass M is obtained equating the centripetal (gravitation) force to the centrifugal force suffered by the orbiting star of mass m:
     ((G M m)/(d^2)) = ((m (v^2))/d)  	  =>
   
   =>	      v = (Sqrt[((G M)/d)])

Clear[G,M,m,v,d]

And  assigning to the gravitation constant G and to the mass M of the galaxy the following values (in MKS units):

G=6.67*(10^(-11));
M=3*(10^42);

And  plotting the orbiting speed of the stars   v = (Sqrt[((G M)/d)])  
as a function of d for a wide interval of astronomical distances, gives:

PvNewton=Plot[(Sqrt[((G M)/d)]),{d,(10^23),1*(10^25)},PlotRange->{{(10^23),1.05*(10^25)},{0,20000}},AxesLabel->{ d,v}]

5)  Rotation curves of galaxies according to the gravitation force FGtrl that accounts for the influence of a Tempo Rate Law.

The orbiting speed v of a star in a circular orbit at a distance d from the nucleus of a galaxy of mass M can again be obtained equating the centripetal force (given here by the gravitation force with TRL)  to the centrifugal force suffered by the orbiting star of mass m.
But, as explained above, having found that the function g[d] = 1 + k d  gives a good fit of the values of the above table of quotients FGtrl[d]/FG[d] then the  gravitation force with TRL can be approximated, as explained above, by 
FGTRL =  ((G M m)/(d^2)) g[d]  =  ((G M m)/(d^2)) (1+ k d)
and therefore equating this centripetal force to the centrifugal force:
 ((G M m)/(d^2))(1+ k d) = ((m (v^2))/d)  	  =>
   
   =>	      v = (Sqrt[((G M (1 + k d))/d)])

A Plot of the star speeds  v = (Sqrt[((G M (1 + k d))/d)]) as a function of their distance d to the center of the galaxy is now shown keeping for  the constants G and M the same values assigned above for the plot PvNewton of the Newtonian case
G=6.67*(10^(-11));
M=3*(10^42);
and  assigning  to the constant  k the value
 k = 5.6 (10^(-26)) (m^(-1))    with which g[d] = 1 + k d  gives a good fit of the values of the above table of quotients FGtrl[d]/FG[d]

Clear[d,G,M,k]
vTRL=(Sqrt[((G M (1+k d))/d)]);

G=6.67*(10^(-11));
M=3*(10^42);
k=5.6*(10^(-26));

PvTRL=Plot[vTRL,{d,(10^23),1.15*(10^25)},PlotRange->{{(10^23),1.15*(10^25)},{0,20000}},AxesLabel->{ d,v},PlotStyle->Hue[0.9]]

Show[PvNewton,PvTRL,PlotRange->{{(10^23),1.15*(10^25)},{0,15000}}]