  In this paper we establish the characterization of the weighted BMO via two
weight commutators in the settings of the Neumann Laplacian $\Delta_{N_+}$ on
the upper half space $\mathbb{R}^n_+$ and the reflection Neumann Laplacian
$\Delta_N$ on $\mathbb{R}^n$ with respect to the weights associated to
$\Delta_{N_+}$ and $\Delta_{N}$ respectively. This in turn yields a weak
factorization for the corresponding weighted Hardy spaces, where in particular,
the weighted class associated to $\Delta_{N}$ is strictly larger than the
Muckenhoupt weighted class and contains non-doubling weights. In our study, we
also make contributions to the classical Muckenhoupt--Wheeden weighted Hardy
space (BMO space respectively) by showing that it can be characterized via area
function (Carleson measure respectively) involving the semigroup generated by
the Laplacian on $\mathbb{R}^n$ and that the duality of these weighted Hardy
and BMO spaces holds for Muckenhoupt $A^p$ weights with $p\in (1,2]$ while the
previously known related results cover only $p\in (1,{n+1\over n}]$. We also
point out that this two weight commutator theorem might not be true in the
setting of general operators $L$, and in particular we show that it is not true
when $L$ is the Dirichlet Laplacian $\Delta_{D_+}$ on $\mathbb{R}^n_+$.
