11institutetext: Department of Physics, University of Chicago, Chicago, IL
60637, USA22institutetext: Theoretical Physics Department, Fermilab, P.O. Box
500, Batavia, IL 60510, USA33institutetext: Department of Physics and
Astronomy, University of California, Irvine, CA 92697, USA

# The impact of neutrino-nucleus interaction modeling on new physics searches

Nina M. Coyle 2,3 Shirley Weishi Li 2 Pedro A. N. Machado<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Accurate neutrino-nucleus interaction modeling is an essential requirement for
the success of the accelerator-based neutrino program. As no satisfactory
description of cross sections exists, experiments tune neutrino-nucleus
interactions to data to mitigate mis-modeling. In this work, we study how the
interplay between _near detector tuning_ and cross section mis-modeling
affects new physics searches. We perform a realistic simulation of neutrino
events and closely follow NOvA’s tuning, the first published of such
procedures in a neutrino experiment. We analyze two illustrative new physics
scenarios, sterile neutrinos and light neutrinophilic scalars, presenting the
relevant experimental signatures and the sensitivity regions with and without
tuning. While the tuning does not wash out sterile neutrino oscillation
patterns, cross section mis-modeling can bias the experimental sensitivity. In
the case of light neutrinophilic scalars, variations in cross section models
completely dominate the sensitivity regardless of any tuning. Our findings
reveal the critical need to improve our theoretical understanding of neutrino-
nucleus interactions, and to estimate the impact of tuning on new physics
searches. We urge neutrino experiments to follow NOvA’s example and publish
the details of their tuning procedure, and to develop strategies to more
robustly account for cross section uncertainties, which will expand the scope
of their physics program.

††preprint: FERMILAB-PUB-22-726-T

##  1 Introduction

The neutrino sector is the least understood of the Standard Model (SM).
Several questions about neutrinos, of both experimental and theoretical
character, remain open to this date, including the conservation of lepton
number, the mass mechanism, the mass ordering, leptonic charge-parity
violation, and puzzling experimental results which remain to be understood,
among others. One of the reasons some of these questions remain unanswered
lies in neutrinos’ small, weak-only interactions. Studying these fermions
requires special experimental setups, with intense sources and typically huge
detectors. As a consequence, the neutrino sector is widely regarded as a
promising portal to beyond the standard model (BSM) physics, while neutrino
experiments are excellent tools to search for novel BSM signatures.

Of particular interest are accelerator neutrino experiments, in which an
intense neutrino beam is produced by impinging protons on a target, and
detectors are placed near and/or far from the neutrino beam origin. While such
a setup provides a rich environment to study standard and beyond standard
physics, it also comes with challenges. These experiments measure neutrinos
through their interactions on nuclei, which are quite complex and not fully
understood processes. In the energy range between a few hundred MeV to a few
GeV, the kinematics span perturbative and non-perturbative regimes, with the
hadronic degrees of freedom being a mixture of nucleons, resonances, and
partons Kronfeld:2019nfb ; Mosel:2019vhx ; Alvarez-Ruso:2020ezu . In addition,
nuclear effects can play an important role in the final states of neutrino-
nuclei interactions Bertini:1963zzc ; Cugnon:1980zz ; Bertsch:1984gb ;
Bertsch:1988ik ; Cugnon:1996xf ; Boudard:2002yn ; Sawada:2012hk ;
Uozumi:2012fm ; Golan:2012wx ; Rocco:2020jlx ; Isaacson:2020wlx ;
Dytman:2021ohr . Although there have been many recent studies on understanding
neutrino-nucleus scattering cross sections Meyer:2016oeg ; MINERvA:2019rhx ;
Ankowski:2019mfd ; MINERvA:2020anu ; MINERvA:2020zzv ; Borah:2020gte ;
T2K:2020jav ; T2K:2021naz ; Tomalak:2021qrg ; MINERvA:2021owq ; NOvA:2021eqi ,
we still do not have a sound theoretical framework to compute the cross
sections or estimate their uncertainties Ankowski:2020qbe ;
electronsforneutrinos:2020tbf ; CLAS:2021neh . On top of that, the neutrino
flux is driven by meson production and decay, which suffer from considerable
QCD hadronization uncertainties MiniBooNE:2008hfu ; T2K:2015sqm ;
MINERvA:2016iqn .

A data-driven approach is typically sought to mitigate the large uncertainties
associated with the neutrino flux and interaction models. Commonly, neutrino
event spectra at near detectors (ND) are used to calibrate the neutrino flux
and interaction cross section models. We will refer to this procedure, which
will be explained in detail later, as ND tuning. This information is then used
in the far detector (FD) to extract oscillation parameters and other neutrino
properties. Although the primary role of NDs is to provide these inputs for
physics searches at the FDs, NDs themselves are very special. Due to the
proximity from the neutrino source and large size, they will typically have a
very high neutrino event rate, accruing millions of events over the lifetime
of an experiment, as opposed to FD statistics of $\mathcal{O}(1000)$ events or
less. Such large statistical data sets have an enormous potential to probe new
physics in the neutrino sector, particularly for models which may affect the
ND more than the FD Batell:2009di ; deNiverville:2011it ; deNiverville:2012ij
; Dharmapalan:2012xp ; Izaguirre:2013uxa ; Batell:2014yra ; Dobrescu:2014ita ;
Coloma:2015pih ; deNiverville:2016rqh ; Magill:2016hgc ; MINOS:2017cae ;
Bertuzzo:2018itn ; Falkowski:2018dmy ; Magill:2018tbb ; Ballett:2018ynz ;
deGouvea:2018cfv ; Ballett:2018uuc ; deNiverville:2018dbu ; Jordan:2018gcd ;
Jordan:2018qiy ; Ballett:2019xoj ; Harnik:2019zee ; ArgoNeuT:2019ckq ;
Altmannshofer:2019zhy ; Tsai:2019mtm ; Batell:2019nwo ; Arguelles:2019xgp ;
DeRomeri:2019kic ; deGouvea:2019wav ; Berryman:2019dme ; DUNE:2020fgq .
Because we do not have the perfect model of neutrino-nucleus interactions, and
the cross section and flux modeling is tuned with the ND data, an obvious
question is raised: _how much does the interplay between ND tuning and cross
section mis-modeling interfere with new physics searches?_

This paper aims to study the interplay between cross section mis-modeling and
searches for BSM physics. To estimate the impact on BSM searches, we analyze a
couple of illustrative scenarios, presenting their signatures before and after
tuning. Because the tuning procedure is not uniquely defined and its impact
depends on the technical details, we follow the procedure described in a
recent, detailed publication of the NOvA collaboration NOvA:2020rbg  as
closely as possible. We choose to follow NOvA’s approach because they are the
only collaboration that published the details of their tuning with codes at
the time of this work; we strongly encourage other experiments to follow
NOvA’s example and publish their tunings. For this work, we generate neutrino
events with state-of-the-art tools and apply NOvA’s tuning on an event-by-
event basis. We identify and discuss how the smoking gun signatures of new
physics models are affected. While some signatures, such as sterile neutrino
oscillations, are robust against tuning, mis-modeling of neutrino interactions
can bias the experimental sensitivity. On the other hand, the sensitivity to
other types of BSM scenarios can be overwhelmed by cross section mis-modeling.

This is the first quantitative study exploring the impact of tuning on new
physics searches. We identify conceptual lessons that should be applicable to
other tuning procedures, detector set ups, analysis details, and BSM
scenarios. For instance, we expect our findings to hold to some extent even in
BSM searches that do not employ a ND tuning but instead rely on a simultaneous
fit to near and far detector data (see, e.g., Ref. MINOS:2017cae ). Our goal
is not to faithfully reproduce any experimental analyses. We hope our results
will encourage theorists and experimentalists to consider the impact of ND
tuning on new physics scenarios and to carefully estimate systematic
uncertainties related to neutrino fluxes and cross sections.

##  2 Near detector tuning

Near detectors in long-baseline neutrino oscillation experiment are designed
to mitigate the uncertainties from both the neutrino beams and neutrino-
nucleus cross sections. Conceptually, this procedure may seem straightforward:

$P(\nu_{\alpha}\rightarrow\nu_{\beta};E_{\nu})=C\frac{\displaystyle\frac{dN^{\mathrm{FD}}_{\beta}}{dE_{\nu}}\bigg{/}\sigma_{\beta}(E_{\nu})}{\displaystyle\frac{dN^{\mathrm{ND}}_{\alpha}}{dE_{\nu}}\bigg{/}\sigma_{\alpha}(E_{\nu})},$
(2.1)

where $P(\nu_{\alpha}\rightarrow\nu_{\beta};E_{\nu})$ is the oscillation
probability the experiment is trying to measure, $C$ is a constant accounting
for detector sizes and distances to the source, the superscripts ND and FD
denote near and far detectors, and $dN_{\alpha,\beta}/dE_{\nu}$ are the event
rates of $\nu_{\alpha,\beta}$ in terms of true neutrino energy. The
differences between $\sigma_{\alpha}$ and $\sigma_{\beta}$, if
$\alpha\neq\beta$, could in principle be computed theoretically
Tomalak:2021qrg . The oscillation probability can be directly derived from the
unoscillated event rate $dN_{\alpha}/dE_{\nu}$ measured at the ND and the
oscillated event rate $dN_{\beta}/dE_{\nu}$ measured at the FD.

In reality, this “logical division” is not practical because one cannot
directly measure neutrino energy, the near and far detectors have different
systematics, and even the unoscillated neutrino fluxes are different at the
near and far detectors solely due to the different solid angles. We can
appreciate this more concretely by studying the following equation,

$\frac{N_{\rm FD}}{N_{\rm ND}}(E_{\mathrm{reco}})=C\frac{\displaystyle\int
dE_{\nu}\frac{d\phi^{\rm
FD}_{\alpha}}{dE_{\nu}}P(\nu_{\alpha}\to\nu_{\beta};E_{\nu})\sigma_{\beta}(E_{\nu})\mathcal{M}^{\mathrm{FD}}_{\beta}(E_{\nu},E_{\mathrm{reco}})}{\displaystyle\int
dE_{\nu}\frac{d\phi^{\rm
ND}_{\alpha}}{dE_{\nu}}\sigma_{\alpha}(E_{\nu})\mathcal{M}^{\mathrm{ND}}_{\alpha}(E_{\nu},E_{\mathrm{reco}})},$
(2.2)

where $N_{\rm FD/ND}$ encode the reconstructed neutrino event spectra,
$d\phi^{\mathrm{FD/ND}}_{\alpha}/dE_{\nu}$ are the fluxes at the far and near
detectors without oscillation, $\sigma_{\alpha,\beta}(E_{\nu})$ are the total
cross sections, and $\mathcal{M}^{\mathrm{FD/ND}}_{\alpha}$ are the migration
matrices. The challenge here is that experiments measure the left-hand side of
Eq. (2.2) and they need to infer the oscillation probability $P$ on the right-
hand side. The main difficulty is encoded in the term
$\sigma_{\alpha}(E_{\nu})\mathcal{M}_{\alpha}(E_{\nu},E_{\mathrm{reco}})$,
where the reconstruction of the true neutrino energy depends on the details of
neutrino-nucleus interaction, as well as detector responses to different
final-state particles. An obvious example is neutrons, for which both the
modeling of neutrons produced by neutrino interactions and the corresponding
detector responses are relevant. Currently, predictions on the number of
outgoing nucleons in a neutrino-nucleus scattering event, as well as their
energy and isospin, differ drastically among generators Ankowski:2019mfd ;
CLAS:2021neh . In addition, neutron detector responses suffer from significant
uncertainties due to neutron propagation and event reconstruction
Friedland:2018vry ; CAPTAIN:2019fxo ; Friedland:2020cdp .

Because of these complications, oscillation experiments adopt a near-detector
tuning procedure. Assuming SM physics, one can predict the measured neutrino
energy spectrum in the ND with an accelerator neutrino beam simulation, an
event generator such as GENIE Andreopoulos:2009rq , NuWro Golan:2012rfa ,
GiBUU Buss:2011mx , or ACHILLES Isaacson:2022cwh  that simulates neutrino-
nucleus interaction cross sections and final states, and a detector simulation
predicting the migration matrix. When the predicted and measured neutrino
spectra disagree, one can modify the cross section simulations until they
match. Because of the nature of such calculations and the complexities of
these simulation packages, there is no unique, agreed-upon way to tune the
models. One can vary the model parameters, adopt alternative models, or take a
model-agnostic approach and add more degrees of freedom.

