Consider $\mathbb{R}^3 \times \mathbb{R}^3$ with standard coordinates $(q_1, q_2, q_3, p_1, p_2, p_3)$. For a fixed $v \in \mathbb{R}^3$, consider the function $f : \mathbb{R}^3 \times \mathbb{R}^3 \to \mathbb{R}$ given by $$f(q, p) = \langle v, q \times p \rangle$$ Writing everything out, it's easy to show that $\nabla f = (- v \times p, v \times q)$. Is there an easier way to see this, that doesn't involve writing out the coordinate-wise formulas for cross product and inner product?
[Math] Gradient of cross product
cross productdifferential-geometrymultivariable-calculus
Related Solutions
The core of the issue is whether the definition of the cross-product is invariant. The definition that Wikipedia provides uses coordinates, and is hence manifestly chart dependent. This is mostly fine for ordinary uses of the cross product, where we are mainly sticking to $\mathbb{R}^3$, but the context here is a bit different. You have a vector space, in your case $T_xM$, and you want to define a cross-product on this vector space in a coordinate invariant manner. It is not at all clear that Wikipedia's definition achieves this.
The question then is really this: given a $3$-dimensional inner product space $V$. How do we define a cross product on $V$ that does not depend on a particular choice of isomorphism with $\mathbb{R}^3$? There is actually a rather natural way of defining an invariant cross product, which is to use the Hodge star.
Definition: Let $(V,g)$ be a $3$-dimensional oriented inner product space. Given $\mathbf{u},\mathbf{v}\in V$, we define their cross-product to be $$\mathbf{u} \times \mathbf{v} \equiv \star(\mathbf{u}\wedge \mathbf{v}),$$ where $\star$ is the Hodge star. This is a manifestly coordinate independent definition, which for $V = \mathbb{R}^3$ reduces to the ordinary cross-product.
Let $(\mathbf{v}_1,\mathbf{v}_2,\mathbf{v}_3)$ be a basis for $V$. Then in components, we have $$\mathbf{u} \wedge \mathbf{v} = \frac{1}{2}(u^iv^j-u^jv^i)\,\mathbf{v}_i\wedge \mathbf{v}_j,$$ where we use Einstein summation notation for repeated indices. Avoiding a rather lengthy calculation, you can show that the Hodge star is given in coordinates by $$(\star (\mathbf{u} \wedge \mathbf{v}))_k = \sqrt{g}\,u^iv^j\epsilon_{ijk},$$ where $\sqrt{g}$ is the square root of the determinant of $g$ in whatever coordinate system you're in. These are components of a covector, so if we raise the index (using the inverse metric, not the metric itself, as you did in your OP) we get $$(\mathbf{u} \times \mathbf{v})^\ell = \sqrt{g}\,u^iv^jg^{k\ell}\epsilon_{ijk},$$ which is very similar to the expression from Wikipedia, albeit with the inclusion of the determinant factor.
Ultimately, that determinant is the source of all your troubles. When you write down an expression in coordinate form, you need to make sure that it is tensorial so that they define genuinely coordinate invariant quantities. The easiest way to make sure your expressions are tensorial (besides remaining manifestly coordinate invariant) is to make sure the constituent objects you use are tensors themselves. This is fine for the (inverse) metric and your vectors, but the Levi-Civita symbol is not a tensor.
As the name suggests, the Levi-Civita symbol $\epsilon_{ijk}$ is a symbol whose components are defined to be completely alternating. The Levi-Civita symbol exists independently of any choice of coordinate system and is not designed to be a coordinate invariant object that transforms appropriately (for example, it doesn't really make sense to raise and lower indices of the Levi-Civita symbol). However, a closely related concept is the Levi-Civita tensor (the Levi-Civita symbol is distinct from the Levi-Civita tensor, although people are unfortunately quite sloppy about the distinction), which I will denote with a tilde as $\tilde{\epsilon}$. This is the object defined as $$\tilde{\epsilon}_{ijk} = \sqrt{g}\,\epsilon_{ijk}.$$ This is a bonafide tensor (in fact, the Riemannian volume form, which underlies the fact that cross-products give volume) which transforms appropriately. You can see then that the cross-product is given by $$(\mathbf{u}\times \mathbf{v})^\ell = u^iv^jg^{k\ell}\tilde{\epsilon}_{ijk},$$ and everything works out correctly if the $\epsilon$ in the Wikipedia definition was intended to be the Levi-Civita tensor and not symbol (again, people are unfortunately quite sloppy about this). Now you can see that all components of the formula above are genuine tensors, and so the expression is guaranteed to be well-defined in any coordinate system. Try your calculation again, but this time include the determinant factor (and use the inverse metric), and hopefully you'll find that everything works out.
Best Answer
Denote the triple vector product, or volume form, by $\epsilon(\cdot,\cdot,\cdot)$. Then $$f(q,p):=\langle v,q\times p\rangle=\epsilon(v,q,p)=\left\{\eqalign{ &\langle q,p\times v\rangle\>,\cr &\langle p,v\times q\rangle\>. \cr}\right.$$ It now follows by inspection that $$\nabla_q f=-v\times p\>,\qquad \nabla_p f=v\times q\ .$$