Yet another integral identity involving a Bessel function and incomplete gamma functions.

Question

Let $a,b > 0 $ and $T > 0$ and $\alpha= 2$. then by combining my answer to this question with this result we noticed that the following identity below holds true.

\begin{eqnarray} &&rhs_{a,b}(T):=\\ &&\underline{\int\limits_0^a \frac{1}{\alpha} (a-s)^{\frac{1}{\alpha}} \Gamma[-\frac{1}{\alpha},(a-s) T^{-\alpha}] \sqrt{b} \frac{I_1(2 \sqrt{b s})}{\sqrt{s}} ds}+\\ &&\int\limits_0^b \frac{1}{\alpha} (b-s)^{-\frac{1}{\alpha}} \gamma[+\frac{1}{\alpha},(b-s) T^{\alpha}] \sqrt{a} \frac{I_1(2 \sqrt{a s})}{\sqrt{s}} ds + \\ &&-\underline{\underline{T I_0(2 \sqrt{a b})}} + \underline{\frac{1}{\alpha} a^{\frac{1}{\alpha}} \Gamma(-\frac{1}{\alpha},a T^{-\alpha})} + \frac{1}{\alpha} b^{-\frac{1}{\alpha}} \gamma(+\frac{1}{\alpha},b T^{\alpha}) = \\ &&\frac{\sqrt{\pi } e^{-2 \sqrt{a b}}}{2 \sqrt{b}}-\frac{\sqrt{\pi } \left(e^{-2 \sqrt{a b}} \text{erfc}\left(\sqrt{b} T-\frac{\sqrt{a}}{T}\right)+e^{2 \sqrt{a b}} \text{erfc}\left(\frac{\sqrt{a}}{T}+\sqrt{b} T\right)\right)}{4 \sqrt{b}}\tag{1} \end{eqnarray}

In[6892]:= {a, b, alpha, T} = RandomReal[{0, 5}, 4]; alpha = 2;
NIntegrate[
  1/alpha (a - s)^(1/alpha) Gamma[-1/
     alpha, (a - s) T^(-alpha)] (b BesselI[1, 2 Sqrt[b s]])/
    Sqrt[b s], {s, 0, a}] + 
 NIntegrate[
  1/alpha (b - s)^(-1/alpha) Gamma[1/alpha, 
    0, (b - s) T^(alpha)] (a BesselI[1, 2 Sqrt[a s]])/Sqrt[a s], {s, 
   0, b}] + -T BesselI[0, 2 Sqrt[b (a)]] + 
 1/alpha a^(1/alpha) Gamma[-1/alpha, a T^(-alpha)] + 
 1/alpha b^(-1/alpha) Gamma[1/alpha, 0, b T^(alpha)]
Sqrt[Pi]/(2 Sqrt[b]) Exp[-2 Sqrt[a b]] - 
 Sqrt[Pi]/(4 Sqrt[b]) (Exp[2 Sqrt[a b]] Erfc[
      Sqrt[b] T + Sqrt[a] 1/T] + 
    Exp[-2 Sqrt[a b]] Erfc[Sqrt[b] T - Sqrt[a] 1/T])
Out[6893]= 1.29874
Out[6894]= 1.29874

My question is twofold. Firstly how would one go about proving $(1)$ in a different way, i.e. without using the approach in the links above. secondly, how does the right hand side look like when $\alpha \neq 2$?

Update:

Note that when $T \rightarrow \infty$ then the identity holds true as I am showing below.

\begin{eqnarray} \lim\limits_{T \rightarrow \infty} rhs_{a,b}(T) &=& \frac{\pi b^{-1/\alpha } \csc \left(\frac{\pi }{\alpha }\right) \left(\, _0\tilde{F}_1\left(;\frac{\alpha -1}{\alpha };a b\right)-a^{\frac{1}{\alpha }} b^{\frac{1}{\alpha }} \, _0\tilde{F}_1\left(;1+\frac{1}{\alpha };a b\right)\right)}{\alpha } \\ &=&\frac{\pi \left(\frac{\cosh \left(2 \sqrt{a} \sqrt{b}\right)}{\sqrt{\pi }}-\frac{\sinh \left(2 \sqrt{a} \sqrt{b}\right)}{\sqrt{\pi }}\right)}{2 \sqrt{b}} \\ &=& \frac{\sqrt{\pi } e^{-2 \sqrt{a b}}}{2 \sqrt{b}} \end{eqnarray}

