I'm having some trouble figuring out the properties with respect to analytical continuation of functions defined using an integral kernel. More particularly, I am working with the electrostatic potential $$ V(x,y,z)=\int \frac{\rho_C(x',y',z')dx'dy'dz'}{\sqrt{(x-x')^2+(y-y')^2+(z-z')^2}} $$ where the charge density $\rho_C$ describes a localized system.

(Scroll down for my specific question.)

My initial expectation is that, if $\rho_C$ decays fast enough and contains no singular, delta-function-like points, then the potential $V$ will be a real analytic function of its variables $x,y,z$, inheriting that behaviour from the kernel. Because of this, I would expect $V$ to extend to an analytical function, probably on the whole complex plane. I additionally expect this analytical continuation to contain either branch cuts or singularities, or to tend to infinity for large complex arguments in some direction, because of Liouville's Theorem.

In my particular case, $\rho_C$ is a molecular transition charge, which means that my knowledge about it must be numerically obtained through quantum chemistry methods, and the integral must be done numerically. I expect the real, physical $\rho_C$ to decay as $e^{-r}$ at large distances. On the other hand, the quantum chemical code approximates the charge density as a superposition of gaussians, $e^{-r^2}$.

Both of these simple cases I can solve for analytically for real coordinates, using their spherical symmetry and Gauss's law. The results are deceptively simple: for an exponential charge density the potential goes as $$V_\text{exp}\propto e^{-r}\left(1+\frac{2}{r}\right)-\frac{2}{r},$$ while for a gaussian distribution $$V_\text{gauss}\propto \frac{\text{erf}(r)}{r}.$$ On the real line, both potentials look essentially identical. However, their behaviour when extended to complex coordinates is wildly different:

- On one hand, $\text{erf}$ is an odd, entire function, so $V_\text{gauss}$ is a function of $r^2=x^2+y^2+z^2$ and has therefore no branch cuts anywhere in any variable. In contrast, $V_\text{exp}$ is a function of the radial distance $r=\sqrt{x^2+y^2+z^2}$, so if one keeps $x,y$ real then $V_\text{exp}$ as a function of $z$ has branch cuts starting at $\pm i\sqrt{x^2+y^2}$, and its Riemann surface has two sheets.
- On top of that, $V_\text{exp}$ is bounded except for the neighbourhood of the branch cut, while $V_\text{gauss}$ blows up as $e^{+|z|^2}$ for large imaginary $z$.

While these differences weird me out slightly, I can accept them as representing the very different origins of the two potentials. The difference in behaviour for $\text{Re}(x^2+y^2+z^2)<0$ is worrisome, as it means the quantum chemical code cannot be trusted there, but I can accept that as inevitable as long as I have warning flags that indicate the code is out of its validity region.

(The above is mostly background, but if I have made conceptual errors or I am working on misconceptions I would dearly like to know. My precise question follows)

What really worries me, though, is the last point on the list: $V_\text{gauss}$ blows up as $e^{+|z|^2}$ for large imaginary $z$. This means that for $x=y=0$ and $z=i\zeta$, the integral $$ V_\text{gauss}(0,0,i\zeta)=\int \frac{e^{-(x^2+y^2+z^2)}dx\ dy\ dz}{\sqrt{x^2+y^2+(z-i\zeta)^2}} $$ (dropping the primes) must go as $e^{+\zeta^2}$. Most importantly, I cannot see how to get this behaviour from a numerical integration procedure, which is necessary to get the right behaviour for a real molecule.

For large $\zeta$, the integrand is bounded except for an integrable ring singularity at $x^2+y^2=\zeta^2$, $z=0$, but the charge density goes *down* as a gaussian on that singularity. I do not see a way of getting a super-exponential growth from that integral that does not involve shifting integration contours and thus evaluating $\rho_C$ on complex coordinates (which does not sit well with the real quantum chemical data on my molecule). On the other hand, Liouville's theorem makes me very wary of solutions to that integral that appear to return bounded entire functions for $V_\text {gauss}$.

Can someone help me see a way out of this conundrum, or point me to appropriate references where this or similar issues are treated in some detail?