# Perturbation theory with Mathematica: Definite integral of polynomial times exponential times hypergeometric function of imaginary argument

I would like to ask also Mathematica users about my question from the math forum. To expand, I'm adding the code which calculates the full double integral for $n=0$ and $\mu=0$ (the second in the post):

f[r_] := r^n E^(-r (1 + μ)) E^(-I k r) Hypergeometric1F1[I/k+1,2,2 I k r]; g = f[r] /. {n -> 0, μ -> 0};

NIntegrate[k^2 Exp[π/k] Abs[Gamma[1 - I/k]]^2/(1 + k^2)Abs[Integrate[g, {r, 0, ∞}]]^2, {k, 0, ∞}]


giving 1.65719 which confirms the result which I obtained by another (perturbative) method.

However, as Jinxed pointed out, if one only evaluates the inner integral (over $r$), it does not converge, if one assumes $k$ to be real:

f[r_] := r^n E^(-r (1 + μ)) E^(-I k r) Hypergeometric1F1[I/k + 1, 2, 2 I k r]; g = f[r] /. {n -> 0, μ -> 0};Integrate[g, {r, 0, ∞}, Assumptions -> {n == 0 || n == 1, μ ∈ Reals, μ >= 0, n ∈ Integers, k > 0, k ∈ Reals}]


Could you help me make sense of this? I appreciate any ideas!

• Have you actually tried solving it in Mathematica ? Mar 25 '15 at 15:44
• He says in the linked question, "Mathematica can handle this integral in about one minute if numerical values of $n$ and $\mu$ are chosen." It would be nice if you could give us your MMA code though. Mar 25 '15 at 16:53
• With f[r_] := r^n E^(-r (1 + \[Mu])) E^(-I k r) Hypergeometric1F1[I/k + 1, 2, 2 I k r]; g = f[r] /. {n -> 0, \[Mu] -> 0};Integrate[g, {r, 0, \[Infinity]}, Assumptions -> {n == 0 || n == 1, \[Mu] \[Element] Reals, \[Mu] >= 0, n \[Element] Integers, k > 0, k \[Element] Reals}] it sure does not converge. Mar 25 '15 at 18:06
• For the assumptions given in the reference, Mathematica claims that the "Integrate does not converge except for mu = -1, and indeed gives error messages for, say, mu = 2. But this makes no sense, because the integrand is smaller at large r for larger mu. Moreover, NIntegrate does converge for numerical values of k and n. Perhaps, Mathematica is mishandling this integral. Mar 25 '15 at 18:55
• I just performed the integral numerically for k = .01 and, even though it did not converge well, it gave the answer 0.432328, which is close to 1/2 - 1/(2 E^2). I have run a few other cases and obtained good agreement too. Mar 25 '15 at 23:44

Here is a rigorous way to deal with the integration problem, making use only of the assumption that the integral is convergent so that we can exchange the integration and the series expansion for the Coulomb function:

f[r_] :=
r^n E^(-r (1 + μ)) E^(-I k r) Hypergeometric1F1[I/k + 1, 2,
2 I k r];
g = f[r] /. {n -> 0, μ -> 0}

(* ==> E^(-r - I k r) Hypergeometric1F1[1 + I/k, 2, 2 I k r] *)


This is the definition you posted. The trouble is caused by the hypergeometric function, so let's simply integrate term by term in its series expansion:

h =
g /. hype_Hypergeometric1F1 :> SeriesCoefficient[hype, {r, 0, ν},
Assumptions -> ν >= 0] r^ν

(*
==> (2^ν E^(-r -
I k r) (I k)^ν r^ν (I/k + ν)!)/((I/
k)! ν! (1 + ν)!)
*)

hype[ν_] =
Integrate[h, {r, 0, Infinity}, Assumptions -> k > 0 && ν >= 0]

(*
==> -((2^ν (k/(-I + k))^(1 + ν) Gamma[1 + I/k + ν])/(
Gamma[I/k] Gamma[2 + ν]))
*)

FullSimplify[Sum[hype[ν], {ν, 0, Infinity}], k > 0]

(* ==> 1/2 - 1/2 ((I + k)/(I - k))^(-(I/k)) *)


This is the exact result. In this way of doing the calculation, Mathematica has no difficulties with the integration.

You can even do this for arbitrary parameters:

h =
f[r] /. hype_Hypergeometric1F1 :>
SeriesCoefficient[hype, {r, 0, ν},
Assumptions -> ν >= 0] r^ν

(*
==> (2^ν E^(-I k r - r (1 + μ)) (I k)^ν r^(
n + ν) (I/k + ν)!)/((I/k)! ν! (1 + ν)!)
*)

