How to compute this integral?

Question

I'm trying to calculate this integral :

$$I(z,k,a)= \displaystyle\int_{1}^{\infty} t^2\, \text{ArcTanh} \left(\sqrt{\frac{t^2-1}{t^2}} \dfrac{k}{z}\right)\, e^{-a\,t} \, dt$$

Where :

• $\textrm{ArcTanh(...)}$ : is the inverse hyperbolic tangent function.

• $z=x+i\,y$ : is a complex with the assymption $y<<<x$.

• $k$ and $a$ are reals $>0$.

Mathematica didn't calculate it. Maybe it is too complicated to be done.

In Mathematica input form :

Assuming[{a>0&&k>0},Integrate[t^2 (ArcTanh[Sqrt[(t^2-1)/t^2](k/z)])Exp[-a t],{t,1,Infinity}]]

If I simply enter that into Mathematica, it instantly returns the same expression.

Please, how do I go about this ? Is there any tricks that can be applied ?

Is there a way to get the symbolic result ?

Thank's.

Answer 1

Preliminary post

For k/z = 1 there is a closed form result:

I[1,1,a]=((4 + a^2) BesselK[0, a] + a (4 BesselK[1, a] + a BesselK[2, a]))/(2 a^3) 

The derivation and the extension to k/z != 1 requires some manual interaction to help Mathematica which we will show in the following.

Solution

Summary

As the integral to be calculated is returned unevaluated by Mathematica we show here to get useful partial results by using a method I like to call "man-machine-interaction". This means we carry out a joint work with paper and pencil on one side and Mathematica on the other side.

The case m = k/z == 1 can be solved symbolically. If m < 1 we show how to get the symbolic soluton to aribtrary precision as a power series in (1-m).

Preliminary remarks

1) The integral to be calculated is

fI := Integrate[
  t^2 (ArcTanh[m Sqrt[(t^2 - 1)/t^2]]) Exp[-a t], {t, 1, Infinity}, 
  Assumptions -> {a > 0, 0 < m < 1}]

The numerical integral presents, however, no problem

fN[a_, m_] := 
 NIntegrate[t^2 (ArcTanh[m Sqrt[(t^2 - 1)/t^2]]) Exp[-a t], {t, 1, Infinity}]

2) Here we have set m = k/z, because k and z enter only as a quotient.

3) Also we have introduced natural restrictions on the parameters "a" and "m". The restriction for m comes from the requirement that fI be a real quantity. This in turn requires the argument of the ArcTanh to be between 0 and 1 if we keep the upper limit of integration of t at infinity.

4) The factor t^2 in the integrand can be generated by differentiating twice with respect to "a". Hence it is sufficient to study the simpler intergral without the factor t^2.

Hence the integrand to be studied becomes

fi = Exp[-a t] ArcTanh[m Sqrt[(-1 + t^2)/t^2]];

Step by step solution

The first step will be a partial integration. This step is not done automatically by Mathematica, hence we do it manually.

fi can be written as

fi = fid + fi1

where

fid1 = -(1/a) Exp[-a t] ArcTanh[m Sqrt[(-1 + t^2)/t^2]];

fid = D[fid1, t];

and

fi1 = Simplify[ (E^(-a t) m (2/t - (2 (-1 + t^2))/t^3))/(
  2 a Sqrt[(-1 + t^2)/t^2] (1 - (m^2 (-1 + t^2))/t^2)), {a > 0, t > 0, 
   0 < m < 1}]

(*
Out[99]= (E^(-a t) m)/(a Sqrt[-1 + t^2] (t^2 - m^2 (-1 + t^2)))
*)

Indeed

Simplify[fi == fid + fi1, {a > 0, t > 0, 0 < m < 1}]

(* Out[101]= True *)

Because

fid1 /. t -> 1

(* Out[102]= 0 *)

Limit[fid1, t -> \[Infinity], Assumptions -> a > 0]

(* Out[103]= 0 *)

the integral becomes

fI2 := Integrate[(E^(-a t) m)/(
  a Sqrt[t^2 - 1] (t^2 + m^2 (1 - t^2))), {t, 1, \[Infinity]}, 
  Assumptions -> {a > 0, 0 < m < 1}]

This integral is also returned unevaluated by Mathematica. Hence we look first at a special case.

The case m = 1

If m = 1 the integrand simplifies and the integral can be done expicitly

