Is $dx\,dy$ really a multiplication of $dx$ and $dy$? On the answers of the question Is $\frac{\textrm{d}y}{\textrm{d}x}$ not a ratio? it was told that $\frac{dy}{dx}$ cannot be seen as a quotient, even though it looks like a fraction. My question is: does $dxdy$ in the double integral represent a multiplication of differentials? The problem than can be generalized for a multiple integral.
 A: 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.
A: $dx$ and $dy$ aren't real numbers; they are things called differential forms. Thus, you can't use the real number multiplication operation to multiply them.
However, $dx \, dy$ is a thing, and it is not terribly unreasonable to define "the multiplication of $dx$ and $dy$ to be $dx \, dy$. The trick is that you have to make the inferences in the opposite direction from what you're used to -- to work out the first properties, it's not because you are understanding $dx \, dy$ in terms of multiplication, it's because you are using your understanding of $dx \, dy$ to figure out what 'multiplication' means.
There are some subtleties in what $dx \, dy$ means that I'm not up to fully explaining at the moment: e.g. it needs to talk about the orientation of a region, so that the fact that $dx \, dy = -dy \, dx$ can be properly explained. (you don't notice this fact when you do ordinary iterated integrals, since you flip the orientation of your region whenever you swap the order of $x$ and $y$, which cancels out the sign change)
A: 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.
A: Let's see what happens to $dx dy$ under a change of variables $x = f(z,w)$, $ y = g(z,w)$.
$$ dx = \frac{df}{dz}dz + \frac{df}{dw}dw $$
$$ dy = \frac{dg}{dz}dz + \frac{dg}{dw}dw $$
So locally this change of variables is a linear transformation and the dilation factor is the determinant:
$$ \left| \begin{array}{cc} 
\frac{df}{dz} & \frac{df}{dw} \\
\frac{dg}{dz} & \frac{dg}{dw} 
\end{array} \right| =\frac{df}{dz}\frac{dg}{dw} -\frac{df}{dw} 
\frac{dg}{dz}   $$
We can derive this using exterior algebra and the wedge product.
$$ dx \wedge dy = 
\left(\frac{df}{dz}dz + \frac{df}{dw}dw \right) \wedge
\left(\frac{dg}{dz}dz + \frac{dg}{dw}dw \right)
$$
Using the identities $dz \wedge dz = dw \wedge dw = 0$ and $dz \wedge dw = - dw \wedge dz$. We recover the Jacobian formula.

In our case, $dx \wedge dy $ behaves like $dx dy$ since they are in a sense perpendicular, "$dx \perp dy$".

A: 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.
A: 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'$.
A: In Cartesian linear coordinates just as we sum up separate rectangle areas as $ \int y. dx $ , we sum up prismatic volumes as $ \int z.dx. dy  $.  Here we consider an infinitesimal(differential) area as a product of two infinitesimal(differential) lengths. 
If we are considering curvilinear coordinates, we need to consider their Jacobian to multiply Cartesian differential area by the Jacobian in order to define relation between the old small area and the  new distorted small area. The way product is looked upon as an area  in the small is exactly the same as in Cartesian before it leads to the integrated in the large situation.While changing from Cartesian to polar,
$ dx dy = r * dr* d\theta $ where $ r$ is Jacobian $ J$ in the general case 
$J = r= \partial(x,y)/\partial ( r,\theta). $
A: No, because the differentials that appear in an integral are just notation; they only signify what variables are to be integrated over.
Consider a sum:
$$\sum_{i=1}^n i^2$$
The $i$ in the bottom of the summation symbol tells you that $i$ is the dummy variable here.
You could, in principle (though no one does this), notate integrals the same way:
$$\int_{x=a}^b x^2$$
And you could do it for a double integral:
$$\int_{y=c}^d \int_{x=a}^b e^{-x^2 - y^2}$$
The differentials being used in integrals are just empty notation.