hype[ν_] =
Integrate[h, {r, 0, Infinity},
Assumptions -> k > 0 && ν >= 0 && μ >= 0 && n >= 0]

(*
==> ((2 I)^ν k^ν (1 + I k + μ)^(-1 -
n - ν) (I/k + ν)! Gamma[1 + n + ν])/((I/
k)! ν! (1 + ν)!)
*)

FullSimplify[Sum[hype[ν], {ν, 0, Infinity}], k > 0]

(*
==> (1 + I k + μ)^(-1 - n)
Gamma[1 + n] Hypergeometric2F1[(I + k)/k, 1 + n, 2, (2 k)/(
k - I (1 + μ))]
*)


This is the general, exact solution for arbitrary values of the physical parameters $n$, $\mu$ and $k$.

Edit to add some more explanation

My goal was to let h represent a single general term in the original integrand, obtained by series expanding the problematic part of that integrand, which is Hypergeometric1F1.

Instead of re-typing this or copying and pasting, I use ReplaceAll (/.) with RuleDelayed to insert a general series term of the hypergeometric function in place of the function itself. For such replacements, one has to specify a pattern that is to be replaced. For purposes of pattern matching, Hypergeometric1F1 is just the Head of an expression whose content is the argument list. To specify a pattern with a desired Head, one uses the _ symbol (Blank) preceding the desired Head. In our case, the replacement is not static in the sense that what we need to insert in place of the hypergeometric function depends on its arguments. So I need to take the expression to be replaced, and first do some manipulation on it (expand it in a series). That means I need a way to refer to this expression, and that's what the hype is for: it's just a name I made up that stands for the expression matching the pattern.

In terms of computer algebra, pattern matching and replacement is the central engine behind Mathematica, and without it we could just as well go back to tables of integrals. Mathematically, what h represents can be illustrated perhaps in one more step by doing (after executing the previous code)

Sum[h, {ν, 0, Infinity}]

(*
==> E^(-I k r - r (1 + μ)) r^n Hypergeometric1F1[1 + I/k, 2,
2 I k r]
*)

% == f[r]

(* ==> True *)


which means that h is indeed represents the original integrand itself, if you sum over the index $\nu$ I introduced as an auxiliary variable above.

• +1 Excellent, thanks so much. As I am a beginning user, could you please explain me the structure of the part of the code where you define h? The problem for me is caused by /. hype_Hypergeometric1F1 :> and how it connects to the surrounding code. Also, what _Hypergeometric1F1 (with the underscore) is exactly used for? Thanks! Mar 27 '15 at 20:33
• P.S. Would it be possible to proceed with the outer integral in a similar fashion? Seems it woudn't, but maybe first I have to understand the hype_Hypergeometric1F1. Mar 27 '15 at 21:18
• I just edited the answer to explain what I did in more detail.
– Jens
Mar 27 '15 at 21:31
• I don't think you can easily use the analogous approach for the $k$ integral because the integrands there involve absolute squares. But in any case, that goes beyond what the original question contains. Since you're doing perturbation theory, it seems reasonable to approximate the $k$ integral first before trying anything with the exact form. Maybe I'll have an idea later...
– Jens
Mar 27 '15 at 21:40
• Oh yes, I forgot that now the formulation of the question is slightly different than it was before my update (the outer integral was the second half of the problem). This means that I can accept your answer and later I can post a new question regarding the full integral - actually the result of the full double integral will be than plugged into an infinite sum, which is the ultimate goal. Mar 27 '15 at 22:04

Define

arg = r^n Exp[-r (1 + mu + I k)] Hypergeometric1F1[1 + I/k, 2, 2 I k r]


Then, performing the inner integral in the Question, complete with Assumptions,

Integrate[arg, {r, 0, Infinity}, Assumptions -> {k > 0, mu ∈ Reals, n ∈ Integers,
n >= 0}]


yields a ConditionalExpression with highly restrictive conditions, n < 1 && mu == -1

(*ConditionalExpression[(Sqrt[Pi]*((Gamma[(-I)/k + n]*
Hypergeometric2F1[I/k, (I + k)/k, 1 + I/k - n, 1/2])/
(8^(I/k)*Gamma[3/2 + I/k]*Gamma[-((2*I + k)/k)]) +
(Gamma[I/k - n]*Gamma[n]*Hypergeometric2F1[n, 1 + n, 1 - I/k + n, 1/2])/
(8^n*Gamma[I/k]*Gamma[-1 - 2*n]*Gamma[3/2 + n])))/(8*(I*k)^n),
n < 1 && mu == -1] *)