In addition, it is well appreciated that neutrino beams have sizable
uncertainties and they interfere with cross section uncertainties. Different
experiments also treat the flux tuning and the cross section tuning
differently. NOvA, for instance, only uses hadronic production data, in-situ
measurements of horn position and current, beam parameters, etc., and MINERvA
neutrino-electron scattering data to tune their flux prediction NOvA:2021nfi .
T2K, on the other hand, uses ND neutrino-nucleus scattering data in addition
to auxiliary data to tune their flux and cross section models at the same time
T2K:2021xwb .

In this work, we follow the tuning procedure outlined by NOvA NOvA:2018gge .
The NOvA experiment is a long-baseline experiment comprised of two
scintillator detectors (CH2) placed along a $\nu_{\mu}/\bar{\nu}_{\mu}$
beamline produced by the NuMI facility at Fermilab: a near and a far detector
placed 1 km and 810 km from the beam source, respectively. The 14 kton FD
observes the muon neutrino beam after long-baseline oscillations, intended to
measure oscillation parameters including $\Delta m_{32}^{2}$ and
$\theta_{23}$. The smaller 0.3 kton ND is nearly identical to the FD to
minimize systematic uncertainties.

The ND tune procedure is detailed in Ref. NOvA:2020rbg . Neutrino interactions
with the material in the ND are first generated using GENIE v2.12.2 with the
default models and parameters. Then, the following changes are applied to the
GENIE simulation based on auxiliary theoretical and experimental studies,
i.e., NOvA ND data is not used for this step:

  * •

adjusting the value of axial mass, $m_{A}$, from 0.99 to 1.04 GeV, based on
recent re-analysis Meyer:2016oeg  of neutrino-deuterium scattering data;

  * •

modifying the momentum distributions of the initial nucleons for quasi-elastic
scattering, based on a MINERvA study Gran:2017psn ;

  * •

lowering the magnitude of neutrino-scattering in non-resonance pion production
regime by 57%, motivated by re-analysis of old bubble chamber data
Rodrigues:2016xjj ;

  * •

suppressing delta resonance production in low-$Q^{2}$ region, motivated by
measurements by MiniBooNE MiniBooNE:2010eis , MINOS MINOS:2014axb , MINERvA
MINERvA:2016sfc ; MINERvA:2017okh , and T2K T2K:2019yqu .

All these changes are to improve the baseline model in GENIE.

After these changes are applied to GENIE, there are still large discrepancies
between the measured neutrino spectrum in NOvA ND and the simulated spectrum.
The last step of NOvA’s ND tune is the crucial step of which we are studying
the effect. To understand how it works, let us first define the kinematics of
neutrino interactions, then explain how kinematic variables are measured in
NOvA. Theoretically, an incoming neutrino with energy $E_{\nu}$ produces an
outgoing lepton with energy $E_{l}$ and a hadronic system. The four-momentum
transfer is $(q_{0},\vec{q})$ with $q_{0}=E_{\nu}-E_{l}$. In NOvA, as in any
experiment, the neutrino energy is not directly observable and is instead
measured through $E_{\nu}^{\mathrm{reco}}\equiv
E_{l}^{\mathrm{reco}}+q_{0}^{\mathrm{reco}}$. For charged-current
interactions, which are the signal channels for oscillation analyses, $E_{l}$
can be measured relatively well either from the total scintillation light
associated with the electron or muon or from the muon range, i.e.,
$E_{l}^{\mathrm{reco}}\simeq E_{l}$. The transferred energy, $q_{0}$, gets
distributed between kinetic energy for knocked out nucleons, total energy for
mesons, and binding energy of the initial nucleus. The sum of nucleon kinetic
energy, $T$, and meson total energy $E$, roughly proportional to the amount of
scintillation light detected, is called the hadronic energy, ${E_{\rm had}}$,
i.e.,

$q_{0}^{\mathrm{reco}}\simeq{E_{\rm
had}}=\sum_{i}^{\textrm{nucleons}}T_{i}+\sum_{j}^{\textrm{mesons}}E_{j}.$
(2.3)

We include a 30% smearing on ${E_{\rm had}}$, which roughly accounts for
NOvA’s detector response. Lastly, NOvA tune uses the reconstructed three-
momentum transfer, $|\vec{q}\,^{\mathrm{reco}}|^{2}=2({E_{\rm
had}}+E_{\mu})(E_{\mu}+p_{\mu}\cos\theta_{\mu})-m_{\mu}^{2}+E_{\mathrm{had}}^{2}$,
where $E_{\mu}$, $p_{\mu}$, $\theta_{\mu}$ are the energy, momentum, and the
angle with respect to the beam of the muon, for $\nu_{\mu}$ charged-current
events. After all the modifications to GENIE described above, in the NOvA ND
tune, one then compares the simulated event distribution to the measured one
in the 2D plane of reconstructed variables
$(|\vec{q}\,^{\mathrm{reco}}|,q_{0}^{\rm reco})$ with 20 bins in
$|\vec{q}\,^{\mathrm{reco}}|$ and 16 bins in $q_{0}^{\rm reco}$. With the
assumption that all the discrepancies between the simulated and measured
distributions are due to mis-modeling of the meson-exchange current (MEC), one
multiplies the MEC event rate predicted by GENIE in every bin by a weight,
i.e., “tune,” such that the tuned prediction matches the measurement. The
total number of free parameters in this tuning is 200. See Fig. 1 for a
schematic illustration.

Figure 1: Schematic illustration of NOvA’s tuning. On the left, we show the
NOvA data versus the output of default GENIE 2.12.2 (colored histograms). The
tuning happens in two steps, first, a few adjustments to GENIE are made, such
as redefining $m_{A}$, and then the meson-exchange current (MEC) component of
the cross section is reweighted to data (see text for details). The result is
the plot on the right, where the output of tuned GENIE matches the data well.
Details have been omitted for clarity. These figures are adapted from Ref.
NOvA:2020rbg .

##  3 Sterile neutrinos

To estimate the impact of the tuning procedure on BSM searches, we choose to
study two illustrative models: eV-scale sterile neutrinos Kainulainen:1990bn
and light neutrinophilic scalars Kelly:2019wow . The reason for choosing these
two models is their different experimental signatures. In this section, we
start with sterile neutrinos with masses around the eV scale, for which
oscillations of active neutrinos to sterile neutrinos will look like wiggles
in the ND neutrino energy spectrum.

###  3.1 Model description

For two decades, the standard three neutrino paradigm has been confronted with
experimental anomalies that could indicate oscillations due to a mass
splitting much larger than the measured solar and atmospheric ones. These
anomalies are fourfold: LSND LSND:2001aii , MiniBooNE MiniBooNE:2020pnu ,
reactor Mueller:2011nm ; Mention:2011rk ; Huber:2011wv , and gallium
GALLEX:1997lja ; SAGE:1998fvr ; Abdurashitov:2005tb ; Kaether:2010ag ;
Giunti:2010zu ; Barinov:2021asz  anomalies. All these anomalies seem
consistent with oscillations at a baseline distance much shorter than
expected, pointing to a new mass splitting $\Delta
m^{2}_{41}\sim\mathcal{O}(1~{}{\rm eV}^{2})$. Nevertheless, direct
ALEPH:1993pqw ; OPAL:1994kgw ; L3:1998uub  and indirect DELPHI:2000wje
measurements of the $Z$ invisible width at LEP would imply that this state
does not interact via weak interactions and thus has no standard model gauge
quantum number whatsoever: a _sterile neutrino_.

Short baseline oscillations may be approximated by

$P(\nu_{\alpha}\to\nu_{\beta})\simeq\delta_{\alpha\beta}-4|U_{\alpha
4}|^{2}(\delta_{\alpha\beta}-|U_{\beta 4}|^{2})\sin^{2}\left(\frac{\Delta
m^{2}_{41}L}{4E_{\nu}}\right),$ (3.4)

where $U_{\alpha i}$ denotes the extended 4$\times$4 PMNS matrix, $\Delta
m^{2}_{41}\equiv m_{4}^{2}-m_{1}^{2}$, and $E_{\nu}$ is the neutrino energy.
Explaining the LSND and MiniBooNE anomalies requires $\nu_{\mu}\to\nu_{e}$
transitions, and thus a nonzero $U_{e4}$ and $U_{\mu 4}$. This, in turn,
implies short baseline $\nu_{\mu}$ and $\nu_{e}$ disappearance. This relation
between appearance and disappearance channels leads to a strong tension
between data sets under the interpretation of sterile neutrino oscillations
Dentler:2018sju ; Moulai:2019gpi ; Giunti:2019aiy . This tension is mainly
driven by the following experiments: MINOS/MINOS+ MINOS:2020iqj  and
IceCubeIceCube:2016rnb ; IceCube:2020tka  for $\nu_{\mu}$ disappearance; Daya
Bay MINOS:2020iqj , as well as solar SNO:2008gqy ; Super-Kamiokande:2016yck ;
BOREXINO:2018ohr , short baseline RENO:2020hva ; DANSS:2018fnn , and
radioactive source experiments Barinov:2021asz  for $\nu_{e}$ disappearance;
and LSND and MiniBooNE for $\nu_{\mu}\to\nu_{e}$ appearance.

Regardless of this tension, testing short baseline oscillations with
disappearance channels has become increasingly important in the last decade.
Here, we are interested in the $\nu_{\mu}$ disappearance channel at NOvA ND
for which

$P(\nu_{\mu}\to\nu_{\mu})\simeq
1-\sin^{2}2\theta_{\mu\mu}\sin^{2}\left(\frac{\Delta
m^{2}_{41}L}{4E_{\nu}}\right),$ (3.5)

where we have defined the effective $\nu_{\mu}$ disappearance angle
$\sin^{2}2\theta_{\mu\mu}\equiv 4|U_{\mu 4}|^{2}(1-|U_{\mu 4}|^{2})$.
Oscillations in NOvA ND, located about 1 km from the beam target, are depicted
for illustrative oscillation parameters in Fig. 2. Clearly, a sterile neutrino
could induce an oscillatory pattern in NOvA’s neutrino spectrum at the ND. The
question we want to ask is how much of this oscillatory pattern will be
affected by NOvA’s tuning procedure, and what the corresponding impact is on
NOvA’s sensitivity to sterile neutrinos.

Figure 2: Oscillation probability in the $\nu_{\mu}$ disappearance channel for
illustrative sterile neutrino oscillation parameters as a function of the true
neutrino energy $E_{\nu}$ (blue line). We also show the shape of the NOvA
$\nu_{\mu}$ flux in neutrino mode for reference (gray histogram) NOvAflux .

###  3.2 Analysis

Let us start with the general approach used to perform analyses in this work.
To mimic a sterile neutrino search in NOvA, we first use an event generator
$\mathcal{G}$ to generate mock data for both ND and FD that contain a sterile
neutrino signal. We then use either the same or a different event generator
$\mathcal{G}^{\prime}$ to tune a SM event set to the mock data in the ND. The
idea behind using a different generator to fit the mock data is to account for
possible mis-modeling of neutrino-nucleus scattering. Lastly, we do a chi-
square analysis using the mock data and the tuned generator
$\mathcal{G}^{\prime}$ to see if we can recover the sterile signals and to
examine how the sensitivity changes because of the tuning.

More concretely, we use the GENIE v2.12.2 event generator to produce mock
data. We produce a set of $\nu_{\mu}$-carbon charged current events and
reweight them according to the NOvA flux and the oscillation probability
$P(\nu_{\mu}\to\nu_{\mu})$ at the ND. We account for the varying decay lengths
of the pions in the decay pipe. Following Ref. NOvAdecay , we take a random
oscillation length $L$ from an exponential distribution that has a maximum at
$L=1$ km and extends to $L=0.3$ km. The reconstructed neutrino energy is
obtained following the procedure discussed in Sec. 2: $E_{\nu}^{\rm
reco}=E_{\mu}+q_{0}^{\rm reco}$, where $q_{0}^{\rm reco}$ is the hadronic
energy $E_{\rm had}$ smeared by 30%.

To fit the data, we study several cases. We use either GENIE v2.12.2 or NuWro
v19.02.2 ***For the MEC component, GENIE adopts the so-called “Empirical MEC”
or “Dytman” model Katori:2013eoa , while NuWro uses the “Nieves” model
Nieves:2011pp . to generate the SM spectra, reweighting the events as above.
We fit mock data both with and without performing the tuning procedure to
understand the extent of the impact of the tuning, especially its interplay
with cross section mis-modeling, on the experimental sensitivity.

