Second gradient electromagnetostatics: electric point charge, electrostatic and magnetostatic dipoles

In this paper, we study the theory of second gradient electromagnetostatics as the static version of second gradient electrodynamics. The theory of second gradient electrodynamics is a linear generalization of higher order of classical Maxwell electrodynamics whose Lagrangian is both Lorentz and U(1)-gauge invariant. Second gradient electromagnetostatics is a gradient field theory with up to second-order derivatives of the electromagnetic field strengths in the Lagrangian. Moreover, it possesses a weak nonlocality in space and gives a regularization based on higher-order partial differential equations. From the group theoretical point of view, in second gradient electromagnetostatics the (isotropic) constitutive relations involve an invariant scalar differential operator of fourth order in addition to scalar constitutive parameters. We investigate the classical static problems of an electric point charge, and electric and magnetic dipoles in the framework of second gradient electromagnetostatics, and we show that all the electromagnetic fields (potential, field strength, interaction energy, interaction force) are singularity-free unlike the corresponding solutions in the classical Maxwell electromagnetism as well as in the Bopp-Podolsky theory. The theory of second gradient electromagnetostatics delivers a singularity-free electromagnetic field theory with weak spatial nonlocality.


Introduction
In recent years, there has been a continuous interest in the so-called Bopp-Podolsky theory [1,2,3,4,5,6,7], a first-order linear gradient theory of electrodynamics. The theory was introduced in the early 1940s by Bopp and Podolsky [8,9,10] as a method to remove singularities in classical Maxwell electrodynamics and to obtain a consistent Lorentz-and gauge-invariant theory of point charges with finite self-energy (see also [11]), thereby proposing an alternative to earlier theories that achieved the same goal through a nonlinear generalization of electrodynamics [12]. While the motivation was ultimately the quantization of the theory (e.g. [13]), it is, first of all, a classical field theory. Through the introduction of an additional gradient term in the Lagrangian the generalized Maxwell equations yield linear partial differential equations of fourth order for the electromagnetic potentials. Along with the additional term a new constant has to be introduced, the Bopp-Podolsky length scale parameter , which by the original idea was supposed to be related to the electron self-energy. Indeed, Iwanenko and Sokolow [14], Kvasnica [15] and Cuzinatto et al. [16] reasoned that the Bopp-Podolsky length scale parameter is of the order of ∼ 10 −15 m which is the order of the classical electron radius, however, Accioly and Mukai [17], and Carley et al. [3] argued that the Bopp-Podolsky length scale parameter should be equal to or smaller than ∼ 10 −18 m. From the mathematical point of view the length scale parameter plays the role of a regularization parameter. The regularization through the higher-order field equations provides a finite self-energy for the point charge, and its electrostatic potential is finite and nonsingular, its electric field finite, but with a directional discontinuity. Also, in Bopp-Podolsky electrodynamics the electromagnetic fields of a non-uniformly moving point charge possess a directional discontinuity on the light cone. In order to obtain electromagnetic fields of a non-uniformly moving point charge with no directional discontinuity on the light cone, the theory of second gradient electrodynamics has recently been proposed and used in [18]. Further important advantages of gradient electrodynamics in comparison to the classical Maxwell electrodynamics are: no infamous "4/3-problem", no unphysical runaway solutions to the equation of motion of a moving charged particle (the analogue of the Lorentz-Dirac equation), no singularities of the electromagnetic fields on the light cone, the self-force of a non-uniformly moving charged particle is regular (see also the discussion in [18]).
Besides the point charge, other important textbook examples of the Maxwell theory are the electrostatic and magnetostatic dipoles possessing typical dipole singularities (1/r 3 -and Dirac delta-singularity) in the electromagnetic fields [19,20,21,22]. In the mathematical literature, there is an interest in the regularization of the dipole singularities arising from the second-order derivatives of 1/r in the sense of generalized functions (see, e.g., [23,24,25]). The ideal magnetostatic dipole in first-order gradient electrodynamics was already studied by Landé and Thomas [26], giving the magnetic fields and a finite selfenergy. The latter result, however, turns out to be erroneous. Using the Bopp-Podolsky theory, the electric and magnetic fields of electrostatic and magnetostatic dipoles are still singular and their self-energy is also infinite, as will be shown in this paper.
Our purpose is to investigate the theory of second gradient electromagnetostatics which is the static version of the theory of second gradient electrodynamics [18]. We will study the textbook examples of electric point charge, electrostatic dipole and magnetostatic dipole in the framework of generalized electrodynamics, and show that second gradient electromagnetostatics yields nonsingular dipole fields and gives a straightforward regularization of the dipole singularities based on higher-order partial differential equations.
In general, in generalized electrodynamics, the electromagnetic fields (electric and magnetic potential, electric and magnetic field strengths) should satisfy the following conditions: • the field must be finite at r = 0, • the field must be everywhere continuous, • the self-energy of the field must be finite.
As mentioned, not all conditions can be satisfied for an electric point charge, an electrostatic dipole and a magnetostatic dipole using the Bopp-Podolsky electrodynamics.
Nowadays, generalized continuum theories and in particular gradient continuum theories are very popular in physics, applied mathematics, material science and engineering science. Gradient continuum theories are continuum theories which might possess characteristic length scales and characteristic time scales in order to describe size effects and memory effects, respectively. In particular, gradient theories are continuum theories valid at small scales unlike classical continuum theories like Maxwell electrodynamics. Because classical continuum theories are not valid at small scales, they lead to unphysical singularities at such scales. Thus, we are forced to regularize at short distances the classical continuum theories by means of generalized continuum theories. Gradient continuum theories provide nonsingular solutions of the field equations and a regularization of classical singularities is achieved. In physics, the most popular gradient continuum theory is the Bopp-Podolsky theory [8,9], which is the first-order gradient version of the theory of electrodynamics as mentioned above. In engineering science, a very popular gradient continuum theory is Mindlin's theory of first strain gradient elasticity [27] . An advantage of gradient elasticity theory is that it can be connected with atomistic theories and all material parameters including the appearing length scale parameters can be determined from ab initio calculations and using atomistic potentials (see, e.g., [28,29,30]). Exciting gradient effects, which are important for applications in material science, exist due to the coupling between gradient elasticity and gradient electricity in gradient electroelasticity [31], like flexoelectricity in solids [32] which is the property of a dielectric material whereby it exhibits a spontaneous electrical polarization induced by an elastic strain gradient. Furthermore, Mindlin [33] introduced the theory of second strain gradient elasticity (see also [34,35,36]). Mindlin's theory of second strain gradient elasticity involves additional material constants, in addition to the elastic constants, which can be determined from atomistic potentials (see [37]). Second strain gradient elasticity provides a better modelling of atomistics than first strain gradient elasticity. Using a simplified version of second gradient elasticity, it was possible to obtain nonsingular solutions for the elastic field produced by point defects which are elastic dipoles in solids [38]. The theory of second gradient electrodynamics has been recently proposed by Lazar [18]. It turns out that second gradient electrodynamics provides a better mathematical modelling of electromagnetic fields at small distances than the Bopp-Podolsky electrodynamics (first gradient electrodynamics). In this paper, we study the static version of it called second gradient electromagnetostatics. Of course, the coupling between second strain gradient elasticity and second gradient electromagnetostatics may lead to many interesting gradient effects of higher order which will be worth to study more in detail in future work. Therefore, gradient continuum theories are very exciting research areas of physics on small scales.
While structurally, as mathematical theories, gradient electrodynamics and gradient elasticity of n-th order are analogous, their physical significance differs slightly. Both can serve the purpose of regularization at small scales, but while gradient elasticity can be interpreted to describe microstructure and can, for example, also be derived as an approximation of lattice theories (see, e.g., [33,39]), in gradient electrodynamics analogous interpretations are not as clear. Moreover, the length scale for gradient elasticity is of the order of 10 −10 m and thus, as mentioned above, can be compared with atomistic simulations, however, in gradient electrodynamics the smallness of the length scale parameter has so far eluded experimental verification. While possible approaches have been suggested (e.g. [16]), so far none has reached the scale of 10 −15 m or smaller for the Bopp-Podolsky parameter and comparisons with quantum mechanical effects have yielded upper estimates for the length scale parameter (e.g. [17,3]). However, gradient electrodynamics remains an interesting subject, as candidate for a consistent classical field theory of electrodynamics including point charges (see also [40]), as candidate for a generalized quantum electrodynamics (see, e.g., [41,42]) and in comparison with other mathematical techniques of regularization.
The outline of this paper is as follows. In Section 2, the theory of second gradient electromagnetostatics is presented. In Section 3, we give the collection of all relevant Green functions and their derivatives. In Section 4, the nonsingular electromagnetic fields of a point charge, an electrostatic dipole and a magnetostatic dipole are computed in the framework of second gradient electromagnetostatics. The limit of those electromagnetic fields to the Bopp-Podolsky theory and to the classical Maxwell theory are given in Section 5 and Section 6, respectively. The conclusions are given in Section 7.