Moreover, if mu ∈ Reals is replaced by mu >= 0 in the Assumptions, Mathematica 10.0.2.0 asserts that the integral does not converge:

Integrate::idiv: Integral of E^(-(1+I k+mu) r) r^n Hypergeometric1F1[(I+k)/k,2,2 I k r] does not converge on {0,∞}. >>


This assertion seems dubious, because arg evaluated at mu == -1 bounds arg evaluated at larger values of mu for r > 0.

An alternative approach is to perform the integral without Assumptions

Integrate[arg, {r, 0, Infinity}]


which yields a ConditionalExpression with conditions such as Im[k] > 0 that ostensibly violate the constraints imposed by the Question.

(* ConditionalExpression[((-I)*k)^(-1 - n)*(-(k/(k - I*(1 + mu))))^(1 + n)*Gamma[1 + n]*
Hypergeometric2F1[(I + k)/k, 1 + n, 2, (2*k)/(k - I*(1 + mu))],
Im[k] > 0 && 1 + Re[n] > 0 && Im[k] <= 1 + Re[mu] &&
(Im[k^(-1)] + Re[n] < 1 || Im[k] < 1 + Re[mu])] *)


So, ignore those conditions and simply assume

ans = %[[1]]
(* ((-I)*k)^(-1 - n)*(-(k/(k - I*(1 + mu))))^(1 + n)*Gamma[1 + n]*
Hypergeometric2F1[(I + k)/k, 1 + n, 2, (2*k)/(k - I*(1 + mu))] *)


This assumption is readily tested by comparing ans with results of NIntegrate for relevant values of n, mu, and k.

Show[Plot[Evaluate[Chop@Table[ans, {n, 0, 1}, {mu, 1, 2}]], {k, .01, 5},
PlotRange -> {0, .4}, ImageSize -> 400, PlotLegends ->
Placed[{"n=0, mu=1", "n=0, mu=2", "n=1, mu=1", "n=1, mu=2"},
Scaled[{{1, 1}, {1, 1}}]]],
ListPlot[Flatten[Table[Chop[NIntegrate[arg, {r, 0, \[Infinity]}], 10^-5],
{n, 0, 1}, {mu, 1, 2}, {k, {.01, 1, 2, 3, 4, 5}}], 1] // Quiet]]


Evidently, ans is a valid answer even for Im[k] = 0, and the assertion by Mathematica cited above that the integral does not converge for mu >= 0 is inaccurate.

For completeness, ans can be simplified for

Table[Simplify[ans], {n, 0, 1}, {mu, 1, 2}]
(* {{1/2 - 1/(2*((2*I + k)/(2*I - k))^(I/k)),
1/2 - 1/(2*((3*I + k)/(3*I - k))^(I/k))},
{1/(((2*I + k)/(2*I - k))^(I/k)*(4 + k^2)),
1/(((3*I + k)/(3*I - k))^(I/k)*(9 + k^2))}} *)


and for larger integer values of n and mu as well.

• +1 Thank you for the nice analysis. Maybe someone on the math forum suggests what could be the reason for such a behavior of the integral (and of Mathematica's treatment of it). Mar 26 '15 at 15:04
• I have another observation, please correct me if I am wrong (I restarted Mathematica to be sure). If I make the replacement $\mu\to\mu-1$ in the integral, then the resulting (equivalent, if the assumptions are properly changed) code does not give the condition Im[k]>0, only \[Mu] > Im[k]. The code: f[r_] := r^n E^(- \[Mu] r) E^(-I k r) Hypergeometric1F1[I/k + 1, 2, 2 I k r]; g = f[r]; Integrate[g, {r, 0, \[Infinity]}, following Assumptions -> {\[Mu] \[Element] Reals, \[Mu] >= 1, n \[Element] Integers, n >= 0}] Could you please check? Did I do any mistake? Mar 26 '15 at 15:13
• @wondering I copied and ran the code in your last comment but still get the usual conditions,  [Mu] >= Im[k] && Im[k] > 0 && ([Mu] != Im[k] || n + Im[k^(-1)] < 1). Mar 26 '15 at 16:12
• I wonder what causes this behavior. I tried a couple of times. Completely restarting Mathematica 10.0.0.0 gives \[Mu] > Im[k] with the above code (with the word "following":-), however, if I restart and add Clear["Global*"]; before the code, I get what you say. Strange. Mar 27 '15 at 20:12
• @wondering That is strange, and I have not good answer. Sorry. Mar 27 '15 at 20:19