  Let $\pi:\mathcal{X}\to S$ be a holomorphic family of canonically polarized
manifolds over a complex manifold $S$, and $f:\mathcal{X}\to N$ a smooth map
into a Riemannian manifold $N$. Consider the energy function $E: S\to
\mathbb{R}$ that assigns $z\in S$ to the Dirichlet energy of the map
$f|_{\mathcal{X}_z}:\mathcal{X}_z\to N$, where $\mathcal{X}_z=\pi^{-1}(z)$ is
the fiber over $z$. In this article, we compute the second variation formula
for the energy function $E$. As a result, we show that if $N$ is of
non-positive complexified sectional curvature and every map
$f|_{\mathcal{X}_z}:\mathcal{X}_z\to N$ is harmonic, then $E:S\to \mathbb{R}$
is plurisubharmonic.
