  The structure of orthogonal polynomials on $\mathbb{R}^{2}$ with the weight
function $| x_{1}^{2}-x_{2}^{2}|^{2k_{0}}|
x_{1}x_{2}|^{2k_{1}}e^{-(x_{1}^{2}+x_{2}^{2})/2}$ is based on the Dunkl
operators of type $B_{2}$. This refers to the full symmetry group of the
square, generated by reflections in the lines $x_{1}=0$ and $x_{1}-x_{2}=0$.
The weight function is integrable if $k_{0},k_{1},k_{0}+k_{1}>-\frac{1}{2}$.
Dunkl operators can be defined for polynomials taking values in a module of the
associated reflection group, that is, a vector space on which the group has an
irreducible representation. The unique 2-dimensional representation of the
group $B_{2}$ is used here. The specific operators for this group and an
analysis of the inner products on the harmonic vector-valued polynomials are
presented in this paper. An orthogonal basis for the harmonic polynomials is
constructed, and is used to define an exponential-type kernel. In contrast to
the ordinary scalar case the inner product structure is positive only when
$(k_{0},k_{1})$ satisfy $-\frac{1}{2}<k_{0}\pm k_{1}<\frac{1}{2}$. For vector
polynomials $(f_{i})_{i=1}^{2}$, $(g_{i})_{i=1}^{2}$ the inner product has the
form $\iint_{\mathbb{R}^{2}}f(x) K(x)
g(x)^{T}e^{-(x_{1}^{2}+x_{2}^{2})/2}dx_{1}dx_{2}$ where the matrix function
$K(x)$ has to satisfy various transformation and boundary conditions. The
matrix $K$ is expressed in terms of hypergeometric functions.
