The Equation of Motion for Multiple Charged Particles with Periodic Motion

G.R.Dixon, July 18, 2007

Reference: Programs at the end of this article.

Let us suppose that two identical particles, of mass m and charge q, have the motions

,

.

We assume that the oscillating pair is non-relativistic, say

.

We define the wavelength to be

,

and stipulate that

.

We should like to determine W, the total work per cycle that must be expended in order to drive the pair of charges.

Let us begin by assuming that Newton’s 2^{nd} law applies. That is, we shall *assume* the net forces that must be applied to the particles are mutually

.

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:

.

Or, rearranging slightly,

.

The rate at which F_{Agent} does work is just F_{Agent}v, and thus

.

The work per cycle done by F_{Agent} is then

.

But

.

Thus the work per cycle done by the driving agent is hypothetically

.

In words, W is hypothetically the work per cycle done by the agent’s *counter*action 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!" *(See References) *There is invariably a shortfall, and the shortfall always has the same formula for varying values of w, A and L (provided wA<<c and L~l). It turns out that the following formula is the one that actually satisfies energy conservation:

,

where **S** is the Poynting vector. It is notable that

.

Evidently, in order to satisfy energy conservation, the driving agent must not only counteract F_{Interact}, but it must also apply the force

.

In previous articles **F**_{Radiative} has been defined to be the agent’s *counter*action to a *radiation reaction force* first suggested by Abraham and Lorentz:

.

Summarizing, the driving agent must evidently apply a 3-part force in order for both Newton and energy conservation to be satisfied:

.

Since mav integrates to zero over any cycle time, we find that

,

and it is this formula for W which satisfies

.

It is notable that, whereas F_{Interact} is the force that one charge feels in the electric field of the other charge, F_{RadReact} is the part of the force a charge feels in its *own* field. The other part of this "self" force, often called F_{InertReact} (for inertial reaction force), is just

.

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 just ma - F_{Interact}, can be remedied by including the radiative component -F_{RadReact}:

.

This, then, is the appropriate form of Newton’s 2^{nd} law when 2 (or more) charges interact. If there is only one isolated charge, then there is no interactive force. Newton 2 then becomes

.

And of course if q=0, then F_{RadReact} = 0 and the formula reduces to the familiar

.

Programs for Computing the Work per Cycle and the Energy Flux per Cycle Through an Enclosing Surface

Option Explicit

Private Sub cmdInteractiveWork_Click()

'*********

'Compute W, the work per cycle that must be expended to drive two

'charged particles sinusoidally along the x-axis.

'*********

'Physical and mathematical constants

Const c As Double = 299792000# 'Speed of light

Const epsilon0 As Double = 0.00000000000885 'Permittivity constant

Const pi As Double = 3.14159265358979

Const Steps As Long = 100 'Number of iterations

Const m As Double = 1 'Mass of each particle

Const q As Double = 1 'Charge of each particle

Const Amp As Double = 1 'Amplitude of each motion

Const omega As Double = 0.01 * c / Amp 'Angular frequency of each motion

Const freq As Double = omega / (2 * pi) 'Frequency of each motion

Const tau As Double = 1 / freq 'Period of each motion

Const lambda As Double = c * tau 'Wavelength

Const deltat As Double = tau / Steps 'Time interval between epochs

Const L As Double = 0.5 * lambda 'Separation of 2 charges

'Variables

Dim FNt1(Steps), FNt2(Steps) As Double 'Inertial force on each particle

Dim vx1, vx2 As Double 'Velocity of each particle

Dim accel1, accel2 As Double 'Acceleration of each particle

Dim Px1, Px2 As Double 'Location of each particle

Dim Index As Long 'Loop index

Dim t(Steps) As Double 'Current time

Dim tr As Double 'Retarded time

Dim dt As Double 't - tr

Dim dtmin As Double 'Minimum value for dt

Dim dtmax As Double 'Maximum value for dt

Dim xr1, xr2 As Double 'Retarded quantities

Dim vrx1, vrx2 As Double

Dim arx1, arx2 As Double

Dim Drx1, Drx2 As Double

Dim Dr1, Dr2 As Double

Dim ux1, ux2 As Double

Dim Ex1(Steps), Ex2(Steps) As Double 'Electric fields of particles

Dim Fx1(Steps), Fx2(Steps) As Double 'Electric forces on each particle

Dim FAgent1(Steps), FAgent2(Steps) As Double 'Agent forces

Dim WInt As Double 'Work per cycle to counteract interactive forces

Dim WSelf As Double 'Work per cycle to counteract Self-force for each particle

Dim W As Double 'Total work per cycle

'Compute the inertial forces using Newton 2.

For Index = 0 To Steps - 1

t(Index) = Index * deltat

accel1 = -omega ^ 2 * Amp * Sin(omega * t(Index))

