Approximating a double sum by a double integral Related to this question, I'm interested in bounding from above the following sum
$$
S:=\sum_{x=0}^\infty \sum_{y=0}^\infty (x+y)^m e^{-\frac{x^2}{2i} - \frac{y^2}{2j}},
$$
which I hope to do by relating it to the integral
$$
I:=\int_0^\infty \int_0^\infty (x+y)^m e^{-\frac{x^2}{2i} - \frac{y^2}{2j}} dx\,dy.
$$
Answers to the previous questions confirmed my expectation that $I = O\left(\exp\left(m\log\sqrt{(i+j)(m)}-\frac{m}{2}\right)\sqrt{ij}\right)$, the intuition for which is probably that the function behaves approximately like a gaussian 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)$.
However I've been unable to show that the difference $|I-S|$ is significantly smaller than this bound. For simple one dimensional integrals, for example with a unique maximum, it's not too hard to bound this difference 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.
I suspect that there should be a relatively standard way to approximate $|I-S|$, and I also wouldn't be surprised if someone more proficient in computing can get a CAS to provide a proof. The former would be more useful, just so that I have a tool for approaching similar questions.
So, very explicitly, I would like to know if
$$
|I-S| = o\left(\exp\left(m\log\sqrt{(i+j)(m)}-\frac{m}{2}\right)\sqrt{ij}\right),
$$
where even big-O would be sufficient for the application I have in mind, and I wouldn't be surprised if the difference is even bounded by a multiple of the maximum of the function. I'm interested in the asymptotics for $i$ and $j$ tending to infinity, $m$ can be fixed or also a function of $i$ and $j$. For the application I have in mind it would probably be sufficient to have such a result for $i = (1+o(1))j$ and $m = o(i)$.
 A: I cannot provide an actual answer , but only some considerations and hint
which might hopefully be helpful.
The function
$$
f(x,y) = \left( {x + y} \right)^{\,m} e^{ - \,{{x^{\,2} } \over {2\,i}}\, - {{y^{\,2} } \over {2\,j}}} 
$$
having a (cutted) bell-shape in the first quadrant, means that it is concave around the maximum and convex
further from it.
This makes quite difficult to relate the integral to the  Riemann sum with a $>, <$, because the sign
of the inequality changes in the two areas.
Furthermore, at increasing $i, \, j$, while the position of the max move $\approx \sqrt{i}$, and so approximately does its spread
the bell peak increases $\approx i^{m/2}$.
Since the $\Delta x , \, \Delta y$ of the sum are fixed at $1$, I doubt that the sum might converge to the integral.
Concerning the integral I would try the following approach
$$
\eqalign{
  & I = \int_{y\, = \,0}^{\,\infty } {\int_{x\, = \,0}^{\,\infty }
 {\left( {x + y} \right)^{\,m} e^{ - \,{{x^{\,2} } \over {2\,i}}\, - {{y^{\,2} } \over {2\,j}}} dxdy} }  =   \cr 
  &  \Rightarrow \left\{ \matrix{
  s = x + y \hfill \cr 
  t = x - y \hfill \cr}  \right.\quad  \Leftrightarrow \quad \left\{ \matrix{
  x = \left( {s + t} \right)/2 \hfill \cr 
  y = \left( {s - t} \right)/2 \hfill \cr}  \right.\quad  \Rightarrow   \cr 
  &  = \int_{y\, = \,0}^{\,\infty } {\int_{x\, = \,0}^{\,\infty }
 {s^{\,m} e^{ - \,{{\left( {s + t} \right)^{\,2} } \over {2\,i}}\,
 - {{\left( {s - t} \right)^{\,2} } \over {2\,j}}} {1 \over 2}dsdt} }  =   \cr 
  &  = \int_{s\, = \,0}^{\,\infty } {\int_{t\, = \, - s}^{\,s}
 {s^{\,m} e^{ - \,{{\left( {s + t} \right)^{\,2} } \over {2\,i}}\,
 - {{\left( {s - t} \right)^{\,2} } \over {2\,j}}} {1 \over 2}dsdt} }  \cr} 
$$
then also consider that
$$
\eqalign{
  &  - \,\left( {{{s^{\,2}  + t^{\,2}  + 2st} \over {2\,i}}\,
 + {{s^{\,2}  + t^{\,2}  - 2st} \over {2\,j}}} \right) =   \cr 
  &  =  - \,{{\left( {i + j} \right)s^{\,2} } \over {2\,i\,j}}
\left( {\left( {{t \over s}} \right)^{\,2}  - 2{{i - j} \over {i + j}}
\left( {{t \over s}} \right) + 1} \right) =   \cr 
  &  =  - \,{{\left( {i + j} \right)s^{\,2} } \over {2i\,j}}\left( {\left( {{t \over s}} \right)^{\,2}
  - 2{{i - j} \over {i + j}}\left( {{t \over s}} \right) + \left( {{{i - j} \over {i + j}}} \right)^{\,2}  + 1 - \left( {{{i - j} \over {i + j}}} \right)^{\,2} } \right) =   \cr 
  &  =  - \,{{\left( {i + j} \right)s^{\,2} } \over {2\,i\,j}}
\left( {\left( {\left( {{t \over s}} \right) - {{i - j} \over {i + j}}} \right)^{\,2}
  + 1 - \left( {{{i - j} \over {i + j}}} \right)^{\,2} } \right) \cr} 
$$
we can change the variables again
$$
\left\{ \matrix{
  s = s \hfill \cr 
  r = t/s \hfill \cr}  \right.\quad J = \left| {\left( {\matrix{
   1 & 0  \cr 
   { - t/s^{\,2} } & {1/s}  \cr 
 } } \right)} \right| = {1 \over s}
$$
and then proceed with approximation or series expansion of the Error Function.
