A doubt regarding change of variables in Double Integrals. So I have recently learnt to calculate "Double Integrals" and I have a doubt on change of variables. Lets say we are evaluating a double integral of the form I = $\int\int _R F(x,y) dx dy $ , where $F(x,y)$ is an integrable smooth function and $R$ is the region of integration in the $xy$ plane. Now I want to convert to polar co-ordinates $x(r,\theta)=r \cos(\theta)$ and $y(r,\theta) = r \sin(\theta)$. Now the region $R$ in the $xy$ palne can be suitably converted to some region $R'$ in the $r \theta$ plane. The function $F(x,y)$ may take the form of some smooth and integrable $G(r,\theta)$. My main doubt is with the change of the area element. According to my textbook, $dx dy = \frac{\partial (x,y)}{\partial(r,\theta)} = rdrd\theta$. So the area element in the $r\theta$ plane is $rdr d\theta$.
Now I want to experiment a little bit with my understanding : I write $dx = d(r\cos\theta) = -r\sin\theta d\theta + \cos \theta dr$  and $dy = d(r\sin\theta)=r\cos\theta d\theta + \sin\theta dr $. And now I multiply the two expression for $dx$ and $dy$ : $dxdy= -r^2\sin\theta\cos\theta d\theta d\theta - r \sin^2 \theta drd\theta + r \cos^2 \theta dr d\theta + \sin \theta \cos\theta dr dr$. Why isn't this expression equal to $rdrd\theta$ ? The $d\theta d\theta$ and $dr dr$ terms look very odd to me. Am I commiting some abuse of notation or something ? Also there are no higher order terms like $dr dr d\theta$ which I can ignore. I am totally confused with all this. Can someone explain this to me a bit clearly? Thanks in Advance :).
EDIT 1 : I know differentials are not the "usual" kind of algebraic quantities and we can't just add and multiply them casually. I am just multiplying $dx$ with $dy$ with the spirit of multiplying some $\delta x$ and $\delta y$ and letting both $\delta x$ and $\delta y$ tend to zero.
EDIT 2 : Is there a branch of Mathematics where addition and multiplication of differentials is very common? Is there something like differential of a differential ?
 A: Side Question

Is there a branch of Mathematics where addition and multiplication of differentials is very common? Is there something like differential of a differential ?

There are lots of different contexts/formalizations with involved manipulation of differentials. Trying to give an overview of all of them would be a big undertaking worthy of a whole separate question. Very briefly: The most common thing that is taught/that I would expect a mathematician to think of is the theory of "differential forms". But there are also many other relevant things like "geometric calculus", "Smooth Infinitesimal Analysis", Bartlett and Khurshudyan's approach to manipulating higher order differentials (arXiv link), debatably: Robinson's approach to nonstandard analysis, and probably a couple other things I'm forgetting or haven't encountered.
Main Question

Why isn't this expression equal to $rdrd\theta$ ?

