  In this paper, the authors study the matrix-valued harmonic functions and
characterize them by the Poisson integral of functions in non-commutative BMO
(bounded mean oscillation) spaces. This provides a very satisfactory
non-commutative analogue of the beautiful result due to Fabes, Johnson and Neri
[Indiana Univ. Math. J. {\bf25} (1976) 159-170; MR0394172].
