In your particular case, it seems clear to sufficiency two conditions: the first is that $\nabla u^*(\omega) = 0,$ i.e, $u^*(\omega)$ is parallel.
The last is quite technical so let me introduce some concepts before:
Recall the Weitzenbock formula: $$\Delta = \nabla^*\nabla + Q^p.$$
Since $\Sigma$ is closed we can intrdocute the $L^2$ inner product of $\Omega^p(\Sigma)$ in the following way
$$(\omega,\beta) := \int_{\Sigma}\omega\wedge\star\beta = \int_{\Sigma} g(\omega,\beta)d\mu_g.$$
Therefore, $\Delta$ is self-adjoint on this inner product and $\nabla^*$ is the dual of $\nabla$ on this inner product.
In the particular case of $n=3$, $Q^p$ is nothing but the Ricci tensor on $2$-forms, i.e, if $R$ denotes the Riemann tensor, then $$R(X,Y)\omega := \nabla_X\nabla_Y\omega - \nabla_Y\nabla_X\omega -\nabla_{[X,Y]}\omega,$$ $$\langle R(X,Y)\omega,\omega\rangle $$ the Ricci tensor corresponds to $$(\mathrm{Ric}~\omega)(X_1,X_2) = \sum_{s=1}^3\left(R(X_s,X_1)\omega(X_s,X_1) + R(X_s,X_2)\omega(X_s,X_2\right)),$$ where $\{X_1,X_2,X_3\}$ is an orthonormal pair.
Hence the Weitzenbock formula reduces to (for $p=2$)
$$\Delta = \nabla^*\nabla + \mathrm{Ric}.$$
Therefore, for any parallel $2$-form $\beta$ one has
$$(\Delta\beta,\beta) =\int_{\Sigma} \langle \Delta\beta,\beta\rangle =\int_{\Sigma} \langle \nabla^*\nabla\beta,\beta\rangle + \int_{\Sigma}\langle\mathrm{Ric}~\beta,\beta\rangle = \int_{\Sigma}|\nabla\beta|^2 +\int_{\Sigma}\langle\mathrm{Ric}~\beta,\beta\rangle = \int_{\Sigma}\langle\mathrm{Ric}~\beta,\beta\rangle.$$
Since $$(\Delta\beta,\beta) = (d\delta\beta,\beta) + (\delta d\beta) = |\delta\beta|^2 + |d\beta|^2,$$ onde thus concludes that $u^*(\omega)$ is closed and co-closed, hence, harmonic, provided if $$\int_{\Sigma}\langle\mathrm{Ric}~u^*(\omega),u^*(\omega)\rangle\leq 0,$$ and this is our second condition.
Note that in general, parallel forms are not harmonic, but in this particular that if $g$ is Ricci flat, we are done, since $$\Delta u^*(d\theta_i) = 0,$$ and for $p=1$, $Q^p$ is also the Ricci tensor, it follows that $$0 = \int_{\Sigma}\langle \mathrm{Ric}(u^*(d\theta_i)),u^*(d\theta_i)\rangle +\int_{\Sigma} |\nabla u^*(\theta_i))|^2 = \int_{\Sigma} |\nabla u^*(\theta_i))|^2,$$ and hence, $u^*(\theta_i)$ is parallel, hence, all of our conditions would be satisfied.