Second gradient electromagnetostatics
In this Section, we provide the theoretical framework of second gradient electromagnetostatics. Second gradient electromagnetostatics is the static version of second gradient electrodynamics 1 given in [18]. In the theory of second gradient electromagnetostatics, the electrostatic and magnetostatic fields are described by the Lagrangian density with the notation ∇∇E Here φ is the electrostatic scalar potential, A is the magnetostatic vector potential, E is 1 For details of second gradient electrodynamics we refer to [18].
the electrostatic field strength vector, B is the magnetostatic field strength vector, ρ is the electric charge density, and J is the electric current density vector. ε 0 is the electric constant and µ 0 is the magnetic constant (also called permittivity of vacuum and permeability of vacuum, respectively). Moreover, 1 and 2 are the two (positive) characteristic length scale parameters in second gradient electrodynamics and ∇ is the vector operator Del (or Nabla). In addition to the classical terms, first and second spatial derivatives of the (static) electromagnetic field strengths (E, B) multiplied by the characteristic lengths 1 and 2 appear in Eq. (1) describing a weak nonlocality in space. While in classical electrodynamics and electromagnetostatics the requirements of isotropy and gauge invariance lead to a unique choice for the Lagrangian, we here only have uniqueness up to null-Lagrangians (cf. [43]). Bopp [8] and Podolsky [9] introduced first-order gradient electrodynamics using different Lagrangians both being equal up to null-Lagrangians and leading to identical field equations. Our choice of Lagrangian is closer to Bopp's convention, using contractions of field gradients rather than divergences. Of course, with the introduction of higher-order terms the number of possible null-Lagrangians increases.
Note that, as is the case in classical electromagnetostatics with linear constitutive relations, the Lagrangian (1) is a sum of two purely electrostatic or magnetostatic terms. In consequence, unlike in Born-Infeld electromagnetostatics [44], it is obvious that electrostatics and magnetostatics are separated: electric currents do not produce electric fields and electric charges do not produce magnetic fields. Also note that the two energy densities are positive definite, which results in positive definite energy functionals and thus in well-posed variational problems. While in the second gradient term the positive sign is both necessary and sufficient for positivity of the energy functional, in the first term it is only sufficient. As long as 4 1 < 4 4 2 a negative sign could be allowed from the mathematical point of view, however, as will be seen below (case (3)), this would be unphysical. While the static theory works formally with this choice of parameters, the dynamic generalization contains serious problems. Also, the Lagrangian with a negative sign in the first gradient term would not be a generalization of the Bopp-Podolsky theory where the positive sign is mandatory.
In electromagnetostatics, the electromagnetic field strengths (E, B) can be expressed in terms of the static electromagnetic potentials (φ, A) because they satisfy the two electromagnetostatic Bianchi identities which are known as the homogeneous Maxwell equations. Eq. (4) states that the electrostatic field E is irrotational, and Eq. (5) states that the magnetostatic field B has no scalar sources. The Euler-Lagrange equations of the Lagrangian (1) with respect to the scalar poten-tial φ and the vector potential A give the electromagnetic field equations respectively, and the scalar differential operator of fourth order is given by where ∆ is the Laplacian. Eqs. (6) and (7) are the generalized inhomogeneous Maxwell equations in second gradient electromagnetostatics which are partial differential equations of fifth order. Eq. (6) represents the generalized Gauss law, and Eq. (7) represents the generalized Ampère law. The electric current density vector fulfills the equation of continuity If we use the variational derivative with respect to the electromagnetic fields (E, B), then we obtain the (isotropic) constitutive relations in second gradient electromagnetostatics for the response quantities (D, H) in vacuum where D is the electric excitation vector and H is the magnetic excitation vector. Therefore, in second gradient electromagnetostatics the (isotropic) constitutive relations (10) and (11) involve an invariant scalar constitutive operator of fourth order, L(∆), in addition to the scalar constitutive parameters ε 0 and 1 µ 0 . The constitutive operator L(∆) is the only linear scalar isotropic operator of fourth order, a fact that is related to the uniqueness up to null-Lagrangians of the Lagrangian for the theory. Constitutive operators of this form already showed up in second strain gradient elasticity (e.g. [33,35,38]). The higher-order terms in Eqs. (10) and (11) describe the polarization of the vacuum present in second gradient electrodynamics (see, e.g., [45,18]). Using the constitutive relations (10) and (11), the Euler-Lagrange equations (6) and (7) can be rewritten in the form of inhomogeneous Maxwell equations From Eqs. (6) and (7), the following inhomogeneous partial differential equations, being partial differential equations of sixth order, can be derived for the static electromagnetic field strengths Using the generalized Coulomb gauge condition 2 (see [46,47,4]) the electromagnetic gauge potentials fulfill the following inhomogeneous partial differential equations of sixth order The differential operator of fourth order (8) can be written in the form as product of two Helmholtz operators with two length scale parameters a 1 and a 2 , which is called bi-Helmholtz operator, and The two length scales a 1 and a 2 may be real or complex. In the theory of second gradient electromagnetostatics, the condition for the character, real or complex, of the two lengths a 1 and a 2 is the condition for the discriminant in Eq. (22), 1 − 4 4 2 / 4 1 , to be positive or negative. Depending on the character of the two length scales a 1 and a 2 one can distinguish between the following cases: (1) 4 1 > 4 4 2 : The length scales a 1 and a 2 are real and distinct and they read with a 1 > a 2 . The limit to the Bopp-Podolsky theory is given by 4 The length scales a 1 and a 2 are real and equal There is no limit to the Bopp-Podolsky theory. This case can lead to Green functions having a time dependence that increases or decreases slowly, which can give rise to unphysical results (e.g. [48,18]).
(3) 4 1 < 4 4 2 : The two length scales a 1 and a 2 are complex conjugate with There is no limit to the Bopp-Podolsky theory. For generalized electrodynamics, this case leads to Green functions having a time dependence that increases exponentially, an acausal propagation and complex mass terms (e.g. [48,49,18]). The dispersion relations of the vacuum, analogous to those computed in [45], have complex coefficients, suggesting instabilities or dissipation in the vacuum. The possible negative sign in the first gradient term of the Lagrangian mentioned above also yields complex a 1 and a 2 and thus has similar consequences.
Therefore, the case (1) is the physical one and is the generalization of the Bopp-Podolsky theory (first gradient electromagnetostatics) towards second gradient electromagnetostatics.