Here the first step follows from the Update in my answer to this question.

Update 1:

Let us take the limit of $ b \rightarrow 0 $ from both sides. Then in the right hand side the two integrals vanish and out of the three terms in the bottom line the first and the third go into $-T$ and $+T$ respectively and we have: We have:

\begin{eqnarray} \lim\limits_{b\rightarrow 0} rhs_{a,b}(T) &=& \underline{\frac{1}{\alpha} a^{\frac{1}{\alpha}} \Gamma(-\frac{1}{\alpha},a T^{-\alpha})} \\ &=&T e^{-\frac{a}{T^2}}-\sqrt{\pi } \sqrt{a} \text{erfc}\left(\frac{\sqrt{a}}{T}\right) \\ &=&\lim\limits_{b\rightarrow 0} lhs_{a,b}(T) \end{eqnarray}

Likewise, let us take the limit of $a \rightarrow 0$. Then, again, the two integrals on the right hand side vanish and out of the the three terms in the bottom line the first two cancel each other and we have:

\begin{eqnarray} \lim\limits_{a\rightarrow 0} rhs_{a,b}(T) &=& \frac{1}{\alpha} b^{-\frac{1}{\alpha}} \gamma(+\frac{1}{\alpha},b T^{\alpha}) \\ &=& \frac{\sqrt{\pi }-\sqrt{\pi } \text{erfc}\left(\sqrt{b} T\right)}{2 \sqrt{b}} \\ &=&\lim\limits_{a\rightarrow 0} lhs_{a,b}(T) \end{eqnarray}

Update 2:

Define ${\hat rhs}_{a,m}(T)$ and ${\hat lhs}_{a,m}(T)$ as the coefficient of the right and the left hand sides at $b^m$. Then we have:

\begin{eqnarray} &&{\hat rhs}_{a,m}(T) := \\ &&\frac{1}{m!}\sum\limits_{m_1=0}^m \left(\frac{(-1)^{m_1} a^{m-m_1} T^{\alpha m_1+1}}{\left(m-m_1\right)! \left(\alpha m_1+1\right) } \right)+\frac{\frac{\frac{1}{\alpha }! (m-1)! \Gamma \left(-\frac{1}{\alpha }\right) a^{\frac{1}{\alpha }+m}}{\alpha \left(\frac{1}{\alpha }+m\right)!}+\frac{T a^m \, _2F_2\left(1,-\frac{1}{\alpha };1-\frac{1}{\alpha },m+1;-a T^{-\alpha }\right)}{m}}{(m-1)! m!}-\frac{T a^m}{(m!)^2} \\ &&= \frac{1}{\alpha m!} \exp\left( -\frac{a}{T^\alpha} \right) \sum\limits_{j=0}^m \frac{a^{m-j} (-1)^j T^{\alpha j+1} }{(\frac{1}{\alpha} +m)_{(m-j+1)}} - \frac{1}{\alpha m!} \frac{a^{m + \frac{1}{\alpha}}}{(\frac{1}{\alpha})^{(m+1)}} \cdot \Gamma\left[1- \frac{1}{\alpha}, \frac{a}{T^\alpha} \right] \tag{1a} \\ &&=\frac{(-1)^m a^{\frac{1}{\alpha }+m} \Gamma \left(-\frac{\alpha m+1}{\alpha },a T^{-\alpha }\right)}{\alpha m!} % &&{\hat lhs}_{a,m}(T) := \\ &&\frac{1}{2} \sqrt{\pi } \left( \right. \\ && \left. % -\frac{2 e^{-\frac{a}{T^2}} }{\sqrt{\pi }} \sum\limits_{n_1=1}^{m} \sum\limits_{m_2=0}^{\left\lfloor \frac{1}{2} \left(2 n_1-1\right)\right\rfloor } \left( \frac{(-1)^{m_2} 2^{2 m-2 m_2-1} a^{m-m_2} T^{2 m_2+1}}{n_1 m_2! \left(2 m-2 n_1+1\right)! \left(-2 m_2+2 n_1-1\right)!} \right)+ \right. \\ && \left. % \frac{2 e^{-\frac{a}{T^2}} }{\sqrt{\pi }} \sum\limits_{n_1=0}^m \sum\limits_{m_2=0}^{\left\lfloor n_1\right\rfloor} \left( \frac{(-1)^{m_2} 2^{2 \left(m-m_2\right)} a^{m-m_2} T^{2 m_2+1}}{\left(2 n_1+1\right) m_2! \left(2 \left(m-n_1\right)\right)! \left(2 n_1-2 m_2\right)!} \right)+ % \right. \\ && \left.\frac{2^{2 m+1} a^{\frac{1}{2} (2 m+1)} \left(\text{erf}\left(\frac{\sqrt{a}}{T}\right)-1\right)}{(2 m+1)!} \right. \\ && \left. % \right) \tag{1b} \end{eqnarray}

