  Let $n\ge 1$ and $\mathfrak{T}_{m}$ be the bilinear square Fourier multiplier
operator associated with a symbol $m$, which is defined by $$
\mathfrak{T}_{m}(f_1,f_2)(x) = \biggl(
\int_{0}^\infty\Big|\int_{(\mathbb{R}^n)^2} e^{2\pi ix\cdot (\xi_1 +\xi_2)
}m(t\xi_1,t\xi_2) \hat{f}_{1}(\xi_1)\hat{f}_{2}(\xi_2)d\xi_1
d\xi_2\Big|^2\frac{dt}{t } \biggr)^{\frac 12}. $$ Let $s$ be an integer with
$s\in[n+1,2n]$ and $p_0$ be a number satisfying $2n/s\le p_0\le 2$. Suppose
that $\nu_{\vec{\omega}}=\prod_{i=1}^2\omega_i^{p/ p_i}$ and each $\omega_i$ is
a nonnegative function on $\mathbb{R}^n$. In this paper, we show that
$\mathfrak{T}_{m}$ is bounded from $L^{p_1}(\omega_1)\times L^{p_2}(\omega_2)$
to $L^p(\nu_{\vec{\omega}})$ if $p_0< p_1, p_2<\infty$ with $1/p=1/p_1+ 1/p_2$.
Moreover, if $p_0>2n/s$ and $p_1=p_0$ or $p_2=p_0$, then $\mathfrak{T}_{m}$ is
bounded from $L^{p_1}(\omega_1)\times L^{p_2}(\omega_2)$ to
$L^{p,\infty}(\nu_{\vec{\omega}})$. The weighted end-point $L\log L$ type
estimate and strong estimate for the commutators of $\mathfrak{T}_{m}$ are also
given. These were done by considering the boundedness of some related
multilinear square functions associated with mild regularity kernels and
essentially improving some basic lemmas which have been used before.