Green functions in second gradient electromagnetostatics
Second gradient electromagnetostatics is a linear theory with partial differential equations of sixth order, and the method of Green functions (fundamental solutions) can be used. The Green function G L∆ of the sixth order differential operator L(∆) ∆ is defined by where R = r − r and δ is the Dirac delta-function. The partial differential equation of sixth order (27) can be written as an equivalent system of partial differential equations of lower order or alternatively where G ∆ is the Green function of the Laplace operator (29) and G L is the Green function of the bi-Helmholtz operator (31).
Using partial fraction decomposition, the inverse differential operators L(∆) −1 and L(∆)∆ −1 with Eq. (19) read in the formal operator notation (see also [50]) and This formal notation directly translates into relations for the Green functions, so that the Green function G L can be written as a linear combination of two Green functions G H (a 1 ) and G H (a 2 ) corresponding to the two length scale parameters a 1 and a 2 and Helmholtz Similarly, the Green function G L∆ can be written as a linear combination of the Green function G ∆ of the Laplace operator and the two Green functions G H (a 1 ) and G H (a 2 ) corresponding to the two length scale parameters a 1 and a 2 Note that the foregoing is essentially an application of proposition 1.4.4 in [51] and could analogously be applied in linear gradient theories of any order. Using Eq. (34), the Green function of the bi-Helmholtz equation might be derived from the Green function of the Helmholtz equation 3 (see, e.g., [52,53,50]). Therefore, the bi-Helmholtz field is a superposition of two Helmholtz fields with length scales a 1 and a 2 . Using Eq. (35), the Green function of the bi-Helmholtz-Laplace equation might be derived by using the expressions of the Green function of the Laplace operator (see, e.g., [55,56,50]) and the Green function of the Helmholtz operator (see, e.g., [52,53,50]). Therefore, the bi-Helmholtz-Laplace field is a superposition of the Laplace field and two Helmholtz fields. On the other hand, resulting from the decomposition into the systems Eqs. (28) and (29), or (30) and (31), the Green function of the bi-Helmholtz-Laplace operator can be written as the convolution of the Green function of the Laplace operator and the Green function of the bi-Helmholtz equation Here, the symbol * denotes the spatial convolution. Therefore, the Green function G L plays the role of the regularization function in second gradient electromagnetostatics. Moreover, the Green function of the bi-Helmholtz equation can be written as convolution of the Green functions of the two Helmholtz operators satisfying Eqs. (31) and (19).

