Polar Coordinate Transformation - Motivation I am trying to work out the reason why the integral int_(-infty)^(infty)int_(-infty)^(infty)(e^-(x^2+y^2))),dx,dy is, in polar coordinates, int_(infty)^(infty)(e^(-r^2)) ,r,dr,d theta As I understand it, a polar coordinate transformation involves the following substitution: (x,y) -> (r cos theta,r sin theta) This would imply that -(x^2+y^2) = -((r cos theta)^2+(r sin theta)^2) = -r^2((cos theta)^2+(sin theta)^2) = -r^2 which gets us this far int_-(infty)^(\infty)int_(-infty)^(infty)e^(-r^2),dx,dy To motivate why dx,dy = r,dr,d theta

No, this mere substitution argument cannot work, and you haven't made an algebraic mistake at all. This goes precisely to show you that these differentials $dx,dy,dr,d\mathrm{\theta <!--\; \theta \; -->}$ are NOT real numbers. So, there is absolutely no reason you should treat them as real numbers and perform the usual arithmetic of real numbers with them. Therefore, the better question to ask yourself is why would you even expect this naive "substitution" to work? Sure, in single variable calculus, many such arguments just happen to "work", but of course this is not at all a good reason to expect it to keep working.
So, how do we treat such integrals? The key result is the multi-dimensional change of variables formula. The proof of this formula is very technical, and not easy at all, but the rough idea is I think straight forward. A very rough way of saying it is that in cartesian coordinates, the function you're integrating is $f(x,y)=e^{-<!--\; -\; -->(x^2+y^2)}$, so when you change to polar coordinates, you have modified your function to
$\begin{array}{r}F(r,\mathrm{\theta <!--\; \theta \; -->})=f(r\cos <!--\; \; -->\mathrm{\theta <!--\; \theta \; -->},r\sin <!--\; \; -->\mathrm{\theta <!--\; \theta \; -->})=e^{-<!--\; -\; -->r^2}.\end{array}$
If we modify the function, then of course we have to modify how we measure the areas. Now, in cartesian coordinates, a little area element is described as dx dy, which you intuitively think of as a small box of size dx and dy. But if you now want to integrate in polar coordinates, you have to think of what does a little area element look like. There are several geometrical justifications available for this on this site and others, and the correct answer is $rdrd\mathrm{\theta <!--\; \theta \; -->}$. That factor of r encodes how the area changes as you measure things in cartesian coordinates vs how you measure things in polar coordinates.
Once again, there is absolutely no reason to think the correct area element should be obtained by plugging in $dx=\cos <!--\; \; -->\mathrm{\theta <!--\; \theta \; -->}dr+\dots <!--\; \dots \; -->$… and likewise for dy, into the "product" (however it is defined) dx dy. Nor should you expect that "because $dA=dxdy$ in cartesian coordinates, therefore $dA=drd\mathrm{\theta <!--\; \theta \; -->}$ in polar coordinates". None of these is a proper way of determining the area element in a different coordinate system. It's not only wrong, but there's also no reason to expect it to work.
In general, the correct scale factor which accounts for the difference in area measurements in different coordinates is the absolute value of the determinant of the derivative of the coordinate transformation. This is a non-trivial issue, and there are several answers on this site which explain in more detail why this is the right "conversion factor" (and also entire chapters of textbooks devoted to this formula and its geometric significance).
Hopefully now that you're convinced that naive substitution is not at all even a reasonable thing to do, let's see how one can arrive at the answer mechanically. It's very simple, and even very memorable:
$\begin{array}{rl}dxdy& =|det{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}(x,y)}{\mathrm{\partial <!--\; \partial \; -->}(r,\mathrm{\theta <!--\; \theta \; -->})}}|drd\mathrm{\theta <!--\; \theta \; -->}\end{array}$
Where that funny derivative is a $2\times <!--\; \times \; -->2$ matrix of partial derivatives, which in this case is:
$\begin{array}{rl}{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}(x,y)}{\mathrm{\partial <!--\; \partial \; -->}(r,\mathrm{\theta <!--\; \theta \; -->})}}& =\left(\begin{array}{cc}\cos <!--\; \; -->\mathrm{\theta <!--\; \theta \; -->}& -<!--\; -\; -->r\sin <!--\; \; -->\mathrm{\theta <!--\; \theta \; -->}\\ \sin <!--\; \; -->\mathrm{\theta <!--\; \theta \; -->}& r\cos <!--\; \; -->\mathrm{\theta <!--\; \theta \; -->}\end{array}\right)\end{array}$
The absolute value of the determinant is |r|, but this is simply r, since $r\ge <!--\; \ge \; -->0$. Therefore, $dxdy=rdrd\mathrm{\theta <!--\; \theta \; -->}$.
Another way to approach the subject is by the means of differential forms. The differential forms machinery pretty much encodes in its definition the determinant of the Jacobian matrix. Now, I won't detail the entire theory of differential forms, but here's the "rules" (again if you want to know more, there's entire books on the subject).
Here, things like $dx,dy,drd\mathrm{\varphi <!--\; \varphi \; -->},d\mathrm{\theta <!--\; \theta \; -->}$ etc are no longer real numbers/infinitesimally small but non-zero numbers (whatever that is supposed to mean in typical analysis). What are they? They're the section of the cotangent bundle. Ok, that's probably just some mumbo jumbo right now, but here's how we calculate with them.
For the purposes of integration, we no longer write dxdy. Instead, we put a little $\wedge <!--\; \wedge \; -->$, like dx$\wedge <!--\; \wedge \; -->$dy; this is called the wedge product. The key property it has is that $dx\wedge <!--\; \wedge \; -->dy=-<!--\; -\; -->dy\wedge <!--\; \wedge \; -->dx$; so it alternates sign whenever you flip two of them. In particular, $dx\wedge <!--\; \wedge \; -->dx=-<!--\; -\; -->dx\wedge <!--\; \wedge \; -->dx$ (I flipped them) and therefore $dx\wedge <!--\; \wedge \; -->dx=0$. In other words whenever you wedge the same thing, it vanishes. Finally, all the other rules, like associativity, distributivity "expanding brackets" etc all work as normal. Now we can compute:
$\begin{array}{rl}dx\wedge <!--\; \wedge \; -->dy& =(\cos <!--\; \; -->\mathrm{\theta <!--\; \theta \; -->}dr-<!--\; -\; -->r\sin <!--\; \; -->\mathrm{\theta <!--\; \theta \; -->}d\mathrm{\theta <!--\; \theta \; -->})\wedge <!--\; \wedge \; -->(\sin <!--\; \; -->\mathrm{\theta <!--\; \theta \; -->}dr+r\cos <!--\; \; -->\mathrm{\theta <!--\; \theta \; -->}d\mathrm{\theta <!--\; \theta \; -->})\\ & =(\cos <!--\; \; -->\mathrm{\theta <!--\; \theta \; -->}\sin <!--\; \; -->\mathrm{\theta <!--\; \theta \; -->})\underset{=0}{\underset{⏟<!--\; ⏟\; -->}{dr\wedge <!--\; \wedge \; -->dr}}+(r{\cos }^2<!--\; \; -->\mathrm{\theta <!--\; \theta \; -->})dr\wedge <!--\; \wedge \; -->d\mathrm{\theta <!--\; \theta \; -->}+(-<!--\; -\; -->r{\sin }^2<!--\; \; -->\mathrm{\theta <!--\; \theta \; -->})\underset{=-<!--\; -\; -->dr\wedge <!--\; \wedge \; -->d\mathrm{\theta <!--\; \theta \; -->}}{\underset{⏟<!--\; ⏟\; -->}{d\mathrm{\theta <!--\; \theta \; -->}\wedge <!--\; \wedge \; -->dr}}+(-<!--\; -\; -->r^2\sin <!--\; \; -->\mathrm{\theta <!--\; \theta \; -->}\cos <!--\; \; -->\mathrm{\theta <!--\; \theta \; -->})\underset{=0}{\underset{⏟<!--\; ⏟\; -->}{d\mathrm{\theta <!--\; \theta \; -->}\wedge <!--\; \wedge \; -->d\mathrm{\theta <!--\; \theta \; -->}}}\\ & =rdr\wedge <!--\; \wedge \; -->d\mathrm{\theta <!--\; \theta \; -->}.\end{array}$
It is a general fact that the factor you get in front will always be the determinant of the Jacobian matrix of the coordinate change. The main reason for this is that there is a deep connection between determinants, alternating objects like differential forms, and volumes of parallelepiped (which you may have learnt in linear algebra).
These relations are very general, and work for any coordinate system. For exmaple, if you consider parabolic coordinates $\mathrm{\xi <!--\; \xi \; -->},\mathrm{\eta <!--\; \eta \; -->}$, whose relation to cartesian coordinates is $x=\mathrm{\xi <!--\; \xi \; -->}\mathrm{\eta <!--\; \eta \; -->}$ and $y={\displaystyle \frac{1}{2}}({\mathrm{\eta <!--\; \eta \; -->}}^2-<!--\; -\; -->{\mathrm{\xi <!--\; \xi \; -->}}^2)$, you can work through a similar mechanical procedure to find
$\begin{array}{rl}dx\wedge <!--\; \wedge \; -->dy& =({\mathrm{\xi <!--\; \xi \; -->}}^2+{\mathrm{\eta <!--\; \eta \; -->}}^2)d\mathrm{\xi <!--\; \xi \; -->}\wedge <!--\; \wedge \; -->d\mathrm{\eta <!--\; \eta \; -->}.\end{array}$
Final Remarks.
There are several subtleties which I have had to gloss over, like the concept of orientation, why differential forms are "suitable objects to integrate", etc. But hopefully this gives you an idea for why naive substitution doesn't (and shouldn't be expected to) work, and that one needs an appropriate scale factor to treat the areas; also hopefully there have been enough "mechanical" rules to help (of course you should eventually try to understand them at a deeper level).
One more thing I'd like to highlight is that while the differential forms approach gives a very quick and easy mechanical approach to the change of volume factor, and while computing determinants algebraically is very straight forward, it is very important to keep in mind the geometry behind these formulas. These Riemann integrals are all about areas of various (simple) figures, so it is absolutely essential to atleast somewhat connect the algebraic significance of the determinant to its geometric significance regarding areas/volumes.