where $m=1,2,3,\cdots$.

By using the Mathematica code snippet below we have checked that those coefficients match up to $m=4$. Here we go:

  a =.; b =.; alpha =.; T =.; M = 4;
rhs := Integrate[
   1/alpha (a - s)^(1/alpha) Gamma[-1/
      alpha, (a - s) T^(-alpha)] (b BesselI[1, 2 Sqrt[b s]])/
     Sqrt[b s], {s, 0, a}] + 
  Integrate[
   1/alpha (b - s)^(-1/alpha) Gamma[1/alpha, 
     0, (b - s) T^(alpha)] (a BesselI[1, 2 Sqrt[a s]])/Sqrt[a s], {s, 
    0, b}] + -T BesselI[0, 2 Sqrt[b (a)]] + 
  1/alpha a^(1/alpha) Gamma[-1/alpha, a T^(-alpha)] + 
  1/alpha b^(-1/alpha) Gamma[1/alpha, 0, b T^(alpha)]
lhs := Sqrt[Pi]/(2 Sqrt[b]) Exp[-2 Sqrt[a b]] - 
  Sqrt[Pi]/(4 Sqrt[b]) (Exp[2 Sqrt[a b]] Erfc[
       Sqrt[b] T + Sqrt[a] 1/T] + 
     Exp[-2 Sqrt[a b]] Erfc[Sqrt[b] T - Sqrt[a] 1/T])
(Coefficients of both sides at b^m for m=1,2,3...)
lhs10 = Drop[
   Simplify[
    CoefficientList[
     Normal[Series[lhs, {b, 0, M}, 
       Assumptions :> b > 0 && a > 0 && T > 0]], b], 
    Assumptions -> a > 0 && T > 0], 1];
lhs1 = 
  Table[
   Sqrt[Pi]/2 (

     (-1 + Erf[Sqrt[a]/T]) (2 Sqrt[a])^(2 m + 1)/(2 m + 1)! +
      E^(-a/T^2)
        2/Sqrt[Pi] Sum[((-1)^m2 2^(2 (m - m2)) a^(m - m2) T^(
          1 + 2 m2))/ 
         m2!  (Sum[
            1/((1 + 2 n1) (2 (m - n1))! (-2 m2 + 2 n1)!), {n1, m2, 
             m}] + Sum[-1/(
            2 n1  (1 + 2 m - 2 n1)! (-1 - 2 m2 + 2 n1)!), {n1, m2 + 1,
              m}]), {m2, 0, m}]
     )
   , {m, 1, M}];
a =.; alpha =.; T =.;
rhs1 = Table[ 
   1/m! 1/alpha Gamma[-1/
       alpha] a^(1/alpha + m) (1/alpha)! 1/(1/alpha + m)! + (
    a^m T  HypergeometricPFQ[{1, -(1/alpha)}, {1 - 1/alpha, 
       1 + m}, -a T^-alpha])/(m! m!) + 
    1/m! Sum[(-1)^
        m1 T^(m1 alpha + 1)/(m1 alpha + 1) a^(m - m1)/(m - m1)! , {m1,
        0, m}] + (-T) a^m/(m! m!)
   , {m, 1, M}];
