{'fraction_non_alphanumeric': 0.03053704178415791, 'fraction_numerical': 0.009644962830043579, 'mean_word_length': 4.57502679528403, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Radiative Damping and Functional Differential Equations

We propose a general technique to solve the classical many-body problem with radiative damping. We modify the short-distance structure of Maxwell electrodynamics. This allows us to avoid runaway solutions as if we had a covariant model of extended particles. The resulting equations of motion are functional differential equations (FDEs) rather than ordinary differential equations. Using recently developed numerical techniques for stiff FDEs, we solve these equations for the one-body central force problem with radiative damping with a view to benchmark our new approach. Our results indicate that locally the magnitude of radiation damping may be well approximated by the standard third-order expression but the global properties of our solutions are dramatically different. We comment on the two body problem and applications to quantum field theory and quantum mechanics.


A. Motivation
It is remarkable how little is known about classical solutions of relativistic many-body electrodynamics with radiative damping. Textbooks often describe simple non-relativistic solutions to the one-body problem with an external force and radiative damping. Studies in accelerators and antennas also consider radiation emitted from charged particles. However, when many bodies are included and the back-reaction of the radiation on the particles involved is considered, the resulting equations are amenable neither to analytic solutions nor to numerics. We explain the reason for this in greater detail below but we emphasize here that the difficulties involved are not merely technical but relate to important questions about the nature of classical electrodynamics.
In fact, until recently, there has been no systematic development of techniques to integrate the equations of relativistic many-body electrodynamics with radiative damping. One reason for this is that, hitherto, there has been no compelling physical reason to consider the intricacies of this problem. For example, in astrophysical processes, the radiation emitted by charged bodies is of interest; however the back-reaction of this radiation on the bodies is small compared to other forces. Hence, the radiation reaction is often ignored or included only heuristically.
In recent years, this situation has changed. Simulations of gamma ray bursts, for example, have reached a level of sophistication that authors have started to investigate the influence of radiation reaction [1]. A more immediate motivation for this study comes from papers of [2,3,4]. In these papers, it was conjectured that some classical orbits of the many-body system have special properties when radiation reaction is considered. However, further progress on this front was prevented by the absence of techniques to integrate the manybody system 1 with radiation damping.
In this paper, we develop such a technique. Our technique may be applied to check the conjectures of [2,3,4]. Moreover, this technique is of independent interest because it is applicable to the investigation of radiative damping effects in astrophysical processes. It would also be interesting to include the effect of radiative damping in molecular-dynamic simulations carried out in codes such as charmm [5,6], amber [7] or gromos [8].
In this paper, we elaborate our technique and explain its advantages. Furthermore, we bench-mark our technique by applying it to the one-body problem. This is feasible since, as we explain further in section I C, some solutions to the one-body Lorentz-Dirac equation have already been obtained. Our technique agrees, excellently, with these solutions. Having established the validity of our methods, we hope to apply them to the concrete physical problems described above in a subsequent paper.

B. Setting
Consider a point particle with mass m, charge q and world-line α µ (τ ) parameterized by its proper time τ . In the presence of an external electromagnetic field strength F ext µν the equation of motion of this particle is given by: The second term on the right may be understood as the regularized action of the particle on itself. We emphasize that the structure of this term is quite delicate and it cannot be modified without modifying the Maxwell theory itself.
This is a third order differential equation and requires 3 initial conditions. However, for a generic choice of initial conditions α(0),α(0),α(0), we obtain a runaway solution where the velocity of the particle increases continuously in the future! This solution is clearly unphysical and the existence of such runaways is an embarrassment for the classical theory of electromagnetism. Several solutions to this problem have been proposed in the literature.
We discuss some of these, in Section I D and Appendix A. However, one common solution to this problem is as follows. Let us say, we assign the initial acceleration to satisfÿ with κ = q 2 6πǫ 0 mc 3 , then the asymptotic solution is free of runaways [9]. A commonly articulated philosophical problem with this is that it requires one to know all forces that act in the future on the particle. Moreover, the particle starts accelerating before any force is actually applied. This is the problem of pre-acceleration.
However, we are more concerned with a different problem that makes Formula (2) completely unsuited to any actual computations. The point is this: if we chooseα(0), to be even infinitesimally different from the value given by (2), the solution ceases to be bounded and turns into a violent runaway. Stated another way, equation (2) specifies a co-dimension 1 surface in the space of all possible initial conditions. The behavior of solutions to equation (1) is discontinuous in this space. On this surface, the solutions are bounded; away from this surface, the solutions comprise runaways.
From a computational point of view, this makes any numerical computations with (2) impossible since we need to specify initial conditions exactly at the critical value prescribed by (2). From a philosophical point of view, this need to 'fine tune' the initial conditions to obtain a meaningful solution is rather unappealing.
Hence, in this paper we do not adopt this point of view but look for another method to tame the runaway solutions of classical electrodynamics.

C. One-Body vs Many-Body
We now pause to explain a critical difference between the many-body system and the one-body system. For the one-body system there exists a simple method that often allows us to sidestep the above complications.
This involves solving equation (1) backward in time. Here, the runaways solutions get rapidly damped, leaving behind the 'physical' solution. To our knowledge this technique was first used by Plass to obtain the solution for a charged particle orbiting an infinitely massive object with the opposite charge [10].
However, for the many-body problem this method fails. This is because the field acting on a particle at a given time depends on the positions of the other particles in the past.
Stated another way, the equations of motion for a relativistic many-body system are functional differential equations. Such equations can be solved forward in time, given past data.
However, in general, they do not have a unique solution backward in time. 2 Hence the trick of solving the equations of motion backward in time to get rid of the runaways cannot be applied to the many-body problem.
In this paper, we will develop a technique that is applicable to the many-body system.
As we have already mentioned, we will check this technique against extant solutions of the one-body problem obtained using the trick above. Thus this trick serves to bench-mark our method; while the trick itself does not generalize to the many-body case, our method does.

D. Other approaches to this problem
The formidable conceptual and calculational problems outlined above are sometimes wished away by appealing to quantum mechanics. This argument goes as follows. Equation (1) is only an approximate equation valid at distances large compared to the Compton wavelength of the particle. At shorter distances we should use the full quantum theory.
Hence, according to some authors, we need not take the difficulties above seriously, since the equation that they arise from is only approximate in the first place.
We do not find this point of view satisfactory. For one, we believe that the classical theory should be well defined in its own right. Second, in the concrete physical situations where we wish to use (1), it is hardly feasible to use the entire quantum theory.
However, the observation that equation (1) is valid only at distances large compared to the Compton wavelength suggests another possible resolution to the problem. This is to modify (1) itself at distances that are short compared to the Compton wavelength of the particle.
This approach has been adopted by many authors. One popular idea has focused on giving the particle itself a structure. For an excellent review of models of extended particles we refer the reader to [12]. However, it is hard to construct a Lorentz invariant theory of extended particles. When this can be done (see for example [13]), one needs to worry about the forces that keep the extended particle stable. In fact, additional terms coming from the energy due to the internal stresses that hold the particle together, must be added to the mass of the particle to obtain covariant equations of motion. This is the famous '4/3' problem (see [13] for a review).
A second approach is to place a momentum cutoff in the theory. Since a momentum cutoff can be only be placed after Wick rotating the momenta to Euclidean momenta, it respects neither Lorentz invariance nor gauge invariance.
We discuss a third approach to this problem in Appendix A.

E. Our Technique
To summarize the discussion above, we are looking for a modification of equation (1).
We would like this equation to be a natural consequence of a modified classical theory.
Furthermore, we would like such a theory to satisfy four properties: 1. It should be Lorentz covariant and gauge invariant.
2. Given generic past data, the theory should possess regular solutions in the future.
3. It should reduce to the usual Maxwell theory at large distances 4. It should be simple enough to make solutions of the many-body problem with radiation damping feasible.
We have seen that the method of assigning the initial acceleration via equation (2) fails on count 2 above while other methods of assigning a structure to the particle fail on counts 1 or 4. Yet other methods discussed in Appendix A are based on ad hoc approximations that are unacceptable because they modify the second term in (1) without modifying the Maxwell theory. In this paper, we will propose a modification of Maxwell electrodynamics that satisfies all four conditions above.
We modify the propagators of classical electrodynamics at small scales in position space, while retaining Lorentz covariance and gauge invariance. The effect is as if we had a covariant model of extended particles. Since, however, it is electrodynamics which is being modified, this allows us to sidestep issues of internal stability of the particle and the equations of motion are covariant without the addition of ad hoc terms. Further, the simplicity of this model makes numerical calculations feasible. At long distances our theory reduces to the usual Maxwell theory.
In the analogy with extended particles, the length-scale at which electrodynamics is modified corresponds intuitively to the 'size' of the extended particle. This finite size introduces a delay in the interaction of the particle with itself and so we obtain a functional differential equation even for the one-body problem with radiation damping. Although the resulting FDEs are numerically stiff, we are able to solve them using numerical techniques recently developed by Guglielmi and Hairer [14,15].
We find the following surprising result. Locally, the force of radiation damping is well approximated by the customary third order expression. However, globally our FDE has dramatically different properties from the solutions of equation (1). In the appropriate regime of parameters, our theory does not exhibit any runaway behavior. The first part of our result ensures that standard calculations of radiation emitted by antennae are not affected noticeably.
As already mentioned, the simplicity of this theory makes numerical solutions of the many-body problem feasible. In this paper, we solve the important test case of a radiating particle in a central potential. The particle gradually spirals in, as is expected on physical grounds. Furthermore, we are able to read off the rate of the decay of this orbit. As already stated for this one-body problem, we can obtain this solution, even in Maxwell electrodynamics, by means of the trick explained in Section I C. The two solutions are in good agreement.
In a subsequent paper [16], we hope to apply the techniques developed here to the full relativistic two-body problem.

II. THE MODEL
Given a 4-current density, j µ (x), the potentials, A µ are given by the formula where the (retarded) Green function is given by [17] 3 χ is an arbitrary scalar function which we will henceforth ignore by choosing Lorenz gauge We emphasize that this is purely for convenience; our results below are independent of this gauge choice. In taking the norm, we use the metric diag(−c 2 , 1, 1, 1) i.e.
Along the lines of [19], in this paper we will modify the short distance structure of Maxwell electrodynamics by modifying G to In the expression above, the null cone has been replaced by a hyperboloid; this regulates the singularities of the theory and leads to a theory analogous to a theory of extended particles with characteristic size d. It is curious that for i (x i − y i ) 2 = 0, equation (6) forces x 0 − y 0 = d. For particles extended in space, one would expect the Green function to have support on a spatial interval at zero time. The behavior above is somewhat different and we may think of our particles as being extended in time! However, we caution the reader against taking this analogy too seriously. The mathematical content of the theory is specified by equation (6) and this is all we will make reference to below.
The delta function in (6) is written in terms of manifestly Lorentz invariant quantities.
Second, for a vector satisfying (x − y) 2 = −d 2 , a Lorentz transformation cannot change the sign of x 0 − y 0 . Hence, the Green function in (6) is Lorentz invariant. Consequently, our modified theory preserves Lorentz invariance.
Consider now a point particle of charge q, with worldline specified by a function α µ (τ ) of the proper time, τ , of the particle. We define the field strength, as usual. 4 Using equations (3), (6), we find that the field produced by this particle at a point x is given by whereα µ ≡ dα µ dτ andα µ ≡ d 2 α µ dτ 2 are evaluated at the retarded proper time τ 0 , which is defined as the solution to that has x 0 > α 0 (τ 0 ). Further, ζ µ = x µ − α µ (τ 0 ). For a slow moving particle, the delay (9) where F µν is the net field strength and includes not only the contribution of other particles but also a contribution from the past history of the particle itself as we explained above.
Despite this, the expression for F µν in (8) is free of divergences. Moreover, the self-interaction in (10) leads to a finite radiative damping term.
Second, notice that even for a free particle the equation of motion (10) is in general a functional differential equation because it relatesα µ (τ ) on the left to α µ (τ 0 ) and its derivatives on the right. To solve this FDE we need the past history of the particle, apart from initial data. For FDEs where the delay is bounded above by some δ 0 , δ < δ 0 , the past history of the particle specified on an interval of length δ 0 (with initial data specified by continuity at the endpoint of this interval) suffices to fix a solution A. An Example Qualitatively, one can understand by a simple linear example how functional differential equations help in getting rid of the runaway solutions of the Lorentz-Dirac equation.
Consider the equation where a, ǫ > 0 are constants. If we expand this equation to first order in ǫ, we find If we discard the higher order terms, then (12) has the solution with κ 1 , κ 2 being constants of integration determined by the initial conditions for f and f ′ .
Note, that for arbitrary κ 2 = 0, these solutions grow exponentially in the future.
On the other hand an exact solution to (11) requires past data over the entire interval [0, ǫ].
Substituting, e mt as a test solution in (11), we find that m must satisfy the characteristic equation It is easy to show that all solutions to (14) have Re(m) < 0. Hence, (11) actually has no runaway solutions of the form (13). The process of expanding FDEs to obtain ODEs, by means of a 'Taylor' series, may lead, as above, to spurious solutions; but this is the process that is used to derive the radiation damping term. If ǫ is very small, then locally it is indeed true that f ′ (t) is well approximated by (12); nevertheless the global properties of the solution (13) are completely incorrect! This leads us to expect that runaway solutions may be avoided by replacing ODEs of the Lorentz-Dirac theory with FDEs, an expectation confirmed by the calculation carried out below.

B. Many Particles
To conclude this section, we would like to remark that the formalism discussed here works for many particles as well as it does for a single particle. If we have many particles, we calculate the field due to each particle by using (8). The equation of motion for each particle is now given by (10), where F is the total field strength at a point, obtained by summing all the individual field strengths.
It is notable, that equation (10) is considerably simpler than (1) which applies in Maxwell electrodynamics. To obtain (1), we need to separate the field strength into two pieces: (a) the field strength due to all other particles and (b) the field strength due to the particle itself. The second term -when regularized -leads to the second term in (1). It is this term that is responsible for the runaway solutions that we have discussed. Equation (10) involves no such artificial separation. The self action of the particle on itself is regular.
Since it is this term that leads to runaways in Maxwell electrodynamics and since we do not find runaways in our one-particle example below, we will not find runaway solutions in a many-body system either.

III. NUMERICAL SOLUTION FOR A CENTRAL FORCE
We now turn to a description of how the equations of motion (10) may be solved numerically for the important test case of a charged particle moving in an attractive Coulomb field.
We represent the effect of this field by the field strength where r = x i x i . The equation of motion for a particle of charge q moving in this field is now given by (10) with where F self is calculated, as explained in the previous section, with the help of formula (8).

A. Choice of Code
To solve the equations (10), we need a numerical solver for stiff functional differential equations. The equations above are stiff for two reasons. First, the time period over which we want to solve the equations is very long compared to the delay. Second, the magnitudes of the Coulomb force and the radiation reaction differ substantially. Two popular earlier codes for FDEs, retard [20] and archi [21] both use the Runge-Kutta coefficients of Dormand and Prince. While these coefficients can be used to construct excellent general-purpose solvers for ODEs, they are not appropriate for stiff equations. (For an example, of how the dopri code fails, see [22].) In contrast, the radar code of Guglielmi and Hairer [14,15]

adapts to FDEs an implicit
Runge-Kutta (collocation) method based on Radau nodes. The corresponding radau code [23] is a long-established code which has successfully been applied to various stiff ODEs.
Another technical issue is that, although we have assumed retarded propagators, the FDE (10) is classified as a 'neutral' equation, in the classification of differential equations with deviating arguments [24], rather than a retarded equation, since the second derivative of α µ (τ 0 ) appears also on the right. Both archi and radar codes may be used also for neutral equations, although the radar code is significantly more efficient(by a factor of about 10 3 in our experience) in view of its better stiffness properties. The efficient solution of such FDEs requires code that can take step sizes much larger than the delay: the radar code allows for this possibility. The successful use of the radar code, however, requires relatively greater preliminary work to evaluate various Jacobians.

B. Rewriting the Equations for Numerical Solution
For the numerical solution it is useful to cast the equations (10) in 3 + 1 form and also rewrite them as first order equations. It is easy to see that with the introduction of an intermediate three vector u i , and an external time t, we may rewrite the equations as where γ is determined by and F is calculated using equations (8), (15), (16).

C. Units, Initial Condition and Parameters
The next task is to choose appropriate units. Since we wish to deal with values that are relevant to the hydrogen atom, it is convenient to choose units so that where For the past history, we used a circular orbit with initial velocity appropriate for the first Bohr orbit. This is given by v in = 0.00730c.
However, as we explain below, the solution is not very sensitive to the exact value of d. For comparison, the classical radius of the electron is d 0 = 1 4πǫ 0 q 2 2mec 2 = 0.417 × 10 −4 in our units. The mass that we actually insert into equation (17), is taken to be This is because, as explained below, the self-action term, locally contains a part that may be thought of as 'renormalizing' the mass back to the value m e .

First, let us briefly investigate the solutions to the Lorentz-Dirac equation in a central
force This ODE is already stiff because of the large difference in the Coulomb term, the Schott term and the radiative damping term. However, it may be solved using the code radau [23].
The solution is shown in Figure (2). It may be seen that the particle shoots off with the a velocity that is very close to the speed of light. This is clearly unphysical. We now turn to the solutions of the equations (17). Recall that to solve the FDE we need to specify the motion of the particle over an appropriate past interval. As stated above, we take this past history to correspond to a circular orbit with the parameters of the first Bohr orbital. The solution obtained by solving (17) using the code radar is then shown in Figure 3. The damping due to radiation causes the curve above to deviate very slightly x from the circular orbit which would have been the solution without radiation. This deviation is plotted separately in Figure 4. Note, that while the particle does spiral in towards the center, this is through a complicated elliptic spiral. The average rate of decrease of distance towards the center may be measured from the graph in Figure 4 and is found to be δr δt ≈ −1.13 × 10 −6 .
This average value is in good agreement with the naive rate of electron decay, δr δt ∼ −4c 3r 2 q 2 4πǫ 0 mec 2 2 , which is derived by assuming that the orbit remains circular while losing energy at a rate given by the Abraham Lorentz formula in each revolution Second, we note that the solution presented above is not very sensitive to the parameter  (20), we find that the two solutions are almost identical and the deviation between them, in the radial coordinate, oscillates with an amplitude of 4 × 10 −11 and does not increase with time. This may be seen in Figure 5 We also note that this solution is in excellent agreement with the backward solution of the Lorentz-Dirac equation discussed above.
Locally, the instantaneous radiation force is described very well by the expression: The term c 1 effectively renormalizes the mass. This is the reason, we choose m as in equation (20). f rad i oscillates as the particle moves in the orbit, and in the Table below  We see that the fit is always very good and although, initially, the fit becomes slightly better as we move to smaller d 2 this gain may be outweighed by numerical errors beyond a certain point.
The stiffness of the equations, and hence the time taken for numerical solution, increases as d decreases. Eventually, the numerical solution becomes unstable for a limiting value of d, that we find empirically to be d ∼ 1.66d 0 , where d 0 is the classical radius of the electron.
V. CONCLUSIONS 1. In this paper we discussed the difficulties in obtaining solutions to the Lorentz-Dirac equations. We described in Section I C, why these difficulties are especially acute for the many-body problem.
2. Next, we proposed a modification of Maxwell electrodynamics via the modification of the Green functions in (6). The modified Green function regularizes the self-action of a particle. This leads to the equation of motion (10). A radiation damping termthat arises from the self-action -is automatically included in (10) since F µν involves a contribution from the particle itself.

5.
As an example, we solved our new equations of motion for the important test case of a single charged particle moving in the Colulomb field defined by (15). Figure 2 shows how an attempt to numerically integrate the Lorentz-Dirac equation in this physical situation leads to violent runaways. Figure 3 shows the solution of our modified theory which is manifestly regular. Figure 4 shows that if one looks at this solution closely, one finds that the particle is slowly spiralling in, due to radiation losses, as one expects on physical grounds. For this simple problem, solutions to the Lorentz-Dirac equation may also be obtained using a trick explained in Section I C. The two solutions are in excellent agreement.
6. Having established and tested our technique in this paper, we hope to apply it to other concrete physical problems in a forthcoming paper.

VI. FUTURE DIRECTIONS
We expect that the technique we have described here will have wide applicability. An immediate application of our technique is to the simulation of relativistic plasmas in astrophysical contexts [25,26].
Another application is to the two body problem that has seen some recent interest. The possibility that the classical solutions corresponding to quantum Bohr orbits may have distinctive features was discussed in [2] and also conjectured from a very different viewpoint in [3,4]. These possibilities can now be investigated numerically and may lead to a better understanding of the link between differential equations with mixed type deviating arguments and quantum mechanics discussed in [27]. (Recently, a link between pure delay differential equations and quantum mechanics was also discussed, from a somewhat different point of view, in [28]).
Long ago, the problem of radiative damping was taken up by Wheeler and Feynman in the hope that clarity about classical electrodynamics might help to resolve the divergences of quantum electrodynamics. In this context, it is noteworthy that the technique of modifying the support of the propagators in position space immediately generalizes to a method of regularization in quantum field theory, in which context it was originally proposed [19]. The regulator here is not merely formal, as in dimensional regularization, and moreover, unlike naive cutoffs it preserves Lorentz invariance. This may present some advantages when one is interested in studying the effect of a UV cutoff. In the quantum theory, the unitarity of the S-matrix places stringent restrictions on the propagator. This leads to additional difficulties that we do not discuss here.
Finally, in the above discussion we have taken for granted that the interactions involved must be retarded. In [27] it was shown that the more general equations with mixed-type deviating arguments could be used to recover much of the mathematical formalism of quantum mechanics. In the present context, the use of mixed-type equations may be expected to decrease the radiative damping from the full value it has in the retarded case. However, solving mixed-type equations involves additional mathematical complexities. We defer this discussion to a future paper.
A covariant version of this equation was derived in [29]. Although the paper initially claimed that the resulting equation was exact, this was shown not to be the case in [31]. Furthermore, since the structure of the Lorentz-Dirac equation is tightly dictated by Maxwell electrodynamics a modification of the form (A2) cannot be performed without modifying the Maxwell theory itself. Unfortunately, no such modification was proposed in [29]. Finally, when relativistic effects are included we note that, once again, triple derivatives appear on the right hand side of (A2) which leads once again to third order equations of motion.

APPENDIX B: NEW FIELD EQUATIONS
Here we briefly discuss the new field equations that result from the modification of propagators in equation (6). We stress however that this is peripheral to the program that we have presented. Given any distribution of particles, we can derive the corresponding field strengths using formula (8). We can then calculate the subsequent motion of all particles using (10). In this formulation of the many-body system in terms of functional differential equations, the fields are only an intermediate crutch.
However, for some formal purposes it may be of interest to analyze the field equations themselves. The modified field equations are easily obtained by means of a Fourier transform of G. Since these field equations are non-local but yet linear they are more conveniently represented in Fourier space and this is what we shall present here.
To lighten the notation below, we shall put c = 1, in this Appendix.
Our task is to find the Fourier transform of the modified Green function (6).
Formally, this integral is identical to the Fourier transform -from momentum space to position space -of a massive propagator. Using the standard techniques, developed for this (see, for example [32]), we find −iG(p) (2π) 3 = 1 4π s(p 0 )δ(p 2 ) + idP −θ(−p 2 ) where J and Y are Bessel Functions of the first and second kind, K is a modified Bessel function [33] and s(p 0 ) is the sign of p 0 . All terms above should be understood to be distributions and P denotes the principal value. Note that the expression (B2) is Lorentz invariant. This is because, for p 2 ≤ 0, we cannot change s(p 0 ) by means of a Lorentz transformation; all other terms in (B2) are written in terms of explicitly Lorentz invariant quantities.
Let us consider this expression near p 2 = 0. In that limit, (B2) may be written as: In the physics literature, the first two terms above are often combined This also tells us that, in the limit d → 0, where A(p) and j(p) are the Fourier transforms of the vector potential and the current density.