The electromagnetic potentials due to charge and current distributions
2026-05-21
What we want to show
To begin this discussion, we state the electromagnetic wave equation for the scalar () and vector () potentials sourced by some time-dependent charge () and current density, respectively. In the Lorenz gauge, the sources and fields are related by the separate wave equations
Supposing that we know the electromagnetic potentials at every point in space, Eq. allows us to compute the charge and current densities that generated them. In the interest of practicality, our goal here is to invert this expression: we want to be able to compute the electromagnetic potentials generated by a given charge and current density. We will show that these relationships are given by
which requires a surprisingly lengthy derivation. The symbol is the retarded time, which indicates that the change in potential due to a disturbance in charge density (or current density) at position is only apparent at a position once the field disturbance propagates the intervening distance at the speed of light.
The impulse response (Green function)
The wave equations in Eq. take the general form
Linear operators for the electromagnetic potentials
The exact expressions for the linear operators are
which acts on the scalar potential by and
which acts on each component of the vector potential for . Since the equations differ at most by a pre-factor, it is easy to swap out this pre-factor to find the other components once one component is solved.
in which the linear operator acting on some profile (in our case, the potentials) results in a quantity equal to the source term (in our case, the charge and current densities). While the form specialized to the case of electromagnetic potentials is accessible in the box above, we will choose to analyze the general form to emphasize that the process of inverting the linear operator to yield from can be done for an arbitrary linear system and that there is nothing special about electromagnetic fields in this sense.
The desired inversion of Eq. is achieved by computing the Green function, which is just a fancy name for the spatiotemporal impulse response. That is to say that the Green function is the field at position and time generated by some instantaneous excitation at position and time . We denote this function as , indicating that is directly a function of and while the parameters and detail the position and time of the excitation. To be more specific, the instantaneous excitations are modelled by delta functions so that the governing equation that defines the Green function is
We can pause here and ask how Eq. is useful to us. The answer lies in the fact that any arbitrary source can be expanded in the delta function basis by
This integral is trivially true, but it is convenient because it allows us to act on each side of Eq. by the integral operator for some arbitrary source distribution . The end result is
Derivation
This is straightforward, but still worth doing carefully. Multiply and integrate:
On the left-hand side, we can pull out of the integral because it does not depend on the integrated (primed) coordinates. The right-hand side is just the integral we introduced in Eq.. Putting it all together, we find that
and all that’s left is to apply our convenient parentheses to make this expression easy to compare to Eq..
Finally comparing Eq. to Eq. , we find that the Green function can be applied to the source term by
to compute the resulting field .
Specialize: The electromagnetic Green function
Now we specialize the Green function analysis to the case of electromagnetic fields. Substituting in the linear operator introduced by the wave equation, the problem we must solve is
where for the scalar potential and for components of the vector potential. Solve Eq. by expanding both sides in the plane-wave (Fourier) basis so that the differential operators reduce to momenta and energies . Doing so, we find a new expression for the Green function
Derivation
To begin, we define as the 3+1-dimensional Fourier transform of the Green function:
Note that the primed coordinate parameters remain spatial because they indicate that we are interested in the response function (in the Fourier domain or otherwise) due to a disturbance at position and time . The corresponding inverse transformation is
which we can substitute into Eq. to find
To simplify the right-hand side, use the fact that
and use the delta function identities to substitute
As a result, we find the right-hand side to be
For the last inequality in Eq., we separated the complex exponentials for primed and unprimed coordinates, allowing us to move everything in Eq. to one side:
It follows from the orthonormality of the plane waves that the term in brackets must equal zero. This equality is solved when the transformed Green function is
Finally, apply the inverse Fourier transformation defined in Eq. to the Green function in Eq. to get the real-space solution
Note that dependence on spatial and temporal coordinates only appear in the differences and . To simplify the notation, we can then define and so that the Green function has the simple form
where and . It then suffices to solve for the Green function in this shifted coordinate system, after which we can just re-insert and (as well as and ) to find a clear expression for .
We can evaluate the integral in Eq. by converting to spherical coordinates. Doing so, it is straightforward to integrate over the azimuthal and polar angles, leaving only an integral over the radial coordinate :
which includes the contribution found by integrating over the angular frequency :
Derivation
Starting with the integral for the Green function in Eq.,
transform into spherical coordinates by the substitution
Since the integral will be taken over all angles and , we may choose a coordinate set in which is parallel to for a given arbitrary input . This choice simplifies the dot product in the complex exponential; in the definition , the angle between vectors and becomes simply the polar angle so that . Finally, insert the defintion of the differential volume element to find
The integral over azimuthal angle in Eq. is trivial because the integrand lacks -dependence. We integrate over the polar angle by substituting and (and accordingly, ). This integral evaluates to
leaving the remaining integration over and ,
The remaining difficulty is in the integration over , which we call . Factoring the denominator , it is clear that there are singularities in the integrand when . Since the integrand is undefined at these points, the best we can do is the principal value integration
which removes an infinitesimally small symmetric interval about each pole so the result stays finite.
A note on the principal value integration
The first thing to point out is that integration is done in the principal value (P.V.) sense, which means that we remove a symmetric region of size centered around each pole .
This must be done to remove the singularities so that the integral remains possible to evaluate. Finally, we take the limit as to remove the error introduced by the reduction in the integral limits. The quick justification for why we do this is that we are interested in this quantity only in the distributional sense (as the limit of a sequence of functions to be integrated against a test function). In this case, we are clearly following by an integral over so that the test function is simply the integrand of Eq..
Most of the remaining effort will be to evaluate the integral in Eq.. We see in the next section that the integral can be calculated by analytic continuation of the integrated variable (and corresponding integrand) into the complex domain and performing a complex line integral whose path includes the integral over the real line that we aim to find.
Principal value integration in the complex plane
We begin to evaluate the integral in Eq. by explicitly imposing the limits dictated by the principal value integration procedure: remove a symmetric interval of length about the poles at . In the sense of distributions, we may then define as the limit of a sequence of functions defined as the integral over these limits as the excluded interval . Restated in mathematical notation (and more clearly),
While this is an integral over the real line, we can evaluate it easily by taking the analytic continuation
and performing a closed line integral in the complex plane. The figure below shows how we can define a closed contour that includes the real-line principal value integration (the red line) and semicircular arcs which avoid the singularities so that is analytic over the path.
Figures 1a) and b) close the contour in the upper and lower-half plane, respectively. As of now, we do not know which direction is corrent. But in either case, we make the choice to exclude the poles from within the contour so that the line integral is zero by the residue theorem. The real-line integral we want is then found from
Derivation
The contours shown in the Figure contain the limits
Though the ordering is swapped, the path traverses the real line from to , a small semicircle about the pole at (), a path along the real line from to , a small semicircle about the pole at (), the real line from to , and a semicircular arc with infinite radius connecting from to ().
If we define the excursions around each singularity so that contains no poles within its boundaries, then is analytic everywhere within the contour so the integral by the residue theorem. Eq. then rearranges to
That is, we can calculate the desired integral by evaluating only the semicircular contours about each pole along with the infinite radius semicircle that connects the positive and negative real axes.
where the notation denotes a semicircular arc of radius at position . We must now impose the physical constraint that , which corresponds to the retarded propagator. In the analysis of the electromagnetic wave equation, this means that the existence of a charge at time is only seen in the electromagnetic fields at later times , and the disturbance does not travel backward in time. However, it is a quirk of the laws of physics that we can also fix and still get a perfectly good mathematical solution. For the retarded propagator, we determine that the integral must follow subfigure a) (closure in the lower half plane). We can further show that the semicircular arc over as is zero, so that we can find strictly from integration over the semicircular contours.
Integration over infinite-radius arc
The constraints on the integration path arise from integration along the arc of radius over an angle . To determine whether this arc should traverse through the upper- or lower-half plane, we need to observe the asymptotic behavior of the numerator (because we know that the contribution from the denominator tends to zero as both the real and imaginary parts of go to zero). Parametrizing the integrand in polar coordinates, the numerator along the arc is
When , must be negative because if , then the term would diverge as . This quantity is negative only when the angle , which correspond to the lower half plane (for strict equality, the sine is zero so the points at and are well-behaved). Consequently, we must choose a contour that closes in the lower half plane where as . In fact, the integral over the arc can be bounded by the triangle inequality to find
Taking the limit as , we find that the exponential dominates so that . In the end, the infinite-radius arch contributes nothing to the contour integral.
Complaint about the way these propagators are commonly calculated
Choosing contours that exclude poles from the integration region is the most convenient choice, but it is not required. We would find the same end result for if we included the poles inside the contour. If (retarded propagator), we always need to close the contour in the lower half plane. Then, included the poles in this contour, the residue theorem tells us that the evaluation of the integral about the entire arc is just the residue from the two poles (instead of zero for the pole-free contour). But now the semicircular arc must go about each pole in the opposite direction compared to the pole-free integration, so half the effect of each pole is negated to give the same result as just excluding them in the first place.
The way we choose to perform an integral does not change its physical interpretation. Though, for example, if we want the solution for , the contour closure will be different so the integral is performed differently. This just makes the function defined piecewise for (retarded) and (advanced), although we don’t have much physical use for the advanced solution. I think things like the Feynman propagator are just a byproduct of the hand-wavy way that the deformation is done around the poles. To me, the idea of imbuing a method for applying math to a problem with physical meaning, much like the insistence that virtual transitions and virtual particles are “real” because we choose to do time-dependent perturbation theory, is lazy and counterproductive. And I still think that’s being charitable to this approach; in reality it’s just bad math.
We can perform the integrations for the semicircular deformations about each pole to find
Derivation
Since the contour is closed in the lower-half plane and we chose to exclude the poles from the line interior, the excursions around the poles must also extend into the lower half plane and proceed in the positive (counter-clockwise diraction). Beginning with the pole at , we evaluate the integral in the neighborhood of the singularity by defining and . Doing so, the integral (over the integrand in Eq.) becomes
The expression in parentheses is analytic on a disk of radius centered at . Since the deformation integral takes place along a semicircle with radius , the analyticity condition holds along the integration path so we can expand the parenthetical term in a power series. Defining the integrand in Eq. as , this power series expansion gives
where we collected the analytic power series terms into the function . Since is a power series, it is bounded by some maximum value along the semicircular contour over which we are integrating. The contour integral about this term is then bounded by the maximum value using the triangle inequality
Since , we conclude from the ineqality that the contour integral of is zero. What remains is the integral over the singular term of . Using the fact that (since is just the power series expansion term when ), we can evaluate the integral directly to find
To evaluate Eq., we defined (and ) to integrate over the semicircular contour in the lower-half plane.
The integral over the other pole goes exactly the same, except for a sign change due to the location of the pole, so that
Applying Eqs. and to the result in Eq., we find
which we may finally use to evaluate the integral in Eq. for the Green function. The result is that
Derivation
For convenience, restate the Green function integral from earlier:
Substituting Eq. for
Defining , we can group these integrals into ones over all values of using
so that
The remaining integrals are the familiar delta function integrals, allowing us to simplify the Green function to
It is tempting to say that Eq., containing two delta functions, includes both the retarded () and advanced () propagators. But we were only able to perform the integrations by fixing the parameter , in which case the advanced delta function condition is never satisfied (since is a magnitude). So we can safely remove this term, and we are left with
If we want the advanced operator, we just need to repeat the process under the condition , which requires closure in the lower-half plane. In the end, we end up with exactly what is suggested by Eq., which is why its possible to be sloppy and still get the “correct” answer.
We can now apply the Green function according to Eq. with source terms and for the scalar potential and th component of the vector potential, respectively. Also applying the definition of our variable (scalar potential) and (vector potential), we finally find
which is the relationship we wanted to show (Eq.).
Derivation
Now that we have the Green function, we can calculate the potentials using Eq.:
Clearly, we need to restore the usual coordinate dependence by substituting back in and to find
It is then easiest to apply the delta function during the integral over , since the spatial coordinates exist within magnitude brackets. Using the delta function identity gives
All that remains is to plug in and to find the scalar potential due to a charge distribution and and to find the vector potential components due to a current distribution.
