  Let $\mathcal{X}$ be a metric space with doubling measure and $L$ a
nonnegative self-adjoint operator in $L^2(\mathcal{X})$ satisfying the
Davies-Gaffney estimates. Let
$\varphi:\,\mathcal{X}\times[0,\infty)\to[0,\infty)$ be a function such that
$\varphi(x,\cdot)$ is an Orlicz function, $\varphi(\cdot,t)\in
A_{\infty}(\mathcal{X})$ (the class of Muckenhoupt weights) and its uniformly
critical lower type index $i(\varphi)\in(0,1]$. In this paper, the authors
introduce a Musielak-Orlicz Hardy space $H_{\varphi,\,L}(\mathcal{X})$ by the
Lusin area function associated with the heat semigroup generated by $L$, and a
Musielak-Orlicz $\mathop\mathrm{BMO}$-type space
$\mathop\mathrm{BMO}_{\varphi,\,L}(\mathcal{X})$ which is further proved to be
the dual space of $H_{\varphi,\,L}(\mathcal{X})$; as a corollary, the authors
obtain the $\varphi$-Carleson measure characterization of
$\mathop\mathrm{BMO}_{\varphi,\,L}(\mathcal{X})$. Characterizations of
$H_{\varphi,\,L}(\mathcal{X})$, including the atom, the molecule and the Lusin
area function associated with the Poisson semigroup of $L$, are presented.
Using the atomic characterization, the authors characterize
$H_{\varphi,\,L}(\mathcal{X})$ in terms of $g^\ast_{\lambda,\,L}$. As further
applications, the authors obtain several equivalent characterizations of the
Musielak-Orlicz Hardy space $H_{\varphi,\,L}(\mathbb{R}^n)$ associated with the
Schr\"odinger operator $L=-\Delta+V$, where $0\le V\in
L^1_{\mathrm{loc}}(\mathbb{R}^n)$ is a nonnegative potential, in terms of the
Lusin-area function, the non-tangential maximal function, the radial maximal
function, the atom and the molecule.
