Solutions of the Wheeler-Feynman equations with discontinuous velocities

We generalize Wheeler-Feynman electrodynamics with a variational boundary-value problem with past and future boundary segments that can include velocity discontinuity points. Critical-point trajectories must satisfy the Euler-Lagrange equations of the action functional, which are neutral-differential delay equations of motion (the Wheeler-Feynman equations of motion). At velocity discontinuity points, critical-point orbits must satisfy the Weierstrass-Erdmann conditions of continuity of partial momenta and partial energies. We study a special class of boundary data having the shortest time-separation between boundary segments, for which case the Wheeler-Feynman equations reduce to a two-point boundary problem for an ordinary differential equation. For this simple case we prove that the extended variational problem has solutions with discontinuous velocities. We construct a numerical method to solve the Wheeler-Feynman equations together with the Weierstrass-Erdmann conditions and calculate some numerical orbits with discontinuous velocities.


I. INTRODUCTION
To overcome the lack of equations of motion for point charges in Maxwell's electrodynamics, J. Wheeler and R. Feynman 1,2 developed an electrodynamics based on the minimization of the Fokker-Schwarzschild-Tetrode [3][4][5] action to replace Maxwell's equations. Besides having sensible equations of motion for point charges, the action functional seemed to point towards a canonical quantization 6,7 of the two-body problem, so far an unfulfilled promise 8 .
In a recent development, Wheeler-Feynman electrodynamics was embedded in a variational boundary-value problem 9 , henceforth variational electrodynamics 10,11 . The development of Ref. 9 was followed by a study of the neutral differential delay equations of motion of the electromagnetic two-body problem 10,11 . Consistent application of variational electrodynamics defines trajectories with discontinuous velocities among the critical points of the variational problem 10,11 , henceforth broken extrema 12 . The conditions for a piecewise C 2 extremal trajectory are (i) to satisfy the Euler-Lagrange equations along C 2 segments, which are neutral differential delay equations, henceforth the Wheeler-Feynman equations of motion, and (ii) at velocity disa) author's email address: daniel.souza@usp.br b) author's email address: jayme.deluca@gmail.com; webpage https://sites.google.com/site/jaymedeluca/ continuity points (henceforth breaking points), a broken extremum must further satisfy the Weierstrass-Erdmann continuity conditions for the partial momenta and partial energies 11,12 , a system of four nonlinear equations at each breaking point 13 . It is known that neutral differential delay equations can generate breaking points 13 which are propagated along trajectories 14 . As regards some physical property of the two-body problem, it can be shown that globally bounded two-body trajectories with vanishing far-fields must have discontinuous velocities 10 .
The electromagnetic variational problem 9 has boundary conditions in past and future, which make it hard to study numerically. Here we study the simplest boundary value problem, having the shortest possible timedistance between past and future boundary segments, as explained in the following. For these, the variational problem reduces to a two-point boundary problem 11 (i.e. a shooting problem for an ODE integrator 15 ). We show well-posedness and calculate some numerical solutions of the variational problem with shortest-length segments near circular orbits of large radii 16,17 . This paper is divided as follows: In Section II we introduce the variational boundary value problem. In Section III we explain the Weierstrass-Erdmann corner conditions for extrema with discontinuous velocities of the boundary-value problem with shortest-length boundary segments. In Section IV we develop the numerical method for shortest-type boundaries and calculate some numerical trajectories. In Section V we put the discussions and conclusion.

