Asymptotic bound for $\int_0^\infty \int_0^\infty (x+y)^m e^{-\frac{x^2}{2i} - \frac{y^2}{2j}} dx\, dy\;$ for $i$ and $j$ large

There is an exact solution for the integral (big surprise for me !).

I try to write the expression for $$-\frac{\sqrt{2} (m+1)}{j}\,I_m=T_1+T_2+T_3+T_4$$ $$T_1=-\frac{\sqrt{\frac{\pi }{2}} (m+1) i^{\frac{m+1}{2}} \Gamma \left(\frac{m+1}{2}\right) \left(\frac{2 j}{i}+2\right)^{m/2}}{\sqrt{j}}$$ $$T_2=2^{\frac{m+1}{2}} j^{m/2} \Gamma \left(\frac{m}{2}+1\right) \, _2F_1\left(\frac{1}{2},1;\frac{m+3}{2};-\frac{j}{i}\right)$$ $$T_3=\frac{2^{\frac{m+1}{2}} (m+1) (i+j) i^{m/2} \Gamma \left(\frac{m}{2}+1\right) \, _2F_1\left(1,\frac{1-m}{2};-\frac{1}{2};-\frac{j}{i}\right)}{j m}$$ $$T_4=-\frac{2^{\frac{m+1}{2}} (m+1) i^{m/2} \Gamma \left(\frac{m}{2}+1\right) (i-j (m-3)) \, _2F_1\left(1,\frac{1-m}{2};\frac{1}{2};-\frac{j}{i}\right)}{j m}$$

I have the feeling that I have been unable to simplify properly.

Edit

Using $$(x+y)^m=\sum_{k=0}^m \binom{m}{k}\, x^{m-k}\,y^k $$ $$\int_0^\infty\int_0^\infty x^{m-k}\, y^k\,e^{-\frac{x^2}{2 i}-\frac{y^2}{2 j}}\,dx\,dy=2^{\frac{m-2}{2}} i^{\frac{m+1-k}{2} }j^{\frac{k+1}{2}} \Gamma \left(\frac{k+1}{2}\right) \Gamma \left(\frac{m+1-k}{2} \right)$$ and then the hypergeometric functions by the summation.

If, as I suggested in a comment, we let $i=p^2$ and $j=a^2p^2$ $$i^{\frac{m+1-k}{2} }j^{\frac{k+1}{2}}=a^{k+1} p^{m+2}$$ which could be more comfortable.

A upper bound

(With the help of Maple)

