Pre-fit MC statistical uncertainties per bin. These values are used to quantify the MC statistical uncertainty contributions to the total uncertainty in $m_t$ based on elements from the covariance matrix obtained in the profile likelihood fit. Each MC statistical contribution is estimated by dividing the relevant covariance matrix entry (see Table 2) by the corresponding pre-fit MC statistical uncertainty for that bin.

The contribution of each source of systematic uncertainty to the total uncertainty in $m_t$ is provided. Each source's contribution is taken directly from the relevant element of the covariance matrix for the profile likelihood fit to $\overline{m_J}$, $m_{jj}$, and $m_{tj}$ using data. The sign of the effect of each uncertainty is included.

The fitted $m_t$ is provided for each replica generated using bootstrapping. Each replica is a 3D histogram of $m_{J}$, $m_{jj}$, and $m_{tj}$, where $m_{J}$ has 50 bins between 145.0 GeV and 205.0 GeV. This method uses the BootstrapGenerator class in ROOT to initialise a pseudo-random number generator for each event, which is seeded by the value of the eventNumber and runNumber variables for a given event in the input dataset. The pseudo-random numbers are drawn from a Poisson distribution with rate parameter equal to 1. These are used to fluctuate the amount by which each replica histogram is filled. There are 1000 replicas, numbered from 0 to 999.

KNO-scaled multiplicity distribution versus the KNO variable $N_{\mathrm{ch}}/\langle N_{\mathrm{ch}}\rangle$ in NSD p-Pb collisions at $\sqrt{s_\mathrm{NN}} = 5020~\mathrm{GeV}$ for the pseudorapidity interval $|\eta|<2.4$.

KNO-scaled multiplicity distribution versus the KNO variable $N_{\mathrm{ch}}/\langle N_{\mathrm{ch}}\rangle$ in NSD p-Pb collisions at $\sqrt{s_\mathrm{NN}} = 5020~\mathrm{GeV}$ for the pseudorapidity interval $|\eta|<3.0$.

KNO-scaled multiplicity distribution versus the KNO variable $N_{\rm ch}/\langle N_{\rm ch}\rangle$ in NSD p-Pb collisions at $\sqrt{s_\mathrm{NN}} = 5020~\mathrm{GeV}$ for the pseudorapidity interval $|\eta|<3.4$.

KNO-scaled multiplicity distribution versus the KNO variable $N_{\mathrm{ch}}/\langle N_{\mathrm{ch}}\rangle$ in NSD p-Pb collisions at $\sqrt{s_\mathrm{NN}} = 5020~\mathrm{GeV}$ for the pseudorapidity interval $-3.4<\eta<5.0$.

KNO-scaled multiplicity distribution versus the KNO variable $N_{\mathrm{ch}}/\langle N_{\mathrm{ch}}\rangle$ in NSD p-Pb collisions at $\sqrt{s_\mathrm{NN}} = 5020~\mathrm{GeV}$ for the pseudorapidity interval $-3.4<\eta<-1.0$.

KNO-scaled multiplicity distribution versus the KNO variable $N_{\mathrm{ch}}/\langle N_{\mathrm{ch}}\rangle$ in NSD p-Pb collisions at $\sqrt{s_\mathrm{NN}} = 5020~\mathrm{GeV}$ for the pseudorapidity interval $2.0<\eta<5.0$.

Average charged-particle multiplicity normalised by the average number of participating nucleon pairs in NSD p-Pb collisions at $\sqrt{s_\mathrm{NN}} = 5020~\mathrm{GeV}$ for the pseudorapidity interval $-3.4<\eta<5.0$.

Suppl. Fig. 7. The reconstructed shower angle for all single photon events. The individual signal and background event type categories added together form the unconstrained prediction.

Suppl. Fig. 7. The constrained covariance matrix for the reconstructed shower angle for all single-photon events. The matrix shows uncertainties and correlations between bins due to flux uncertainties, cross-section uncertainties, hadron reinteraction uncertainties, detector systematic uncertainties, Monte-Carlo statistical uncertainties, and dirt (outside cryostat) uncertainties. Data statistical uncertainties are not included. An example of how to add Pearson data statistical uncertainties can be found in the example code repository.

Suppl. Fig. 7. The unconstrained covariance matrix for reconstructed shower angle for all single-photon events. The matrix shows uncertainties and correlations between bins due to flux uncertainties, cross-section uncertainties, hadron reinteraction uncertainties, detector systematic uncertainties, Monte-Carlo statistical uncertainties, and dirt (outside cryostat) uncertainties. Data statistical uncertainties are not included. An example of how to add Pearson data statistical uncertainties can be found in the example code repository.

Fig. 3. The reconstructed number of protons, with a proton KE energy threshold of 35 MeV, for all single-photon events. The individual signal and background event type categories added together form the unconstrained prediction.

Fig. 4, Suppl. Fig. 8. The reconstructed shower energy for zero proton events binned from 0 to 600 MeV. The individual signal and background event type categories added together form the unconstrained prediction.

Fig. 4. The constrained covariance matrix for the reconstructed shower energy for zero proton events binned from 0 to 600 MeV. The matrix shows uncertainties and correlations between bins due to flux uncertainties, cross-section uncertainties, hadron reinteraction uncertainties, detector systematic uncertainties, Monte-Carlo statistical uncertainties, and dirt (outside cryostat) uncertainties. Data statistical uncertainties are not included. An example of how to add Pearson data statistical uncertainties can be found in the example code repository.

Fig. 4. The unconstrained covariance matrix for reconstructed shower energy for zero proton events binned from 0 to 600 MeV. The matrix shows uncertainties and correlations between bins due to flux uncertainties, cross-section uncertainties, hadron reinteraction uncertainties, detector systematic uncertainties, Monte-Carlo statistical uncertainties, and dirt (outside cryostat) uncertainties. Data statistical uncertainties are not included. An example of how to add Pearson data statistical uncertainties can be found in the example code repository.

Fig. 6. The reconstructed shower energy for zero proton events binned from 150 to 1250 MeV, with LEE prediction. The individual signal and background event type categories added together form the unconstrained prediction. The "LEE" category is the unfolded MiniBooNE LEE model prediction scaled under the neutrino-target nucleon interaction assumption.

Fig. 6. The constrained covariance matrix for the reconstructed shower energy for zero proton events binned from 150 to 1250 MeV. Both with and without the LEE model added. The matrix shows uncertainties and correlations between bins due to flux uncertainties, cross-section uncertainties, hadron reinteraction uncertainties, detector systematic uncertainties, Monte-Carlo statistical uncertainties, and dirt (outside cryostat) uncertainties. Data statistical uncertainties are not included. An example of how to add Pearson data statistical uncertainties can be found in the example code repository.

