  Let $L$ be a homogeneous divergence form higher order elliptic operator with
complex bounded measurable coefficients on $\mathbb{R}^n$ and $X$ a ball
quasi-Banach function space on $\mathbb{R}^n$ satisfying some mild assumptions.
Denote by $H_{X,\, L}(\mathbb{R}^n)$ the Hardy space, associated with both $L$
and $X$, which is defined via the Lusin area function related to the semigroup
generated by $L$. In this article, the authors establish both the maximal
function and the Riesz transform characterizations of $H_{X,\,
L}(\mathbb{R}^n)$. The results obtained in this article have a wide range of
generality and can be applied to the weighted Hardy space, the variable Hardy
space, the mixed-norm Hardy space, the Orlicz--Hardy space, the Orlicz-slice
Hardy space, and the Morrey--Hardy space, associated with $L$. In particular,
even when $L$ is a second order divergence form elliptic operator, both the
maximal function and the Riesz transform characterizations of the mixed-norm
Hardy space, the Orlicz-slice Hardy space, and the Morrey--Hardy space,
associated with $L$, obtained in this article, are totally new.
