Obstructions for the wedge of coordinate differentials to be harmonic

Update (1 June 2018): I have now figured out the 'little linear algebra lemma' in all dimensions $d = 2n$ and can give a complete answer to the OP's question: A metric $(M^{2n},g)$ possesses coordinate charts of the kind that the OP desires if and only if it is 'locally conformally unimodular Hessian', i.e., every point $p\in M$ has an open neighborhood $V$ on which there exist coordinates $x^1,\ldots,x^{2n}$ and a function $u$ with positive definite $x$-Hessian satisfying the Monge-Ampère equation $\det(\mathrm{Hess}_x(u)) = 1$ so that $g$ is a multiple of the Hessian metric $$ h = \frac{\partial^2u}{\partial x^i\partial x^j}\,\mathrm{d}x^i\mathrm{d}x^j\,. \tag0 $$ When $d = 2n > 2$, the generic metric $g$ is not locally conformally unimodular Hessian, so such coordinate charts do not exist.

Because the proof is a bit clearer in the case $d=4$, I'm going to leave that in and simply indicate how the fundamental lemma (equation (2) below) changes in higher dimensions after the $d=4$ argument. What follows up until then is what I had written before:

I now have a complete answer in the case $d=4$ (the case $d=2$ being trivial). (I suspect that the answer holds for all (even) $d$, but that would require proving a little linear algebra lemma that I don't see an immediate proof of, but maybe it will come to me.)

Here is the answer: When $d=2n>2$ there are definitely obstructions (though, what they are explicitly in terms of curvature, I haven't a clue at this point), as my original argument (retained below) shows.

Meanwhile, in the case $d=4$, one has the following necessary and sufficient condition: $(M^4,g)$ has the desired local coordinate systems if and only if $(M^4,g)$ is locally conformally unimodular Hessian, i.e., every point $p\in M$ has a neighborhood $V$ on which there exist coordinates $x^1,x^2,x^3,x^4$ and a function $u$ with positive definite $x$-Hessian satisfying the Monge-Ampère equation $\det(\mathrm{Hess}_x(u)) = 1$ so that $g$ is a multiple of the Hessian metric $$ h = \frac{\partial^2u}{\partial x^i\partial x^j}\,\mathrm{d}x^i\mathrm{d}x^j\,. \tag1 $$ In fact, for such coordinates, we have $\mathrm{d}\bigl(\ast_g(\mathrm{d}x^i\wedge\mathrm{d}x^j)\bigr) = 0$ for all $1\le i<j\le 4$.

Note: The set of such metrics properly contains the locally conformally flat metrics in dimension $4$.

The argument: Assume that $x^1,x^2,x^3,x^4$ are local coordinates on a simply connected set $V\subset M$ that satisfy the condition $\mathrm{d}\bigl(\ast_g(\mathrm{d}x^i\wedge\mathrm{d}x^j)\bigr) = 0$ for all $1\le i<j\le 4$. Replacing $g$ by a scalar multiple $h$, we can assume that the volume form of $h$ is $\mathrm{d}x^1\wedge\mathrm{d}x^2\wedge\mathrm{d}x^3\wedge\mathrm{d}x^4$. By linear algebra (special to dimension $4$), there exists a basis $\alpha_1,\alpha_2,\alpha_3,\alpha_4$ of $1$-forms on $V$ such that the following identities hold: $$ \begin{aligned} \ast_h(\mathrm{d}x^1\wedge\mathrm{d}x^2)&= \alpha_3\wedge\alpha_4\\ \ast_h(\mathrm{d}x^1\wedge\mathrm{d}x^3)&= \alpha_4\wedge\alpha_2\\ \ast_h(\mathrm{d}x^1\wedge\mathrm{d}x^4)&= \alpha_2\wedge\alpha_3\\ \end{aligned}\qquad \begin{aligned} \ast_h(\mathrm{d}x^3\wedge\mathrm{d}x^4)&= \alpha_1\wedge\alpha_2\\ \ast_h(\mathrm{d}x^4\wedge\mathrm{d}x^2)&= \alpha_1\wedge\alpha_3\\ \ast_h(\mathrm{d}x^2\wedge\mathrm{d}x^3)&= \alpha_1\wedge\alpha_4\\ \end{aligned} \tag2 $$ In fact, the solution of these equations (unique up to replacing each $\alpha_i$ by $-\alpha_i$) is given by $$ \alpha_i = h_{ij}\,\mathrm{d}x^j\qquad\qquad\text{where}\quad h = h_{ij}\,\mathrm{d}x^i\mathrm{d}x^j. $$ Now, the closure assumption is equivalent to assuming that $\mathrm{d}(\alpha_i\wedge\alpha_j) = 0$ for $1\le i<j\le 4$, and, as remarked above, this implies that the $\alpha_i$ themselves are closed. In particular, there exist functions $p_i$ such that $\alpha_i = \mathrm{d}p_i$ where, because of the symmetry of $h_{ij}=h_{ji}$, it follows that $$ \frac{\partial p_i}{\partial x^j} = h_{ij} = \frac{\partial p_j}{\partial x^i}. $$ Hence there exists a function $u$ such that $\mathrm{d}u = p_j\,\mathrm{d}x^j$. Whence we have
$$ h_{ij} = \frac{\partial^2u}{\partial x^i\partial x^j}\,, $$ so that $h$ has the Hessian form above, where, because of the volume form normalization, $u$ satisfies the specified Monge-Ampère equation in the $x$-coordinate system.

Higher dimensions: For general $d=2n$ the version of $(2)$ that has to be proved is the following one: If $\xi^1,\ldots,\xi^{2n}$ is any positively oriented basis of $1$-forms on $V$ such that $h = h_{ij}\xi^i\xi^j$ with $\det(h_{ij})= 1$, then, for any multi-index $I = i_1i_2\cdots i_n$ with $1\le i_1<i_2<\cdots<i_n\le 2n$, if we set $\xi^I = \xi^{i_1}{\wedge}\cdots\wedge\xi^{i_n}$, then we have $$ \ast_h(\xi^I) = \sigma^{IJ}\alpha_{J} $$ where $J=j_1j_2\cdots j_n$ with $1\le j_1<j_2<\cdots<j_n\le 2n$ is the multi-index of length $n$ that is disjoint from $I$, $\sigma^{IJ}=\pm1$ is the sign of the permutation that $IJ$ represents of $12\cdots(2n)$, and $\alpha_J = \alpha_{j_1}{\wedge}\cdots\wedge\alpha_{j_n}$ where $$ \alpha_i = h_{ij}\,\xi^j. $$ (This lemma about the Hodge star operator can be proved using techniques very similar to the case $n=4$.)

Once you have this lemma, taking $\xi^i = \mathrm{d}x^i$, the proof proceeds as in the case $d=4$: We show that $\mathrm{d}\alpha_j = 0$, so that $\alpha_j = \mathrm{d}p_j$ for some function $p_j$, etc.

Original answer (of the OP's side question):

Since the OP asked, I can say something about the simpler case of being able to do this for two complementary pairs of $\tfrac12d$-sets of indices. (This has some bearing on the case of all co-closure for all $\tfrac12d$-sets of indices, since that is asking for much more.)

The basic point is this: For the generic metric $(M^{2n},g)$ in dimension $d=2n$, there will not exist, even locally, $n$ independent functions $x^1,\ldots,x^n$ such that the nonvanishing $n$-form $\alpha = \mathrm{d}x^1\wedge\cdots\wedge\mathrm{d}x^n$ be co-closed.

The reason is this: If $\mathrm{d}(\ast\alpha)=0$, then, since $\ast\alpha$ will be decomposable (and nonzero), there would have to exist, in a neighborhood of any point $p$ in the domain of the $x^i$, $n$ local functions $y^1,\ldots,y^n$ such that $\ast\alpha = \mathrm{d}y^1\wedge\cdots\wedge\mathrm{d}y^n$.

Since $\alpha\wedge(\ast\alpha)\not=0$, it follows that $x^1,\ldots,x^n,y^1,\ldots,y^n$ is a local coordinate system in a neighborhood of $p$. Moreover, the $g$-duality relation shows that the span of $\mathrm{d}x^i$ is $g$-orthogonal to the span of $\mathrm{d}y^i$. Thus, in this local coordinate system, $g$ must take the form $$ g = g_{ij}(x,y)\,\mathrm{d}x^i\mathrm{d}x^j + h_{ij}(x,y)\,\mathrm{d}y^i\mathrm{d}y^j.\tag1 $$ Moreover, setting $G = \sqrt{\det(g_{ij})}$ and $H = \sqrt{\det(h_{ij})}$, we clearly have $$ G\,(\ast\alpha)=\ast(G\alpha) = \ast\bigl(G\,\mathrm{d}x^1\wedge\cdots\wedge\mathrm{d}x^n\bigr) = H\,\mathrm{d}y^1\wedge\cdots\wedge\mathrm{d}y^n = H\,(\ast\alpha), $$ so $G = H$.

Thus, a metric $g$ for which a nontrivial solution to $\mathrm{d}\bigl(\ast(\mathrm{d}x^1\wedge\cdots\wedge\mathrm{d}x^n)\bigr)=0$ exists can always be put locally in the form $(1)$ for some functions $g_{ij}=g_{ji}$ and $h_{ij}=h_{ji}$ with $\sqrt{\det(g_{ij})}=\sqrt{\det(h_{ij})}$. Thus, modulo diffeomorphisms (i.e., changes of coordinates), such metrics depend on at most $n(n{+}1)-1$ functions of $2n$ variables.

However, the general metric in dimension $2n$ depends, modulo diffeomorphisms, on $n(2n{+}1)-2n = n(2n{-}1)$ functions of $2n$ variables, and, by Cartan's count, there is no normal form that can reduce this minimum number of arbitrary functions.

Thus, when $n>1$, the generic metric $(M^{2n},g)$ does not admit any local coordinate system satisfying even the single condition $$ \mathrm{d}\bigl(\ast(\mathrm{d}x^1\wedge\cdots\wedge\mathrm{d}x^n)\bigr)=0. \tag2 $$

In the OP's formulation, one is asking for the conditions that $g$ admit a coordinate system such that $$ \mathrm{d}\bigl(\ast(\mathrm{d}x^{i_1}\wedge\cdots\wedge\mathrm{d}x^{i_n})\bigr) =0.\tag3 $$ for each choice of $1\le i_1<i_2<\cdots<i_n\le 2n$, so this is very highly overdetermined. Update (in view of the added information above): It is now clear that this set of metrics contains the locally conformally unimodular Hessian metrics, which, of course, properly contain the conformally flat metrics.