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.

Self-Normalized $J/\psi$ yields in the same arms before and after dimuon subtraction and the comparison with PYTHIA 8 simulations

Self-normalized $\psi(2S)$ to $J/\psi$ yields ratio as a function of Self-normalized $N_{ch}$ in the same arms

Self-normalized $\psi(2S)$ to $J/\psi$ yields ratio as a function of Self-normalized $N_{ch}$ in opposite arms

$\xi$ distributions for different jet $p_T$ bins.

$\xi$ distributions for different jet $p_T$ bins.

The $j_T/p_{T}^{\rm jet}$ distributions for different jet $p_T$ bins.

The $j_T/p_{T}^{\rm jet}$ distributions for different jet $p_T$ bins.

$r$ distributions for different jet $p_T$ bins.

The 95% CL limits on $|V_{Ne}|^2$ as a function of the HNL mass for a Dirac HNL. Values of $-1$ indicate that no limit is available for the mass point.

The 95% CL limits on $|V_{N\mu}|^2$ as a function of the HNL mass for a Dirac HNL. Values of $-1$ indicate that no limit is available for the mass point.

The 95% CL limits on mixed coupling as a function of the HNL mass for a Dirac HNL. Values of $-1$ indicate that no limit is available for the mass point.

The transverse-mass dependence of the correlation-strength parameter $\lambda$ in 30-40% centrality bin obtained from Lévy fits with Eq. (9).

The transverse-mass dependence of the correlation-strength parameter $\lambda$ in 40-50% centrality bin obtained from Lévy fits with Eq. (9).

The transverse-mass dependence of the correlation-strength parameter $\lambda$ in 50-60% centrality bin obtained from Lévy fits with Eq. (9).

The transverse-mass dependence of the Lévy-scale parameter $R$ in 0-10% centrality bin obtained from Lévy fits with Eq. (9).

The transverse-mass dependence of the Lévy-scale parameter $R$ in 10-20% centrality bin obtained from Lévy fits with Eq. (9).

The transverse-mass dependence of the Lévy-scale parameter $R$ in 20-30% centrality bin obtained from Lévy fits with Eq. (9).

The transverse-mass dependence of the Lévy-scale parameter $R$ in 30-40% centrality bin obtained from Lévy fits with Eq. (9).

The transverse-mass dependence of the Lévy-scale parameter $R$ in 40-50% centrality bin obtained from Lévy fits with Eq. (9).

The transverse-mass dependence of the Lévy-scale parameter $R$ in 50-60% centrality bin obtained from Lévy fits with Eq. (9).

The transverse-mass dependence of the Lévy-index of stability parameter $\alpha$ in 0-10% centrality bin obtained from Lévy fits with Eq. (9).

The transverse-mass dependence of the Lévy-index of stability parameter $\alpha$ in 10-20% centrality bin obtained from Lévy fits with Eq. (9).

The transverse-mass dependence of the Lévy-index of stability parameter $\alpha$ in 20-30% centrality bin obtained from Lévy fits with Eq. (9).

The transverse-mass dependence of the Lévy-index of stability parameter $\alpha$ in 30-40% centrality bin obtained from Lévy fits with Eq. (9).

The transverse-mass dependence of the Lévy-index of stability parameter $\alpha$ in 40-50% centrality bin obtained from Lévy fits with Eq. (9).

The transverse-mass dependence of the Lévy-index of stability parameter $\alpha$ in 50-60% centrality bin obtained from Lévy fits with Eq. (9).

The transverse-mass dependence of the normalized correlation-strength parameter $\lambda/\lambda_{\rm max}$ for 0-10% centrality intervals. For each centrality bin the data are fitted with the Gaussian function of Eq. (15). This parameterization can describe the data and is discussed in detail in Section VI B. The fit parameters with the statistic and systematic uncertainties are shown in Fig. 7.

The transverse-mass dependence of the normalized correlation-strength parameter $\lambda/\lambda_{\rm max}$ for 10-20% centrality intervals. For each centrality bin the data are fitted with the Gaussian function of Eq. (15). This parameterization can describe the data and is discussed in detail in Section VI B. The fit parameters with the statistic and systematic uncertainties are shown in Fig. 7.

The transverse-mass dependence of the normalized correlation-strength parameter $\lambda/\lambda_{\rm max}$ for 20-30% centrality intervals. For each centrality bin the data are fitted with the Gaussian function of Eq. (15). This parameterization can describe the data and is discussed in detail in Section VI B. The fit parameters with the statistic and systematic uncertainties are shown in Fig. 7.