For the tuning, we follow NOvA’s procedure as discussed in Sec. 2; however,
because we are working with two sets of simulated events, we perform the tune
entirely in the true $(|\vec{q}|,q_{0})$ plane. When considering sterile
neutrinos in the fit generator $\mathcal{G}^{\prime}$, one could either first
reweight according to the oscillation probability and then tune, or the other
way around. Performing the analysis by first reweighting ensures us that when
we use the same generator for both mock data and fit model, and we compare the
same parameters for sterile neutrinos, the tuning does not make any changes to
the generated events because, by design, the mock data and the fit generator
predict the same ND spectra.

To appreciate the impact of the tuning, we present Fig. 3. In the left panel,
we compare the mock data (black points) against the fit spectra using GENIE
without (blue) and with (red) tuning. Mock data was generated with a sterile
neutrino signal, assuming $\Delta m^{2}_{41}=5$ eV2 and
$\sin^{2}2\theta_{\mu\mu}=0.2$, while the fit spectra have no sterile signal.
The top panel shows the spectra and the bottom panel shows the ratios of the
data to fit spectra. The blue lines reflect the situation that one would
intuitively expect in a sterile neutrino analysis. In the top panel, we see
that the data rate is lower than the SM rate because of the existence of a
sterile neutrino; in the bottom panel, the blue line shows the characteristic
oscillatory feature of a sterile neutrino signal. The wiggles look slightly
different from Fig. 2 due to energy and baseline smearing. The red lines,
meanwhile, illustrate the effect of tuning. The top panel shows that the
tuning does not “fix” all the discrepancies between the fit generator
prediction and the data. This might be surprising because the tuning has many
degrees of freedom; however, the tuning is performed in the $(\vec{q},q_{0})$
plane, and these observables do not correlate with $E_{\nu}$ in a
straightforward way. In fact, the red line in the bottom panel shows that if
our modeling of the SM neutrino-carbon interactions is perfect (when we use
GENIE to model both the data and the fit spectra), the tuning almost perfectly
preserves the sterile neutrino feature. For more detailed discussions, see
Appendix A.1.

The right panel of Fig. 3 shows the same mock data but with the fit generator
being NuWro. This is a more realistic situation, where our modeling of
neutrino-carbon interactions, i.e., by the fit generator NuWro, is different
from the true interactions, i.e., by the mock data generator GENIE. To mimic
the situation in NOvA ND, in which the pre-tune model (NuWro) predicts fewer
overall events than the data (GENIE), we rescale the model weights relative to
the data by an overall factor of 0.9. For reference, we plot in the bottom
panel the ratio between SM GENIE and SM NuWro spectra (grey dashed line). We
can see from the red and blue lines in the bottom panel that some oscillatory
features are still present in the ratio between mock data and fits, but less
prominent than in the left panel. This can be partially attributed to the
intrinsic differences between baseline generators, without the presence of any
new physics, which can be seen in the grey dashed line. We can already see
that this may cause a bias in the experimental sensitivity to sterile
neutrinos.

Figure 3: ND distributions of $E_{\nu}^{\mathrm{reco}}$ for sterile data and
SM expectation for $\Delta m_{41}^{2}=5$ eV2 and
$\sin^{2}(2\theta_{\mu\mu})=0.2$. In the left panel, the data and model are
both generated using GENIE; in the right panel, the model is generated using
NuWro. The top panel shows a histogram of events for the SM (blue), mock data
(black points), and tuned SM (red) at the ND, while the lower panel shows the
ratio ND data / ND SM for both the tuned (red) and untuned (blue) SM
expectation.

To estimate NOvA’s sensitivity, we build a covariance matrix, accounting for
the following systematics: 20% overall normalization, 4% normalization of ND
relative to FD, 2% correlated near-to-far spectral uncertainty (bin-by-bin),
and fully uncorrelated 2% bin-to-bin uncertainty on both near and far
detectors. Our covariance matrix is meant to capture the general features of
NOvA’s sensitivity, and not to reproduce precisely the experimental results.
The sensitivity is estimated by a chi-square function $\chi^{2}=(D-F)\cdot
C^{-1}\cdot(D-F)$, where $D$ and $F$ are the data and fit spectra, and $C$ is
the covariance matrix including statistical uncertainties.

We perform a sensitivity analysis for the tuned and untuned case after adding
a sterile neutrino signal with $\sin^{2}2\theta_{\mu\mu}=0.1$ and two choices
of mass splittings, $\Delta m_{41}^{2}=5.0$ eV2 and 2.0 eV2, as shown in the
left and right panels of Fig. 4. Again, we generate the mock data with GENIE
and estimate the sensitivities for four cases: GENIE without tuning (shaded
region), GENIE with tuning (thin line), NuWro without tuning (blue), and NuWro
with tuning (red). Let us start with the left panel. The grey shaded region
indicates an ideal analysis, where our fit generator and the data generator
are the same, and there is no tuning. In this case, we recover the true input
parameter, and the systematic uncertainties determine the region size. Note
that because the FD has very few events, it contributes very little to the
$\chi^{2}$, and the grey region would look almost the same if we did a ND-only
analysis. We also show the impact of the tuning procedure if one had the
correct model of neutrino-carbon interaction (thin line). We can see that
tuning indeed enlarges the region slightly, but the changes are small. The
blue region shows that when we mis-model neutrino-carbon interaction and we do
not tune our generator, we would still identify a sterile neutrino signal, but
the fit parameters would be biased by more than 2$\sigma$. Lastly, the red
region shows that if we adopt the ND tuning, we improve the accuracy of the
fitted parameters. This could seem to indicate that the tuning procedure makes
the sensitivity more robust to mis-modeling.

The right panel of Fig. 4 tells a different story. With just a different
$\Delta m_{41}^{2}$, we see that neither the untuned nor the tuned regions
include the real input values at 2$\sigma$. The untuned region is consistent
with no sterile neutrinos, and the tuned region prefers large mixing angles.
This illustrates that cross section mis-modeling can significantly bias the
experimental sensitivities and that the tuning procedure does not fix this
issue.

Figure 4: Sensitivities to sterile neutrinos for different choices of event
generators and tuning, assuming sterile neutrino parameters
$\sin^{2}2\theta_{\mu\mu}=0.1$ and $\Delta m^{2}_{41}=5$ eV2 (left panel) and
$\Delta m^{2}_{41}=2$ eV2 (right panel).

##  4 Light neutrinophilic scalars

The second scenario we analyze is the light neutrinophilic scalar model
proposed in Refs. Berryman:2018ogk ; Kelly:2019wow . Its experimental
signature consists of an excess of missing transverse momentum. This is a good
example of a class of models in which the experimental signature depends on
certain aspects of the interaction kinematics on an event-by-event basis. The
correlations between these kinematic variables and the tuning variables
$q_{0}$ and $|\vec{q}|$ can be nontrivial. By studying a representative
signature, we hope to uncover common lessons for these kinematic searches.

###  4.1 Model description

Light neutrinophilic scalars could mediate interactions between neutrinos and
dark matter candidates. These scalars would couple to neutrinos via the
effective operator,

$\mathcal{O}=\frac{(L_{\alpha}H)(L_{\beta}H)}{\Lambda_{\alpha\beta}^{2}}\phi+{\rm
h.c.}\to\frac{1}{2}\lambda_{\alpha\beta}\nu_{\alpha}\nu_{\beta}\phi+{\rm
h.c.},$ (4.6)

where $\phi$ is the neutrinophilic scalar itself; $H$ and $L_{\alpha}$ are the
Higgs and lepton doublets; $\alpha,\beta=e,\mu,\tau$; $\Lambda_{\alpha\beta}$
is the scale of the dimension-6 operator;
$\lambda_{\alpha\beta}=v^{2}/\Lambda_{\alpha\beta}^{2}$ is an effective
coupling between the $\phi$ and neutrinos after electroweak symmetry breaking;
and lastly, $v=246$ GeV is the Higgs vacuum expectation value. A sum over
flavor indices is implicit. In addition to this neutrino interaction term,
$\phi$ could also couple to dark matter, so it serves as a portal between
neutrinos and the dark sector.

In this scenario, when a neutrino scatters via charged current interactions,
it may radiate a $\phi$, which would either leave the detector or decay
invisibly to neutrinos or dark matter, thus constituting missing momentum.
From Eq. (4.6), we note that $\phi$ carries lepton number; thus, when emitted,
$\phi$ changes a neutrino to an antineutrino and vice versa. We depict the
Feynman diagram inducing this transition in Fig. 5. Note that we have not been
very explicit on the hadronic current, as the neutrino may scatter on
nucleons, quarks, and even two-body currents. This signature has been named
the _mono-neutrino_ signal, in tandem with collider mono-$X$ searches
Abercrombie:2015wmb .

$W$$\phi$$\overline{\nu}_{\alpha}$$\nu_{\beta}$$\ell_{\beta}^{-}$p, n, …
Figure 5: Feynman diagram relevant to the mono-neutrino signal. Note that the
hadronic current can go beyond nucleon level (e.g. two-body currents, deep
inelastic scattering), so we leave it general.

The typical experimental signature of this scenario will be an excess of
missing transverse momentum, $\not{p}_{T}$, compared to usual neutrino
scattering events. The missing transverse momentum is defined as the magnitude
of the sum of the transverse momenta of the visible particles,

$\not{p}_{T}=\left|\sum_{\mathrm{i,vis}}(\vec{p}_{T})_{i}\right|.$ (4.7)

While the $\phi$ emission will necessarily lead to a missing momentum,
standard neutrino-nucleus interactions could also display relatively large
missing transverse momenta. There are three main reasons for that. If a
neutrino interacts with a nucleon inside a nucleus, the struck nucleon is not
at rest due to Fermi motion. The Fermi momentum is typically of the order of
250 MeV or so, leading to an unavoidable spectrum of $\not{p}_{T}$. Moreover,
neutrinos scatter with particles inside a dense nuclear medium. As particles
propagate throughout the nucleus, secondary scatterings may occur, referred to
as final state interactions. This effect is fairly common and may lead to
large $\not{p}_{T}$. Lastly, outgoing hadrons, particularly neutrons, could
escape detection and thus lead to $\not{p}_{T}$.

###  4.2 Analysis

To simulate a set of mono-neutrino events, we use MadGraph5 Alwall:2011uj  to
generate $\overline{\nu}_{\mu}+n\to\mu^{-}+\phi+p$ events, treating the proton
and neutron as elementary particles. Although this neglects Fermi motion and
final state interactions, the emission of a $\phi$ dominates the missing
transverse momentum distribution. We choose to work with an antineutrino beam
because the final-state protons are visible and can be used to reduce SM
backgrounds, as we will discuss later. We generate events at fixed neutrino
energies, then weight the events by the convolution of the mono-neutrino cross
section and the NOvA antineutrino flux at the given energy for each event. In
this analysis, we use NuWro to simulate SM interactions, which accounts for
Fermi motion and final state interactions. Here we treat neutrons as missing
energy and include a 30% smearing on hadronic momenta. Mock data is obtained
by combining the signal and background events, producing a simulated expected
missing momentum spectrum at the ND for this model. To account for mis-
modeling of neutrino-nucleus interactions, we use GENIE to simulate our SM fit
background.

In the left panel of Fig. 6, we show the missing momentum spectra for the mock
data background (black points), the fit model (blue), and the mono-neutrino
signal (red) for coupling $\lambda_{\mu\mu}=0.5$ and scalar mass
$m_{\phi}=0.5$ GeV. As we can see, the mono-neutrino signal is more than two
orders of magnitude below the background. This is partially due to the
multiplicity of the final states of the signal. The extra phase space factor
of three- compared to two-body final states leads to a suppression of about
two orders of magnitude on the cross section. This low signal rate forces us
to impose experimental cuts to improve the signal-to-background ratio.

In antineutrino mode, the signal events are given by
$\overline{\nu}_{\mu}+n\rightarrow\mu^{-}+\phi+p$. The visible final states in
the detector are the muon and the proton. In contrast, we can think of the
background events naively as given by
$\overline{\nu}_{\mu}+p\rightarrow\mu^{+}+n$. The presence of a neutron could
lead to large $\not{p}_{T}$, which would mimic the signal. While neutrons are
hard to reconstruct, those above 100 MeV are likely to have hard interactions
with nuclear matter, leading to hadronic activity that allows us to identify
its presence. Given this, we cut events with neutrons above 100 MeV. †††This
is consistent with setting all the neutron energy as missing energy. We assume
that, even though NOvA can identify the presence of a neutron above 100 MeV,
they cannot reconstruct its energy, and therefore its energy still goes to
missing energy. Moreover, we only select events with exactly one visible
proton, that is, protons with kinetic energy above 100 MeV, and no pions at
all. We include a 30% smearing on hadronic momenta.

