Two-neutrino $\beta\beta$ decay of $^{136}$Xe to the first excited $0^+$ state in $^{136}$Ba

We calculate the nuclear matrix element for the two-neutrino $\beta\beta$ decay of $^{136}$Xe into the first excited $0^+$ state of $^{136}$Ba. We use different many-body methods: the quasiparticle random-phase approximation (QRPA) framework, the nuclear shell model, the interacting boson model (IBM-2), and an effective field theory (EFT) for $\beta$ and $\beta\beta$ decays. While the QRPA suggests a decay rate at the edge of current experimental limits, the shell model points to a half-life about two orders of magnitude longer. The predictions of the IBM-2 and the EFT lie in between, and the latter provides systematic uncertainties at leading order. An analysis of the running sum of the nuclear matrix element indicates that subtle cancellations between the contributions of intermediate states can explain the different theoretical predictions. For the EFT, we also present results for two-neutrino $\beta\beta$ decays to the first excited $0^+$ state in other nuclei.