Fig. 6. The unconstrained covariance matrix for reconstructed shower energy for zero proton events binned from 150 to 1250 MeV. Both with and without the LEE model added. The matrix shows uncertainties and correlations between bins due to flux uncertainties, cross-section uncertainties, hadron reinteraction uncertainties, detector systematic uncertainties, Monte-Carlo statistical uncertainties, and dirt (outside cryostat) uncertainties. Data statistical uncertainties are not included. An example of how to add Pearson data statistical uncertainties can be found in the example code repository.

Suppl. Fig. 10. The reconstructed shower backwards projected distance for zero proton events. The individual signal and background event type categories added together form the unconstrained prediction.

Suppl. Fig. 9. The reconstructed shower angle for zero proton events. Includes LEE model prediction. The individual signal and background event type categories added together form the unconstrained prediction. The "LEE" category is the unfolded MiniBooNE LEE model prediction scaled under the neutrino-target nucleon interaction assumption.

Suppl. Fig. 9. The constrained covariance matrix for the reconstructed shower angle for zero proton events. The matrix shows uncertainties and correlations between bins due to flux uncertainties, cross-section uncertainties, hadron reinteraction uncertainties, detector systematic uncertainties, Monte-Carlo statistical uncertainties, and dirt (outside cryostat) uncertainties. Data statistical uncertainties are not included. An example of how to add Pearson data statistical uncertainties can be found in the example code repository.

Suppl. Fig. 9. The unconstrained covariance matrix for reconstructed shower angle for zero proton events. The matrix shows uncertainties and correlations between bins due to flux uncertainties, cross-section uncertainties, hadron reinteraction uncertainties, detector systematic uncertainties, Monte-Carlo statistical uncertainties, and dirt (outside cryostat) uncertainties. Data statistical uncertainties are not included. An example of how to add Pearson data statistical uncertainties can be found in the example code repository.

1$\gamma$Xp shower energy and angle selection efficiencies for true 1$\gamma$Xp events. NaN values indicate a "divide by zero" case where zero events were present in the bin before any selection was applied. Energy and Angle values represent bin lower edges.

1$\gamma$0p shower energy and angle selection efficiencies for true single-photon events with no other true particles of any energy. NaN values indicate a "divide by zero" case where zero events were present in the bin before any selection was applied. Energy and Angle values represent bin lower edges.

Covariance matrix showing uncertainties and correlations between bins due to flux uncertainties, cross-section uncertainties, hadron reinteraction uncertainties, detector systematic uncertainties, Monte-Carlo statistical uncertainties, and dirt (outside cryostat) uncertainties. Pearson data statistical uncertainties have been included, and include small correlations due to events which can be selected by both WC and Pandora. The first four bins are the WC $1\gamma Np$, WC $1\gamma 0p$, Pandora $1\gamma 1p$, and Pandora $1\gamma 0p$ channels, and the next four sets of 16 bins are the NC $\pi^0 Np$, NC $\pi^0 0p$, $\nu_\mu$CC $Np$, and $\nu_\mu$CC $0p$ channels. This corresponds to Fig. 1 of the paper and Fig. 6 of the Supplemental Material.

Constrained WC $1\gamma Np$ energy distribution. There are 6 bins from 0 to 600 MeV of reconstructed shower energy, with an additional overflow bin. Constrained systematic uncertainties on the predictions are illustrated, and a more detailed covariance matrix is included in the Constrained Energy Distributions Covariance Matrix tab. This corresponds to Fig. 2 of the paper.

Constrained WC $1\gamma 0p$ energy distribution. There are 14 bins from 0 to 700 MeV of reconstructed shower energy, with an additional overflow bin. Constrained systematic uncertainties on the predictions are illustrated, and a more detailed covariance matrix is included in the Constrained Energy Distributions Covariance Matrix tab. This corresponds to Fig. 2 of the paper.

Constrained Pandora $1\gamma 1p$ energy distribution. There are 6 bins from 0 to 600 MeV of reconstructed shower energy, with an additional overflow bin. Constrained systematic uncertainties on the predictions are illustrated, and a more detailed covariance matrix is included in the Constrained Energy Distributions Covariance Matrix tab. This corresponds to Fig. 2 of the paper.

Constrained Pandora $1\gamma 0p$ energy distribution. There are 14 bins from 0 to 700 MeV of reconstructed shower energy, with an additional overflow bin. Constrained systematic uncertainties on the predictions are illustrated, and a more detailed covariance matrix is included in the Constrained Energy Distributions Covariance Matrix tab. This corresponds to Fig. 2 of the paper.

Constrained covariance matrix showing constrained uncertainties and correlations between bins due to flux uncertainties, cross-section uncertainties, hadron reinteraction uncertainties, detector systematic uncertainties, Monte-Carlo statistical uncertainties, and dirt (outside cryostat) uncertainties. Pearson data statistical uncertainties have been included, and include small correlations due to events which can be selected by both WC and Pandora. The first seven bins are the WC $1\gamma Np$ reconstructed shower energy from 0-600 MeV with an overflow bin, the next 15 bins are the WC $1\gamma 0p$ reconstructed shower energy from 0-700 MeV with an overflow bin, the next seven bins are the Pandora $1\gamma 1p$ reconstructed shower energy from 0-600 MeV with an overflow bin, and the last 15 bins are the Pandora $1\gamma 0p$ reconstructed shower energy from 0-700 MeV with an overflow bin. This corresponds to Fig. 2 of the paper.

2D efficiency matrix showing the efficiency of the WC $1\gamma Np$ selection with respect to all true NC Delta Radiative Decay events which interact in the fiducial volume. The efficiency is shown as a function of true photon energy and angle. There are 10 bins from -1 to 1 in true photon $\cos \theta$ and 15 bins from 0 to 1500 MeV in true photon energy. NaN values indicate a "divide by zero" case where zero events were present in the bin before any selection was applied. This corresponds to Fig. 9 of the Supplemental Material.

2D efficiency matrix showing the efficiency of the WC $1\gamma 0p$ selection with respect to all true NC Delta Radiative Decay events which interact in the fiducial volume. The efficiency is shown as a function of true photon energy and angle. There are 10 bins from -1 to 1 in true photon $\cos \theta$ and 30 bins from 0 to 1500 MeV in true photon energy. NaN values indicate a "divide by zero" case where zero events were present in the bin before any selection was applied. This corresponds to Fig. 9 of the Supplemental Material.

2D efficiency matrix showing the efficiency of the Pandora $1\gamma 1p$ selection with respect to all true NC Delta Radiative Decay events which interact in the fiducial volume. The efficiency is shown as a function of true photon energy and angle. There are 10 bins from -1 to 1 in true photon $\cos \theta$ and 15 bins from 0 to 1500 MeV in true photon energy, with lower bin edges indicated in the table. NaN values indicate a "divide by zero" case where zero events were present in the bin before any selection was applied. This corresponds to Fig. 9 of the Supplemental Material.