The result of these cuts can be seen in the right panel of Fig. 6. While the
signal-to-background ratio has improved considerably, a new issue arises: the
background-only mock data is significantly different from the fit model
(blue). Under further consideration, this should not be surprising because the
more exclusive we look at cross section predictions, the larger the expected
theoretical uncertainties and discrepancies. To understand the extent to which
this large discrepancy could be addressed by ND tuning, we perform the tuning
before (gray) and after (black) experimental cuts. We can clearly see that
neither options can reproduce the mock data spectrum. But, more importantly,
the tuned spectra depend strongly on the cuts used, so even if a tune is
performed with one event sample and measured with another, this dependence may
introduce biases in the experimental analyses. On top of that, the differences
between tuned model and data are larger than the signal rate itself. For more
details on discrepancy between generators and the impact of the tuning to
mono-neutrino signal, see Appendix A.2.

Figure 6:  Missing $p_{T}$ distributions before (left) and after (right) cuts.
Mock SM data (black points) was generated with NuWro, while GENIE was used as
fit model (blue). The mono-neutrino signal is given in red. In the right panel
we also show the predicted $/\\!\\!\\!p_{T}$ spectra for tuned fit model when
the tuning is performed before cuts (gray) or after cuts (black).

This large discrepancy between mock data and tuned fit model already shows us
that the mis-modeling of neutrino-nucleus interactions may render an
experiment completely unable to probe certain BSM physics scenarios.
Regardless, we can perform a statistical analysis of the experimental
sensitivity to appreciate quantitatively how relevant cross section modeling
is to this BSM search. We perform the analysis by employing a similar
covariance matrix to that built in the previous section, but excluding the FD
data, and with the uncertainties now applicable to bins of $\not{p}_{T}$. For
systematics we take an overall normalization uncertainty of 20% and a bin-to-
bin uncorrelated uncertainty of 5%; our chi-square is $\chi^{2}=(D-F)\cdot
C^{-1}\cdot(D-F)$, where $D$ and $F$ are the data and fit model binned in
$\not{p}_{T}$. While this is not a realistic treatment of experimental
sensitivity to missing momentum signals, it suffices for our purposes.

Figure 7 shows our results, assuming the presence of a mono-neutrino signal
with coupling $\lambda_{\mu\mu}=0.5$ and $m_{\phi}=0.5$ GeV. In grey, we
present the allowed region when using the same generator for the signal and
the fit model, which would correspond to the case of perfect modeling of
neutrino-nucleus interactions. We also show the allowed region for data and
fit model using different generators, which would correspond to cross section
mis-modeling. We see that without tuning (blue), the fit prefers no new
physics, as the allowed region includes $\lambda_{\mu\mu}=0$. In fact, the fit
excludes the true parameter point at more than 5$\sigma$. With tuning (red),
the results actually get worse, as the fit strongly prefers no new physics. We
can conclude that the discrepancy between mock data and the tuned fit model
dominates the experimental sensitivity.

This shows that while near detectors have the potential to probe new physics,
without proper modeling of neutrino-nucleus interactions, we may lose this
potential. Tuning is not the solution.

Figure 7:  Sensitivity for a mono-neutrino signal with $\lambda_{\mu\mu}=0.5$
and $m_{\phi}=0.5$ GeV for different choices of event generators and tuning.
Because the $2\sigma$ tuned region is not visible on the figure, we show the
$5\sigma$ region (in the lower left corner).

##  5 Discussions and conclusions

In this paper, we have estimated the impact of the ND neutrino-nucleus cross
section tuning procedure on BSM signatures at neutrino experiments. To perform
a realistic study, we follow very closely the NOvA collaboration tuning
procedure, and we simulate neutrino scattering events with state-of-the-art
event generators GENIE and NuWro. We study two illustrative BSM scenarios: eV-
scale sterile neutrinos, which would manifest via wiggles in the $\nu_{\mu}$
survival oscillation probability; and neutrinophilic scalars, which could be
emitted by neutrinos in charged current interactions leading to large amounts
of missing transverse momentum.

Our results show that the interplay between cross section mis-modeling, ND
tuning, and sensitivity to new physics is far from trivial. Concretely, the
wiggles in the ND neutrino energy spectrum induced by eV-scale sterile
neutrino oscillations largely remain after tuning. That is, the tuning of the
neutrino-nucleus cross section cannot mimic sterile neutrino signatures.
Nevertheless, mis-modeling neutrino-nucleus interactions can lead to
discrepancies between theory and data that may introduce a bias on the
experimental sensitivity.

The situation for neutrinophilic scalars is quite different. First, the beyond
standard model signature is significantly lower than the standard model
background. This forces us to impose experimental cuts to improve the signal-
to-background ratio. In doing so, we go to regions of parameter space in which
the background is suppressed but most of the signal remains; however, in these
regions the tuning of the fit model is inefficient at reproducing the data,
leading to large discrepancies between data and theory. We have seen that for
the case of neutrinophilic scalars, this discrepancy dominates the
experimental sensitivity, rendering it unrealistic.

In both BSM scenarios studied here, the mis-modeling of neutrino-nucleus
interactions can lead to significant biases on the experimental sensitivity,
regardless of the tuning procedure. One may wonder if this is a widespread
issue in BSM searches at neutrino experiments, such as T2K, HK, and DUNE. It
is reasonable to presume that a tuning based on the correct physics would be
able to resolve mis-modeling; however, ad hoc tunings will most likely fail to
capture the correct physics of neutrino-nucleus interactions, e.g. by lacking
important correlations or by being arbitrarily limited to specific processes.
Without the correct physics, a tuning will only be able to reproduce
experimental data _by chance_. Agreement between an ad hoc tuning and data in
one observable does not necessarily imply that any other observable will be
well described by the tuning, as can be seen for instance in recent MINERvA
analyses MINERvA:2020zzv . It should therefore not come as a surprise that a
tuning that is not correlated with a BSM signature and thus does not wash it
out is also not capable of addressing mis-modeling.

One may think that practically, a conservative approach would be to assign a
large systematic error that includes the possibility that there is indeed
gross mis-modeling of the cross sections. However, this does not change the
fact that significant cross section mis-modeling and therefore large
systematic errors would render the ND blind to new physics. Our findings
reveal the importance of properly accounting for the tuning procedure, as well
as properly estimating cross section uncertainties for BSM searches This will
be even more relevant for the future experiments HK and DUNE due to their high
statistics. We hope that this will motivate both theorists and
experimentalists to carefully consider the interplay between ND tuning and BSM
searches.

###### Acknowledgements.

We would like to thank André De Gouvêa, Andrew Furmanski, Joachim Kopp, Ivan
Martinez-Soler, Jonathan Paley, Afroditi Papadopoulou, Bryan Ramson, Noemi
Rocco, Alexandre Sousa, and Jeremy Wolcott for discussions, and especially
Kevin Kelly for providing the mono-neutrino UFO files. Fermilab is operated by
the Fermi Research Alliance, LLC under contract No. DE-AC02-07CH11359 with the
United States Department of Energy. This project has received support from the
European Union’s Horizon 2020 research and innovation programme under the
Marie Skłodowska-Curie grant agreement No 860881-HIDDeN. This research was
supported in part by the National Science Foundation under Grant No. NSF
PHY-1748958. The work of N.C. at Fermilab was supported by the U.S. Department
of Energy, Office of Science, Office of Workforce Development for Teachers and
Scientists, Office of Science Graduate Student Research (SCGSR) program, which
is administered by the Oak Ridge Institute for Science and Education (ORISE)
for the DOE. ORISE is managed by ORAU under contract number DE-SC0014664. The
work of N.C. at the University of Chicago has been supported by the DOE grant
DE-SC0013642.

##  Appendix A Tuning details

In this appendix, we provide further details on the tuning procedure and its
impact on the two observables in which we are interested: reconstructed
neutrino energy and missing transverse momentum. As discussed in Sec. 2, the
NOvA tuning is performed in the $(|\vec{q}\,|,\,q^{0})$ plane, only for MEC
events. While this tuning has many parameters (200 bins), it has two main
limitations. First, $|\vec{q}\,|$ and $q^{0}$ may not be fully correlated with
the relevant experimental observables in an analysis. Second, since the tuning
is applied only to the MEC component of the cross section, it does not have
full coverage of the $(|\vec{q}\,|,\,q^{0})$ plane. The latter point can be
appreciated in Fig. A1, where we show the distributions of $\nu_{\mu}$
charged-current events in NOvA ND for all events (left panel) and for MEC
events only (right panel), according to the GENIE event generator. The NOvA ND
flux was taken from Ref. NOvAflux .

In the left panel, the brightest strip is dominated by quasi-elastic
scattering events, while the lighter strip above has mostly resonance
production and some deep-inelastic scattering events. The top triangular area
is not kinematically accessible simply because charged current scattering on
nuclei proceeds via $t$-channel, and thus
$q^{2}=(q^{0})^{2}-|\vec{q}\,|^{2}<0$. Note also that the weights applied in
the tuning of the MEC cross section cannot be negative, further limiting the
coverage of the tuning procedure.

Figure A1: Distributions of all $\nu_{\mu}$-carbon charged current scattering
events (left) and MEC events (right) in NOvA ND in the true
$(|\vec{q}\,|,\,q^{0})$ plane. In the left panel, the lower strip is dominated
by quasi-elastic events, while the upper strip comes mainly from resonant and
deep inelastic events.

###  A.1 Sterile neutrinos

In Sec. 3.2 we have shown that the ND tuning essentially does not change the
shape of the reconstructed neutrino energy, and therefore it does not wash out
the spectral wiggles coming from short baseline oscillations. The main reason
for this is the lack of a clear correlation between the tuning parameters
$(|\vec{q}\,|,\,q^{0})$ and the neutrino energy. We can appreciate this in
Fig. A2, where we present the normalized distribution of events in the NOvA
near detector, restricting the true neutrino energy to two disparate regions,
$E_{\nu}\in[0.5,\,0.6]$ GeV (left panel) and $E_{\nu}\in[2.9,\,3.0]$ GeV
(right panel).

While there are visible differences in these distributions, the large overlap
between them shows the lack of correlation among tuning variables and neutrino
energy. The main goal of the experiment is to measure the true neutrino energy
spectra in both ND and FD so as to extract the oscillation probability.
Therefore, the tuning is performed in parameters that are not fully correlated
with neutrino energy. Regardless, while the lack of correlation prevents the
tuning of washing out ND oscillations, it also makes the tuning unable to fix
arbitrary cross section mis-modelings, which as discussed in Sec. 3.2 can bias
the experimental sensitivity to sterile neutrinos.

Figure A2: Normalized distribution of charged current $\nu_{\mu}$-carbon
events at NOvA for $E_{\nu}\in[0.5,\,0.6]$ GeV (left panel) and
$E_{\nu}\in[2.9,\,3]$ GeV (right panel). The large overlap for these two
disparate energies shows the lack of correlation between neutrino energy and
NOvA’s tuning variables.

###  A.2 Light neutrinophilic scalars

In Sec. 4, we have seen a large discrepancy between the missing transverse
momentum spectra generated by GENIE and NuWro, particularly after experimental
cuts. The goal of this section is to understand the origin of this
discrepancy. Let us first recall that we select events in antineutrino beam
mode with exactly one proton above 100 MeV, no neutrons above 100 MeV and no
pions. Naively, only a small fraction of MEC events, in which a low energy
neutron and a proton are emitted, contributes to this sample. Nevertheless,
final-state interactions will give rise to contributions to this sample coming
from other cross section processes. For example, although quasi-elastic
scattering of antineutrinos on nucleons produce a final-state neutron,
$\overline{\nu}_{\mu}+p\to\mu^{+}+n$, as this neutron propagates throughout
the nuclear medium it can kick out a proton, leading to an event that passes
experimental cuts. Even if at the NOvA energy scale these effects are
subleading, we call attention to the fact that cuts imposed here suppress
standard model backgrounds by two orders of magnitude, see Fig. 6. Therefore,
the experimental search is susceptible to percent-level differences between
generators.

In Fig. A3, we show the missing transverse momentum spectra broken down by
neutrino-nucleus cross section process for NuWro (left) and GENIE (right),
after experimental cuts. We can clearly see that the two generators strongly
disagree on the physics behind the events that pass selection: while in NuWro
the spectrum is dominated by quasi-elastic interactions (blue), in GENIE the
MEC events (red) are much more prominent. Since the magnitude of the quasi-
elastic component of the cross section in the two generators are quite
similar, this difference can only come from final-state interactions. This
already puts in question the focus on MEC-only tuning. More importantly, given
this wildly different physics in selected events between generators, there is
no reason to expect that an unphysical tuning would mitigate the observed
discrepancy.