Disclaimer
This is almost entirely copied from my answer to the similar question Why can't we convert the area element $dA$ to polar by multiplying the polar expressions for $dx$ and $dy$?.
Intro
There are two main types of ways to think about things like $\mathrm{d}x\mathrm{d}y$ in multivariable calculus, and we often switch between them depending on the context. (This clarification was inspired in part by Terry Tao's preprint on “differential forms and integration”.) $\mathrm{d}x$ can either act kind of like a number, or act kind of like a vector.
For the “number” interpretation, there are things like limit arguments or infinitesimals in nonstandard analysis to make things rigorous. For the “vector” interpretation, there are things like “differential forms” or “geometric calculus” to make things rigorous. But I'm going to gloss over those details because there are many ways to make things formal, and the exact choices don't affect the intuition here.
Numbers
One way to think about things is that $\mathrm{d}x$ and $\mathrm{d}y$ are in some way like tiny positive numbers representing the width and length of a tiny rectangle, so that $\mathrm{d}x\mathrm{d}y$ is the area of a tiny rectangle. Then when we write something like $\iint f\left(x,y\right)\,\mathrm{d}x\mathrm{d}y$ or $\iint g\left(r,\theta\right)\,\mathrm{d}r\mathrm{d}\theta$, we just add up the signed volumes (in case $f$ or $g$ is negative) of thin rectangular prisms with cross-sectional area represented by $\mathrm{d}x\mathrm{d}y$ or $\mathrm{d}r\mathrm{d}\theta$.
Under this interpretation, $\mathrm{d}x=\mathrm d(r\cos\theta)=\cos\theta\mathrm{d}r-r\sin\theta\mathrm{d}\theta$ doesn't make too much sense. For example, if $\theta=\pi/2$, then we would have $\mathrm{d}x=-r\mathrm{d}\theta$, so that $\mathrm{d}\theta$ and $\mathrm{d}x$ couldn't both represent positive lengths. But we can still understand the relationship between the areas $\mathrm{d}x\mathrm{d}y$ and $\mathrm{d}r\mathrm{d}\theta$ with arguments like the geometric one in this answer by Mike Spivey.
Vectors
The other way to think about things is that $\mathrm{d}x$ and $\mathrm{d}y$ are in some way like tiny vectors whose direction we care about, and this leads to a slightly different discussion. To emphasize this vector idea, I will use some nonstandard notation. Let's write $\overrightarrow{\mathrm{d}x}=\left\langle \Delta x,0,0\right\rangle$  for some positive $\Delta x$, and $\overrightarrow{\mathrm{d}y}=\left\langle 0,\Delta y,0\right\rangle$  for some positive $\Delta y$. So $\overrightarrow{\mathrm{d}x}$ points to the right in the $xy$-plane and $\overrightarrow{\mathrm{d}y}$ points “up” in the $xy$-plane. Then the area of the little rectangle they make is $\left\Vert \overrightarrow{\mathrm{d}x}\times\overrightarrow{\mathrm{d}y}\right\Vert =\Delta x\Delta y$.
However, now that we have vectors, we could choose to care about the orientation. When we think about a usual integral like $\int_{\left[a,b\right]}f\left(x\right)\,\mathrm{d}x$ when $f$ is negative, we decide to count that area in a negative way. Similarly, we could count the “signed area” of the oriented rectangle with sides $\overrightarrow{\mathrm{d}x},\overrightarrow{\mathrm{d}y}$ as the positive value $\Delta x\Delta y$ (which is just the $z$-component of $\overrightarrow{\mathrm{d}x}\times\overrightarrow{\mathrm{d}y}$). But the signed area of the oriented rectangle with sides $\overrightarrow{\mathrm{d}y},\overrightarrow{\mathrm{d}x}$ would be the negative value $-\Delta x\Delta y$ (which is just the $z$-component of $\overrightarrow{\mathrm{d}y}\times\overrightarrow{\mathrm{d}x}$).
Then when we write something like $\iint f\left(x,y\right)\,\mathrm{d}x\mathrm{d}y$ or $\iint g\left(r,\theta\right)\,\mathrm{d}r\mathrm{d}\theta$, we could (if we chose) care about the order of things, with something like an implied $\overrightarrow{\mathrm{d}x}\times\overrightarrow{\mathrm{d}y}$ or $\overrightarrow{\mathrm{d}r}\times\overrightarrow{\mathrm{d}\theta}$ in our heads, if not always our notation.
When using things like vectors, then something like $\overrightarrow{\mathrm{d}x}=\cos\theta\overrightarrow{\mathrm{d}r}-r\sin\theta\overrightarrow{\mathrm{d}\theta}$ makes a lot of sense. $\overrightarrow{\mathrm{d}r}$ points away from the origin, and $\overrightarrow{\mathrm{d}\theta}$ points perpendicularly to $\overrightarrow{\mathrm{d}r}$ in the counterclockwise way (so that $\overrightarrow{\mathrm{d}r}\times\overrightarrow{\mathrm{d}\theta}$ points in the same positive-$z$ direction as $\overrightarrow{\mathrm{d}x}\times\overrightarrow{\mathrm{d}y}$).
Finally, we can correct your calculation:
$$\overrightarrow{\mathrm{d}x}\times\overrightarrow{\mathrm{d}y}=\left(\cos\theta\overrightarrow{\mathrm{d}r}-r\sin\theta\overrightarrow{\mathrm{d}\theta}\right)\times\left(\sin\theta\overrightarrow{\mathrm{d}r}+r\cos\theta\overrightarrow{\mathrm{d}\theta}\right)$$
$$=\cos\theta\sin\theta\overrightarrow{\mathrm{d}r}\times\overrightarrow{\mathrm{d}r}-r^{2}\sin\theta\cos\theta\overrightarrow{\mathrm{d}\theta}\times\overrightarrow{\mathrm{d}\theta}+r\cos^{2}\theta\overrightarrow{\mathrm{d}r}\times\overrightarrow{\mathrm{d}\theta}-r\sin^{2}\theta\overrightarrow{\mathrm{d}\theta}\times\overrightarrow{\mathrm{d}r}$$
$$=\overrightarrow{0}-\overrightarrow{0}+r\cos^{2}\theta\overrightarrow{\mathrm{d}r}\times\overrightarrow{\mathrm{d}\theta}-r\sin^{2}\theta\left(-\overrightarrow{\mathrm{d}r}\times\overrightarrow{\mathrm{d}\theta}\right)=\boxed{r\,\overrightarrow{\mathrm{d}r}\times\overrightarrow{\mathrm{d}\theta}}$$
A: My understanding of this is that in integration the area element which we denote $dA = dxdy$ is in fact a differential form as Miguel mentioned. A differential 2-form measures a certain amount of area, and is made up of differential 1-forms (which measure length) through the wedge product $\wedge$. The wedge product provides natural measures for a space through what's called the exterior algebra of the space. For example, given two vectors $\vec{a}$ and $\vec{b}$ the area of the parallelopiped they span is given by the scalar part of $\vec{a}\wedge \vec{b}$.
That is, when integrating over an oriented surface the area 2-form is $dA = dx\wedge dy.$ Because it describes oriented areas the wedge product is anticommutative, i.e. $dx\wedge dy = -dy\wedge dx$. Then when changing variables, you encounter the quadratic forms $d\theta d\theta$ for example. Interpreting the differentials here as measures (like in differential geometry), when you substitute your change of variables into $dx\wedge dy$ (called a pull-back if you want to get fancy), the terms with repeated differential 1-forms must vanish: $$d\theta \wedge d\theta = -d\theta\wedge d\theta = 0.$$
Differential forms in integration should be distinguished from Riemannian notation like $ds^2$, which really denotes a tensor product. A good elementary book on differential forms is "Differential Forms: Theory and Practice" by Weintraub. Also check out the page on exterior algebra on Wikipedia.
A: The change of variables transforms a small rectangular region $R$ with sides $\Delta x$ and $\Delta y$ in a region $R’$, which is not necessarily rectangular, so an area of $R’$ can differ from a product $\Delta r\cdot\Delta\theta$ of its “dimensions” $\Delta r$ and $\Delta\theta$.
For the case of polar coordinates, $R’$ is bounded by circles of radii $r$ and $r+\Delta r$ and radial lines directed from the origin at angles $\theta$ and $\theta+\Delta\theta $, so the area of $R’$ equals $$\frac 12\Delta\theta((r+\Delta r)^2-r^2)= r\Delta r\Delta\theta+\frac 12(\Delta r)^2\Delta\theta.$$