With the substitution $u = x+y, v = y$, we have \begin{align} I &= \int_0^\infty \int_0^u u^m \mathrm{e}^{-(u-v)^2/(2i) - v^2/(2j)} \mathrm{d} v \mathrm{d}u\\ &= \int_0^\infty \sqrt{\frac{\pi ij}{2i+2j}}\, u^m \mathrm{e}^{-\frac{u^2}{2i+2j}} \left[\mathrm{erf}\Big(\tfrac{u}{i}\sqrt{\tfrac{ij}{2i+2j}}\Big) + \mathrm{erf}\Big(\tfrac{u}{j}\sqrt{\tfrac{ij}{2i+2j}}\Big) \right] \mathrm{d}u\\ &= \int_0^\infty \sqrt{\frac{\pi ij}{2i+2j}}\, u^m \mathrm{e}^{-\frac{u^2}{2i+2j}} \mathrm{erf}\Big(\tfrac{u}{i}\sqrt{\tfrac{ij}{2i+2j}}\Big)\mathrm{d}u \\ &\qquad + \int_0^\infty \sqrt{\frac{\pi ij}{2i+2j}}\, u^m \mathrm{e}^{-\frac{u^2}{2i+2j}} \mathrm{erf}\Big(\tfrac{u}{j}\sqrt{\tfrac{ij}{2i+2j}}\Big) \mathrm{d}u\\ &= \sqrt{\pi}(\tfrac{2i+2j}{ij})^{m/2}i^{m+1}\int_0^\infty w^m \mathrm{erf}(w)\mathrm{e}^{-w^2i/j} \mathrm{d} w\\ &\qquad + \sqrt{\pi}(\tfrac{2i+2j}{ij})^{m/2}j^{m+1}\int_0^\infty w^m \mathrm{erf}(w)\mathrm{e}^{-w^2j/i}\mathrm{d} w\\ &= \sqrt{\pi}(\tfrac{2i+2j}{ij})^{m/2}i^{m+1} \Big(\int_0^\infty w^m \mathrm{e}^{-w^2i/j} \mathrm{d} w - \int_0^\infty w^m (1 - \mathrm{erf}(w))\mathrm{e}^{-w^2i/j} \mathrm{d} w\Big)\\ &\qquad + \sqrt{\pi}(\tfrac{2i+2j}{ij})^{m/2}j^{m+1} \Big(\int_0^\infty w^m \mathrm{e}^{-w^2j/i}\mathrm{d} w - \int_0^\infty w^m (1-\mathrm{erf}(w))\mathrm{e}^{-w^2j/i}\mathrm{d} w\Big)\\ &= 2\sqrt{\pi}2^{m/2-1}(i+j)^{m/2}\sqrt{ij}\, \Gamma(\tfrac{m+1}{2})\\ &\qquad - \sqrt{\pi}(\tfrac{2i+2j}{ij})^{m/2}i^{m+1} \int_0^\infty w^m (1 - \mathrm{erf}(w))\mathrm{e}^{-w^2i/j} \mathrm{d} w\\ &\qquad - \sqrt{\pi}(\tfrac{2i+2j}{ij})^{m/2}j^{m+1} \int_0^\infty w^m (1-\mathrm{erf}(w))\mathrm{e}^{-w^2j/i}\mathrm{d} w\\ &\le 2\sqrt{\pi}2^{m/2-1}(i+j)^{m/2}\sqrt{ij}\, \Gamma(\tfrac{m+1}{2})\\ &\qquad - \sqrt{\pi}(\tfrac{2i+2j}{ij})^{m/2}i^{m+1} \int_0^\infty w^m \Big(\sqrt{\frac{2\mathrm{e}}{\pi}}\frac{\sqrt{\beta-1}}{\beta}\mathrm{e}^{-\beta w^2}\Big)\mathrm{e}^{-w^2i/j} \mathrm{d} w\\ &\qquad - \sqrt{\pi}(\tfrac{2i+2j}{ij})^{m/2}j^{m+1} \int_0^\infty w^m \Big(\sqrt{\frac{2\mathrm{e}}{\pi}}\frac{\sqrt{\beta-1}}{\beta}\mathrm{e}^{-\beta w^2}\Big)\mathrm{e}^{-w^2j/i} \mathrm{d} w\\ &= 2\sqrt{\pi}2^{m/2-1}(i+j)^{m/2}\sqrt{ij}\, \Gamma(\tfrac{m+1}{2})\\ &\qquad - \sqrt{\pi}(\tfrac{2i+2j}{ij})^{m/2}i^{m+1} \sqrt{\frac{\mathrm{e}}{2\pi}}\frac{\sqrt{\beta-1}}{\beta} (\beta +\tfrac{i}{j})^{-(m+1)/2}\Gamma(\frac{m+1}{2})\\ &\qquad - \sqrt{\pi}(\tfrac{2i+2j}{ij})^{m/2}j^{m+1} \sqrt{\frac{\mathrm{e}}{2\pi}}\frac{\sqrt{\beta-1}}{\beta} (\beta +\tfrac{j}{i})^{-(m+1)/2}\Gamma(\frac{m+1}{2}) \end{align} where $\mathrm{erf}(w) = \frac{2}{\sqrt{\pi}}\int_0^w \mathrm{e}^{-t^2}\mathrm{d} t$ is the error function, and we have used $1 - \mathrm{erf}(w) \ge \sqrt{\frac{2\mathrm{e}}{\pi}}\frac{\sqrt{\beta-1}}{\beta}\mathrm{e}^{-\beta w^2}$ (for $w\ge 0$, $\beta > 1$; see https://en.wikipedia.org/wiki/Error_function). We may choose $\beta = \frac{5}{4}$.

Wondering if I made or not a mistake somewhere in my previous answer, I restarted using the binomial expansion of $(x+y)^m$ and ended with something apparently simpler (but also apparently different). The end result write $$I=\frac{(2i)^{\frac m2}}{4j} \Big[\cdots\Big] $$ with $$\Big[\cdots\Big]=\Gamma \left(\frac{m-2}{2}\right) \left(i (i-j(m-3))-(i+j)^2 \, _2F_1\left(1,\frac{1-m}{2};-\frac{1}{2};-\frac{j}{i}\right)\right)+$$ $$2j \sqrt{\pi ij}\, \left(\frac{i+j}{i}\right)^{\frac m2} \Gamma \left(\frac{m+1}{2}\right)$$

If the first term can be neglected (or similar to the second one), then $$I\sim \sqrt \pi \,2^{\frac{m-2}{2}} \Gamma \left(\frac{m+1}{2}\right) (i+j)^{\frac m2}\sqrt{ij} $$ which looks like what you wrote.

Edit

In order to check, I made $j=i$ which makes $$\frac I{2^{\frac{m-4}{2}} i^{\frac{m+2}{2}}}=$$ $$\sqrt{\pi }\, 2^{\frac{m}{2}+1} \Gamma \left(\frac{m+1}{2}\right)-\left(4 \, _2F_1\left(1,\frac{1-m}{2};-\frac{1}{2};-1\right)+m-4\right) \Gamma \left(\frac{m-2}{2}\right)$$

I rewrote it as $$I=\sqrt{\pi }\, 2^{m-1}\, i^{\frac{m+2}{2}} \Gamma \left(\frac{m+1}{2}\right)\, (1-K)$$ with $$K=\frac {\left(4 \, _2F_1\left(1,\frac{1-m}{2};-\frac{1}{2};-1\right)+m-4\right) \Gamma \left(\frac{m}{2}-1\right) } {\sqrt{\pi }\, 2^{\frac{m+2}{2}} \Gamma \left(\frac{m+1}{2}\right) }$$

As shown below, when $m$ increases, factor $K$ tends very fact to $-1$ $$\left( \begin{array}{cc} m & K \\ 3 & -0.883883 \\ 4 & -0.924413 \\ 5 & -0.950175 \\ 6 & -0.966854 \\ 7 & -0.977796 \\ 8 & -0.985044 \\ 9 & -0.989880 \\ 10 & -0.993128 \\ 15 & -0.998968 \\ 20 & -0.999839 \end{array} \right)$$

In other words, at least for $j=i$, for large $m$, an asymptotics is $$I \sim \sqrt{\pi }\, 2^m\, i^{\frac{m+2}{2}} \,\Gamma \left(\frac{m+1}{2}\right)$$

Using Stirling approximation $$\log(I) =m\log \left(\frac{2 i m}{e}\right)+\log \left(\sqrt{2} \pi i\right)-\frac{1}{12 m}+O\left(\frac{1}{m^3}\right)$$ which seems to be very close to what you wrote.