Is $dx\,dy$ really a multiplication of $dx$ and $dy$?
Solution 1:
In a double integral, you are actually integrating a differential two-form:
$$\int_R \mathrm{f}(x,y) \ \mathrm{d}x \wedge \mathrm{d}y$$
Here, $\mathrm{d}x$ and $\mathrm{d}y$ are the basis differential one-forms and $\mathrm{d}x \wedge \mathrm{d}y$ is their exterior product.
Solution 2:
Just as one can think of the derivative in Robinson's framework as a true ratio $\frac{\Delta y}{\Delta x}$ modulo an infinitesimal error (eliminated by applying the shadow), so also one can think of a single-variable integral as an infinite sum of infinitesimal terms of type $dx$ (again up to applying shadow). Double integrals can naturally be viewed as double (infinite) sums, where $dx\,dy$ is most decidedly an ordinary product. And of course this generalizes to multiple integrals as the OP suggested. If one is in Euclidean space, talking about differential forms is an unnecessary obfuscation.
Edit 1: For finite Riemann sums approxiating the double integral, it is obvious that the term $\Delta x \Delta y$ is a product; it seems nobody in his right mind would deny this. The difference is that one cannot deduce the value of the integral from a finite Riemann sum. On the other hand, with an infinite Riemann sum when $\Delta x$ is replaced by $dx$, etc., the value of the integral is deduced from the value of the Riemann sum by taking the shadow (see above). That's the advantage of having the richer syntax of the hyperreal approach.
Edit 2: the OP's question is in fact equivalent to a question about single-variable integrals, namely: does $f(x)dx$ denote multiplication of $f(x)$ by $dx$? Perhaps the right answer is that it denotes a memory of multiplication. Namely, the multiplication is still there at the level of the hyperfinite Riemann sum. To pass from this to the integral one applies the standard part function, after which we have only a "memory" left. Similarly, one can form the differential quotient Δy/Δx which is still a ratio, but one doesn't get the derivative until one applies the standard part function. Here also there is only a memory of a division left. The advantage of the hyperreal framework is that one has a direct procedure for passing from the ratio to the derivative which isn't the case in the traditional real-based framework where one must appeal to an indirect notion of an epsilon, delta limit.
A survey of various approaches to Robinson's framework is due to appear in Real Analysis Exchange.
Solution 3:
As others have said $dx\, dy$ does not represent a product of differentials. But it represents a product of measures. We have the "natural" Lebesgue measure $\lambda$ on the $x$-axis, and integration with respect to this measure is signalled by writing ${\rm d}x$ as right parenthesis of the integral. Similarly we have the "natural" Lebesgue measure $\lambda$ on the $y$-axis, and integration with respect to this measure is signalled by writing ${\rm d}y$ as right parenthesis of an integral involving the variable $y$. The individual measures $\lambda$ on the $x$-axis ${\mathbb R}$ and the $y$-axis ${\mathbb R}$ define a product measure $\lambda\otimes\lambda$ on the cartesian product ${\mathbb R}^2$, again called Lebesgue measure on ${\mathbb R}^2$. Integration with respect to this product measure is is signalled by writing ${\rm d}(x,y)$, ${\rm d}x\otimes {\rm d}y$, or simply $dx\,dy$, as right parenthesis of an integral over some subset $A\subset{\mathbb R}^2$. Fubini's theorem then tells us that $$\int\nolimits_A f(x,y)\>{\rm d}(x,y)=\int\nolimits_{A'}\left(\int\nolimits_{A_x} f(x,y)\> {\rm d}y\right)\ {\rm d}x\ ,$$ where $A'$ denotes the projection of $A$ onto the $x$-axis and $A_x:=\{y\mid (x,y)\in A\}$ collects the $y$-values to be weighted in for given $x\in A'$.
Solution 4:
Let's say $f(x,y)$ is in $\dfrac{\mathbf{kg}\cdot \mathbf{m}}{\mathbf{sec}\cdot\mathbf{dollar}}$ and $x$ is in $\mathbf{sec}$ and $\mathbf{dollar}$. Then $f(x,y)\,dx\,dy$ is in $\mathbf{kg}\cdot \mathbf{m}$, just as if we are multiplying.
If an infinitely small rectangle has length and width respectively $dx$ and $dy$, then its area is $dx\,dy$; if $f(x,y)$ is density of something (mass, probability, energy$\ldots$) with respect to area, then $f(x,y)\,dx\,dy$ is measure in the same units as that "something". Why does $dx\,dy$ become $r\,dr\,d\theta$. Sometimes people say "because you multiply by a Jacobian". That's Ok as far as it goes, maybe. I regard it as $(dr)(r\,d\theta)$. If $r$ is in meters and $\theta$ is dimensionless and in radians, then $r\,d\theta$ is in meters. A length of an arc of a circle is the radius times the radian measure of the arc. The radius is $r$; the radian measure of the arc is $d\theta$, so $r\,d\theta$ is the length, and $dr$ is disance in a direction orthogonal to that so $(dr)(r\,d\theta)$ is area of that rectangle. It's a rectangle because an infinitely small arc is a straight line.
If you say that this is not rigorous, I agree.
If you object that this is not rigorous, I disagree.
Intuitive ideas can be made rigorous [comment inspired by comments below: The following should be obvious, but apparently there was one person to whom it wasn't, so maybe there are others. I do not condone making infinitesimals rigorous in most first-year calculus courses.]. What is the right way to do that may be subject to philosophical disagreements. But one should not assert that intuition to be made rigorous is the rigorous end-product. Just which way of making something rigorous is the right one depends on the context. Some other way of making something rigorous that will be discovered 100 years from now may have its place. But the idea that is to be made rigorous exists independently of the ways of making it rigorous.