How to calculate infinitesimal analytic continuation?

Question

Many problems in science and engineering are related to the analytic continuation and in particular infinitesimal analytic continuation to the upper or lower complex plane, i.e., a generic complex function $f(\omega)$ with real $\omega$ is changed to $$F(\omega)=\lim_{\eta\rightarrow0^+}f(\omega\pm i\eta).$$ In many applications, only the part introduced by $\eta$ is necessary, which is $$A(\omega)=\lim_{\eta\rightarrow0^+}[f(\omega+i\eta)-f(\omega-i\eta)].$$

Question: How to symbolically obtain $A(\omega)$?

A related question here deals with the simplest case (the first formula below), but I am asking about more general cases when symbolic Fourier transform or integral does not work.

Simple Limit does not work, because the Sokhotski–Plemelj formula $$\lim_{\eta\to0^{+}} \frac{1}{x\pm i\eta}= \mp i\pi\delta(x) + {\mathcal{P}} {\Big(\frac{1}{x}\Big)}$$ plays a role in the calculation. Here $\delta$ denotes Dirac delta function (DiracDelta), $\mathcal{P}$ denotes the Cauchy principal value. A more general version of this formula is $$\lim_{\eta\to0^{+}} \frac{1}{(x\pm i\eta)^{n+1}}= \mp i\pi(-1)^n\frac{\delta^{(n)}(x)}{n!} + {\mathcal{P}} {\Big(\frac{1}{x^{n+1}}\Big)}$$ where $\delta^{(n)}$ denotes the $n$-th derivative of $\delta$-function. But note that all these will just appear as convenient symbols, not to be evaluated anywhere.

---

Two examples are given below, where all parameters are real. Therefore only the imaginary part is purely generated by $\eta$. For $f_1(\omega)=\frac{a\omega+\sqrt{c+\omega^2}}{e}$, we have \begin{equation} \begin{split} A(\omega)&=2i\lim_{\eta\rightarrow0^+}\Im f(\omega+i\eta)\\ &=2i\lim_{\eta\rightarrow0^+}\Im\frac{a(\omega+i\eta)+\sqrt{c+(\omega+i\eta)^2}}{e}\\ &\approx2i\lim_{\eta\rightarrow0^+}\Im\frac{a(\omega+i\eta)+c+\omega^2+i\eta\frac{\omega}{c+\omega^2}}{e}\\ &=2i\lim_{\eta\rightarrow0^+}\frac{ai\eta+i\eta\frac{\omega}{c+\omega^2}}{e}\\ &=0. \end{split} \end{equation} For $f_2(\omega)=\frac{a+b\omega}{\omega-c}$ we have \begin{equation} \begin{split} A(\omega)&=\lim_{\eta\rightarrow0^+} f(\omega+i\eta)-f(\omega-i\eta)\\ &=\lim_{\eta\rightarrow0^+}(a + b c)(\frac{1}{\omega-c+i\eta}-\frac{1}{\omega-c-i\eta})\\ &=[-i\pi\delta(\omega-c)+\mathcal{P}(\frac{1}{\omega-c})]-[+i\pi\delta(\omega-c)+\mathcal{P}(\frac{1}{\omega-c})]\\ &=-2i\pi(a+bc)\delta(\omega-c). \end{split} \end{equation} Here, simple Limit does not work, which merely gives 0 as shown below.

f[ω_] := (a + b ω)/(ω - c);
Limit[f[ω + I η] - f[ω - I η], η -> 0, Direction -> "FromAbove"]

Edit

We aim at generic functions $f(\omega)$ that can be more complicated than the above minimal examples, but let's assume the denominator, if x is involved, can always be factorized. For instance (use x as the variable),

f[x_] := (a - b x)/(c^3 + 3 c^2 d + 3 c d^2 + d^3) + Sqrt[x + b]/(
   x - e) + (
   b x + Sqrt[x + b^2])/((x - c)^2 (x - Sqrt[d + c^2]) (x^2 + g^2));

In the result, there should be $\delta(x-e),\delta(x-\sqrt{d+c^2}),\delta(x-c),\delta^{(1)}(x-c)$. Also similar

f[x_] := (
       b x + Sqrt[x + b^2])/((x - c)^2 (x - Sqrt[d + c^2]+a b^2) (x^2 + g^2))

which should give $\delta(x-\sqrt{d+c^2}+a b^2),\delta(x-c),\delta^{(1)}(x-c)$.

I thought about this for some time, but due to my very limited MMA skill, I don't see how to realize such a calculation.

Answer 1

The manual transform can be automated with pattern matching in a straightforward manner:

generalLimit[expr_, η_ -> 0] := 
 Limit[expr /. 
   a_. (x_ + c_. Complex[0, i_] η)^n_ :> 
    a (-Complex[0, Sign@i] Pi (-1)^(-1 - n)
         Derivative[-1 - n][δ][x]/(-1 - n)! + \[ScriptCapitalP][
        1/x^-n]), η -> 0]

Since you've mentioned that "all these will just appear as convenient symbols, not to be evaluated anywhere", I use δ instead of DiracDelta so Simplify, etc. won't have any effect on it.

I dare not say the generalLimit will handle every case, but it does work on the examples in the question:

generalLimit[1/(x + I η)^(n + 1), η -> 0]

generalLimit[1/(x + I η), η -> 0]

generalLimit[ccc/(x + w - 2 y I η), η -> 0]

f[ω_] := (a + b ω)/(ω - c);

generalLimit[
  f[ω + I η] - f[ω - I η], η -> 0] // Simplify

And once you've fully understood the definition of generalLimit, you should know how to patch it for more complicated cases. To learn more about pattern matching, consider reading posts under semanic-matching, see also tips and tricks of Ted Ersek (Notebook edition is here).