Radiation from charges in the continuum limit

It is known that an accelerating charge radiates according to Larmor formula. On the other hand, any DC current following a curvilinear path, e.g. a circular loop, consists of accelerating charges, but in such case the radiated power is 0. The scope of this paper is to analyze and quantify how the radiation vanishes when one goes to the continuum DC limit.


I. INTRODUCTION
There are many physical configurations in which discreet charges are in acceleration, but in spite of that they behave almost like steady state DC. This situation is encountered in any DC or low frequency electrical circuit, because in spite of the fact that the charge distribution is always discreet and each charge accelerates in the influence of the electric field, one may consider the charge distribution as almost continuous. For example the model used for conducting materials is of an average drift velocity for positive and negative charges which is proportional to the electric field v ± = ±µ ± E, µ ± being the mobility of the positive and negative charges, respectively. Hence the current density is J = ρ + v + +ρ − v − = (ρ + µ + − ρ − µ − )E, where ρ ± is the charge density of the positive and negative charge carriers, ρ − being negative, so that the total charge density is ρ = ρ + + ρ − = ǫ 0 ∇ · E. In case the charges are protons and electrons, in a conductor, the protons mobility is µ + = 0 and ρ + = −ρ − , so that usually ∇ · E = 0 like in free space. In addition, the model considers ρ − to be almost uniform, so that one defines the conductivity σ = −ρ − µ − , resulting in Ohm's law J = σE.
Another situation, for which charges are somehow "more discreet", are DC or low frequency ion drift devices, and for those we usually have one type of charge carriers, say positive ions. Here we use J = ρv = ρµE, and one may not assume the charge density is uniform, but rather has to use Gauss's law ǫ 0 ∇ · E = ρ, resulting in a set of nonlinear equations. In principle, such problem is time dependent, because discreet ions move in the space, but it appears that the DC approximation ∇ · J = 0 works very well for such cases [11][12][13][14].
In either of the above situations, currents may follow curvilinear paths, in which case, the charges clearly accelerate, also if the magnitude of their velocity is constant, and still the DC model works well for most practical cases for which the density of the charge carriers is big.
The purpose of this work is to understand how the radiation vanishes when the density of the charge carriers approaches the continuum steady state. Certainly, the radiation has to disappear gradually, so that this vanishing may be quantified. Some preliminary work [10] has been done in this direction, but because this work has been done a priori in a non relativistic approach, its results are inaccurate and incomplete.
The calculation is done in a canonical configuration of charges in circular motion at constant speed. The configuration and the formulation are explained in section 2. In section 3 we calculate the fields, explain their behavior and derive an exact expression for the radiated power. For completeness, we also calculate the radiation reaction and show that the power needed to support the radiation equals the radiated power. In section 4 we derive an asymptotic result for the case the number of charges is big (continuum limit). As a special case, we also derive the limit for slow charges, and this case is of importance because it represents the typical situation of DC currents in devices as discussed before. The paper is ended with some concluding remarks.

II. CONFIGURATION AND FORMULATION
A total amount of charge q is rotating in a circle of radius d at constant speed v so that the angular velocity is ω = v/d. The charge q is "split" into N charges of value q/N, uniformly distributed around the circle, so that the charge number k is at the angle ωt + 2πk/N, where The configuration is shown in Figure 1, for N = 3.
The location of the charge k as function of time is given by: The fields propagate with the speed of light c. Hence the fields at the observer location r at time t are influenced by the motion of each charge, at an earlier (retarded) time.
Specifically, the fields are influenced by the motion of the charge k at time t ′ k so that At large distance from the charges, one may approximate: where hence the retarded time t ′ k may be calculated from the following implicit equation: which may be solved numerically by setting a "1st guess" t ′ k = t − r/c in the right side of eq. (5) and recalculate t ′ k until convergence is obtained. Figure 2 emphasizes the meaning of retarded positions of the charges.
Our purpose is to calculate the power radiated by those rotating charges. Of course, for N = 1 the result is known by Larmor formula (and will be confirmed later on).
The general form of the Larmor formula [1][2][3][4][5][6][7][8][9] is: FIG. 2: (color online) The dark colored spheres represent the charges at the current position at time t and the light colored spheres represent the charges at the retarded positions at times t ′ k . The retarded positions are connected with dashed lines to the observer location, and those distances are called R 0 , R 1 and R 2 . This emphasizes that the field at observer is determined by the velocities and accelerations of the charges at the retarded times.
where β ≡ v/c. and If one defines the angle between the velocity and the acceleration as α, one may express (β ×β) 2 = (ββ sin α) 2 an rewrite In our case of circular motion, the velocity is perpendicular to the acceleration, so that sin 2 α = 1. Therefore by replacing 1 − β 2 = 1/γ 2 , the Larmor formula simplifies for our case to where a is the acceleration. We will calculate how this power decreases when the number of charges N increases.

