  The $k$-Cauchy-Fueter operators, $k=0,1,\ldots$, are quaternionic
counterparts of the Cauchy-Riemann operator in the theory of several complex
variables. The weighted $L^2$ method to solve Cauchy-Riemann equation is
applied to find the canonical solution to the non-homogeneous $k$-Cauchy-Fueter
equation in a weighted $L^2$-space, by establishing the weighted $L^2$
estimate. The weighted $k$-Bergman space is the space of weighted $L^2$
integrable functions annihilated by the $k$-Cauchy-Fueter operator, as the
counterpart of the Fock space of weighted $L^2$-holomorphic functions on
$\mathbb{C}^n$. We introduce the $k$-Bergman orthogonal projection to this
closed subspace, which can be nicely expressed in terms of the canonical
solution operator, and its matrix kernel function. We also find the asymptotic
decay for this matrix kernel function.
