Neutral scalar pair productions through $W$-boson fusion at multi--TeV muon colliders

Full one-loop electroweak radiative corrections to $\mu^- \mu^+ \to W^\pm W^\mp \to hh$ in the Standard Model are computed for the first time in this work. We then evaluate neutral scalar pair production through vector boson fusion at multi--TeV muon colliders within Two-Higgs-Doublet Model (THDM). In the phenomenological analysis, the enhancement factor, defined as the ratio of the cross sections for SM-like Higgs pair production in the THDM with respect to the corresponding ones in the SM, is examined over the viable regions of the parameter space in the Type-X and Type-Y THDMs. Our findings show that this ratio can reach a factor of $3$ in several regions within the valid parameter space of the Type-X THDM, whereas it ranges from $0.91$ to $0.96$ over the entire parameter space of the Type-Y THDM. Finally, we scan the cross sections for double CP-odd Higgs production over the updated parameter space of the Type-X and Type-Y THDMs. In the Type-Y case at $\sqrt{s} = 10$~TeV with an integrated luminosity of $\mathcal{L} = 10000~\text{fb}^{-1}$, CP-odd Higgs pair production in the $t\bar{t}b\bar{b}$ final state, with subsequent top-quark decays into leptons and bottom quarks taken into account, can be tested with a statistical significance exceeding the $2\sigma$ level at several viable parameter points.