How to construct a Riemann integral from a Riemann sum that includes $(\Delta x)^2$ Suppose, for argument's sake, that we have a function $f: \mathbb{R}^2 \to \mathbb{R}$ that we want to integrate over a circular surface. We can use polar coordinates and cut up the $r$ and $\theta$ intervals into $N$ and $M$ pieces respectively, and then use that:
$$
\begin{aligned}
\hat S &= \sum_i^{N - 1}\sum_j^{M-1} f(r_i, \theta_j)\Delta S \\
&= \sum_i^{N - 1}\sum_j^{M-1} \frac{1}{2}f(r_i, \theta_j)(\theta_{j+1} - \theta_j)(r_{i+1}^2 - r_i^2) \\
&= \sum_i^{N-1}\sum_j^{M-1} \frac{1}{2}f(r_i, \theta_j)\Delta\theta((r_i + \Delta r)^2 - r_i^2) \\
&= \sum_i^{N-1}\sum_j^{M-1} \frac{1}{2}f(r_i, \theta_j)\Delta\theta(2r_i\Delta r + (\Delta r)^2) \\
&= \sum_i^{N-1}\sum_j^{M-1} f(r_i, \theta_i)r_i\Delta\theta\Delta r + \frac{1}{2}f(r_i, \theta_j)\Delta\theta(\Delta r)^2 \\
\end{aligned}
$$
Breaking this up, we can take the limit w.r.t.  $N$ and $M$ and get:
$$
\begin{aligned}
S &= \lim_{N, M \to \infty} \Big[\sum_i^{N-1}\sum_j^{M-1} f(r_i, \theta_i)r_i\Delta\theta\Delta r\Big] + \lim_{N, M \to \infty}\Big[\sum_i^{N-1}\sum_j^{M-1}\frac{1}{2}f(r_i, \theta_j)\Delta\theta(\Delta r)^2\Big] \\
&= \int_0^R\int_0^{2\pi} f(r, \theta)rdrd\theta + \lim_{N, M \to \infty}\Big[\sum_i^{N-1}\sum_j^{M-1}\frac{1}{2}f(r_i, \theta_j)\Delta\theta(\Delta r)^2\Big]
\end{aligned}
$$
Now, by the standard definition of the surface integral, the result should be:
$$
S = \int_0^R\int_0^{2\pi} f(r, \theta)rdrd\theta
$$
however that leaves the question: Why is it that
$$
\lim_{N, M \to \infty}\Big[\sum_i^{N-1}\sum_j^{M-1}\frac{1}{2}f(r_i, \theta_j)\Delta\theta(\Delta r)^2\Big] = 0
$$
And in general, how am I supposed to approach Riemann sums with non-linear, arbitrary $\Delta x$ terms, e.g. $(\Delta x)^2$, $\sqrt{\Delta x}$, etc.?
 A: I think I have found a way to prove this, although it involves working with the definition of a Riemann integral directly. I would be extremely happy if someone showed a simpler way to do this (if it exists).
Since the above example is a double integral where $\theta$ basically doesn't matter, I'm going to prove this for the general case of a function $f(x)$ over an interval $(a, b)$. Here it goes.
We want to show that for all $\epsilon$ there exists a given partition $\Pi = \{x_1, x_2, ...\}$ such that for any refinement $\Pi_1$ of $\Pi$
$$
\Bigg\lvert \sum_i^{\Pi_1} f(x_i)(\Delta x_i)^2 \Bigg\rvert < \epsilon
$$
A refinement is simply the original partition with more cuts. Therefore, we only need to show that splitting one sub-interval into two sub-intervals gives a smaller sum than the original. First, consider that if $\lvert f(x) \rvert < M$ for all $x\in (a, b)$, then
$$
\Bigg\lvert \sum_i^{\Pi} f(x_i)(\Delta x_i)^2 \Bigg\rvert \leq \sum_i^{\Pi} M(\Delta x_i)^2
$$
Now let us introduce another splitting point to $\Pi$, say $x_j$, between two other existing points $x_p < x_q$. The above sum changes to
$$
\sum_i^{\Pi_1} M(\Delta x_i)^2 = \sum_i^{\Pi}M(\Delta x_i)^2 - M(x_q - x_p)^2 + M(x_q - x_j)^2 + M(x_j - x_p)^2
$$
Setting $\alpha = x_q - x_p$, $\beta = x_q - x_j$ and $\gamma = x_j - x_p$, we know that:
$$
\alpha =  \beta + \gamma \implies \alpha^2 > \beta^2 + \gamma^2
$$
Therefore, we can say that:
$$
\Bigg\lvert \sum_i^{\Pi_1}f(x_i)(\Delta x_i)^2 \Bigg\rvert \leq \sum_i^{\Pi_1} M(\Delta x_i)^2 < \sum_i^{\Pi} M(\Delta x_i)^2
$$
Now we show that any $\epsilon$ can be expressed as $\sum_i^{\Pi} M_1(\Delta x_i)^2$ given the right $\Pi$. Essentially, we have to solve the following two equations:
$$
\sum_i^\Pi \Delta x_i = b - a \\
\sum_i^\Pi (\Delta x_i)^2 = \frac{\epsilon}{M_1}
$$
If we take $\Delta x_i = \frac{b-a}{N}$ for some $N$, then
$$
\sum_i^\Pi (\Delta x_i)^2 = N\frac{(b - a)^2}{N^2} = \frac{(b - a)^2}{N} = \frac{\epsilon}{M_1}
$$
Rewriting this, we have
$$
\epsilon = \frac{M_1(b - a)^2}{N}
$$
We need to show that we can choose a natural number $N$ and a real number $M_1$ such that $\lvert f(x)\rvert \leq M \leq M_1$ and the above is possible. If $\frac{\epsilon}{(b - a)^2} \geq M$, then taking $M_1 = \frac{\epsilon}{(b - a^2)}$ and $N = 1$ is sufficient. If $\frac{\epsilon}{(b - a)^2} < M$, then we take $N = \Big\lceil \frac{M(b-a)^2}{\epsilon} \Big\rceil$ and $M_1 = \frac{N\epsilon}{(b - a)^2}$, which gives us
$$
\epsilon = \frac{\frac{N\epsilon}{(b - a)^2}(b-a)^2}{N} = \epsilon
$$
and
$$
M_1 = \frac{\lceil \frac{M(b-a)^2}{\epsilon}\rceil\epsilon}{(b-a)^2} \geq \frac{\frac{M(b-a)^2}{\epsilon}\epsilon}{(b-a)^2} = M \implies M_1 \geq M
$$

Finally, we can say that for any $\epsilon$, we can construct a partition $\Pi$ such that:
$$
\sum_i^{\Pi}f(x_i)(\Delta x_i)^2 \leq \sum_i^{\Pi}M_1(\Delta x_i)^2 = \epsilon
$$
Above we showed that whenever we refine a partition, the resulting sum is always less than the $M_1$-bound of the original partition, therefore for every refinement $\Pi_1$ of $\Pi$
$$
\Bigg\lvert \sum_i^{\Pi_1}f(x_i)(\Delta x_i)^2 \Bigg\rvert < \epsilon
$$
Thus we conclude that the above sum converges to $0$, or in other words if the mesh of the partition is $P$
$$
\lim_{P\to 0} \sum_i^{\Pi} f(x_i)(\Delta x_i)^2 = 0
$$
Q.E.D.

On the question about having $\sqrt(\Delta x)$ etc., I believe that the above method would yield the same results.