The transverse-mass dependence of the normalized correlation-strength parameter $\lambda/\lambda_{\rm max}$ for 30-40% centrality intervals. For each centrality bin the data are fitted with the Gaussian function of Eq. (15). This parameterization can describe the data and is discussed in detail in Section VI B. The fit parameters with the statistic and systematic uncertainties are shown in Fig. 7.

The transverse-mass dependence of the normalized correlation-strength parameter $\lambda/\lambda_{\rm max}$ for 40-50% centrality interval. For each centrality bin the data are fitted with the Gaussian function of Eq. (15). This parameterization can describe the data and is discussed in detail in Section VI B. The fit parameters with the statistic and systematic uncertainties are shown in Fig. 7.

The transverse-mass dependence of the normalized correlation-strength parameter $\lambda/\lambda_{\rm max}$ for 50-60% centrality interval. For each centrality bin the data are fitted with the Gaussian function of Eq. (15). This parameterization can describe the data and is discussed in detail in Section VI B. The fit parameters with the statistic and systematic uncertainties are shown in Fig. 7.

The transverse-mass dependence of the inverse square of the Lévy-scale parameter $R$, shown in 0-10% centrality bin. The values and uncertainties of the A and B linear fit parameters with the statistic and systematic uncertainties are shown in Fig. 8.

The transverse-mass dependence of the inverse square of the Lévy-scale parameter $R$, shown in 10-20% centrality bin. The values and uncertainties of the A and B linear fit parameters with the statistic and systematic uncertainties are shown in Fig. 8.

The transverse-mass dependence of the inverse square of the Lévy-scale parameter $R$, shown in 20-30% centrality bin. The values and uncertainties of the A and B linear fit parameters with the statistic and systematic uncertainties are shown in Fig. 8.

The transverse-mass dependence of the inverse square of the Lévy-scale parameter $R$, shown in 30-40% centrality bin. The values and uncertainties of the A and B linear fit parameters with the statistic and systematic uncertainties are shown in Fig. 8.

The transverse-mass dependence of the inverse square of the Lévy-scale parameter $R$, shown in 40-50% centrality bin. The values and uncertainties of the A and B linear fit parameters with the statistic and systematic uncertainties are shown in Fig. 8.

The transverse-mass dependence of the inverse square of the Lévy-scale parameter $R$, shown in 50-60% centrality bin. The values and uncertainties of the A and B linear fit parameters with the statistic and systematic uncertainties are shown in Fig. 8.

The transverse-mass dependence of the $\widehat{R}$ parameter, shown in 0-10% centrality bin. The values and uncertainties of the linear fit parameters $\widehat{A}$ and $\widehat{B}$ with the statistic and systematic uncertainties are shown in Fig. 9

The transverse-mass dependence of the $\widehat{R}$ parameter, shown in 10-20% centrality bin. The values and uncertainties of the linear fit parameters $\widehat{A}$ and $\widehat{B}$ with the statistic and systematic uncertainties are shown in Fig. 9

The transverse-mass dependence of the $\widehat{R}$ parameter, shown in 20-30% centrality bin. The values and uncertainties of the linear fit parameters $\widehat{A}$ and $\widehat{B}$ with the statistic and systematic uncertainties are shown in Fig. 9

The transverse-mass dependence of the $\widehat{R}$ parameter, shown in 30-40% centrality bin. The values and uncertainties of the linear fit parameters $\widehat{A}$ and $\widehat{B}$ with the statistic and systematic uncertainties are shown in Fig. 9

The transverse-mass dependence of the $\widehat{R}$ parameter, shown in 40-50% centrality bin. The values and uncertainties of the linear fit parameters $\widehat{A}$ and $\widehat{B}$ with the statistic and systematic uncertainties are shown in Fig. 9

The transverse-mass dependence of the $\widehat{R}$ parameter, shown in 50-60% centrality bin. The values and uncertainties of the linear fit parameters $\widehat{A}$ and $\widehat{B}$ with the statistic and systematic uncertainties are shown in Fig. 9

The $H$ parameter that characterize the magnitude of the $\lambda/\lambda_{max}(m_T)$ and is defined in Eq. (15)., shown as functions of $N_{part}$.

The $\sigma$ parameter that characterize the width of the $\lambda/\lambda_{max}(m_T)$ and is defined in Eq. (15)., shown as functions of $N_{part}$.

The parameter $A$ of the affine linear fit to the inverse square of the Lévy-scale parameter that are defined in Eq.(12) as function of $N_{part}$

