The model-independent 95\% \CL expected and observed upper limits set on ${\sigma(\PP\to 2\Pa+\PX)\mathcal{B}^2(\Pa\to 2\PGm)\alphaGen}$ over the range $0.21 < \MPa < 9\GeV$ for the combined 2016, 2017, and 2018 analyses. Mass ranges that overlap with \JPsi and \PgU resonances are excluded from the search

The 95\% \CL observed upper limits on the effective coupling $\CahLambda$ of the ALP to the SM Higgs assuming the branching fraction of the ALP to muons is 1 (blue) and 0.1 (orange), for both the expected (dashed) and observed (solid) limits

The 95\% \CL observed upper limits on the effective coupling $\cll/\Lambda$ of the ALP to the SM leptons. The orange shaded region represents the parameter space excluded by this search under three choices of $\CahLambda$ $1\TeV^{-2}$ (solid), $0.1\TeV^{-2}$ (dashed), and $0.01\TeV^{-2}$ (dotted)

The 95\% \CL observed upper limits on $\varepsilon^{2} \mathcal{B}(\ZDark \to \sDark \overline{\sDark}) \mathcal{B}^{2}(\sDark \rightarrow 2\PGm)$ for the vector portal model as a function of the dark scalar mass $\MsDark$ and dark vector boson mass $\MZDark$. The yellow and green bands indicate one and two standard deviation values of the limit, respectively. Because the model-independent limits are calculated only up to a dimuon mass of 60\GeV, the \sDark mass considered for these limits is below 60\GeV

The 95\% \CL observed upper limits for $\sigma(\PP\to\PhOneTwo\to 2\PaOne) \mathcal{B}^2(\PaOne\to 2\PGm)$ for the NMSSM as a function of $\MPaOne$ for $\MPhOneTwo$ = 90~GEV. The data presented here reflect the results of the 2017 and 2018 datasets combined with the previously published results of the 2016 dataset in \cite{CMS:2018jid}

The 95\% \CL observed upper limits for $\sigma(\PP\to\PhOneTwo\to 2\PaOne) \mathcal{B}^2(\PaOne\to 2\PGm)$ for the NMSSM as a function of $\MPaOne$ for $\MPhOneTwo$ = 125~GEV. The data presented here reflect the results of the 2017 and 2018 datasets combined with the previously published results of the 2016 dataset in \cite{CMS:2018jid}

The 95\% \CL observed upper limits for $\sigma(\PP\to\PhOneTwo\to 2\PaOne) \mathcal{B}^2(\PaOne\to 2\PGm)$ for the NMSSM as a function of $\MPaOne$ for $\MPhOneTwo$ = 150~GEV. The data presented here reflect the results of the 2017 and 2018 datasets combined with the previously published results of the 2016 dataset in \cite{CMS:2018jid}

The 95\% \CL observed upper limits (black solid curves) from this search as interpreted in the dark SUSY scenario for the process $\PP \to \Ph \to 2\nOne \to 2\gammaDark + 2\nDark \to 4\PGm + \PX$, with \mbox{$\MnOne = 60\GeV$} and $\MnDark = 1\GeV$. The limits are presented in the plane of $\varepsilon$ and $\MgammaDark$. The color gradient represents different branching fraction assumptions for \mbox{$\mathcal{B}(\Ph \to 2\nOne \to 2\gammaDark + 2\nDark)$}, ranging from dark orange (0.05\%) to light orange (10\%). The degradation of the limit around 1\GeV is attributed to the drop in the dimuon branching fraction $\mathcal{B}(\gammaDark \to 2\mu)$ due to the dimuon resonance of the $\phi$ meson~\cite{Batell:2009yf}

Distribution of the missing transverse energy in CRTop. The Z+HF and $t\bar{t}$ scale factors, described in the text, have been applied to the simulated samples. The distribution labeled "Signal" is for the model with ($m_a$, $m_{\tilde{\chi}_{1}^{0}}$, $m_{\tilde{\chi}_{2}^{0}}$) = (45 GeV, 10 GeV, 80 GeV).

Distribution of the dijet invariant mass in the signal region, shown together with the parameterized background model (labelled "Bkg Model"). For reference, the MC prediction for the SM background is also shown (labelled "SM MC"). The Z+HF and $t\bar{t}$ scale factors, described in the text, have been applied to the simulated samples. The signal region is defined to have dijet invariant mass > 20 GeV. The distribution labeled "Signal" is for the model with ($m_a$, $m_{\tilde{\chi}_{1}^{0}}$, $m_{\tilde{\chi}_{2}^{0}}$) = (45 GeV, 10 GeV, 80 GeV).

