Let
n
≥
2
n\ge 2
and
Ω
\Omega
be a bounded non-tangentially accessible domain (for short, NTA domain) of
R
n
\mathbb {R}^n
. Assume that
L
D
L_D
is a second-order divergence form elliptic operator having real-valued, bounded, measurable coefficients on
L
2
(
Ω
)
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
\{K_t^{L_D}\}_{t>0}
generated by
L
D
L_D
are Hölder continuous. Second, for any
p
∈
(
0
,
1
]
p\in (0,1]
, the authors introduce the ‘geometrical’ Hardy space
H
r
p
(
Ω
)
H^p_r(\Omega )
by restricting any element of the Hardy space
H
p
(
R
n
)
H^p(\mathbb {R}^n)
to
Ω
\Omega
and the authors show that, when
p
∈
(
n
n
+
δ
0
,
1
]
p\in (\frac {n}{n+\delta _0},1]
,
H
r
p
(
Ω
)
=
H
p
(
Ω
)
=
H
L
D
p
(
Ω
)
H^p_r(\Omega )=H^p(\Omega )=H^p_{L_D}(\Omega )
with equivalent quasi-norms, where
H
p
(
Ω
)
H^p(\Omega )
and
H
L
D
p
(
Ω
)
H^p_{L_D}(\Omega )
respectively denote the Hardy space on
Ω
\Omega
and the Hardy space associated with
L
D
L_D
and where
δ
0
∈
(
0
,
1
]
\delta _0\in (0,1]
is the critical index of the Hölder continuity for the kernels
{
K
t
L
D
}
t
>
0
\{K_t^{L_D}\}_{t>0}
. Third, as applications, the authors obtain the global gradient estimates in both
L
p
(
Ω
)
L^p(\Omega )
, with
p
∈
(
1
,
p
0
)
p\in (1,p_0)
, and
H
z
p
(
Ω
)
H^p_z(\Omega )
, with
p
∈
(
n
n
+
1
,
1
]
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
∈
(
2
,
∞
)
p_0\in (2,\infty )
is a constant depending only on
n
n
,
Ω
\Omega
, and the coefficient matrix of
L
D
L_D
. Here, the ‘geometrical’ Hardy space
H
z
p
(
Ω
)
H^p_z(\Omega )
is defined by restricting any element of the Hardy space
H
p
(
R
n
)
H^p(\mathbb {R}^n)
supported in
Ω
¯
\overline {\Omega }
to
Ω
\Omega
, where
Ω
¯
\overline {\Omega }
denotes the closure of
Ω
\Omega
in
R
n
\mathbb {R}^n
. It is worth pointing out that the range
p
∈
(
1
,
p
0
)
p\in (1,p_0)
for the global gradient estimate in the scale of Lebesgue spaces
L
p
(
Ω
)
L^p(\Omega )
is sharp and the above results are established without any additional assumptions on both the coefficient matrix of
L
D
L_D
and the domain
Ω
\Omega
.