A search is reported for pairs of light Higgs bosons ($\mathrm{H}_1$) produced in supersymmetric cascade decays in final states with small missing transverse momentum. A data set of LHC pp collisions collected with the CMS detector at $\sqrt s=13\;\mathrm{TeV}$ and corresponding to an integrated luminosity of $137\;\mathrm{fb}^{-1}$ is used. The search targets events where both $\mathrm{H}_1$ bosons decay into $b\bar{b}$ pairs that are reconstructed as large-radius jets using substructure techniques. No evidence is found for an excess of events beyond the background expectations of the standard model (SM). Results from the search are interpreted in the next-to-minimal supersymmetric extension of the SM, where a "singlino" of small mass leads to squark and gluino cascade decays that can predominantly end in a highly Lorentz-boosted singlet-like $\mathrm{H}_1$ and a singlino-like neutralino of small transverse momentum. Upper limits are set on the product of the squark or gluino pair production cross section and the square of the $b\bar{b}$ branching fraction ($\mathcal{B}$) of the $\mathrm{H}_1$ in a benchmark model containing almost mass-degenerate gluinos and light-flavour squarks at the mass $m_{\mathrm{SUSY}}$. Under the assumption of an SM-like $\mathrm{H}_1\to b\bar{b}$ branching fraction, $\mathrm{H}_1$ bosons with masses in the range $40$-$120\;\mathrm{GeV}$ arising from the decays of squarks or gluinos with a mass of $1200$ to $2500\;\mathrm{GeV}$ are excluded at 95% confidence level.
Simplified analysis for reinterpretation
To aid reinterpretation of the search, a simplified analysis is performed using only the 10 search regions in the upper $H_T$ category of the nominal analysis. The results, which are approximately model-independent, are presented here. The value $A_{\mathrm{kin}}$ is defined as the product of acceptance and efficiency for a signal event to satisfy the kinematic selection and the $H_T>3500\;\mathrm{GeV}$ requirement of the nominal analysis. The value of $A_{\mathrm{kin}}$ is common among all 10 search regions in the simplified analysis, and is quoted for the benchmark signal model in Table 3. Upper limits on the product $\sigma \times \mathcal{B}^2 \times A_{\mathrm{kin}}$ as a function of $m_{\mathrm{H_1}}$ are set in Figure 11, from which $\sigma \times \mathcal{B}^2$ limits for different signal models can be derived through division by the appropriate value of $A_{\mathrm{kin}}$. Since the upper $H_T$ category provides most of the sensitivity for $m_{\mathrm{SUSY}}>2000\;\mathrm{GeV}$ in the nominal analysis, the $\sigma\mathcal{B}^2$ upper limits in this region are not much weaker in the simplified analysis. This is not the case in the region $m_{\mathrm{SUSY}}<2000\;\mathrm{GeV}$, where the lower $H_T$ categories become important. The double-b tag and mass region selections are not considered in $A_{\mathrm{kin}}$. This is done for simplicity, and because the fraction of events satisfying these selections is not found to be strongly model-dependent (except for the dependence on $m_{\mathrm{H_1}}$, which is accounted for explicitly in Figure 11). For the benchmark model, this fraction is found to be independent of $m_{\mathrm{SUSY}}$ within 10% in the region $1600<m_{\mathrm{SUSY}}<2800\;\mathrm{GeV}$ and $35<m_{\mathrm{H_1}}<125\;\mathrm{GeV}$. This approximate independence does not hold for $m_{\mathrm{SUSY}}<1600\;\mathrm{GeV}$, where the $\mathrm{H}_1$ $p_T$ distribution has substantial contributions below the $p_T$ necessary for the $\mathrm{H}_1\rightarrow b\bar{b}$ decay products to be merged in a single AK8 jet. Only models with typical $b\bar{b}$ angular separation $\Delta R<0.8$ should be considered for reinterpretation.