Figure A3: Missing transverse momentum spectrum broken down by neutrino-
nucleus cross section process for NuWro (left) and GENIE (right), after
experimental cuts. The cuts select events with exactly one proton above 100
MeV, no neutrons above 100 MeV and no pions.

## References

  * (1) USQCD collaboration, “Lattice QCD and Neutrino-Nucleus Scattering”, _Eur. Phys. J. A_ 55 (2019) 196 [1904.09931].
  * (2) U. Mosel, “Neutrino event generators: foundation, status and future”, _J. Phys. G_ 46 (2019) 113001 [1904.11506].
  * (3) L. Alvarez-Ruso et al., “Snowmass 2021 LoI: Neutrino-induced Shallow- and Deep-Inelastic Scattering”, in _2022 Snowmass Summer Study_ , 9, 2020, 2009.04285.
  * (4) H. W. Bertini, “Low-Energy Intranuclear Cascade Calculation”, _Phys. Rev._ 131 (1963) 1801.
  * (5) J. Cugnon, “Monte Carlo calculation of high-energy heavy-ion interactions”, _Phys. Rev. C_ 22 (1980) 1885.
  * (6) G. F. Bertsch, H. Kruse and S. D. Gupta, “BOLTZMANN EQUATION FOR HEAVY ION COLLISIONS”, _Phys. Rev. C_ 29 (1984) 673.
  * (7) G. F. Bertsch and S. Das Gupta, “A Guide to microscopic models for intermediate-energy heavy ion collisions”, _Phys. Rept._ 160 (1988) 189.
  * (8) J. Cugnon, C. Volant and S. Vuillier, “Improved intranuclear cascade model for nucleon - nucleus interactions”, _Nucl. Phys. A_ 620 (1997) 475.
  * (9) A. Boudard, J. Cugnon, S. Leray and C. Volant, “Intranuclear cascade model for a comprehensive description of spallation reaction data”, _Phys. Rev. C_ 66 (2002) 044615.
  * (10) Y. Sawada, Y. Uozumi, S. Nogamine, T. Yamada, Y. Iwamoto, T. Sato et al., “Intranuclear cascade with emission of light fragment code implemented in the transport code system PHITS”, _Nucl. Instrum. Meth. B_ 291 (2012) 38.
  * (11) Y. Uozumi, T. Yamada, S. Nogamine and M. Nakano, “Intranuclear cascade model including collective excitations and trajectory deflections for (p, p-prime x) reactions around 50 MeV”, _Phys. Rev. C_ 86 (2012) 034610.
  * (12) T. Golan, C. Juszczak and J. T. Sobczyk, “Final State Interactions Effects in Neutrino-Nucleus Interactions”, _Phys. Rev. C_ 86 (2012) 015505 [1202.4197].
  * (13) N. Rocco, “Ab initio Calculations of Lepton-Nucleus Scattering”, _Front. in Phys._ 8 (2020) 116.
  * (14) J. Isaacson, W. I. Jay, A. Lovato, P. A. N. Machado and N. Rocco, “New approach to intranuclear cascades with quantum Monte Carlo configurations”, _Phys. Rev. C_ 103 (2021) 015502 [2007.15570].
  * (15) S. Dytman, Y. Hayato, R. Raboanary, J. T. Sobczyk, J. Tena Vidal and N. Vololoniaina, “Comparison of validation methods of simulations for final state interactions in hadron production experiments”, _Phys. Rev. D_ 104 (2021) 053006 [2103.07535].
  * (16) A. S. Meyer, M. Betancourt, R. Gran and R. J. Hill, “Deuterium target data for precision neutrino-nucleus cross sections”, _Phys. Rev. D_ 93 (2016) 113015 [1603.03048].
  * (17) MINERvA collaboration, “Measurement of $\bar{\nu}_{\mu}$ Charged-Current Single $\pi^{-}$ Production on Hydrocarbon in the Few-GeV Region using MINERvA”, _Phys. Rev. D_ 100 (2019) 052008 [1906.08300].
  * (18) A. M. Ankowski, A. Friedland, S. W. Li, O. Moreno, P. Schuster, N. Toro et al., “Lepton-Nucleus Cross Section Measurements for DUNE with the LDMX Detector”, _Phys. Rev. D_ 101 (2020) 053004 [1912.06140].
  * (19) MINERvA collaboration, “Probing nuclear effects with neutrino-induced charged-current neutral pion production”, _Phys. Rev. D_ 102 (2020) 072007 [2002.05812].
  * (20) MINERvA collaboration, “Double-differential inclusive charged-current $\nu_{\mu}$ cross sections on hydrocarbon in MINERvA at $\langle E_{\nu}\rangle\sim$ 3.5 GeV”, _Phys. Rev. D_ 101 (2020) 112007 [2002.12496].
  * (21) K. Borah, R. J. Hill, G. Lee and O. Tomalak, “Parametrization and applications of the low-$Q^{2}$ nucleon vector form factors”, _Phys. Rev. D_ 102 (2020) 074012 [2003.13640].
  * (22) T2K collaboration, “Simultaneous measurement of the muon neutrino charged-current cross section on oxygen and carbon without pions in the final state at T2K”, _Phys. Rev. D_ 101 (2020) 112004 [2004.05434].
  * (23) T2K collaboration, “First T2K measurement of transverse kinematic imbalance in the muon-neutrino charged-current single-$\pi^{+}$ production channel containing at least one proton”, _Phys. Rev. D_ 103 (2021) 112009 [2102.03346].
  * (24) O. Tomalak, Q. Chen, R. J. Hill and K. S. McFarland, “QED radiative corrections to neutrino-nucleon elastic scattering”, 2105.07939.
  * (25) MINERvA collaboration, “Measurement of inclusive charged-current $\nu_{\mu}$ cross sections as a function of muon kinematics at $<E_{\nu}>\sim 6~{}GeV$ on hydrocarbon”, _Phys. Rev. D_ 104 (2021) 092007 [2106.16210].
  * (26) NOvA collaboration, “Measurement of the Double-Differential Muon-neutrino Charged-Current Inclusive Cross Section in the NOvA Near Detector”, 2109.12220.
  * (27) A. M. Ankowski and A. Friedland, “Assessing the accuracy of the GENIE event generator with electron-scattering data”, _Phys. Rev. D_ 102 (2020) 053001 [2006.11944].
  * (28) electrons for neutrinos collaboration, “Inclusive Electron Scattering And The GENIE Neutrino Event Generator”, _Phys. Rev. D_ 103 (2021) 113003 [2009.07228].
  * (29) CLAS, e4v collaboration, “Electron-beam energy reconstruction for neutrino oscillation measurements”, _Nature_ 599 (2021) 565.
  * (30) MiniBooNE collaboration, “The Neutrino Flux prediction at MiniBooNE”, _Phys. Rev. D_ 79 (2009) 072002 [0806.1449].
  * (31) T2K collaboration, “Measurements of neutrino oscillation in appearance and disappearance channels by the T2K experiment with 6.6×1020 protons on target”, _Phys. Rev. D_ 91 (2015) 072010 [1502.01550].
  * (32) MINERvA collaboration, “Neutrino Flux Predictions for the NuMI Beam”, _Phys. Rev. D_ 94 (2016) 092005 [1607.00704].
  * (33) B. Batell, M. Pospelov and A. Ritz, “Exploring Portals to a Hidden Sector Through Fixed Targets”, _Phys. Rev. D_ 80 (2009) 095024 [0906.5614].
  * (34) P. deNiverville, M. Pospelov and A. Ritz, “Observing a light dark matter beam with neutrino experiments”, _Phys. Rev. D_ 84 (2011) 075020 [1107.4580].
  * (35) P. deNiverville, D. McKeen and A. Ritz, “Signatures of sub-GeV dark matter beams at neutrino experiments”, _Phys. Rev. D_ 86 (2012) 035022 [1205.3499].
  * (36) MiniBooNE collaboration, “Low Mass WIMP Searches with a Neutrino Experiment: A Proposal for Further MiniBooNE Running”, 1211.2258.
  * (37) E. Izaguirre, G. Krnjaic, P. Schuster and N. Toro, “New Electron Beam-Dump Experiments to Search for MeV to few-GeV Dark Matter”, _Phys. Rev. D_ 88 (2013) 114015 [1307.6554].
  * (38) B. Batell, P. deNiverville, D. McKeen, M. Pospelov and A. Ritz, “Leptophobic Dark Matter at Neutrino Factories”, _Phys. Rev. D_ 90 (2014) 115014 [1405.7049].
  * (39) B. A. Dobrescu and C. Frugiuele, “GeV-Scale Dark Matter: Production at the Main Injector”, _JHEP_ 02 (2015) 019 [1410.1566].
  * (40) P. Coloma, B. A. Dobrescu, C. Frugiuele and R. Harnik, “Dark matter beams at LBNF”, _JHEP_ 04 (2016) 047 [1512.03852].
  * (41) P. deNiverville, C.-Y. Chen, M. Pospelov and A. Ritz, “Light dark matter in neutrino beams: production modelling and scattering signatures at MiniBooNE, T2K and SHiP”, _Phys. Rev. D_ 95 (2017) 035006 [1609.01770].
  * (42) G. Magill and R. Plestid, “Neutrino Trident Production at the Intensity Frontier”, _Phys. Rev. D_ 95 (2017) 073004 [1612.05642].
  * (43) MINOS+ collaboration, “Search for sterile neutrinos in MINOS and MINOS+ using a two-detector fit”, _Phys. Rev. Lett._ 122 (2019) 091803 [1710.06488].
  * (44) E. Bertuzzo, S. Jana, P. A. N. Machado and R. Zukanovich Funchal, “Dark Neutrino Portal to Explain MiniBooNE excess”, _Phys. Rev. Lett._ 121 (2018) 241801 [1807.09877].
  * (45) A. Falkowski, G. Grilli di Cortona and Z. Tabrizi, “Future DUNE constraints on EFT”, _JHEP_ 04 (2018) 101 [1802.08296].
  * (46) G. Magill, R. Plestid, M. Pospelov and Y.-D. Tsai, “Millicharged particles in neutrino experiments”, _Phys. Rev. Lett._ 122 (2019) 071801 [1806.03310].
  * (47) P. Ballett, S. Pascoli and M. Ross-Lonergan, “U(1)’ mediated decays of heavy sterile neutrinos in MiniBooNE”, _Phys. Rev. D_ 99 (2019) 071701 [1808.02915].
  * (48) A. de Gouvêa, P. J. Fox, R. Harnik, K. J. Kelly and Y. Zhang, “Dark Tridents at Off-Axis Liquid Argon Neutrino Detectors”, _JHEP_ 01 (2019) 001 [1809.06388].
  * (49) P. Ballett, M. Hostert, S. Pascoli, Y. F. Perez-Gonzalez, Z. Tabrizi and R. Zukanovich Funchal, “Neutrino Trident Scattering at Near Detectors”, _JHEP_ 01 (2019) 119 [1807.10973].
  * (50) P. deNiverville and C. Frugiuele, “Hunting sub-GeV dark matter with the NO$\nu$A near detector”, _Phys. Rev. D_ 99 (2019) 051701 [1807.06501].
  * (51) J. R. Jordan, Y. Kahn, G. Krnjaic, M. Moschella and J. Spitz, “Signatures of Pseudo-Dirac Dark Matter at High-Intensity Neutrino Experiments”, _Phys. Rev. D_ 98 (2018) 075020 [1806.05185].
  * (52) J. R. Jordan, Y. Kahn, G. Krnjaic, M. Moschella and J. Spitz, “Severe Constraints on New Physics Explanations of the MiniBooNE Excess”, _Phys. Rev. Lett._ 122 (2019) 081801 [1810.07185].
  * (53) P. Ballett, M. Hostert, S. Pascoli, Y. F. Perez-Gonzalez, Z. Tabrizi and R. Zukanovich Funchal, “$Z^{\prime}$s in neutrino scattering at DUNE”, _Phys. Rev. D_ 100 (2019) 055012 [1902.08579].
  * (54) R. Harnik, Z. Liu and O. Palamara, “Millicharged Particles in Liquid Argon Neutrino Experiments”, _JHEP_ 07 (2019) 170 [1902.03246].
  * (55) ArgoNeuT collaboration, “Improved Limits on Millicharged Particles Using the ArgoNeuT Experiment at Fermilab”, _Phys. Rev. Lett._ 124 (2020) 131801 [1911.07996].
  * (56) W. Altmannshofer, S. Gori, J. Martín-Albo, A. Sousa and M. Wallbank, “Neutrino Tridents at DUNE”, _Phys. Rev. D_ 100 (2019) 115029 [1902.06765].
  * (57) Y.-D. Tsai, P. deNiverville and M. X. Liu, “Dark Photon and Muon $g-2$ Inspired Inelastic Dark Matter Models at the High-Energy Intensity Frontier”, _Phys. Rev. Lett._ 126 (2021) 181801 [1908.07525].
  * (58) B. Batell, J. Berger and A. Ismail, “Probing the Higgs Portal at the Fermilab Short-Baseline Neutrino Experiments”, _Phys. Rev. D_ 100 (2019) 115039 [1909.11670].
  * (59) C. A. Argüelles et al., “New opportunities at the next-generation neutrino experiments I: BSM neutrino physics and dark matter”, _Rept. Prog. Phys._ 83 (2020) 124201 [1907.08311].
  * (60) V. De Romeri, K. J. Kelly and P. A. N. Machado, “DUNE-PRISM Sensitivity to Light Dark Matter”, _Phys. Rev. D_ 100 (2019) 095010 [1903.10505].
  * (61) A. de Gouvea, P. A. N. Machado, Y. F. Perez-Gonzalez and Z. Tabrizi, “Measuring the Weak Mixing Angle in the DUNE Near Detector Complex”, _Phys. Rev. Lett._ 125 (2020) 051803 [1912.06658].
  * (62) J. M. Berryman, A. de Gouvea, P. J. Fox, B. J. Kayser, K. J. Kelly and J. L. Raaf, “Searches for Decays of New Particles in the DUNE Multi-Purpose Near Detector”, _JHEP_ 02 (2020) 174 [1912.07622].
  * (63) DUNE collaboration, “Prospects for beyond the Standard Model physics searches at the Deep Underground Neutrino Experiment”, _Eur. Phys. J. C_ 81 (2021) 322 [2008.12769].
  * (64) NOvA, R. Group collaboration, “Adjusting neutrino interaction models and evaluating uncertainties using NOvA near detector data”, _Eur. Phys. J. C_ 80 (2020) 1119 [2006.08727].
  * (65) A. Friedland and S. W. Li, “Understanding the energy resolution of liquid argon neutrino detectors”, _Phys. Rev. D_ 99 (2019) 036009 [1811.06159].
  * (66) CAPTAIN collaboration, “First Measurement of the Total Neutron Cross Section on Argon Between 100 and 800 MeV”, _Phys. Rev. Lett._ 123 (2019) 042502 [1903.05276].
  * (67) A. Friedland and S. W. Li, “Simulating hadron test beams in liquid argon”, _Phys. Rev. D_ 102 (2020) 096005 [2007.13336].
  * (68) C. Andreopoulos et al., “The GENIE Neutrino Monte Carlo Generator”, _Nucl. Instrum. Meth. A_ 614 (2010) 87 [0905.2517].
  * (69) T. Golan, J. T. Sobczyk and J. Zmuda, “NuWro: the Wroclaw Monte Carlo Generator of Neutrino Interactions”, _Nucl. Phys. B Proc. Suppl._ 229-232 (2012) 499.
  * (70) O. Buss, T. Gaitanos, K. Gallmeister, H. van Hees, M. Kaskulov, O. Lalakulich et al., “Transport-theoretical Description of Nuclear Reactions”, _Phys. Rept._ 512 (2012) 1 [1106.1344].
  * (71) J. Isaacson, W. I. Jay, A. Lovato, P. A. N. Machado and N. Rocco, “ACHILLES: A novel event generator for electron- and neutrino-nucleus scattering”, 2205.06378.
  * (72) NOvA collaboration, “Improved measurement of neutrino oscillation parameters by the NOvA experiment”, _Phys. Rev. D_ 106 (2022) 032004 [2108.08219].
  * (73) T2K collaboration, “Improved constraints on neutrino mixing from the T2K experiment with $\mathbf{3.13\times 10^{21}}$ protons on target”, _Phys. Rev. D_ 103 (2021) 112008 [2101.03779].
  * (74) NOvA collaboration, “New constraints on oscillation parameters from $\nu_{e}$ appearance and $\nu_{\mu}$ disappearance in the NOvA experiment”, _Phys. Rev. D_ 98 (2018) 032012 [1806.00096].
  * (75) R. Gran, “Model Uncertainties for Valencia RPA Effect for MINERvA”, 1705.02932.
  * (76) P. Rodrigues, C. Wilkinson and K. McFarland, “Constraining the GENIE model of neutrino-induced single pion production using reanalyzed bubble chamber data”, _Eur. Phys. J. C_ 76 (2016) 474 [1601.01888].
  * (77) MiniBooNE collaboration, “Measurement of Neutrino-Induced Charged-Current Charged Pion Production Cross Sections on Mineral Oil at E${}_{\nu}\sim 1~{}\textrm{GeV}$”, _Phys. Rev. D_ 83 (2011) 052007 [1011.3572].
  * (78) MINOS collaboration, “Study of quasielastic scattering using charged-current $\nu$$\mu$ -iron interactions in the MINOS near detector”, _Phys. Rev. D_ 91 (2015) 012005 [1410.8613].
  * (79) MINERvA collaboration, “Cross sections for $\nu_{\mu}$ and $\bar{\nu}_{\mu}$ induced pion production on hydrocarbon in the few-GeV region using MINERvA”, _Phys. Rev. D_ 94 (2016) 052005 [1606.07127].
  * (80) MINERvA collaboration, “Measurement of $\nu_{\mu}$ charged-current single $\pi^{0}$ production on hydrocarbon in the few-GeV region using MINERvA”, _Phys. Rev. D_ 96 (2017) 072003 [1708.03723].
  * (81) T2K collaboration, “Measurement of the muon neutrino charged-current single $\pi^{+}$ production on hydrocarbon using the T2K off-axis near detector ND280”, _Phys. Rev. D_ 101 (2020) 012007 [1909.03936].
  * (82) K. Kainulainen, J. Maalampi and J. T. Peltoniemi, “Inert neutrinos in supernovae”, _Nucl. Phys. B_ 358 (1991) 435.
  * (83) K. J. Kelly and Y. Zhang, “Mononeutrino at DUNE: New Signals from Neutrinophilic Thermal Dark Matter”, _Phys. Rev. D_ 99 (2019) 055034 [1901.01259].
  * (84) LSND collaboration, “Evidence for neutrino oscillations from the observation of $\bar{\nu}_{e}$ appearance in a $\bar{\nu}_{\mu}$ beam”, _Phys. Rev. D_ 64 (2001) 112007 [hep-ex/0104049].
  * (85) MiniBooNE collaboration, “Updated MiniBooNE neutrino oscillation results with increased data and new background studies”, _Phys. Rev. D_ 103 (2021) 052002 [2006.16883].
  * (86) T. A. Mueller et al., “Improved Predictions of Reactor Antineutrino Spectra”, _Phys. Rev. C_ 83 (2011) 054615 [1101.2663].
  * (87) G. Mention, M. Fechner, T. Lasserre, T. A. Mueller, D. Lhuillier, M. Cribier et al., “The Reactor Antineutrino Anomaly”, _Phys. Rev. D_ 83 (2011) 073006 [1101.2755].
  * (88) P. Huber, “On the determination of anti-neutrino spectra from nuclear reactors”, _Phys. Rev. C_ 84 (2011) 024617 [1106.0687].
  * (89) GALLEX collaboration, “Final results of the Cr-51 neutrino source experiments in GALLEX”, _Phys. Lett. B_ 420 (1998) 114.
  * (90) SAGE collaboration, “Measurement of the response of the Russian-American gallium experiment to neutrinos from a Cr-51 source”, _Phys. Rev. C_ 59 (1999) 2246 [hep-ph/9803418].
  * (91) J. N. Abdurashitov et al., “Measurement of the response of a Ga solar neutrino experiment to neutrinos from an Ar-37 source”, _Phys. Rev. C_ 73 (2006) 045805 [nucl-ex/0512041].
  * (92) F. Kaether, W. Hampel, G. Heusser, J. Kiko and T. Kirsten, “Reanalysis of the GALLEX solar neutrino flux and source experiments”, _Phys. Lett. B_ 685 (2010) 47 [1001.2731].
  * (93) C. Giunti and M. Laveder, “Statistical Significance of the Gallium Anomaly”, _Phys. Rev. C_ 83 (2011) 065504 [1006.3244].
  * (94) V. V. Barinov et al., “Results from the Baksan Experiment on Sterile Transitions (BEST)”, _Phys. Rev. Lett._ 128 (2022) 232501 [2109.11482].
  * (95) ALEPH collaboration, “A Direct measurement of the invisible width of the Z from single photon counting”, _Phys. Lett. B_ 313 (1993) 520.
  * (96) OPAL collaboration, “Measurement of single photon production in e+ e- collisions near the Z0 resonance”, _Z. Phys. C_ 65 (1995) 47.
  * (97) L3 collaboration, “Determination of the number of light neutrino species from single photon production at LEP”, _Phys. Lett. B_ 431 (1998) 199.
  * (98) DELPHI collaboration, “Cross-sections and leptonic forward backward asymmetries from the Z0 running of LEP”, _Eur. Phys. J. C_ 16 (2000) 371.
  * (99) M. Dentler, A. Hernández-Cabezudo, J. Kopp, P. A. N. Machado, M. Maltoni, I. Martinez-Soler et al., “Updated Global Analysis of Neutrino Oscillations in the Presence of eV-Scale Sterile Neutrinos”, _JHEP_ 08 (2018) 010 [1803.10661].
  * (100) M. H. Moulai, C. A. Argüelles, G. H. Collin, J. M. Conrad, A. Diaz and M. H. Shaevitz, “Combining Sterile Neutrino Fits to Short Baseline Data with IceCube Data”, _Phys. Rev. D_ 101 (2020) 055020 [1910.13456].
  * (101) C. Giunti and T. Lasserre, “eV-scale Sterile Neutrinos”, _Ann. Rev. Nucl. Part. Sci._ 69 (2019) 163 [1901.08330].
  * (102) MINOS+, Daya Bay collaboration, “Improved Constraints on Sterile Neutrino Mixing from Disappearance Searches in the MINOS, MINOS+, Daya Bay, and Bugey-3 Experiments”, _Phys. Rev. Lett._ 125 (2020) 071801 [2002.00301].
  * (103) IceCube collaboration, “Searches for Sterile Neutrinos with the IceCube Detector”, _Phys. Rev. Lett._ 117 (2016) 071801 [1605.01990].
  * (104) IceCube collaboration, “Searching for eV-scale sterile neutrinos with eight years of atmospheric neutrinos at the IceCube Neutrino Telescope”, _Phys. Rev. D_ 102 (2020) 052009 [2005.12943].
  * (105) SNO collaboration, “An Independent Measurement of the Total Active B-8 Solar Neutrino Flux Using an Array of He-3 Proportional Counters at the Sudbury Neutrino Observatory”, _Phys. Rev. Lett._ 101 (2008) 111301 [0806.0989].
  * (106) Super-Kamiokande collaboration, “Solar Neutrino Measurements in Super-Kamiokande-IV”, _Phys. Rev. D_ 94 (2016) 052010 [1606.07538].
  * (107) BOREXINO collaboration, “Comprehensive measurement of $pp$-chain solar neutrinos”, _Nature_ 562 (2018) 505.
  * (108) RENO, NEOS collaboration, “Search for sterile neutrino oscillations using RENO and NEOS data”, _Phys. Rev. D_ 105 (2022) L111101 [2011.00896].
  * (109) DANSS collaboration, “Search for sterile neutrinos at the DANSS experiment”, _Phys. Lett. B_ 787 (2018) 56 [1804.04046].
  * (110) NOvA Collaboration, “Flux Histograms for 2017 Analysis.” http://publicdocs.fnal.gov/cgi-bin/ShowDocument?docid=8.
  * (111) Keloth, Rijeesh, “NOvA Short-Baseline Tau-Neutrino Appearance Search.” https://absuploads.aps.org/presentation.cfm?pid=12831.
  * (112) T. Katori, “Meson Exchange Current (MEC) Models in Neutrino Interaction Generators”, _AIP Conf. Proc._ 1663 (2015) 030001 [1304.6014].
  * (113) J. Nieves, I. Ruiz Simo and M. J. Vicente Vacas, “Inclusive Charged–Current Neutrino–Nucleus Reactions”, _Phys. Rev. C_ 83 (2011) 045501 [1102.2777].
  * (114) J. M. Berryman, A. De Gouvêa, K. J. Kelly and Y. Zhang, “Lepton-Number-Charged Scalars and Neutrino Beamstrahlung”, _Phys. Rev. D_ 97 (2018) 075030 [1802.00009].
  * (115) D. Abercrombie et al., “Dark Matter benchmark models for early LHC Run-2 Searches: Report of the ATLAS/CMS Dark Matter Forum”, _Phys. Dark Univ._ 27 (2020) 100371 [1507.00966].
  * (116) J. Alwall, M. Herquet, F. Maltoni, O. Mattelaer and T. Stelzer, “MadGraph 5 : Going Beyond”, _JHEP_ 06 (2011) 128 [1106.0522].


