Why is the graph of the implicit relation $\sin(x^2+y^2) = \cos(xy)$ so cool?

Here are some observations to try to explain what is going on.

First of all, why does there appear to be concentric repetition of the image expanding outwards? This is caused by the periodicity of sine and cosine. Each pair of crossing ovals happens periodically, but not at a constant period. Since the input is quadratic in the variables, the period is shrinking as we move outward, making the repetition more frequent. I tackle this later on by looking at radial slices and equations of the ovals.

Secondly, why are there pairs of intersecting oval shapes? As one quickly guesses, the graph has 90 degree rotational symmetry. The transformation $(x,y)\to(y,-x)$ carries solutions of the equation to other solutions. Obviously $x^2+y^2=y^2+(-x)^2$, and $\cos(xy)=\cos(y(-x))$ since cosine is an even function.

One can try to get a grip on the solutions closest to the origin and the subsequent repetition by considering radial slices of the plane. The first two to investigate would be the one with $x=0$ (Which covers the y=0 case, after symmetry.) you would be solving the equation $\sin(y^2)=1$, which happens whenever $y^2$ hits $\pi/2 + n2\pi$.

The solutions on the line $x=y$, for example, would have to satisfy $\sin(2x^2)=\cos(x^2)$. Actually, you an work out what the equation should be for $y=kx$ for any constant k.

Added for completeness:

Are they actually ellipeses? Yes! Achille Hui provided a derivation of why the ovals are really ellipses in comments to the original post. I'm echoing part of it here to make it a little more visible and durable.

Using these equations:

$$\begin{align} \sin(A) = \cos(B) \iff & B = \pm ( \frac{\pi}{2} - A ) + 2m\pi \quad\text{ for some }\; m \in \mathbb{Z}\\ \iff & A \pm B = (2n+\frac{1}{2})\pi\quad\text{ for some }\; n \in \mathbb{Z} \end{align}$$

substitution gives the form of the ellipses:

$$x^2+y^2 \pm xy = (2n + \frac12)\pi$$

Work of my own expanding on this:

If we had a family of equations of concentric circles $x^2+y^2=an+b$ for constants $a,b$ and natural number $n$, then we would recognize that the square of the radius varies with $n$, so the radius should vary as $\sqrt{n}$, and we would expect the radii to be "slowing down" since the values of $\sqrt{n}$ get closer together as $n$ increases.

Rotating this family of equations by $\pi/4$, we can eliminate the $xy$ term and think about the ellipses as basic ones we do in calculus. Here's the rotated version using $+$ for the $\pm$ part of the equation above:

$$ \frac{3x^2}{2}+\frac{y^2}{2}=2\pi n +\frac{\pi}{2}\\ 3x^2+y^2=4\pi n +\pi\\ \frac{x^2}{(\frac{4\pi n+\pi}{3})}+\frac{y^2}{4\pi n+\pi}=1 $$

So here we can see that the major and minor principal axes vary similarly to $\sqrt{n}$, and moreover that the major principal axis is $\sqrt{3}$ times that of the minor axis. If you do the same trick with the "$-$" for $\pm$, you get the other family of ellipses rotated $\pi/2$ from these ellipses.