Obstructions for the wedge of coordinate differentials to be harmonic Let $(M,g)$ be a smooth $d$-dimensional Riemannian manifold, $d$ even. Are there obstructions (I guess in terms of curvature) for $g$ to have the following property:

For every $p \in M$ there exist a coordinate system around $p$, such that the co-frame associated with it satisfy:
$$ \delta(dx^{i_1} \wedge dx^{i_2} \dots \wedge dx^{i_{\frac{d}{2}}})=0 \, \, \text{for every choice of indices } 1 \le i_1 < i_2 < \dots < i_{\frac{d}{2}} \le d.$$

Here $\delta=d^*$ is the adjoint of the exterior derivative.
In fact, for my purposes it suffices that $dx^{i_1} \wedge dx^{i_2} \dots \wedge dx^{i_{\frac{d}{2}}}$ would be co-closed only for two complementing sets of indices, but I am not sure this problem is easier to analyse in practice. (though any ideas on that would be welcomed).

For the full problem, we have $\binom{d}{\frac{d}{2}} \cdot \binom{d}{\frac{d}{2}-1}$ equations, while the metric has $\frac{d(d+1)}{2}$ degress of freedom, so this problem is probably overdetermined. (Can we prove that a generic metric admits no solutions?).
Given a coordinate system, we can write $dx^i=a^i_je^j$ where $e^j$ is some (positive) orthonormal coframe. Writing $A=a^i_j$, we get that $A^TA=G^{-1}$, where $G=g_{ij}$ is the coordinate representation of the metric. This means we can assume that $A= \sqrt{G^{-1}}$. Then
$$ dx^{i_1} \wedge dx^{i_2} \dots \wedge dx^{i_{k}}=a^{i_1}_{j_1}\dots a^{i_k}_{j_k} e^{j_1} \wedge \dots \wedge e^{j_k}=$$
$$\sum_{1 \le j_1 < j_2 < \dots < j_k \le d} \sum_{\sigma \in S^k} a^{i_1}_{j_{\sigma(1)}}\dots a^{i_k}_{j_{\sigma(k)}} \text{sgn}(\sigma) e^{j_1} \wedge \dots \wedge e^{j_k}= \sum_J A^I_Je^J,$$ where $A^I_J$ is the $k$-minor of the matrix $A= \sqrt{G^{-1}}$ corresponding to columns $I=(i_1,\dots,i_k)$ and rows $J=(j_1,\dots,j_k)$. Here $e^J:=e^{j_1} \wedge \dots \wedge e^{j_k}$. Taking the Hodge dual, we obtain
$$\star dx^{i_1} \wedge dx^{i_2} \dots \wedge dx^{i_{k}}=\sum_J A^I_Je^{J^c},$$
so the final equation is
$$d \star dx^{i_1} \wedge dx^{i_2} \dots \wedge dx^{i_{k}}=\sum_J dA^I_J \wedge e^{J^c}+A^I_J \wedge de^{J^c}.$$
I am not sure how to proceed from here. My idea was to expand $dA^I_J ,de^{J^c}$ and get some first order equation on $G$ and its minors. Then, I guess that second differentiation might give us something which is related to the curvature. However, it doesn't seem easy to do so.
Finally, we might simplify things a bit by assuming $e^{J^c} $ is closed. This raises the question what are the obstructions for such an orthonormal co-frame $e^j$ to exist. Note that if we want the $e^i$ themselves to be closed (not just their wedge product) then this forces the metric to be flat. I am not sure what is the obstruction when $|J|=k>1$.
 A: 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.
A: I don't have an answer to the question about obstructions. However, I can suggest some simplifications to the condition you are trying to satisfy.
Let $J_a^i = \nabla_a x^i$, where $x^i$ are your sought after coordinates. Let us restrict our attention to the case when $J_a^i$ is invertible. In the sequel, the Latin letters $a,b,c,\cdots$ may be freely raised and lowered using the metric $g_{ab}$. The equation you want to satisfy is
$$
  \nabla^a (\wedge^n J)_{a_1\cdots a_n}^{i_i\cdots i_n} = 0 ,
$$
where the dimension is $d=2n$ and the $\wedge$ operation is carried out over the $a_1,\ldots,a_n$ indices. Though as a consequence the above expression is also antisymmetric in the $i_1,\ldots,i_n$ indices. This equation can be simplified by a judicious application of the Leibniz rule:
\begin{align*}
  \nabla^{a_1} (\wedge^n J)_{a_1\cdots a_n}^{i_i\cdots i_n}
  &= n g^{ba_1} (\nabla_b J \wedge (\wedge^{n-1} J))_{a_1\cdots a_n}^{i_1\cdots i_n} \\
  &= n (\nabla^2 x)^{[i_1} (\wedge^{n-1} J)_{a_2\cdots a_n}^{i_2\cdots i_n]}
    + n(n-1) (J^{b[i_2} \nabla_b J^{i_1}) \wedge (\wedge^{n-2} J)^{i_3\cdots i_n]})_{a_2\cdots a_n} \\
  &= n[((\nabla^2 x)^{[i_1} J^{i_2} + (n-1) J^{b[i_2} \nabla_b J^{i_1}) \wedge (\wedge^{n-2} J)^{i_3\cdots i_n]}]_{a_2\cdots a_n}
\end{align*}
I don't have a slick argument for it, but I suspect that the $\wedge^{n-2} J$ factor may be canceled, since we are assuming the non-degeneracy of $J$. Extracting the non-trivial factor, we get the new equation
\begin{align*}
  0
  &= J^{ka}(g_{ab} \nabla^2 x + (n-1)\nabla_b \nabla_a x)^{[i} J^{j]b}
  =: E^{ijk} ,
\end{align*}
where I've added the extra $J^{ka}$ factor for symmetry. In this form, the expression $E^{ijk}$ has the symmetries $E^{[ij]k} = E^{ijk}$ and $E^{[ijk]} = 0$, which corresponds to a Young diagram of the form
$$
  \begin{array}{|c|c|}
    \hline
    i & k \\ \hline
    j \\
    \hline
  \end{array} ,
$$
where, to get the right shape, you should ignore any box not filled with an index. This is an overdetermined system for the unknown coordinate functions $x^i$. Unfortunately, from this point on, I don't see an easy way to extract integrability conditions on the $x^i$, other than taking more covariant derivatives and trying to see what that gives you.