Green functions
The (three-dimensional) Green functions (or fundamental solutions) of the Laplace operator (29), the Helmholtz operator with length parameter a 1 , the bi-Helmholtz operator (31) and the bi-Helmholtz-Laplace operator are given by Eq. (40) is obtained by substituting Eq. (39) into Eq. (34), and Eq. (41) is obtained by substituting Eqs. (38) and (39) into Eq. (35). Moreover, the Green function (41) may be written as with the auxiliary function The series expansion (near field) of the auxiliary function (43) reads as Therefore, the function f 0 (R, a 1 , a 2 ) regularizes up to a 1/R-singularity towards a nonsingular field expression. Indeed, the Green function (41) is nonsingular and finite at R = 0, namely .
On the other hand, the Green function (40) is nonsingular and possesses a maximum value at R = 0, namely (see Fig. 1) .
Moreover, the Green function (40) is a Dirac-delta sequence with parametric dependence a 1 and a 2 lim with lim a 2 →0 where the limit is to be understood in the weak sense for distributions. The Green functions (39) and (40) are plotted in Fig. 1.

Derivatives of the Green function G L∆
Now, we calculate the first, second and third gradients of the Green function G L∆ . The first, second and third gradients of the Green function (41) are obtained as 4 ∇∇G L∆ (R) = 1 4π 1 and (cf. Eq. (30)) ∆G L∆ (R) = 3 4π 1 with the auxiliary functions 4 The expression sym(1R) means δ ij R k + δ jk R i + δ ki R j . Note that the auxiliary functions for the gradients of a Green function in the form (42) obey for i = 0, 1, 2. 5 This relation has direct consequences for the series expansions (the f i are analytic and can be differentiated term by term): f 1 has no linear term, f 2 neither linear nor cubic, f 3 no R, R 3 and R 5 -term, and so forth. Therefore, the first non-vanishing term of even order in f 0 , which here is the fourth-order term, determines the strength of the regularization.
The relevant series expansions of the auxiliary functions (54)-(56) (near fields) read as In Eqs. (58)-(60), it can be seen that the function f 1 (R, a 1 , a 2 ) regularizes up to a 1/R 3 -singularity and the functions f 2 (R, a 1 , a 2 ) and f 3 (R, a 1 , a 2 ) regularize up to a 1/R 4singularity towards nonsingular field expressions. The auxiliary functions (43), (54)-(56) are plotted in Fig. 2. In the far field, the auxiliary functions (43), (54)-(56) approach the value of 1 and in the near field, they are modified due to gradient parts and approach the value of 0 at the position R = 0 (see Fig. 2).