# Baryogenesis in $f(R,L_{m})$ gravity

Lakhan V. Jaybhaye 0000-0003-1497-276X<EMAIL_ADDRESS>Department of
Mathematics, Birla Institute of Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India. Snehasish Bhattacharjee
0000-0002-7350-7043<EMAIL_ADDRESS>Institute of Astronomy,
National Central University, 32001 Taoyuan, Taiwan Department of Particle
Physics and Astrophysics, Weizmann Institute of Science, 76100 Rehovot, Israel
P.K. Sahoo 0000-0003-2130-8832<EMAIL_ADDRESS>Department
of Mathematics, Birla Institute of Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India.

###### Abstract

This paper aims to recreate the gravitational baryogenesis epoch in the
framework of the $f(R,L_{m})$ theory of gravity, where $R$ and $L_{m}$ are the
curvature scalar and the matter Langragian, respectively. In particular, we
examine the model, $f(R,L_{m})=\frac{R}{2}+L_{m}^{\alpha}+\zeta$, under the
supposition that the universe is saturated with dark energy and perfect fluid,
with a non-zero baryon to entropy ratio during a radiation dominance era. We
confine the model with the gravitational baryogenesis scenario, emphasizing
the appropriate values of model parameters compatible with the baryon-to-
entropy ratio observation data. Our study found that $f(R,L_{m})$ gravity can
considerably and steadily make a contribution to the phenomenon of
gravitational baryogenesis.

