Krdytkyu

For a given function $f(x,y)$, its double integral is defined as:

$$iint_R f(x,y);mathrm{d}x;mathrm{d}y=lim_{ntoinfty}; sum_{i=1}^{n }f(x_{i},y_{i})Delta A_{i}$$

Where $Delta A_{i}=Delta x_{i}Delta y_{i}$.

Later I was told it's $mathrm{d}x wedge mathrm{d}y$ not $mathrm{d}x;mathrm{d}y$. $wedge$ being the wedge product.

This implies for instance that in polar, the differential with respect one integrates is given by:

$mathrm{d}x wedge mathrm{d}y=rmathrm{d}r wedge mathrm{d} theta $

Whereas if it were the ordinary product $dxdy$ then we would have

begin{align}
dxdy & = left(dr cos theta - r sin theta dtheta right) left( dr sin theta + r cos theta dthetaright)\
& = dr^2 cos theta sin theta - r^2 dtheta^2 cos theta sin theta + r dr dtheta (cos^2 theta - sin^2theta )
end{align}

But nowhere in the above definition

$$iint_R f(x,y);mathrm{d}x;mathrm{d}y=lim_{ntoinfty}; sum_{i=1}^{n }f(x_{i},y_{i})Delta A_{i}$$

I see the wedge product. It Seems to me that $$lim_{ntoinfty}; sum_{i=1}^{n }Delta A_{i}=iint mathrm{d}A=iint dxdy$$

Then where does the wedge product arise in the definition?

Sorry for being incredibly late to the party, but I recently stumbled upon this question after falling into a similar rabbit hole with my research, and I had to piece together my understanding from several sources. After a lot of thinking, I feel as though I have a solid answer to this question, and even though it is late, hopefully it may aid those who may stumble upon this question in future as I did.

There are three main motivations for favouring the wedge product $text dxwedge text dy$ over $text dx text dy$. The first is that the latter is simply a special case of the former--only applicable in the case of orthogonal axes. The second reason is that the latter does not generalise orientation in a way that is consistent with $1$D integration. The third (which was covered by Thomas' answer) is that it makes the algebra for transformations self-consistent.

$textbf{Infinitesimal Area}$

To expand on the first motivation, consider the question: "just what is $text dx text dy$"? Well, intuitively it should be a generalisation of the $1$D infinitesimal line element $text dx$. Extending this, we expect $text dxtext dy$ to be an infinitesimal area element. However, there are many ways in which we can extend a $1$D line to a $2$D area. We focus primarily on infinitesimal parallelograms, as parallelograms are obtained via linear transformations on squares and any transformation can be approximated locally by a linear transformation (this is essentially where the Jacobian comes into play for changing of integration variables, as it is a first order approximation to a coordinate transformation).

In the special scenario where $hat x$ and $hat y$ are orthogonal directions, then the parallelogram $text dx text dy$ corresponds to a square and its area is given exactly by the product $text dx text dy$. In contrast, if we have two directions $hat u$ and $hat v$ which are not orthogonal (but still linearly independent), then the area of the parallelogram would not be given by their product, but rather by the parallelogram area formula. Since $text du$ and $text dv$ are difference vectors, they lie in the plane, and their cross product $text dutimes text dv$ gives the area of the parallelogram that we seek. It is quite clear from this analysis that the cross product (and in higher dimensions, the wedge product) should be the default operation, and the fact that the standard product works is only a lucky coincidence afforded to us by a convenient choice of coordinates.

$textbf{Orientation}$

We know that integration is oriented, that is $int_a^b f(x) text dx = - int_b^a f(x) text dx$. However, this is not captured in $2$D by the standard area product $text dxtext dy$. If we imagine tracing out the square $text dxtext dy$, then we expect these four equations to be the same orientation (and thus equivalent)
begin{align*}
int_c^dint_a^b f(x,y) text dxtext dy\
int_d^cint_b^a f(x,y) text dxtext dy\
int_a^bint_d^c f(x,y) text dytext dx\
int_b^aint_c^d f(x,y) text dytext dx
end{align*}
The first two correspond to tracing out the $hat x$ direction first, and the latter two correspond to tracing out $hat y$ direction first. It seems I don't have the reputation to post an image yet, but here is the link to an album which shows all the four orders of integration. The remaining four equations that we can form by permutation of boundaries and integration order will represent the reverse orientation, and hence be the negative of the above equations. Focusing on one relation in particular, we expect
begin{equation}
int_c^dint_a^b f(x,y) text dxtext dy = -int_c^dint_a^b f(x,y) text dytext dx
end{equation}
In the case for the simple product $text dxtext dy$, only half (by combinatoric exhaustion) of the above mentioned relations are consistent. To fix this, we simply enforce anticommutivity by replacing $text dxtext dy$ with $text dxwedgetext dy$, which automatically carries the rule $text dxwedgetext dy = -text dywedgetext dx$.

$textbf{Conclusion}$

The takeaway message from these two examples (as well as the transformation rule given in Thomas' answer) is that the naïve area element $text dxtext dy$ just isn't well equipped to handle any of the subtler points that arise when generating a good definition of multi-variable integration. $text dxtext dy$ is really just an idealised case of $text dx timestext dy$, and the latter is what should've been in your integration formula from the start, as it is the more general formula for the area $text{dA}$. Then, when generalising to more dimensions still, the cross product should be replaced with the wedge product.

Short answer: it's the transformation formula.

Slightly longer answer, with some handwaving:
The $n$-fold wedge product in $n$ dimensions is basically the same as the determinant (the determinant and the wedge product are both alternating n-Forms and that space is $1$ dimensional on $n$-dimensional vector spaces, so it's clear that the first is a multiple of the second).

Volume (area in two dimensions) transforms under linar transformations by multiplication with the determinant of the map associated with the transformation. If you look at a rectangle and apply a linear map to it the volume change is given by the determinant of the linear map. For differentiable transformation the determinant is replaced by the determinant of the differential of the transformation (up to sign, but that's a technical thing I'll ignore here). If you already came across the transformation formula you will have seen the Jacobian pop up:

$$ int_Omega fcirc phi(x) ,|detDphi|(x) ,dx= int_{phi(Omega)} f(y) , dy $$
here, $dx$ is actually an abbreviation for $dx^1 dx^2cdots dx^n$, indicating the variables to integrate over. So the 'infinitesimal' area change is given by the determinant of the linearization of the transformation.

This is just the transformation behaviour of the wedge product, which is no surprise if you are aware of the remark I started with. So the wedge product comes in naturally via the transformation behaviour of volume. In your question it would enter on transforming the $Delta A$ term (like in the proof of the transformation formula). In the polar coordinates example $phi$ from the transformation formula corresponds to the coordinate change.