rhs1a = Table[
   1/(alpha m!)
      Exp[-a/T^alpha] Sum[(a^(m - j) (-1)^j T^(alpha j + 1))/
      Pochhammer[1/alpha + m - (m - j + 1) + 1, m - j + 1], {j, 0, 
       m}] + 1/( alpha m!) a^(m + 1/alpha)/
     Pochhammer[1/alpha, 
      m + 1] (-Gamma[(-1 + alpha)/alpha, a T^(-alpha)])
   , {m, 1, M}];
(lhs1 - lhs10) // FullSimplify
(lhs1 - (rhs1 /. alpha :> 2)) // FullSimplify
FullSimplify[(rhs1 - rhs1a)] // PowerExpand

Answer 1

This is an answer to the first question.

Take $j \in {\mathbb N}$ and $0 \le j \le m$. By comparing expressions $(1a)$ and $(1b)$ all what we need to do is to prove a following identity:

\begin{eqnarray} &&rhs_j:=\\ &&\left( \sum\limits_{n_1=j+0}^m \frac{1}{(1+2 n_1)(2m-2 n_1)! (-2 j+2 n_1)!} + \sum\limits_{n_1=j+1}^m \frac{-1}{(0+2 n_1)(1+2m-2 n_1)! (-1-2 j+2 n_1)!} \right) \frac{2^{2(m-j)}}{j!} = \frac{1}{\alpha m!} \frac{1}{(1/\alpha + m)_{(m-j+1)}} =: lhs_j \end{eqnarray}

But now we have:

\begin{eqnarray} &&rhs_j=\\ &&\left(\frac{1}{(1+2 j)(2m-2 j)!} + \sum\limits_{n_1=2 j+2}^{2 m+1} \frac{(-1)^{n_1-1}}{n_1 (1+ 2 m-n_1)! (-1-2 j+n_1)!}\right) \frac{2^{2(m-j)}}{j!}= \\ && \left( \sum\limits_{n_1=0}^{2 m-2 j} \frac{(-1)^{n_1}}{(n_1+2 j+1)(2 m-2 j-n_1)! n_1 !} \right) \frac{2^{2(m-j)}}{j!}= \\ &&\frac{(2 j)!}{(2m+1)!} \frac{2^{2(m-j)}}{j!}= \\ &&\frac{2^j j! (2j-1)!!}{(2m+1)!} \frac{2^{2(m-j)}}{j!}=\\ && \frac{2^{2 m}}{(2m+1)!} \frac{(2j-1)!!}{2^j} =\\ lhs_j \end{eqnarray}

In[1145]:= m = RandomInteger[{1, 10}]; j =.; eta =.;
{alpha, a, T} = RandomReal[{1, 3}, 3]; alpha = 2;
l1 = Table[(Sum[
       1/((1 + 2 n1) (2 (m - n1))! (-2 j + 2 n1)!), {n1, j, m}] + 
      Sum[-1/(2 n1  (1 + 2 m - 2 n1)! (-1 - 2 j + 2 n1)!), {n1, j + 1,
         m}])  (2^(2 (m - j))) /j!, {j, 0, m}];

l2 = Table[(1/((1 + 2 j) (2 (m - j))! ) + 
      Sum[(-1)^(n1 - 1)/(
       n1  (1 + 2 m - n1)! (-1 - 2 j + n1)!), {n1, 2 j + 2, 
        2 m + 1}])  (2^(2 (m - j))) /j!, {j, 0, m}];
l3 = Table[(Sum[(-1)^(n1)/((n1 + 2 j + 1)  (2 m - 2 j - n1)! ( 
         n1)!), {n1, 0, 2 m - 2 j}])  (2^(2 (m - j))) /j!, {j, 0, m}];
l4 = Table[(2 j)!/(2 m + 1)!  (2^(2 (m - j))) /j!, {j, 0, m}];
l4a = Table[(2^j  j! (2 j - 1)!!)/(2 m + 1)!  (2^(2 (m - j))) /j!, {j,
     0, m}];