Upper limits at 95% CL on the $pp \rightarrow ZH$ cross section times the branching ratio for $Z \rightarrow \ell^{+}\ell^{-}$ (where $\ell = e, \mu \;\mathrm{or}\; \tau$) and $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0} \rightarrow a \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0} \rightarrow b\bar{b} \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ as a function of $m_{a}$ for $m_{\tilde{\chi}_{1}^{0}} = 10$ GeV and $m_{\tilde{\chi}_{2}^{0}} = 65$ GeV in the NMSSM scenario described in the text. All branching ratios in the Higgs boson decay chain after the decay $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0}$ are set to 100%.

Upper limits at 95% CL on the $pp \rightarrow ZH$ cross section times the branching ratio for $Z \rightarrow \ell^{+}\ell^{-}$ (where $\ell = e, \mu \;\mathrm{or}\; \tau$) and $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0} \rightarrow a \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0} \rightarrow b\bar{b} \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ as a function of $m_{a}$ for $m_{\tilde{\chi}_{1}^{0}} = 10$ GeV and $m_{\tilde{\chi}_{2}^{0}} = 80$ GeV in the NMSSM scenario described in the text. All branching ratios in the Higgs boson decay chain after the decay $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0}$ are set to 100%.

Upper limits at 95% CL on the $pp \rightarrow ZH$ cross section times the branching ratio for $Z \rightarrow \ell^{+}\ell^{-}$ (where $\ell = e, \mu\; \mathrm{or}\; \tau$) and $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0} \rightarrow a \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0} \rightarrow b\bar{b} \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ as a function of $m_{a}$ for $m_{\tilde{\chi}_{1}^{0}} = 10$ GeV and $m_{\tilde{\chi}_{2}^{0}} = 95$ GeV in the NMSSM scenario described in the text. All branching ratios in the Higgs boson decay chain after the decay $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0}$ are set to 100%.

Upper limits at 95% CL on the $pp \rightarrow ZH$ cross section times the branching ratio for $Z \rightarrow \ell^{+}\ell^{-}$ (where $\ell = e, \mu\; \mathrm{or}\; \tau$) and $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0} \rightarrow a \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0} \rightarrow b\bar{b} \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ as a function of $m_{a}$ for $m_{\tilde{\chi}_{1}^{0}} = 10$ GeV and $m_{\tilde{\chi}_{2}^{0}} = 110$ GeV in the NMSSM scenario described in the text. All branching ratios in the Higgs boson decay chain after the decay $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0}$ are set to 100%.

Upper limits at 95% CL on the $pp \rightarrow ZH$ cross section times the branching ratio for $Z \rightarrow \ell^{+}\ell^{-}$ (where $\ell = e, \mu\; \mathrm{or}\; \tau$) and $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0} \rightarrow a \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0} \rightarrow b\bar{b} \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ as a function of $m_{a}$ for $m_{\tilde{\chi}_{1}^{0}} = 20$ GeV and $m_{\tilde{\chi}_{2}^{0}} = 80$ GeV in the NMSSM scenario described in the text. All branching ratios in the Higgs boson decay chain after the decay $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0}$ are set to 100%.

Upper limits at 95% CL on the $pp \rightarrow ZH$ cross section times the branching ratio for $Z \rightarrow \ell^{+}\ell^{-}$ (where $\ell = e, \mu \;\mathrm{or}\; \tau$) and $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0} \rightarrow a \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0} \rightarrow b\bar{b} \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ as a function of $m_{a}$ for $m_{\tilde{\chi}_{1}^{0}} = 30$ GeV and $m_{\tilde{\chi}_{2}^{0}} = 80$ GeV in the NMSSM scenario described in the text. All branching ratios in the Higgs boson decay chain after the decay $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0}$ are set to 100%.

Upper limits at 95% CL on the branching ratio $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0} \rightarrow a \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0} \rightarrow b\bar{b} \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ as a function of $m_{a}$ for $m_{\tilde{\chi}_{1}^{0}} = 10$ GeV and $m_{\tilde{\chi}_{2}^{0}} = 65$ GeV in the NMSSM scenario described in the text. The SM $ZH$ cross section is assumed.

Upper limits at 95% CL on the branching ratio $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0} \rightarrow a \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0} \rightarrow b\bar{b} \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ as a function of $m_{a}$ for $m_{\tilde{\chi}_{1}^{0}} = 10$ GeV and $m_{\tilde{\chi}_{2}^{0}} = 80$ GeV in the NMSSM scenario described in the text. The SM $ZH$ cross section is assumed.