Electromagnetic fields in second gradient electromagnetostatics
The electromagnetic potentials are the solutions of the inhomogeneous partial differential equations of sixth order (17) and (18) for given charge and current densities (ρ, J )

Electric point charge
The charge density of an electric point charge located at the position r is given by where q denotes the electric charge. This means that r is the position vector of the point charge and r is the field vector. Substituting Eq. (63) into Eq. (61) and performing the convolution, the electrostatic potential of a point charge reads as If we insert the Green function (42) into Eq. (64), the explicit expression of the electrostatic potential of a point charge reads in terms of the auxiliary function (43) Using the near field of f 0 , Eq. (44), it can be seen that the electrostatic potential (65) is finite at R = 0, namely (see Fig. 3a) .
The electric field strength (2) of a point charge is given by the negative gradient of the electrostatic potential (64) and reads as Using Eq. (49), Eq. (67) reduces to Using the near field of f 1 , Eq. (58), it can be seen that the electric field (68) is zero at R = 0. In general, it is nonsingular and possesses an extremum value near the origin (see Fig. 3b). It does not have a directional discontinuity at the origin unlike the electric field  The electric field energy can be written as where we have used integration by parts and the surface terms vanish at infinity. Substituting Eq. (63) into Eq. (69), the interaction energy of two point charges reads which is finite in the whole space. The electrostatic self-energy of a point charge is obtained as which is finite. The electrostatic part of the Lorentz force reads as [18] Substituting Eqs. (63) and (68) into Eq. (72), the electrostatic interaction force between two point charges q at r and q at r is obtained as which is zero at R = 0 and nonsingular. Eq. (73) is the force exerted by one charge q at r on the other charge q at r. It holds F qq = −F q q .

