I have difficulty understanding "trace of map" in the definition of harmonic map.
Let $\phi: (M,g)\to (N,h)$ is map between two Riemannian manifolds, the energy density is defined as $$e(\phi)=\frac12trace_g\phi^*h$$
From my understanding, the trace of $\phi^*h$ is the trace of its matrix representation with basis $\{ \frac{\partial}{\partial x^i}\}$, that is $$trace_g\phi^*h=\sum_{i=1}^m\phi^*h(\frac{\partial}{\partial x^i},\frac{\partial}{\partial x^i})=\sum_{i=1}^mh(\phi_*\frac{\partial}{\partial x^i},\phi_*\frac{\partial}{\partial x^i})=\sum_{i=1}^mh(\frac{\partial\phi^{\alpha}}{\partial x^i}\frac{\partial}{\partial y^{\alpha}},\frac{\partial\phi^{\beta}}{\partial x^i}\frac{\partial}{\partial y^{\beta}})=\sum_{i=1}^m\frac{\partial\phi^{\alpha}}{\partial x^i}\frac{\partial\phi^{\beta}}{\partial x^i}h_{\alpha \beta}$$
One can see the reference here, that it explains that in local coordinates, the energy density can be written as: $$e(\phi)=\frac12g^{ij}h_{\alpha \beta}\frac{\partial\phi^{\alpha}}{\partial x^i} \frac{\partial\phi^{\beta}}{\partial x^j}$$
which is different with my understanding! So the trace certainly involves metric $g$, it seems as the interpretation of $trace_g$ is to raise it up by $g^{ij}$ and regard it as a (1,1)-tensor, instead of a (0,2)-tensor. But why?
Any help would be appreciated!
Abstractly, when you have a linear operator $f: V\to V$, one can take the trace of this linear map as
$$\text{tr} f = \sum_{i=1}^nf_{ii},$$
where we write $f = (f_{ij})$ as a matrix once we fix a basis $\{v_1, \cdots, v_n\}$ on $V$. One can check that $\text{tr}f$ is independent of the choice of basis.
In your case, $G = \phi^* h$ is not a linear operator $V\to V$, but instead a bilinear form $G : V\times V \to \mathbb R$. So one cannot define a trace unless you identify $G$ as $\tilde G :V \to V^*$, where $$\tilde G(v) (\omega) = G(v, \omega)$$ and then use the metric $g$ on $V$ to identify $V^*$ with $V$. That is, by definition,
$$\text{tr}_g (G):= \text{tr} (\tilde G : V\to V^*\overset{\cong}{\to} V).$$
Now let $\{v_1, \cdots, v_n\}$ be a basis on $V$, $\{v^1, \cdots, v^n\}$ its dual basis in $V^*$, and write $G_{ij} = G(v_i, v_j)$. Then $\tilde G : V\to V^*$ in matrix form is given by
$$\tilde G (v_i) =G_{ij} v^j.$$
The identification $V^* \to V$ is given by (you need to check that)
$$ v^j \mapsto g^{ij} v_i.$$
Thus the composition $V \overset{\tilde G}{\to} V^* \to V$ is given by
$$ v_i \mapsto G_{ij} g^{jk} v_k.$$
Lastly, we obtain
$$\text{tr}_g (G) = g^{ij} G_{ij}.$$
Going back to your question. Since $G = \phi^* h$. If we use the coordinate basis on $V = T_xM$, then as you calculated,
$$G_{ij} = \phi^* h\left( \frac{\partial }{\partial x^i}, \frac{\partial }{\partial x^j}\right) = h_{\alpha \beta}\frac{\partial\phi^{\alpha}}{\partial x^i} \frac{\partial\phi^{\beta}}{\partial x^j}$$
and that explains the expression you provided.