We provide an analytical description of the galaxy bispectrum covariance and the power spectrum-bispectrum cross-covariance in redshift space that captures the dominant non-Gaussian contributions at mildly non-linear scales, particularly those arising from the squeezed configurations. The Gaussian prediction for the variance of the halo bispectrum monopole significantly underestimates numerical estimates particularly for squeezed triangles, that is bispectrum triangular configurations where one side is much smaller than the other two, whereas the effect is relatively less important when considering the quadrupole. We propose an expression for the missing non-Gaussian contribution valid at mildly non-linear scales in the squeezed limit that requires an accurate modeling of the bispectrum alone. We validate our model against the numerical covariance estimated from a large suite of mock catalogs and find that it accurately predicts the variance as well as the dominant off-diagonal terms. We also present an expression for the cross-covariance between power spectrum and bispectrum multipoles and likewise find it to provide a good description of the numerical results.
