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$.