# Numerically Recovering Imaginary Part of Analytic Continuation from Real Part

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$.

• Is numerical analytic continuation definitely helpful here? Could you instead evaluate it along $\mathbb{R}+0.99\mathrm{i}$, where the integrand of $\bar f$ would be exponentially decaying but $f$ is available directly? It is also somewhat surprising that the oscillatory integral could not be evaluated, because usually the specialized methods for oscillatory integrals would be able to handle a quickly polynomially-decaying function like $f$. I say this because after playing with the integral in Mathematica I'm concerned that numerical analytic continuation might happen to be a blind alley. – Kirill Jun 22 '18 at 12:41
• I tried to implement the integral directly in Mathematica, and I got it to evaluate $\bar f$ in 20s on my laptop: not great, but also not hours. Would it help you if I wrote up an answer for how to evaluate it directly? – Kirill Jun 22 '18 at 16:57
• @Kirill After many failed attempts at numerical analytic continuation I am inclined to wholeheartedly agree with your first comment. Please, if you were able to evaluate $\bar{f}$ in 20s, I would be very grateful for a write up. By the way, it might help to add that actually it wasn't just one single evaluation of $\bar{f}$ that was taking hours, but approximately 30 evaluations (evaluating for $\alpha\in [-1,2]$ with interval size $0.1$). A single evaluation though was taking me about 14 minutes. – Arturo don Juan Jun 22 '18 at 17:05
• I wrote it up, but I discovered an issue with my code, so I'm no longer sure if what I computed is at all valid. Do you have any known-valid reference values? – Kirill Jun 23 '18 at 12:23

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.

• Thanks, I hadn't really thought about fiddling with the NIntegrate options/strategies. Oh and btw the integral over $k_4$ goes from $-\infty$ to $\infty$, not from $0$ to $\infty$. The result should be a symmetric/even function, and so there should be no imaginary part in the Fourier Transform. Also, I don't really understand why you defined g[t,e0] the way you did (at the end of your code). What's all that stuff starting with the first "EvaluationMonitor" for? – Arturo don Juan Jun 24 '18 at 20:27
• @ArturodonJuan Are you sure $f$ is symmetric in $E$? The powers in the numerator are different and $p_1$ and $p_2$ are not interchangeable under $E\leftrightarrow -E$. Your function doesn't look symmetric to me, unless I did something wrong. I changed the integral over $k_4$ to $2\times$ the integral over $(0,\infty)$ for simplicity. EvaluationMonitor is left over from debugging it. – Kirill Jun 24 '18 at 21:22
• (1) Oh yeah, sorry you're right. My actual expression/integral of interest is the one I put up plus $p_1\leftrightarrow p_2$, which makes it symmetric in $E$. I should have clarified that in my original post. (2) I made a mistake in my definition of $p_{1,2}$. There should only be one power of $k_4$ sitting in those parentheses, not two. I just changed that. – Arturo don Juan Jun 24 '18 at 23:22
• @ArturodonJuan I think it makes no real difference to how the answer works, only the numbers would change. – Kirill Jun 25 '18 at 7:38