I'm interested in the asympotic behaviour (assuming an exact solution is intractable) of $$I(m,d)=\int\_0^{\infty} \left[ Q(m,x)\right]^d dx$$


for fixed $d \in \mathbb{N}$ (in particular, for $d=3$) and $m\to +\infty$.


Here $Q(m,x) = \frac{\Gamma(m,x)}{\Gamma(m)} $ is the upper regularized gamma function. 


[Empirically](https://math.stackexchange.com/questions/3370609/integral-involving-gamma-incomplete-function?noredirect=1&lq=1), it seems that $I(m,3) = m - a \sqrt{m} +O(1)$ with $a \approx 
0.835$



Asympotics of $\int_0^{\infty} \left[ Q(m,x)\right]^d dx$
Short:


The answer for $d=3$ is $a=\frac{3}{2\sqrt{\pi}}=0.846283 \dots$.


Long:


The integral is splitted in two at $x = m$.
$$
\int\_{0}^{m} \ Q(m,x)^d \ dx + \int\_{m}^{\infty} \ Q(m,x)^d \ dx
$$
For the asymptotic calculation of the integrals two approximations of the reguralized gamma function, $Q(m,x)$, for large $m$ are needed:


1. for $0\le x \le m$ (see, e.g., <http://dlmf.nist.gov/8.12.E4> ):
$$
Q(m, x)\sim \frac{1}{2} \text{erfc} \left(-\sqrt{m}\ \sqrt{\frac{x}{m}-1- \ln \frac{x}{m}}\right),
$$
2. for $m \le x$ (see, e.g., <http://dlmf.nist.gov/8.12.E18> )
$$
Q(m, x)\sim \frac{x^{m-\frac{1}{2}} e^{-x}}{\Gamma(m)} \ e^{\frac{(x-m)^2}{2\ x}} \text{erfc} \left(\frac{x-m}{\sqrt{2\ x}}\right).
$$


ad 1.


First we derive the following asymptotic integral for $m \rightarrow \infty$ ($d \in \mathbb{N}$) using an integral representation of the error function (see, e.g., <http://dlmf.nist.gov/7.7.E1> ). Please, note that there is no minus sign in the argument.
$$
\int\_{0}^{1} dy \ \text{erfc}(\sqrt{m}\ \sqrt{y-1-\ln y})^{n}
$$


$$
 = \left( \frac{2}{\pi}\right)^{n} \int\_{0}^{1} dy \ e^{-n\ m (y-1-\ln\ y)} \ \prod\_{i=1}^{n} \left( \ \int\_{0}^{\infty} dt\_{i}\ (1 + t\_{i}^{2} )^{-1} \ e^{- m \ t\_{i}^{2} (y-1-\ln y) }\right)
$$


$$
= \left( \frac{2}{\pi} \right)^{n} \left( \prod\_{i=1}^{n}\int\_{0}^{\infty} dt\_{i} \right)\ \left( \prod\_{i=1}^{n} (1 + t\_{i}^{2})^{-1} \right) \int\_{0}^{1} dy \ e^{-m\ (n + \sum\_{i}t\_{i}^{2})\ (y-1-\ln y)}.
$$
where the summation in the exponent is understood from $i=1$ to $n$. We assume, that all integrals converge to justify the exchange of integrations.


We use Laplace' method for asymptotic expansion of the $y$ integral: The exponent as a function of $y$ is expanded around its minimum $y=1$. With $\eta= y-1$
$$
y-1-\ln y \sim \frac{1}{2}\eta^2.
$$
After changing the integral variable to $\eta$ and extending the integration limit neglecting an exponentially small error, we get for the inner integral
$$
\int\_{0}^{1} dy \ e^{-m\ (d+\sum\_{i}t\_{i}^{2})\ (y-1-\ln y)} 
$$


$$
\sim \int\_{0}^{\infty} d\eta \ e^{-m\ (d+\sum\_{i}t\_{i}^{2})\frac{\eta^2}{2}\ } 
$$


$$
= \sqrt{\frac{\pi}{2\ m}} \left(d+\sum\_{i=1}^{d}t^{2}\_{i}\right)^{-1/2}.
$$


Together
$$
\int\_{0}^{1} dy \ \text{erfc}(\sqrt{m}\ \sqrt{y-1-\ln y})^{n} \sim
\left(\frac{2}{\pi}\right)^{n-\frac{1}{2}} \ m^{-1/2}\ I\_{n},
$$


with
$$
I\_{n}:= \left(\prod\_{i=1}^{n}\ \int\_{0}^{\infty}dt\_{i}\right)\
\left(n + \sum\_{i=1}^{n} t\_{i}^{2}\right)^{-1/2}\prod\_{i=1}^{n}\ \left(1 + t^{2}\_{i} \right)^{-1}
$$


Mathematica 11 gives the following results for $I\_{n}$
$$
I\_{1} = 1, \ I\_{2} = \pi\left(1-\frac{1}{\sqrt{2}}\right), \ I\_{3} = 1.0356625\dots, 
$$


$$
I\_{4} = 1.273085\dots, \ I\_{5} = 1.6458\dots,
$$


(I did not invest much effort to calculate the integrals symbolically. This might be worth a question on MO ?)


Altogether and by using the identity for the error function,
$$
\text{erfc} (-z)=2-\text{erfc}(z),
$$


one gets for the first part of the integral (a change to variable $y=x/m$ is included)
$$
\int\_{0}^{m}dx \ Q(m,x)^d \sim 2^{-d} \int\_{0}^{m}dx \ \text{erfc} \left(-\sqrt{m}\ \sqrt{\frac{x}{m}-1-\ln \frac{x}{m}}\right)^{d} 
$$


$$
= m \int\_{0}^{1}dy\ \left[1-\frac{1}{2}\text{erfc} (\sqrt{m}\ \sqrt{y-1-\ln y})\right]^{d} 
$$


$$
= m - m \sum\_{i=1}^{d}(-2)^{-i} {d \choose i} \int\_{0}^{1}dy\ \text{erfc}(\sqrt{m}\ \sqrt{y-1-\ln y})^{i}
$$


$$
= m - m^{-1/2} \sqrt{\frac{\pi}{2}}\sum\_{i=1}^{d}(-\pi)^{-i} {d \choose i}\ I\_{i}
$$


ad 2.


For the second part transform integration variable to $y=x/m$, use the same integral representation of the error function (see, e.g., <http://dlmf.nist.gov/7.7.E1> ) and interchange integrations
$$
\int\_{m}^{\infty}dx \ Q(m,x)^{d} \sim \int\_{m}^{\infty}dx\ \left(\sqrt{\frac{\pi}{2}} \frac{x^{m-\frac{1}{2}} \ e^{-x}}{\Gamma(m)} e^{-\frac{(x-m)^2}{2x}}\ \text{erfc}\left(\frac{x-m}{\sqrt{2x}}\right)\right)^{d} 
$$


$$
= \frac{m^{d(m-\frac{1}{2})+1}}{\Gamma(m)^d} \left(\frac{\pi}{2}\right)^{d/2} \int^{\infty}\_{1} dy\ e^{-d m \left(y+\frac{(y-1)^2}{2y}\right)}\ y^{d \left(m-\frac{1}{2}\right)} \ \text{erfc} \left(\sqrt{m}\frac{y-1}{\sqrt{2y}}\right)^{d}
$$


$$
 = \frac{m^{d \left(m-\frac{1}{2}\right) + 1}}{\Gamma(m)^d} \left(\frac{2}{\pi}\right)^{d/2} \left( \prod\_{i=1}^{d}\int\_{0}^{\infty}dt\_{i} \ (1+t\_{i}^{2})^{-1} \right) \int\_{1}^{\infty} dy \ e^{-m \left(d+\sum\_{i}t\_{i}^{2}\right) \frac{(y-1)^{2}}{2y} }\ y^{d \left( m - \frac{1}{2} \right)}
$$


The $y$ integral is approximated by Laplace' method. The integrand is written as
$$
\exp \left[ - m \left( \left( d + \sum\_{i} t^{2}\_{i} \right) \frac{ (y-1)^{2}}{2y} + d \ln y \right) \right] y^{-d/2}.
$$


The term in the square brackets is extremal at the lower integration limit $y=1$. Expanding to second order in $\eta=y-1$ gives
$$
-m d - m \ \eta^2 \ \frac{\sum\_{i}t\_{i}^{2}+d}{2}.
$$


Integrating from $\eta=0$ to infinity, again neglecting errors by exponentially small contributions, gives for the asymptotic approximation of the second part of the integral, $I(m,d)$,
$$
\int\_{m}^{\infty}dx \ Q(m,x)^d \sim \frac{m^{d \left(m-\frac{1}{2}\right) + 1}}{\Gamma(m)^d} \left(\frac{2}{\pi}\right)^{(d-1)/2} e^{- d m} m^{-1/2} I\_{d}
$$


The term $ y^{-d/2}$ in the integrand is not dependent on $m$ and can be set to 1 in this approximation.


After approximating the $\Gamma$ function by Stirling this gives for the second part of the integral, $I(d,m)$,
$$
\int\_{m}^{\infty}dx \ Q(m,x)^d \sim \frac{\pi^{-d+\frac{1}{2}}}{\sqrt{2}} m^{-1/2} I\_{d} 
$$


Taken all together we have asymptotically for large $m$
$$
I(m,d) \sim m - m^{-1/2} \sqrt{\frac{\pi}{2}}\left[ \sum\_{i=1}^{d-1}\left( \pi^{-i} (-1)^{i+1} {d \choose i}\ I\_{i}\right) - (1 + (-1)^{d})\ \pi^{-d} I\_{d} \right]
$$


and finally
$$
I(m,3) \sim m - m^{1/2}\frac{3}{2\sqrt{\pi}}
$$


