$\hat{e}_R$ is the unit radial vector. This is simply $(x,y,z)$ divided by its length, $r$: $\hat{e}_R=\frac{(x,y,z)}{r}$.
$\hat{e}_\theta$ is the unit vector tangent to the sphere, thus perpendicular to $(x,y,z)$, which is also perpendicular to $(0,0,1)$ since changing $\theta$ does not change $z$. Therefore, we simply need to take the cross product of $(x,y,z)$ and $(0,0,1)$ and normalize to get $\hat{e}_\theta=\frac{(-y,x,0)}{\sqrt{r^2-z^2}}$. The sign is chosen so that the vector points counter-clockwise.
$\hat{e}_\phi$ is perpendicular to the other two, so we take the cross product and normalize: $\hat{e}_\phi=\frac{(-xz,-yz,r^2-z^2)}{r\sqrt{r^2-z^2}}$. Again, the sign is chosen to have positive $z$ component.
Thus, we get
$$
\begin{align}
\hat{e}_R&=\frac{(x,y,z)}{r}\\
\hat{e}_\theta&=\frac{(-y,x,0)}{\sqrt{r^2-z^2}}\\
\hat{e}_\phi&=\frac{(-xz,-yz,r^2-z^2)}{r\sqrt{r^2-z^2}}
\end{align}
$$
References
- Source of the polar-coordinate image here from Wikipedia.
Solving (*)
We can solve $(\ast)$ without changing to polar coordinates
$$
\frac{\mathrm{d}y}{\mathrm{d}x}=\frac{2xy}{x^2-y^2}\tag{1}
$$
which can be manipulated to
$$
\frac{\mathrm{d}}{\mathrm{d}y}\frac{x^2}{y}=\frac{2x}{y}\frac{\mathrm{d}x}{\mathrm{d}y}-\frac{x^2}{y^2}=-1\tag{2}
$$
Integrating $(2)$ yields
$$
\frac{x^2}{y}=2c-y\tag{3}
$$
For some $c$. Therefore,
$$
x^2+y^2=2cy\tag{4}
$$
Answer to the Question
A change of coordinate may make the solution more apparent, but the equation should be solvable in either coordinate system.
Best Answer
As you go once around the unit circle ($\theta$ runs from $0$ to $2\pi$), the vector $(\cos2\theta, \sin2\theta)$ makes two full turns counterclockwise, therefore the index is $2$.