Edit: Fly by Night suggests that differentials mean something because they transform.  There is a transformation involved when changing coordinates, but it is absolutely incorrect to say this has to do with the differential.
This can be seen by considering the integral abstractly as a map, with its internals being the limit of a Riemann sum.  Let $\mathcal I(f, \ell, d)$ be the integral of a function $f$ over a region described by the function $\ell$, with an interval $d$ describing the parameter region of $\ell$.  That is, if $d = [a,b)$, then we integrate $f$ from $\ell(a)$ to $\ell(b)$.
The integral can then be written as a limit of a Riemann sum:
$$\mathcal I(f, \ell, d) = \lim_{N \to \infty} \sum_{n=0}^{N-1} [f \circ \ell]\left (a + n \frac{b-a}{N} \right) \ell'\left(a + n \frac{b-a}{N}\right) \frac{b-a}{N} $$
Consider the case $\ell(x) = x$.  Then $\ell'(x) = 1$, and you get the usual, familiar form of a Riemannian integral.
Now instead, consider a reparameterization.  Let $x = g(u) = u^2$, as Fly by Night suggested.  Then there is a new paremeterization function $m$:
$$m(u) = (\ell \circ g)(u) = u^2$$
There is also a transformed interval:
$$e = [g^{-1}(a), g^{-1}(b)) = [p, q)$$
Now then, since $m([p, q)) = \ell([a, b))$, the integral of $f$ should be the same for both.  That is, $\mathcal I(f, \ldots)$ only really cares about the region $f$ is actually being integrated over, not how that region is parameterized.  So we should conclude that
$$\mathcal I(f, \ell, d) = \mathcal I(f, m, e)$$
But we know that $m'(u) = {\color{red}{2u}}$, and we get
$$\mathcal I(f, m, e) = \lim_{N \to \infty} \sum_{n=0}^{N-1} [f \circ m]\left( p + n \frac{q - p}{N} \right) \color{red}{2\left (p + n \frac{q-p}{N} \right)} \frac{q-p}{N}$$
What does this look like in traditional integral notation?  We started with
$$\int_a^b [f \circ \ell](t) \ell'(t) \, dt = \int_a^b f(t) (1) \, dt$$
And then we changed to
$$\int_{p = \sqrt{a}}^{q = \sqrt{b}} [f \circ m](t) m'(t) \, dt = \int_p^q [f \circ m](t) 2t \, dt$$
All of this follows from the transformation of the tangent vector along the parameterized curve.  Writing it as a change of differentials is nothing more than misleading voodoo.  It is a fundamental mistake to say that the differentials are 1-forms.  They are not; they are totally vacuous and without meaning, and thus my answer to the question is that it's meaningless to talk about multiplying differentials because they are themselves meaningless.  All of the supposed properties of differentials actually come from other, better defined geometric principles and concepts.
A: I find it very helpful in this case, and with multi-dimensional integrals in general, to write the $dx$ and $dy$ after the $\int$ sign. I.e.
$$ \int dy \int dx\, f(x,y) \;\;\;\text{ rather than }\;\;\; \int\int f(x,y)\, dx\,dy $$
This makes it clear that the $dx$ and $dy$, far from being infinitessimally small $\Delta x$ and $\Delta y$s, are part of the integral operators. Mathematical physicists often follow this convention, not just because it makes multi-dimensional integrals easier to understand, but also because it emphasizes  integration as an operator -- and a linear operator at that.
$\int dy$ and $\int dx$ each operate on $f(x,y)$ as "atomic units" that can't be split up. $\int dy$ and $\int dx$ apply to $y$ and $x$ respectively just as $\sum_y$ and $\sum_x$ are units that apply to $y$ and $x$. It still makes perfect sense to see that when $f(x,y)$ factorizes to  $f(x)f(y)$ you can rearrange
$$\int dy \int dx\, f(x)f(y) \;\;\;\text{ to  }\;\;\;  \int dy \,f(y) \int dx\, f(x)$$
just as you can
$$\sum_y \sum_x\, f(x)f(y) \;\;\;\text{ to  }\;\;\;  \sum_y \,f(y) \sum_x f(x)$$
However if you write $dxdy$ following the integrand it appears as though you are also factorizing $dxdy$ into $dx\times dy$ when in fact the $dx$ and $dy$ are part of the $\int$ operators.
One reason why $dx$ is often treated as if it is a standalone object to be manipulated separately is the convenience of doing just that when substituting variables. For example $\int \sin^2 x \cos x \, dx = \int u^2 \,du$, substituting $u=\sin x$ and $dx = du / \cos x$. This is very convenient and everyone does it but strictly speaking there is a rule being applied which transforms the integration operation from one to another: $\int dx\, \frac{du}{dx} = \int du$. Similarly in differentiation it is very convenient to think about cancelling the numerator of one derivative with the denominator of another, e.g. when applying the chain rule, but this is not strictly speaking correct.
