The main purpose of this paper is to establish, using the bi-parameter Littlewood–Paley–Stein theory (in particular, the bi-parameter Littlewood–Paley–Stein square functions), a Hörmander–Mihlin type theorem for the following bi-parameter Fourier multipliers on bi-parameter Hardy spaces H^p(\mathbb{R}^{n_1}\!\times \mathbb{R}^{n_2}) ( 0 < p\le 1 ) with optimal smoothness: T_mf(x_1, x_2)=\frac{1}{(2\pi)^{n_1+n_2}}\int_{\mathbb{R}^n\times \mathbb{R}^m}m(\xi, \eta)\hat{f}(\xi, \eta)e^{2\pi(x_1\xi+x_2\eta)}\,d\xi \,d\eta. One of our results (Theorem 1.7) is the following: assume that m(\xi,\eta) is a function on \mathbb{R}^{n_1}\!\times \mathbb{R}^{n_1} satisfying \sup_{j,k\in \mathbb{Z}}\|m_{j,k}\|_{W^{(s_1,s_2)}}<\infty, with s_1 > n_1({1}/{p}-{1}/{2}) , s_2 > n_2({1}/{p}-{1}/{2}) . Then T_m is bounded from H^p(\mathbb{R}^{n_1}\!\times \mathbb{R}^{n_2}) to H^p(\mathbb{R}^{n_1}\!\times \mathbb{R}^{n_2}) for all 0 < p\le 1 , and \|T_m\|_{H^p\rightarrow H^p}\lesssim \sup_{j,k\in \mathbb{Z}}\|m_{j,k}\|_{W^{(s_1,s_2)}}. Moreover, the smoothness assumption on s_1 and s_2 is optimal. Here, m_{j,k}(\xi,\eta)=m(2^j\xi,2^k\eta)\Psi(\xi)\Psi(\eta) , where \Psi(\xi) and \Psi(\eta) are suitable cut-off functions on \mathbb{R}^{n_1} and \mathbb{R}^{n_2} , respectively, and W^{(s_1, s_2)} is a two-parameter Sobolev space on \mathbb{R}^{n_1}\!\times \mathbb{R}^{n_2} . We also establish that under the same smoothness assumption on the multiplier m , \|T_m\|_{H^p\rightarrow L^p}\!\lesssim\! \sup_{j,k\in \mathbb{Z}}\|m_{j,k}\|_{W^{(s_1,s_2)}} and \|T_m\|_{{\rm CMO}_p\rightarrow {\rm {CMO}}_p}\lesssim \sup_{j,k\in \mathbb{Z}}\|m_{j,k}\|_{W^{(s_1,s_2)}} for all 0 < p\le 1 . Moreover, \|T_m\|_{L^p\rightarrow L^p}\lesssim \sup_{j,k\in \mathbb{Z}}\|m_{j,k}\|_{W^{(s_1,s_2)}} for all 1 < p < \infty under the assumption s_1 > n_1/2 and s_2 > n_2/2 .