Assume that \(L_{1}\) and \(L_{2}\) are nonnegative self-adjoint operators whose associated semigroups \(e^{-tL_{1}}\) and \(e^{-tL_{2}}\) satisfy the Davies–Gaffney estimates. In this work, we define the product Hardy spaces \(H^{p}_{L_{1},L_{2}}(\mathbb {R}^{n_{1}}\times \mathbb {R}^{n_{2}})\) and \(H^{p}_{\sqrt{L_{1}},\sqrt{L_{2}}}(\mathbb {R}^{n_{1}}\times \mathbb {R}^{n_{2}})\) associated to operators \(L_{1}\) and \(L_{2}\) via the area integrals \(S_{L_{1},L_{2}}\) and \(S_{\sqrt{L_{1}},\sqrt{L_{2}}}\) for \(0< p\le 1\) , respectively. We prove that these Hardy spaces defined by the square functions can be characterized in terms of the product atoms.
