## 3 Integro-Dipolar Expansion  
We divide the \(\mathbf{y}\) integration region using cubes centered at the origin. When we compute the contributions of the regions with \(0\leq y_{i}\leq 10^{-6}\), then \(10^{-6}\leq y_{i}\leq 10^{-5}\) etc. (values in cm), we obtain precise results up to approximately \(\,10^{-5}\); then the fluctuations become large, even in long runs (\(10^{11}\) to \(10^{12}\) sampling points). This happens because the two opposite monopolar contributions in the integral are large and when the sampling points are spread over bigger volumes, their cancellation is affected by large casual errors. We therefore make recourse to a dipolar expansion in the integration region far from the primary source. This is a non-standard expansion because of the presence of the double retarded integration, so it needs special care and must be cross-checked numerically by comparing its results to those of the full integral in the intermediate integration region where \(|{\bf y}|\) is small enough that the fluctuations are still under control but large enough that the assumption \(|{\bf y}|\gg a\) for the dipolar expansion is valid. Let us first consider the case of dipole charges having spherical symmetry, so that \(d=\varepsilon\) in the definition of \(f({\bf z})\). We rewrite the integral for \(E_{i}^{\phi}\) as the sum of two integrals \(E_{i}^{\phi,{\bf a}}\) and \(E_{i}^{\phi,-{\bf a}}\) for the sources at \({\bf a}\) and \(-{\bf a}\), in which the \({\bf z}\) variable is shifted by \(-{\bf a}\) and \({\bf a}\), respectively:  
\[E_{i}^{\phi}=E_{i}^{\phi,{\bf a}}+E_{i}^{\phi,-{\bf a}} \tag{24}\]  
For the first integral, with shift \(-{\bf a}\), we define a new variable \({\bf u}={\bf z}-{\bf a}\). Define a regularized \(\delta\)-function for a source centered at the origin:  
\[F({\bf u})=\frac{1}{\sqrt{(2\pi)^{3}}\varepsilon^{2}}e^{-\frac{1}{2}\frac{{ \bf a}^{2}}{\varepsilon^{2}}} \tag{25}\]  
The electric field generated by the scalar potential of the source at \({\bf a}\) can be written as  
\[E_{i}^{\phi,{\bf a}}=-qK^{2}\int_{M-{\bf a}}d^{3}u\int d^{3}yH_{i}^{\phi}({\bf x },{\bf y},{\bf u})G^{\phi}(t,{\bf x},{\bf y},{\bf u}+{\bf a}) \tag{26}\]  
where  
\[H_{i}^{\phi}({\bf x},{\bf y},{\bf u})=\frac{F({\bf u})(x_{i}-y_{i})}{|{\bf x}- {\bf y}|^{3}} \tag{27}\]  
\[G^{\phi}(t,{\bf x},{\bf y},{\bf u}+{\bf a})=\frac{K|{\bf x}-{\bf y}|\sin\Omega -\cos\Omega}{|{\bf y}-{\bf u}-{\bf a}|} \tag{28}\]  
\[\Omega=\omega(t-k|{\bf x}-{\bf y}|-k|{\bf y}-{\bf u}-{\bf a}|) \tag{29}\]  
We have symbolically denoted the integration range of \({\bf u}\) as \"\(M-{\bf a}\)\", meaning that it is equal to the integration range \(M\) of \(z\) (\(-R_{z}\leq z_{i}\leq R_{z}\)) shifted by a quantity \(-{\bf a}\). The function \(G^{\phi}(t,{\bf x},{\bf y},{\bf u}+{\bf a})\) can be expanded as a term of order zero in \(a=|{\bf a}|\) and a term of order 1:  
\[G^{\phi}(t,{\bf x},{\bf y},{\bf u}+{\bf a})\simeq G_{0}^{\phi}(t,{\bf x},{\bf y },{\bf u})+a\cdot G_{1}^{\phi}(t,{\bf x},{\bf y},{\bf u}) \tag{30}\]  
Actually, the small quantity in which we make the expansion is \(a/|{\bf y}|\) and we therefore expect that the expansion is accurate where \(|{\bf y}|\gg a\), which is what we need, as explained above. Let us expand the factor \(1/|{\bf y}-{\bf u}-{\bf a}|\) to first order in \(a\). Define \({\bf v}={\bf y}-{\bf u}\). \(|{\bf v}|\) is of order \(|{\bf y}|\), because \(F({\bf u})\) has range \(\simeq\varepsilon<a\); therefore \(|{\bf v}|\gg a\). In the following we denote \(v=|{\bf v}|\).