II. VARIATIONAL BOUNDARY VALUE PROBLEM
Here we consider the natural electronic units where the speed of light, electronic charge and electronic mass are c = −e 1 = m 1 ≡ 1. We henceforth use the index i = 1 to indicate electronic quantities and the index i = 2 to denote quantities of the positive charge 2.
The variational problem is to find (extended) trajectories given by continuous and piecewise C 2 functions of a (real) parameter s → Henceforth a dot over vectors of R 4 denotes derivative respect to the parameter s, i.e.,ẋ i ≡ dx i ds = (ṫ i (s),ẋ i (s)) for i = 1, 2. On regular segments the particle velocities are recovered by the chain rule, dx i dt i =ẋ i ṫ i . The space R 4 is equipped with the usual Euclidean norm of R 4 (double bars), defined by the inner product, i.e., . A form that appears often is the (real) Minkowski bilinear form between two vectors v 1 ≡ (a 1 , b 1 ) and v 2 ≡ (a 2 , b 2 ), defined by the scalar product with the other vector's dual, i.e., Central to the construction of the variational problem are the light-cone conditions for j ≡ 3 − i and i = 1, 2, which are implicit conditions for the trajectories. In Eq. (1) the plus sign defines the future light-cone condition for particle i, while the minus sign defines the past light-cone condition for particle i. In Ref. 19 it is shown that if trajectories are sub-luminal, i.e., for i = 1, 2, the light-cone conditions (1) have unique solutions s j± (s). Henceforth a ± sign after the particle index j ≡ 3 − i indicates a quantity of particle j evaluated at the advanced/retarded argument s j± (s). The variational two-body problem 9 is defined by the critical-point conditions of the action functional where and the interaction energy with i = 1, 2 and j ≡ 3 − i, which has a non-zero denominator along non-collisional trajectories 19 . Henceforth we specify for the attractive problem by setting e i e j = −1 in Eq. (5). The first variation of (3) naturally decomposes in a sum of two partial variations, δS = δS 1 + δS 2 , as follows; (i) for variation δS 1 trajectory 1 varies for while trajectory 2 is kept fixed, and therefore the first term on the right-hand-side of (3) is a constant term, I 1 . The remaining three integrals are over ds 1 , thus defining the partial Lagrangian L 1 as the integrand, and (ii) along variation δS 2 trajectory 2 varies for s 2 ∈ [s O + , s L B ] while trajectory 1 is kept fixed. To calculate variation δS 2 it is convenient to express functional (6) with three integrals over ds 2 , plus a constant term I 2 , as obtained using a change of variable on the last two integrals of Eq. (3) to the other particle's parameter in light-cone condition (1) 11 . Since the problem thus defined is totally symmetric, bellow we explain the critical point conditions and equations of motion for particle 1 only. Form (6) with its integral over partial Lagrangian L 1 is used to calculate the partial variation δS 1 and the Euler-Lagrange equations of motion of particle 1. The boundary conditions in past and future are described in Fig. 1 Since past and future boundary segments are supposed to be independent of each other, we must have that the past boundary-segment of particle 2 does not interact with the future boundary-segment of particle 1, in which minimal case point L − is in the forward light-cone condition of point O + .
Next we discuss the variational problem for piecewise C 2 continuous trajectories x 1 (s), x 2 (s) having monotonically increasing time-components and satisfying the above boundaries conditions. To calculate the first variation δS 1 we assume trajectory x i (s) and its continuous perturbation b 1 (s) are C 2 inside the intervals s i ∈ (s µ−1 , s µ ) defined by the grid of possible discontinuities s µ with µ = 1, ..., N . The perturbed trajectory is defined by and outside the grid of discontinuity points, s ≠ s µ , and for i = 1 satisfy the fixed-ends boundary conditions The first variation δS i induced by a trajectory variation (7) about a non-collisional sub-luminal trajectory is 19 where b i is the norm of piecewise C 2 perturbations, Integrating (10) by parts over each C 2 segment to eliminate the integral containingḃ i (s) yields Since the b i (s) are continuous and vanish at the endpoints, the second term of Eq. (12) can be re-arranged as where is the (possible) discontinuity of partial momentum i. For a critical point we must have δS i = 0 for arbitrary b i (s), and since the first term on the right-handside of (13) is an integral and the second term depends on the discrete values of b i (s µ ) on the finite number of grid points, each must vanish independently, yielding (a) Euler-Lagrange equations piecewise, henceforth the Wheeler-Feynman equations of motion, which can be written for the spatial components as 9 where the Liénard-Wiechert fields are defined by with Still in Eq. (17), γ j± ∶= (1 − v 2 j± ) −1 2 while r ij± is the distance in light-cone, r ij± ∶= x i −x j± , and the unit vector n ij± ∶= (x i − x j± ) r ij± points from the advanced/retarded position x j± to the position x i . The vanishing of the second term on the right-hand-side of (13) imposes four continuity conditions at each grid point, henceforth the Weierstrass-Erdmann corner conditions 12 of continuity of partial momenta and partial energies Using definition (19) with Eq. (3) yields the partial momentum and the partial energy Defining the right and left limits of the partial momenta/energies at each breaking point by P r 1 , P l 1 , E r 2 and E l 2 , respectively, the Weierstrass-Erdmann corner conditions can be expressed at each grid point s = s µ by