2D efficiency matrix showing the efficiency of the Pandora $1\gamma 0p$ selection with respect to all true NC Delta Radiative Decay events which interact in the fiducial volume. The efficiency is shown as a function of true photon energy and angle. There are 10 bins from -1 to 1 in true photon $\cos \theta$ and 30 bins from 0 to 1500 MeV in true photon energy, with lower bin edges indicated in the table. NaN values indicate a "divide by zero" case where zero events were present in the bin before any selection was applied. This corresponds to Fig. 9 of the Supplemental Material.

2D efficiency matrix showing the efficiency of the WC $1\gamma 0p$ selection with respect to all true zero-primary-final-state-proton NC Delta Radiative Decay events which interact in the fiducial volume. The efficiency is shown as a function of true photon energy and angle. There are 10 bins from -1 to 1 in true photon $\cos \theta$ and 30 bins from 0 to 1500 MeV in true photon energy, with lower bin edges indicated in the table. NaN values indicate a "divide by zero" case where zero events were present in the bin before any selection was applied.

2D efficiency matrix showing the efficiency of the Pandora $1\gamma 0p$ selection with respect to all true zero-primary-final-state-proton NC Delta Radiative Decay events which interact in the fiducial volume. The efficiency is shown as a function of true photon energy and angle. There are 10 bins from -1 to 1 in true photon $\cos \theta$ and 30 bins from 0 to 1500 MeV in true photon energy, with lower bin edges indicated in the table. NaN values indicate a "divide by zero" case where zero events were present in the bin before any selection was applied.