Electric dipole
The charge density of an ideal electric dipole is given by where p is the electric dipole moment. Substituting Eq. (74) into Eq. (61) and employing the convolution, we obtain for the electrostatic potential of an electric dipole, or the electric dipole potential, Now, using Eq. (49), Eq. (75) becomes (R, a 1 , a 2 ) . (76) It is zero at R = 0 and possesses an extremum value near the origin (see Fig. 3b). The electric field strength (2) of an electric dipole is given by the negative gradient of the electric dipole potential (75) and, using Eq (50), it reads as Moreover, using the near fields of f 1 and f 2 , Eqs. (58) and (59), it can be seen that the electric dipole field (78) is finite at R = 0 due to the f 1 -term, namely (see Figs. 3c and 3d) and it does not have a directional discontinuity. See Fig. 5 for a comparison of dipole potentials and fields in second gradient electromagnetostatics and the limits to Bopp-Podolsky theory and classical Maxwell theory. Along with the regularization directional changes in the electric field of the dipole in the gradient theory are introduced, due to the appearance of two different auxiliary functions in the two terms of the dipole field (78). These modifications are stronger in second gradient electromagnetostatics than in the Bopp-Podolsky theory and even lead to qualitatively new behavior: There exist two zeros of the dipole field. Substituting Eq. (74) into Eq. (69), the interaction energy of two electric dipoles p and p , a distance R = r − r apart, becomes The electrostatic self-energy of an electric dipole reduces to which is finite.
Substituting Eqs. (74) and (77) into Eq. (72) and using Eq. (51), the electrostatic interaction force between two electric dipoles reduces to which is finite at R = 0 and nonsingular (see Figs. 3e and 3f), but it has a directional discontinuity at R = 0. Eq. (73) is the force exerted by one electric dipole p at r on the other electric dipole p at r. It holds F pp = −F p p .

Magnetic dipole
The electric current density vector of a magnetic dipole is given by where m is the magnetic dipole moment. Substituting Eq. (83) into Eq. (62) and performing the convolution, we obtain for the magnetic vector potential of a magnetic dipole, or the magnetic dipole potential, Using Eq. (49), Eq. (84) becomes which is zero at R = 0 and possesses an extremum value (see Fig. 3b). The magnetic field strength (3) of a magnetic dipole is given by the curl of the magnetic dipole potential (84) and, using Eq. (50), it becomes Using the near fields of f 1 and f 2 , Eqs. (58) and (59), it can be seen that the magnetic dipole field (87) is finite at R = 0 due to the f 1 -term, namely (see Figs. 3c and 3d) and it does not have a directional discontinuity. A comparison of vector potentials and magnetic fields in second gradient electromagnetostatics, the Bopp-Podolsky theory and classical electromagnetostatics is given in Fig. 6. Again, similarly to the electric dipole fields, the regularization modifies the field direction in the near field. The magnetic field energy can be written as where we have used integration by parts and the surface terms vanish at infinity. Substituting Eq. (83) into Eq. (89), the interaction energy of two magnetic dipoles m and m , a distance R = r − r apart, is given by The magnetostatic self-energy of a magnetic dipole reduces to which is finite. The magnetostatic part of the Lorentz force reads as [18] Substituting Eqs. (83) and (86) into Eq. (92) and using Eqs. (51) and (53), the magneto- static interaction force between two magnetic dipoles becomes which is finite at R = 0 and nonsingular (see Figs. 3e and 3f), but it possesses a directional discontinuity at R = 0.