III. FIELDS AND POWER CALCULATION
To calculate the power radiated from the collection of charges in Figure 1, one needs only the far fields, i.e. those who behave like 1/R. The far electric and magnetic fields due to the moving charge k are given by [2,3]: and where R k is the unit vector pointing from the position of the charge to observer and R k is the distance between the charge and the observer, as defined in eq. (2) -see Figure 2.
is the velocity relative to c and a k =v k is the acceleration of the charge. All the dynamical variables are evaluated at the retarded time (defined in eq. (5)).
Defining r as the unit vector pointing from the coordinates origin to the observer, we may calculate in the far field the difference: Hence one may use r instead of R k in eqs. (11) and (12) with an error of order 1/R 2 , which does not affect the calculations of the radiated power. Also, in the denominator of eq. (11) we may set R k = r, as always done for far field. So we express the electric and magnetic fields as the sum of the contribution from all the charges: and Now using eq. (1), we evaluate r × [( r − β k ) × a k ] and express it in spherical coordinates: and β k · r evaluates to Setting those results into eq. (14), we obtain and where the functions Fc and Fs are defined as: and and the functions fc and fs are defined as and and p is defined as and is a parameter which controls the behavior of φ k , as will be soon shown.
The power per unit of normal area (or Poynting vector) is given by, The total power is calculated via which results in Here we factored out the Larmor formula for the radiation of a single charge -see eq. (10), so that G(t, β, N) is dimensionless and represents the decay of the power. The function G(t, β, N) is given by hence G = 1 for N = 1, for any t or β and the function F is Now we eliminate t ′ k from eq. (5) and rewrite the implicit eqs. (5) and (4) in terms of φ k This allows us to change variable ϕ ′ = ϕ − ω(t − r/c) in (28) obtaining: so that eq. (30) is rewritten as We see that if φ k and ϕ ′ satisfy eq. (32), also φ k − 2π and ϕ ′ + 2π satisfy it, hence cos φ k and sin φ k are periodic functions of ϕ ′ , with a periodicity of 2π. Therefore, the dϕ ′ integral in eq. (31) may be evaluated over any period of 2π, showing that G (and therefore also the radiated power P ) does not depend on time, so that we may simplify eq. (31) to We redefined G to be time independent, so that fc, fs, Fc and Fs in eqs.  Fc elements Fs elements   Fc elements Fs elements  It is interesting to remark that the average of the functions fc and fs is always 0, although this is not always visible for the fs functions. This may be proved by calculating where for brevity we called fc,s the functions fc or fs, and we called their average fc,s .
We also use the abbreviation hc,s for the functions hc and hs defined as Fc elements Fs elements for the fc and fs average calculation, respectively. We change variable from ϕ ′ to φ k , and we find from eq. (32) that getting: where φ k (0) is the value of φ k at ϕ ′ = 0. Because the integrand has a periodicity of 2π, one may integrate over any period of 2π, obtaining: This integral may be solved by the residue method on the complex plane. After changing variable z = exp(iφ k ), one obtains fc,s = 1 2π 2i where C is the counterclockwise unit circle integration contour shown in Figure (9) and lc,s are abbreviations for and where b is a pure imaginary number defined by and the 2nd order poles are (expressed in terms of p or b) where indices 1 and 2 refer to upper and lower signs respectively -see Figure (9).
In the last expression of the poles in terms of b, the magnitude under the square root is real and negative, so the square root is understood to be positive pure imaginary.
Also we named the integrand in eq. (40) Lc,s, defined as to ease manipulations. We remark that only z 1 is inside the integration contour (see Figure (9)), and its residue is Using the relations z 1 z 2 = −1 and z 1 + z 2 = −2i/p, one finds that for both fc and fs cases the nominator of eq. (46) is 0, and hence The poles z 1,2 are negative pure imaginary, and |z 1 | < 1 and |z 2 | > 1, the integration contour C is on the unit circle and the integration contour C 1 is on a circle of radius smaller than |z 1 |.
Now we continue with the calculation of G in eq. (33). Rewriting eq. (32) for the charge which may be rewritten as showing that if we could explicitly express φ k (ϕ ′ ) from eq. (32), φ m would be the same function of ϕ ′ only shifted: as also evident from Figures 3 -8. Therefore both Fs and Fc functions remain unchanged for a shift of multiples of 2π/N, say: and the sum being on all charges (and k is modulo N), we are left with the same result.
We may therefore rewrite in eq. (33) and after changing variable ϕ ′′ = ϕ ′ − 2πn/N, we are left with N identical integrals over the period 0 to 2π/N. For simplicity we rename ϕ ′′ back to ϕ ′ and rewrite the function G as: and this will significantly reduce the time of a numerical integration. Now looking at the θ dependence of Fs and Fc, we remark from eq. (32) that φ k depends on p = β sin θ, hence it is invariant under replacing θ by π − θ, and so is cos 2 θ. We may therefore replace the integration from 0 to π by twice the integration from 0 to π/2, getting Now we change to the variable p defined in eq. (24) and rewrite eq. (54) obtaining: For the case of N = 1 we know G must be 1, but for consistency we shall prove it: where Fc 2 and Fs 2 reduced here to a single term. By changing integration order, we may perform the dϕ ′ integration by the change of variable dϕ ′ /dφ 0 defined in eq. (37), obtaining All the functions have a 2π periodicity on φ 0 , so one may use any limits of interval 2π for φ 0 . By changing variable z = exp(iφ 0 ) and using the residue method on the complex plane we obtain hence which results in We perform now the calculation in eq. (55) numerically. Knowing that G(β, 1) = 1, this calculation actually shows the power radiated by N charges, divided by the power radiated For big N , G goes asymptotically to 0 and as smaller β is, G goes faster to 0.