III. SHOOTING PROBLEM AND WEIERSTRASS-ERDMANN CORNER CONDITIONS
The shortest-length boundary value problem occurs when event L − is in the future light-cone of event O + , as illustrated in Fig. 2. Again, for boundaries with a smaller than the minimum time-separation illustrated in Fig. 2, the boundary-segments would interact in lightcone, an absurd.
For shortest-length boundary conditions, points x 2− and x 1+ fall each on a past/future boundary segment (illustrated in red in Fig. 2), and therefore are given functions of the running positions x 2 and x 1 , thus reducing the Wheeler-Feynman equations to a two-point boundary problem for an ODE, as explained in the following.
The Wheeler-Feynman equations of motion (15) can be expressed as 11  where the fields E ± j can be re-arranged using (18) as The linear dependence of the E ± j on the accelerations a j± allows the use of a (simpler) matrix form where (26) Further using (18), Eq. (23) thus becomes where We can isolate the accelerations a 1 and a 2 from Eqs. (23) with i = 1, 2, yielding an algebraic-differential equation where we have defined the matrices and M i ∶= m i γ i , A + 12 ∶= K +12 B 12+ and A − 21 ∶= K −21 B 21− . If one can invert the matrices on the left-hand-side of (29) locally, the algebraic differential equation (29) reduces to an ODE.
In order to be able to invert the matrices on the lefthand-side of (29), we restrict to boundary-segments satisfying For example, along small perturbations of segments of circular orbits with large radii, r ij± ≫ 1 16 , condition (32) holds.
Notice that the running accelerations in (29) are each defined respect to a different independent time variable, i.e., a 1 ≡ d 2 x 1 dt 2 1 and a 2 ≡ d 2 x 2 dt 2 2+ , where t 2+ (t 1 ) is defined by the implicit-function theorem and the future light-cone relation (1). The derivative dt 2+ dt 1 is obtained by taking a derivative of (1) piecewise, yielding Last, if (32) holds we can define the vector fields and transform (29) into the following non-autonomous ODE with the two-point boundary conditions thus defining a two-point boundary problem. The factor λ 12+ in Eq. (35) insures the running positions satisfy the t 1 -explict condition (1), thus arriving at the end-point with x 2 (t L B ) in the future light-cone of x 1 (t L − ) (the second column of boundary condition (36)). We solve the two-point boundary problem (35) and (36) with a shooting method that searches initial veloci- x O − such that the initial value problem terminates at the specified end-points To define the map for the shooting method, we consider that position X ≡ [x 1 , , and linearize about some reference initial velocity, V 0 ≡ [v 10 , v 20 ] ⊺ , yielding Restricted to small perturbations of circular boundary segments, the velocity is a constant at O(1 r 12 ). The shooting method is used with (37) to numerically calculate matrix J s0 and vector X 0 for the perturbed boundary data. For that we solve seven initial value problems (35) for the shooting map (37). The seven initial velocities V are used with the same small perturbation of boundary segments from a large radius circular orbit, defined as follows; (a) we start with the velocity of the unperturbed circular orbit, for which the second term on the right-hand-side of (37) vanishes, thus calculating the constant X 0 within the numerical precision, and (b) we perturb each of the six components of the initial velocity away from the circular orbit's velocity V 0 , one component at a time. We further solve Eq. (37) for V, substitute V 0 by V and solve again the seven initial value problems (35) to find the new X 0 and J s0 . This iterative process results in the following map with k = 1, 2, ... and X k = [x L − , x L B ] ⊺ . If for each iteration k the matrix J sk is well conditioned and the map (38) converges, these velocities solve the two-point boundary problem given by (35) and (36) within the numerical error. The condition number of the shooting matrix depends on the stability of the initial value problem 15 and with generic boundary segments one may not be able either to invert matrix J sk or to find a unique solution or any solution for map (38). As we show next, for small perturbations of circular orbits of large radii the boundary value problem is well-posed in a local subspace of boundary segments, matrix J sk is well conditioned and iteration (38) converges. In general, we expect the orbital velocities at points O + and L − to be different from the velocity on the boundary segments, and thus discontinuous.
In Theorem 1 we analyze the condition of matrix J sk and the convergence of map (38) for boundary-segments near segments of circular orbits of large radii. Theorem 1 Let x S i (t) ∈ R 3 and v S i (t) ∈ R 3 , denote the positions and velocities along a doubly circular orbit 16 , and x h i (t) ∈ R 3 be the boundary-segments for i = 1, 2. Assume the x h i (t) are C 2 and in a neighborhood of circular orbits with large radii, r 12 ≫ 1. Then the two-point boundary problem posed by ODE (35) with boundary conditions (36) has a unique solution that can be found using the shooting method.
Proof. The positions are given by the integral and for boundary segments near circular orbits with large radii, the velocity V(t) is almost constant because the acceleration falls at the most as 1 r 12 , as can be found by inspecting (23) and (24). Notice that for near-circular boundary data of the shortest type the time span of the circular flight is so short that the trajectories are almost straight lines.
Using the above we have the approximation and Eqs. (40) and (41) yield Comparing (42) with (37) yields (37) with J s0 = ∆t1, where ∆t is the time span and 1 the 6×6 identity matrix. Therefore matrix J s0 has a well-conditioned inverse and the unique solution depends continuously on the boundary segments and initial velocities, i.e., a well-posed shooting method (38). ∎ When the continuous boundary-segments are only piecewise C 2 and have velocity discontinuities, one has to stop the integration of ODE (35) at all breaking points to satisfy the Weierstrass-Erdmann corner conditions (22), as follows.
For boundary segments with pre-specified jumps, the Weierstrass-Erdmann corner conditions (22) form an overdetermined system of 8 equations for the 6 velocities v r 1 , v r 2 to continue the orbit on the right-hand-side of each discontinuity point. The former shows that for generic boundary segments the solution may not even exist.
In order to describe boundary segments that do have solutions with nontrivial discontinuous velocities, we include as variables two components of each right-velocity along each boundary segment, v r 2− , v r 1+ , which have not been used until the corner point. Next we show that such augment of the set of variables must be made very carefully.
The above described augment generates an underdetermined nonlinear system having 8 equations and 12 variables, where the superscripts l , r denote the left-limit and rightlimit of the velocities at the breaking points. Defining the velocity jumps ∆v i ∶= v r i − v l i , the vector of jumps v ∶= (∆v 1 , ∆v 2 , ∆v 1+ , ∆v 2− ) ⊺ and the local value of the Eq. (43) as Equations (20) and (21) Definingf 0 ≡ d−F(v 0 ), we obtain at lowest order a linear system of 8 equations for 12 variables, Next we analyze the linear system (46) for near-circular boundary segments of large radii. The idea is to choose four of the twelve variables on the left-hand-side of (46) as independent variables to be placed on the right-hand side of (46). The system thus generated should have an invertible linear 8 × 8 matrix with maximum row rank on the left-hand-side, yielding a unique solution for the eight "slave variables" on the left-hand-side, i.e., Vector f 0 ≡ d−F(v 0 ) on the right-hand-side of (47) should depend on the four independent variables, and for a nontrivial discontinuity, f 0 must be nonzero. We solve (46) for v using the following iterative process. Starting from the solution v of linear system (47), we replace v 0 by v and recalculate J 0 (v 0 ) and f to find the next iterate v by (47), thus generating the following map with k = 0, 1, 2, . . . and f k ≡ d − F(v k ). Theorem 2 gives an example where linear system (47) has a well conditioned matrix J 0 on the left-hand-side and a nonzero f 0 , thus determining a unique solution to (43) by iterating map (48).
, v S i (t)) ∈ R 3 denote positions and velocities along a doubly circular orbit 16 , and (x h i (t), v h i (t)) ∈ R 3 be the C 2 boundary-segments in a neighbourhood of circular orbits with large radii, r 12 ≫ 1, for i = 1, 2 and with v h 2 (t) along thex axis and v h 1 (t)along the −x axis. Then we can choose the x and z components of ∆v h 1 (t) and ∆v h 2 (t) as independent variables. The reduced linearized problem for the remaining (slave) variables v ≡ (∆v 1 , ∆v 2 , ∆v y 1+ , ∆v y 2− ) ⊺ is an inhomogeneous linear system given by (48) with a well-conditioned matrix J k , yielding a unique solution to Eq. (43).
Proof. For the above described near-circular orbits with a large inter-particle separation, r 12± ≡ r ≫ 1, and low velocities, v r The linearized expansion (47) evaluated in the limit of large radii where trajectories are approximated by straight lines near v 0 = 0 yields f 0 = − 1 2r (∆v x 2− , 0, ∆v z 2− , ∆v x 1+ , 0, ∆v z 1+ , 0, 0) ⊺ and Matrix J 0 is well-conditioned and one can calculate the next iterate v k+1 using map (48) with k = 1. For a continuous dependence on boundary data, velocity discontinuities in histories can be chosen arbitrarily in the subspace of independent discontinuity variables, (∆v x 1+ , ∆v x 2− , ∆v z 1+ , ∆v z 2− ). Thus restricted, matrix J k is invertible and the unique solution for the slave variables v ≡ (∆v 1 , ∆v 2 , ∆v y 1+ , ∆v y 2− ) ⊺ depends continuously on the boundary data and on (∆v x 1+ , ∆v x 2− , ∆v z 1+ , ∆v z 2− ) and the unique solution to (43) is given by the fixed point of map (48). ∎

