In this paper, we investigate the bilinear Stein’s square functions associated with the bilinear Bochner-Riesz means, defined as \( \mathcal {G}^\alpha (f, g)(x):=\left( \int _0^{\infty }\left| \frac{\partial }{\partial R} \mathcal {B}_R^{\alpha +1}(f, g)(x)\right| ^2 R d R\right) ^{\frac{1}{2}}, \) where the bilinear Bochner-Riesz means is given by \(\begin{aligned} \mathcal {B}_R^\alpha (f, g)(x)=\int _{\mathbb {R}^n} \int _{\mathbb {R}^n}\left( 1-\frac{|\xi |^2+|\eta |^2}{R^2}\right) _{+}^\alpha \hat{f}(\xi ) \hat{g}(\eta ) e^{2 \pi i x \cdot (\xi +\eta )} d \xi d \eta , \quad R>0. \end{aligned}\) The weighted strong and weak type estimates for the operator \(\mathcal {G}^\alpha \) and its associated commutators are given. We also show that the commutators \([\vec {b},\mathcal {G}^\alpha ]\) are compact from \(L^{p_1}(w_1)\times L^{p_2}(w_2)\) to \(L^p(v_{\vec {w}})\) for \(1<p_1,p_2<\infty ,\frac{1}{p}=\frac{1}{p_1}+\frac{1}{p_2}\) and \(\alpha >n-\frac{1}{2}\) , where \(\vec {b}=(b_1,b_2)\) , \(b_1,b_2\in {\text{ CMO }}(\mathbb {R}^n)\) , \(\vec {w}=\left( w_1, w_2\right) \in {A}_{\vec {p}}\) and \(v_{\vec {w}}=w_1^{{p}/{p_1}} w_2^{{p}/{ p_2}}.\) Here \(\textrm{CMO}(\mathbb {R}^n)\) is the closure of \({C}_c^\infty (\mathbb {R}^n)\) in the \(\textrm{BMO}(\mathbb {R}^n)\) topology.
