AbstractA search for new heavy resonances decaying to a pair of Higgs bosons (HH) in proton-proton collisions at a center-of-mass energy of 13 TeV is presented. Data were collected with the CMS detector at the LHC in 2016–2018, corresponding to an integrated luminosity of 138 fb−1. Resonances with a mass between 0.8 and 4.5 TeV are considered using events in which one Higgs boson decays into a bottom quark pair and the other into final states with either one or two charged leptons. Specifically, the single-lepton decay channel$$ \mathrm{HH}\to \mathrm{b}\overline{\mathrm{b}}{\mathrm{WW}}^{\ast}\to \mathrm{b}\overline{\mathrm{b}}\ell v\mathrm{q}{\overline{\mathrm{q}}}^{\prime } $$HH→bb¯WW∗→bb¯ℓvqq¯′and the dilepton decay channels$$ \mathrm{HH}\to \mathrm{b}\overline{\mathrm{b}}{\mathrm{WW}}^{\ast}\to \mathrm{b}\overline{\mathrm{b}}\ell v\ell v $$HH→bb¯WW∗→bb¯ℓvℓvand$$ \mathrm{HH}\to \mathrm{b}\overline{\mathrm{b}}\uptau \uptau \to \mathrm{b}\overline{\mathrm{b}}\ell vv\ell vv $$HH→bb¯ττ→bb¯ℓvvℓvvare examined, whereℓin the final state corresponds to an electron or muon. The signal is extracted using a two-dimensional maximum likelihood fit of the$$ \mathrm{H}\to \mathrm{b}\overline{\mathrm{b}} $$H→bb¯jet mass and HH invariant mass distributions. No significant excess above the standard model expectation is observed in data. Model-independent exclusion limits are placed on the product of the cross section and branching fraction for narrow spin-0 and spin-2 massive bosons decaying to HH. The results are also interpreted in the context of radion and bulk graviton production in models with a warped extra spatial dimension. The results provide the most stringent limits to date for X→HH signatures with final-state leptons and at some masses provide the most sensitive limits of all X→HH searches.