Next-to-leading order unitarity fits in Two-Higgs-Doublet models with soft $\mathbb{Z}_2$ breaking

We fit the next-to-leading order unitarity conditions to the Two-Higgs-Doublet model with a softly broken $\mathbb{Z}_2$ symmetry. In doing so, we alleviate the existing uncertainty on how to treat higher order corrections to quartic couplings of its Higgs potential. A simplified approach to implementing the next-to-leading order unitarity conditions is presented. These new bounds are then combined with all other relevant constraints, including the complete set of LHC Run I data. The upper $95\%$ bounds we find are $4.2$ on the absolute values of the quartic couplings, and $235$ GeV ($100$ GeV) for the mass degeneracies between the heavy Higgs particles in the type I (type II) scenario. In type II, we exclude an unbroken $\mathbb{Z}_2$ symmetry with a probability of $95\%$. All fits are performed using the open-source code HEPfit.