##  I  Introduction

Even before cosmology became an independent research division, one of the
unanswered questions was the abundance of matter over antimatter in our
Universe. Significant observational evidence, like Big-Bang Nucleosynthesis
(BBN) [1] and the Cosmic Microwave Background (CMB) [2], has strongly
indicated that matter predominates over antimatter in the universe. This
superiority is known as baryogenesis. Recent research suggests that the
asymmetrical correlation between matter and antimatter originally came at the
beginning of the cosmos. However, the true origin of baryon asymmetry (BA)
remains a mystery that requires further investigation.
Numerous theories have emerged to solve this mystery of BA by evaluating
interactions in the primitive universe that go beyond the standard model, some
of which are Affleck-Dine baryogenesis [3, 4, 5], spontaneous baryogenesis [6,
7, 8] electroweak baryogenesis [9, 10], grand unified theories (GUTs)[11],
baryogenesis of thermal and black hole evaporation [12], these multiple
baryogenesis contexts discuss how this universe could have more matter than
antimatter during the matter or radiation epoch. The gravitational
baryogenesis process employs one of the Sakharov criteria [13] proposed that
the baryon asymmetry can be produced by three necessary conditions: (1)
processes that violate baryon number, (2) Charge (C) and ChargeParity (CP)
violation, and (3) Out of thermal equilibrium interactions. The key component
is a CP-violating interaction stipulated by coupling between the baryon matter
current $J^{\mu}$ and the derivative of the Ricci scalar curvature $R$, in the
form

$\frac{1}{M_{*}^{2}}\int{\sqrt{-g}J^{\mu}\partial_{\mu}(R)d^{4}x}$ (1)

In Eq. (1), the parameter $M_{*}$ represents the cutoff scale of the effective
theory, whereas $R$, $J^{\mu}$, and $g$ stand for the Ricci scalar, baryonic
matter current, and the trace of the metric tensor respectively. The baryon to
entropy ratio $\frac{n_{B}}{s}\propto\dot{R}$, in the case of flat FRW
Universe, where the dot represents cosmic time derivative. In the situation of
a radiation-dominated era with $\omega=\frac{1}{3}$, the net baryon asymmetry
generated by Eq.(1) is zero. In the past few years, a number of authors have
studied the mysterious concept of Baryogenesis within the context of modified
gravity. In [14, 15], $f(R)$ gravity theories are addressed in terms of
gravitational baryogenesis, whereas in [16] Gauss-Bonnet gravity, in [17]
Gauss-Bonnet braneworld cosmology, in [18] $f(T)$ gravity, in [19] $f(P)$
gravity, teleparallel gravity [26], in [20, 21, 22] $f(R,T)$ gravity, in [23]
$f(R,T,X)$ gravity, in [24] $f(T,B)$ gravity, and in [25] $f(Q,T)$ gravity,
etc. The objective of this work is to look into the framework of gravitational
baryogenesis in the $f(R,L_{m})$ gravity theory, which is an extension of the
$f(R)$ gravity theory [27], that includes an explicit coupling of Rcicci
scalar $R$ and matter Lagrangian $L_{m}$. This Theory was proposed by Harko
and Lobo [28] as well as several works with interesting results have been
found [29, 30, 31, 32, 33]. The following is how the present manuscript is
organized: Section II includes an overview of $f(R,L_{m})$ gravity. In Section
III, we will provide some essential baryogenesis factors before investigating
gravitational baryogenesis in $f(R,L_{m})$ gravity and then we explain
baryogenesis in $f(R,L_{m})$ gravity in detail and imply generating
observationally acceptable baryon to entropy ratios for $f(R,L_{m})$ gravity
model and we will also go over the generalized form of baryogenesis that
applies to assuming gravity model. In the final section IV, we will discuss
the conclusions of the present work.

##  II Overview of $f(R,L_{m})$ Gravity

Consideration of the overall action in modified $f(R,L_{m})$ gravity [28]

$S=\int{f(R,L_{m})\sqrt{-g}d^{4}x}$ (2)

Where $R$, $L_{m}$ stand for the Ricci scalar curvature and the matter
Lagrangian. By varying the action (2) with respect to the metric tensor
$g_{\mu\nu}$, we were able to derive the field equation governing the dynamics
of gravitational interactions.

$f_{R}R_{\mu\nu}+(g_{\mu\nu}\square-\nabla_{\mu}\nabla_{\nu})f_{R}-\frac{1}{2}(f-f_{L_{m}}L_{m})g_{\mu\nu}=\frac{1}{2}f_{L_{m}}T_{\mu\nu}$
(3)

Here $f_{R}\equiv\frac{\partial f}{\partial R}$,
$f_{L_{m}}\equiv\frac{\partial f}{\partial L_{m}}$, and $T_{\mu\nu}$
represents the stress-energy tensor for the ideal fluid, which is described by

$T_{\mu\nu}=\frac{-2}{\sqrt{-g}}\frac{\delta(\sqrt{-g}L_{m})}{\delta
g^{\mu\nu}}$ (4)

It further takes the covariant derivative in equation (3) to yield the
following result:

$\nabla^{\mu}T_{\mu\nu}=2\nabla^{\mu}ln(f_{L_{m}})\frac{\partial
L_{m}}{\partial g^{\mu\nu}}$ (5)

Now, we think of a spatially flat FLRW metric as

$ds^{2}=-dt^{2}+a^{2}(t)[dx^{2}+dy^{2}+dz^{2}]$ (6)

where, $a(t)$ is the cosmic scale factor. So the Ricci scalar produced
pertaining to the line element (6) is

$R=6\frac{\ddot{a}}{a}+6\bigl{(}\frac{\dot{a}}{a}\bigr{)}^{2}=6(\dot{H}+2H^{2})$
(7)

Overall, we get the modified Friedman equations that characterize the dynamics
of the universe in $f(R,L_{m})$ gravity by employing (6) in (3) [30],

$3H^{2}f_{R}+\frac{1}{2}\left(f-f_{R}R-f_{L_{m}}L_{m}\right)+3H\dot{f_{R}}=\frac{1}{2}f_{L_{m}}\rho$
(8)

and

$\dot{H}f_{R}+3H^{2}f_{R}-\ddot{f_{R}}-3H\dot{f_{R}}+\frac{1}{2}\left(f_{L_{m}}L_{m}-f\right)=\frac{1}{2}f_{L_{m}}p$
(9)

where $H=\frac{\dot{a}}{a}$, the Hubble parameter, and the dot stand for the
derivative with respect to cosmic time $t$. The matter density and pressure
are denoted by $\rho$ and $p$ in the equations above.

##  III  Baryogenesis in $f(R,L_{m})$

We shall demonstrate how $f(R,L_{m})$ gravity addresses the gravitational
baryogenesis difficulty in the cosmos in this part. A crucial parameter to
figure out asymmetry is known as BAF given by

$\eta_{B}=\frac{n_{B}-\bar{n}_{B}}{s}$ (10)

where $n_{B}$ is the baryon number and $\bar{n}_{B}$ is the antibaryon number,
and s is the entropy of the universe. Observational evidence like the BBN [1]
and CMB [2] confirm that the restriction on this BAF is $\frac{n_{B}}{s}\simeq
9\times 10^{-11}$. For the $f(R,L_{m})$ gravity, we take into consideration an
interaction term that violates CP and is produced by the baryon asymmetry of
the universe of the form

$\frac{1}{M_{*}^{2}}\int{\sqrt{-g}J^{\mu}\partial_{\mu}(R+L_{m})d^{4}x}$ (11)

As a result, for an interaction of (11) that violates CP, the resulting baryon
to entropy ratio in $f(R,L_{m})$ gravity is as follows:

$\frac{n_{B}}{s}\simeq-\frac{15g_{B}}{4\pi^{2}g_{*s}}\frac{(\dot{R}+\dot{L_{m}})}{M^{2}_{*}T_{D}}$
(12)

In (12), $g_{B}$ is the total number of intrinsic degrees of freedom of
baryons, $g_{*}s$ is the total number of degrees of freedom of massless
particles, and the critical temperature is $T_{D}$ is the temperature of the
universe when all interactions that cause baryon asymmetry to begin. We will
assume that a thermal equilibrium exists, with energy density being
proportional to temperature $T$ as

$\rho(T)=\frac{\pi^{2}}{30}g_{*s}T^{4}$ (13)

Inside the framework of Einstein’s general theory of relativity, suppose the
universe’s matter content is a perfect fluid with the constant equation of
state parameter $\omega=\frac{p}{\rho}$ and the Ricci scalar $R$, as

$R=-8\pi G(1-3\omega)\rho$ (14)

In General Relativity (GR), if the universe is filled with radiation, the
baryon number to entropy ratio equals zero. For the other content of the
matter, this result differs from zero. However, in $f(R,L_{m})$ gravity
theories, a net baryon asymmetry may be generated during the radiation-
dominated era. To accomplish this, we focus on a specific $f(R,L_{m})$ model
to describe how we can recover the baryogenesis epoch with this model. We
calculate the baryon-to-entropy ratio for the model by imagining a universe
full of dark energy and perfect fluid with the constant equation of state
parameter $\omega=\frac{P}{\rho}$.

### 3.1 The perfect fluid with $f(R,L_{m})$ gravity

The recently proposed model we are using in this section was presented in Ref.
[30]. In our work, we borrowed this model and its correlating field equations
from this paper and utilized them to explain baryogenesis. That is

$f(R,L_{m})=\frac{R}{2}+L_{m}^{\alpha}+\zeta$ (15)

where $\alpha$ and $\zeta$ are model parameters. Then, for this specific
$f(R,L_{m})$ model with $L_{m}=\rho$ [29], the Friedmann equations (8) and (9)
are transformed into

$3H^{2}=(2\alpha-1)\rho^{\alpha}-\zeta$ (16)

and

$-2\dot{H}-3H^{2}=(\alpha p+(1-\alpha)\rho)\rho^{\alpha-1}+\zeta$ (17)

One can obtain the following energy balance equation by taking the trace of
the field equations

$\alpha\dot{\rho}+3H\rho=0$ (18)

To continue, we will presume that the scale factor adapts as a power-law
$a(t)=Bt^{\beta}$, where $\beta=\frac{2}{3(1+\omega)}$, and $B$ is a constant
parameter [16]. Then the Hubble parameter $H(t)$ and the energy density
$\rho(t)$ expressions for this model are as follows

$H(t)=\frac{\beta}{t}$ (19)

$\rho(t)=\big{[}\frac{3\beta^{2}+\zeta
t^{2}}{(2\alpha-1)t^{2}}\big{]}^{\frac{1}{\alpha}}$ (20)

Equating Eqs. (20) and (13), we obtain the decoupling time $t_{D}$ as a
function of the decoupling temperature $T_{D}$ written as

$t_{D}=\big{[}\frac{3\beta^{2}}{(2\alpha-1)(\frac{\pi^{2}}{30}g_{*s}T_{D}^{4})^{\alpha}-\zeta}\big{]}^{\frac{1}{2}}$
(21)

