Another strong constraint on the quartic couplings of potential is the requirement for the theory to be unitary. For the SM, this implies an upper bound on the mass of the Higgs boson [1]}, [2]},
and similar studies have been applied to the 2HDM (general unitarity conditions are presented in Ref. [3]}) and other models with a larger scalar content. Essentially, the method consists in computing all scalar-scalar \(J = 0\)  scattering amplitudes (usually denoted \(a_0\) ) and requiring that they respect unitarity in the high energy regime. This translates into an upper bound on those amplitudes, \(|a_0| < 1/2\) . Theories with many scalars complicate the calculation somewhat due to a growing multiplicity of such scattering amplitudes and bounds must be imposed upon eigenvalues of the S-matrix. General unitarity bounds for a 3HDM with a \(\mathbb {Z}_2\times \mathbb {Z}_2\) -symmetric potential are presented in Ref. [4]},
of which our \({\rm U}(1)\times \mathbb {Z}_2\)  symmetry is a special case. Since our model has a larger symmetry, it has less parameters than that of Ref. [4]} and one could read off from Eqs. (91) to (102) of that article the unitarity constraints imposed upon the quartic couplings. However, as explained below, we will use the SARAH/SPheno machinery instead to take into account these bounds.
