  Let $n\ge2$ and $\Omega$ be a bounded non-tangentially accessible domain (for
short, NTA domain) of $\mathbb{R}^n$. Assume that $L_D$ is a second-order
divergence form elliptic operator having real-valued, bounded, measurable
coefficients on $L^2(\Omega)$ with the Dirichlet boundary condition. The main
aim of this article is threefold. First, the authors prove that the heat
kernels $\{K_t^{L_D}\}_{t>0}$ generated by $L_D$ are H\"older continuous.
Second, for any $p\in(0,1]$, the authors introduce the `geometrical' Hardy
space $H^p_r(\Omega)$ by restricting any element of the Hardy space
$H^p(\mathbb{R}^n)$ to $\Omega$, and show that, when
$p\in(\frac{n}{n+\delta_0},1]$, $H^p_r(\Omega)=H^p(\Omega)=H^p_{L_D}(\Omega)$
with equivalent quasi-norms, where $H^p(\Omega)$ and $H^p_{L_D}(\Omega)$
respectively denote the Hardy space on $\Omega$ and the Hardy space associated
with $L_D$, and $\delta_0\in(0,1]$ is the critical index of the H\"older
continuity for the kernels $\{K_t^{L_D}\}_{t>0}$. Third, as applications, the
authors obtain the global gradient estimates in both $L^p(\Omega)$, with
$p\in(1,p_0)$, and $H^p_z(\Omega)$, with $p\in(\frac{n}{n+1},1]$, for the
inhomogeneous Dirichlet problem of second-order divergence form elliptic
equations on bounded NTA domains, where $p_0\in(2,\infty)$ is a constant
depending only on $n$, $\Omega$, and the coefficient matrix of $L_D$. It is
worth pointing out that the range $p\in(1,p_0)$ for the global gradient
estimate in the scale of Lebesgue spaces $L^p(\Omega)$ is sharp and the above
results are established without any additional assumptions on both the
coefficient matrix of $L_D$, and the domain $\Omega$.
