We write, for geometric index values, a probabilistic proof of the product formula for spherical Bessel functions. Our proof has the merit to carry over without any further effort to Bessel-type hypergeometric functions of one matrix argument. Moreover, the representative probability distribution involved in the matrix setting is shown to be closely related to matrix-variate normal distributions and to the symmetrization of upper-left corners of Haar distributed orthogonal matrices. Once we did, we use the latter relation to perform a detailed analysis of this probability distribution. In case it is absolutely continuous with respect to Lebesgue measure on the space of real symmetric matrices, the product formula for Bessel-type hypergeometric functions of two matrix arguments is obtained from Weyl integration formula.