The parameter $B$ of the affine linear fit to the inverse square of the Lévy-scale parameter that are defined in Eq.(12) as function of $N_{part}$

The parameter $\widehat{A}$ of the affine linear fit to the inverse of the $\widehat{R}$ parameter that are defined in Eq.(14) as function of $N_{part}$.

The parameter $\widehat{B}$ of the affine linear fit to the inverse of the $\widehat{R}$ parameter that are defined in Eq.(14) as function of $N_{part}$.

The centrality dependence of $\alpha_0$, the $m_T$-averaged value of $\alpha$ as a function of $N_{part}$.

The Lévy scale parameter $R$ as a function of $N^{1/3}_{part}$ at $m_T=0.326$ GeV.

The Lévy scale parameter $R$ as a function of $N^{1/3}_{part}$ at $m_T=0.438$ GeV.

The Lévy scale parameter $R$ as a function of $N^{1/3}_{part}$ at $m_T=0.553$ GeV.

The Lévy scale parameter $R$ as a function of $N^{1/3}_{part}$ at $m_T=0.670$ GeV.

The Lévy scale parameter $R$ as a function of $N^{1/3}_{part}$ at $m_T=0.787$ GeV.

Monte Carlo calcuations with no mass drop assumed in 0-10%.

Monte Carlo calcuations with no mass drop assumed in 10-20%.

Monte Carlo calcuations with no mass drop assumed in 20-30%.

Monte Carlo calcuations with no mass drop assumed in 30-40%.

Monte Carlo calcuations with no mass drop assumed in 40-50%.

Monte Carlo calcuations with no mass drop assumed in 50-60%.

Monte Carlo calcuations with the best mass drop assumed in 0-10%.

Monte Carlo calcuations with the best mass drop assumed in 10-20%.

Monte Carlo calcuations with the best mass drop assumed in 20-30%.

Monte Carlo calcuations with the best mass drop assumed in 30-40%.

Monte Carlo calcuations with the best mass drop assumed in 40-50%.

Monte Carlo calcuations with the best mass drop assumed in 50-60%.

The best values of $B^{-1}_{\eta\prime}$ spectrum of the $\eta\prime$ condensate in the six centrality classes.

The centrality dependence of the best values of the in-medium mass of the $m^{*}_{\eta\prime}$.

The CL maps of the comparison of the MC to data in 0-10% (see Fig. 11).

The CL maps of the comparison of the MC to data in 10-20% (see Fig. 11).

The CL maps of the comparison of the MC to data in 20-30% (see Fig. 11).

The CL maps of the comparison of the MC to data in 30-40% (see Fig. 11).

The CL maps of the comparison of the MC to data in 40-50% (see Fig. 11).

The CL maps of the comparison of the MC to data in 50-60% (see Fig. 11).

Expected (dashed line) value of $-2\ln\Lambda$ as a function of $\kappa_{\lambda}$ for the combined fit, when all other coupling modifiers are fixed to their SM predictions.

Observed (solid line) value of $-2\ln\Lambda$ as a function of $\kappa_{2V}$ for the combined fit, when all other coupling modifiers are fixed to their SM predictions.

Expected (dashed line) value of $-2\ln\Lambda$ as a function of $\kappa_{2V}$ for the combined fit, when all other coupling modifiers are fixed to their SM predictions.

Likelihood observed contours at 68% CL (solid line) and 95% CL (dashed line) in the ($\kappa_{\lambda}$, $\kappa_{2V}$) parameter space, when all other coupling modifiers are fixed to their SM predictions. The best-fit value, denoted by a cross in the plot, corresponds to the very first entry in the table.

Likelihood expected contours at 68% CL (teal shaded region) and 95% CL (yellow shaded region) in the ($\kappa_{\lambda}$, $\kappa_{2V}$) parameter space, when all other coupling modifiers are fixed to their SM predictions. In the plot, the SM prediction is indicated by the star.

Likelihood contours at 68% (solid line) and 95% (dashed line) CL in the ($c_{hhh}$, $c_{gghh}$) parameter space, when all other Wilson coefficients are fixed to their SM values. The corresponding expected contours are shown by the inner and outer shaded regions. The SM prediction is indicated by the star, while the observed best-fit value is denoted by the black cross. The best-fit value corresponds to the very first entry in the table.

Likelihood contours at 68% (solid line) and 95% (dashed line) CL in the ($c_{hhh}$, $c_{gghh}$) parameter space, when all other Wilson coefficients are fixed to their SM values. The corresponding expected contours are shown by the inner and outer shaded regions. The SM prediction is indicated by the star, while the observed best-fit value is denoted by the black cross.

