  Let $\Omega$ be either $\mathbb{R}^n$ or an unbounded strongly Lipschitz
domain of $\mathbb{R}^n$, and $\Phi$ be a continuous, strictly increasing,
subadditive and positive function on $(0,\infty)$ of upper type 1 and of
strictly critical lower type $p_{\Phi}\in(n/(n+1),1]$. Let $L$ be a divergence
form elliptic operator on $L^2 (\Omega)$ with the Neumann boundary condition
and the heat semigroup generated by $L$ have the Gaussian property
$(G_{\infty})$. In this paper, the authors introduce the Orlicz-Hardy space
$H_{\Phi,\,L}(\Omega)$ via the nontangential maximal function associated with
$\{e^{-t\sqrt{L}}\}_{t\ge0}$, and establish its equivalent characterization in
terms of the Lusin area function associated with $\{e^{-t\sqrt{L}}\}_{t\ge0}$.
The authors also introduce the "geometrical" Orlicz-Hardy space
$H_{\Phi,\,z}(\Omega)$ via the classical Orlicz-Hardy space
$H_{\Phi}(\mathbb{R}^n)$, and prove that the spaces $H_{\Phi,\,L}(\Omega)$ and
$H_{\Phi,\,z}(\Omega)$ coincide with equivalent norms, from which,
characterizations of $H_{\Phi,\,L}(\Omega)$, including the vertical and the
nontangential maximal function characterizations associated with
$\{e^{-tL}\}_{t\ge0}$, and the Lusin area function characterization associated
with $\{e^{-tL}\}_{t\ge0}$, are deduced. All the above results generalize the
well-known results of P. Auscher and E. Russ by taking $\Phi(t)\equiv t$ for
all $t\in(0,\infty)$.
