  Products of Siegel upper half spaces are Siegel domains, whose Silov
boundaries have the structure of products $\mathscr H_1\times\mathscr H_2$ of
Heisenberg groups. By the reproducing formula of bi-parameter heat kernel
associated to sub-Laplacians, we show that a function in holomorphic Hardy
space $H^1$ on such a domain has boundary value belonging to bi-parameter Hardy
space $ H^1 (\mathscr H_1\times \mathscr H_2)$. With the help of atomic
decomposition of $ H^1 (\mathscr H_1\times \mathscr H_2)$ and bi-paramete
rharmonic analysis, we show that the Cauchy-Szeg\H o projection is a bounded
operator from $ H^1 (\mathscr H_1\times \mathscr H_2)$ to holomorphic Hardy
space $H^1$, and any holomorphic $H^1$ function can be decomposed as a sum of
holomorphic atoms. Bi-parameter atoms on $\mathscr H_1\times\mathscr H_2$ are
more complicated than $1$-parameter ones, and so are holomorphic atoms.