FNt1(Index) = m * accel1

accel2 = -omega ^ 2 * Amp * Sin(omega * t(Index))

FNt2(Index) = m * accel2

Next Index

'Compute the work per cycle to counteract the radiation reaction forces.

WSelf = 2 * (q ^ 2 * omega ^ 3 * Amp ^ 2 / (6 * epsilon0 * c ^ 3))

'Zero out the interactive work variable.

WInt = 0

'Compute the work per cycle to counteract the interactive forces.

For Index = 0 To Steps - 1

'Compute the current positions of q1 and q2.

t(Index) = Index * deltat

Px1 = -L / 2 + Amp * Sin(omega * t(Index))

Px2 = L / 2 + Amp * 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

Do

dt = (dtmax + dtmin) / 2

tr = t(Index) - dt

xr1 = -L / 2 + Amp * Sin(omega * tr)

Drx1 = Px2 - xr1

Dr1 = Abs(Drx1)

If Abs(c * dt - Dr1) < 2 ^ (-30) Then Exit Do

If c * dt - Dr1 > 0 Then

dtmax = dt

Else

dtmin = dt

End If

Loop

vrx1 = omega * Amp * Cos(omega * tr)

arx1 = -(omega ^ 2) * Amp * Sin(omega * tr)

'Then compute the needed retarded quantities for q2.

dtmin = 0

dtmax = 5 * L / c

Do

dt = (dtmax + dtmin) / 2

tr = t(Index) - dt

xr2 = L / 2 + Amp * Sin(omega * tr)

Drx2 = Px1 - xr2

Dr2 = Abs(Drx2)

If Abs(c * dt - Dr2) < 2 ^ (-30) Then Exit Do

If c * dt - Dr2 > 0 Then

dtmax = dt

Else

dtmin = dt

End If

Loop

vrx2 = omega * Amp * Cos(omega * tr)

arx2 = -(omega ^ 2) * Amp * 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(Index) = q / (4 * pi * epsilon0) * Dr1 / ((Drx1 * ux1) ^ 3) * (ux1 * (c ^ 2 - vrx1 ^ 2))

Ex2(Index) = q / (4 * pi * epsilon0) * Dr2 / ((Drx2 * ux2) ^ 3) * (ux2 * (c ^ 2 - vrx2 ^ 2))

Next Index

'Now compute the electric (interactive) forces on q1 and q2,

'and the residue agent forces that must be applied.

For Index = 0 To Steps - 1

Fx1(Index) = q * Ex2(Index) 'Electric force

FAgent1(Index) = FNt1(Index) - Fx1(Index) 'Inertial minus electric force

Fx2(Index) = q * Ex1(Index)

FAgent2(Index) = FNt2(Index) - Fx2(Index)

Next Index

'Compute the work per cycle expended by the residue agent forces.

For Index = 0 To Steps - 1

vx1 = omega * Amp * Cos(omega * t(Index))

vx2 = omega * Amp * Cos(omega * t(Index))

WInt = WInt + FAgent1(Index) * vx1 * deltat + FAgent2(Index) * vx2 * deltat

Next Index

'Now add in the work per cycle that must be expended to counteract the

'radiation reaction force.

W = WInt + WSelf

'Display the result for comparison with the field energy flux per cycle.

MsgBox ("W = " & W)

Stop

End Sub

Private Sub cmdTotalWork_Click()

'**********

'Compute the energy flux per cycle through a spherical surface surrounding

'an oscillating pair of charges, separated by a distance L.

'**********

'Physical and mathematical constants

Const c As Double = 299792000# 'Speed of light

Const epsilon0 As Double = 0.00000000000885 'Permittivity constant

Const pi As Double = 3.14159265358979

Const Steps As Long = 100 'Iterations

Const q As Double = 1 'Charge of each particle

Const Amp As Double = 1 'Amplitude of oscillation

Const omega As Double = 0.01 * c / Amp 'Angular frequency

Const freq As Double = omega / (2 * pi) 'Frequency

Const tau As Double = 1 / freq 'Period

Const deltat As Double = tau / Steps 'Time interval between epochs

Const lambda As Double = c * tau 'Wavelength

Const L As Double = 0.5 * lambda 'Separation

Const dtheta As Double = pi / Steps 'Interval between angles over sphere

Const Radius As Double = 3 * lambda 'Spherical radius

'Variables

Dim theta As Double 'Angle between x-axis and point on sphere (in xy-plane)

Dim Py, Px As Double 'Coordinates of point on sphere

Dim TimeIndex, ThetaIndex As Long 'Loop indexes

Dim t(Steps) As Double 'Current time

Dim tr As Double 'Retarded time

Dim dt As Double 't - tr

Dim dtmin As Double 'Minimum value for dt

Dim dtmax As Double 'Maximum value for dt

