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)$.
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 $&gt;, &lt;$, 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{
  &amp; 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 
  &amp;  \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 
  &amp;  = \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 
  &amp;  = \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{
  &amp;  - \,\left( {{{s^{\,2}  + t^{\,2}  + 2st} \over {2\,i}}\,
 + {{s^{\,2}  + t^{\,2}  - 2st} \over {2\,j}}} \right) =   \cr 
  &amp;  =  - \,{{\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 
  &amp;  =  - \,{{\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 
  &amp;  =  - \,{{\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 &amp; 0  \cr 
   { - t/s^{\,2} } &amp; {1/s}  \cr 
 } } \right)} \right| = {1 \over s}
$$
and then proceed with approximation or series expansion of the Error Function.
%94 // FullSimplify

(1/(m (1 + m)))2^(
 m/2) (Sqrt[i j] (i + j)^(m/2) m Sqrt[\[Pi]] Gamma[(3 + m)/2] - 
   Gamma[1 + m/
      2] (j^(1 + m/2)
        m Hypergeometric2F1[1/2, 1, (3 + m)/2, -(j/i)] + (
      i^(m/2) (1 + m) (-i^2 + 
         i j (-3 + m) + (i + j)^2 Hypergeometric2F1[1, (1 - m)/
           2, -(1/2), -(j/i)]))/(j (-2 + m))))

This shows a professional grade solution. It is meant to be exact and not an approximation of some kind. It sets the requirements ot the question properly first and then integrate for all values that appear in the integrant and integral correctly. To type this in MathML is a really tedious task. So this is Mathematica, Wolfram Language code and can for some extent be verified on Wolfram Alpha. Have a try on the free Maple app too.
The names of the functions appear in the reference question partly and are not reduced to somehow more elementary functions in this calculation. This respects Mathematica great internal convention in favour to some universitary or college simplicity.
This is still a function of {i,j,m} and very lengthy. Only two terms have this $\sqrt{i j}$. Others are dependent on m in the exponent.
This can be calculated to numbers. For example for $j=10=i$ and $m=1,2,3$:
{
{79.2665},
{514.159},
{3963.33},
{34849.6},
{340846.},
{3.6449610^6},
{4.2090510^7},
{5.1989410^8},
{6.8177210^9},
{9.43489*10^10}
}
That looks pretty divergent for bigger ${i,j}$.
For example for $i=10$ and $j$ free and $m=1,2,3,...0$ this looks this way:
$
 5 \sqrt{j} \sqrt{2 \pi} + j \sqrt{5 \pi},$
$20 j + \sqrt{5/2} \sqrt{j} (10 + j) \pi,$
$ \sqrt{j} (100 \sqrt{2} + 30 \sqrt{5} \sqrt{j} + 15 \sqrt{2} j + 
 2 \sqrt{5} j^(3/2)) \sqrt{\pi},$
$ (80 j (10 + j) + 3 \sqrt{5/2} \sqrt{j} (10 + j)^2 \pi,$
$ \sqrt{j} (4000 \sqrt{2} + 1500 \sqrt{5} \sqrt{j} + 1000 \sqrt{2} j + 
 200 \sqrt{5} j^(3/2) + 75 Sqrt{2} j^2 + 
 8 \sqrt{5} j^(5/2)) \sqrt{\pi}, ...$
These are polynomials in j and $\sqrt{j}$.
If there would be another condition for the solution from the Mathematica engine for the input it would show that in the output or on the message cue. But there is none. OK. I am lazy because this is a tedious expression. I added an answer to the reference question too to earn more seriosity for my answer. This is not an expansion of any form and high mathematical intuition to take and offer that solution path for this highly complex problem and the restrictions on it.
$I$ is a function independent on $x$ and $y$. $I(i,j,m)$!
S rises strongly. If the summation is trunkated for fixed $m$ than the function may be plotted. The higher the terms that are taken into the sum the higher the function value. The plot is taken over $l$ and $k$.
I is a type of convolution integral. It acts in information theory and related as s low pass filter.
The binomial polynom is dominating for lower summation parts. The summands to not converge fast enough. This is what concerns one factor: Convolution. This requires that the other function is finite on the reals. That is not the case for $(x+y)^m$ unless it is limited to an finite interval. That is not the case in the question. Therefore the infinite summation fails. $S$ is returned unevaluated by Mathematica and other CAS.
The cognitive dissonance is the the summation is over discretized $x$ and $y$ with out and exhaustion. Integers are not dense in the reals. Approximation of convolution integrals require the finiteness of the convoluted function. This binomial is unlimited.
I selected $m=1$.
p1 = Plot3D[
  With[{m = 1}, 
   1/(m (1 + m)) 2^(
    m/2) (Sqrt[i j] (i + j)^(m/2) m Sqrt[\[Pi]] Gamma[(3 + m)/2] - 
      Gamma[1 + m/
         2] (j^(1 + m/2)
           m Hypergeometric2F1[1/2, 1, (3 + m)/2, -(j/i)] + (
         i^(m/2) (1 + m) (-i^2 + 
            i j (-3 + m) + (i + j)^2 Hypergeometric2F1[1, (1 - m)/
              2, -(1/2), -(j/i)]))/(j (-2 + m))))], {i, 1, 20}, {j, 1,
    20}]

p2 = Plot3D[\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(x = 0\), \(100\)]\(
\*UnderoverscriptBox[\(\[Sum]\), \(y = 0\), \(100\)]
\*SuperscriptBox[\((x + y)\), \(1\)] Exp[\(-1\)\ \((
\*SuperscriptBox[\((
\*FractionBox[\(x\), \(k\)])\), \(2\)] + 
\*SuperscriptBox[\((
\*FractionBox[\(y\), \(l\)])\), \(2\)])\)]\)\), {k, 1, 20}, {l, 1, 
   20}]



Apparently not very much alike. I seems to have something of a Puiseux_series.
1/(m (1 + m)) 2^(
   m/2) (Sqrt[i j] (i + j)^(m/2) m Sqrt[\[Pi]] Gamma[(3 + m)/2] - 
     Gamma[1 + m/
        2] (j^(1 + m/2)
          m Hypergeometric2F1[1/2, 1, (3 + m)/2, -(j/i)] + (
        i^(m/2) (1 + m) (-i^2 + 
           i j (-3 + m) + (i + j)^2 Hypergeometric2F1[1, (1 - m)/
             2, -(1/2), -(j/i)]))/(j (-2 + m)))) /. 
  m -&gt; 1 // FullSimplify

$(\sqrt{i}*j + \sqrt{i*j}*\sqrt{i + j} - 
   i*\sqrt{j}*(-1 + \sqrt{\frac{(i + j)}{i}}))*\sqrt{\frac{\pi}{2}}$
On the diagonal $i=j$ or $k=l$ for $m=1$
Plot[{\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(x = 0\), \(100\)]\(
\*SuperscriptBox[\((2  x)\), \(1\)] Exp[\(-1\)\ 2 
\*SuperscriptBox[\((
\*FractionBox[\(x\), \(k\)])\), \(2\)]]\)\), (k^(
     3/2) - (-1 + Sqrt[2]) k^(3/2) + 
     Sqrt[2] Sqrt[k] Sqrt[k^2]) Sqrt[\[Pi]/2]}, {k, 1, 1000}, 
 PlotRange -&gt; Full, AxesLabel -&gt; Automatic]


The sum S exits for all values, since exponential is stronger than polynomial and polynomial is kept fix while the argument of the exponential drops to zero with inverse square order. Blue is S and I is goldish. The infinite summation takes endless on my machine. The FullSimplify might get complex and needs extra care.
It is in general very complicated or delicate to approximate convolution integrals. The Gaussian function is problematic in representation and integration. Such profile function play an important roles in all measurements of modern technologies. So they are studied with great care. Since we know this function we just have to handle it with care. Mathematica has extra feature for dealing with the representation of very large and very small numbers. This both are needed unless the solution is symbolic.
There some complexity function around on this community. Mixing polynomials with exponentials is always very hard and complicated. Help might get easier if the wordings would be presented.
There might be impulses for a different approach by interpretation as a Kampe de Feriet Function. The given sum is related to the answer of this question: evaluation-of-sum-x-0-infty-e-x2. Mathematica can do this:
Sum[E^(-r^2), {r, 0, Infinity}]
1/2 (1 + EllipticTheta[3, 0, 1/E])

If we do this:
Sum[D[D[E^(-r^2), r], r], {r, 0, Infinity}]
Sum[-2/E^r^2 + (4*r^2)/E^r^2, {r, 0, Infinity}]

Nice for it is this Jacobi Theta Functions. So S is the second derivative of the EllipticTheta with first argument 3 with respect to the second argument evaluated at 0 for the second argument and 1/e for the third up to $2^3$ as a factor following the formula for the corresponding derivative on EllipticThetaPrime4. Now we have a representation for $i=j=1$ and $m=1$. This is a very simple cases of $S$.
We can plug our $e{-(\frac{x^2}{i^2}+\frac{y^2}{j^2})}$ into this too. That is only representation. But for $m&gt;1$ there is still no help. The $i$ and $j$ make the exponent still elliptical. The $m$ add circularity.
So the next step is hard again. Find relationships between EllipticTheta and Hypergeometric2F1. There are some for derivatives but this needs identities. Then generalize. Candidates are Elliptic hypergeometric series. Perhaps there is need for a revision of the hypothesis.