Where does the wedge product arise in the definition of an integral? 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_{n\to\infty}\; \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=r\mathrm{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\ d\theta \right) \left( dr \sin \theta + r \cos \theta\ d\theta\right)\\
& = dr^2 \cos \theta \sin \theta  - r^2 d\theta^2 \cos \theta\ \sin\ \theta + r\ dr\ d\theta\ (\cos^2 \theta\ - \sin^2\theta )
\end{align}
But nowhere in the above definition
$$\iint_R f(x,y)\;\mathrm{d}x\;\mathrm{d}y=\lim_{n\to\infty}\; \sum_{i=1}^{n }f(x_{i},y_{i})\Delta A_{i}$$
I see the wedge product. It Seems to me that $$\lim_{n\to\infty}\; \sum_{i=1}^{n }\Delta A_{i}=\iint \mathrm{d}A=\iint dxdy$$
Then where does the wedge product arise in the definition?
 A: 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 f\circ \phi(x) \,|detD\phi|(x) \,dx= \int_{\phi(\Omega)} f(y) \, dy $$
here, $dx$ is actually an abbreviation for $dx^1 dx^2\cdots 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.
A: 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 dx\wedge \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 dx\text 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 du\times \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 dx\text dy$. If we imagine tracing out the square $\text dx\text dy$, then we expect these four equations to be the same orientation (and thus equivalent)
\begin{align*}
    \int_c^d\int_a^b f(x,y) \text dx\text dy\\
    \int_d^c\int_b^a f(x,y) \text dx\text dy\\
    \int_a^b\int_d^c f(x,y) \text dy\text dx\\
    \int_b^a\int_c^d f(x,y) \text dy\text 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^d\int_a^b f(x,y) \text dx\text dy = -\int_a^b\int_c^d f(x,y) \text dy\text dx
\end{equation}
In the case for the simple product $\text dx\text dy$, only half (by combinatoric exhaustion) of the above mentioned relations are consistent. To fix this, we simply enforce anticommutivity by replacing $\text dx\text dy$ with $\text dx\wedge\text dy$, which automatically carries the rule $\text dx\wedge\text dy = -\text dy\wedge\text 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 dx\text 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 dx\text dy$ is really just an idealised case of $\text dx \times\text 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.
