This paper presents the characterization of the covariance matrix function of a Gaussian or second-order elliptically contoured vector random field on the sphere which is stationary, isotropic, and mean square continuous. This characterization involves an infinite sum of the products of positive definite matrices and Gegenbauer's polynomials, and may not be available for other non-Gaussian vector random fields on spheres such as a chi (2) or log-Gaussian vector random field. We also offer two simple but efficient constructing approaches, and derive some parametric covariance matrix structures on spheres.