Orthogonal group synchronization aims to recover orthogonal group elements from their noisy pairwise measurements. It has found numerous applications including computer vision, imaging science, and community detection. Due to the orthogonal constraints, it is often challenging to find the least squares estimator in presence of noise. In the recent years, semidefinite relaxation (SDR) and spectral methods have proven to be powerful tools in recovering the group elements. In particular, under additive Gaussian noise, the SDR exactly produces the maximum likelihood estimator (MLE), and both MLE and spectral methods are able to achieve near-optimal statistical error. In this work, we take one step further to quantify the uncertainty of the MLE and spectral estimators by considering their distributions. By leveraging the orthogonality constraints in the likelihood function, we obtain a second-order expansion of the MLE and spectral estimator with the leading terms as an anti-symmetric Gaussian random matrix that is on the tangent space of the orthogonal matrix. This also implies state-of-the-art min-max risk bounds and a confidence region of each group element as a by-product. Our works provide a general theoretical framework that is potentially useful to find an approximate distribution of the estimators arising from many statistical inference problems with manifold constraints. The numerical experiments confirm our theoretical contribution.