(E^(-a t) m)/(a Sqrt[t^2 - 1] (t^2 + m^2 (1 - t^2))) /. m -> 1

(* Out[125]= E^(-a t)/(a Sqrt[-1 + t^2]) *)

fI3 = Integrate[%, {t, 1, \[Infinity]}, Assumptions -> a > 0]

(* Out[126]= BesselK[0, a]/a *)

In order to find the original integral we need to differentiate this expression twice

fI4 = D[fI3, {a, 2}] // Simplify

(*
Out[127]= ((4 + a^2) BesselK[0, a] + a (4 BesselK[1, a] + a BesselK[2, a]))/(2 a^3)
*)

The result in the original designation

$$I(1,1,a)=\frac{\left(a^2+4\right) K_0(a)+a (4 K_1(a)+a K_2(a))}{2 a^3}$$

The numerical agreement is excellent:

Plot[fN[a, 1]/fI4, {a, 0.1, 2}, PlotRange -> {{0, 2}, {0.99, 1.01}}, 
 PlotLabel -> 
  "Integral for m = 1\nComparison of symbolical and numerical solution", 
 AxesLabel -> {"a", "fN/fI4"}]

Series expansion

Now we can proceed further developing the integrand of fi1 into a power series about m = 1.

Remark: Although it leads to cumbersome expressions I took four terms in order to possibly find a rule in the result (but, alas, I didn't)

Series[fi1, {m, 1, 4}] // Normal

(* 
Out[122]= E^(-a t)/(a Sqrt[-1 + t^2]) + (E^(-a t) (-1 + m) (-1 + 2 t^2))/(
 a Sqrt[-1 + t^2]) + (E^(-a t) (-1 + m)^2 (1 - 5 t^2 + 4 t^4))/(
 a Sqrt[-1 + t^2]) + (
 E^(-a t) (-1 + m)^4 (-1 + t^2)^(3/2) (1 - 12 t^2 + 16 t^4))/a + (
 E^(-a t) (-1 + m)^3 (-1 + 9 t^2 - 16 t^4 + 8 t^6))/(a Sqrt[-1 + t^2])
*)

fis = List @@ %;

Timing[Integrate[fis, {t, 1, \[Infinity]}, Assumptions -> {a > 0}]]

(*
{1233.4843069`, {BesselK[0, a]/a, ((-1 + m) BesselK[2, a])/a, (
  3 (-1 + m)^2 (4 a BesselK[0, a] + (8 + a^2) BesselK[1, a]))/a^4, (
  15 (-1 + m)^4 (a (112 + a^2) BesselK[2, a] + 
     4 (168 + 5 a^2) BesselK[3, a]))/
  a^6, ((-1 + 
     m)^3 (24 a (20 + a^2) BesselK[0, a] + (960 + 168 a^2 + a^4) BesselK[1, 
       a]))/a^6}}
*)

Finally, we need the second derivative with respect to a:

D[%[[2]], {a, 2}] // Simplify

(*
{((4 + a^2) BesselK[0, a] + a (4 BesselK[1, a] + a BesselK[2, a]))/(
 2 a^3), ((-1 + m) (a^2 BesselK[0, a] + 4 a BesselK[1, a] + 
    8 BesselK[2, a] + 2 a^2 BesselK[2, a] + 4 a BesselK[3, a] + 
    a^2 BesselK[4, a]))/(4 a^3), (
 3 (-1 + m)^2 (16 a (20 + a^2) BesselK[0, 
      a] + (640 + 144 a^2 + 3 a^4) BesselK[1, a] + 
    a (8 + a^2) (16 BesselK[2, a] + a BesselK[3, a])))/(4 a^6), (1/(
 4 a^8))15 (-1 + m)^4 (a^3 (112 + a^2) BesselK[0, a] + 
    32 a^2 (91 + a^2) BesselK[1, a] + 29568 a BesselK[2, a] + 
    592 a^3 BesselK[2, a] + 2 a^5 BesselK[2, a] + 
    112896 BesselK[3, a] + 5184 a^2 BesselK[3, a] + 
    52 a^4 BesselK[3, a] + 16128 a BesselK[4, a] + 
    432 a^3 BesselK[4, a] + a^5 BesselK[4, a] + 
    672 a^2 BesselK[5, a] + 20 a^4 BesselK[5, a]), (1/(
 4 a^8))(-1 + m)^3 (8 a (10080 + 600 a^2 + 7 a^4) BesselK[0, a] + 
    3 (53760 + 11840 a^2 + 368 a^4 + a^6) BesselK[1, a] + 
    a (8 (2880 + 456 a^2 + 7 a^4) BesselK[2, a] + 
       a (960 + 168 a^2 + a^4) BesselK[3, a]))}
*)

Observation

Related simpler example where Mathematica needs some help

This integral is returned unevaluated

Integrate[Exp[-a Cosh[u]], {u, 0, \[Infinity]}, Assumptions -> a > 0]

(*
Out[124]= Integrate[E^(-a Cosh[u]), {u, 0, \[Infinity]}, Assumptions -> a > 0] *)

But with the substitution Cosh[u] -> t we obtain an eqivalent expression of the integral which now is deone by Mathematica:

Integrate[Exp[-a t]/Sqrt[t^2 - 1], {t, 1, \[Infinity]}, Assumptions -> a > 0]

(* Out[142]= BesselK[0, a] *)

Answer 2

For k/z = 1, and integrating by parts:

ClearAll["Global`*"]
Iv[t_] := Exp[-a*t]*t^2;
u[t_] := ArcTanh[Sqrt[(t^2 - 1)/t^2]];
v = Integrate[Iv[t], t];
Du = Simplify@D[u[t], t];
Int == u[t]*v - Integrate[Du*v, t]

HoldForm[Integrate[Exp[-a*t]*t^2*ArcTanh[Sqrt[(t^2 - 1)/t^2]], {t, 1, Infinity}] == 
Limit[-((E^(-a t) (2 + 2 a t + a^2 t^2) ArcTanh[Sqrt[(-1 + t^2)/
   t^2]])/a^3), t -> Infinity] - 
Limit[-((
 E^(-a t) (2 + 2 a t + a^2 t^2) ArcTanh[Sqrt[(-1 + t^2)/t^2]])/
 a^3), t -> 1] + 
Integrate[(E^(-a t) (2 + 2 a t + a^2 t^2))/(
Sqrt[1 - 1/t^2] t), {t, 1, Infinity}]/a^3]

$\int_1^{\infty } \exp (-a t) t^2 \tanh ^{-1}\left(\sqrt{\frac{t^2-1}{t^2}}\right) \, dt=\lim_{t\to \infty } \, -\frac{e^{-a t} \left(2+2 a t+a^2 t^2\right) \tanh ^{-1}\left(\sqrt{\frac{-1+t^2}{t^2}}\right)}{a^3}-\lim_{t\to 1} \, -\frac{e^{-a t} \left(2+2 a t+a^2 t^2\right) \tanh ^{-1}\left(\sqrt{\frac{-1+t^2}{t^2}}\right)}{a^3}+\frac{\int_1^{\infty } \frac{e^{-a t} \left(2+2 a t+a^2 t^2\right)}{\sqrt{1-\frac{1}{t^2}} t} \, dt}{a^3}$

Limit[u[t]*v, t -> Infinity, Assumptions -> a > 0] - 
Limit[u[t]*v, t -> 1, Assumptions -> a > 0]

$0$

ExpandAll@Simplify[Du*v, t > 0]

in1 = Assuming[{a > 0, t > 0}, Integrate[-((2 E^(-a t))/(a^3 Sqrt[-1 + t^2])), {t, 1, Infinity}]];
in2 = Assuming[{a > 0, t > 0}, Integrate[-((2 E^(-a t) t)/(a^2 Sqrt[-1 + t^2])), {t, 1, 
 Infinity}]];
in3 = Assuming[{a > 0, t > 0}, Integrate[(-E^(-a t) t^2)/(a Sqrt[-1 + t^2]), {t, 1, Infinity}]];
HoldForm[Integrate[Exp[-a*t]*t^2*ArcTanh[Sqrt[(t^2 - 1)/t^2]], {t, 1, 
Infinity}]] == -(in1 + in2 + in3)

$\int_1^{\infty } \exp (-a t) t^2 \tanh ^{-1}\left(\sqrt{\frac{t^2-1}{t^2}}\right) \, dt=\frac{2 K_0(a)}{a^3}+\frac{2 K_1(a)}{a^2}+\frac{a K_0(a)+K_1(a)}{a^2}$

where: $K_0(a)$ and $K_1(a)$ the modified Bessel function of the second kind.