l4b = Table[ (2 j - 1)!!/(2 m + 1)!  (2^(2 (m) - j)) /1, {j, 0, m}];
l5 = 1/(alpha m!)
    Table[1/
    Pochhammer[1/alpha + m - (m - j + 1) + 1, m - j + 1], {j, 0, m}];
(l1 - l2)
(l1 - l3)
(l1 - l4)
(l1 - l4a)
(l1 - l4b)
(l1 - l5)


Out[1154]= {0, 0, 0, 0, 0, 0, 0, 0, 0}

Out[1155]= {0, 0, 0, 0, 0, 0, 0, 0, 0}

Out[1156]= {0, 0, 0, 0, 0, 0, 0, 0, 0}

Out[1157]= {0, 0, 0, 0, 0, 0, 0, 0, 0}

Out[1158]= {0, 0, 0, 0, 0, 0, 0, 0, 0}

Out[1159]= {0, 0, 0, 0, 0, 0, 0, 0, 0}

Answer 2

This is a partial answer to the second question. By using the representation -- second from the top in $(1a)$-- we found a following result for our integral:

\begin{eqnarray} &&rhs_{a,b}(T):= \int\limits_0^T \exp\left( -\frac{a}{u^\alpha} - b u^\alpha \right) du =\\ &&e^{-a T^{-\alpha }} \cdot \sum\limits_{j=0}^\infty \frac{(-1)^j b^j T^{\alpha j+1} \, _1F_2\left(1;j+1,j+\frac{1}{\alpha }+1;a b\right)}{j! (\alpha j+1)} + \\ &&-a^{\frac{1}{\alpha }} \, _0F_1\left(;1+\frac{1}{\alpha };a b\right) \Gamma \left(\frac{\alpha -1}{\alpha },a T^{-\alpha }\right) \tag{3} \end{eqnarray}

From numerical tests it seems to follow that the series have a infinite radius of convergence-- which is good. But now the question still remains what is the limit of this result when $T \rightarrow \infty$.

In[1710]:= {a, b, alpha, T} = 
 RandomReal[{0, 3}, 4, WorkingPrecision -> 50]; M = 50;
(*T=1/T;*)
l1 =
 Take[Accumulate[
    Flatten@(Table[
       Exp[-a/T^alpha] HypergeometricPFQ[{1}, {1 + j, 
          1 + 1/alpha + j}, a b]/( 
        j! (1 + j alpha)) ( (-1)^j b^j T^(alpha j + 1)), {j, 0, 
        M}])], -5] + 
  a^(1/alpha)
    Hypergeometric0F1[1 + 1/alpha, 
    a b] (-Gamma[(-1 + alpha)/alpha, a T^(-alpha)])
NIntegrate[Exp[-a/u^alpha - b u^alpha], {u, 0, T}]
Out[1711]= {0.02133704081131646175505079569936099802131653, 

0.02133704080942688858103605226282512528146042, 

0.02133704080985795756643855713611553196628419, 

0.02133704080976158349817294852610096499615789, 

0.02133704080978270756025453626977940452914675}
Out[1712]= 0.021337

Update:

From my answer to this question we know that:

\begin{eqnarray} &&rhs_{a,b}(\infty) = \\ && b^{-\frac{1}{\alpha}} \, _0F_1 \left(; 1- \frac{1}{\alpha}; a b \right) \Gamma\left( \frac{1}{\alpha} + 1\right) +\\ && -a^{\frac{1}{\alpha }} \, _0F_1\left(;1+\frac{1}{\alpha };a b\right) \Gamma \left(\frac{\alpha -1}{\alpha }\right) \end{eqnarray}

In[1979]:= {alpha} = RandomReal[{1, 2}, 1]; M = 20; {a, b} = 
 RandomReal[{1, 2}, 2];
NIntegrate[ Exp[-b t^alpha] Exp[-a/t^alpha], {t, 0, Infinity}]
(

Gamma[-1/alpha + 1] (a)^(+1/alpha) Hypergeometric0F1[1 + 1/alpha, 
 a b] +

Gamma[+1/alpha + 1] (b)^(-1/alpha) Hypergeometric0F1[1 - 1/alpha, 
    a b]
 )
Out[1980]= 0.0196757
Out[1981]= 0.0196757