Polyharmonic Representation of the Electromagnetic Field Generated by an Oscillating Particle near a Dispersive Bulk

: This study introduces a polyharmonic framework for analyzing the electromagnetic (EM) field generated by an oscillating point charge near a dispersive bulk of size comparable to the wavelength under study. We critically evaluate traditional approaches such as Liénard-Wiechert, Landau


Introduction
In this study, we investigate the electromagnetic (EM) field generated by a stiff oscillating charged particle in proximity to a nanosphere by utilizing the diffracted field formalism.Our primary objective is to elucidate the EM field produced by an oscillating particle in vacuo.Despite its seemingly academic nature, this problem has practical applications in various configurations, such as Hertz and Sommerfeld antenna studies [1][2][3][4].Another example where this kind of source is used lies in the investigation of the potential generated by a charged point that is fixed in positionand oscillating in magnitude (as provided by ρ(x, t) = q exp(iω 0 t)δ(x)) in an inhomogenous cold plasma with no magnetic field [5].
In the context of nanoemitters and their interaction with plasmonic structures, Klimov et al. [6,7] treated an atom near a dielectric microsphere as a non-relativistic oscillator consisting of a stationary charge −q and an oscillating charge q.They highlighted the necessity of a quantum electrodynamics approach to ensure consistent calculation of photon emission, yet acknowledge that a classical electric dipole model can serve as a useful approximation [7].This perspective was echoed by Lobanov et al. in their numerical study of a point dipole with a periodic array of Yagi-Uda nanoantennas, where they stated that "The exact description of photon radiation from quantum nanoemitters is very complicated.A convenient approximation is the model of an oscillating point dipole."Similarly, Lassalle et al. modeled an excited two-level atom near a nanosphere as a harmonically oscillating point dipole [8].This approach was adopted in [9] as well, where the study of a point electric dipole in a flat-slab composite structure was conducted.Thus, as described in various foundational texts, it is common practice to employ a far-field dipole approximation [2,[10][11][12].
Remark 1.However, the dipole approximation is only valid when dealing with the far field and for non-relativistic particles.
Thus, it is necessary to find a better way to describe the field generated by an oscillating charge ρ(x, t) = qδ(x − a cos(ω 0 t)e z ).The immediate approach is to use the Liénard-Wiechert fields, which describe the fields produced by a moving charge in quite a compact form; however, their complexity renders them impractical from a numerical standpoint.Moreover, the methods proposed by Landau et al. involving a harmonic expansion of the delta and J.M. Raimond's Taylor series approximation of the delta distribution [13] are limited both analytically and numerically, as discussed in Section 2.2.2 of this work.
To address these limitations, we introduce a harmonic decomposition of the source terms ρ and ȷ.After establishing charge conservation, this method enables us to derive a harmonic representation of the electromagnetic fields.Our approach is distinguished by its analytical clarity and directness, allowing us to examine electromagnetic radiation without resorting to the common simplifications found in the existing literature, particularly the dipole approximation.
When the sources have been obtained in harmonic form, we use them as the incident field in the diffraction problem to obtain the diffracted field numerically based on a sphere and using the Finite Element Method (FEM) [14].This choice is motivated by the flexibility of FEM in handling various geometries of the dispersive bulk, although in this case we have chosen a sphere for simplicity.Furthermore, our solution is presented in the time domain, moving away from the simplifying assumption of working in the harmonic regime.
The rest of this paper is structured as follows.Section 2 lays the groundwork for our study, establishing the mathematical formalism of the diffracted field problem in Section 2.1.We then explore various methodologies for deriving the EM field produced by an oscillating charge in Section 2.2.This includes an analysis of the EM field generated by an arbitrary charge distribution and a critical evaluation of the Liénard-Wiechert, Landau, and Raimond approaches, highlighting their respective advantages and limitations.
In Section 2.3, we focus on developing a harmonic representation of the sources and demonstrating the conservation of charge.Section 2.4 is particularly significant, as it presents the polyharmonic construction of the electric and magnetic fields, complemented by numerical illustrations.
Section 3 presents the outcomes derived from employing the polyharmonic representation of the fields.Sections 3.1 and 3.2 detail the computation of the far-field radiated power from the polyharmonic fields, showing its alignment with the numerical integration of the Poynting vector around the trajectory of the oscillating charge.Furthermore, Section 3.3 establishes that the oscillating charge cannot exceed the speed of light, thereby upholding the principles of Einstein's theory of relativity.Building on these findings, Section 3.4 demonstrates how our polyharmonic solution encompasses the scenarios of a static monopole and a point dipole within a multipolar expansion framework.Section 3.5 is dedicated to applying the polyharmonic representation of the fields as the incident field in the diffracted field problem.This leads to a formulation involving spatial coefficients obtainable through the Finite Elements Method (FEM).A numerical illustration of the total field is provided in Section 3.6.
Finally, Section 4 summarizes the main contributions of this work, discussing its limitations and outlining potential avenues for future research.

Mathematical Description of the Problem
The mathematical description of the interaction of an EM field with a time-dispersive bulk is provided by the Maxwell equations: where ξ = {I, I I} represents the restriction of the fields, I corresponds to the total field outside the bulk, and I I corresponds to to the field inside the bulk.The constitutive relations are then H ξ = µ −1 0 B ξ , with µ 0 the magnetic permeability in vacuo (i.e., the dispersive bulk has no magnetization) and with ϵ 0 being the vacuum electric permittivity, δ(t) the Dirac's delta distribution, and χ(t) the causal electric susceptibility with its support within the bulk.The convolution is the one corresponding to the Fourier Transform convention in Appendix A. Thus, our new system of equations in terms of the (E ξ , B ξ ) fields reads: where c = (µ 0 ϵ 0 ) −1/2 is the speed of light in vacuo .The sources ρ and ȷ represent the charge density and current corresponding to an oscillating particle; their explicit expression is provided in Section 2.3.The fields generated by these sources in vacuo are: ∇ where the constitutive relations D 0 = ϵ 0 E 0 and H 0 = µ −1 0 B 0 have been used.Let us remark here that (E 0 ,B 0 ) is the incident field of our diffraction problem and that it is defined everywhere in space.The next step is to consider a new set of fields, defined as which, in the present linear case, satisfy the system of so-called diffracted fields: the sources of which are defined as and From these definitions, it is easy to see that charge conservation is satisfied.Note that the support of these new sources is within the bulk, and depends on the electromagnetic field generated by the oscillating particle.However, obtaining a handy expression of such a field requires very careful crafting, as shown in the next subsection.

The Search for the EM Field Generated by an Oscillating Charge In Vacuo
The problem of obtaining the electromagnetic field generated by a charged oscillating particle is a very important problem.This subsection begins with the very general approach of finding the (E, B) fields generated by an arbitrary pair of charge and current distributions (ρ, ȷ).Usually, these solutions are provided by the so-called Jefimenko's equations.However, we show that by following a similar method to the one in [15] for the far field approximations, it is possible to express the electric field mostly in terms of the current density.Next, we discuss a number of the efforts that have been made towards representing the EM field created by an oscillating charge.

The EM Field Generated by an Arbitrary Charge Distribution
First, we consider the system of Equations ( 1)-(4) (in the sequelae, the superscript 0 denoting the incident field radiated by the oscillating particle in free space has been removed to alleviate notations).By considering the Lorentz gauge , we can obtain the electric and magnetic fields (E, B) using the expressions where ϕ and A are the retarded potentials [10,11] with the retarded time t r = t − R c and R = |R| with R = x − x ′ .The next step consists of plugging the retarded potentials from Equations ( 13) and (14) into Equations ( 11) and (12) and then taking the appropriate derivatives.Nevertheless, this implies taking the time derivatives with respect to a function which is in terms of the retarded time t r .In order to avoid this difficulty and carry our calculations further, we propose using the Fourier transform in time (see Appendix A, f denoting the Fourier transform of f ) and then considering only the spatial derivatives.Starting by transforming the potential ϕ, we have where k 0 = ω c .In a similar way the vector potential in the angular frequency domain reads Before proceeding, the following identities are necessary: and by combining Equations ( 17) and ( 18) we obtain Equipped with these tools, it is quite easy to obtain B by simply taking the curl of Equation ( 15): after which, by taking the inverse Fourier transform, we arrive at the expression where B int is the intermediate field and B rad is the radiated field, which are respectively provided by the integrals In order to obtain the electric field, it is necessary to consider the Fourier transform of Equation (11): The first term on the right-hand side of Equation ( 21) is quite easy to calculate while the second term is much trickier.Its derivation is as follows.
Up to this point, most textbooks, e.g., [10,11,16], simply take the inverse Fourier transform of Equation ( 23) and provide the electric field in terms of the derivative of the density of charges, obtaining the so-called Jefimenko's equations [11,16].Despite the straightforward nature of this derivation, we are going to carry the calculations further.The reason for doing this is laid out in the following subsection.For our purposes, the second integral in Equation ( 23) is already in optimum form, and we focus our attention on the integral defined within a bounded volume Ω: Notice that the integral we are looking for is the limit case when Ω → R 3 .Next, we make use of the continuity equation in the angular frequency domain:

By defining the function
the integral I Ω in Equation ( 24) can be written in a more compact way (omitting the x ′ , R, and ω dependencies) as follows: Remark 2. Notice here that the divergence is being taken with respect to the primed coordinates (hence, the prime superscript in ∇ ′ ).
Due to the fact that we are working with Cartesian coordinates, it is possible to write where (e x , e z , e z ) denote the Cartesian unit vectors; then, we have Each one of these η integrals can be evaluated by means of the identity (η F and the Green-Ostrogradsky theorem [17,18], as follows: Taking the limit Ω → R 3 and following the assumption in [15] that ȷ is spatially bounded, the surface integral vanishes, leading to The gradient (with respect to the primed coordinates) can be computed explicitly: and we have Plugging Equation (26) into Equation (25), we obtain Before taking the inverse Fourier transform, we will try to express the term between square brackets (which we call ȷ) in a more illuminating way.First, we rearrange ȷ as Now, we consider the vector identity [18], and from this, we have Then, ȷ can be seen as Substituting this result into Equation (27) and then plugging that new integral into Equation (23), we finally arrive at the following expression.
From Equation (22), we can recognize the last integral as −iω Â; then, after taking the inverse Fourier transform, we obtain where E c is the Coulomb field, E int the intermediate field, and E rad the radiated field provided by the integrals At this point, the reader may be wondering about the reason we have used all these extra steps, when Jefimenko's equations already provide an explicit expression for the electric field.When dealing with Jefimenko's equations, the magnetic field is in terms of the electric current ȷ, and the electric field is in terms of the distribution of charge ρ [11,16].This point of view, albeit intuitively very clear, makes it difficult to compare the terms corresponding to the radiation field.By carrying out the manipulations described above, we ensure that the intermediate and radiated electric and magnetic fields are all expressed solely in terms of the electric current.The utility of this approach is shown in Section 3.1.

•
The Liénard-Wiechert's Field approach The academic problem of describing the E and B fields generated by a charged particle moving along a given trajectory u(t), called the Liénard-Wiechert fields, has been studied in many books [10,11,15,[19][20][21].The basic idea is to consider a charge density ρ and an electric current j provided by where δ is a Delta distribution and v(t) = du dt .From here, there are many ways to tackle the problem of obtaining the E and B generated fields: Jackson [10] and Landau [19] considered an elegant formalism using quadrivector approach, while Panofsky [15] and Heald [20] used the so-called Liénard-Wiechert potentials, which can be obtained by direct substitution in Equations ( 13) and (14) and then carrying all the necessary derivatives.The deduction of the fields E and B following this procedure can be seen in [11].For this subsection, we have decided not to follow any of these approaches, instead proceeding by direct substitution of the sources (31) and (32) into Equations ( 19), ( 20) and ( 28)- (30), that is: where we have introduced the usual following shorthand conventions: The procedure that is shown in Appendix C follows the ideas expressed by Heald and Marion in [21]; the main difference is that while Heald and Marion considered Jefimenko's equations, we use Equations ( 33) and (37).
Remark 3. We have decided to include the full deduction of the Liénard-Wiechert fields because, as far as we have seen, while this result has been quoted.the steps towards its obtention have not been shown.Griffiths simply states that the deduction is very difficult, and Heald and Marion say that it is necessary to perform heroic algebra .Therefore, we believe that it is important to show, as best we can, how to obtain one of the main results in classical electrodynamics.
When the fields E and B have been provided by Equation (A25) and Equation (A19), that is, and , with a = Kc .
β, it might be easy to think that the fields produced by an oscillating particle could be retrieved by considering the specific trajectory where a and ω 0 are the oscillation amplitude and angular frequency, respectively; nevertheless, as pointed out by Spohn in [20], the Liénard-Wiechert fields are less explicit than they appear to be.This is due to the fact that Equation (A25) and Equation (A19) depend on the retarded time, which is itself a solution of a transcendental equation (in general, a non-trivial one), namely, Let us note here that if the particle is at constant speed with a straight trajectory, then the Liénard-Wiechert fields can be found almost straightforwardly [11].However, the solution to this problem when dealing with an oscillating charge, in this case the retarded time, is a function of the present time and the position (t r := t r (t, x)).

•
The Landau Spectral Resolution Approach In the The Classical Theory of Fields , Landau et al. [19] considered that the fields produced by moving charges can be expanded as a superposition of monochromatic waves.Assuming that ρ(x, t) and ȷ(x, t) have a Fourier integral representation, we can write According to Landau, it is clear that each Fourier component of ρ(x, t) and ȷ(x, t) is responsible for the creation of the corresponding monochromatic component of the field.Thus, it is natural to consider the following Fourier integral representations of the potentials ϕ(x, t) and A(x, t): A(x, t) = ω∈R Â(x, ω)e −iωt dω.
In the sequelae, we only work with ϕ(x, t), as all the discussion applies to A(x, t) as well; substituting Equation (38) into Equation ( 13), we obtain Comparing this last equality with Equation ( 39), we find that Remembering that ρ(x, ω) is the Fourier transform of ρ(x, τ), after some manipulations, φ(x, ω) can be written as Now, we consider the singular charge distribution, as in Equation ( 31), to obtain Upon substitution of this expression into Equation (39), we arrive at and similarly, for A(x, t), we have From the above expressions, it is evident that the right-hand side of the equation explicitly depends on the present time.This feature avoids issues related to the retarded time, which is a notable problem when using the Liénard-Wiechert fields.However, difficulties emerge when attempting to explicitly compute these integrals.This arises due to the fact that the particle's trajectory, denoted by u, is in principle unrestricted in its choice of any argument τ.Moreover, the exponential function in Equations ( 40) and ( 41) necessitates a sweep over all the (ω, τ) values in R 2 .• The Raimond's Taylor Series Expansion Approach In [13] J.M. Raimond proposes another way to represent the charge density ρ(x, t) = qδ(x − a cos(ω 0 t)e z ) using a Taylor series expansion of the Dirac delta in the sense of distributions: where g(t) = a cos(ω 0 t).In this way, the charge density can be rewritten: where the ρ n (x, t) is the n-th charge density and ϱ n (x) is provided by The reader may identify ϱ n (x) as the singular charge distribution of a 2 n -pole (i.e., n = 0 monopole, n = 1 dipole, n = 2 quadrupole, etc.) [22,23].Now, by considering the charge conservation, we can write: From this last expression, it can be seen that the current distribution ȷ(x, t) can be written as where ȷ 0 (x, t) = 0 and ȷ n (x, t) = g ′ (t)ρ n−1 (x)e z (n > 0).From the above expressions, it is easy to see that there is a conservation of charge between the n-th current density ȷ n (x, t) and the n-th charge density ρ n (x, t) for n ≥ 0. Thus, it is possible to consider the charge distribution pairs (ρ n (x, t), ȷ n (x, t)) in order to obtain the fields generated by an oscillating charge.In this case, the retarded potentials are In the same fashion, we can write the vectorial magnetic potential as The potentials in Equations ( 42) and ( 43) are expressed in terms of the present time t.However, the multiple derivatives that must be computed render the expressions impractical.One could argue that for a large value of ||x||, only a few terms are necessary to obtain a good approximation, as demonstrated in [13] by retaining up to the dipole term.However, this approach would essentially involve another far-field approximation [10].In the following subsection, we present a more practical and elegant method for representing an oscillating charge.

Harmonic Decomposition of the Sources
As seen in the previous subsection, the Liénard-Wiechert fields are not the best way to obtain the fields produced by an oscillating particle.The main problem is that the source terms depend on the trajectory that describes the charge.For this reason, it is convenient to find a way to decompose the source terms in a polyharmonic way.This idea has been previously considered by Landau [19] and, as evident from Equations ( 42) and (43), by Raimond [13], albeit indirectly.Remark 4.However, as demonstrated in the previous subsection, the expressions derived from Landau's and Raimond's ideas are challenging to implement in practice.Thus, a new approach for describing the sources is necessary.
In this subsection, we propose another approach inspired in quantum mechanics, which can be summarized as follows: The superposition of waves spread in a certain domain can be seen as a particle .Physically, this means that a very localized source can be seen as the interference of a certain kind of waves.Mathematically speaking, we are looking to find a sequence such that we can have convergence in the sense of distributions to a Dirac delta [24][25][26].It is worth noting that a similar concept was employed in [27,28] for the case of a charged particle in uniform motion.

Two Fourier Expansions for the Sources
Let us start our analysis by providing the charge density ρ and its current density ȷ simply by with ϱ(z, t) := δ(z − a cos(ω 0 t)) and and with ȷ(z, t) := −aω 0 sin(ω 0 t)ϱ(z, t) out of charge conservation.Notice that ϱ and ȷ are not multiplied by the charge.It then turns out that the charge density is not harmonic, despite the harmonic motion of the particle.In other words, no complex function ϱ(z) can be found such that ϱ(z, t) = Re{ϱ(z)e iωt }.
Remark 5.This is quite important, as in references such as [10,22] this is the starting point when representing an oscillating dipole.
Nevertheless, it is appropriate to notice that the distribution ϱ(z, •) is a 2π ω 0 -periodic distribution.Thus, ϱ and ȷ(z, t) can be expanded as a Fourier series (See Appendix B): where T l are the Chebyshev polynomials of the first kind [29] and w(z) is the weight function with χ [−a,a] (z) a characteristic function.

Continuity Equation for the Harmonic Components of the Sources
In the previous paragraph, an expansion for ρ and ȷ was linked by the so-called charge conservation, namely, ∇ • ȷ + .ρ = 0. Notice that for moving point particles this equation has to be understood in the sense of distributions [17,[24][25][26].However, what about the different components ϱ F l and ȷ F l ?In other words, is there any transference of energy between the different waves oscillating with the different frequencies at stake ω 0 , 2ω 0 , etc. . .?To answer this question, we have to care much more about the notation referring to l.While ϱ F l oscillates with angular frequency lω 0 , the scalar function ȷ F l is a mix of two oscillations with different frequencies, namely, (l − 1)ω 0 and (l + 1)ω 0 , due to the presence of the sin(ω 0 t) term in Equation (46).Making use of the complex representation of sin(ω 0 t), it follows that )[e i(l−1)ω 0 t − e i(l+1)ω 0 t ].
As can be seen, this representation is not in a convenient form; instead, we would like to see each term in the series oscillating at angular frequency lω 0 (where l is a dummy index), namely, This can be easily achieved after renaming indices (l + 1 → l and l − 1 → l) for the corresponding terms e i(l−1)ω 0 t and e i(l+1)ω 0 t , which allows us to obtain where ξ l (z) = w(z)T l ( z a ).Analogously, for ρ(z, t) we obtain Thus, the arcane meaning of the superscripts T and F is clear: T (resp.F) means true (resp.false) in the sense that each spatial coefficient corresponds to only one l ∈ Z.
Correspondingly, the conservation of charge can be now formulated for each multiple of the angular frequency ω 0 .From the definition of ρ(x, t) and ȷ(x, t), the conservation of charge implies that and by plugging in the harmonic expansions of ϱ(z, t) and ȷ(z, t) we obtain In light of the fact that e +ilω 0 t with l ∈ Z is a basis [24,25] for our polyharmonic decomposition, all the terms between square brackets are equal to zero.Thus, the following identity is obtained: and the conservation of charge for each l-th term in the expressions Equations ( 47) and (49) follows.Upon demonstrating that there is no transfer among the distinct harmonic components of the charge density ρ(x, t) and the current density ȷ(x, t), we may proceed to employ these Fourier expansions to derive the electric and magnetic induction fields for an oscillating charged particle.This is shown in the following subsection.

The Building of E and B via Polyharmonic Computations
The main consequence of the term-by-term harmonic decomposition of the sources and the conservation of charge is that the electric and magnetic induction fields, respectively E and B, can be seen as a superposition of the fields E l (x, t) and B l (x, t), that is, fields generated by the harmonic densities of the charge ϱ T l (z)e +ilω 0 t and current ȷ T l (z)e +ilω 0 t .Each one of these fields oscillates in terms of multiples of the fundamental angular frequency ω 0 .This subsection is devoted to this issue, starting work with Equations ( 19), ( 20) and ( 28)- (30) and considering R = x − x ′ , R = |R|, and t r = t − R c the retarded time.
Remark 6.In this case, the retarded time is not in terms of a transcendental equation, and is instead provided explicitly in terms of the present time t.

Fourier Expansion of the Fields E and B
We start by applying the delta distribution δ(x ′ ⊥ ) from Equations ( 44) and ( 45) to the expressions (19) and ( 20) and ( 28)- (30), which implies that E and B are provided by Note that in the previous equations and the sequelae, the explicit dependencies on x, z ′ and t are omitted.The next step is to define the functions allowing the electric and magnetic fields to be written in a more compact way as From the definitions of ϱ(z, t) in Equation (49) and ȷ(z, t) in Equation ( 47), we know that Equations ( 50) and (51) read where Notice that there is a phase shift e −ilk 0 R in these expressions due to the retarded time t r .Therefore, E(x, t) and B(x, t) can be seen as a superposition of elementary harmonic terms, i.e., with spatially dependent coefficients provided by where we have used the fact that the support of ϱ T l (z ′ ) and ȷ T l (z ′ ) lies within the interval [−a, a].These coefficients can be computed numerically; thus, the building of the E and B fields is complete.

A Geometrical Description of the Fields
Despite the complicated appearance of Equations ( 52) and (53) it is possible to extract a priori information about the geometric behaviour of the (E,B) fields.Starting with the vector identity (e z × R) × R = (e z • R) R − R2 e z , we can rewrite the spatial coefficients as where Φ sph l and Φ f lat l are defined by Finally, the spatial coefficients for the magnetic induction field read From these representations, it is easy to see that the first integral term of E, which is called E sph , is a field with spherical symmetry.However, due to the action of the second integral term, in the sequel E f lat the total field is flattened out in the direction perpendicular to the motion.On the other hand, the field lines of B circle around the z-axis, and as expected are perpendicular to the field lines of the electric field.
The only drawback we have encountered is that E 0 l might be singular between (0, 0, −a) and (0, 0, +a). Figure 1 shows cuts along the planes x = 0 and y = 0 of the imaginary part of the harmonic electric field components E l for l = 1, 2, 3, 4. For l = 1 (cf.Figure 1a), it is possible to recognize the typical shape of the Green function (or point dipole along z).The real part of the harmonic magnetic field components H l = H l /µ 0 is shown in Figure 2, which, as expected, circles around the z-axis.Finally, the projection of the Poynting vector field on the canonical planes is shown in Figure 3. Again, it is possible to recognize the radiation pattern of the point dipole in Figure 3a.The computation of the electric and magnetic fields has been carried out by means of Equations ( 52) and (53).In order to compute these integrals numerically, the following change of variables is used: z ′ = a sin(θ).This allows the term 1 √ a 2 −z ′2 , which comes from the definition of w(z ′ ) (see (A3)), to be eliminated.With this change of variable, the formulae were coded into GetDP by P. Dular and C. Geuzaine (Liege, Belgium) open-source software [30], from the source available at https://gitlab.onelab.info/getdp/getdp/-/tree/oscillating_multiharm(accessed on 9 November 2014), using a simple trapezoidal rule with 300 integration points.Visualization was then performed using Gmsh v4.12.1 open-source software [31].

Far-Field Radiated Power
As was mentioned before, the main approach towards the study of the EM field generated by an oscillating charge requires consideration of the far field approximation .In this subsection, we demonstrate how it is possible, using the polyharmonic expression from the previous subsections, to retrieve the same results reported in the literature [10].

The Polyharmonic Representation of the EM Far Field Approximations
Remembering Equations ( 20) and (30), the radiated electric and magnetic induction fields are provided by where n = R R .Using the fact that the electric current is defined as ȷ(x ′ , t r ) := qδ(x ⊥ ) ⊗ ȷ(z, t r )e z , the fields read and, as before, R = x − z ′ e z , ñ = R R , and tr = t − R c .Now, in order to compute the radiated fields ad infinitum, we make the following approximations for when |x ′ | ≪ |x|: By plugging these approximations into Equations ( 54) and (55), we arrive at the expressions From this last equation, it is clear that we can focus our attention on just B ∞ .Therefore, the next step is to consider its polyharmonic representation.Remembering the definitions of ȷ(z ′ , τ r ) in Equation ( 47) and τ r in (56), we know that with the spatial coefficient B l ∞ defined as Here, it is important to emphasize that B 0 ∞ = 0; thus, the results obtained in this subsection are for l ̸ = 0. Notice that the integral term resembles a spatial Fourier transform that goes from z ′ → η l , with the new variable This Fourier transform can be easily computed using the definition in Equation ( 48): Fortunately, the Fourier transform of ξ l is provided in Equation (A5), and after using identity (A27) we obtain where we have defined the dipole moment vector p z = qae z .Notice that this last expression is very similar to the one in [10] for the magnetic dipole field in the radiation zone.Calling A l = i l−1 µ 0 4π c(lk 0 ) 2 e −ilk 0 |x| , we obtain and its square norm is provided by where we have used the fact that η l = lk 0 z |x| .This last result is used later in the study of the radiated power.

Radiated Power
The expression for the far-field radiated power for an oscillating particle is now derived, first for the relativistic case and later for the non-relativistic one.

•
Relativistic case.We know from the definition of the Poynting vector and (57) that Remembering the Fourier expansion of B ∞ in Equation (58), we have Now, we can consider the time-averaged Poynting vector; in light of the fact that B ∞ is 2π ω 0 -periodic, we integrate Equation (60) from − π ω 0 to π ω 0 to obtain where the quantities between the brackets are time-averaged.From the orthogonality of the complex exponential, we have When obtaining the power, it is customary to calculate the flux of the averaged Poynting vector through a spherical surface of radius R 0 (while this could be any surface encompassing the EM sources, the spherical surface is the one that allows the simplest calculations), that is, Next, we perform the change of variable: after which the differential du reads du = −alk 0 sin θ dθ.
After these calculations, the flux of the averaged Poynting vector reads Calling where Z 0 = µ 0 c is the impedance of free space, because D −l = −D l we can obtain the averaged energy flux as follows: with P l as the l-th contribution provided by where, by means of the identities obtained in Appendix D, we can see that the integrals are provided by where the function I R 0 is defined trough an integral (see (A37)).Employing Equations ( 65)-( 67), it is possible to compute the radiated power flux in a semianalytical manner.In particular, if we focus our attention on the case where l = 1.
or in more explicit fashion, Remark 7. It is important to note that in this case we have made no restrictions regarding the velocity of the charged particle, albeit it cannot be faster than c (this is studied in greater depth in Section 3.3).Thus, the expressions derived here are valid even for charges with speeds near c.
Let us recall that the position of the charged particle follows the law a cos(ω 0 t), which means that the maximum speed attainable by this particle is aω 0 .In the context of the non-relativistic regime, this necessitates that aω 0 ≪ c, or equivalently, alk 0 ≪ l.For l = 1, we have 0 < u < ak 0 ≪ 1.In this case, the behaviour of J 1 near the origin is well known, namely, J 1 (u) ∼ u 2 ; as a result, we have and obtain The reader may recognize this last equality as the expression used in [10] for the total radiated power.Finally, Figure 4 shows a comparison between the power P l (for l = 1, 2, 3) obtained analytically by means of (64) when using the far field approximation and the power computed numerically by considering the Poynting vector flux across the surface of a pill box that encloses the trajectory of the charged particle, as depicted in Figure 5b.The Poynting vector is obtained via Equations ( 52) and (53), and the numerical integrations were performed using the GetDP solver [30].As can be seen, both approaches fit perfectly.

The Particle Cannot Be Supraluminal
It is common sense that for the l-th component of the magnetic induction field in Equation ( 58 57)), and the power P l (provided by Equations ( 61)-( 64)), the value converges with growing l.A divergence would result in a total value of infinity.Due to the connection between B l ∞ , ⟨S l ∞ ⟩ and P l , it follows that if the magnetic field diverges, then the two others diverge as well.
The l th magnetic field depends on l only in For large l, J l (lβ) can be written as follows [32]: Combining this result with Equation (59), it is possible to see that B l ∞ behaves as This only converges with l → ∞ if the argument of the exponential function is negative, which, in the form of Equation ( 69), leads to the following inequalities.
With v = aω 0 being the maximum velocity of our point charge in this sinusoidal movement, the condition for convergence is matched in the physical sense of Einstein's postulate that nothing is faster than light [10,11].
Remark 8.This result holds irrespective of whether the particle has a mass.

Multipolar Expansion of the Fields
The results presented in the preceding subsection indicate that the radiated power calculated using the polyharmonic representation of the EM field aligns with the results derived from the oscillating point dipole approximation, when the far field approximation is used .This subsection aims to elucidate the reasons behind this phenomenon and to provide a method for estimating the error associated with this approximation.In order to illustrate this method, we examine the Coulomb field outlined in Equation ( 28), applying the polyharmonic representation of the charge density (ρ(x, t)) as defined in Equations ( 44) and ( 49), that is, To simplify our analysis, it is convenient to define the following function: where r = |x|.Upon substitution of Equation (71) into Equation (70), we obtain The reader can identify the Type 1 integral as defined in Appendix F; however, our analysis requires summing only for l ≥ 0. To facilitate this, we note the following relationships: where the overline symbol denotes the complex conjugate.Consequently, Equation (72) can be reformulated to incorporate these observations, as follows: The integrals presented in Equation (73) can be reformulated using a vector Taylor series (multipole expansion), as outlined in Equation (A58): where N represents the highest term in the Taylor expansion and d l,n is defined as follows.
Equation ( 74) establishes a selection rule that links the l-th multiple of the angular frequency ω 0 with the n-th term of the Taylor expansion (2n-pole).This relationship is more clearly demonstrated in Table 1, where the truncation term is considered with N = 3.
Table 1.Compatible pairs of l and n, indicated by • symbol, when N = 3.

This table exemplifies a broader principle.
Remark 9. When employing multipole expansion to represent the electromagnetic field generated by an oscillating charge, it becomes apparent that the 2n-pole order of our expansion limits the observable oscillating frequencies lω 0 .These angular frequencies must adhere to the rule requiring that n ≥ l and n + l be even.
Therefore, the Coulomb field E c (x, t) reads Our next objective is to conduct a similar analysis, this time focusing on the radiated field as described in Equation (30).For this purpose, we use the polyharmonic representation of the current density, ȷ(x, t), as detailed in Equations ( 45), (47) and (48).The expression for the radiated field is provided by By defining the function we can reformulate Equation ( 76) as To ensure that the summation includes only terms with l ≥ 0, we observe the following relationships: .
With these relationships in mind, Equation (78) can be rewritten as By considering Equation (A62), we can reformulate Equation (79) as follows: Here, d l+1,n and d l−1,n follow the same rule established in Equation ( 74).After conducting elementary manipulations, Equation (80) can be further simplified to The methodology applied to the E rad (x, t) field can be similarly employed for the E int (x, t), B int (x, t), and B rad (x, t) fields.This is achieved by considering the following functions: By adopting these formulations, we achieve a more lucid representation of the electromagnetic field generated by an oscillating charge.This approach enhances our understanding of the field's representation as derived in this work as well as of the multipole expansion method.In the following subsection, we demonstrate how the point dipole approximation is inherently encompassed within our more general expression.
Retrieving the Oscillating Point Dipole: Taking N ≤ 1 Our aim here is to derive the classical time-dependent point dipole expressions found in standard texts such as [33], specifically under the context of harmonic motion.Additionally, we provide expressions of the truncation error for both the electric and magnetic fields.To illustrate our approach, let us start with the electric field E c (x, t) as described in Equation (75), focusing on the case where the truncation term N is set to 1.This choice simplifies our analysis and allows us to directly compare our results with the established point dipole model.
=F 0 E,c (x) + 2Re e iω 0 t a 2 D (1) For this last equality, we have applied the selection rule as outlined in Equation ( 74), particularly focusing on its specific cases where d 0,0 = 1 and d 1,1 = 1 2 .Remembering that F 0 E,c (x) is defined in Equation ( 71), we have e z [e ik 0 r x r 3 ] + V E,c , where is the truncation error.Next, we apply the differential operator D e z , as defined in Equation (A45): with n = x r .Substituting this result into Equation (85), after some elementary manipulations we obtain where p z = qae z .For the E rad (x, t) field, considering a truncation with N = 0, Equation (81) simplifies to It is important to note that the double sums in Equation (81) vanish due to our choice to set N = 0. Recalling the definition of F 1 Erad in Equation (77), we obtain where Applying the same methodology to the E int (x, t), B int (x, t), and B rad (x, t) fields and setting the truncation term to N = 0 leads us to derive the following concise representations of these fields: B,int , (90) with the corresponding truncation errors: B,rad =2Re Synthesizing our findings, we can express the electric and magnetic fields with truncation term set to N ≤ 1 as follows: and The expressions derived here may be familiar to the reader as those corresponding to a point oscillating dipole [10,33]; however, it is important to note the inclusion of an additional term representing the Coulomb field of a static point charge located at the origin.This inclusion aligns with our expectations, considering that the zeroth order in the Taylor expansion of qδ(x − a cos(ω 0 t)e z ) is qδ(x), a component already encompassed within our polyharmonic representation of the electromagnetic field.Furthermore, our approach allows us to articulate the truncation error associated with the point oscillating dipole approximation.This highlights the comprehensive nature of the polyharmonic representation of electromagnetic fields developed in this work.

Obtaining the Diffracted Field
After all this work, we find that the fields E 0 and B 0 of the system of Equations ( 1)-( 4) can be regarded as a superposition of waves in the form of Equations ( 52) and (53).Even more, the analysis of the radiated energy from the previous subsection shows that it is possible to retrieve the classical results from our expressions of E 0 and B 0 when the nonrelativistic far field approximations are considered.The second part of this work, which is going to be considerably shorter that the first part, can be tackled in straightforward fashion.First, the sources in Equation ( 9) can be easily retrieved by noting that It is important to notice that εr,ξ (lω 0 ) is the complex conjugate of the Fourier transform of ϵ r,ξ (t); thus, It is then natural to propose the following as solutions of the system (1)-( 4): Plugging these solutions into Equations ( 5)-( 8) and recalling that E 1 ξ,l := E ξ,l − E 0 , we arrive at the following system which must be satisfied for each l ∈ Z: Taking the curl on Faraday's Law, we obtain where the right-hand side of this expression is a source term.Moreover, the equation can be solved, for instance, using the Finite Element Method, as explained in [34,35].

Numerical Illustration
In order to illustrate the above discussion, we present the following numerical results for the case of an oscillating charged particle close to a sphere (this method can be readily applied to more complicated geometries, as finite elements are used).For the purposes of this example, we consider that the particle oscillates at an angular frequency ω 0 = 2πc/λ 0 with λ 0 = 900 nm.The particle oscillates between points of the coordinates (0, 0, −a) and (0, 0, −a) with a = λ 0 /7.The radius of the sphere is λ 0 /6 and its center is located at (0, λ 0 /5, a).The relative permittivity of the sphere is set to 9 + i.It is important to note that for the sake of simplicity we have used a non-dispersive sphere here; however, it would have been simple to consider temporal dispersion with the proposed method.Finally, all geometries and conformal meshes were obtained using Gmsh software [31], and all finite element formulations in this article were implemented thanks to the flexibility of the GetDP finite element software [30].Perfectly Matched Layers (PMLs) [36] were used to truncate the surrounding free space.The incident field was hard-coded in GetDP as well, as detailed in Section 2.4.1.
The mesh size (see Figure 6) was set to λ 0 /20, which is very fine at ω 0 and reasonably coarse at ω 4 = 4ω 0 ↔ λ 4 = λ 0 /4.The 3D scattering problem uses high order Webb hierarchical edge elements [14,37,38] with twenty unknowns per tetrahedron (two unknowns per edge, two unknowns per face).The direct problem corresponding to Equation (97) was solved using the MUMPS direct solver [39], which natively interfaces with GetDP.The number of unknowns of the final discrete system in this example is about 1.8 million, which uses 130 Gb of RAM memory for a resolution time of 20 min per harmonic l on a workstation equipped with of Intel Xeon Platinum processors (2.50 GHz) and with a multithreaded version of MUMPS v5.5.1 running on 32 threads.Figure 7 shows the harmonic components of the diffracted field E 1 l for l = 1, 2, 3, 4. Notice that this field was obtained as a numerical solution of (97) using FEM.Moreover, due to the fact that the support of the source term on the right-hand side of Equation ( 97) is within the sphere, the possible singularity of E 0 l does not affect our numerical results.Figure 8 shows the total electric field (E l = E 1 l + E 0 ) and its interaction with the sphere, whereas the total Poynting vector has been computed as S l = 1 2µ 0 E l × B l and can be seen in Figure 9 for l = 1, 2, 3, 4. It is important in the case of the Poynting vector to see how this is somewhat pulled by the sphere.This is due to the passivity of the material; remember that the permittivity of the sphere is 9 + i.Finally, all of our results have been corroborated by considering an energy balance that measures the total energy flux that crosses a pill-box that surrounds the sphere.This is shown in Figure 5.

Conclusions
In this study, the electromagnetic field generated by an oscillating point charge has been thoroughly examined.We have explored and analyzed the methodologies of Liénard-Wiechert, Landau, and Raimond, discussing their respective merits and limitations.Our proposed Fourier representation of the sources has led to the development of a polyharmonic framework for constructing electromagnetic (EM) fields.This approach offers an analytical and practical representation of the EM fields, effectively circumventing the complexities inherent in the previously discussed methods.
Our polyharmonic technique is applicable to both relativistic and non-relativistic scenarios, yielding more comprehensive formulas for the radiated power.Under specific simplifying conditions, these formulas align with the well-known far-field time-dependent dipole radiation power equations frequently referenced in the literature.Furthermore, our multipolar decomposition demonstrates that our solution includes both the Coulomb field and the oscillating point dipole fields, which are typically employed to represent the EM field of an oscillating charge.We elucidate why such approximations fail to capture the more intricate polyharmonic fields with angular frequencies.
The polyharmonic representation is particularly advantageous for studying the interaction between a stiff oscillating charge and a dispersive bulk in the time domain.To illustrate this, we have chosen to focus on a dispersive nanosphere, employing the Finite Element Method (FEM) for our numerical solutions.This choice is justified by the flexibility and adaptability of FEM, which accommodates dispersive bulks with more complex geometries.Our results are validated through an energy balance, ensuring their accuracy and reliability.
It is important to note, however, that our approach has certain limitations.For instance, the impact of the diffracted field on the charge trajectory and the Abraham-Lorentz force are not considered in this study.Future work will address these limitations, along with exploring different configurations of materials, oscillation frequencies, and geometries.
where t l are the zeros of the function f (t), i.e., cos(ω 0 t) = z a .
In light of the fact that the absolute value of the cosine is bounded by 1, there are only two solutions within the range [−a, a].Then, by considering the principal branch of the arc cosine function and denoting t 0 (z) := 1 ω 0 arccos( z a ), these zeros are t + l (z) = t 0 (z) + 2lπ/ω 0 , t − l (z) = −t 0 (z) + 2lπ/ω 0 .
Taking the Fourier transform of Equation (A4) and keeping in mind that these two Fourier series are equal term-by-term, we arrive at this beautiful and unexpected expression: which is used in Section 3.1.
Now, it is necessary to compute the derivatives with respect to the present time t.Assuming the convention that the doted quantities denote partial derivation with respect to time t, the following identities will be useful: From Equations (A10) and (A11), we can solve for and upon substitution of Equation (A12) into Equations (A10) and (A11) we obtain .
From the expression Rn = R, we can perform derivation with respect to t; after some manipulations, we obtain In addition, we have .
where a = .
v denotes the particle acceleration.Finally, Using these identities, we first deal with B in Equation (A8): as defined in Equation (A35).Next, we obtain an expression for I R 1 (A) by considering Equation (A30) with l = 0, d dr r d dr J 0 + J 0 r = 0.
Remembering that d dr J 0 = −J 1 , we obtain − d dr rJ 1 + J 0 r = 0, and after multiplying by dJ 0 dr = −J 1 , we have J 1 d dr rJ 1 + dJ 0 dr J 0 r = 0, which implies After integration by parts, we have This expression can be arranged in a more illuminating way: Therefore, at the end Equation (A35) depends only on An analytical expression for this integral is provided in Section 2.1.3 of [46], as follows: where x = A in this case.Alternatively, this integral can be accurately evaluated numerically using the periodisation method detailed in [47].

Figure 1 .
Figure 1.Harmonic components of the imaginary part of the electric field (in V/m) generated by an oscillating particle, E l for l = 1, 2, 3, 4 (cuts along the planes x = 0 and y = 0).The size of the arrows is redundant with the colorscale.

Figure 2 .
Figure 2. Harmonic components of the real part of the magnetic field (in A/m) generated by an oscillating particle, H l = B l /µ 0 for l = 1, 2, 3, 4 (cuts along the planes x = 0 and y = 0).

Figure 3 .
Figure 3. Harmonic components of the real part of the Poynting vector field (in W/m 2 ) generated by an oscillating particle, S l for l = 1, 2, 3, 4 (cuts along the planes x = 0 and y = 0).

Figure 4 .
Figure 4. Comparison of the powers P l (l = 1, 2, 3) obtained analytically (thick lines) and numerically (crosses), normalized by the quantity P ref 0 obtained analytically for l = 1 and λ 0 /a = 1.The analytical data were obtained by numerical integration of the formula in Equation (64).The numerical data were obtained by 3D numerical intgration of the Poynting vector shown in Figure3.

Figure 5 .
Figure 5. (a) Pill-box surface (green lines) enclosing the sphere (blue lines) on which is plotted the Poynting vector.The trajectory of the particle is indicated by the red line.(b) Zoom on the pill-box and sphere.The yellow-green map represents the real part of the total Poynting vector for l = 0.The blue-yellow map represents the square norm of the total electric field for l = 0.

Figure 6 .
(a) Mesh showing the different regions of integration and (b) zoom enclosing the sphere and the region where the particle oscillates.

Figure 7 .
Figure 7. Harmonic components of the imaginary part of the diffracted electric field (in V/m) E 1 l for l = 1, 2, 3, 4 (cut in the plane x = 0).For clarity, data in close vicinity to the particle trajectory (elongated rectangle below the sphere) are not displayed.

Figure 8 .
Figure 8. Harmonic components of the imaginary part of the total electric field (in V/m) E l for l = 1, 2, 3, 4 (cut in the plane x = 0).for clarity, data in close vicinity to the particle trajectory (elongated rectangle below the sphere) are not displayed.

Figure 9 .
Figure 9. Harmonic components of the real part of the total Poynting vector field (in W/m 2 ) S l for l = 1, 2, 3, 4 (cut in the plane x = 0).For clarity, data in close vicinity to the particle trajectory (elongated rectangle below the sphere) are not displayed.