  In this paper, we will investigate the boundedness of the bi-parameter
Fourier integral operators (or FIOs for short) of the following form:
$$T(f)(x)=\frac{1}{(2\pi)^{2n}}\int_{\mathbb{R}^{2n}}e^{i\varphi(x,\xi,\eta)}\cdot
a(x,\xi,\eta)\cdot\widehat{f}(\xi,\eta)d\xi d\eta,$$ where for $x=(x_1,x_2)\in
\mathbb{R}^{n}\times \mathbb{R}^{n}$ and $\xi,\eta \in
\mathbb{R}^{n}\setminus\{0\}$, the amplitude $a(x,\xi,\eta)\in L^\infty
BS^m_\rho$ and the phase function is of the form $
\varphi(x,\xi,\eta)=\varphi_1(x_1,\xi)+\varphi_2(x_2,\eta)$ with $\quad
\varphi_1,\varphi_2 \in L^\infty \Phi^2
(\mathbb{R}^{n}\times\mathbb{R}^{n}\setminus\{0\})$ and $\varphi(x, \xi, \eta)$
satisfies a certain rough non-degeneracy condition.
  The study of these operators are motivated by the $L^p$ estimates for
one-parameter FIOs and bi-parameter Fourier multipliers and pseudo-differential
operators. We will first define the bi-parameter FIOs and then study the $L^p$
boundedness of such operators when their phase functions have compact support
in frequency variables with certain necessary non-degeneracy conditions. We
will then establish the $L^p$ boundedness of the more general FIOs with
amplitude $a(x,\xi,\eta)\in L^\infty BS^m_\rho$ and non-smooth phase function
$\varphi(x,\xi,\eta)$ on $x$ satisfying a rough non-degeneracy condition.
