The derivation of zonal polynomials involves evaluating the integral \[ \exp\left( - \frac{1}{2} \operatorname{tr} D_β Q D_{l} Q \right) \] with respect to orthogonal matrices \(Q\), where \(D_β\) and \(D_{l}\) are diagonal matrices. The integral is expressed through a polynomial expansion in terms of the traces of these matrices, leading to the identification of zonal polynomials as symmetric, homogeneous functions of the variables \(l_1, l_2, \ldots, l_n\). The coefficients of these polynomials are derived systematically from the structure of the integrals, revealing relationships between them and illustrating the significance of symmetry in their formulation. Furthermore, properties such as the uniqueness up to normalization are established, reinforcing the foundational role of zonal polynomials in statistical and mathematical applications involving orthogonal matrices.