Dim xr1, xr2 As Double 'Retarded positions (on x-axis)

Dim vrx1, vrx2 As Double 'Retarded velocities

Dim arx1, arx2 As Double 'Retarded accelerations

Dim Drx1, Drx2 As Double 'x-components of vectors Dr

Dim Dry1, Dry2 As Double 'y-components of vectors Dr

Dim Dr1, Dr2 As Double 'magnitudes of vectors Dr

Dim ux1, ux2 As Double 'x-components of vectors u

Dim uy1, uy2 As Double 'y-components of vectors u

Dim Ex1, Ex2, Ex As Double 'x-components of electric field vectors

Dim Ey1, Ey2, Ey As Double 'y-components of electric field vectors

Dim Bz1, Bz2, Bz As Double 'z-components of magnetic field vectors

Dim Sx As Double 'x-component of Poynting vector at point on sphere

Dim Sy As Double 'y-component of Poynting vector

Dim unitx, unity, dA As Double 'unit vectors and spherical area increment

Dim EFluxRate(Steps), EFlux(Steps) As Double 'Enery flux at different angles

Dim Snormal As Double 'normal component of Poynting vector

Dim Erad As Double 'energy flux per cycle at separation L

Erad = 0

'Compute Erad

'Zero out the energy fluxes for each epoch.

For TimeIndex = 0 To Steps - 1

EFlux(TimeIndex) = 0

Next TimeIndex

'Compute the energy flux through the spherical surface, for each time

'interval.

For TimeIndex = 0 To Steps - 1 'time epochs

Debug.Print TimeIndex

t(TimeIndex) = TimeIndex * deltat

'Zero out the energy flux rates for each theta interval.

For ThetaIndex = 0 To Steps - 1

EFluxRate(ThetaIndex) = 0

Next ThetaIndex

'Compute the energy flux rate through each theta interval.

For ThetaIndex = 0 To Steps - 1

theta = ThetaIndex * dtheta + dtheta / 2

unitx = Cos(theta)

unity = Sin(theta)

Px = Radius * Cos(theta)

Py = Radius * Sin(theta)

dA = 2 * pi * Py * Radius * dtheta

'First do q1.

dtmin = 0

dtmax = 2 * Radius / c

Do

dt = (dtmax + dtmin) / 2

tr = t(TimeIndex) - dt

xr1 = -L / 2 + Amp * Sin(omega * tr)

Drx1 = Px - xr1

Dry1 = Py

Dr1 = Sqr(Drx1 ^ 2 + Dry1 ^ 2)

If Abs(c * dt - Dr1) < 2 ^ (-30) Then Exit Do

If c * dt - Dr1 > 0 Then

dtmax = dt

Else

dtmin = dt

End If

Loop

vrx1 = omega * Amp * Cos(omega * tr)

arx1 = -(omega ^ 2) * Amp * Sin(omega * tr)

'Then do q2.

dtmin = 0

dtmax = 2 * Radius / c

Do

dt = (dtmax + dtmin) / 2

tr = t(TimeIndex) - dt

xr2 = L / 2 + Amp * Sin(omega * tr)

Drx2 = Px - xr2

Dry2 = Py

Dr2 = Sqr(Drx2 ^ 2 + Dry2 ^ 2)

If Abs(c * dt - Dr2) < 2 ^ (-30) Then Exit Do

If c * dt - Dr2 > 0 Then

dtmax = dt

Else

dtmin = dt

End If

Loop

vrx2 = omega * Amp * Cos(omega * tr)

arx2 = -(omega ^ 2) * Amp * 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 * pi * epsilon0) * Dr1 / (Drx1 * ux1 + Dry1 * uy1) ^ 3 * (ux1 * (c ^ 2 - vrx1 ^ 2) + Dry1 * (-uy1 * arx1))

Ey1 = q / (4 * 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 * pi * epsilon0) * Dr2 / (Drx2 * ux2 + Dry2 * uy2) ^ 3 * (ux2 * (c ^ 2 - vrx2 ^ 2) + Dry2 * (-uy2 * arx2))

Ey2 = q / (4 * 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.

EFluxRate(ThetaIndex) = EFluxRate(ThetaIndex) + Snormal * dA

'Repeat for the next angle theta.

Next ThetaIndex

'Now update the running sum of energy fluxes.

For ThetaIndex = 0 To Steps - 1

EFlux(TimeIndex) = EFlux(TimeIndex) + EFluxRate(ThetaIndex) * deltat

Next ThetaIndex

Next TimeIndex

'Once the energy fluxes for each time interval have been computed, add

'the energy fluxes to Erad for the current separation.

For TimeIndex = 0 To Steps - 1

Erad = Erad + EFlux(TimeIndex)

Next TimeIndex

'Display the results for comparison with W.

MsgBox ("Erad = " & Erad)

End Sub