A. The radiation reaction
We calculate now the Lorentz force on the charges and the resulting radiation resistance power for comparing with the radiated power. We will need the electric field on charge n at time t due to charge m at its retarded position at the earlier time t ′ , so that: Using the expression for r ′ k in eq. (1) we obtain which is exact for any m and n. By definition, t − t ′ > 0, so to take the correct square root from the left side of this equation, we need to know the connection between m and n.
To simplify, we restrict: for which we obtain where Φ nm is defined by Now we isolate t − t ′ from eq. (66) and set it in eq. (65), obtaining which is an implicit equation, that can be solved by setting a 1st guess Φ nm = π(n−m)/N in the right side of the equation and recalculate Φ nm till convergence is obtained.  The retarded distance c(t − t ′ ) is the big (red) segment on the 2Φ mn arc, hence equal to 2d sin Φ mn according to eq. (65). We also see that the big angle 2Φ mn (marked in red) equals the sum of the angles ω(t − t ′ ) and 2π(n − m)/N according to eq. (66). Two orthogonal unit vectors (green) r(t) and ϕ(t) are drawn near charge n, representing the radial and tangential directions of the moving charge.
It is to be mentioned that the solution Φ mn of eq. Now the electric field on charge n due to charge m at its retarded position is given by where the 2nd part is the far field which we used in eq. (11) (after replacing a byβc and 1/ √ ǫ 0 µ 0 by c) and the first part is the near field which behaves like 1/R 2 .
The quantities appearing in eq. (68) are: β m is the retarded velocity of charge m (relative to c), R mn is the distance between the retarded position of charge m and the current position of charge n (and equals to c(t − t ′ ) -see Figure (11)), and R mn is the unit vector pointing from the retarded position of charge m to the current position of charge n.
We need the field at charge n in local components r(t) and ϕ(t) (see Figure (11)) and we calculate now all the needed quantities to express the electric field E mn . So we obtain and which is easily derived also from Figure (11), because R mn = c(t − t ′ ), which equals to 2d sin Φ nm according to eq. (65). Putting eqs. (69)-(72) in eq. (68) we obtain With the aid of this field we will calculate the total force acted on a charge by the other charges, but we also need the "self" radiation reaction force. For the most general case, this is given by [2,3] In our case, of circular motion, the velocity is perpendicular to the acceleration so β ·β = 0, we therefore remain with 2 terms For the circular motion we know thatβ = −ω 2 β, and using ω = v/d = βc/d, eq. (75) reduces to showing that the self reaction force is in the direction opposite to the velocity of the charge.
Now we chose n to be the "last" charge, i.e. n = N − 1, hence the restriction in eq. (64) holds, and calculate the total force on it, given by the self force plus the force acted by all other charges (i.e. the Lorentz force): where v N −1 (t) is the velocity of charge N − 1 and B m,N −1 is the retarded magnetic field on charge N − 1 due to charge m.
The dumping power on charge N − 1 is given by v N −1 (t) · F N −1 , therefore we do not need the magnetic part, obtaining The velocity of the given charge N −1 is in the tangential direction, i.e. v N −1 (t) = cβ ϕ(t), so only the tangential part of eq. (73) affects the power. We obtain Because the dumping power on one charge does not depend on time and by symmetry is the same for all charges, the total dumping power on the whole system of charges is just N times the above: The dumping power must be identical with the radiated power, with a minus sign, so to compare them, we may factor out − P | N =1 from eq. (10): where which clearly shows that for N = 1, the sum is 0, so that G dump (β, 1) = 1 for any β.
The numerical calculation of eq. (82) shows identical results with those of G calculated from eq. (55), as shown in Figure 10.

