  Let $L$ be a homogeneous divergence form higher order elliptic operator with
complex bounded measurable coefficients and $(p_-(L),\, p_+(L))$ be the maximal
interval of exponents $q\in[1,\,\infty]$ such that the semigroup
$\{e^{-tL}\}_{t>0}$ is bounded on $L^q(\mathbb{R}^n)$. In this article, the
authors establish the non-tangential maximal function characterizations of the
associated Hardy spaces $H_L^p(\mathbb{R}^n)$ for all $p\in(0,\,p_+(L))$,
which, when $p=1$, answers a question asked by Deng et al. in [J. Funct. Anal.
263 (2012), 604-674]. Moreover, the authors characterize $H_L^p(\mathbb{R}^n)$
via various versions of square functions and Lusin-area functions associated to
the operator $L$.