Upper limits at 95% CL on the branching ratio $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0} \rightarrow a \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0} \rightarrow b\bar{b} \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ as a function of $m_{a}$ for $m_{\tilde{\chi}_{1}^{0}} = 10$ GeV and $m_{\tilde{\chi}_{2}^{0}} = 95$ GeV in the NMSSM scenario described in the text. The SM $ZH$ cross section is assumed.

Upper limits at 95% CL on the branching ratio $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0} \rightarrow a \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0} \rightarrow b\bar{b} \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ as a function of $m_{a}$ for $m_{\tilde{\chi}_{1}^{0}} = 10$ GeV and $m_{\tilde{\chi}_{2}^{0}} = 110$ GeV in the NMSSM scenario described in the text. The SM $ZH$ cross section is assumed.

Upper limits at 95% CL on the branching ratio $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0} \rightarrow a \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0} \rightarrow b\bar{b} \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ as a function of $m_{a}$ for $m_{\tilde{\chi}_{1}^{0}} = 20$ GeV and $m_{\tilde{\chi}_{2}^{0}} = 80$ GeV in the NMSSM scenario described in the text. The SM $ZH$ cross section is assumed.

Upper limits at 95% CL on the branching ratio $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0} \rightarrow a \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0} \rightarrow b\bar{b} \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ as a function of $m_{a}$ for $m_{\tilde{\chi}_{1}^{0}} = 30$ GeV and $m_{\tilde{\chi}_{2}^{0}} = 80$ GeV in the NMSSM scenario described in the text. The SM $ZH$ cross section is assumed.

Unweighted and weighted number of events after each stage of selection for the NMSSM scenario with $pp \rightarrow ZH$, $Z \rightarrow \ell^{+}\ell^{-}$, $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0} \rightarrow a \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0} \rightarrow b\bar{b} \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ where $(m_{a}, m_{\tilde{\chi}_{1}^{0}}, m_{\tilde{\chi}_{2}^{0}}) = (45,10,80)$ GeV and the Z boson decaying to $e^{+}e^{-}, \;\mu^{+}\mu^{-} \;\mathrm{or}\; \tau^{+}\tau^{-}$. All branching ratios in the Higgs boson decay chain are set to 100%. The weighted number of events corresponds to an integrated luminosity of $139 \;\mathrm{fb}^{-1}$. The "skimming selection" required either a single electron or muon with $p_{T} > 100\;$ GeV or a pair of electrons or muons each with $p_{T} > 20\;$ GeV. The preselection requirements include the trigger, absence of a bad jet or bad muon, and two leptons, where a lepton is either an electron or a muon.

Acceptance and efficiency of this analysis for the signal models considered in this paper. The signal is an NMSSM scenario with $pp \rightarrow ZH$, $Z \rightarrow \ell^{+}\ell^{-}$, $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0} \rightarrow a \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0} \rightarrow b\bar{b} \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ where the values of $(m_{a}, m_{\tilde{\chi}_{1}^{0}}, m_{\tilde{\chi}_{2}^{0}})$ are varied. The Z boson decays to $e^{+}e^{-}, \mu^{+}\mu^{-} \;\mathrm{or}\; \tau^{+}\tau^{-}$. All branching ratios in the Higgs boson decay chain are set to 100%. The product of acceptance times efficiency is defined by the fraction of simulated events that pass all the selection criteria of this analysis. The acceptance is defined by the fraction of simulated events that pass all the selection criteria as applied to Monte Carlo truth-level quantities. The efficiency is then defined as the acceptance times efficiency, divided by the acceptance.

Acceptance and efficiency of this analysis for the signal models considered in this paper. The signal is an NMSSM scenario with $pp \rightarrow ZH$, $Z \rightarrow \ell^{+}\ell^{-}$, $H \rightarrow \tilde{\chi}_{2}^{0}\tilde{\chi}_{1}^{0} \rightarrow a \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0} \rightarrow b\bar{b} \tilde{\chi}_{1}^{0}\tilde{\chi}_{1}^{0}$ where the values of $(m_{a}, m_{\tilde{\chi}_{1}^{0}}, m_{\tilde{\chi}_{2}^{0}})$ are varied. The Z boson decays to $e^{+}e^{-}, \mu^{+}\mu^{-} \;\mathrm{or}\; \tau^{+}\tau^{-}$. All branching ratios in the Higgs boson decay chain are set to 100%. The product of acceptance times efficiency is defined by the fraction of simulated events that pass all the selection criteria of this analysis. The acceptance is defined by the fraction of simulated events that pass all the selection criteria as applied to Monte Carlo truth-level quantities. The efficiency is then defined as the acceptance times efficiency, divided by the acceptance.