  The aim of this article is to develop the theory of product Hardy spaces
associated with operators which possess the weak assumption of Davies--Gaffney
heat kernel estimates, in the setting of spaces of homogeneous type. We also
establish a Calder\'on--Zygmund decomposition on product spaces, which is of
independent interest, and use it to study the interpolation of these product
Hardy spaces. We then show that under the assumption of generalized Gaussian
estimates, the product Hardy spaces coincide with the Lebesgue spaces, for an
appropriate range of~$p$.
