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 Whilst trying to count certain types of bipartite graphs, I'm lead to try to bound the following quantity
$$
I:=\int_0^\infty \int_0^\infty (x+y)^m e^{-\frac{x^2}{2i} - \frac{y^2}{2j}} dx\,dy
$$
where $i,j$ and $m$ are integers, and I'm interested in the asymptotics for large $i$ and $j$ and potentially $m$ (although it would suffice to have a good upper bound when $i \approx j$ and $m=o(i)$).
One can derive an exact expression for the integral by multiplying out the terms and using known identities for the quantities $\int_0^\infty x^k e^{-\frac{x^2}{2i}} dx$, however the asymptotics of this sum is unclear to me.
It would seem more natural to use a type of `saddle-point' method here, approximate the logarithm of the function around its maximum at $(x_0,y_0) = \left(i \sqrt{\frac{m}{i+j}},j \sqrt{\frac{m}{i+j}} \right)$ using the first two terms of the Taylor series, and so evaluate the integral in this region as a standard Gaussian, and then show that the contribution from outside this region is negligible.
This would lead to the following bound, which I would guess is in fact the correct asymptotic order
$$
I \approx \exp\left(m\log\sqrt{(i+j)(m)}-\frac{m}{2}\right)\pi\sqrt{2ij}.
$$
However, I can't get the regions in which the approximation is correct and the region in which the integral is negligible to overlap.
I suspect that this integral will have been considered somewhere in the literature, or at the very least will be susceptible to standard techniques in a field I'm not familiar with.
 A: 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: 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}$.
A: 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.
A: Not an answer, but this is too long for a comment. Perhaps one can expand $(x+y)^m$ into a sum (which will probably work out nicely if $m\in\Bbb{N}$, otherwise perhaps not) and use an identity given by Mathematica:
$$\int_0^\infty x^n\exp\left(\frac{-ax^2}{2}\right)\exp(-bx)\mathrm{d}x$$
$$=2^{\frac{n-1}{2}} a^{-\left(\frac{n+2}{2}\right)}\left( -b\sqrt{2} \ \Gamma \left(\frac{n+2}{2}\right) \ _{1} F_{1}\left(\left[\frac{n+2}{2} ,\frac{3}{2}\right] ;\frac{b^{2}}{2a}\right) +\sqrt{a} \ \Gamma \left(\frac{n+1}{2}\right) \ _{1} F_{1}\left(\left[\frac{n+1}{2} ,\frac{1}{2}\right] ;\frac{b^{2}}{2a}\right)\right)$$
For $\operatorname{Re}(a)>0$ and $\operatorname{Re}(n)>-1$. Here ${}_1F_1$ is the Kummer confluent hypergeometric function of the first kind, defined by the power series
$${}_1 F_1\left([\alpha,\beta];z\right)=\sum_{k=0}^\infty \frac{\Gamma(\alpha+k)\Gamma(\beta)}{\Gamma(\beta+k)\Gamma(\alpha)}z^k$$
For $\alpha,\beta,z \in\Bbb{C}$.
A: $$I=$$

For $i\inℕ$ and $j\inℕ$ always true. But $y\inℝ$ and $x\inℝ$ without restrictions.
If the restriction are presented first to the integral then the solutions looks this way:
$$I(x/;ℝ,y/;ℝ,i/;ℕ,j/;ℕ,m/;ℕ)=$$

This is done with Mathematica and Wolfram Language. I used the built-ins $Assuming$ and $Integrate$. This should be working with CAS that have not the high-level mathmatical knwoledge implemented.
