The baryon acoustic oscillation (BAO) reconstruction plays a crucial role in cosmological analysis for spectroscopic galaxy surveys because it can make the density field effectively more linear and more Gaussian. The combination of the power spectra before and after the BAO reconstruction helps break degeneracies among parameters, then improve the constraints on cosmological parameters. It is therefore important to estimate the covariance matrix between pre- and post-reconstructed power spectra. In this work, we use perturbation theory to estimate the covariance matrix of the related power spectra multipoles, and check the accuracy of the derived covariance model using a large suite of dark matter halo catalogs at $z=0.5$. We find that the diagonal part of the auto covariance is well described by the Gaussian prediction, while the cross covariance deviates from the Gaussian prediction quickly when $k > 0.1\,h\,\mathrm{Mpc}^{-1}$. Additionally, we find the non-Gaussian effect in the non-diagonal part of the cross covariance is comparable to, or even stronger than the pre-reconstruction covariance. By adding the non-Gaussian contribution, we obtain good agreement between analytical and numerical covariance matrices in the non-diagonal part up to $k \simeq 0.15\,h\,\mathrm{Mpc}^{-1}$. The agreement in the diagonal part is also improved, but still under-predicts the correlation in the cross covariance block.