Likelihood contours at 68% (solid line) and 95% (dashed line) CL in the ($c_{hhh}$, $c_{tthh}$) parameter space, when all other Wilson coefficients are fixed to their SM values. The corresponding expected contours are shown by the inner and outer shaded regions. The SM prediction is indicated by the star, while the observed best-fit value is denoted by the black cross. The best-fit value corresponds to the very first entry in the table.

Likelihood contours at 68% (solid line) and 95% (dashed line) CL in the ($c_{hhh}$, $c_{tthh}$) parameter space, when all other Wilson coefficients are fixed to their SM values. The corresponding expected contours are shown by the inner and outer shaded regions. The SM prediction is indicated by the star, while the observed best-fit value is denoted by the black cross.

Likelihood contours at 68% (solid line) and 95% (dashed line) CL in the ($c_{H}$, $c_{H◻}$) parameter space, when all other Wilson coefficients are fixed to their SM values. The corresponding expected contours are shown by the inner and outer shaded regions. The SM prediction is indicated by the star, while the observed best-fit value is denoted by the black cross. The best-fit value corresponds to the very first entry in the table.

Likelihood contours at 68% (solid line) and 95% (dashed line) CL in the ($c_{H}$, $c_{H◻}$) parameter space, when all other Wilson coefficients are fixed to their SM values. The corresponding expected contours are shown by the inner and outer shaded regions. The SM prediction is indicated by the star, while the observed best-fit value is denoted by the black cross.

The acceptance times efficiency [in %] for the signal ggF $HH$ process as a function of the coupling modifier $\kappa_{\lambda}$ in each analysis category for the $\tau_{\mathrm{had}}\tau_{\mathrm{had}}$ channel. The other coupling modifiers affecting ggF $HH$ production are fixed to their SM predictions.

The acceptance times efficiency [in %] for the signal ggF $HH$ process as a function of the coupling modifier $\kappa_{\lambda}$ in each analysis category for the $\tau_{\mathrm{lep}}\tau_{\mathrm{had}}$ SLT channel. The other coupling modifiers affecting ggF $HH$ production are fixed to their SM predictions.

The acceptance times efficiency [in %] for the signal ggF $HH$ process as a function of the coupling modifier $\kappa_{\lambda}$ in each analysis category for the $\tau_{\mathrm{lep}}\tau_{\mathrm{had}}$ LTT channel. The other coupling modifiers affecting ggF $HH$ production are fixed to their SM predictions.

The acceptance times efficiency [in %] for the signal VBF $HH$ process as a function of the coupling modifier $\kappa_{\lambda}$ in each analysis category for the $\tau_{\mathrm{had}}\tau_{\mathrm{had}}$ channel. The other coupling modifiers affecting VBF $HH$ production are fixed to their SM predictions.

The acceptance times efficiency [in %] for the signal VBF $HH$ process as a function of the coupling modifier $\kappa_{\lambda}$ in each analysis category for the $\tau_{\mathrm{lep}}\tau_{\mathrm{had}}$ SLT channel. The other coupling modifiers affecting VBF $HH$ production are fixed to their SM predictions.

The acceptance times efficiency [in %] for the signal VBF $HH$ process as a function of the coupling modifier $\kappa_{\lambda}$ in each analysis category for the $\tau_{\mathrm{lep}}\tau_{\mathrm{had}}$ LTT channel. The other coupling modifiers affecting VBF $HH$ production are fixed to their SM predictions.

The acceptance times efficiency [in %] for the signal VBF $HH$ process as a function of the coupling modifier $\kappa_{2V}$ in each analysis category for the $\tau_{\mathrm{had}}\tau_{\mathrm{had}}$ channel. The other coupling modifiers affecting VBF $HH$ production are fixed to their SM predictions.

The acceptance times efficiency [in %] for the signal VBF $HH$ process as a function of the coupling modifier $\kappa_{2V}$ in each analysis category for the $\tau_{\mathrm{lep}}\tau_{\mathrm{had}}$ SLT channel. The other coupling modifiers affecting VBF $HH$ production are fixed to their SM predictions.

The acceptance times efficiency [in %] for the signal VBF $HH$ process as a function of the coupling modifier $\kappa_{2V}$ in each analysis category for the $\tau_{\mathrm{lep}}\tau_{\mathrm{had}}$ LTT channel. The other coupling modifiers affecting VBF $HH$ production are fixed to their SM predictions.

