Recently, it was established, via lower order moments, that the univariate q-normal distribution, which is the weighting function for q-Hermite polynomials, describes the ensemble-averaged eigenvalue density from many-particle random matrix ensembles generated by k-body interactions (Vyas and Kota 2019 J. Stat. Mech. 103103). These ensembles are generically called embedded ensembles of k-body interactions [EE(k)] and their Gaussian orthogonal ensemble (GOE) and Gaussian unitary ensemble (GUE) versions are called embedded Gaussian orthogonal ensemble (EGOE), EGOE(k), and embedded Gaussian unitary ensemble (EGUE), EGUE(k), respectively. Going beyond the earlier work, this paper examines the lower-order bivariate reduced moments of the distribution of the square of the transition matrix elements (usually called transition strengths) generated by the action of a transition operator O (such as a dipole operator in atoms and molecules, a quadrupole operator in atomic nucleic, etc) on the eigenstates generated by a k-body Hamiltonian H(k). To this end, H is represented by EGOE(k) [or EGUE(t)] and O, that is to say, a t-body, by an independent EGOE(t) [or EGUE(t)]. Given this, it is shown that the ensembleaveraged distribution of transition strengths follows a bivariate q-normal, giving a first application of the bivariate q-normal distribution in the statistical physics of quantum systems. Also presented are the formulas for the bivariate correlation coefficient ρ and the q value, that define a bivariate q-normal as a function of