For $n=1$ or $2$, there is no $u\in C^1(D_\rho)$ for any $\rho>0$ that satisfies $|\nabla u|^2 = (xy)^{2n}$. (Note that $f=(xy)^2$ has a degenerate minimum at $(0,0)$, so $n=1$ should be allowed in this discussion.) Meanwhile, for $n\ge 3$, there do exist $u\in C^1(\mathbb{R}^2)$ that satisfy $|\nabla u|^2 = (xy)^{2n}$.
The above non-existence result is special to the case $f = (xy)^{2n}$. For example, in the case $f=(xy)^{2n}(x^2{+}y^2)$, which has an even more degenerate minimum at $(0,0)$, there is a real-analytic, global solution $u(x,y) = (xy)^{n+1}/(n{+}1)$.
A few preliminaries before sketching the main argument are in order.
First, note that the disk radius $\rho>0$ actually plays no role in the problem. If $u\in C^1(D_\rho)$ were to satisfy $|\nabla u|^2 = (xy)^{2n}$, then for any $r>0$ the scaled function $\tilde u(x,y) = r^{-(2n+1)}\,u(rx,ry)$ would satisfy $\tilde u\in C^1(D_{\rho/r})$ and $|\nabla \tilde u|^2 = (xy)^{2n}$. Thus, one can assume that $\rho$ be arbitrarily large.
Second, one can assume, by adding a constant to $u$, that $u(0,0)=0$, so I will assume this normalization made henceforth. Then, the obvious integral inequality arising from $|\nabla u| = |xy|^n$ and the Cauchy-Schwartz inequality would imply that
$$
|u(x,y)| \le \frac{|xy|^n\sqrt{x^2+y^2}}{(2n{+}1)}.
$$
In particular, $u$ would vanish to order $2n{+}1$ at $(0,0)$ and would satisfy $u(x,0) = u(0,y) = 0$. It follows from this that $u$ could not be of differentiability class $C^{2n+2}$, since, if it were, the limit function
$$
p(x,y) = \lim_{r\to0} \frac{u(rx,ry)}{r^{2n+1}}
$$
would exist and be a polynomial homogeneous of degree $2n{+}1$ that satisfied $|\nabla p|^2 = (xy)^{2n}$, and it is easy to show that there is no such polynomial. However, as will be seen, when $n\ge3$, there exists a $u\in C^{n-1}(\mathbb{R}^2)$ satisfying $|\nabla u|^2 = (xy)^{2n}$ and the homogeneity condition $u(rx,ry) = |r|^{2n+1}\,u(x,y)$ for all $r$. This $u$ is real-analytic away from the lines $x\pm y = 0$ but fails to be $C^n$ on these two lines.
From now on, I will assume that $u\in C^1(D_\rho)$ (with $\rho>>0$ as large as necessary for the argument) satisfies $u(0,0)=0$
and $|\nabla u|^2 = (xy)^{2n}$. In particular, as the OP points out, $u$ satisfies the eikonal equation $|\nabla^g u|^2 = 1$ for the singular 'metric' $g = (xy)^{2n}(\mathrm{d}x^2+\mathrm{d}y^2)$, so that the gradient flow lines of $\nabla^g u$ are $g$-geodesics in the four quadrants of the $xy$-plane where $xy\not=0$.
Now, the metric $g$ has some interesting properties: First, it is homogeneous of degree $2n{+}2$, so that its family of geodesics is preserved under the scaling homothety, and, moreover, it is invariant under the discrete symmetries $(x,y)\to(-x,y)$, $(x,y)\to(x,-y)$, and $(x,y)\to(y,x)$. Consequently, it suffices to study the behavior of its geodesics in the 'first' quadrant, where $x>0$ and $y>0$, and it is immediate that the lines $y\pm x = 0$ are geodesics in the quadrants. Let the ray $\{(x,x)\ |\ x>0\}$ in the first quadrant be known as the fundamental geodesic. Of course, $g$ is not complete in the first quadrant, as the two boundary rays can be reached from, say, $(1,1)$, by curves of finite $g$-length.
Now, computation shows that the Gauss curvature of $g$ is $K = n(x^2{+}y^2)/(xy)^{2n+2}>0$, which suggests that nearly all of the geodesics of $g$ will avoid going to the singular boundary where $xy=0$, and, indeed, this turns out to be the case (see below). To parametrize the geodesics, it turns out to be convenient to use a parameter $t$ other than arc length. A curve $\bigl(x(t),y(t)\bigr)$ in the first quadrant parametrizes a $g$-geodesic when there is a function $\phi(t)$ satisfying the ODE system
$$
\dot x = xy\cos\phi,
\quad \dot y = xy\sin\phi,
\quad \dot\phi = n\,(x\cos\phi-y\sin\phi),\tag1
$$
and every $g$-geodesic in the first quadrant has such a parametrization, unique up to replacing $t$ by $t+t_0$ for some constant $t_0$.
In this case, arclength $s(t)$ along the geodesic satisfies $\dot s = (xy)^{n+1}$. [The advantage of writing the geodesic equations this way is that they extend smoothly across the singular locus $xy=0$.] Note that these equations are invariant under the homothetical scaling $(t,x,y,\phi)\to(t/r,rx,ry,\phi)$. Because of the scaling symmetry of the equations, one can extract a 2D phase portrait that makes clear the behavior of the geodesics as follows: Let $x+iy = \mathrm{e}^{u+iv}$. Then the above equations become (after a change of independent variable)
$$
u' = \cos(\phi{-}v)\,\cos v\sin v,\qquad
v' = \sin(\phi{-}v)\,\cos v\sin v,\qquad
\phi' = n\,\cos(\phi{+}v).\tag2
$$
One can now draw the $v\phi$-phase portrait, concentrating on the strip $0\le v\le \pi/2$, which represents the first quadrant in the $xy$-plane, and bearing in mind that these equations are invariant under the involution $(v,\phi)\to(\tfrac12\pi{-}v,\tfrac12\pi{-}\phi)$ and reverse under $(v,\phi)\to (v,\phi{+}\pi)$.
There are a sink at $S_- = (v,\phi)=(\pi/4,\pi/4)$, a source at $S_+ = (v,\phi)=(\pi/4,-3\pi/4)$, and saddles at $S_1 = (v,\phi)=(0,\pi/2)$, $S_2 = (v,\phi)=(\pi/2,0)$, $S_3=(v,\phi)=(0,-\pi/2)$, and $S_4=(\pi/2,-\pi)$. In addition to the 'trivial' separatrices that make up the boundary lines $v=0$ and $v=\tfrac12\pi$, there are four 'non-trivial' separatrices: $L_1$ leaving $S_1$ and going to $S_-$, $L_2$ leaving $S_2$ and going to $S_-$, $L_3$ leaving $S_+$ and going to $S_3$, and $L_4$ leaving $S_+$ and going to $S_4$.
Here is where the difference between the cases $n=1,2$ and the cases $n\ge 3$ becomes evident. When $n\le 2$, the fixed points $S_\pm$ are spiral, i.e., the linearization of the flow at these two points has eigenvalues that are non-real (and complex conjuate), while, when $n\ge 3$, the linearizations of the flow at these two points have distinct real eigenvalues whose ratio is a real number strictly between $n{-}2$ and $n{-}1$.
One then finds that, when $n\ge 3$, the union of the two separatrices $L_1$ and $L_2$, together with their endpoints $S_1$, $S_2$, and $S_-$ is the graph $\phi = f_n(v)$ of a function $f_n:[0,\tfrac12\pi]\to[0,\tfrac12\pi]$ that is real-analytic except at the midpoint $v=\tfrac14\pi$, where it is $C^{n-2}$. One then shows that the corresponding $g$-geodesics in the first quadrant define a $C^{n-2}$ foliation by geodesics that meet the boundary rays $x=0$ and $y=0$ orthogonally. Moreover, by reflecting this foliation across the lines $x=0$ and $y=0$, one can construct a foliation $\mathcal{F}$ of $\mathbb{R}^2$ minus the origin by $g$-geodesics that is real-analytic except along the lines $x\pm y=0$, where it is $C^{n-2}$. It then follows easily that there is a unique $C^{n-1}$ function $u$ that vanishes on the axes $x=0$ and $y=0$, satisfies $|\nabla u|^2=(xy)^{2n}$ and, away from the axes, the gradient lines of $u$ are the leaves of the foliation $\mathcal{F}$.
Meanwhile, when $n\le 2$, the spiral nature of the two fixed points $S_\pm$ leads to an analysis that shows that there is no foliation of the 'first quadrant' quarter of a disk $D_\rho$ by $g$-geodesics, which leads to the conclusion that there is no $C^1$ solution $u$ on any $D_\rho$ to the equation $|\nabla u|^2=(xy)^{2n}$.
If there is interest, I can supply details of these arguments when I get the time.