Using (21), we retrieve a final expression for the baryon-to-entropy ratio for
the current $f(R,L_{m})$ model

$\frac{n_{B}}{s}\simeq\left(\frac{-15g_{B}}{4\pi^{2}g_{*s}M^{2}_{*}T_{D}}\right)\left(\frac{\left[(2\alpha-1)(\frac{\pi^{2}}{30}g_{*s}T_{D}^{4})^{\alpha}-\zeta\right]^{\frac{3}{2}}}{\sqrt{27}\beta^{3}}\right)\left(12(\beta-2\beta^{2})-\frac{6\beta\left(\frac{\pi^{2}}{30}g_{*s}T_{D}^{4}\right)^{1-\alpha}}{\alpha(2\alpha-1)}\right)$
(22)

In the radiation-dominated phase, $\beta=\frac{1}{2}$. Hence, (22) reduces to

$\frac{n_{B}}{s}\simeq\left(\frac{-15g_{B}}{4\pi^{2}g_{*s}M^{2}_{*}T_{D}}\right)\left(\frac{8\left[(2\alpha-1)(\frac{\pi^{2}}{30}g_{*s}T_{D}^{4})^{\alpha}-\zeta\right]^{\frac{3}{2}}}{\sqrt{27}}\right)\left(-\frac{3\left(\frac{\pi^{2}}{30}g_{*s}T_{D}^{4}\right)^{1-\alpha}}{\alpha(2\alpha-1)}\right)$
(23)

As shown in Eq.(23), the resulting baryon to entropy ratio is non-zero. The
ratio (23) can be adjusted to satisfy the observational constraints depending
on the matter content, but the most interesting feature of a perfect fluid
dominated Universe for the case of $f(R,L_{m})$ gravity baryogenesis is that
the ratio is non-zero in the radiation dominated case.

Figure 1: The baryon to entropy ratio for the model
$f(R,L_{m})=\frac{R}{2}+L_{m}^{\alpha}+\zeta$. The graphs are plotted for
$\alpha$, for the varying values of $\beta$, $\beta=0.55$ (Red), $\beta=0.60$
(Blue), $\beta=0.65$ (Green), $\beta=0.70$ (Orange), $\beta=0.75$ (Purple),
and $\zeta=2$. The dashed line represents the observational value.

Substituting $g_{*s}=106$, $g_{B}=1$, $T_{D}=2\times 10^{12}GeV$ and
$M_{*}=2\times 10^{16}GeV$ [34], with model parameters $\alpha=0.79$ and
$\zeta=2$ in Eq. (23) the resultant baryon to entropy ratio reads
$\frac{n_{B}}{s}\simeq 7.28749\times 10^{-11}$ which is in excellent agreement
with observations. The intersections of the curves in Fig 1 with the curve
that reveals the observational value (dashed line) are in fine contract value
of baryon to entropy ratio for particular values of $\alpha$ about including
$0.75$ and $0.85$. Interestingly, when the parameter $0.5\leq\alpha\leq 0.69$,
thus every curve tends to zero, which is consistent with theoretical results.

### 3.2 Generalized Gravitational Baryogenesis

We will now attempt to investigate the effects of a more comprehensive and
generalized CP-violating interaction proportional to
$\partial_{\mu}(f(R,L_{m}))$ rather than $\partial_{\mu}(R+L_{m})$ in
attempting to address the baryon asymmetry of the Universe for the chosen
$f(R,L_{m})=\frac{R}{2}+L_{m}^{\alpha}+\zeta$ model. In $f(R,L_{m})$ gravity,
we can express the generalized CP-violating interaction as [20]

$\frac{1}{M_{*}^{2}}\int{\sqrt{-g}J^{\mu}\partial_{\mu}(f(R,L_{m}))d^{4}x}$
(24)

The resulting baryon-to-entropy ratio for (24) is

$\frac{n_{B}}{s}\simeq\frac{-15g_{B}(\dot{R}f_{R}+\dot{L_{m}}f_{L_{m}})}{4g_{*s}M^{2}_{*}\pi^{2}T_{D}}$
(25)

Substituting (14), (15), (19), and (21) in (25), we get the baryon to entropy
ratio as

$\frac{n_{B}}{s}\simeq\left(\frac{-15g_{B}}{4\pi^{2}g_{*s}M^{2}_{*}T_{D}}\right)\left(\frac{\left[(2\alpha-1)(\frac{\pi^{2}}{30}g_{*s}T_{D}^{4})^{\alpha}-\zeta\right]^{\frac{3}{2}}}{\sqrt{27}\beta^{3}}\right)\left(6(\beta-2\beta^{2})-\frac{6\beta}{(2\alpha-1)}\right)$
(26)

We finally obtain, as discussed in the previous section for a radiation-
dominated universe by setting $\beta=\frac{1}{2}$, in Eq. (26)

$\frac{n_{B}}{s}\simeq\left(\frac{-15g_{B}}{4\pi^{2}g_{*s}M^{2}_{*}T_{D}}\right)\left(\frac{8\left[(2\alpha-1)(\frac{\pi^{2}}{30}g_{*s}T_{D}^{4})^{\alpha}-\zeta\right]^{\frac{3}{2}}}{\sqrt{27}}\right)\left(-\frac{3}{4(2\alpha-1)}\right)$
(27)

By substituting $M_{*}$, $g_{*s}$, $g_{B}$, $T_{D}$ as before, $\alpha=0.93$,
and $\zeta=2$, the obtained baryon to entropy ratio is $\frac{n_{B}}{s}\simeq
7.01\times 10^{-11}$, which is also in excellent agreement with observations.

Figure 2: The baryon to entropy ratio for the model
$f(R,L_{m})=\frac{R}{2}+L_{m}^{\alpha}+\zeta$. The graphs are plotted for
$\alpha$, for the varying values of $\beta$, $\beta=0.018$ (Red),
$\beta=0.026$ (Blue), $\beta=0.036$ (Green), $\beta=0.046$ (Orange),
$\beta=0.056$ (Purple), and $\zeta=2$. The dashed line represents the
observational value.

We show $\frac{n_{B}}{s}$ for the generalised baryogenesis interaction as a
function of $\alpha$ in Fig 2. The intersections of the curves in Fig 2 with
the curve that reveals the observational value (dashed line) are in fine
contract value of baryon to entropy ratio for particular values of $\alpha$
about including $0.79$ and $0.86$. Interestingly, when the parameter
$0.6\leq\alpha\leq 0.79$, thus every curve tends to zero, which is consistent
with theoretical results. Acceptable baryon to entropy values was obtained by
$f(R,L_{m})=\frac{R}{2}+L_{m}^{\alpha}+\zeta$, which resulted in physically
acceptable baryon to entropy ratios. Thus, in $f(R,L_{m})$ gravity, the
problem of baryogenesis can be rectified.

##  IV  Conclusions

The paper investigates the methodology of gravitational baryogenesis from the
perspective of $f(R,L_{m})$ gravity theory. We evaluate the baryon to entropy
ratio for specific $f(R,L_{m})=\frac{R}{2}+L_{m}^{\alpha}+\zeta$ model based
on the CP-violating interaction that will produce the Universe’s baryon
asymmetry and assessing the Universe’s matter content as a perfect fluid with
a constant equation of state parameter $\omega$. In comparison to GR, we
demonstrate that its baryon-to-entropy ratio in a radiation-dominated era is
nonzero for our model. Then we find the baryon-to-entropy ratio for a
radiation-dominated universe $\beta=\frac{1}{2}$, and when the Universe is
filled with perfect fluid, and cosmic dynamics are governed by the $f(R,Lm)$
theory of gravity. We determined the baryon-to-entropy ratio with specific
model $\frac{n_{B}}{s}\simeq 7.28749\times 10^{-11}$, which is in excellent
agreement with the observational value of $9\times 10^{-11}$. We assumed a
scale factor of the form $a(t)=Bt^{\beta}$, of $\beta=\frac{2}{3(1+\omega)}$,
and $B$ is a constant parameter for this work. Because $\beta=\frac{1}{2}$
were the appropriate values required to obtain a viable baryon-to-entropy
ratio, we assert that for such a scale factor, the power law part predominated
at early times, which is consistent with observations [35, 36]. Finally, we
conclude our research by investigating a more complete and generalized
baryogenesis interaction that is proportional to $\partial_{\mu}f(R,L_{m})$.
In this kind of interaction, our model produced a theoretical value in leading
order of $\frac{n_{B}}{s}\simeq 7.01\times 10^{-11}$, which is close to the
observed value as well as the value obtained in the first case.

## Data Availability Statement

There are no new data associated with this article.

## Acknowledgements

L.V.J. acknowledges the University Grant Commission (UGC), Govt. of India, New
Delhi, for awarding JRF (NTA Ref. No.: 191620024300). PKS acknowledges Science
and Engineering Research Board, Department of Science and Technology,
Government of India for financial support to carry out Research project No.:
CRG/2022/001847 and IUCAA, Pune, India for providing support through the
visiting Associateship program. We are very much grateful to the honorable
referee and to the editor for the illuminating suggestions that have
significantly improved our work in terms of research quality, and
presentation.

## References

  * [1] S. Burles et al.,  Phys. Rev. D 63, 063512 (2001).
  * [2] C. Bennett,  Ap. J. S. 148, 1 (2003)
  * [3] E. D. Stewart et al., Phys. Rev. D  54, 6032 (1996)
  * [4] M. Yamada,  Phys. Rev. D 93, 083516 (2016)
  * [5] K. Akita et al.,  J. Cosmol. Astropart. Phys. 4, 42 (2017)
  * [6] F. Takahashi et al.,  Phys. Rev. D 4, 083506 (2004)
  * [7] R. H. Brandenberge et al.,  Phys. Rev. D  68, 023505 (2003)
  * [8] A. D. Simone et al.,  J. Cosmol. Astropart. Phys. 8, 52 (2016)
  * [9] M. Trodden, Rev. Modern Phys. 71, 1463 (1999)
  * [10] D. E. Morrissey et al., New J. Phys. 14, 125003 (2012)
  * [11] E. W. Kolb et al., Phys. Rev. Lett. 77, 4290 (1996)
  * [12] A. D. Dolgov,  Phys. Rev. D 24, 1042 (1981)
  * [13] A. D. Sakharov,  JETP Lett. 2, 24 (1967).
  * [14] G. Lambiase at el., Phys. Rev. D  74, 087504 (2006).
  * [15] M. P. L. P. Ramos at el., Phys. Rev. D  96, 104024 (2017).
  * [16] S. D. Odintsov at el.,  Phys. Lett. B 760, 259-262 (2016).
  * [17] M. C. Bento at el.,  Phys. Rev. D 71, 123517 (2005).
  * [18] V. K. Oikonomou at el.,  Phys. Rev. D 94, 124005 (2016).
  * [19] S. Bhattacharjee, Int. J. Modern phys. A 36, 2150200 (2021).
  * [20] K. Nozari at el., Commun. Theor. phys. 70, 451 (2018).
  * [21] E. H. Baffou at el., Eur. Phys. J. C 112, 112 (2019).
  * [22] P. K. Sahoo & S.Bhattacharjee, Int. J. Theor. Phys. 59, 1451-1459 (2020).
  * [23] R. Saleem at el., Eur. Phys. J. C 137, 961 (2022).
  * [24] N. Azhar at el., Phys. Dark Uni. 30, 100724 (2020).
  * [25] S. Bhattacharjee at el., Eur. Phys. J. C  80, 289 (2020).
  * [26] S. Bhattacharjee at el., Phys.Dark Univ. 30, 100612 (2020).
  * [27] A. A. Starobinsky,  Phys. Lett. B 91, 99-102 (1980).
  * [28] T. Harko at el.,  Eur. Phys. J. C 70, 373-379 (2010).
  * [29] T. Harko at el.,  Eur. Phys. J. C 75, 386 (2015).
  * [30] L. V. Jaybhaye at el.,  Phys. Lett. B 831, 137148 (2022).
  * [31] R. V. Labato at el.,  Eur. Phys. J. C 82, 540 (2022).
  * [32] R. V. Labato at el.,  Eur. Phys. J. C 81, 1013 (2021).
  * [33] L. V. Jaybhaye at el.,  Int. J. Geom. Methods Mod. 19, 2250050 (2022).
  * [34] G. Lambiase at el.,  Int. J. Mod. Phys. D 22, 1330030 (2013).
  * [35] A. G. Riess at el., Astron. J. 116, 1009 (1998)
  * [36] G. Perlmutter at el., Astrophys. J. 517, 565 (1999)