Expected (dashed line) value of $-2\ln\Lambda$ as a function of $\kappa_{\lambda}$ for the fit to the $\tau_{\mathrm{had}}\tau_{\mathrm{had}}$ channel, when all other coupling modifiers are fixed to their SM predictions.

Expected (dashed line) value of $-2\ln\Lambda$ as a function of $\kappa_{\lambda}$ for the fit to the $\tau_{\mathrm{lep}}\tau_{\mathrm{had}}$ channel, when all other coupling modifiers are fixed to their SM predictions.

Observed (solid line) value of $-2\ln\Lambda$ as a function of $\kappa_{\lambda}$ for the fit to the $\tau_{\mathrm{had}}\tau_{\mathrm{had}}$ channel, when all other coupling modifiers are fixed to their SM predictions.

Observed (solid line) value of $-2\ln\Lambda$ as a function of $\kappa_{\lambda}$ for the fit to the $\tau_{\mathrm{lep}}\tau_{\mathrm{had}}$ channel, when all other coupling modifiers are fixed to their SM predictions.

Expected (dashed line) value of $-2\ln\Lambda$ as a function of $\kappa_{2V}$ for the fit to the $\tau_{\mathrm{had}}\tau_{\mathrm{had}}$ channel, when all other coupling modifiers are fixed to their SM predictions.

Expected (dashed line) value of $-2\ln\Lambda$ as a function of $\kappa_{2V}$ for the fit to the $\tau_{\mathrm{lep}}\tau_{\mathrm{had}}$ channel, when all other coupling modifiers are fixed to their SM predictions.

Observed (solid line) value of $-2\ln\Lambda$ as a function of $\kappa_{2V}$ for the fit to the $\tau_{\mathrm{had}}\tau_{\mathrm{had}}$ channel, when all other coupling modifiers are fixed to their SM predictions.

Observed (solid line) value of $-2\ln\Lambda$ as a function of $\kappa_{2V}$ for the fit to the $\tau_{\mathrm{lep}}\tau_{\mathrm{had}}$ channel, when all other coupling modifiers are fixed to their SM predictions.

Summary of the expected (blue) and observed (orange) one-dimensional constraints at 68% (solid error bars) and 95% (dashed error bars) CL on the coefficients of selected HEFT (top) and SMEFT (bottom) operators describing Higgs boson interactions affecting $HH$ production through gluon-gluon fusion. The results are obtained from one dimensional scans of the profile log-likelihood assuming that all other Wilson coefficients are fixed to their SM values. The contribution from VBF $HH$ production is neglected.

Values of $N^{EXP}_{coll}$ versus $p_{T}$ as defined in Eq. 2 for three d+Au event classes, (a) $0\%–100\%$, (b) $0\%–5\%$ and (c) $60\%–88\%$. Also shown are fits to the data (solid lines) and the corresponding value $N^{GL}_{coll}$ (dashed lines).

Panels (d) to (f) show the nuclear modification factors $R^{\pi^{0}}_{dAu,EXP}$, calculated with Eq. 3, for the same event selections as in (a) to (c) together with fits to the data, (d) $0\%–100\%$, (e) $0\%–5\%$ and (f) $60\%–88\%$.

The Ratio $N^{EXP}_{coll}/N^{GL}_{coll}$ as a function of $N^{EXP}_{coll}$. Horizontal and vertical bars are the statistical uncertainties. The values for $0\%–100\%$ d+Au collisions are represented by a solid blue line, with the statistical uncertainty given as band. The scale uncertainties that are common to all data points are shown for the $0\%–100\%$ value.

Average $R^{\pi^{0}}_{dAu,EXP}$ as a function of $N^{EXP}_{coll}$. Horizontal and vertical bars are the statistical uncertainties. The values for $0\%–100\%$ d+Au collisions are represented by a solid blue line, with the statistical uncertainty given as band. The scale uncertainties(16.5%) that are common to all data points are shown for the $0\%–100\%$ value.

Average $R^{\pi^{0}}_{dAu,EXP}$ as a function of $N^{EXP}_{coll}$. Horizontal and vertical bars are the statistical uncertainties. The values for $0\%–100\%$ d+Au collisions are represented by a solid blue line, with the statistical uncertainty given as band. The scale uncertainties($16.5\%$) that are common to all data points are shown for the $0\%–100\%$ value. The systematic uncertainties on $N^{EXP}_{coll}$, $\approx16\%$, are dominated by uncertainties on the p+p data set and thus are a common scale uncertainty for all d+Au event classes.