Asymptotic bound for $\sum_{x=0}^\infty \sum_{y=0}^\infty (x+y)^m e^{-\frac{x^2}{2i} - \frac{y^2}{2j}}$ for $i$ and $j$ large Note: This question relates to two previous questions on math.stackexchange (1 and 2), neither of which had satisfactory answers after posting bounties.
Whilst trying to count certain types of bipartite graphs, I'm lead to try to bound the following quantity
$$
S:=\sum_{x=0}^\infty \sum_{y=0}^\infty (x+y)^m e^{-\frac{x^2}{2i} - \frac{y^2}{2j}}
$$
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)$).
A natural strategy seems to be to consider the integral
$$
I:=\int_0^\infty \int_0^\infty (x+y)^m e^{-\frac{x^2}{2i} - \frac{y^2}{2j}} dx\,dy
$$
and show that $S \approx I$, and hopefully bound $I$ by analytic means. Unfortunately both parts of this strategy have been causing me difficulty.
One can derive an exact expression for $I$ 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)$, where the function takes the value $\exp\left(m\log\sqrt{(i+j)(m)}-\frac{m}{2}\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 believe 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 wasn't able to get this approach to work, there is a gap between where I can show the approximation to the Gaussian holds and where the integral is negligible. With the help of a CAS this quantity can be shown to be an upper bound for $I$, which would be sufficient for my application, but it would be nice to know of a more natural method to approach such problems (and I suspect there are standard tools for this sort of thing).
With this, it would be sufficient to show that
$$
|I-S| = o\left(\exp\left(m\log\sqrt{(i+j)(m)}-\frac{m}{2}\right)\sqrt{ij}\right),
$$
and I wouldn't be surprised if the difference is even bounded by a multiple of the maximum of the function. Indeed, for functions of a single variable with a unique maximum it's possible to bound the difference between the integral and the sum in terms of the maximum by considering appropriate telescoping sums. However, a naive analogue of this argument doesn't seem to work in two dimensions, and trying to apply this argument to each `slice' of the integral led to some pretty horrendous calculations. I also looked into using the Euler-Maclaurin formula but it's a bit out of my area of expertise.
Of course, there may be a way to bound or approximate $S$ without considering $I$ at all, which I would be happy to hear about!
 A: Let $a:=\sqrt i$ and $b:=\sqrt j$. Let us find the asymptotics of $I$ (leaving the proof of $S\sim I$ for a hopefully near future). Passing, if necessary, to a subsequence, we see that without loss of generality (wlog) one of the following two cases occurs:
Case 1: $m\to\infty$ or Case 2: $m=m_0$ is constant.
Also, wlog
\begin{equation}
    t_{a,b}:=\arccos\frac a{\sqrt{a^2+b^2}}\to t_0
\end{equation}
for some $t_0\in[0,\pi/2]$.
By substitutions $x=ar\cos t$ and $y=br\sin t$,
\begin{equation}
    I=ab \int_0^\infty r\,dr\,e^{-r^2/2}r^m\ J_m=ab\, 2^{m/2} \Gamma \left(\frac{m}{2}+1\right) J_m,
\end{equation}
where
\begin{equation}
    J_m:=\int_0^{\pi/2} dt\,e^{mg(t)},\quad g(t):=\ln(a\cos t+b\sin t); 
\end{equation}
this is the key observation.
A straightforward saddlepoint approximation (noting, in particular, that $g''(t)=-\frac{a^2+b^2}{(a \cos t+b \sin t)^2}\le-1$) shows that
\begin{equation}
    J_m\sim\sqrt{2\pi/m}\,(a^2+b^2)^{m/2}  
\end{equation}
if $m\to\infty$, that is, if Case 1 occurs, whence
\begin{equation}
    I\sim ab\, 2^{m/2} \Gamma \left(\frac{m}{2}+1\right) \sqrt{2\pi/m}\,(a^2+b^2)^{m/2},
\end{equation}
with the asymptotics of $ \Gamma \left(\frac{m}{2}+1\right)$ obtained by Stirling's formula, if desired.
In Case 2,
\begin{equation}
    J_m\sim (a^2+b^2)^{m/2} \int_0^{\pi/2} dt\,(\cos t_0\,\cos t+\sin t_0\sin t)^m
    =(a^2+b^2)^{m/2} \int_0^{\pi/2} dt\,\cos^m(t-t_0), 
\end{equation}
whence
\begin{equation}
    I\sim ab\, 2^{m_0/2} \Gamma \left(\frac{m_0}{2}+1\right) (a^2+b^2)^{m_0/2} \int_0^{\pi/2} dt\,\cos^{m_0}(t-t_0). 
\end{equation}