IV. NUMERICAL EXPERIMENTS
The family of circular orbits 16 with angular velocity ω and radius r i = v i ω, with i = 1, 2, can be parametrized by the retardation angle 11 θ. The light-cone time τ for light to travel the inter-particle distance is related to the constant retardation angle of the circular orbit by 11 τ = θ ω. In the limit of small θ, the circular radii are given by r i = 1 (m i θ 2 ) for i = 1, 2 and the constant angular frequency is ω = m 1 m 2 θ 3 (m 1 + m 2 ) 11 .
Our first numerical experiment uses boundary data given by a continuous perturbation of circular orbit's segments with θ = 0.77 for m 1 = 1 and m 2 = 2, henceforth boundary data (I). The perturbed boundary segments include a velocity discontinuity in the middle of the histories, i.e., are only piecewise C 2 . The numerical solution is calculated with a shooting method that solves each initial value problem (35) with a fourth-order Runge-Kutta method, as described in Section (III). When the integration reaches the breaking point the Runge-Kutta integrator is halted and we solve the Weierstrass-Erdmann corner conditions using the linear solution (47) as initial guess for the function fsolve of MatLab R2011a.
In Fig. 3 we show the trajectories of the particles, history segments in red and numerically calculated trajectories in black and blue lines. In Fig. 4 we show the components of the velocity of particle 1 and Fig. 5 shows the components of the velocity of particle 2. Notice that the numerically calculated solutions have discontinuous velocities at points O + and L − (which is a generic feature of the shortest-lenght boundaries) and at one extra pair of points in light-cone and along each trajectory as caused by the breaking point in histories. In our second numerical experiment the boundary segments are given by a C 2 perturbation of circular orbit's segments with θ = 0.077 for m 1 = 1 and m 2 = 10, henceforth boundary data (II). The perturbed boundary segments have no velocity discontinuity in the history segments. The numerical solution in again calculated with a shooting method that solves each initial value problem (35) with a fourth-order Runge-Kutta method, as described in Section (III). In Fig. 6 we show the trajectories of the particles, history segments in red and numerically calculated trajectories in black and blue lines. Last, Fig.  7 shows the components of the velocity of particle 1 and Fig. 8 shows the components of the velocity of particle 2. Notice that for our second experiment the numerically calculated solutions have velocity discontinuities only at points O + and L − .

V. DISCUSSIONS AND CONCLUSION
We studied the variational boundary value problem of the electromagnetic two-body problem with shortestlength boundary-segments in a neighbourhood of circular orbits with large inter-particle separations, r ij ≫ 1. The Wheeler-Feynman equations reduce to a two-point boundary problem, (35) and (36). For this case the initial value problem given by (35) is well-posed and the solution is unique (Theorem 1). We observed that the shooting method converged even for some perturbations of circular orbits with radii r ij ≈ 0.5 and having relativistic velocities v i ≈ 0.9, provided that n ij± ⋅ v j± ≈ 0.3. The shooting problem uses up all the initial-velocity freedoms and the occurrence of discontinuous velocities at points O + and L − is expected even for C 2 perturbations of circular-orbit segments (i.e., without breaking points in histories). We have also shown existence of solutions with discontinuous velocities for nearcircular boundary-segments having discontinuous velocities in histories. For boundary-segments with continuous velocities, trajectories may still have discontinuous velocities satisfying Eq. (43) for inter-particle separation in the nuclear magnitude r 12 ≃ 1 √ m 1 m 2 , a case described by the algebraic-differential equation (29). Situations where velocity discontinuities necessarily occur are (a) one of the boundary-segments has discon-