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, it seems that $I(m,3) = m - a \sqrt{m} +O(1)$ with $a \approx 0.835$
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:
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}} $$
