The Equation of Motion for Multiple Charged Particles with Periodic Motion Abstract. An argument is made for the reality of the radiation traction force. Let us suppose that two identical particles, of mass m and charge q, have the motions ,    (1) .    (2) We assume that the oscillating pair is non-relativistic, say .    (3) We define the wavelength to be ,    (4) and stipulate that .    (5) We should like (1) to determine W, the total work per cycle that must be expended by a driving agent to maintain the oscillation, and (2) to compare W to Erad, the total radiant energy per cycle that fluxes through a spherical surface containing the oscillating pair. Let us begin by assuming that Newton’s 2nd law applies. That is, we shall assume that the net forces which must be applied to the particles are mutually .    (6) Part of this force will be supplied by the interactive force that one charge experiences in the electric field of the other. And the remainder must presumably be supplied by the driving agent. For either charge, .    (7) Or, rearranging slightly, .    (8) The rate at which FAgent does work is just FAgentv, and thus .    (9) The work per cycle done by FAgent is then .    (10) But .    (11) Thus the work per cycle done by the driving agent is hypothetically .    (12) In words, W is hypothetically the work per cycle done by the agent’s counteraction to the interactive force. A question suggested by Energy Conservation is: "Does this formula for W equal the field energy flux per cycle through an enclosing surface?" It turns out that the answer is "No!" There is a shortfall, and the shortfall always has the same formula for varying values of w, A and L. It turns out that the following formula is the one that actually satisfies energy conservation: ,    (13) where S is the Poynting vector. It is notable that .    (14) Evidently, in order to satisfy energy conservation, the driving agent must not only counteract FInteract, but it must also apply the force .    (15) Now FRadiative is defined to be the agent’s counteraction to a radiation reaction force first suggested by Abraham and Lorentz: .    (16) Summarizing, the driving agent must evidently apply a 3-part force in order for both Newton and energy conservation to be satisfied: .    (17) Since mav integrates to zero over any cycle time, we find that ,    (18) and it is this formula for W which satisfies .    (19) It is notable that, whereas FInteract is the force that one charge feels in the electric field of the other charge, FRadReact is the part of the force a charge feels in its own field. The other part of this "self" force, often called FInertReact (for inertial reaction force), is just .    (20) This part of the self-force does not contribute anything when integrated over a complete cycle time. Evidently the shortfall, incurred by assuming that the net agent force is simply ma - FInteract, can be remedied by including the radiatiation component -FRadReact: .    (21) This, then, is the appropriate form of Newton’s 2nd law when two (or more) charges interact. If there is only one isolated charge, then there is no interactive force. Newton 2 then becomes .    (22) And of course if q=0, then FRadReact = 0 and the formula reduces to the familiar .    (23) Python programs for computing (a) W, the Work per cycle, and (b) Erad, the energy flux per cycle through an enclosing surface, follow. The computed results are W = 49112 joules, Erad = 49115 joules. These practically equal results appear to confirm that the driving agent must counteract a radiation reaction force, in addition to an interactive force, in order to maintain the oscillatory motion.  Program A. #********* #Compute W, the work per cycle that must be expended to drive two #charged particles sinusoidally along the x-axis. #********* import math c = 299792000# #Speed of light epsilon0 = 0.00000000000885 #Permittivity constant Steps = 500 #Number of iterations q = 1. #Charge of each particle Amp = 1. #Amplitude of each motion omega=.01*c freq=omega/(2*math.pi) tau = 1 / freq #Period of each motion wavelngth = c * tau #Wavelength deltat = tau / Steps #Time interval between epochs L=.5*wavelngth t=[] #Current time m=1. Ex1=[] Ex2=[] #Electric fields of particles Fx1=[] Fx2=[] #Electric forces on each particle FAgent1=[] FAgent2=[] #Agent forces FNt1=[] FNt2=[] #Compute the work per cycle to counteract the radiation reaction forces. WSelf = 2. * (q ** 2 * omega ** 3 * Amp ** 2 / (6 * epsilon0 * c ** 3)) print('WSelf=',WSelf) input('OK?') #Zero out the interactive work variable. WInt = 0. #Compute the work per cycle to counteract the interactive forces. t.clear() for Index in range(Steps):     #Compute the current positions of q1 and q2.     t.append(Index * deltat)     Px1 = -L / 2 + Amp * math.sin(omega * t[Index])     Px2 = L / 2 + Amp * math.sin(omega * t[Index])     #Compute the electric field of q1, at the position of q2, and     #the electric field of q2 at the position of q1.     #First compute the needed retarded quantities for q1.     dtmin = 0     dtmax = 5 * L / c     while(1==1):         dt = (dtmax + dtmin) / 2         tr = t[Index] - dt         xr1 = -L / 2 + Amp * math.sin(omega * tr)         Drx1 = Px2 - xr1         Dr1 = abs(Drx1)         if(abs(c * dt - Dr1) < 2 ** (-30)):             break         if((c * dt - Dr1) > 0.):             dtmax = dt         else:             dtmin = dt     vrx1 = omega * Amp * math.cos(omega * tr)     arx1 = -(omega ** 2) * Amp * math.sin(omega * tr)     #Then compute the needed retarded quantities for q2.     dtmin = 0.     dtmax = 5 * L / c     while(1==1):         dt = (dtmax + dtmin) / 2         tr = t[Index] - dt         xr2 = L / 2 + Amp * math.sin(omega * tr)         Drx2 = Px1 - xr2         Dr2 = abs(Drx2)         if(abs(c * dt - Dr2) < 2 ** (-30)):             break         if((c * dt - Dr2) > 0):             dtmax = dt         else:             dtmin = dt     vrx2 = omega * Amp * math.cos(omega * tr)     arx2 = -(omega ** 2) * Amp * math.sin(omega * tr)     #Compute the components of vectors u.     ux1 = c * Drx1 / Dr1 - vrx1     ux2 = c * Drx2 / Dr2 - vrx2     #Compute the electric field components     Ex1.append(q / (4 * math.pi * epsilon0) * Dr1 / ((Drx1 * ux1) ** 3) * (ux1 * (c ** 2 - vrx1 ** 2)))     Ex2.append(q / (4 * math.pi * epsilon0) * Dr2 / ((Drx2 * ux2) ** 3) * (ux2 * (c ** 2 - vrx2 ** 2)))     #Now compute the electric (interactive) forces on q1 and q2,     #and the residue agent forces that must be applied. for Index in range(Steps):     Fx1.append(q * Ex2[Index]) #Electric force     FAgent1.append(- Fx1[Index])     Fx2.append(q * Ex1[Index])     FAgent2.append(- Fx2[Index])     #Compute the work per cycle expended by the residue agent forces. for Index in range(Steps):     vx1 = omega * Amp * math.cos(omega * t[Index])     vx2 = omega * Amp * math.cos(omega * t[Index])     WInt += (FAgent1[Index] * vx1 * deltat) + (FAgent2[Index] * vx2 * deltat) #Now add in the work per cycle that must be expended to counteract the #radiation reaction force. W = WInt + WSelf print('WInt=',WInt,', WSelf=',WSelf) #Display the result for comparison with the field energy flux per cycle. print('W = ',W) print('End of Program') Program B. #********** #Compute the energy flux per cycle through a spherical surface surrounding #an oscillating pair of charges, separated by a distance L. #********** import math c = 299792000# #Speed of light epsilon0 = 0.00000000000885 #Permittivity constant Steps = 500 #Iterations q = 1. #Charge of each particle Amp = 1. #Amplitude of oscillation omega=.01*c freq=omega/(2*math.pi) tau = 1. / freq #Period deltat = tau / Steps #Time interval between epochs wavelngth = c * tau #Wavelength L = .5*wavelngth #Separation dtheta = math.pi / Steps #Interval between angles over sphere Radius = 3 * wavelngth #Spherical radius t=[] #Current time EFluxRate=[] EFlux=[] #Enery flux at different angles Erad = 0. EFlux.clear() #Compute the energy flux through the spherical surface, for each time #interval. for TimeIndex in range(Steps):     #EFlux.append(TimeIndex)     t.append(TimeIndex * deltat) for TimeIndex in range(Steps): #time epochs     #Compute the energy flux rate through each theta interval.     EFluxRateTemp=0.     for ThetaIndex in range(Steps):         theta = ThetaIndex * dtheta + dtheta / 2         unitx = math.cos(theta)         unity = math.sin(theta)         Px = Radius * math.cos(theta)         Py = Radius * math.sin(theta)         dA = 2 * math.pi * Py * Radius * dtheta         #First do q1.         dtmin = 0.         dtmax = 2 * Radius / c         while(1==1):             dt = (dtmax + dtmin) / 2             tr = t[TimeIndex] - dt             xr1 = -L / 2 + Amp * math.sin(omega * tr)             Drx1 = Px - xr1             Dry1 = Py             Dr1 = math.sqrt(Drx1 ** 2 + Dry1 ** 2)             if( abs(c * dt - Dr1) < 2 ** (-30)):                 break             if( c * dt - Dr1 > 0):                 dtmax = dt             else:                 dtmin = dt         vrx1 = omega * Amp * math.cos(omega * tr)         arx1 = -(omega ** 2) * Amp * math.sin(omega * tr)         #Then do q2.         dtmin = 0.         dtmax = 2 * Radius / c         while(1==1):             dt = (dtmax + dtmin) / 2             tr = t[TimeIndex] - dt             xr2 = L / 2 + Amp * math.sin(omega * tr)             Drx2 = Px - xr2             Dry2 = Py             Dr2 = math.sqrt(Drx2 ** 2 + Dry2 ** 2)             if( abs(c * dt - Dr2) < 2 ** (-30)):                 break             if( c * dt - Dr2 > 0):                 dtmax = dt             else:                 dtmin = dt         vrx2 = omega * Amp * math.cos(omega * tr)         arx2 = -(omega ** 2) * Amp * math.sin(omega * tr)         #Compute the components of vectors u.         ux1 = c * Drx1 / Dr1 - vrx1         uy1 = c * Dry1 / Dr1         ux2 = c * Drx2 / Dr2 - vrx2         uy2 = c * Dry2 / Dr2         #Compute the electric and magnetic field components         Ex1 = q / (4 * math.pi * epsilon0) * Dr1 / (Drx1 * ux1 + Dry1 * uy1) ** 3 * (ux1 * (c ** 2 - vrx1 ** 2) + Dry1 * (-uy1 * arx1))         Ey1 = q / (4 * math.pi * epsilon0) * Dr1 / (Drx1 * ux1 + Dry1 * uy1) ** 3 * (uy1 * (c ** 2 - vrx1 ** 2) - Drx1 * (-uy1 * arx1))         Bz1 = 1 / (c * Dr1) * (Drx1 * Ey1 - Dry1 * Ex1)         Ex2 = q / (4 * math.pi * epsilon0) * Dr2 / (Drx2 * ux2 + Dry2 * uy2) ** 3 * (ux2 * (c ** 2 - vrx2 ** 2) + Dry2 * (-uy2 * arx2))         Ey2 = q / (4 * math.pi * epsilon0) * Dr2 / (Drx2 * ux2 + Dry2 * uy2) ** 3 * (uy2 * (c ** 2 - vrx2 ** 2) - Drx2 * (-uy2 * arx2))         Bz2 = 1 / (c * Dr2) * (Drx2 * Ey2 - Dry2 * Ex2)         #Sum the field components for the 2 charges.         Ex = Ex1 + Ex2         Ey = Ey1 + Ey2         Bz = Bz1 + Bz2         #Compute the Poynting vector components, using the net E and B         #field components.         Sx = epsilon0 * c ** 2 * (Ey * Bz)         Sy = epsilon0 * c ** 2 * (-Ex * Bz)         #Find the Poynting vector component that is normal to the         #spherical surface.         Snormal = Sx * unitx + Sy * unity         #Use it to update the energy flux from this value of theta.         EFluxRateTemp+=Snormal*dA     EFlux.append(EFluxRateTemp*deltat #Once the energy fluxes for each time interval have been computed, add #the energy fluxes to Erad for the current separation. for TimeIndex in range(Steps):     Erad += EFlux[TimeIndex] #Display the results for comparison with W. print('Erad = ',Erad) print('End of program')