A conceptual reason for why the Jacobian of a rotation by a changing angle is $1$? Consider
$$f(x,y)=(x\cos r^2+y\sin r^2, y\cos r^2-x\sin r^2)\qquad\text{with }r=\sqrt{x^2+y^2},$$
as a map $\mathbb{R}^2 \to \mathbb{R}^2$. Geometrically, $f(x,y)$ is obtained from the vector $(x,y)$ by rotating it by angle $r^2$. (we rotate the more we move away from the origin).

Wolframe Alpha claim that $\det (df)=1$ holds identically. Is there an elegant rigorous way of seeing this, without too much computations?

Should this be "obvious" in retrospect? I was a bit surprised by this result... and direct computation is not so nice to do by hand (although tractable).
I thought using the chain rule (treating the angle as a function of $x,y$) but got nowhere.
Edit:
I agree that roughly speaking, near a point $p$, this function is like a standard rotation by a fixed angle $|p|^2$. However, I do not consider this a rigorous explanation. Indeed, this vague intuition is still with us when we replace $r^2=x^2+y^2$ by $x^4+y^4$, but then the Jacobian is non-constant. 

So, this property doesn't even hold for smooth radially symmetric angle function $\theta(x,y)$.

Is this just a coincidence then? Can we charavterize the angle functions (or at least radially symmetric ones) which satisfy this?
(Wolframe says the Jacobian remains $1$ when we replace $r^2$ by $r$ or $r^4$. Perhaps this remains true for any power of $r$?)
 A: Any symplectomorphism has Jacobian with determinant one. It is relatively easy to see that the Hamiltonian $H_n(x)=\Vert x\Vert^n$ for $n \geq 2$ has as flows the "rotations" you mention. (I.e., by powers of $r$.) 
Indeed, $\nabla H_n(x)=n\Vert x\Vert^{n-2}x,$ thus letting $\Phi_t^n$ denote the Hamiltonian flow of $H_n$ at time $t$, we have that $\Phi^4_{1/4}$ is the rotation you mainly want in the question. For the generalized question of rotations by powers $r^n$ of $r$, just take $\Phi_{1/(n+2)}^{n+2}$.
This also generalizes for the case that Jyrki mentions, i.e. $\theta \circ r$. For this case, take $F$ to be an antiderivative for $f(r)=\theta(r) \cdot r$ and pick the Hamiltonian $H=F \circ \Vert \cdot \Vert$. Then 
$$\nabla H(x)=F'(\Vert x \Vert)\cdot \frac{x}{\Vert x \Vert}=\theta(\Vert x \Vert) \cdot \Vert x \Vert \cdot \frac{x}{\Vert x \Vert}=\theta(\Vert x \Vert) \cdot x.$$
Then the flow at time $1$ is the desired rotation.
It may not be clear why this is conceptual and natural without some background in symplectic geometry, so some words about that:


*

*It is a known fact disguised in several different ways, among which is the Liouville's theorem (phrased in Wikipedia in a physics POV), that Hamiltonian flows preserve volume. This is of course intrinsically related to the determinant of the Jacobian being one. So we search for a way to construct your map as a Hamiltonian flow.

*Hamiltonian flows stay confined in the level-sets and have velocity equal to the norm of the gradient. So whatever Hamiltonian we will use to guess, since your map is a rotation, it must be constant in circles. Luckily, this also tells us that if we can find such Hamiltonian that is constant in circles, then the flows will be rotations by some angle in each circle. This is because the gradient will have constant norm in each level set. (A fact which is intuitively clear since the level sets are homogeneously spreading out, but backed up by the computation I gave above.) 

*So we just have to find a Hamiltonian which has its gradient at every point $x$ with norm equal to $\theta\cdot \Vert x\Vert $, where $\theta$ is the angle we want to rotate in that radius value for $x$. In your case we want it to be equal to $\Vert x \Vert ^2 \cdot \Vert x \Vert$. Thus it is natural to choose $H(x)=\frac{1}{4}\Vert x\Vert ^4$, or simply $H(x)=\Vert x \Vert^4$ and pick the time $1/4$ of the flow.

A: 
The claim holds whenever $\theta(x,y)=\theta(r)$, $r=\sqrt{x^2+y^2}$.

One way to think about this would be to view your function as a composition of two differentiable functions, $f=g\circ h$: 
$$
\begin{aligned}
h:\Bbb{R}^2\to\Bbb{R}^3,\ &(x,y)\mapsto (x,y,\theta(x,y))\\
g:\Bbb{R}^3\to\Bbb{R}^2,\ &(x,y,\theta)\mapsto R(\theta)(x,y)=(x\cos\theta -y\sin\theta ,x\sin\theta +y\cos\theta ),
\end{aligned}
$$
where $R(\theta)$ is the rotation of the plane about the origin by the angle $\theta$.
The Jacobian (or the linear/Frechét derivative) of $h$ is
$$
Dh(x,y)=\left(
\begin{array}{cc}
1&0\\
0&1\\
\partial\theta/\partial x&\partial\theta/\partial y
\end{array}
\right)
$$
whereas the Jacobian of $g$ is
$$
Dg(x,y,\theta)=\left(
\begin{array}{ccc}
\cos\theta&-\sin\theta&-x\sin\theta-y\cos\theta\\
\sin\theta&\cos\theta&x\cos\theta-y\sin\theta
\end{array}
\right).
$$
You are interested in the determinant of the product of these two matrices that a calculation (heavily using $\cos^2\theta+\sin^2\theta=1$) reveals to be
$$
\det(Dg\cdot Dh)=1+x\frac{\partial\theta}{\partial y}-y\frac{\partial\theta}{\partial x}.
$$
But we have $\partial r/\partial x=x/r$, $\partial r/\partial y=y/r$, so by the chain rule $\theta_x=(x/r)\theta_r$, $\theta_y=(y/r)\theta_r$, implying that
$$
x\frac{\partial\theta}{\partial y}-y\frac{\partial\theta}{\partial x}=\frac{\partial\theta}{\partial r}\left(x\cdot\frac yr-y\cdot \frac xr\right)=0.
$$
Therefore $\det Dg\cdot Dh=1$.
A: I would say it is not completely obvious, but it is believable. The value of the Jacobian at $x$ is the factor by which area is distorted locally near $x$. Roughly speaking, near $x$, this function is like a standard rotation by a fixed angle $|x|$. The fact that this isn't true globally seems to be a negligible error, but I don't think it's obvious it was going to work out like that.
A: Partition the plane into very thin annular sets $A_k: \>r_{k-1}\leq\sqrt{x^2+y^2}< r_k$ $(k\in {\mathbb N}_{\geq1})$. The map $f$ (approximately) rotates each $A_k$ by an angle $\phi_k$ around the origin, hence is (approximately) an isometry on $A_k$. This implies $\mu\bigl(f(X)\bigr)=\mu(X)$ for  any $X\subset A_k$. For an arbitrary set $B\subset{\mathbb R}^2$ consider the annular subsets $B_k:=B\cap A_k\subset A_k$ and obtain
$$\mu\bigl(f(B)\bigr)=\sum_{k=1}^\infty \mu\bigl(f(B_k)\bigr)=\sum_{k=1}^\infty \mu(B_k)=\mu(B)\ .$$
This shows that $f$ is globally area preserving. (But I'd say that using the chain rule on $df$ is simpler.)