IV. ASYMPTOTIC RESULT FOR MANY CHARGES
To understand how the radiation goes to 0 when the number of charges N goes to infinity, one may try to approximate either G from eq. (55) or G dump from eq. (82), for big N.
Although G dump seems more compact, it is more difficult to handle, and we shall develop G for large N.
As mentioned in the previous section (see , for large N the functions Fc and Fs tend to be harmonic, hence we may approximate them by the first term of their Fourier series. For a function a(x) with periodicity X, we specify the Fourier coefficients by and a(x) is represented by its Fourier series For brevity, to refer to the functions Fc and Fs we call them Fc,s (as in eq. (34)). We know those functions have a periodicity of 2π/N in ϕ ′ (see eq. (51)), so we define their Fourier coefficients: and the functions Fc and Fs are expressed as: Now we calculate the Fourier coefficients Ac,s n in eq. (85). The integrand is periodic in 2π/N, so increasing the integration interval to 2π multiplies the result by N, hence we may express Ac, and by using the definitions of Fc,s (eqs. (20) and (21)) and the property of φ k from eq. (50), we get We interchange the sum and the integral and change variable ϕ ′′ = ϕ ′ + 2πk/N, obtaining In the above integral, fc,s has a periodicity of 2π and e −iN nϕ ′′ has a periodicity of 2π/N, therefore the integrand is periodic by 2π. We may therefore move the integration range to be between 0 and 2π, showing that the integral does not depend on k. After renaming ϕ ′′ to ϕ ′ we obtain because e −in2πk = 1 for any k. We see that the n Fourier coefficient of Fc,s is actually the Therefore, for large N we get the asymptotic Fc,s from its first Fourier coefficient (i.e. coefficients 1 and −1), so that we get from eq. (86): Because Fc,s are real, Ac,s −1 = Ac,s * 1 and we obtain Fc,s(ϕ ′ ) → N π |Ac,s 1 | cos [Nϕ ′ + arg(Ac,s 1 )] .
So we have to calculate the first Fourier coefficients of the Fc and Fs functions, Ac,s 1 .
The integrand being periodic on 2π, we may shift the limits by any value getting We change variable z = exp(iφ 0 ) to solve this integral on the complex plane, obtaining: where b is the pure imaginary number defined in eq. (43), C is the counterclockwise unit circle integration contour (see Figure (9 The integrand has two 2nd order poles, z 1,2 , defined in eq. (44) and only z 1 lies inside the integration contour C -see Figure (9). In addition there is an essential singularity at z = 0, because of the z −1 in the exponent.
We first calculate the residue at z = z1. Rewriting f (z) where Lc,s(z) is defined in eq. (45), we obtain We already showed that the first part is 0 -(see eq. (46)). This is because the integral in eq. (95) reduces for N = 0 to the integral in eq. (39) (up to 2π). So we are left with which comes out We start with Lc. To handle this derivative we express it as The last rational function is the Fibonacci polynomials generating function, with argument b −1 . Hence the result is (m + 1)! multiplied by the m + 1 Fibonacci polynomial. This may be directly calculated by factorizing the Fibonacci generating function to obtain or explicitly The Ls case is handled similarly. We express: The last rational function is the Lucas polynomials generating function, with argument b −1 . Hence the result is (m + 1)! multiplied by the m + 1 Lucas polynomial. This may be directly calculated by factorizing the Lucas generating function to obtain which results in By using eqs. (108), (104) and (102) and replacing b = −0.5ip we obtain a closed form expression for Ac 1 and for As 1 We remark that if N is a multiple of 4, the angle of Ac 1 is 90 o , and each increment of N removes 90 o from the angle of Ac 1 . Also we see that the angle of As 1 is always bigger by 90 o than the angle of Ac 1 and this is visible in Figures (3), (5) and (8). We shall name those angles ϕ c and ϕ s in the following calculations.
So by setting eqs. (92) in (55), and inverting the integration order between dp and dϕ ′ we obtain: The dϕ ′ integrals result in half the integration interval, i.e. π/N, so after simplifying we obtain G(β, N → ∞) = 3 4πNγ 4 β 2 N π β 0 dp p The asymptotic result is much easier calculable than the exact one, and does not require to solve for each step the implicit equation (32), but is still not given by a simple formula.
We will calculate in the next subsection an asymptotic expression for small β, i.e. G(β → 0, N → ∞), and for this case one arrives to a simple formula, as we shall see below.
A. Asymptotic result for many charges and low velocity For small β, |p| ≪ 1, hence 1 − p 2 in the denominator of eq. (113) may be set to 1, and 1 + 1 − p 2 ≈ 2, neglecting 1 − 1 − p 2 in eqs. (113) and (114). So we obtain (117) and Rearranging eq. (117) we obtain and defining K ≡ −4/(Np 2 ), one may express the sum over m in eq. (119) as a derivative with respect to K as follows and by changing the summation variable m ′ = n − m, this is written as We may sum exactly the last sum over m ′ to obtain n m ′ =0 where Γ with 2 arguments is the incomplete gamma function. We are interested in small |p|, so for |1/K| = Np 2 /4 ≪ 1, we obtain the approximated result given in eq. (122), which means that for small argument, the exponential series in eq. (122) needs very few terms to converge to an exponent. So continuing the calculation started in eq. (121) we obtain ∂ K K n+1 exp(1/K)] = K n exp(1/K)[n + 1 − 1/K] ≈ (n + 1)K n = (n + 1) where the last approximation used again the fact that |1/K| ≪ 1. Now using the result from eq. (123) in eq. (119), we obtain n + 1 (n + N + 1)! N n .
For large N, Γ(N + 3/2)/Γ(N) ≈ N 3/2 , so we obtain (130) Figure 13 shows the asymptotic results calculated with eq. (130) versus exact results according to eq. (55). It is to be mentioned that the case of small β and large N analyzed here fits the situations of currents in conducting materials or ion drift currents, mentioned in the introduction. In a conducting loop, the number of charges may be of order of 10 23 , and β may be of order 10 −12 , so that β 2N results in completely unmeasurable radiated power. In a ion drift device, the number of charges may be of order of 10 10 and β may be of order 10 −6 , so that the radiated power is somehow bigger than for the conducting loop, but still unmeasurable.

V. CONCLUSIONS
The purpose of this work was to learn how the radiation from discreet charges vanishes in the continuum steady state limit. We used a canonical configuration of charges in uniform circular motion, uniformly spread around a circle.
We found that the log of the power decreases almost linearly with the increase in the number of charges, and arrived to a close form solution to calculate this power if the number of charges is big -see Figure (