Abstract Let L be a homogeneous divergence form higher order elliptic operator with complex bounded measurable coefficients and ( p - ( L ) , p + ( L ) ) ${(p_{-}(L),p_{+}(L))}$ be the maximal interval of exponents q ∈ [ 1 , ∞ ] ${q\in[1,\infty]}$ such that the semigroup { e - t L } t > 0 ${\{e^{-tL}\}_{t>0}}$ is bounded on L q ( ℝ n ) ${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 ( ℝ n ) ${H_{L}^{p}(\mathbb{R}^{n})}$ for all p ∈ ( 0 , p + ( L ) ) ${p\in(0,p_{+}(L))}$ , which when p = 1 ${p=1}$ , answers a question asked by Deng, Ding and Yao in [21]. Moreover, the authors characterize H L p ( ℝ n ) ${H_{L}^{p}(\mathbb{R}^{n})}$ via various versions of square functions and Lusin-area functions associated to the operator L.