Numerically Recovering Imaginary Part of Analytic Continuation from Real Part

Question

My situation.

I have a function of a complex variable $f(z)$ defined through a complicated integral. What I am interested in is the value of this function on the imaginary axis. I have numerical access to this function on the following ribbon: $z=(x,y)\in (-\infty,\infty)\times[-1,1]$. Formally the integral expression is divergent outside this domain, and therefore I need an analytic continuation. To sum up my situation in a picture,

Here is what I know about $f(z)$ on this ribbon from numerics:

1. It is simultaneously symmetric about the imaginary and real axes.

2. It decays to zero at $Re(z)\rightarrow\infty$.

3. It blows up near $z=\pm i$. It could be pole or a branch point, I don't know. I suspect the nature of this singularity (and possibly all other isolated singularities of the analytic continuation) depends on the specific parameterization $\xi$ of this function (see integral below for details)

In fact it looks very similar to a $\text{sech}^2(z)$ or a $1/(1+z^2)^{2n}$ when plotted. Here is a plot of the real part:

My question is, given the sheer amount of information I have about the function (total numerical access to it on that ribbon), is there some way for me to numerically calculate an approximation to this function along the imaginary axis? I am using Mathematica by the way.

The reason I am interested in the values along the imaginary axis is because I need to evaluate the following Fourier transform of this function:

$$\bar{f}(t)=\int_{-\infty}^{\infty} dx\, e^{itx}\frac{1}{x^2+x^2_0} f(x) \tag{1}$$

for large values of $t$, which in my case is actually on the order of $10$. Although I know the integrand well, this Fourier transform is formidably oscillatory, so the only other way I know how to calculate this is by a Contour integration.

---

What I have Tried.

1. I have actually tried to calculate the ultimate highly oscillatory integral, eq. (1). Evaluating eq. (1) for a single value of 't' takes a few hours to compute. I have carried out a few of these integrals already and the results actually make sense, but I would like an alternative approach.

2. I have tried analytically continuing with Pade approximants, but this is also computationally expensive, but not as much as direct evaluation. More importantly, I could not establish convergence with increasing order of the approximants (nor the average of their partial sums!), which is in contrast with how my tests with simple functions like $\text{sech}^2(z)$ went (I could easily get very quick convergence on wide ranges of the complex $z$-plane with simple test functions).

3. I have tried symbolic integration to no avail. I have tried massaging the integrand into a more digestible form for Mathematica, but my attempts have not succeeded.

---

The offending integral.

Let $k_4$,$k_{\perp}$, $\xi$, and $\alpha$ be positive real numbers while $E$ is the complex number we're interested in (plays the role of $z$ in the previous discussion). Define:

$$\begin{align} p_1^2&=\left(k_4+\frac{1}{2}E\right)^2+k_{\perp}^2+\alpha^2\\ p_2^2&=\left(k_4-\frac{1}{2}E\right)^2+k_{\perp}^2+(1-\alpha)^2 \end{align}$$

The integral I'm interested in is the following:

$$\begin{align} f(E;\,\alpha,\xi)&=\int_{-\infty}^{\infty} dk_4 \int_{0}^{\infty} d(k_{\perp}^2)\left[\frac{\alpha (1+p_1^2)^{3\xi/2}(1+p_2^2)^{\xi/2}}{(1+p_1^2(1+p_1^2)^{2\xi})(1+p_2^2(1+p_2^2)^{2\xi})}+\\ \,\,\,\,\,\,\,\,\,\,\,\,\,+(p_1\leftrightarrow p_2)\right] \end{align}$$

where I have suppressed functional dependence notation in the integrand for brevity. I'm particularly interested in the values $\xi=1,2,3$, the range $0<\alpha<1$, and (as stated above) the Fourier transform (1) for $t~10$.

Answer 1

Note: I'm somewhat worried at this point that the integral values Mathematica gives me are bogus. I thought it was working because it gave a sensible-looking result in a short time, but it might be the case that the method it tries to use is buggy or that I did something wrong. So it might be that the code below isn't working at all, I don't know, sorry.

Note 2: It bothered me, so I wrote another version (code here, sorry about the code quality) using Julia and GSL, and it evaluates g in 2 seconds to the same answer that Mathematica gives below. So I think the code is probably okay.

The main reason I was suspicious about numerical analytic continuation is that your integral actually looked kind of nice in my limited testing. In particular, the integrands of both $f$ and $\bar f$ decay polynomially and rather quickly, and this is exactly the sort of thing that conventional quadrature routines are designed to handle well. There are also no tricky singularities.

My past experience with numerical integration leads me to believe that the fancier mathematical methods can sometimes be spectacularly helpful, but also that evaluating numerical Fourier transforms and integrating rational and algebraic functions are the bread-and-butter of numerical integration algorithms, so one can often make easy progress by picking algorithms carefully and playing with their parameters. This is usually the easier option if it's hard to see how to make the mathematical technique work right.

ClearAll[ξ, α, p1, p2, fi, f, g];
ξ = 1;
α = 1/2;
fi[e_, k4_, kp_] := Module[{
   p1 = (k4^2 + e/2)^2 + kp^2 + α^2,
   p2 = (k4^2 - e/2)^2 + kp^2 + (1 - α)^2},
  2 * (* because integrate k4 over (0,∞) *)
   2 kp * (* because d(kp^2) *)
   (α (1 + p1)^(3 ξ/2) (1 + p2)^(ξ/2)) /
     ((1 + p1 (1 + p1)^(2 ξ)) (1 + p2 (1 + p2)^(2 ξ)))
  ]
f[e_?NumericQ] := NIntegrate[
   fi[e, k4, kp], {k4, 0, ∞}, {kp, 0, ∞},
   Method -> {Automatic, "SymbolicProcessing" -> 0}];

(* !!! This gives a bogus result: *)
gBogus[t_?NumericQ, e0_?NumericQ] := 
 NIntegrate[
  Exp[I t e]/(e^2 + e0^2) f[e], {e, -∞, ∞}, 
  Method -> {"DoubleExponentialOscillatory", 
    "SymbolicProcessing" -> 0}]

(* This gives *a* result, different from above despite being equivalent *)
g[t_?NumericQ, e0_?NumericQ] := 
 NIntegrate[Exp[I t e]/(e^2 + e0^2) f[e], {e, -\[Infinity], 0}, 
   Method -> {"DoubleExponentialOscillatory", 
     "SymbolicProcessing" -> 0},
   EvaluationMonitor :> Print["e=", e]] +
  NIntegrate[Exp[I t e]/(e^2 + e0^2) f[e], {e, 0, \[Infinity]}, 
   Method -> {"DoubleExponentialOscillatory", 
     "SymbolicProcessing" -> 0},
   EvaluationMonitor :> Print["e=", e]]

Result:

In[18]:= Timing@g[10,1]
Out[18]= {78.0828, 0.0000704303 + 9.78009*10^-6 I}

In[338]:= Timing@g[1,1]
Out[338]= {14.3125,0.389542 +0.024758 I}

I made Mathematica spend zero time on symbolically preprocessing the integrands, because in this case it wouldn't have been able to figure out anything useful about it anyway. I also told it to specifically use an oscillatory quadrature method for the second integral.

My guess for why randomly fiddling with integration strategies (see NIntegrateIntegrationStrategies) works at all is that sometimes Mathematica might accidentally pick a bad strategy automatically, killing performance, whereas anything I ask it to do is at least a little bit meaningful even if suboptimal. You could also consider getting help at https://mathematica.stackexchange.com, they might know more about Mathematica's internals over there.