Some qualitative side-effects of regularization
As we have seen, the gradient theory for large distances approaches the classical results asymptotically, while on small length scales Green functions and therefore electric and magnetic fields are modified. These modifications introduce new qualitative behavior that the reader should be aware of. For example, apart from weakening or removal of singularities, the field of the electric dipole is reversed in the direct vicinity of the dipole and two singular points appear, zeros of the electric field, where a point charge would be in equilibrium (see Fig. 5f). These singular points also exist in the fields of real dipoles, i.e., two point charges a certain distance apart, and thus motivate the following.
Consider three point charges on a straight line, for simplicity we take two charges q at the same distance r to a charge q in the middle. Due to the symmetry of the situation the force on the charge q is zero, and through Eq. (73) the force on the other two charges is For q q > 0 this expression is positive, the force is therefore repulsive for all r, while for q q < 0 we find asymptotically for large r and for small r. Now, for a charge ratio with 1 2 |q | < |q| < 4|q |, we can have F negative at large distances and positive at sufficiently small distances, by continuity we thus find a zero. This means there exists a situation where three point charges at rest can be in equilibrium just through the electromagnetic forces in vacuum in second gradient electromagnetostatics. This situation has no analogue in classical or Bopp-Podolsky electromagnetostatics. A linear stability analysis shows that the state is stable with respect to small perturbations in r.

Reinterpretation of second gradient electromagnetostatics as electromagnetostatics with extended charge and current densities
Now, we demonstrate that the solutions, Eqs. (61) where we have introduced the following extended charge and current densities as convolution of the Green function G L and the singular, classical charge and current densities ρ and J . Because the Green function G L of the bi-Helmholtz operator is finite, also the extended charge and current densities are finite unlike the Bopp-Podolsky case with singular Green function G H of the Helmholtz operator (see Fig. 1). The interpretation of extended, singular charge distribution was proposed in the Bopp-Podolsky theory in [15,5]. From the physical point of view, the Green function G L , Eq. (40), plays the role of a "form factor" in Eqs. (99) and (100) since it characterizes the shape of the charge and the current. So, the electromagnetostatic potentials (97) and (98) satisfy the following Poisson equations with the extended charge and current densities (99) and (100) as sources. Moreover, the convolution of the Euler-Lagrange equations (6) and (7) with the Green function G L and the use of Eqs. (31), (99) and (100) lead to "classical" inhomogeneous Maxwell equations with the extended charge and current densities (99) and (100) as sources. Of course, the extended charge and current densities (99) and (100) satisfy inhomogeneous bi-Helmholtz equations This proves that Eqs. (103) and (104) are equivalent to Eqs. (6) and (7).

Electromagnetic fields in first gradient electromagnetostatics (Bopp-Podolsky electromagnetostatics)
Now we perform the limit from second gradient electromagnetostatics to the Bopp-Podolsky electromagnetostatics for the electromagnetic fields of a point charge, an electric dipole and a magnetic dipole. In the limit a 2 → 0 and a 1 → (Bopp-Podolsky limit), the auxiliary functions (43), (54)-(56) simplify to The auxiliary functions (107)-(110) are plotted in Fig. 7. The relevant series expansion (near field behavior) of the auxiliary functions (108)-(110) (near fields) reads In Eq. (111), it can be seen that the function f 0 (R, ) regularizes up to a 1/R-singularity and gives a nonsingular field expression. In Eqs. (112)-(114), it can be seen that the functions f 1 (R, ), f 2 (R, ) and f 3 (R, ) regularize up to a 1/R 2 -singularity and give nonsingular fields. At R = 0 the auxiliary functions (107)-(110) are zero and in the far field they approach the value of 1.

Electric point charge
In the Bopp-Podolsky limit, the electrostatic potential (65) reduces to which is finite at R = 0, namely (see Fig. 3a) Of course, Eq. (115) is in agreement with the original expression given by Bopp [8] and Podolsky [9]. Moreover, the electric field strength (68) reduces to which is in agreement with the expression given by Bopp [8] (see also [26]). The electric field strength of a point charge is nonsingular but has a directional discontinuity at the origin, namely whereR = R/R. Therefore, the projection of the electric field (117) onto a curve passing trough the location of the charge jumps from q/(8πε 0 2 ) to −q/(8πε 0 2 ) at R = 0. The self-energy (71) of a point charge becomes which is positive and finite for > 0.

