\documentclass[a4paper,11pt]{article}
\parskip=1.5ex
\usepackage{verbatim,a4margins,amsmath,bm,hyperref}

\newcommand{\dd}[2]{\frac{\mathrm d #1}{\mathrm d #2}}
\newcommand{\dds}[2]{\frac{\mathrm d^2 #1}{{\mathrm d #2}^2}}
\newcommand{\pdd}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\pdds}[2]{\frac{\partial^2 #1}{\partial{#2}^2}}
\newcommand{\pddd}[3]{\frac{\partial^2 #1}{\partial #2\partial #3}}

\newcommand{\rowfour}[4]{\left[\begin{array}{cccc}#1&#2&#3&#4\end{array}\right]}
\newcommand{\vectwo}[2]{\left[\begin{array}{c}#1\\#2\end{array}\right]}
\newcommand{\mxtwo}[4]{\left[\begin{array}{cc}#1&#2\\#3&#4\end{array}\right]}

\newcommand{\sq}[1]{\left(#1\right)^2}
\newcommand{\floor}[1]{\left\lfloor#1\right\rfloor}
\newcommand{\abs}[1]{\left|#1\right|}

\newcommand{\rem}[1]{}

\newcommand{\A}{\mathbf A}
\newcommand{\E}{\mathbf E}
\newcommand{\B}{\mathbf B}
\newcommand{\J}{\mathbf J}

\newcommand{\x}{\mathbf x}
\newcommand{\y}{\mathbf y}
\newcommand{\z}{\mathbf z}

\newcommand{\s}{\mathbf s}
\newcommand{\rv}{\mathbf r}

\newcommand{\e}{\mathbf e}

\newcommand{\de}{\,\mathrm d}

\begin{document}

\title{Smoothed Retarded Potentials and Fields}
\author{Stephen Brooks}
\maketitle

% https://en.wikipedia.org/wiki/Li%C3%A9nard%E2%80%93Wiechert_potential
% https://en.wikipedia.org/wiki/Retarded_potential

\section{Aims}
This note will follow the derivation of the Li\'enard--Wiechert potentials for a moving charge, which give the entire electromagnetic field generated by that charge including time delays, while concurrently deriving a version for a smoothed charge.  The smoothed charge replaces the delta function charge density $\rho\propto\delta^3(\x)$ by a spherically symmetric function $f(\abs{\x})$, so instead of a potential $\phi\propto1/\abs{\x}$ for a charge at rest, the potential will be the 3D convolution of $f$ with $1/\abs{\x}$.

The smoothing will be defined in the `laboratory frame' rather than the rest frame of the particle.  If a spherically-symmetric distribution is required in a different rest frame, for example the averaged rest frame of a whole particle beam, the particle histories may be Lorentz boosted and the resulting fields Lorentz boosted back.  Another motivation for using the lab frame is that if a beam is uniformly accelerated by an electric field, it maintains its lab frame size, rather than shrinking as a rigid object would (if the beam enters an electric field from the side, the particles maintain their time separation and the length actually \emph{increases} with the beam velocity).  Properly defining a charge distribution of fixed size in the changing reference frame of an accelerating particle is somewhat complicated and relates to the concept of Born rigidity.

\section{Notation}
The potentials and fields are being evaluated at a spacetime point $(\x,t)$ and sometimes will be expressed as integrals over other spacetime points $(\y,u)$.  The source has charge $q$ and position history $\s(t)$; for the smoothed charge this is the centre of the charge distribution in the lab frame.  The `retarded time' at which a signal leaving a point $(\y,t_r)$ would arrive at our evaluation point $(\x,t)$ is $t_r(\x,t,\y)=t-\frac1c\abs{\x-\y}$.
\subsection{Four-vectors}
The four-current density $J^\mu=(c\rho,\J)$ and four-potential $A^\mu=(\phi/c,\A)$ allow some of the formulae to be expressed more compactly.

\section{Potentials}
\subsection{Maxwell's Equations}\label{maxpot}
Maxwell's equations with sources in the Lorenz gauge may be expressed compactly in terms of the four-vectors as
\[ \partial_\mu\partial^\mu A^\nu = \mu_0J^\nu, \]
where the repeated index summation convention has been used.  A general solution with free boundary conditions in space can be written
\[ A^\mu(\x,t) = \frac{\mu_0}{4\pi} \int \frac{J^\mu(\y,t_r)}{\abs{\x-\y}}\de^3\y. \]
This is known as the retarded solution, where fields depend on the particle position in the past where they were emitted.  It assumes there were no freely-propagating waves at $t=-\infty$ and all radiation is created by the particle.  If $t_r$ is replaced by $t_a(\x,t,\y)=t+\frac1c\abs{\x-\y}$, this becomes the advanced solution, where propagating waves from the past converge on and are absorbed by the particle, so mathematically it `depends on' particle positions in the future where the fields are absorbed.

The solution may be re-expressed as a 4D integral
\[ A^\mu(\x,t) = \mu_0 \iint J^\mu(\y,u)\frac{\delta(u-t_r)}{4\pi\abs{\x-\y}}\de^3\y\de u, \]
showing that the Green's function is shaped like a past light cone as $t_r$ depends on $\y$.

\subsection{Source Definitions}
The four-current density for our moving point particle is
\[ J^\mu(\x,t) = qV^\mu(t)\delta^3(\x-\s(t)) \qquad \mbox{where} \qquad V^\mu(t)=(c,\dot\s(t)). \]
This $V^\mu$ is not the usual four-velocity $U^\mu=(\gamma c,\gamma\dot\s)$ that would appear if a charge density were being transformed to a boosted frame, in fact $V^\mu=U^\mu/\gamma$, so $V^\mu$ does not transform as a Lorentz vector.  The reason for this is that the $\delta^3$ function only cares about the integral of the charge density over a small volume, and the relativistic boost would compress volume to $1/\gamma$ while the charge density increases by $\gamma$.  The two effects cancel, leaving the total charge $q$ constant under accelerating the particle, as observed in reality.

The smoothed charge source has an analogous form
\[ J_s^\mu(\x,t) = qV^\mu(t)f(\abs{\x-\s(t)}), \]
with the delta function replaced by a spherically-symmetric distribution that should also integrate to one:
$ \int_0^\infty f(r)4\pi r^2 \de r = 1. $
Again, $V^\mu$ appears rather than $U^\mu$ because the smoothing distribution is defined in the lab frame and is not Lorentz transformed.

\subsection{Derivation for Point Particle}\label{pointpot}
Substituting the point particle source into the solution of Maxwell's equations gives
\[ A^\mu(\x,t) = \frac{\mu_0q}{4\pi} \iint V^\mu(u)\delta^3(\y-\s(u)) \frac{\delta(u-t_r)}{\abs{\x-\y}}\de^3\y\de u. \]
Noting that $\s(u)$ does not depend on $\y$, the 3D integral and delta function can be cancelled, replacing $\y$ by $\s(u)$:
\[ A^\mu(\x,t) = \frac{\mu_0q}{4\pi} \int V^\mu(u) \frac{\delta(u-t_r(\x,t,\s(u)))}{\abs{\x-\s(u)}} \de u. \]
The remaining delta function is not so simple, as its argument varies non-linearly with $u$ and the rule $\delta(f(x))=\sum_{f(x_n)=0}\frac1{\abs{f'(x_n)}}\delta(x-x_n)$ needs to be used before integrating it.
\begin{eqnarray*}
\pdd{}{u}t_r(\x,t,\s(u)) &=& \pdd{}{u}\left(t-\frac1c\abs{\x-\s(u)}\right) \\
 &=& -\frac1c \pdd{}{u}\abs{\x-\s(u)} \\
 &=& -\frac1c\frac{\s(u)-\x}{\abs{\s(u)-\x}}\cdot\dot\s(u).
\end{eqnarray*}
The delta function is active when $u=t_r(\x,t,\s(u))$, which is when the particle intersects the past light cone of $(\x,t)$.  This does not have a closed form solution but can be computed numerically given the particle's position history $\s(u)$.  From now on this solution $u$ will be referred to as $t_{r,\s}(\x,t)$; the solution is unique for slower-than-light particles.

Making the required substitution and putting the gradient of the delta function on the denominator gives:
\[ A^\mu(\x,t) = \frac{\mu_0q}{4\pi} V^\mu(t_{r,\s}) \frac{1}{\abs{\x-\s(t_{r,\s})} \abs{1+\frac1c\frac{\s(t_{r,\s})-\x}{\abs{\s(t_{r,\s})-\x}}\cdot\dot\s(t_{r,\s})} }. \]
The second term on the denominator is always positive for slower-than-light particles, so the absolute value may be removed.  Introducing $\rv=\x-\s$ and evaluating $\s$, $\dot\s$ and $\rv$ at $t_{r,\s}$ implicitly gives a recognisable form of the Li\'enard--Wiechert potential:
\[ A^\mu(\x,t) = \frac{\mu_0q}{4\pi} {\left.\frac{(c,\dot\s)^\mu}{\abs{\rv} (1-\frac{\rv}{\abs{\rv}}\cdot\frac{\dot\s}{c}) }\right|}_{t_{r,\s}(\x,t)}. \]

\subsection{Derivation for Smoothed Charge}
Substituting the smoothed charge source into the solution of Maxwell's equations gives
\begin{eqnarray*}
 A^\mu(\x,t) &=& \frac{\mu_0q}{4\pi} \iint V^\mu(u)f(\abs{\y-\s(u)}) \frac{\delta(u-t_r)}{\abs{\x-\y}}\de^3\y\de u \\
 &=& \frac{\mu_0q}{4\pi} \int V^\mu(t_r) \frac{f(\abs{\y-\s(t_r)})}{\abs{\x-\y}}\de^3\y.
\end{eqnarray*}
This time, the $\y$ integral cannot be removed because there is no 3D delta function.  There are two routes to proceed: one is to convolve the point particle solution from section \ref{pointpot} with the charge distribution $f(|\x|)$, which requires integrating $t_{r,\s}(\x,t)$ and numerically solving an equation with the arbitrary function $\s$ for every value of $\x$.  Instead, here the \emph{time} integral has been removed before $\s$ can complicate the process.
\rem
{
The integral will be evaluated using a series solution centred around the point where $\y=\s(t_r(\x,t,\y))$, which corresponds to $f(0)$, where we assume most of the charge will be.  Trying $\y=\s(t_{r,\s})$, the definition of $t_{r,\s}$ gives $\y=\s(t_r(\x,t,\s(t_{r,\s})))=\s(t_r(\x,t,\y))$, so this is the point in question.  Note that calculating $t_{r,\s}(\x,t)$ here only requires one numerical solve because $\x$ is not being integrated over.
% Work out dt_r/dy_i, then ds_j(t_r)/dy_i, then change variables to z=y-s(t_r) but using the Jacobian
}

% Another idea: change variables to t_r and two spherical angles.  Integrate across the relevant t_r range numerically.
% The spherical integral should then have a simpler analytic form and be possible via fast-converging series around the peak
This 3D integral still includes the arbitrary function $\s$ and will have to be integrated numerically.  However, $\s$ itself only has one parameter, so it will be faster to evaluate if all terms involving $\s$ can be moved to a 1D integral.  Noting that $t_r=t-\frac1c\abs{\x-\y}$, it seems like a good idea to reparameterise so that $\abs{\x-\y}$ is one of the variables.  This can be done by using spherical polar coordinates centered around $\x$, with $\de^2\Omega$ representing solid angle:
\[ \de^3\y = \abs{\x-\y}^2\de\abs{\x-\y}\de^2\Omega \]
\begin{eqnarray*}
\Rightarrow \qquad 
 A^\mu(\x,t) &=& \frac{\mu_0q}{4\pi} \iint V^\mu(t_r) \frac{f(\abs{\y-\s(t_r)})}{\abs{\x-\y}} \abs{\x-\y}^2\de\abs{\x-\y}\de^2\Omega \\
  &=& \frac{\mu_0q}{4\pi} \int V^\mu(t_r)\abs{\x-\y} \left(\int f(\abs{\y-\s(t_r)}) \de^2\Omega\right) \de\abs{\x-\y}.
\end{eqnarray*}
Within the inner integral, $\s(t_r)$ is a constant since $t_r$ only depends on $\abs{\x-\y}$.  If $S(\x,r)$ is the sphere of radius $r$ centred on $\x$, the inner integral can be written
\begin{eqnarray*}
\int_{\y\in S(\x,\abs{\x-\y})}  f(\abs{\y-\s(t_r)}) \de^2\Omega 
 &=& 
\int_{\y\in S(\x-\s(t_r),\abs{\x-\y})}  f(\abs{\y}) \de^2\Omega \\
 &=& 
 I(\x-\s(t_r),\abs{\x-\y}),
\end{eqnarray*}
where $I(\x,r)$ is defined to integrate $f$ over solid angle on a sphere.

\subsubsection{Spherical Integral}\label{sphint}
As $f$ only depends on the modulus of its argument, the spherical integral can be reduced to a single variable if the $\theta=0$ pole of the spherical polar coordinates is chosen to point away from the origin:
\begin{eqnarray*}
 I(\x,r) &=& \int_{\y\in S(\x,r)} f(\abs{\y}) \de^2\Omega \\
 &=& \int_{\theta=0}^\pi \int_{\phi=0}^{2\pi} f\left(\sqrt{(\abs{\x}+r\cos\theta)^2+(r\sin\theta)^2}\right) \sin\theta \de\theta\de\phi  \\
 &=& 2\pi \int_0^\pi f\left(\sqrt{|\x|^2+2\abs{\x}r\cos\theta+r^2}\right) \sin\theta \de\theta.
\end{eqnarray*}
With the substitution $p=\cos\theta$, the differential is $\de p = -\sin\theta \de\theta$, simplifying things further:
\begin{eqnarray*}
 I(\x,r) &=& 2\pi \int_1^{-1}  f\left(\sqrt{|\x|^2+2\abs{\x}rp+r^2}\right) (-\de p) \\
 &=& 2\pi \int_{-1}^1 f\left(\sqrt{a+bp}\right) \de p,
\end{eqnarray*}
where $a=\abs{\x}^2+r^2$ and $b=2\abs{\x}r$.  Note that $a\ge 0$ and $\abs{b}\le a$ so the $a+bp$ under the square root never becomes negative.

A common form for the smoothed charge density would be a 3D spherical Gaussian
\begin{eqnarray*}
f_{3,\sigma}(\abs{\x}) &=& \frac{1}{(\sigma\sqrt{2\pi})^3}e^{-|\x|^2/2\sigma^2} \\
\Rightarrow \qquad 
f_{3,\sigma}\left(\sqrt{a+bp}\right) &=& \frac{1}{(\sigma\sqrt{2\pi})^3}e^{-(a+bp)/2\sigma^2} \\
\Rightarrow \qquad 
 I(\x,r) &=& \frac{2\pi}{(\sigma\sqrt{2\pi})^3} e^{-a/2\sigma^2} \int_{-1}^1 e^{-bp/2\sigma^2} \de p \\
 &=& \frac{1}{\sigma^3\sqrt{2\pi}} e^{-a/2\sigma^2} \frac{2\sigma^2}{-b}\left( e^{-b/2\sigma^2} - e^{b/2\sigma^2} \right) \\
 &=& \frac{2}{\sigma\sqrt{2\pi}b} \left( e^{-(a-b)/2\sigma^2} - e^{-(a+b)/2\sigma^2} \right).
\end{eqnarray*}
Noting that $a+b=(\abs{\x}+r)^2$ and $a-b=(\abs{\x}-r)^2$,
\[ I(\x,r) = \frac{e^{-(\abs{\x}-r)^2/2\sigma^2} - e^{-(\abs{\x}+r)^2/2\sigma^2} }{\sigma\sqrt{2\pi}\abs{\x}r}. \]

More generally, if $g(x)=f(\sqrt{x})$ and this function has an integral $G$ such that $G'(x)=g(x)$,
\begin{eqnarray*}
 I(\x,r) &=& 2\pi \int_{-1}^1 g(a+bp) \de p \\
 &=& \frac{2\pi}{b} \left(G(a+b)-G(a-b)\right) \\
 &=& \frac{\pi}{\abs{\x}r} \left(G((\abs{\x}+r)^2)-G((\abs{\x}-r)^2)\right).
\end{eqnarray*}


\subsubsection{Single-Parameter Integral}\label{onedint}
The remaining 1D integral will have to be integrated numerically over the range where $I$ is non-negligible:
\begin{eqnarray*}
A^\mu(\x,t) &=& \frac{\mu_0q}{4\pi} \int V^\mu(t_r)\abs{\x-\y} I(\x-\s(t_r),\abs{\x-\y}) \de\abs{\x-\y} \\
&=& \frac{\mu_0q}{4\pi} \int_0^\infty \left(c,\dot\s\left(t-\frac{r}{c}\right)\right)^\mu r I\left(\x-\s\left(t-\frac{r}{c}\right),r\right) \de r,
\end{eqnarray*}
where for compactness $r=\abs{\x-\y}$ has been defined, making $t_r=t-\frac{r}{c}$.  In fact, the variable may now be changed to $t_r$:
\begin{eqnarray*}
A^\mu(\x,t) &=& \frac{\mu_0q}{4\pi} \int_{-\infty}^t (c,\dot\s(t_r))^\mu c(t-t_r) I(\x-\s(t_r),c(t-t_r)) c \de t_r \\
 &=& \frac{q}{4\pi\epsilon_0} \int_{-\infty}^t (c,\dot\s(t_r))^\mu (t-t_r) I(\x-\s(t_r),c(t-t_r)) \de t_r.
\end{eqnarray*}

\subsubsection{Integration Range}\label{intrange}
Suppose the smoothed charge satisfies $f(\abs{\y})\simeq0$ for $\abs{\y} > R$.  The general form of $I(\x,r)$ integrates it over the range $\abs{\abs{\x}-r} \le \abs{\y} \le \abs{\x}+r$, noting that there are two cases depending if the sphere $S(\x,r)$ encloses the origin.  For $I$ to include only `negligible' values of $f$, the lower end of this range must be greater than $R$, or $\abs{\abs{\x}-r} > R$.

In the integral for $A^\mu(\x,t)$, this means the range of non-negligible values is where
\[ \abs{\abs{\x-\s(t_r)}-c(t-t_r)} \le R. \]
The distance $\abs{\x-\s(t_r)}$ changes at less than $c$ for a slower-than-light particle, so $\abs{\x-\s(t_r)}-c(t-t_r)$ is always increasing with $t_r$.  This means the integral is a single range between $t_{r-}$ and $t_{r+}$ defined by 
\begin{eqnarray*}
\abs{\x-\s(t_{r\pm})}-c(t-t_{r\pm}) &=& \pm R \\
\Leftrightarrow \qquad
\frac1c\abs{\x-\s(t_{r\pm})}-(t-t_{r\pm}) &=& \pm \frac{R}{c} \\
\Leftrightarrow \qquad
t_{r\pm} &=& t\pm\frac{R}{c}-\frac1c\abs{\x-\s(t_{r\pm})}\\
 &=& t_r\left(\x,t\pm\frac{R}{c},\s(t_{r\pm})\right) \\
\Leftrightarrow \qquad t_{r\pm} &=& t_{r,\s}\left(\x,t\pm\frac{R}{c}\right).
\end{eqnarray*}
Here $t_{r,\s}(\x,t)$ is the retarded time calculated by root finding on the $\s$ trajectory as used in the formula for a point particle.  For $t_{r+}$, the particle might not have left the past light cone of $(\x,t+\frac{R}{c})$ yet, so if a root is not found for this reason, set $t_{r+}=t$.

\subsubsection{Check: Potential of a Stationary Gaussian Smoothed Charge}
The formula in section \ref{onedint} may be used to derive potentials for a smoothed charge at the origin that is not moving ($\s=\mathbf 0$).  The spatial elements are obviously zero, leaving
\begin{eqnarray*}
A^0(\x,t) = \phi/c &=& \frac{q}{4\pi\epsilon_0} \int_{-\infty}^t c (t-t_r) I(\x,c(t-t_r)) \de t_r \\
 &=& \frac{q}{4\pi\epsilon_0} \int_{-\infty}^t c (t-t_r) \frac{e^{-(\abs{\x}-c(t-t_r))^2/2\sigma^2} - e^{-(\abs{\x}+c(t-t_r))^2/2\sigma^2} }{\sigma\sqrt{2\pi}\abs{\x}c(t-t_r)} \de t_r \\
 &=& \frac{q}{4\pi\epsilon_0\sigma\sqrt{2\pi}\abs{\x}} \int_{-\infty}^t e^{-(\abs{\x}-c(t-t_r))^2/2\sigma^2} - e^{-(\abs{\x}+c(t-t_r))^2/2\sigma^2} \de t_r \\
 &=& \frac{q}{4\pi\epsilon_0\sigma\sqrt{2\pi}\abs{\x}} \int_0^\infty e^{-(\abs{\x}-c\Delta t)^2/2\sigma^2} - e^{-(\abs{\x}+c\Delta t)^2/2\sigma^2} \de \Delta t \\
 &=& \frac{q}{4\pi\epsilon_0\sigma\sqrt{2\pi}\abs{\x}} \left( \frac1c \int_{-\infty}^{\abs{\x}} e^{-x^2/2\sigma^2} \de x - \frac1c \int_{\abs{\x}}^\infty  e^{-x^2/2\sigma^2} \de x \right).
\end{eqnarray*}
 At this point, it is useful to define the cumulative function of the standard normal distribution:
\begin{eqnarray*}
\Phi(z) &=& \frac1{\sqrt{2\pi}}\int_{-\infty}^z e^{-Z^2/2}\de Z \\
 &=& \frac1{\sigma\sqrt{2\pi}}\int_{-\infty}^{z\sigma} e^{-x^2/2\sigma^2}\de x
\end{eqnarray*}
\begin{eqnarray*}
\Rightarrow \qquad
A^0(\x,t) = \phi/c &=& \frac{q}{4\pi\epsilon_0\abs{\x}c} \left( \Phi\left(\frac{\abs{\x}}{\sigma}\right) - \left( 1 - \Phi\left(\frac{\abs{\x}}{\sigma}\right) \right) \right) \\
 &=& \frac{q}{4\pi\epsilon_0\abs{\x}c} \left( 2\Phi\left(\frac{\abs{\x}}{\sigma}\right) - 1 \right). 
\end{eqnarray*}
For a spherically symmetrical potential, we expect
\[ \nabla^2 \phi = \phi''(r)+\frac2{r}\phi'(r) = -\frac{\rho}{\epsilon_0}. \]
The radial derivatives of $\phi$ are:
\begin{eqnarray*}
\phi(r) &=&  \frac{q}{4\pi\epsilon_0}\frac1r\left(2\Phi\left(\frac{r}{\sigma}\right) - 1 \right) \\
\phi'(r) &=&  \frac{q}{4\pi\epsilon_0}\left( \frac{-1}{r^2}\left(2\Phi\left(\frac{r}{\sigma}\right) - 1 \right) + \frac2{r\sigma}\Phi'\left(\frac{r}{\sigma}\right) \right) \\
\phi''(r) &=&  \frac{q}{4\pi\epsilon_0}\left( \frac{2}{r^3}\left(2\Phi\left(\frac{r}{\sigma}\right) - 1 \right) -  \frac{2}{r^2\sigma}\Phi'\left(\frac{r}{\sigma}\right) - \frac2{r^2\sigma}\Phi'\left(\frac{r}{\sigma}\right) + \frac{2}{r\sigma^2}\Phi''\left(\frac{r}{\sigma}\right) \right) \\
\end{eqnarray*}
\begin{eqnarray*}
\Rightarrow \qquad
\phi''(r)+\frac2{r}\phi'(r) &=& \frac{q}{4\pi\epsilon_0} \frac{2}{r\sigma^2}\Phi''\left(\frac{r}{\sigma}\right) \\
&=& \frac{q}{4\pi\epsilon_0} \frac{2}{r\sigma^2} \frac1{\sqrt{2\pi}} \left(-\frac{r}{\sigma}\right) e^{-(r/\sigma)^2/2}
\end{eqnarray*}
\begin{eqnarray*}
\Rightarrow \qquad
\rho &=& \frac{q}{4\pi} \frac{2}{r\sigma^2} \frac1{\sqrt{2\pi}} \frac{r}{\sigma} e^{-(r/\sigma)^2/2} \\
&=& \frac{q}{(\sigma\sqrt{2\pi})^3}  e^{-r^2/2\sigma^2} = qf_{3,\sigma}(r),
\end{eqnarray*}
as required.

\section{Fields (Smoothed Charge)}
The Li\'enard--Wiechert fields for a moving point charge are well-known and do not seem to have a strong analogy with the derivation for a smoothed charge, so only the latter is given here.

\subsection{Maxwell's Equations}
The electromagnetic fields are defined from the potential via
\[ F_{\mu\nu} = \partial_\mu A_\nu - \partial_\nu A_\mu,
\qquad \mbox{where} \qquad
F_{0i}=E_i/c,  \qquad F_{32,13,21}=B_{x,y,z} \]
and the index-lowered four-potential is $A_\mu=(\phi/c,-\A)$.  This index lowering operation is merely multiplication by the $4\times4$ matrix $\mathrm{diag}(1,-1,-1,-1)$.  For example, there is a lowered form of the retarded solution of Maxwell's equations in the Lorenz gauge from section \ref{maxpot}:
\[ A_\mu(\x,t) = \frac{\mu_0}{4\pi} \int \frac{J_\mu(\y,t_r)}{\abs{\x-\y}}\de^3\y, \]
where $J_\mu=(c\rho,-\J)$.  Taking the derivative with respect to $x^\mu=(ct,\x)$ gives
\begin{eqnarray*}
 \partial_\mu A_\nu(\x,t)
  &=& \frac{\mu_0}{4\pi} \int \partial_\mu\left(\frac1{\abs{\x-\y}}\right)J_\nu(\y,t_r) + \frac{\dot J_\nu(\y,t_r)\partial_\mu t_r(\x,t,\y)}{\abs{\x-\y}}\de^3\y \\
  &=& \frac{\mu_0}{4\pi} \int \frac{(0,\y-\x)_\mu}{\abs{\x-\y}^3}J_\nu(\y,t_r) + \frac{(\abs{\x-\y},\y-\x)_\mu}{c\abs{\x-\y}^2}\dot J_\nu(\y,t_r)\de^3\y. 
\end{eqnarray*}
Evaluating $\partial_\mu A_\nu - \partial_\nu A_\mu$ with the formula above and identifying the $E$ and $B$ components gives Jefimenko's equations, which give the fields directly from $J$ and $\dot J$.  However, for the smoothed charge case, the potential has already been simplified into a 1D integral in section \ref{onedint}, so it is simpler to differentiate that rather than go all the way back to the source terms.

\rem
{
\subsection{Time Derivative of Source}
The index-lowered form of the four-current density for the smoothed charge is
\[ J_{s\mu}(\x,t) = qV_\mu(t)f(\abs{\x-\s(t)}) \qquad \mbox{where} \qquad V_\mu(t)=(c,-\dot\s(t)). \]
Its time derivative is
\[ \dot J_{s\mu}(\x,t) = q\dot V_\mu(t)f(\abs{\x-\s(t)}) + qV_\mu(t)f'(\abs{\x-\s(t)})\frac{\s(t)-\x}{\abs{\s(t)-\x}}\cdot\dot\s(t), \]
where $\dot V_\mu(t)=(0,-\ddot\s(t))$.
}

%\subsection{Derivation for Point Particle}

\subsection{Derivation for Smoothed Charge}
\rem
{
Substituting the smoothed four-current density into the expression for $\partial_\mu A_\nu$ gives
\begin{eqnarray*}
 \partial_\mu A_\nu(\x,t)
  &=& \frac{\mu_0q}{4\pi} \int V_\nu(t_r) \frac{(0,\y-\x)_\mu}{\abs{\x-\y}^3} f(\abs{\y-\s(t_r)}) \de^3\y \\
  &+& \frac{\mu_0q}{4\pi} \int \dot V_\nu(t_r) \frac{(\abs{\x-\y},\y-\x)_\mu}{c\abs{\x-\y}^2} f(\abs{\y-\s(t_r)}) \de^3\y \\
  &+& \frac{\mu_0q}{4\pi} \int V_\nu(t_r) \frac{(\abs{\x-\y},\y-\x)_\mu}{c\abs{\x-\y}^2} f'(\abs{\y-\s(t_r)})\frac{\s(t_r)-\y}{\abs{\s(t_r)-\y}}\cdot\dot\s(t_r) \de^3\y
\end{eqnarray*}
}

The smoothed charge potential derived in section \ref{onedint} can have its indices lowered:
\[ A_\mu(\x,t) = \frac{q}{4\pi\epsilon_0} \int_{-\infty}^t (c,-\dot\s(t_r))_\mu (t-t_r) I(\x-\s(t_r),c(t-t_r)) \de t_r. \]
Remembering that $t_r$ is now an integration variable, it has no dependence on $\x$ and $t$.  Taking the derivative with respect to $x^\mu=(ct,\x)$ will also change the upper integration limit, which would normally produce an extra term, but the integrand is zero there because of the factor $(t-t_r)$ in the product.  Derivatives of other parts are
\[ \partial_\mu(t-t_r) = \left(\frac1c,\mathbf0\right)_\mu \qquad \qquad \partial_\mu(\x-\s(t_r)) = (\mathbf 0,\e_x,\e_y,\e_z)_\mu \]
where $\e_i$ are the unit vectors in each direction.  If the gradient of the function $I(\x,r)$ with respect to its first argument is written $\nabla I$ and its partial derivative with respect to the second argument as $\partial_r I$, then
\begin{eqnarray*}
 \partial_\mu I(\x-\s(t_r),c(t-t_r)) &=& \nabla I\cdot\partial_\mu(\x-\s(t_r)) + \partial_r Ic\partial_\mu(t-t_r) \\
 &=& (0,\nabla I)_\mu + \partial_r I (1,0)_\mu \\
 &=& (\partial_r I,\nabla I)_\mu,
\end{eqnarray*}
where $\nabla I$ and $\partial_r I$ are evaluated at the same point as $I$.  Finally, the derivative of the potential can be written
\[ \partial_\mu A_\nu(\x,t) = \frac{q}{4\pi\epsilon_0} \int_{-\infty}^t (c,-\dot\s(t_r))_\nu \left( \left(\frac1c,\mathbf0\right)_\mu I(\x-\s(t_r),c(t-t_r)) + (t-t_r)(\partial_r I,\nabla I)_\mu \right) \de t_r. \]
% &=& \frac{q}{4\pi\epsilon_0} \int_{-\infty}^t \left(\frac{I(\x-\s(t_r),c(t-t_r))}{c}+ (t-t_r)\partial_r I,(t-t_r)\nabla I\right)_\mu (c,-\dot\s(t_r))_\nu \de t_r.
To calculate the fields on a computer, this integral of a $4\times4$ matrix is evaluated over the range found in section \ref{intrange} and the $F_{\mu\nu} = \partial_\mu A_\nu - \partial_\nu A_\mu$ combinations formed afterwards to get the $E$ and $B$ fields. % Could do the E and B explicitly here like I do in the code

% Can separate I's derivatives into contact/near-field/far-field parts.  The far-field part can be expressed as Gaussian'(f(t_r)).  Change of variable, move deriv, change back gives a more well-conditioned version of this integral

A disadvantage of using this integral directly is that in the far field region, the derivatives of $I$ contain a part that is mostly an odd function, which has large positive and then almost equally large negative values as the smoothed charge goes through the past light cone.  This leads to loss of accuracy on the computer when the large numbers are cancelled.  A solution to this will be explored in section \ref{farfieldint}.

\rem
{ % Probably doesn't work except under special symmetry conditions
\subsection{Integral Trick for Finite Timesteps}
In a particle simulation, it may be more appropriate to apply the average field over a timestep so that brief spikes in field intensity are not missed.  If the timestep runs from $t_1$ to $t_2$, this average is given by:
\begin{eqnarray*}
\bar F_{\mu\nu} = \frac1{t_2-t_1}\int_{t_1}^{t_2} F_{\mu\nu}(\s(t),t) \de t
\end{eqnarray*}
}

\subsection{Derivatives of $I$ for a 3D Gaussian}
It was derived in section \ref{sphint} that a 3D Gaussian smoothed charge has
\[ I(\x,r) = \frac{1}{\sigma\sqrt{2\pi}\abs{\x}r}\left( e^{-(\abs{\x}-r)^2/2\sigma^2} - e^{-(\abs{\x}+r)^2/2\sigma^2} \right). \]
The gradient with respect to $\x$ is simplified by noting $I$ is a function of $\abs{\x}$ and $r$ only:
\begin{eqnarray*}
 \nabla I(\x,r) &=& \partial_{\abs{\x}} I(\abs{\x},r) \frac{\x}{\abs{\x}} \\
 &=& \frac{-1}{\sigma\sqrt{2\pi}\abs{\x}^2r}\left( e^{-(\abs{\x}-r)^2/2\sigma^2} - e^{-(\abs{\x}+r)^2/2\sigma^2} \right) \frac{\x}{\abs{\x}} \\
 &+& \frac{1}{\sigma\sqrt{2\pi}\abs{\x}r}\left( \frac{-2(\abs{\x}-r)}{2\sigma^2} e^{-(\abs{\x}-r)^2/2\sigma^2} -  \frac{-2(\abs{\x}+r)}{2\sigma^2} e^{-(\abs{\x}+r)^2/2\sigma^2} \right) \frac{\x}{\abs{\x}} \\
 &=& \frac{\x}{\sigma\sqrt{2\pi}\abs{\x}^2r}
 \left(
 \left(-\frac{1}{\abs{\x}} - \frac{\abs{\x}-r}{\sigma^2} \right)e^{-(\abs{\x}-r)^2/2\sigma^2} + 
 \left(\frac{1}{\abs{\x}} + \frac{\abs{\x}+r}{\sigma^2} \right)e^{-(\abs{\x}+r)^2/2\sigma^2}
 \right).
\end{eqnarray*}
The partial derivative with respect to $r$ is
\begin{eqnarray*}
\partial_r I(\x,r) &=& \frac{-1}{\sigma\sqrt{2\pi}\abs{\x}r^2}\left( e^{-(\abs{\x}-r)^2/2\sigma^2} - e^{-(\abs{\x}+r)^2/2\sigma^2} \right) \\
&+&  \frac{1}{\sigma\sqrt{2\pi}\abs{\x}r}\left( \frac{2(\abs{\x}-r)}{2\sigma^2}e^{-(\abs{\x}-r)^2/2\sigma^2} - \frac{-2(\abs{\x}+r)}{2\sigma^2}e^{-(\abs{\x}+r)^2/2\sigma^2} \right) \\
&=& \frac{1}{\sigma\sqrt{2\pi}\abs{\x}r}
\left( 
\left(-\frac1r + \frac{\abs{\x}-r}{\sigma^2} \right)e^{-(\abs{\x}-r)^2/2\sigma^2} + 
\left(\frac1r + \frac{\abs{\x}+r}{\sigma^2} \right)e^{-(\abs{\x}+r)^2/2\sigma^2} 
\right).
\end{eqnarray*}
Since all of these quantities need to be calculated together, on a computer it is efficient to first calculate
\[ k = \frac{1}{\sigma\sqrt{2\pi}\abs{\x}r} \qquad 
e_- =  e^{-(\abs{\x}-r)^2/2\sigma^2} \qquad
e_+ =  e^{-(\abs{\x}+r)^2/2\sigma^2} \]
and then
\begin{eqnarray*}
I &=& k(e_- - e_+) \\
\nabla I &=& k
 \left(
 \left(-\frac{1}{\abs{\x}} - \frac{\abs{\x}-r}{\sigma^2} \right)e_- + 
 \left(\frac{1}{\abs{\x}} + \frac{\abs{\x}+r}{\sigma^2} \right)e_+
 \right) \frac{\x}{\abs{\x}} \\
 \partial_r I &=& k
\left( 
\left(-\frac1r + \frac{\abs{\x}-r}{\sigma^2} \right)e_- + 
\left(\frac1r + \frac{\abs{\x}+r}{\sigma^2} \right)e_+ 
\right).
\end{eqnarray*}

\subsection{Derivatives of $I$ for General Case}
Using the definitions of $g(x)=f(\sqrt{x})$ and $G'=g$ from section \ref{sphint},
\[ I(\x,r) = \frac{\pi}{\abs{\x}r} \left(G((\abs{\x}+r)^2)-G((\abs{\x}-r)^2)\right). \]
The partial derivative with respect to $r$ is
\begin{eqnarray*}
\partial_r I(\x,r) &=& \frac{-1}{r}I + \frac{\pi}{\abs{\x}r} \left(2(\abs{\x}+r)g((\abs{\x}+r)^2) + 2(\abs{\x}-r)g((\abs{\x}-r)^2)\right) \\
&=& \frac{-1}{r}I + \frac{2\pi}{\abs{\x}r} \left((\abs{\x}+r)f(\abs{\x}+r) + (\abs{\x}-r)f(\abs{\x}-r)\right),
\end{eqnarray*}
where negative arguments of $f$ are defined by $f(-x)=f(x)$.  Interpreting $I(\abs{\x},r)$ with scalar arguments, note that $I(x,y)=I(y,x)$ so that the other partial derivative is the same with the arguments swapped over:
\begin{eqnarray*}
 \nabla I(\x,r) &=& \partial_{\abs{\x}} I(\abs{\x},r) \frac{\x}{\abs{\x}} \\
&=& \left(
 \frac{-1}{\abs{\x}}I + \frac{2\pi}{\abs{\x}r} \left((\abs{\x}+r)f(\abs{\x}+r) - (\abs{\x}-r)f(\abs{\x}-r)\right)
\right) \frac{\x}{\abs{\x}}.
\end{eqnarray*}

\subsection{Numerically Stable Integration of Far Field Terms}\label{farfieldint}
The (additive) part of the fields coming from derivatives of $I$ is 
\[ (\partial_\mu A_\nu)_\mathrm{derivs}(\x,t) = \frac{q}{4\pi\epsilon_0} \int_{-\infty}^t (c,-\dot\s(t_r))_\nu (t-t_r)(\partial_r I,\nabla I)_\mu \de t_r, \]
where $I$ and its derivatives are evaluated at $(\y,r)=(\x-\s(t_r),c(t-t_r))$.  Breaking this down further to just the term including the odd function $(\abs{\y}-r)f(\abs{\y}-r)$ gives the far field part % When I break dI into parts, are the parts still well-behaved like I is?
\begin{eqnarray*}
(\partial_\mu A_\nu)_\mathrm{ff}(\x,t)
 &=& \frac{q}{4\pi\epsilon_0} \int_{-\infty}^t (c,-\dot\s(t_r))_\nu (t-t_r)\frac{2\pi}{\abs{\y}r}\left(1,\frac{-\y}{\abs{\y}}\right)_\mu (\abs{\y}-r)f(\abs{\y}-r) \de t_r \\
 &=& \frac{q}{4\pi\epsilon_0} \int_{-\infty}^t \frac{2\pi}{c} (c,-\dot\s(t_r))_\nu \frac{1}{\abs{\y}}\left(1,\frac{-\y}{\abs{\y}}\right)_\mu (\abs{\y}-r)f(\abs{\y}-r) \de t_r.
\end{eqnarray*}
Note that the distance $z=\abs{\y}-r=\abs{\x-\s(t_r)}-c(t-t_r)$ is a monotonically increasing function of $t_r$ for slower-than-light particles.  The odd function $zf(z)$ will have an even integral $H(z)$ such that $H'(z)=zf(z)$ and $H\rightarrow0$ as $\abs{z}\rightarrow\infty$.  Finally collect the various front terms together by defining
\[ k_{\mu\nu}(t_r) =\frac{2\pi}{c}(c,-\dot\s(t_r))_\nu \left(\frac{1}{\abs{\y}},\frac{-\y}{\abs{\y}^2}\right)_\mu. \]
The far field now has the form
\begin{eqnarray*}
\frac{4\pi\epsilon_0}{q}(\partial_\mu A_\nu)_\mathrm{ff}(\x,t)
 &=& \int_{-\infty}^t k_{\mu\nu}(t_r) H'(z(t_r)) \de t_r \\
 &=& \int_{-\infty}^{z=\abs{\x-\s(t)}} k_{\mu\nu} H'(z) \dd{t_r}{z} \de z \\
 &=& \left[ H(z)  k_{\mu\nu} \dd{t_r}{z} \right]_{-\infty}^{z=\abs{\x-\s(t)}} - \int_{-\infty}^{z=\abs{\x-\s(t)}} H(z) \dd{}{z}\left(k_{\mu\nu} \dd{t_r}{z}\right) \de z \\
 &=& \left. H(\abs{\x-\s(t)})  k_{\mu\nu} \dd{t_r}{z} \right|_t - \int_{-\infty}^t H(z(t_r)) \dd{}{t_r}\left(k_{\mu\nu} \dd{t_r}{z}\right) \de t_r.
\end{eqnarray*}
The function that contains detail on the scale of the charge size (e.g. $f$, $g$, $G$, $H$) within the integral is now even.  The rest is algebraically slightly messy but varies with distance from the source and its velocity and acceleration:
\begin{eqnarray*}
\dd{}{t_r}\left(k_{\mu\nu} \dd{t_r}{z}\right)
 &=& \dd{}{t_r}\left(k_{\mu\nu} \left(\dd{z}{t_r}\right)^{-1} \right) \\
 &=& \dd{k_{\mu\nu}}{t_r}\left(\dd{z}{t_r}\right)^{-1} - k_{\mu\nu} \left(\dd{z}{t_r}\right)^{-2} \dds{z}{t_r}; \\
\dd{z}{t_r} &=& -\frac{\y}{\abs{\y}}\cdot\dot\s(t_r) + c; \\
\dds{z}{t_r} &=& \frac{\abs{\dot\s(t_r)}^2 - (\frac{\y}{\abs{\y}}\cdot\dot\s(t_r))^2 - \y\cdot\ddot\s(t_r)}{\abs{\y}}; \\
\dd{k_{\mu\nu}}{t_r} &=&
\frac{2\pi}{c}
\left(
(0,-\ddot\s(t_r))_\nu \left(\frac{1}{\abs{\y}},\frac{-\y}{\abs{\y}^2}\right)_\mu
\right.\\&+&\left.
(c,-\dot\s(t_r))_\nu \left(\frac{\y\cdot\dot\s(t_r)}{\abs{\y}^3},\frac{\dot\s(t_r)}{\abs{\y}^2} - \frac{2(\y\cdot\dot\s(t_r))\y}{\abs{\y}^4}\right)_\mu
\right).
\end{eqnarray*}

For the 3D Gaussian, the $H$ function is the integral of 
\begin{eqnarray*}
xf_{3,\sigma}(x) &=& \frac{1}{(\sigma\sqrt{2\pi})^3}xe^{-x^2/2\sigma^2} \\
\Rightarrow \qquad H(x) &=& \frac{1}{(\sigma\sqrt{2\pi})^3}(-\sigma^2)e^{-x^2/2\sigma^2}
 = -\sigma^2 f_{3,\sigma}(x).
\end{eqnarray*}

%\subsection{Integration Step Sizes}
% See what I did in CSRTrack_gauss_field 

\end{document}