The main purpose of this paper is to develop the theory of product Hardy spaces built on Banach lattices on \(\mathbb R^n\times {\mathbb {R}}^m\) . First we introduce new product Hardy spaces \({H}_X({\mathbb {R}}^n\times {\mathbb {R}}^m)\) associated with ball quasi-Banach function spaces \(X({\mathbb {R}}^n\times {\mathbb {R}}^m)\) via applying the Littlewood–Paley–Stein theory. Then we establish a decomposition theorem for \({H}_X({\mathbb {R}}^n\times {\mathbb {R}}^m)\) in terms of the discrete Calderón’s identity. Moreover, we explore some useful and general extrapolation theorems of Rubio de Francia on \(X({\mathbb {R}}^n\times {\mathbb {R}}^m)\) and give some applications to boundedness of operators. Finally, we conclude that the bi-parameter singular integral operators \({{\widetilde{T}}}\) are bounded from \({H}_X({\mathbb {R}}^n\times {\mathbb {R}}^m)\) to itself and bounded from \({H}_X({\mathbb {R}}^n\times {\mathbb {R}}^m)\) to \(X(\mathbb R^n\times {\mathbb {R}}^m)\) via extrapolation. The main results obtained in this paper have a wide range of generality. Especially, we can apply these results to many concrete examples of ball quasi-Banach function spaces, including product Herz spaces, weighted product Morrey spaces and product Orlicz spaces.