Electric dipole
In the Bopp-Podolsky limit, the electrostatic potential (76) reduces to which is in agreement with the expression given by Bonin et al. [4]. The electrostatic potential of an electric dipole is finite and possesses a directional discontinuity at the origin, namely The electrostatic field strength (78) becomes which is singular because it has a 1/R-singularity at the origin. Therefore, the selfenergy (81) of an electric dipole and the interaction force (82) between two electric dipoles become infinite and singular in the Bopp-Podolsky theory.

Magnetic dipole
In the Bopp-Podolsky limit, the magnetic vector potential (85) reduces to which is in agreement with the expression given by Bonin et al. [4]. The magnetostatic potential of a magnetic dipole is finite and possesses a directional discontinuity at the origin The magnetic field strength (87) becomes which is singular since it possesses a 1/R-singularity at the origin. Therefore, the selfenergy (89) of a magnetic dipole and the interaction force (93) between two magnetic dipoles become infinite and singular in the Bopp-Podolsky theory. Note that Eq. (125) is in agreement with the expression given by Bonin et al. [4].

Electromagnetic fields in classical Maxwell electromagnetostatics
The limit from the Bopp-Podolsky electromagnetostatics to the classical Maxwell electromagnetostatics is → 0.

Electric point charge
respectively, are recovered.

Electric dipole
From Eqs. (120) and (122), the classical electrostatic potential and the electric field strength of an electric dipole are obtained as and respectively. Note that Eqs. (128) and (129) are in agreement with the expressions 6 given in the literature [19,20,21,22].

Magnetic dipole
From Eqs. (123) and (125), the classical magnetic vector potential and the magnetic field strength of a magnetic dipole are obtained as and respectively. Note that Eqs. (130) and (131) are in agreement with the expressions given in the literature [19,20,21,22]. 6 In the classical Maxwell electromagnetostatics of dipoles, the term proportional to the Dirac delta function follows from the second order derivatives of 1/R in the sense of generalized functions [52] and is known as the Frahm formula [20]

Conclusion
We have presented second gradient electromagnetostatics, which is the static version of second gradient electrodynamics, a generalization of Bopp-Podolsky electrodynamics, as singularity-free field theory of electromagnetic fields, in particular, for an electrostatic point charge, an electrostatic dipole and a magnetostatic dipole. Through linear field equations of sixth order, second gradient electromagnetostatics yields an even stronger regularization than the Bopp-Podolsky theory: the field of a point charge is nonsingular and zero at the origin, the fields of electric and magnetic dipoles are nonsingular as well and dipoles have a finite self-energy, which diverges in the limit to Bopp-Podolsky electrodynamics (see Table 1). The regularization stems from the fact that the Green function of the bi-Helmholtz-Laplace operator appearing in second gradient electromagnetostatics and its first, second, and third gradients are singularity-free: In every single problem studied, we were able to recover the Bopp-Podolsky results and the Maxwell results in the proper limits, when the length scale parameter 2 (or a 2 ) goes to zero and the Bopp-Podolsky length scale parameter goes to zero, respectively. We have also demonstrated that, much like in Bopp-Podolsky electrodynamics, the above results can be obtained through a special ansatz for extended charge distributions in classical electrostatics instead of an invariant generalized field theory. This is due to the structure of the Green function, a convolution of multiple Green functions to second-order operators, some of which can be attributed to a charge distribution in this interpretation.
Another interesting detail in second gradient electromagnetostatics is the existence of stable equilibria of three point charges. While in classical electrodynamics the maximum principle for the Poisson equation forbids such states (this is sometimes referred to as Earnshaw's theorem), the higher-order generalized analogue to the Poisson equation has no maximum principle. The same is true for the fourth-order equation for the potential in Bopp-Podolsky electrostatics, although a stable configuration of point charges might be harder to construct there.
Finally, we conclude that the theory of second gradient electromagnetostatics provides a singularity-free generalized continuum theory of electromagnetostatics with generalized Gauss law and generalized Ampère law valid down to short distances. The covariant form of second gradient electrodynamics and its meaning as a nonsingular relativistic field theory will be given in a forthcoming publication.