  The main aim of this article is to establish boundedness of singular
integrals with non-smooth kernels on product spaces. Let $L_1$ and $L_2$ be
non-negative self-adjoint operators on $L^2(\mathbb{R}^{n_1})$ and
$L^2(\mathbb{R}^{n_2})$, respectively, whose heat kernels satisfy Gaussian
upper bounds. First, we obtain an atomic decomposition for functions in
$H^1_{L_1,L_2}(\mathbb{R}^{n_1}\times\mathbb{R}^{n_2})$ where the Hardy space
$H^1_{L_1,L_2}(\mathbb{R}^{n_1}\times\mathbb{R}^{n_2})$ associated with $L_1$
and $L_2$ is defined by square function norms, then prove an interpolation
property for this space. Next, we establish sufficient conditions for certain
singular integral operators to be bounded on the Hardy space
$H^1_{L_1,L_2}(\mathbb{R}^{n_1}\times\mathbb{R}^{n_2})$ when the associated
kernels of these singular integrals only satisfy regularity conditions
significantly weaker than those of the standard Calder\'on--Zygmund kernels. As
applications, we obtain endpoint estimates of the double Riesz transforms
associated to Schr\"dingier operators and a Marcinkiewicz-type spectral
multiplier theorem for non-negative self-adjoint operators on product spaces.
