Analytic continuation of a data set from the upper complex plane to the lower complex plane?

Question

Context

I am interested in identifying damped modes such as those in self gravitating galaxies:

This requires extending to the lower complex plane a dispersion relation which is computed numerically and known over a grid in the upper complex plane.

The (interpolated) complex data in the upper plane looks like this:

Question

How to I extend this function below the real axis?

I.e. how can I sample the analytic continuation of my data set in the lower complex plane?

Toy model

For the sake of constructing a toy model, let me assume that I know numerically the (shifted) Zeta function above the plane as

ComplexPlot[Zeta[z + I], {z, 0, 4 + 2 I}, 
 ColorFunction -> "CyclicReImLogAbs"]

My goal to be able to make a plot like this

ComplexPlot[Zeta[z + I], {z, -2 I, 4 + 2 I}, 
 ColorFunction -> "CyclicReImLogAbs"]

I am told this is possible to some extend (?)

I am aware it is not an easy question in general and that the analytic continuation is only valid on some limited interval which depends on the accuracy of my sampling in the upper plane.

I hope having a genetic tool to do this would be useful to others?

Pade Approximation seems an obvious venue, but I did not understand how to feed it with data rather than functions?

I believe an implementation in fortran is discussed in this paper (Appendix D).

Thank you for your help!

Answer 1

Some colleagues have provided me with a solution to this question which I post for reference. All the credits goes to them.

Let us try it on the Plasma function:

Z[\[Xi]_] = 
 1/Sqrt[Pi] Integrate[
D[Exp[-z^2], z]/(z - \[Xi]) // Evaluate, {z, -Infinity, Infinity}, 
Assumptions -> Im[\[Xi]] > 0];
fz[z_] = 1/(Z[z + I] + 1)

It looks like this:

pl0 = ComplexPlot[fz[z], {z, -6 - 7 I, 6 + 2 I},
PlotRange -> 5 {-1, 1}, PlotPoints -> 5, 
ColorFunction -> "CyclicLogAbsArg", ImageSize -> Medium, 
WorkingPrecision -> 20];

the purpose of the analytic continuation is here to identify all the extrema in the lower complex plane while having only a knowledge of the behaviour of the function on a very finite set of points in the upper complex plane.

---

Let us define a region in the upper complex plane over which the function is sampled.

rg = {{-12, 12}, {0, 1}};
region = 
DiscretizeRegion[(Rectangle @@ Transpose[rg]), 
MaxCellMeasure -> 0.15];

and sample (note the high precision requested in the sampling):

\[Omega]s = 
N[Rationalize[#, 10^-20] & /@ 
MapApply[Complex, (MeshCoordinates@region)], 2500]; n = 
Length@\[Omega]s

Now let us define the actual continuation algorithm:

Clear[\[Rho]];
  \[Rho][i_, k_] := \[Rho][i, k] = 
  If[k == 0, fz[\[Omega]s[[i]]], 
   If[k == 1, (\[Omega]s[[i]] - \[Omega]s[[
        i + 1]])/(fz[\[Omega]s[[i]]] - fz[\[Omega]s[[i + 1]]]), 
    (\[Omega]s[[i]] - \[Omega]s[[i + k]])/(\[Rho][i, k - 1] - \[Rho][
         i + 1, k - 1]) + \[Rho][i + 1, k - 2]
    ]]
   as = Join[{\[Rho][1, 0], \[Rho][1, 1]}, 
   Table[\[Rho][1, k] - \[Rho][1, k - 2], {k, 2, n - 1}]];
   Clear[g];
   g[i_, \[Omega]_] := 
   If[i == n, as[[n]], 
    as[[i]] + (\[Omega] - \[Omega]s[[i]])/g[i + 1, \[Omega]]]
    g[\[Omega]_] := g[1, \[Omega]]

Let us finally compare this analytic continuation to the (known) complex function in the lower half plane:

pl1 = ComplexPlot[g[z], {z, -6 - 7 I, 6 + 2 I}, PlotPoints -> 5, 
   PlotRange -> 3 {-1, 1}, AspectRatio -> Automatic, 
   ColorFunction -> "CyclicLogAbsArg", ImageSize -> Medium];
pl1 = Show[pl1, Graphics@(Point[ReIm[#]] & /@ \[Omega]s), 
   Graphics@InfiniteLine[{-1, 0}, {1, 0}]
   ];
{pl0, Show[pl1, PlotRange -> {{-6, 6}, Automatic}]}

Given that it only knows about the few black dots at the top, it's doing a rather good extrapolation job!