Fig. 3 bottom figure - Distributions of MC simulation compared with data for reconstructed shower energy in the 1$e$0$p$0$\pi$ signal channel, along with the LEE Signal Model 2. The signal and background event categories are summed to form the unconstrained prediction (excluding LEE). Signal events correspond to $\nu_e$ CC events. Background events include $\nu$ with $\pi^0$ events, $\nu$ other events, and cosmic ray events. In Fig. 3, the LEE component is plotted on top of the constrained prediction (excluding LEE) for illustrative purposes. In all statistical tests (results summarized in Table I), the prediction under an LEE hypothesis corresponds to a constrained prediction including LEE. The statistical uncertainties of data use a combined Neyman-Pearson (CNP) version (Eq.(19) in https://doi.org/10.1016/j.nima.2020.163677).

Fig. 4 top figure - Distributions of MC simulation compared with data for reconstructed shower cos($\theta$) in the 1$e$N$p$0$\pi$ signal channel, along with the LEE Signal Model 2. The signal and background event categories are summed to form the unconstrained prediction (excluding LEE). Signal events correspond to $\nu_e$ CC events. Background events include $\nu$ with $\pi^0$ events, $\nu$ other events, and cosmic ray events. In Fig. 4, the LEE component is plotted on top of the constrained prediction (excluding LEE) for illustrative purposes. In all statistical tests (results summarized in Table I), the prediction under an LEE hypothesis corresponds to a constrained prediction including LEE. The statistical uncertainties of data use a combined Neyman-Pearson (CNP) version (Eq.(19) in https://doi.org/10.1016/j.nima.2020.163677).

Fig. 4 bottom figure - Distributions of MC simulation compared with data for reconstructed shower cos($\theta$) in the 1$e$0$p$0$\pi$ signal channel, along with the LEE Signal Model 2. The signal and background event categories are summed to form the unconstrained prediction (excluding LEE). Signal events correspond to $\nu_e$ CC events. Background events include $\nu$ with $\pi^0$ events, $\nu$ other events, and cosmic ray events. In Fig. 4, the LEE component is plotted on top of the constrained prediction (excluding LEE) for illustrative purposes. In all statistical tests (results summarized in Table I), the prediction under an LEE hypothesis corresponds to a constrained prediction including LEE. The statistical uncertainties of data use a combined Neyman-Pearson (CNP) version (Eq.(19) in https://doi.org/10.1016/j.nima.2020.163677).

Fig. 1 top figure - Reconstructed neutrino energy distributions of the 1$\mu$N$p$0$\pi$ control sample (sideband). The signal and background event categories are summed to form the unconstrained prediction. Signal events correspond to $\nu_\mu$ CC events. Background events include $\nu_e$ CC events, NC without $\pi^0$ events, NC with $\pi^0$ events, and cosmic ray events. The statistical uncertainties of data use a combined Neyman-Pearson (CNP) version (Eq.(19) in https://doi.org/10.1016/j.nima.2020.163677).

Fig. 1 middle figure - Reconstructed neutrino energy distributions of the 1$\mu$0$p$0$\pi$ control sample (sideband). The signal and background event categories are summed to form the unconstrained prediction. Signal events correspond to $\nu_\mu$ CC events. Background events include $\nu_e$ CC events, NC without $\pi^0$ events, NC with $\pi^0$ events, and cosmic ray events. The statistical uncertainties of data use a combined Neyman-Pearson (CNP) version (Eq.(19) in https://doi.org/10.1016/j.nima.2020.163677).

Fig. 1 bottom figure - Reconstructed neutrino energy distributions of the NC$\pi^{0}$ control sample (sideband). The signal and background event categories are summed to form the unconstrained prediction. Signal events correspond to NC with $\pi^0$ events. Background events include $\nu_\mu$ CC events, $\nu_e$ CC events, NC without $\pi^0$ events, and cosmic ray events. Only bins up to 1 GeV in this sideband are released and are used in the constraint procedure due to low statistics above this energy. The statistical uncertainties of data use a combined Neyman-Pearson (CNP) version (Eq.(19) in https://doi.org/10.1016/j.nima.2020.163677).

The covariance matrix includes both signal channels and sideband channels, with the LEE component either excluded or included in the signal channels. All channels are binned in reconstructed neutrino energy. Fig. 2 in Supplemental Material shows the corresponding correlation matrix. The constraint procedure uses the correlations with sideband channels to update the prediction and its uncertainty in the signal region. The matrix accounts for uncertainties and correlations between bins due to flux, cross-section, hadron reinteraction, detector systematics, MC statistics, and cosmic ray statistics. Data statistical uncertainties are not included.

The covariance matrix includes both signal channels and sideband channels, with the LEE component either excluded or included in the signal channels. Signal channels are binned shower energy and sideband channels are binned in reconstructed neutrino energy. The constraint procedure uses the correlations with sideband channels to update the prediction and its uncertainty in the signal region. The matrix accounts for uncertainties and correlations between bins due to flux, cross-section, hadron reinteraction, detector systematics, MC statistics, and cosmic ray statistics. Data statistical uncertainties are not included.

The covariance matrix includes both signal channels and sideband channels, with the LEE component either excluded or included in the signal channels. Signal channels are binned shower cos($\theta$) and sideband channels are binned in reconstructed neutrino energy. The constraint procedure uses the correlations with sideband channels to update the prediction and its uncertainty in the signal region. The matrix accounts for uncertainties and correlations between bins due to flux, cross-section, hadron reinteraction, detector systematics, MC statistics, and cosmic ray statistics. Data statistical uncertainties are not included.

The covariance matrices of the sideband channels, binned in reconstructed neutrino energy. The systematic covariance matrix accounts for uncertainties and correlations between bins due to flux, cross-section, hadron reinteraction, detector systematics, MC statistics, and cosmic ray statistics. The total covariance matrix additionally includes data statistical uncertainties (CNP version). The $\chi^2$ values shown in Fig. 1 are calculated using this total covariance matrix, together with the data and MC predictions from the released sideband distribtions.

The unconstrained covariance matrices of the signal channels, binned in reconstructed neutrino energy, with the LEE component either excluded or included. The systematic covariance matrix accounts for uncertainties and correlations between bins due to flux, cross-section, hadron reinteraction, detector systematics, MC statistics, and cosmic ray statistics. The total covariance matrix additionally includes data statistical uncertainties (CNP version). The unconstrained systematic uncertainty in neutrino energy, shown in Supplemental Material Fig. 19, is derived from this systematic covariance matrix (excluding LEE).

The constrained covariance matrices of the signal channels, binned in reconstructed neutrino energy, with the LEE component either excluded or included. The constrained systematic uncertainty in neutrino energy, shown in Supplemental Material Fig. 19, is derived from the systematic covariance matrix (excluding LEE). The $\chi^2$ values shown in Fig. 2 and Table I are calculated using these total covariance matrices, together with the data and constrained predictions from the released signal distributions under the $H_0$ hypothesis (LEE excluded) and $H_1$ hypothesis (LEE included).

The unconstrained covariance matrices of the signal channels, binned in reconstructed shower energy, with the LEE component either excluded or included. The systematic covariance matrix accounts for uncertainties and correlations between bins due to flux, cross-section, hadron reinteraction, detector systematics, MC statistics, and cosmic ray statistics. The total covariance matrix additionally includes data statistical uncertainties (CNP version). The unconstrained systematic uncertainty in shower energy, shown in Supplemental Material Fig. 20, is derived from this systematic covariance matrix (excluding LEE).

The constrained covariance matrices of the signal channels, binned in reconstructed shower energy, with the LEE component either excluded or included. The constrained systematic uncertainty in shower energy, shown in Supplemental Material Fig. 20, is derived from the systematic covariance matrix (excluding LEE). The $\chi^2$ values shown in Fig. 3 and Table I are calculated using these total covariance matrices, together with the data and constrained predictions from the released signal distributions under the $H_0$ hypothesis (LEE excluded) and $H_1$ hypothesis (LEE included).

The unconstrained covariance matrices of the signal channels, binned in reconstructed shower cos($\theta$), with the LEE component either excluded or included. The systematic covariance matrix accounts for uncertainties and correlations between bins due to flux, cross-section, hadron reinteraction, detector systematics, MC statistics, and cosmic ray statistics. The total covariance matrix additionally includes data statistical uncertainties (CNP version). The unconstrained systematic uncertainty in shower cos($\theta$), shown in Supplemental Material Fig. 21, is derived from this systematic covariance matrix (excluding LEE).

The constrained covariance matrices of the signal channels, binned in reconstructed shower cos($\theta$), with the LEE component either excluded or included. The constrained systematic uncertainty in shower cos($\theta$), shown in Supplemental Material Fig. 21, is derived from the systematic covariance matrix (excluding LEE). The $\chi^2$ values shown in Fig. 4 and Table I are calculated using these total covariance matrices, together with the data and constrained predictions from the released signal distributions under the $H_0$ hypothesis (LEE excluded) and $H_1$ hypothesis (LEE included).

Normalized transverse momentum correlator, $\sqrt{ \langle\langle \Delta p_{{\rm T}1}\Delta p_{{\rm T}2} \rangle\rangle }/\langle\langle p_{\rm T} \rangle\rangle $, as a function of the average number of participant nucleons, $\langle{ N_{\rm part}}\rangle$, in Pb--Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV

Normalized transverse momentum correlator, $\sqrt{ \langle\langle \Delta p_{{\rm T}1}\Delta p_{{\rm T}2} \rangle\rangle }/\langle\langle p_{\rm T} \rangle\rangle $, as a function of the average number of participant nucleons, $\langle{ N_{\rm part}}\rangle$, in Xe--Xe collisions at $\sqrt{s_{\rm NN}}$ = 5.44 TeV

Normalized transverse momentum correlator, $\sqrt{ \langle\langle \Delta p_{{\rm T}1}\Delta p_{{\rm T}2} \rangle\rangle }/\langle\langle p_{\rm T} \rangle\rangle $, as a function of the average number of participant nucleons, $\langle{ N_{\rm part}}\rangle$, in Pb--Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV

Normalized transverse momentum correlator, $\sqrt{ \langle\langle \Delta p_{{\rm T}1}\Delta p_{{\rm T}2} \rangle\rangle }/\langle\langle p_{\rm T} \rangle\rangle $ as a function of the charged-particle multiplicity density for spherocity-integrated events in pp collisions at $\sqrt{s}=5.02$ TeV

Normalized transverse momentum correlator, $\sqrt{ \langle\langle \Delta p_{{\rm T}1}\Delta p_{{\rm T}2} \rangle\rangle }/\langle\langle p_{\rm T} \rangle\rangle $ as a function of the charged-particle multiplicity density for jetty events in pp collisions at $\sqrt{s}=5.02$ TeV

Normalized transverse momentum correlator, $\sqrt{ \langle\langle \Delta p_{{\rm T}1}\Delta p_{{\rm T}2} \rangle\rangle }/\langle\langle p_{\rm T} \rangle\rangle $ as a function of the charged-particle multiplicity density for isotropic events in pp collisions at $\sqrt{s}=5.02$ TeV

Normalized transverse momentum correlator, $\sqrt{ \langle\langle \Delta p_{{\rm T}1}\Delta p_{{\rm T}2} \rangle\rangle }/\langle\langle p_{\rm T} \rangle\rangle $ as a function of the charged-particle multiplicity density for spherocity-integrated events in pp collisions at $\sqrt{s}=13$ TeV

Normalized transverse momentum correlator, $\sqrt{ \langle\langle \Delta p_{{\rm T}1}\Delta p_{{\rm T}2} \rangle\rangle }/\langle\langle p_{\rm T} \rangle\rangle $ as a function of the charged-particle multiplicity density for jetty events in pp collisions at $\sqrt{s}=13$ TeV

Normalized transverse momentum correlator, $\sqrt{ \langle\langle \Delta p_{{\rm T}1}\Delta p_{{\rm T}2} \rangle\rangle }/\langle\langle p_{\rm T} \rangle\rangle $ as a function of the charged-particle multiplicity density for isotropic events in pp collisions at $\sqrt{s}=13$ TeV

Ratio of the $\eta/\pi^{0}$ ratio to the $m_{\rm T}$-scaling prediction as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV.

$x_{\rm T}$ scaled $\pi^{0}$ spectrum for pp collisions at $\sqrt{s}$ = 0.9 TeV, with $n = 5.01$.

$x_{\rm T}$ scaled $\pi^{0}$ spectrum for pp collisions at $\sqrt{s}$ = 2.76 TeV, with $n = 5.01$.

$x_{\rm T}$ scaled $\pi^{0}$ spectrum for pp collisions at $\sqrt{s}$ = 7 TeV, with $n = 5.01$.

$x_{\rm T}$ scaled $\pi^{0}$ spectrum for pp collisions at $\sqrt{s}$ = 8 TeV, with $n = 5.01$.

$x_{\rm T}$ scaled $\pi^{0}$ spectrum for pp collisions at $\sqrt{s}$ = 13 TeV, with $n = 5.01$.

$x_{\rm T}$ scaled $\eta$ spectrum for pp collisions at $\sqrt{s}$ = 2.76 TeV, with $n = 4.91$.

$x_{\rm T}$ scaled $\eta$ spectrum for pp collisions at $\sqrt{s}$ = 7 TeV, with $n = 4.91$.

$x_{\rm T}$ scaled $\eta$ spectrum for pp collisions at $\sqrt{s}$ = 8 TeV, with $n = 4.91$.

$x_{\rm T}$ scaled $\eta$ spectrum for pp collisions at $\sqrt{s}$ = 13 TeV, with $n = 4.91$.

Invariant differential yields of the $\pi^{0}$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 0-100%. (Scaled for visibility in the figure)

Invariant differential yields of the $\pi^{0}$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 70-100%. (Scaled for visibility in the figure)

Invariant differential yields of the $\pi^{0}$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 50-70%. (Scaled for visibility in the figure)

Invariant differential yields of the $\pi^{0}$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 30-50%. (Scaled for visibility in the figure)

Invariant differential yields of the $\pi^{0}$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 20-30%. (Scaled for visibility in the figure)

Invariant differential yields of the $\pi^{0}$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 10-20%. (Scaled for visibility in the figure)

Invariant differential yields of the $\pi^{0}$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 5-10%. (Scaled for visibility in the figure)

Invariant differential yields of the $\pi^{0}$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 1-5%. (Scaled for visibility in the figure)

Invariant differential yields of the $\pi^{0}$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 0-1%. (Scaled for visibility in the figure)

Invariant differential yields of the $\pi^{0}$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 0.05-0.1%. (Scaled for visibility in the figure)

Invariant differential yields of the $\pi^{0}$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 0-0.1%. (Scaled for visibility in the figure)

Invariant differential yields of the $\pi^{0}$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 0.01-0.05%. (Scaled for visibility in the figure)

Invariant differential yields of the $\pi^{0}$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 0-0.01%. (Scaled for visibility in the figure)

Invariant differential yields of the $\eta$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 0-100%. (Scaled for visibility in the figure)

Invariant differential yields of the $\eta$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 70-100%. (Scaled for visibility in the figure)

Invariant differential yields of the $\eta$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 50-70%. (Scaled for visibility in the figure)

Invariant differential yields of the $\eta$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 30-50%. (Scaled for visibility in the figure)

Invariant differential yields of the $\eta$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 20-30%. (Scaled for visibility in the figure)

Invariant differential yields of the $\eta$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 10-20%. (Scaled for visibility in the figure)

Invariant differential yields of the $\eta$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 5-10%. (Scaled for visibility in the figure)

Invariant differential yields of the $\eta$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 1-5%. (Scaled for visibility in the figure)

Invariant differential yields of the $\eta$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 0-1%. (Scaled for visibility in the figure)

Invariant differential yields of the $\eta$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 0.05-0.1%. (Scaled for visibility in the figure)

Invariant differential yields of the $\eta$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 0-0.1%. (Scaled for visibility in the figure)

Invariant differential yields of the $\eta$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 0.01-0.05%. (Scaled for visibility in the figure)

Invariant differential yields of the $\eta$ meson versus transverse momentum for pp collisions at $\sqrt{s}$ = 13 TeV in the INEL>0 event class in the V0M multiplicity class 0-0.01%. (Scaled for visibility in the figure)

Ratios of the invariant differential yields of $\pi^{0}$ meson in the V0M multiplicity class 70-100% to the inclusive spectrum in the INEL>0 event class.

Ratios of the invariant differential yields of $\pi^{0}$ meson in the V0M multiplicity class 50-70% to the inclusive spectrum in the INEL>0 event class.

Ratios of the invariant differential yields of $\pi^{0}$ meson in the V0M multiplicity class 30-50% to the inclusive spectrum in the INEL>0 event class.

Ratios of the invariant differential yields of $\pi^{0}$ meson in the V0M multiplicity class 20-30% to the inclusive spectrum in the INEL>0 event class (data points not shown in figure 11).

Ratios of the invariant differential yields of $\pi^{0}$ meson in the V0M multiplicity class 10-20% to the inclusive spectrum in the INEL>0 event class.

Ratios of the invariant differential yields of $\pi^{0}$ meson in the V0M multiplicity class 5-10% to the inclusive spectrum in the INEL>0 event class (data points not shown in figure 11).

Ratios of the invariant differential yields of $\pi^{0}$ meson in the V0M multiplicity class 1-5% to the inclusive spectrum in the INEL>0 event class.

Ratios of the invariant differential yields of $\pi^{0}$ meson in the V0M multiplicity class 0-1% to the inclusive spectrum in the INEL>0 event class (data points not shown in figure 11).

Ratios of the invariant differential yields of $\pi^{0}$ meson in the V0M multiplicity class 0.05-0.1% to the inclusive spectrum in the INEL>0 event class (data points not shown in figure 11).

Ratios of the invariant differential yields of $\pi^{0}$ meson in the V0M multiplicity class 0-0.1% to the inclusive spectrum in the INEL>0 event class (data points not shown in figure 11).

Ratios of the invariant differential yields of $\pi^{0}$ meson in the V0M multiplicity class 0.01-0.05% to the inclusive spectrum in the INEL>0 event class (data points not shown in figure 11).

Ratios of the invariant differential yields of $\pi^{0}$ meson in the V0M multiplicity class 0-0.01% to the inclusive spectrum in the INEL>0 event class.

Ratios of the invariant differential yields of $\eta$ meson in the V0M multiplicity class 70-100% to the inclusive spectrum in the INEL>0 event class.

Ratios of the invariant differential yields of $\eta$ meson in the V0M multiplicity class 50-70% to the inclusive spectrum in the INEL>0 event class.

Ratios of the invariant differential yields of $\eta$ meson in the V0M multiplicity class 30-50% to the inclusive spectrum in the INEL>0 event class.

Ratios of the invariant differential yields of $\eta$ meson in the V0M multiplicity class 20-30% to the inclusive spectrum in the INEL>0 event class (data points not shown in figure 11).

Ratios of the invariant differential yields of $\eta$ meson in the V0M multiplicity class 10-20% to the inclusive spectrum in the INEL>0 event class.

Ratios of the invariant differential yields of $\eta$ meson in the V0M multiplicity class 5-10% to the inclusive spectrum in the INEL>0 event class (data points not shown in figure 11).

Ratios of the invariant differential yields of $\eta$ meson in the V0M multiplicity class 1-5% to the inclusive spectrum in the INEL>0 event class.

Ratios of the invariant differential yields of $\eta$ meson in the V0M multiplicity class 0-1% to the inclusive spectrum in the INEL>0 event class (data points not shown in figure 11).

Ratios of the invariant differential yields of $\eta$ meson in the V0M multiplicity class 0.05-0.1% to the inclusive spectrum in the INEL>0 event class (data points not shown in figure 11).

Ratios of the invariant differential yields of $\eta$ meson in the V0M multiplicity class 0-0.1% to the inclusive spectrum in the INEL>0 event class (data points not shown in figure 11).

Ratios of the invariant differential yields of $\eta$ meson in the V0M multiplicity class 0.01-0.05% to the inclusive spectrum in the INEL>0 event class (data points not shown in figure 11).

Ratios of the invariant differential yields of $\eta$ meson in the V0M multiplicity class 0-0.01% to the inclusive spectrum in the INEL>0 event class.

Slope parameter (n) of a power-law fit ($f(p_{\rm T}) = \alpha p_{\rm T}^{-n}$) for $5 < p_{\rm T} < 15$ GeV/$c$ as a function the charged-particle multiplicity density in units of the average normalized multiplicity density for the INEL>0 event class for the neutral pion.

Slope parameter (k) of an exponential fit ($f(p_{\rm T}) = \alpha e^{-p_{\rm T}/k}$) for $0.4 < p_{\rm T} < 2$ GeV/$c$ as a function the charged-particle multiplicity density in units of the average normalized multiplicity density for the INEL>0 event class for the neutral pion.

Slope parameter (n) of a power-law fit ($f(p_{\rm T}) = \alpha p_{\rm T}^{-n}$) for $5 < p_{\rm T} < 15$ GeV/$c$ as a function the charged-particle multiplicity density in units of the average normalized multiplicity density for the INEL>0 event class for the eta meson.

Slope parameter (k) of an exponential fit ($f(p_{\rm T}) = \alpha e^{-p_{\rm T}/k}$) for $1.5 < p_{\rm T} < 4$ GeV/$c$ as a function the charged-particle multiplicity density in units of the average normalized multiplicity density for the INEL>0 event class for the eta meson.

$\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV in INEL>0 collisions in the V0M multiplicity class 0-100%.

$\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV in INEL>0 collisions in the V0M multiplicity class 70-100% (data points not shown in figure 13).

$\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV in INEL>0 collisions in the V0M multiplicity class 50-70%.

$\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV in INEL>0 collisions in the V0M multiplicity class 30-50% (data points not shown in figure 13).

$\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV in INEL>0 collisions in the V0M multiplicity class 20-30% (data points not shown in figure 13).

$\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV in INEL>0 collisions in the V0M multiplicity class 10-20% (data points not shown in figure 13).

$\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV in INEL>0 collisions in the V0M multiplicity class 5-10% (data points not shown in figure 13).

$\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV in INEL>0 collisions in the V0M multiplicity class 1-5% (data points not shown in figure 13).

$\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV in INEL>0 collisions in the V0M multiplicity class 0-1% (data points not shown in figure 13).

$\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV in INEL>0 collisions in the V0M multiplicity class 0.05-0.1% (data points not shown in figure 13).

$\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV in INEL>0 collisions in the V0M multiplicity class 0-0.1%.

$\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV in INEL>0 collisions in the V0M multiplicity class 0.01-0.05% (data points not shown in figure 13).

$\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV in INEL>0 collisions in the V0M multiplicity class 0-0.01% (data points not shown in figure 13).

Ratio of $\eta/\pi^{0}$ ratio in INEL>0 collisions in the V0M multiplicity class 70-100% to the inclusive $\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV (data points not shown in figure 13).

Ratio of $\eta/\pi^{0}$ ratio in INEL>0 collisions in the V0M multiplicity class 50-70% to the inclusive $\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV.

Ratio of $\eta/\pi^{0}$ ratio in INEL>0 collisions in the V0M multiplicity class 30-50% to the inclusive $\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV (data points not shown in figure 13).

Ratio of $\eta/\pi^{0}$ ratio in INEL>0 collisions in the V0M multiplicity class 20-30% to the inclusive $\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV (data points not shown in figure 13).

Ratio of $\eta/\pi^{0}$ ratio in INEL>0 collisions in the V0M multiplicity class 10-20% to the inclusive $\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV (data points not shown in figure 13).

Ratio of $\eta/\pi^{0}$ ratio in INEL>0 collisions in the V0M multiplicity class 5-10% to the inclusive $\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV (data points not shown in figure 13).

Ratio of $\eta/\pi^{0}$ ratio in INEL>0 collisions in the V0M multiplicity class 1-5% to the inclusive $\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV (data points not shown in figure 13).

Ratio of $\eta/\pi^{0}$ ratio in INEL>0 collisions in the V0M multiplicity class 0-1% to the inclusive $\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV (data points not shown in figure 13).

Ratio of $\eta/\pi^{0}$ ratio in INEL>0 collisions in the V0M multiplicity class 0.05-0.1% to the inclusive $\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV (data points not shown in figure 13).

Ratio of $\eta/\pi^{0}$ ratio in INEL>0 collisions in the V0M multiplicity class 0-0.1% to the inclusive $\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV.

Ratio of $\eta/\pi^{0}$ ratio in INEL>0 collisions in the V0M multiplicity class 0.01-0.05% to the inclusive $\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV (data points not shown in figure 13).

Ratio of $\eta/\pi^{0}$ ratio in INEL>0 collisions in the V0M multiplicity class 0-0.01% to the inclusive $\eta/\pi^{0}$ ratio as a function of $p_{\rm T}$ for pp collisions at $\sqrt{s}$ = 13 TeV (data points not shown in figure 13).

Mean values of the ratio of the $\eta/\pi^{0}$ ratio in a certain multiplicity interval to the $\eta/\pi^{0}$ ratio in the INEL>0 event class for $1 < p_{\rm{ T}} < 4 \rm{~GeV}/c$ as a function of the normalized mean charged-particle multiplicity.

Mean values of the ratio of the $\eta/\pi^{0}$ ratio in a certain multiplicity interval to the $\eta/\pi^{0}$ ratio in the INEL>0 event class for $5 < p_{\rm{ T}} < 15 \rm{~GeV}/c$ as a function of the normalized mean charged-particle multiplicity.

Correlation functions $P_2^{\rm LS}$ of charged hadrons measured in minimum bias pp collisions at $\sqrt{s}=13\;\text{TeV}$. Charged hadrons are selected in the range $0.2 < p_{\rm T} < 2.0$ GeV/$c$ and with pseudorapidity $|\eta| < 0.8$.

Correlation functions $R_2^{\rm CI}$ of charged hadrons measured in minimum bias pp collisions at $\sqrt{s}=13\;\text{TeV}$. Charged hadrons are selected in the range $0.2 < p_{\rm T} < 2.0$ GeV/$c$ and with pseudorapidity $|\eta| < 0.8$.

Correlation functions $P_2^{\rm CI}$ of charged hadrons measured in minimum bias pp collisions at $\sqrt{s}=13\;\text{TeV}$. Charged hadrons are selected in the range $0.2 < p_{\rm T} < 2.0$ GeV/$c$ and with pseudorapidity $|\eta| < 0.8$.

Projections onto $\Delta\eta$ of $R_2^{\rm CI}$ correlation functions of charged hadrons calculated in the selected $p_{\rm T}$ range in pp collisions at $\sqrt{s}=13\;\text{TeV}$. The $\Delta\eta$ projection is calculated as averages of the two-dimensional correlations in the range $|\Delta\varphi| \leq \pi$ and $|\Delta\eta| \leq 1.6$, respectively. Simulations using PYTHIA8 with the Monash 2013 tune and color reconnection (CR) enabled are added.

Projections onto $\Delta\varphi$ of $R_2^{\rm CI}$ correlation functions of charged hadrons calculated in the selected $p_{\rm T}$ range in pp collisions at $\sqrt{s}=13\;\text{TeV}$. The $\Delta\varphi$ projection is calculated as averages of the two-dimensional correlations in the range $|\Delta\varphi| \leq \pi$ and $|\Delta\eta| \leq 1.6$, respectively. Simulations using PYTHIA8 with the Monash 2013 tune and color reconnection (CR) enabled are added.

Projections onto $\Delta\eta$ of $P_2^{\rm CI}$ correlation functions of charged hadrons calculated in the selected $p_{\rm T}$ range in pp collisions at $\sqrt{s}=13\;\text{TeV}$. The $\Delta\eta$ projection is calculated as averages of the two-dimensional correlations in the range $|\Delta\varphi| \leq \pi$ and $|\Delta\eta| \leq 1.6$, respectively. Simulations using PYTHIA8 with the Monash 2013 tune and color reconnection (CR) enabled are added.

Projections onto $\Delta\varphi$ of $P_2^{\rm CI}$ correlation functions of charged hadrons calculated in the selected $p_{\rm T}$ range in pp collisions at $\sqrt{s}=13\;\text{TeV}$. The $\Delta\varphi$ projection is calculated as averages of the two-dimensional correlations in the range $|\Delta\varphi| \leq \pi$ and $|\Delta\eta| \leq 1.6$, respectively. Simulations using PYTHIA8 with the Monash 2013 tune and color reconnection (CR) enabled are added.

Correlation functions $R_2^{\rm CD}$ of charged hadrons measured in minimum bias pp collisions at $\sqrt{s}=13\;\text{TeV}$. Charged hadrons are selected in the range $0.2 < p_{\rm T} < 2.0$ GeV/$c$ and with pseudorapidity $|\eta| < 0.8$.

Correlation functions $P_2^{\rm CD}$ of charged hadrons measured in minimum bias pp collisions at $\sqrt{s}=13\;\text{TeV}$. Charged hadrons are selected in the range $0.2 < p_{\rm T} < 2.0$ GeV/$c$ and with pseudorapidity $|\eta| < 0.8$.

Projections onto $\Delta\eta$ of $R_2^{\rm CD}$ correlation functions of charged hadrons calculated in the selected $p_{\rm T}$ range in pp collisions at $\sqrt{s}=13\;\text{TeV}$. The $\Delta\eta$ projection is calculated as averages of the two-dimensional correlations in the range $|\Delta\varphi| \leq \pi$ and $|\Delta\eta| \leq 1.6$, respectively. Simulations using PYTHIA8 with the Monash 2013 tune and color reconnection (CR) enabled are added.

Projections onto $\Delta\varphi$ of $R_2^{\rm CD}$ correlation functions of charged hadrons calculated in the selected $p_{\rm T}$ range in pp collisions at $\sqrt{s}=13\;\text{TeV}$. The $\Delta\varphi$ projection is calculated as averages of the two-dimensional correlations in the range $|\Delta\varphi| \leq \pi$ and $|\Delta\eta| \leq 1.6$, respectively. Simulations using PYTHIA8 with the Monash 2013 tune and color reconnection (CR) enabled are added.

Projections onto $\Delta\eta$ of $P_2^{\rm CD}$ correlation functions of charged hadrons calculated in the selected $p_{\rm T}$ range in pp collisions at $\sqrt{s}=13\;\text{TeV}$. The $\Delta\eta$ projection is calculated as averages of the two-dimensional correlations in the range $|\Delta\varphi| \leq \pi$ and $|\Delta\eta| \leq 1.6$, respectively. Simulations using PYTHIA8 with the Monash 2013 tune and color reconnection (CR) enabled are added.

Projections onto $\Delta\varphi$ of $P_2^{\rm CD}$ correlation functions of charged hadrons calculated in the selected $p_{\rm T}$ range in pp collisions at $\sqrt{s}=13\;\text{TeV}$. The $\Delta\varphi$ projection is calculated as averages of the two-dimensional correlations in the range $|\Delta\varphi| \leq \pi$ and $|\Delta\eta| \leq 1.6$, respectively. Simulations using PYTHIA8 with the Monash 2013 tune and color reconnection (CR) enabled are added.

Width of $R_2^{\rm CI}(\Delta\eta)$ and $P_2^{\rm CI}(\Delta\eta)$ correlation functions along $\Delta\eta$ measured within $|\Delta\varphi| < \pi$ in pp collisions and compared with results from $p-Pb$ and $Pb-Pb$ collisions as a function of $\langle \mathrm{d}N_\mathrm{ch}/\mathrm{d}\eta \rangle_{|\eta|<0.5}$.

Width of $R_2^{\rm CD}(\Delta\eta)$ and $P_2^{\rm CD}(\Delta\eta)$ correlation functions along $\Delta\eta$ measured within $|\Delta\varphi| < \pi$ in pp collisions and compared with results from p$-$Pb and Pb$-$Pb collisions as a function of $\langle \mathrm{d}N_\mathrm{ch}/\mathrm{d}\eta \rangle_{|\eta|<0.5}$.

Width of $R_2^{\rm CD}(\Delta\varphi)$ and $P_2^{\rm CD}(\Delta\varphi)$ correlation functions along $\Delta\varphi$ measured within $|\Delta\eta| < 1.6$ in pp collisions and compared with results from p$-$Pb and Pb$-$Pb collisions as a function of $\langle \mathrm{d}N_\mathrm{ch}/\mathrm{d}\eta \rangle_{|\eta|<0.5}$.

Balance function of charged hadrons, in the selected $p_{\rm T}$ range, obtained using ALICE data in pp collisions at $\sqrt{s}=13\;\text{TeV}$.

Integral of $B$ results for charged hadrons in pp collisions at $\sqrt{s}=13\;\text{TeV}$ compared with the integral of $B$ for identified hadron ($\pi \pi$, $\pi$K, and $\pi$p) pairs in Pb$-$Pb collisions at $\sqrt{s_{\rm NN}}=2.76\;\text{TeV}$.

Total invariant mass distribution $m_{4\ell}$ in Signal Region 2.

Selected events in the ($m_{4\ell}$, $\left\langle m_{\ell\ell}\right\rangle$) plane, SR1

Selected events in the ($m_{4\ell}$, $\left\langle m_{\ell\ell}\right\rangle$) plane, SR2

Expected background events in the $(m_{4\ell}, \langle m_{\ell\ell}\rangle)$ plane, SR1

Expected signal shape for $(m_{S}, m_{Z_d}) = (110~\text{GeV},30~\text{GeV}$) in the $(m_{4\ell}, \langle m_{\ell\ell}\rangle)$ plane, SR1. The signal distributions are normalized to unit volume.

Expected background events in the $(m_{4\ell}, \langle m_{\ell\ell}\rangle)$ plane, SR2. Events that satisfy $|m_{ab,cd} - m_Z| < 8~\text{GeV}$ are removed by the Z-veto requirement, however not all events lying behind the mask in the plot, corresponding to $\langle m_{\ell\ell}\rangle$ values between 83 and 99 GeV, are deleted.

Expected signal shape for $(m_{S}, m_{Z_{d}}) = (750~\text{GeV},150~\text{GeV}$) in the $(m_{4\ell}, \langle m_{\ell\ell}\rangle)$ plane, SR2. The signal distributions are normalized to unit volume. Events that satisfy $|m_{ab,cd} - m_Z| < 8~\text{GeV}$ are removed by the Z-veto requirement, however not all events lying behind the mask in the plot, corresponding to $\langle m_{\ell\ell}\rangle$ values between 83 and 99 GeV, are deleted.

Expected 95% CL limit on $\sigma \mathcal{B}(S\rightarrow Z_{d}Z_{d} \rightarrow 4\ell)$, SR1

Observed 95% CL limit on $\sigma \mathcal{B}(S\rightarrow Z_{d}Z_{d} \rightarrow 4\ell)$, SR1

Observed / Expected 95% CL limit ratio in SR1

Expected 95% CL limit on $\mathcal{B}(S\rightarrow Z_{d}Z_{d} \rightarrow 4\ell)$ with Z-veto, SR2

Observed 95% CL limit on $\sigma \mathcal{B}(S\rightarrow Z_{d}Z_{d} \rightarrow 4\ell)$ with Z-veto, SR2

Observed / Expected 95% CL limit ratio with Z-veto, SR2

Local $p_0$-values in the ($m_{S}, m_{Z_{d}}$) plane, SR1

Local $p_0$-values in the ($m_{S}, m_{Z_{d}}$) plane, SR2 with Z-veto applied

Observed data in Signal Region 1 in the $(m_{4\ell}, \langle m_{\ell\ell}\rangle)$ plane.

Observed data in Signal Region 2 in the $(m_{4\ell}, \langle m_{\ell\ell}\rangle)$ plane. Events that satisfy $|m_{ab,cd} - m_Z| < 8~\text{GeV}$ are removed by the Z-veto requirement, however not all events lying behind the mask in the plot, corresponding to <m_ll> values between 83 and 99 GeV, are deleted.

Distributions of the average dilepton mass $\left\langle m_{\ell\ell}\right\rangle = \frac{1}{2}\left(m_{ab} + m_{cd}\right)$ for Signal Region 1 in the $4e$ final state.

Distributions of the average dilepton mass $\left\langle m_{\ell\ell}\right\rangle = \frac{1}{2}\left(m_{ab} + m_{cd}\right)$ for Signal Region 1 in the $2e2\mu$ final state.

Distributions of the average dilepton mass $\left\langle m_{\ell\ell}\right\rangle = \frac{1}{2}\left(m_{ab} + m_{cd}\right)$ for Signal Region 1 in the $4\mu$ final state.

Distributions of the average dilepton mass $\left\langle m_{\ell\ell}\right\rangle = \frac{1}{2}\left(m_{ab} + m_{cd}\right)$ for Signal Region 2 in the $4e$ final state.

Distributions of the average dilepton mass $\left\langle m_{\ell\ell}\right\rangle = \frac{1}{2}\left(m_{ab} + m_{cd}\right)$ for Signal Region 2 in the $2e2\mu$ final state.

Distributions of the average dilepton mass $\left\langle m_{\ell\ell}\right\rangle = \frac{1}{2}\left(m_{ab} + m_{cd}\right)$ for Signal Region 2 in the $4\mu$ final state.

Distributions of the total invariant mass $m_{4\ell}$ for Signal Region 1 in the $4e$ final state. The SR1 requirement $m_{4\ell} < 115~\text{GeV}$ is not applied here.

Distributions of the total invariant mass $m_{4\ell}$ for Signal Region 1 in the $2e2\mu$ final state. The SR1 requirement $m_{4\ell} < 115~\text{GeV}$ is not applied here.

Distributions of the total invariant mass $m_{4\ell}$ for Signal Region 1 in the $4\mu$ final state. The SR1 requirement $m_{4\ell} < 115~\text{GeV}$ is not applied here.

Distributions of the total invariant mass $m_{4\ell}$ for Signal Region 2 in the $4e$ final state.

Distributions of the total invariant mass $m_{4\ell}$ for Signal Region 2 in the $2e2\mu$ final state.

Distributions of the total invariant mass $m_{4\ell}$ for Signal Region 2 in the $4\mu$ final state.