Abstract: Assume that $L$ is a non-negative self-adjoint operator on $L^2(\mathbb{R}^n)$ with its heat kernels satisfying the so-called Gaussian upper bound estimate and that $X$ is a ball quasi-Banach function space on $\mathbb{R}^n$ satisfying some mild assumptions. Let $H_{X,\,L}(\mathbb{R}^n)$ be the Hardy space associated with both $X$ and $L,$ which is defined by the Lusin area function related to the semigroup generated by $L$. In this article, the authors establish various maximal function characterizations of the Hardy space $H_{X,\,L}(\mathbb{R}^n)$ and then apply these characterizations to obtain the solvability of the related Cauchy problem. These results have a wide range of generality and, in particular, the specific spaces $X$ to which these results can be applied include the weighted space, the variable space, the mixed-norm space, the Orlicz space, the Orlicz-slice space, and the Morrey space. Moreover, the obtained maximal function characterizations of the mixed-norm Hardy space, the Orlicz-slice Hardy space, and the Morrey-Hardy space associated with $L$ are completely new.

Key words: Hardy space, ball quasi-Banach function space, Gaussian upper bound estimate, non-negative self-adjoint operator, maximal function