††thanks: now at Lawrence Berkeley National Laboratory, Berkeley, CA 94720,
USA††thanks: now at University of Michigan, Ann Arbor, MI 48109, USA††thanks:
now at Tufts University, Medford, MA 02155, USA

The MINER$\nu$A Collaboration

# Tuning the GENIE Pion Production Model with MINERvA Data

P. Stowell  L. Pickering  C. Wilkinson  C. Wret  F. Akbar  AMU Campus,
Aligarh, Uttar Pradesh 202001, India  D. A. Andrade  Campus León y Campus
Guanajuato, Universidad de Guanajuato, Lascurain de Retana No. 5, Colonia
Centro, Guanajuato 36000, Guanajuato México.  M. V. Ascencio  Sección Física,
Departamento de Ciencias, Pontificia Universidad Católica del Perú, Apartado
1761, Lima, Perú  L. Bellantoni  Fermi National Accelerator Laboratory,
Batavia, Illinois 60510, USA  A. Bercellie  M. Betancourt  Fermi National
Accelerator Laboratory, Batavia, Illinois 60510, USA  A. Bodek  A. Bravar
University of Geneva, 1211 Geneva 4, Switzerland  H. Budd  G. Caceres  Centro
Brasileiro de Pesquisas Físicas, Rua Dr. Xavier Sigaud 150, Urca, Rio de
Janeiro, Rio de Janeiro, 22290-180, Brazil  T. Cai  M. F. Carneiro  Dept. of
Physics, Oregon State University, Corvallis, Oregon 97331, USA  J. Chaves
Dept. of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA
19104  H. da Motta  Centro Brasileiro de Pesquisas Físicas, Rua Dr. Xavier
Sigaud 150, Urca, Rio de Janeiro, Rio de Janeiro, 22290-180, Brazil  S. A.
Dytman  Dept. of Physics and Astronomy, University of Pittsburgh, Pittsburgh,
Pennsylvania 15260, USA  G.A. Díaz  Sección Física, Departamento de Ciencias,
Pontificia Universidad Católica del Perú, Apartado 1761, Lima, Perú  J. Felix
Campus León y Campus Guanajuato, Universidad de Guanajuato, Lascurain de
Retana No. 5, Colonia Centro, Guanajuato 36000, Guanajuato México.  L. Fields
Fermi National Accelerator Laboratory, Batavia, Illinois 60510, USA
Northwestern University, Evanston, Illinois 60208  A. Filkins  Dept. of
Physics, College of William & Mary, Williamsburg, Virginia 23187, USA  R. Fine
N. Fiza  IISER, Mohali, Knowledge city, Sector 81, Manauli PO 140306  R.
Galindo  Departamento de Física, Universidad Técnica Federico Santa María,
Avenida España 1680 Casilla 110-V, Valparaíso, Chile  H. Gallagher  Physics
Dept., Tufts University, Medford, Massachusetts 02155, USA  A. Ghosh
Departamento de Física, Universidad Técnica Federico Santa María, Avenida
España 1680 Casilla 110-V, Valparaíso, Chile  Centro Brasileiro de Pesquisas
Físicas, Rua Dr. Xavier Sigaud 150, Urca, Rio de Janeiro, Rio de Janeiro,
22290-180, Brazil  R. Gran  Dept. of Physics, University of Minnesota –
Duluth, Duluth, Minnesota 55812, USA  D. A. Harris  Fermi National Accelerator
Laboratory, Batavia, Illinois 60510, USA  S. Henry  S. Jena  IISER, Mohali,
Knowledge city, Sector 81, Manauli PO 140306  D. Jena  Fermi National
Accelerator Laboratory, Batavia, Illinois 60510, USA  J. Kleykamp  M. Kordosky
Dept. of Physics, College of William & Mary, Williamsburg, Virginia 23187, USA
D. Last  Dept. of Physics and Astronomy, University of Pennsylvania,
Philadelphia, PA 19104  T. Le  Physics Dept., Tufts University, Medford,
Massachusetts 02155, USA  Rutgers, The State University of New Jersey,
Piscataway, New Jersey 08854, USA  X.-G. Lu  Oxford University, Dept. of
Physics, Oxford, United Kingdom  E. Maher  Massachusetts College of Liberal
Arts, 375 Church Street, North Adams, MA 01247  S. Manly  W. A. Mann  Physics
Dept., Tufts University, Medford, Massachusetts 02155, USA  C. M. Marshall  K.
S. McFarland  Fermi National Accelerator Laboratory, Batavia, Illinois 60510,
USA  B. Messerly  Dept. of Physics and Astronomy, University of Pittsburgh,
Pittsburgh, Pennsylvania 15260, USA  J. Miller  Departamento de Física,
Universidad Técnica Federico Santa María, Avenida España 1680 Casilla 110-V,
Valparaíso, Chile  J. G. Morfín  Fermi National Accelerator Laboratory,
Batavia, Illinois 60510, USA  J. Mousseau  University of Florida, Dept. of
Physics, Gainesville, FL 32611  D. Naples  Dept. of Physics and Astronomy,
University of Pittsburgh, Pittsburgh, Pennsylvania 15260, USA  J. K. Nelson
Dept. of Physics, College of William & Mary, Williamsburg, Virginia 23187, USA
C. Nguyen  University of Florida, Dept. of Physics, Gainesville, FL 32611  A.
Norrick  Dept. of Physics, College of William & Mary, Williamsburg, Virginia
23187, USA  Nuruzzaman  Rutgers, The State University of New Jersey,
Piscataway, New Jersey 08854, USA  Departamento de Física, Universidad Técnica
Federico Santa María, Avenida España 1680 Casilla 110-V, Valparaíso, Chile  V.
Paolone  Dept. of Physics and Astronomy, University of Pittsburgh, Pittsburgh,
Pennsylvania 15260, USA  G. N. Perdue  Fermi National Accelerator Laboratory,
Batavia, Illinois 60510, USA  M. A. Ramírez  Campus León y Campus Guanajuato,
Universidad de Guanajuato, Lascurain de Retana No. 5, Colonia Centro,
Guanajuato 36000, Guanajuato México.  R. D. Ransome  Rutgers, The State
University of New Jersey, Piscataway, New Jersey 08854, USA  H. Ray
University of Florida, Dept. of Physics, Gainesville, FL 32611  D. Rimal
University of Florida, Dept. of Physics, Gainesville, FL 32611  P. A.
Rodrigues  University of Mississippi, Oxford, Mississippi 38677, USA  D.
Ruterbories  H. Schellman  Dept. of Physics, Oregon State University,
Corvallis, Oregon 97331, USA  Northwestern University, Evanston, Illinois
60208  C. J. Solano Salinas  Universidad Nacional de Ingeniería, Apartado
31139, Lima, Perú  H. Su  Dept. of Physics and Astronomy, University of
Pittsburgh, Pittsburgh, Pennsylvania 15260, USA  M. Sultana  V. S. Syrotenko
Physics Dept., Tufts University, Medford, Massachusetts 02155, USA  E.
Valencia  Dept. of Physics, College of William & Mary, Williamsburg, Virginia
23187, USA  Campus León y Campus Guanajuato, Universidad de Guanajuato,
Lascurain de Retana No. 5, Colonia Centro, Guanajuato 36000, Guanajuato
México.  J. Wolcott  M. Wospakrik  University of Florida, Dept. of Physics,
Gainesville, FL 32611  B. Yaeggy  Departamento de Física, Universidad Técnica
Federico Santa María, Avenida España 1680 Casilla 110-V, Valparaíso, Chile  L.
Zazueta  Dept. of Physics, College of William & Mary, Williamsburg, Virginia
23187, USA  D. Zhang  Dept. of Physics, College of William & Mary,
Williamsburg, Virginia 23187, USA

###### Abstract

Faced with unresolved tensions between neutrino interaction measurements at
few-GeV neutrino energies, current experiments are forced to accept large
systematic uncertainties to cover discrepancies between their data and model
predictions. In this paper, the widely used pion production model in GENIE is
compared to four MINERvA charged current pion production measurements using
NUISANCE. Tunings, i.e., adjustments of model parameters, to help match GENIE
to MINERvA and older bubble chamber data are presented here. We find that
scattering off nuclear targets as measured in MINERvA is not in good agreement
with expectations based upon scattering off nucleon (hydrogen or deuterium)
targets in existing bubble chamber data. An additional ad hoc correction for
the low-$Q^{2}$ region, where collective nuclear effects are expected to be
large, is presented. While these tunings and corrections improve the agreement
of GENIE with the data, the modeling is imperfect. The development of these
tunings within the NUISANCE framework allows for straightforward extensions to
other neutrino event generators and models, and allows omitting and including
new data sets as they become available.

##  I Introduction

In recent years, experimental groups have started to publish neutrino
interaction cross-section measurements on nuclear targets in terms of
measurable final state particle content, instead of inferred initial
interaction channels. This avoids the problem of correcting for complex
nuclear effects to make a measurement in terms of the initial interaction
channels. For example, events with only a single pion can be produced by the
decay of hadronic resonances formed at the primary neutrino interaction,
followed by loss of a nucleon from the resonance’s decay as a result of final
state interactions (FSI) within the nuclear medium. Such events can also be
produced by other sequences of interactions, such as a deep inelastic
collision where only a single pion is produced after FSI. A measurement of
charged current events with one identified pion in the final state is a
benchmark for models, independent of the details of how each model assesses
any particular interaction channel’s contributions to that final state. The
limitation of giving results in terms of final state particle content is that
FSI are important, and result in the contribution of many different
interaction channels into a specific final state.

There are tensions between published data from the T2K, MiniBooNE, and MINERvA
experiments Mosel and Gallmeister (2017); Alvarez-Ruso  _et al._ (2018); R.
González-Jiménez and Jachowicz (2018); Wilkinson  _et al._ (2016); Mahn  _et
al._ (2018). These tensions exist in the charged current production of both
zero and one pion final states, and a model has yet to emerge that can
reliably simulate all experiments at once. This is troubling, as current and
future neutrino oscillation experiments require a cross section model which is
predictive across the range of energies covered by these experiments and for a
variety of targets.

The differences in neutrino fluxes, scattering targets, available phase space
and signal definitions between experiments make it difficult to diagnose the
exact causes of disagreement within the global data set. In particular, as
results must be averaged over the neutrino flux distribution of each
experiment, it is difficult to disentangle the energy dependence of an
observed deficiency in a particular model, and decide how uncertainties should
be propagated in neutrino energy. Tensions between measurements from a single
experiment can uncover fundamental problems with a model which should be
addressed, before considering the more difficult issue of developing, or
empirically tuning a model which fits data from multiple experiments. In this
work, we employ published MINERvA pion production data. The cross-section
measurements utilized in this effort have not been reanalyzed or modified in
any way.

NUISANCE Stowell  _et al._ (2017) was developed to provide the neutrino
scattering community with a flexible framework in which various neutrino
interaction generators can be validated and empirically tuned to data. Its
structure allows for generator tunings to be easily adapted to account for
changes in the underlying model or data. In this work, the default pion
production model in the GENIE Andreopoulos  _et al._ (2009, 2015) neutrino
interaction generator is tuned to MINERvA data. Although more sophisticated
pion production models exist (e.g. Kabirnezhad (2018); Leitner  _et al._
(2009); Nieves  _et al._ (2011); Hernadez and Nieves (2017)), GENIE is widely
used by the neutrino scattering community, and its model uncertainties have a
central importance to the field. Although the work is only directly applicable
to one generator, the methods developed in this paper are easily adaptable to
different generators. All the data and methods are publicly available and
integrated into the open source NUISANCE framework, facilitating similar
studies using other generators and models.

In Section II, the data are reviewed and the goodness-of-fit test statistic is
defined for the tuning process. Section III describes the default GENIE pion
production model, and reviews comparisons of this model to data. In Section
IV, the parameter reweighting package in GENIE is discussed along with the
specific parameters tuned therewith. We also discuss other corrections to the
GENIE model made to improve agreement with bubble chamber data Rodrigues  _et
al._ (2016); Wilkinson  _et al._ (2014). In Section V, we tune additional
systematic parameters in GENIE to improve agreement with the MINERvA data in
combination with the bubble chamber data. In Section VI, additional
low-$Q^{2}$ ad hoc corrections are added to the model to resolve observed
tensions, motivated by the need for similar corrections observed at both MINOS
Adamson  _et al._ (2015) and MiniBooNE Aguilar-Arevalo  _et al._ (2011).
Finally, in Section VII we present our conclusions.

##  II Data included in the fits

We tune to four of MINERvA’s published charged current pion production
measurements taken on a polystyrene scintillator target:
$\nu_{\mu}$CC1$\pi^{\pm}$ Eberly  _et al._ (2015), $\nu_{\mu}$CC$N\pi^{\pm}$
McGivern  _et al._ (2016), $\nu_{\mu}$CC$1\pi^{0}$ Altinok  _et al._ (2017)
and $\bar{\nu}_{\mu}$CC$1\pi^{0}$ McGivern  _et al._ (2016), summarized in
Table 1111In “$\nu_{\mu}$CC$N\pi^{\pm}$ ”, the $N$ indicates one or more
identified pions and does not refer to a nucleon.. The MINERvA detector Aliaga
_et al._ (2014) does not determine the polarity of charged pions. The fraction
of $\pi^{-}$ in $\nu_{\mu}$CC1$\pi^{\pm}$ sample is small ($\sim$ 1%).
Furthermore, the $\nu_{\mu}$CC1$\pi^{\pm}$ and $\nu_{\mu}$CC$N\pi^{\pm}$
signal definition allows for any number of neutral pions. Approximately 3% of
the MINERvA $\nu_{\mu}$CC1$\pi^{\pm}$ signal events have at least one neutral
pion in the final state. All four analyses include signal definition cuts on
the true “reconstructed” mass of the hadronic system assuming the struck
nucleon is at rest, $W_{\mathrm{rec}}$, and the true neutrino energy
$E_{\nu}$.

Channel |  $\nu_{\mu}$CC1$\pi^{\pm}$ Eberly  _et al._ (2015) |  $\nu_{\mu}$CC$N\pi^{\pm}$ McGivern  _et al._ (2016) |  $\nu_{\mu}$CC$1\pi^{0}$ Altinok  _et al._ (2017) |  $\bar{\nu}_{\mu}$CC$1\pi^{0}$ McGivern  _et al._ (2016)
---|---|---|---|---
$\mathrm{N_{bins}}$ $\;{p_{\mu}}$ | 8 | 9 | 8 | 9
$\mathrm{N_{bins}}$ $\;{\theta_{\mu}}$ | 9 | 9 | 9 | 9
$\mathrm{N_{bins}}$ $\;{T_{\pi}}$ | 7 | 7 | 7 | 7
$\mathrm{N_{bins}}$ $\;{\theta_{\pi}}$ | 14 | 14 | 11 | 11
$\mathrm{N_{bins}}$ total | 38 | 39 | 35 | 36
Signal definition | 1$\pi^{\pm}$, $\geq$ 0$\pi^{0}$ |  $>0$$\pi^{\pm}$, $\geq$ 0$\pi^{0}$ | 1$\pi^{0}$, 0$\pi^{\pm}$ | 1$\pi^{0}$, 0$\pi^{\pm}$
| 1$\mu^{-}$ | 1$\mu^{-}$ | 1$\mu^{-}$ | 1$\mu^{+}$
|  $W_{\mathrm{rec}}$$<1.4~{}\mbox{GeV}$ |  $W_{\mathrm{rec}}$$<1.8~{}\mbox{GeV}$ |  $W_{\mathrm{rec}}$$<1.8~{}\mbox{GeV}$ |  $W_{\mathrm{rec}}$$<1.8~{}\mbox{GeV}$
| — | — | $\theta_{\mu}<25^{\circ}$ | —
Table 1: Summary of the measurements used in this analysis. $W_{\mathrm{rec}}$
is the true reconstructed hadronic mass assuming the struck nucleon is at
rest. None of the measurements veto on activity other than the $\mu$ and $\pi$
in their signal definition, and all selections require
$1.5<E_{\nu}<20~{}\mbox{GeV}$.

The kinematic variable distributions used in this work are the momentum and
angle of the outgoing muon with respect to the incoming neutrino beam,
$p_{\mu}$ and $\theta_{\mu}$, and the kinetic energy and angle of the outgoing
pion with respect to the incoming neutrino beam, $T_{\pi}$ and $\theta_{\pi}$.
In the $\nu_{\mu}$CC$N\pi^{\pm}$ channel, where there is at least one
$\pi^{\pm}$ in the final state, there is one entry in the distributions of
$\theta_{\pi}$ and $T_{\pi}$ for each $\pi^{\pm}$ in an event. The data are
reported as efficiency corrected results unfolded to true kinematic variables,
which may introduce model dependence. This is notably problematic in regions
of low efficiency—present in the charged pion channels at $\theta_{\pi}\sim
90^{\circ}$, $T_{\pi}<50~{}\mbox{MeV}$ and $T_{\pi}>350~{}\mbox{MeV}$, where
the signal efficiency is zero Eberly  _et al._ (2015). The pion selection
cuts, not present in the signal definition, remove about 50% of signal events,
with little dependence upon the muon variables, but a clear impact on the
shape of the pion kinematic variables.

The published cross-sections are one dimensional with correlations provided
between the bins within each distribution. No correlations are provided
between measurements of different final states, or between different one-
dimensional projections of the same measurement. These correlations are
expected to be large, coming predominantly from flux and detector
uncertainties. Additionally, the $\nu_{\mu}$CC1$\pi^{\pm}$ event sample is a
subset ($\sim$64%) of the $\nu_{\mu}$CC$N\pi^{\pm}$ event sample, and
including both channels introduces a statistical correlation. Not assessing
correlations between the distributions, while common practice in this field,
is a limitation when tuning models to multiple data sets. It introduces a bias
in the $\chi^{2}$ statistic that is difficult to quantify, and requires
imposing ad hoc uncertainties Wilkinson  _et al._ (2016) as the test-statistic
is not expected to follow a $\chi^{2}$ distribution for the given degrees of
freedom.

The covariance matrices contain a flux-dominated normalization component which
we expect to be fully correlated across all distributions. To account for the
correlated uncertainty, we use the full covariance matrix, $M_{ij}$, for the
$p_{\mu}$ distribution and shape-only covariance matrices, $S_{ij}$, for the
other three distributions in each of the topologies. Whilst any distribution
could set the normalization constraint, the shape of the $p_{\mu}$
distribution for each channel was chosen since it was found to be relatively
insensitive to model variations and had good shape agreement with the data.
The joint $\chi^{2}$ is therefore defined as the sum of the full $p_{\mu}$
$\chi^{2}$ and shape-only $\theta_{\mu}$, $T_{\pi}$ and $\theta_{\pi}$
$\chi^{2}$’s:

$\displaystyle\chi^{2}$
$\displaystyle=\sum_{ij}^{N_{p_{\mu}}}\Delta_{i}(M^{-1})_{ij}\Delta_{j}$
$\displaystyle+\sum_{kij}^{N_{k}}\Delta^{S}_{k,i}(S^{-1})_{ij}\Delta^{S}_{k,j}$
(1)

where $i$ and $j$ are bin indices,

$\displaystyle\Delta_{i}$ $\displaystyle=d_{p_{\mu},i}-m_{p_{\mu},i}$ (2)
$\displaystyle\Delta^{S}_{k,i}$
$\displaystyle=d_{k,i}-\left(m_{k,i}\times\frac{\sum_{j}d_{k,j}}{\sum_{j}m_{k,j}}\right),$
(3)

and $d_{k,i}$ and $m_{k,i}$ are the data and MC values respectively for the
$i^{th}$ bin in the $k^{th}$ distribution. The shape-only covariance matrices
are provided in the public data release for the $\nu_{\mu}$CC1$\pi^{\pm}$ and
$\nu_{\mu}$CC$1\pi^{0}$ measurements, and the method of Ref. Katori (2008)
(section 10.6.3) was used to extract them for the $\nu_{\mu}$CC$N\pi^{\pm}$
and $\bar{\nu}_{\mu}$CC$1\pi^{0}$ channels.

##  III Pion Production in GENIE

This analysis begins with version 2.12.6 of GENIE, which is close to what is
used by MINERvA, T2K, NOvA and MicroBooNE. We use the Smith-Moniz relativistic
Fermi gas (RFG) model Smith and Moniz (1972) with an added high momentum tail
as per Bodek and Ritchie Bodek and Ritchie (1981). The Valencia random phase
approximation screening Nieves  _et al._ (2004) is applied as a weight to
quasielastic events. The two-particle two-hole process is simulated using the
Valencia model Nieves  _et al._ (2011); Gran  _et al._ (2013). MINERvA
currently uses a modification of v2.8.4 Altinok  _et al._ (2017); Patrick  _et
al._ (2018); Mislivec  _et al._ (2018) with an increased rate for the Valencia
two-particle two-hole process; that modification is not used here. An
important difference in single pion production between v2.8.x and v2.12.x is
the angular distributions of single pion events in the Rein-Sehgal model,
discussed below. A sample of 2.5 million events were generated using the
MINERvA flux predictions Aliaga  _et al._ (2016), a polystyrene target and the
official GENIE 2.12.6 splines Andreopoulos  _et al._ (2017).

To simulate pion production, GENIE uses the Rein-Sehgal (RS) model Rein and
Sehgal (1981) with a hadronic invariant mass cut of $W\leq 1.7~{}\mbox{GeV}$.
Of the 18 resonances in the RS model, the $\Delta(1600)$ and $N(1990)$ were
not included due to their unclear experimental status at the time of
implementation. Resonance-resonance and resonance-nonresonance interference
terms are not included. Lepton mass terms are only included in calculating
phase space limits and are neglected when calculating the cross sections. A
discussion of the limitations of this simplification can be found in Ref.
Graczyk and Sobczyk (2008). In earlier versions—including v2.8.4—the pion-
nucleon distribution was isotropic in the resonance rest frame, but was
changed in 2.12.x. Here we use the non-isotropic model as our default and
reweight to the isotropic distribution, explained later. The RS nonresonant
background is not used by GENIE; rather, a deep inelastic scattering (DIS)
model is extended to cover that invariant mass region. The DIS model uses the
Bodek-Yang parametrization Bodek and Yang (2003), and the AGKY model to
describe hadronization Yang  _et al._ (2009). In the AGKY model, the KNO model
Koba  _et al._ (1972) is used for $W\leq 2.3~{}\mbox{GeV}$ and PYTHIA T.
Sjostrand and Skands (2008) is used for $W\geq 3.0~{}\mbox{GeV}$, with a
smooth transition in between the two, implemented by randomly selecting the
results of one model or the other for each event.

In addition to pion production on a single nucleon, it is also possible for a
neutrino to produce a pion by scattering coherently off the nucleus. GENIE
uses the Rein-Sehgal coherent pion production model Rein and Sehgal (1983,
2007), including the effect of lepton masses in the cross-section calculation.
MINERvA has found that the RS coherent pion production model needs to be
suppressed by $\sim 50\%$ at $T_{\pi}$ $<500~{}\mbox{MeV}$ to agree with data
Mislivec  _et al._ (2018). This correction also moves the shape of the
$T_{\pi}$ spectrum closer to the predictions of the Berger-Sehgal coherent
model Berger and Sehgal (2009). The $\nu_{\mu}$CC1$\pi^{\pm}$ channel has a
small contribution from coherent production in the lowest $Q^{2}$ bins but the
inclusion of this suppression has only a small effect on the MC predictions.
To maintain a model similar to that currently being used by MINERvA, this
suppression is included in the analysis presented in Sections IV and onwards.

The “hA Intranuke” effective cascade model Dytman (2009) is used to model pion
and nucleon FSI. In this model, the effect of intranuclear scattering is
parameterized as a single cascade step applied to each particle emanating from
the primary interaction. This model steps hadrons through a nucleus of radius
$r\sim A^{1/3}$ and a nuclear density function derived from electron
scattering data. The hadron’s mean free path is determined from tabulated
hadron-proton and hadron-neutron cross sections Arndt  _et al._ (2006). The
probability to interact with the nucleus is high; it is e.g. $\sim$73% for a
pion from an $E_{\nu}$ = 3 GeV quasielastic event in carbon. When a FSI
occurs, the possible interactions (absorption, pion production, knockout,
charge exchange, elastic scatter) are chosen according to their proportions
for iron.

Distribution | Channel | $\mathrm{N_{bins}}$ | Default | ANL/BNL | FrAbs Tune | FrInel Tune
---|---|---|---|---|---|---
$p_{\mu}$ (Rate) |  $\nu_{\mu}$CC1$\pi^{\pm}$ | 8 | 19.1 | 13.8 | 12.0 | 12.3
|  $\nu_{\mu}$CC$N\pi^{\pm}$ | 9 | 35.4 | 19.5 | 26.1 | 26.8
|  $\nu_{\mu}$CC$1\pi^{0}$ | 8 | 11.1 | 19.6 | 19.0 | 19.3
|  $\bar{\nu}_{\mu}$CC$1\pi^{0}$ | 9 | 7.4 | 6.4 | 6.2 | 6.3
$\theta_{\mu}$ (Shape) |  $\nu_{\mu}$CC1$\pi^{\pm}$ | 9 | 7.1 | 12.4 | 7.5 | 7.4
|  $\nu_{\mu}$CC$N\pi^{\pm}$ | 9 | 4.5 | 10.4 | 4.0 | 4.1
|  $\nu_{\mu}$CC$1\pi^{0}$ | 9 | 35.1 | 71.5 | 44.5 | 45.6
|  $\bar{\nu}_{\mu}$CC$1\pi^{0}$ | 9 | 9.3 | 14.0 | 10.2 | 10.3
$T_{\pi}$ (Shape) |  $\nu_{\mu}$CC1$\pi^{\pm}$ | 7 | 2.9 | 2.6 | 2.5 | 2.3
|  $\nu_{\mu}$CC$N\pi^{\pm}$ | 7 | 39.8 | 34.7 | 31.2 | 29.4
|  $\nu_{\mu}$CC$1\pi^{0}$ | 7 | 28.3 | 31.4 | 30.9 | 29.9
|  $\bar{\nu}_{\mu}$CC$1\pi^{0}$ | 7 | 19.3 | 17.9 | 16.6 | 16.0
$\theta_{\pi}$ (Shape) |  $\nu_{\mu}$CC1$\pi^{\pm}$ | 14 | 25.4 | 26.5 | 13.0 | 12.6
|  $\nu_{\mu}$CC$N\pi^{\pm}$ | 14 | 11.7 | 11.1 | 6.9 | 6.2
|  $\nu_{\mu}$CC$1\pi^{0}$ | 11 | 13.5 | 15.0 | 8.3 | 8.9
|  $\bar{\nu}_{\mu}$CC$1\pi^{0}$ | 11 | 5.7 | 5.9 | 3.4 | 3.5
Total $\chi^{2}$ |  | 148 | 275.6 | 312.7 | 242.3 | 240.7
Table 2: Channel by channel contributions to the $\chi^{2}$ at different
stages of the tuning process.

Default GENIE predictions separated by nucleon level interaction channels for
the MINERvA data are shown in Fig. 1. The shape of the $p_{\mu}$ distributions
agree well with the data for all four measurements. However, the model
overestimates the cross section for $\pi^{\pm}$ production and as a result the
$\chi^{2}$ for the $\nu_{\mu}$CC1$\pi^{\pm}$ and $\nu_{\mu}$CC$N\pi^{\pm}$,
given in the fourth column (“Default”) of Table 2, are large. The model
overestimates $\theta_{\mu}$ below $<5^{\circ}$ in the $\pi^{0}$ channels,
although it does correctly predict the shape of the $\theta_{\mu}$
distribution in the $\pi^{\pm}$ channels. The model underestimates the
production rate at large $\theta_{\mu}$ in $\nu_{\mu}$CC$1\pi^{0}$. The shape
of the $T_{\pi}$ distribution is in larger disagreement for
$\nu_{\mu}$CC$N\pi^{\pm}$ data than for $\nu_{\mu}$CC1$\pi^{\pm}$. Since the
$\nu_{\mu}$CC$N\pi^{\pm}$ distributions summed over all identified
$\pi^{\pm}$, redistributing kinetic energy between $\pi^{\pm}$ in events with
more than one $\pi^{\pm}$ could resolve some of this tension. The $\pi^{0}$
channels are under-predicted at low $T_{\pi}$. Finally, GENIE predictions are
too high in magnitude at $\theta_{\pi}\approx 50^{\circ}$ in both the
$\nu_{\mu}$CC1$\pi^{\pm}$ and $\nu_{\mu}$CC$N\pi^{\pm}$ channels, and the
prediction has the wrong shape in the $\nu_{\mu}$CC1$\pi^{\pm}$ channel.
Comparisons using the transport theory based GiBUU model Buss  _et al._ (2012)
show similar shape disagreements despite GiBUU’s use of an advanced
semiclassical cascade model to simulate FSI Mosel and Gallmeister (2017).

Figure 1: GENIE 2.12.6 Default model predictions compared to MINERvA data.
Colors correspond to particle content at the nucleon interaction. “Other” is
dominated by coherent pion production. “MC Shape” shows the total MC
prediction after it has been normalized to match the total data normalization.
In the case of the shape-only distributions ($\theta_{\mu}$, $T_{\pi}$,
$\theta_{\pi}$) the shape-only $\chi^{2}$/$\mathrm{N_{bins}}$ values are
shown. All cross sections are per nucleon.

Each of the measurements are shown as MC/data ratio distributions in Fig. 2.
Similar comparisons between the MiniBooNE and MINERvA experiments are found in
Ref. Mahn  _et al._ (2018). The shape-only data sets ($\theta_{\mu}$,
$\theta_{\pi}$, $T_{\pi}$) were normalized to match the data before the ratio
was taken and the error bars in Fig. 2 reflect the extracted shape-only
uncertainties on the data, so that the distributions reflect their
contributions to the total $\chi^{2}$.

Figure 2: MC/data ratios for the default GENIE predictions. The $p_{\mu}$
distribution provides a rate comparison in the $\chi^{2}$ calculation; the
other distributions are treated as shape-only, i.e., the MC is normalized to
match the data and the uncertainties are from the shape-only covariance
matrix.

##  IV Tunable parameters in the GENIE Model

The GENIE event generator allows assessment of systematic uncertainties
through the GENIE reweighting package. A large number of event weighting
“dials” are included to allow model uncertainties to be evaluated. The dials
adjusted in this note are summarized in Table 3 and are chosen because of
their connection to the kinematic variables and interaction modes studied
herein.

Parameter | Default Value | GENIE parameter name
---|---|---
CC Resonant Axial Mass ($M_{\mathrm{A}}^{\mathrm{res}}$) | $1.12\pm 0.22~{}\mbox{GeV}$ | MaCCRES
CC Resonant Normalization ($\mathrm{NormRes}$) | $100\pm 20~{}\mbox{\%}$ | NormCCRES
CC1$\pi$ Nonresonant Normalization (NonRes$1\pi$) | $100\pm 50~{}\mbox{\%}$ | NonRESBGvnCC1pi
|  | NonRESBGvpCC1pi
|  | NonRESBGvbarnCC1pi
|  | NonRESBGvbarpCC1pi
CC2$\pi$ Nonresonant Normalization (NonRes$2\pi$) | $100\pm 50~{}\mbox{\%}$ | NonRESBGvnCC2pi
|  | NonRESBGvpCC1pi
|  | NonRESBGvbarnCC1pi
|  | NonRESBGvbarpCC1pi
Pion Angular Emission ($\pi$-iso) | 0 (RS) | Theta_Delta2Npi
Pion Absorption FSI Fraction (FrAbs) | $100\pm 30~{}\mbox{\%}$ | FrAbs_pi
Pion Inelastic FSI Fraction (FrInel) | $100\pm 40~{}\mbox{\%}$ | FrInel_pi
Table 3: Summary of the GENIE dials optimized in this note, their default
values and the uncertainties recommended by the GENIE collaboration. We do not
use the defaults for $M_{\mathrm{A}}^{\mathrm{res}}$, $\mathrm{NormRes}$ and
NonRes$1\pi$ and instead impose central values and uncertainties from tunings
to ANL and BNL data as described in the text.

Experiments often use variations in the charged-current resonant axial mass,
$M_{\mathrm{A}}^{\mathrm{res}}$, as a systematic uncertainty which varies both
the normalization and $Q^{2}$ shape of resonant interactions along with
variations in a total resonant cross-section normalization dial,
$\mathrm{NormRes}$. Variations in $\mathrm{NormRes}$ approximates the
behaviour of varying $F_{\mathrm{A}}(0)$ in the axial form factor in the Rein-
Sehgal model. Since low $\theta_{\mu}$ correlates with low $Q^{2}$, variations
in $M_{\mathrm{A}}^{\mathrm{res}}$ have the largest effect on the shape of the
muon angular distributions as shown in Fig. 3, and have a small effect on the
$\theta_{\pi}$ spectrum.

Figure 3: The effect varying the $M_{\mathrm{A}}^{\mathrm{res}}$ dial on the
default GENIE prediction for $\theta_{\mu}$. The red bands show the variation
to the total rate and shape. The green bands are obtained by normalizing the
reweighted curves to the default predicted rate to highlight the smaller
effect the dial has on the shape of the distributions.

Dials are available to vary the normalization of the nonresonant $1\pi$
production channels in GENIE (e.g. NonRESBGvnCC1pi, NonRESBGvpCC1pi) but each
dial introduces similar modifications to the predictions. To reduce the number
of free parameters in the fit described in Section V, these dials were grouped
into a single background scaling for nonresonant $1\pi$ production,
NonRes$1\pi$, following the approach in Ref. Rodrigues  _et al._ (2016);
Wilkinson  _et al._ (2014). A similar treatment was also applied to
nonresonant $2\pi$ production, NonRes$2\pi$, with the neutrino and
antineutrino related parameters assumed to be 100% correlated in both cases.
The effects of varying the nonresonant contributions are shown in Fig. 4.
Variations in the NonRes$2\pi$ dial introduce a large change in normalization
for the $\nu_{\mu}$CC$N\pi^{\pm}$ channel and has a minor effect in the other
single pion channels as the fraction of multi-$\pi$ events is small.

Figure 4: The effect of varying the NonRes$1\pi$ dial on the default GENIE
prediction for $T_{\pi}$ (top) and of varying the NonRes$2\pi$ dial on the
prediction for $\theta_{\pi}$ (bottom).

Reanalysis of data from ANL and BNL bubble chambers has provided a tuning of
GENIE’s single pion production model on free nucleons. The work showed that a
small shift in $M_{\mathrm{A}}^{\mathrm{res}}$ was required to model the
low-$Q^{2}$ region and a large suppression of the nonresonant $\pi$ production
($-54\%$) was required to match the observed cross sections of $\pi^{+}$ and
$\pi^{0}$ production. The reanalysis used the measured ratios of the rates of
single $\pi$ production to CCQE measurements to cancel errors in the flux. We
note that by using CCQE data multiple times, they introduce hidden
correlations which may have a small effect on the postfit uncertainties.
However, as the single pion statistical uncertainties at ANL Radecky and
others (ANL) and BNL Kitagaki and others (BNL) were magnitudes higher than the
CCQE statistical uncertainty Barish and others (ANL); Kitagaki and others
(BNL), the effect was neglected in that work, and is also neglected here. The
resulting parameter tunes shown in Table 4 and Fig. 5 have been partially
adopted by MINERvA and NOvA which both apply the nonresonant rescaling of 43%
but leave the other parameters unchanged.

Parameter | GENIE default | ANL/BNL tune
---|---|---
$M_{\mathrm{A}}^{\mathrm{res}}$ [GeV] | $1.12\pm 0.22$ | $0.94\pm 0.05$
$\mathrm{NormRes}$ [%] | $100\pm 20$ | $115\pm 7$
NonRes$1\pi$ [%] | $100\pm 50$ | $43\pm 4$
Table 4: Prefit and best fit central values and uncertainties from tuning
GENIE to the ANL/BNL pion production measurements. The prefit uncertainties
are those recommended by the GENIE collaboration. The tuned values are used as
penalty terms with supplied covariance matrix of Fig 5.

Figure 5: Correlation matrix from tuning GENIE to reproduce the ANL/BNL pion
production measurements included in our $\chi^{2}$ penalty term.

.

Fig. 6 shows MINERvA data and the predictions of GENIE when its output has
been reweighted to reflect the parameter changes of Table 4. The channel-by-
channel contributions to the $\chi^{2}$ are given in the fifth column
(“ANL/BNL”) of Table 2. Incorporating the parameter changes improves the total
normalization agreement in the $p_{\mu}$ distributions for
$\nu_{\mu}$CC1$\pi^{\pm}$ and $\nu_{\mu}$CC$N\pi^{\pm}$. The $\chi^{2}$ for
the $p_{\mu}$ distribution is also improved in the
$\bar{\nu}_{\mu}$CC$1\pi^{0}$ channel, even though the ANL/BNL data is from
neutrino interactions only. The $\chi^{2}$ for the $p_{\mu}$ distribution in
the $\nu_{\mu}$CC$1\pi^{0}$ channel is somewhat worse as the parameter tunes
reduce the predicted nucleon $\nu_{\mu}$CC$1\pi^{0}$ cross section. The
modification of $M_{\mathrm{A}}^{\mathrm{res}}$ shifts the $\theta_{\mu}$
predictions to lower values, increasing the $\chi^{2}$ contributions. The
$T_{\pi}$ and $\theta_{\pi}$ distributions change mostly by normalization,
having a smaller effect on the $\chi^{2}$. The overall agreement of GENIE with
MINERvA data is not improved by incorporating the ANL/BNL information. Indeed,
the total $\chi^{2}$ increases, largely because of the $\chi^{2}$
contributions from the $\theta_{\mu}$ distributions.

Figure 6: GENIE ANL/BNL single pion tuning model predictions compared to
MINERvA data. The distributions have been weighted to the ANL/BNL tuning
parameter set, and have had the coherent pion correction applied. Colors
correspond to particle content at the nucleon interaction. “Other” is
dominated by coherent pion production. “MC Shape” shows the total MC
prediction after it has been normalized to match the total data normalization.
In the case of the shape-only distributions ($\theta_{\mu}$, $T_{\pi}$,
$\theta_{\pi}$) the shape-only $\chi^{2}$/$\mathrm{N_{bins}}$ values are
shown. All cross sections are per nucleon.

GENIE provides a dial that influences the resonances’ decay into the pion-
nucleon system in the resonance rest frame, $\pi$-iso, and allows events to be
reweighted continuously between the default anisotropic distribution
($\pi$-iso = 0) and the isotropic distribution ($\pi$-iso = 1). The Adler
angle222The Adler angle is the angle between the pion and the three momentum
transfer in the resonance rest frame Adler (1968). is highly sensitive to the
$\pi$-iso parameter and has been measured by neutrino induced pion production
experiments on single nucleons, such as ANL Radecky and others (ANL), BNL
Kitagaki and others (BNL), BEBC Allen  _et al._ (1980); Allasia  _et al._
(1983) and FNAL Bell.  _et al._ (1978). Nucleon data strongly prefers an
anisotropic process, as shown in Fig. 7. Nonetheless, $\pi$-iso has some
impact, albeit one that does depend on how FSI are modelled, on the shape of
MINERvA $\theta_{\pi}$ and $T_{\pi}$ distributions, seen in the bottom of Fig.
7, and was therefore included in this work.

Figure 7: Effect of varying the $\pi$-iso dial on pion angular distributions
for BNL $\nu_{\mu}$CC1$\pi^{+}$1p data (top) and MINERvA
$\nu_{\mu}$CC1$\pi^{+}$ (bottom).

The GENIE hA model for FSI has uncertainties from the $\pi-A$ cross-section
data to which the model was tuned. The total $\pi-A$ cross section has a
stronger constraint than each of the individual interaction cross sections, so
GENIE provides dials to vary the fractional contribution of each component.
The available fractional dials are pion absorption (FrAbs), pion inelastic
scattering (FrInel), pion elastic scattering, pion charge exchange and pion
production.

##  V Tuning the GENIE Model

Fig. 6 and Table 2 show the unsatisfactory agreement of the GENIE prediction
against MINERvA data. The disagreement worsens after incorporating the prior
constraint from ANL and BNL data; this correction, based on nucleon data, is
inadequate. This section describes fits that improve the agreement with
MINERvA data. The parameters $M_{\mathrm{A}}^{\mathrm{res}}$,
$\mathrm{NormRes}$ and NonRes$1\pi$ are included in the fits with a penalty
term added to the $\chi^{2}$ from the ANL and BNL data. The penalty term uses
the covariance, $M$, shown in Fig. 5:

$\displaystyle\chi^{2}_{\textrm{pen}}=\sum_{i,j}^{N=3}(x_{i}-f_{i})\left(M^{-1}\right)_{ij}(x_{j}-f_{j}),$
(4)

where $x_{i}$ are the parameter values $i$ at each iteration of the fit, and
$f_{i}$ are the parameter values from the fit to ANL and BNL data. The GENIE
default model is strongly disfavoured with $\chi^{2}_{\textrm{pen}}$ = 299.3,
but changing NonRes$1\pi$ to 43% while leaving all the other parameters at
their default values reduces the $\chi^{2}_{\textrm{pen}}$ to 21.8, showing
that the largest tension is due to the NonRes$1\pi$ parameter.

The $\pi$-iso dial is allowed to vary in the range 0–1 in the fit,
corresponding to a continuous variation between an RS angular distribution and
an isotropic distribution for $\Delta(1232)$ decay. To avoid the normalization
of the $\nu_{\mu}$CC$N\pi^{\pm}$ measurement pulling parameters in the
$\nu_{\mu}$CC1$\pi^{\pm}$ model, the NonRes$2\pi$ dial was allowed to vary
between 0-300% of the nominal value.

When varying one of the five hA pion FSI dials, GENIE automatically adjusts
the remaining parameters to preserve the total pion cross-section and maintain
agreement with pion-nucleus scattering data. This “cushion” technique
introduced instabilities in the $\chi^{2}$ surface, so it was not possible to
include multiple pion FSI parameters in a simultaneous fit. Instead we
performed fits with only one of the FSI parameters floating. No $\chi^{2}$
penalty terms were added for the FSI dial in either tuning: the parameters
were driven solely by MINERvA data. The charge exchange and pion production
dials had small contributions to the overall $\chi^{2}$ for the selected data,
forcing the parameters to be inflated beyond +3$\sigma$ of GENIE’s
recommendation, with large post-fit uncertainties. Furthermore, the pion
elastic scattering parameter is strongly constrained by external data, so its
1$\sigma$ variation has a small impact on the MINERvA distributions. The non-
FSI fit parameters’ (e.g. $M_{\mathrm{A}}^{\mathrm{res}}$) central values and
uncertainties all agreed for the five fits. Here we present the results from
the FrAbs and FrInel fits.

The NUISANCE interface to MINUIT2 James and Roos (1975) was used to perform
the fits. At each iteration, the GENIE-ReWeight package was used on an event-
by-event basis to update the MC predictions before the total $\chi^{2}$ was
calculated. The uncertainties in the fitted parameters were determined using
the HESSE routine in MINUIT2. The best fit results from the joint tuning are
shown in Table 5. Fig. 8 shows the ratios of the best fit prediction to the
data for all four kinematic variables of interest when the pion absorption FSI
parameter (FrAbs) is floated in the fit; Fig. 9 is the same, but when the
inelastic scattering FSI parameter (FrInel) is floated. Notably, the two FSI
fits are very similar in both minimum $\chi^{2}$ and best-fit parameter
values.

Comparing to the results of the ANL and BNL reanalysis, larger values of
$M_{\mathrm{A}}^{\mathrm{res}}$ and smaller values of $\mathrm{NormRes}$ were
found by the fit, pulling the parameters closer to GENIE nominal. The
NonRes$1\pi$ parameter is strongly bound by the bubble chamber data and the
MINERvA data did little to improve on this constraint. The penalty term
contributed to the $\chi^{2}$ by 9.3 for the FrAbs fit and 11.1 for the FrInel
fit. This is a significant improvement over the default, but indicates that
there is mild tension between the nucleon and nuclear data. The post-fit
correlation matrices are provided in Fig 10. The ANL/BNL input correlations
are largely maintained in our fit.

Parameter | Default Value | ANL/BNL Value | FrAbs Fit Result | FrInel Result
---|---|---|---|---
$M_{\mathrm{A}}^{\mathrm{res}}$ (GeV) | $1.12\pm 0.22$ | $0.94\pm 0.05$ | $1.07\pm 0.04$ | $1.08\pm 0.04$
$\mathrm{NormRes}$ (%) | $100\pm 20$ | $115\pm 7$ | $94\pm 6$ | $92\pm 6$
NonRes$1\pi$ (%) | $100\pm 50$ | $43\pm 4$ | $44\pm 4$ | $44\pm 4$
NonRes$2\pi$ (%) | $100\pm 50$ | - | $166\pm 32$ | $161\pm 33$
$\pi$-iso | 0 = RS | - | 1 = Iso (limit) | 1 = Iso (limit)
FrAbs (%) | 100 $\pm$ 30 | - | $109\pm 16$ | -
FrInel (%) | 100 $\pm$ 40 | - | - | $109\pm 24$
MINERvA $\chi^{2}$ | 275.6 | 312.7 | 242.3 | 240.7
$\chi^{2}_{\textrm{pen}}$ | 299.3 | 0.0 | 9.3 | 11.1
Total $\chi^{2}$ | 574.8 | 312.7 | 251.6 | 251.8
$\mathrm{N_{DoF}}$ | 148 | 148 | 145 | 145
Table 5: Fit results from tuning GENIE parameters in NUISANCE. The “ANL/BNL
Value” column shows the contributions when parameters are fixed at values of
Table 4.

Tables 6 and 7 show the results when individual MINERvA data were tuned in
separate fits. Since three of the four channels were removed in these fits the
constraint from data is weakened and the total $\chi^{2}$ is steered by the
bubble chamber $\chi^{2}$ penalty. The individual channel fits also found
values at the 300% limit for NonRes$2\pi$ dial, except in the
$\nu_{\mu}$CC$N\pi^{\pm}$ channel, where the result was unchanged by the fit.
Only the $\nu_{\mu}$CC$N\pi^{\pm}$ channel has a significant contribution from
nonresonant 2$\pi^{\pm}$ production. In the other fits, the parameter is
largely unconstrained and has little impact on the fitted distributions. The
$\chi^{2}$ per degree of freedom is indicative of a poor fit in the
$\nu_{\mu}$CC$N\pi^{\pm}$ and $\nu_{\mu}$CC$1\pi^{0}$ channels, but not in the
$\nu_{\mu}$CC1$\pi^{\pm}$ or $\bar{\nu}_{\mu}$CC$1\pi^{0}$ channels.
Furthermore, the $\nu_{\mu}$CC$1\pi^{0}$ shows the strongest $\chi^{2}$
penalty, indicating tension with the ANL/BNL prior. Given the different
kinematic regions covered by the channels (see Table 1) and the different
physics (e.g. fraction of coherent pion production) it is difficult to infer
what combination of effects are at work. Isotropic emission was preferred in
all fits, driven by the $\theta_{\pi}$ distributions. Disagreements in the
$\theta_{\pi}$ spectrum are clearly seen in the data/MC ratios of Fig. 8 and
9, and the large $\chi^{2}$ values observed for the $\nu_{\mu}$CC$N\pi^{\pm}$
and $\nu_{\mu}$CC$1\pi^{0}$ channels.

The individual $\chi^{2}$ contributions in the joint tuning best fit, shown in
sixth and seventh columns (“FrAbs Tune” and “FrInel Tune”) of Table 2, show
that not all distributions in all channels benefit from the model variations,
as the default GENIE fits have a better $\chi^{2}$ for some distributions. In
particular, the $\nu_{\mu}$CC$1\pi^{0}$ channel distributions have worse
agreement after the tuning, with only the $\theta_{\pi}$ distribution
improving in $\chi^{2}$, whereas all channels benefit from the shift to
isotropic emission. While there is an overall improvement over the ANL/BNL
tune when comparing the combined $\chi^{2}$ results, Figs. 8 and 9 show that
there are still unresolved shape disagreements in both the $T_{\pi}$ and
$\theta_{\mu}$ kinematics.

Parameter |  $\nu_{\mu}$CC1$\pi^{\pm}$ |  $\nu_{\mu}$CC$N\pi^{\pm}$ |  $\nu_{\mu}$CC$1\pi^{0}$ |  $\bar{\nu}_{\mu}$CC$1\pi^{0}$
---|---|---|---|---
$M_{\mathrm{A}}^{\mathrm{res}}$ (GeV) | 0.97 $\pm$ 0.05 | $0.97\pm 0.05$ | $1.02\pm 0.05$ | $0.96\pm 0.05$
$\mathrm{NormRes}$ (%) | 110 $\pm$ 7 | $110\pm 7$ | $104\pm 7$ | $111\pm 7$
NonRes$1\pi$ (%) | 43 $\pm$ 4 | $42\pm 4$ | $44\pm 4$ | $43\pm 4$
NonRes$2\pi$ (%) |  $300$ (limit) | 99$\pm$30 |  $300$ (limit) |  $300$ (limit)
$\pi$-iso | 1 = Iso (limit) | 1 = Iso (limit) | 1 = Iso (limit) | 1 = Iso (limit)
FrAbs (%) | 156 $\pm$ 53 | $128\pm 34$ | $126\pm 17$ | $82\pm 31$
MINERvA $\chi^{2}$ | 36.6 | 64.1 | 92.3 | 34.6
$\chi^{2}_{\textrm{pen}}$ | 0.5 | 0.7 | 3.2 | 0.3
Total $\chi^{2}$ | 37.1 | 64.8 | 95.5 | 34.9
$\mathrm{N_{DoF}}$ | 35 | 36 | 32 | 33
Table 6: Individual channel tuning results when the FrAbs dial is treated as the free FSI parameter.  Parameter |  $\nu_{\mu}$CC1$\pi^{\pm}$ |  $\nu_{\mu}$CC$N\pi^{\pm}$ |  $\nu_{\mu}$CC$1\pi^{0}$ |  $\bar{\nu}_{\mu}$CC$1\pi^{0}$
---|---|---|---|---
$M_{\mathrm{A}}^{\mathrm{res}}$ (GeV) | 0.97 $\pm$ 0.05 | $0.97\pm 0.05$ | $1.03\pm 0.05$ | $0.96\pm 0.05$
$\mathrm{NormRes}$ (%) | 109 $\pm$ 7 | $108\pm 7$ | $103\pm 7$ | $112\pm 7$
NonRes$1\pi$ (%) | 42 $\pm$ 4 | $42\pm 4$ | $43\pm 4$ | $43\pm 4$
NonRes$2\pi$ (%) |  $300$ (limit) | 110$\pm$30 |  $300$ (limit) |  $300$ (limit)
$\pi$-iso | 1 = Iso (limit) | 1 = Iso (limit) | 1 = Iso (limit) | 1 = Iso (limit)
FrInel (%) | 117 $\pm$ 54 | $127\pm 33$ | 0 (limit) | $80\pm 59$
MINERvA $\chi^{2}$ | 37.1 | 63.4 | 86.9 | 34.9
$\chi^{2}_{\textrm{pen}}$ | 0.7 | 1.3 | 3.4 | 0.2
Total $\chi^{2}$ | 37.8 | 64.7 | 90.3 | 35.1
$\mathrm{N_{DoF}}$ | 35 | 36 | 32 | 33
Table 7: Individual channel tuning results when the FrInel dial is treated as
the free FSI parameter.

Figure 8: MC/data ratios at the best fit point for the FrAbs joint tuning.

Figure 9: MC/data ratios at the best fit point for the FrInel joint tuning.

Figure 10: Correlation matrix from from tuning GENIE parameters in NUISANCE
with FrAbs included as a fit parameter (left) and with FrInel included as a
fit parameter (right).

The tension between MINERvA’s nuclear data and the constraints from ANL and
BNL nucleon data is difficult to confidently pinpoint; the lack of lepton mass
effects Berger and Sehgal (2007), modification to the resonance propagator in
the nucleus Singh  _et al._ (1998); Oset and Salcedo (1987), missing diagrams
describing the non-resonant background contributions Kabirnezhad (2018),
dynamical coupled channels Nakamura  _et al._ (2017), interactions on
correlated initial states, and the pion FSI model Mosel and Gallmeister (2017)
are all part of an incomplete list of possible culprits.

##  VI Ad hoc $Q^{2}$ Supression

Further modifications beyond the standard GENIE dials are required to resolve
the observed tensions. Fig. 11 (not used for any tuning) shows the $Q^{2}$
distributions observed at MINERvA in for our tunes. The data is below the
predictions of the tunes of Section V at low values of $Q^{2}$. There are also
also differences at low $\theta_{\mu}$, as shown in Fig. 8 and 9. Measurements
of $\nu_{\mu}$CC1$\pi^{\pm}$ and $\nu_{\mu}$CC$1\pi^{0}$ interactions on
mineral oil at MiniBooNE have shown a data/MC shape discrepancy for the RS
implementation in the NUANCE model Casper (2002); Przewlocki (2009) in both
$Q^{2}$ and $\cos\theta_{\mu}$ distributions Aguilar-Arevalo  _et al._ (2011);
Aquilar-Arevalo  _et al._ (2011). In the MINOS quasielastic analysis Adamson
_et al._ (2015) on iron, which used NeuGen Gallagher (2002), a similar
disagreement was observed when studying pion production dominated sidebands.
Indeed, concerns about low-$Q^{2}$ modeling date back almost a decade Graczyk
and Sobczyk (2008). The data from MINOS and MiniBooNE experiments and the
MINERvA data on CH studied herein suggest that the RS implementation common to
each of the generators needs to be suppressed at low $Q^{2}$. Collective
effects, which are usually modeled in the random phase approximation, are
known to affect the $Q^{2}$ distribution of neutrino-nucleus reactions at low
$Q^{2}$. Motivated by these considerations, we attempted to improve the
$\theta_{\mu}$ modeling by introducing a $Q^{2}$-dependent correction to the
model.

The MINOS collaboration suppression was expressed as

$\displaystyle R=\frac{A}{1+\exp\\{1-\sqrt{Q^{2}}/Q_{0}\\}},$ (5)

where the free parameters $A=1.010$ and $Q_{0}=0.156~{}\mbox{GeV}$ were
empirically extracted from bin-by-bin fits in $Q^{2}$ to the data, and a hard
cut-off at $Q^{2}<0.7~{}\mbox{GeV}^{2}$ was imposed.

We chose an empirical function so that the shape of the suppression preferred
by each of the MINERvA channels could be extracted. The empirical correction
function is applied to events with a resonance decay inside the nucleus giving
rise to a pion. Our suppression term is defined by choosing 3 points
$(x_{i},R_{i})_{i=1,2,3}$ between $0.0<x<1.0$ and $0.0<R<1.0$, where $x\equiv
Q^{2}$. Motivated the ANL/BNL curves in Fig. 11, the correction is assumed to
approach unity as $Q^{2}$ approaches $0.7~{}\mbox{GeV}^{2}$, providing the
constraint $(x_{3},R_{3})=(0.7~{}\mbox{GeV}^{2},1.0)$. Lagrange interpolation
is used to derive a curvature from $R_{2}$ by assuming a simple interpolation
between the points $(x_{1},0.0)$, $(x_{2},R_{2})$, and
$(0.7~{}\mbox{GeV}^{2},1.0)$:

$\displaystyle R(Q^{2}<x_{3})$
$\displaystyle=\frac{R_{2}(Q^{2}-x_{1})(Q^{2}-x_{3})}{(x_{2}-x_{1})(x_{2}-x_{3})}$
$\displaystyle+\frac{(Q^{2}-x_{1})(Q^{2}-x_{2})}{(x_{3}-x_{1})(x_{3}-x_{2})}.$
(6)

This interpolation function is then used to calculate the correction for each
event as

$\displaystyle w(Q^{2})=1-(1-R_{1})(1-R(Q^{2}))^{2}.$ (7)

where $R_{1}$ defines the magnitude of the correction function at the
intercept, $x_{1}=0.0$. $x_{2}$ is chosen to be $Q^{2}=0.35~{}\mbox{GeV}^{2}$
so that $R_{2}$ describes the curvature at the centre point of the correction.
Expressing the weights with Equations 6 and 7 ensures that the magnitude at
$x_{2}$ always lies between $R_{1}$ and 1.0, avoiding parameter sets with
large unphysical peaks in the correction function. Additionally, the squared
term in Equation 7 ensures that $w(Q^{2})\rightarrow 1.0$ as $x\rightarrow
x_{3}$, avoiding discontinuous steps in the weighting function at $x_{3}$. The
fitted parameters $R_{1}$ and $R_{2}$ were limited to $0.0<R_{1}<1.0$ and
$0.5<R_{2}<1.0$ to avoid extraneous solutions, e.g. double peaks.

Figure 11: Comparisons of the nominal and tuned models to MINERvA
$\nu_{\mu}$CC1$\pi^{\pm}$ (left top), $\nu_{\mu}$CC$N\pi^{\pm}$ (right top),
$\nu_{\mu}$CC$1\pi^{0}$ (left bottom) and $\bar{\nu}_{\mu}$CC$1\pi^{0}$ (right
bottom) distributions in $Q^{2}$. The $\chi^{2}$ is computed using the full
covariance matrices. The distributions were not explicitly used in the tuning
procedure.

The fit results are shown in Table 8. The correction from the fit with FrAbs
taken as a free parameter are compared to the MINOS low-$Q^{2}$ correction in
Fig. 12. Our fits obtain a suppression factor that is similar to the MINOS
one, with almost identical suppression at $Q^{2}$ = 0, albeit with less
curvature, particularly in the $\nu_{\mu}$CC1$\pi^{\pm}$ and
$\nu_{\mu}$CC$N\pi^{\pm}$ channels. The correction factors from the fit with
FrInel or FrAbs as free parameters give similar results.

The correlation matrices for the fits including a $Q^{2}$ dependent
suppression are provided in Fig 13. Again, the ANL/BNL input prior covariance
is maintained. The parameters largely correlate in the same way for the FrAbs
and FrInel fit, and for the FrInel fit the $R_{1}$ and $R_{2}$ parameters are
negatively correlated.

Parameter | FrAbs Tune | FrAbs + low-$Q^{2}$ Tune | FrInel Tune | FrInel + low-$Q^{2}$ Tune
---|---|---|---|---
$M_{\mathrm{A}}^{\mathrm{res}}$ $(GeV)$ | $1.07\pm 0.04$ | $0.92\pm 0.02$ | $1.08\pm 0.04$ | $0.93\pm 0.05$
$\mathrm{NormRes}$ (%) | $94\pm 6$ | $116\pm 3$ | $92\pm 6$ | $116\pm 7$
NonRes$1\pi$ (%) | $43\pm 4$ | $46\pm 4$ | $44\pm 4$ | $46\pm 4$
NonRes$2\pi$ (%) | $166\pm 32$ | 99$\pm$31 | $161\pm 33$ | $120\pm 32$
$\pi$-iso |  $1.0$ (limit) | 1.0 (limit) | 1.0 (limit) | 1.0 (limit)
FrAbs (%) | $109\pm 16$ | $48\pm 21$ | - | -
FrInel (%) | - | - | $109\pm 24$ | $132\pm 27$
Lag. $R_{1}$ | - | $0.32\pm 0.06$ | - | $0.37\pm 0.09$
Lag. $R_{2}$ | - |  $0.5$ (limit) | - | $0.60\pm 0.16$
MINERvA $\chi^{2}$ | 242.3 | 212.2 | 240.7 | 215.7
$\chi^{2}_{\textrm{pen}}$ | 9.3 | 0.7 | 11.1 | 0.5
Total$\chi^{2}$ | 251.6 | 212.9 | 251.8 | 216.2
$\mathrm{N_{DoF}}$ | 145 | 143 | 145 | 143
Table 8: Ad hoc low-$Q^{2}$ suppression model tuning results compared to the
tuning results without the low-$Q^{2}$ suppression.

Figure 12: Extracted low-$Q^{2}$ suppression factors from the FrAbs +
low-$Q^{2}$ tuning to each channel. The left and right plots compare the
results for the charged and neutral pion production channels respectively.
Shown in red is the uncertainty band extracted from the joint fit to all 4
channels simultaneously.

Figure 13: Correlation matrix from tuning GENIE parameters with an ad hoc
low-$Q^{2}$ supression with FrAbs included as a fit parameter (left) and with
FrInel included as a fit parameter (right).

Figure 14 (Fig. 15) shows the ratio of the resulting fits to the MINERvA data
when FrAbs (FrInel) is taken as a free parameter. As anticipated, the
predictions now have better agreement with the data in regards to the
$\theta_{\mu}$ distribution, and the $\chi^{2}$ values are improved by the
introduction of our ad hoc low-$Q^{2}$ correction. Other fit parameters are
for the most part unchanged by the introduction of the low-$Q^{2}$ correction.
Furthermore, $M_{\mathrm{A}}^{\mathrm{res}}$ and $\mathrm{NormRes}$ are closer
to their values when fitting ANL and BNL data, indicating the $Q^{2}$
correction alleviates the tension between nucleon and nuclear modeling. Fig.
11 shows the comparison of all our models directly against MINERvA data in
$Q^{2}$. Although the tuning sees improvement in the $\chi^{2}$ for the
$\nu_{\mu}$CC$1\pi^{0}$ and $\bar{\nu}_{\mu}$CC$1\pi^{0}$ distributions, the
$\nu_{\mu}$CC1$\pi^{\pm}$ and $\nu_{\mu}$CC$N\pi^{\pm}$ distributions get
worse, hinting at tensions in the charged and neutral pion production
channels.

Figure 14: MC/data ratios at the best fit points from the FrAbs tuning with
low-$Q^{2}$ suppression included.

Figure 15: MC/data ratios at the best fit points from the FrInel tuning with
low-$Q^{2}$ suppression included.

Tables 9 and 10 show the results of the fits to individual channels, and Table
11 shows the breakdown of contributions to the $\chi^{2}$ from the individual
channels. The best fit $\chi^{2}$ value was significantly improved for each
channel tuning when using a low-$Q^{2}$ suppression and the extracted
parameters were consistent with the ANL/BNL tunings. Pion kinematic
distributions are not improved, and in some cases are slightly worse, as a
result of including the low-$Q^{2}$ suppression. It is clear from Table 8 (or
by comparing Tables 6 and 9) that the low-$Q^{2}$ suppression has a similar
effect in the fit to the FrAbs parameter. When the low-$Q^{2}$ suppression is
introduced, FrAbs tends to consistently lower values. It is also clear that
the $\bar{\nu}_{\mu}$CC$1\pi^{0}$ channel favors stronger low-$Q^{2}$
suppression than the other channels.

Parameter |  $\nu_{\mu}$CC1$\pi^{\pm}$ |  $\nu_{\mu}$CC$N\pi^{\pm}$ |  $\nu_{\mu}$CC$1\pi^{0}$ |  $\bar{\nu}_{\mu}$CC$1\pi^{0}$
---|---|---|---|---
$M_{\mathrm{A}}^{\mathrm{res}}$ (GeV) | $0.93\pm 0.02$ | $0.92\pm 0.02$ | $0.96\pm 0.05$ | $0.94\pm 0.05$
$\mathrm{NormRes}$ (%) | $115\pm 3$ | $117\pm 3$ | $114\pm 7$ | $115\pm 7$
NonRes$1\pi$ (%) | $43\pm 4$ | $43\pm 4$ | $45\pm 4$ | $43\pm 4$
NonRes$2\pi$ (%) |  $300$ (limit) | $70\pm 28$ |  $300$ (limit) |  $300$ (limit)
$\pi$-iso | 1 = Iso (limit) | 1 = Iso (limit) | 1 = Iso (limit) | 1 = Iso (limit)
FrAbs (%) | 92 $\pm$ 65 | $79\pm 40$ | $74\pm 22$ | $34\pm 35$
Lag. $R_{1}$ | $0.53\pm 0.16$ | $0.43\pm 0.13$ | $0.21\pm 0.14$ | $0.14\pm 0.22$
Lag. $R_{2}$ |  $0.50$ (limit) |  $0.50$ (limit) | $0.63\pm 0.31$ |  $1.00$ (limit)
MINERvA $\chi^{2}$ | 32.2 | 55.7 | 71.2 | 27.7
$\chi^{2}_{\textrm{pen}}$ | 0.1 | 0.4 | 0.5 | 0.0
Total $\chi^{2}$ | 32.3 | 56.1 | 71.7 | 27.7
$\mathrm{N_{DoF}}$ | 33 | 34 | 30 | 31
Table 9: Individual channel FrAbs + low-$Q^{2}$ tuning results.  Parameter |  $\nu_{\mu}$CC1$\pi^{\pm}$ |  $\nu_{\mu}$CC$N\pi^{\pm}$ |  $\nu_{\mu}$CC$1\pi^{0}$ |  $\bar{\nu}_{\mu}$CC$1\pi^{0}$
---|---|---|---|---
$M_{\mathrm{A}}^{\mathrm{res}}$ (GeV) | $0.93\pm 0.02$ | $0.91\pm 0.02$ | $0.95\pm 0.05$ | $0.94\pm 0.05$
$\mathrm{NormRes}$ (%) | $116\pm 3$ | $117\pm 3$ | $114\pm 7$ | $115\pm 7$
NonRes$1\pi$ (%) | $43\pm 4$ | $43\pm 4$ | $44\pm 4$ | $43\pm 4$
NonRes$2\pi$ (%) |  $300$ (limit) | $78\pm 28$ |  $300$ (limit) | 300 (limit)
$\pi$-iso | 1 = Iso (limit) | 1 = Iso (limit) | 1 = Iso (limit) | 1 = Iso (limit)
FrInel (%) | $179\pm 63$ | $173\pm 37$ | $8\pm 125$ | $103\pm 57$
Lag. $R_{1}$ | $0.49\pm 0.14$ | $0.38\pm 0.13$ | $0.25\pm 0.17$ | $0.31\pm 0.26$
Lag. $R_{2}$ |  $0.50$ (limit) |  $0.50$ (limit) | $0.76\pm 0.37$ |  $1.00$ (limit)
MINERvA $\chi^{2}$ | 30.8 | 52.1 | 69.5 | 30.9
$\chi^{2}_{\textrm{pen}}$ | 0.1 | 0.6 | 0.2 | 0.0
Total $\chi^{2}$ | 30.9 | 52.7 | 69.7 | 30.9
$\mathrm{N_{DoF}}$ | 33 | 34 | 30 | 31
Table 10: Individual channel FrInel + low-$Q^{2}$ tuning results.  Distribution | Channel | $\mathrm{N_{bins}}$ | FrAbs Tune |  FrAbs + low-$Q^{2}$ Tune | FrInel Tune |  FrInel + low-$Q^{2}$ Tune
---|---|---|---|---|---|---
$p_{\mu}$ (Rate) |  $\nu_{\mu}$CC1$\pi^{\pm}$ | 8 | 12.0 | 10.8 | 12.3 | 10.9
|  $\nu_{\mu}$CC$N\pi^{\pm}$ | 9 | 26.1 | 16.2 | 26.8 | 17.9
|  $\nu_{\mu}$CC$1\pi^{0}$ | 8 | 19.0 | 26.2 | 19.3 | 26.9
|  $\bar{\nu}_{\mu}$CC$1\pi^{0}$ | 9 | 6.2 | 7.1 | 6.3 | 7.2
$\theta_{\mu}$ (Shape) |  $\nu_{\mu}$CC1$\pi^{\pm}$ | 9 | 7.5 | 7.4 | 7.4 | 7.1
|  $\nu_{\mu}$CC$N\pi^{\pm}$ | 9 | 4.0 | 6.3 | 4.1 | 5.6
|  $\nu_{\mu}$CC$1\pi^{0}$ | 9 | 44.5 | 20.0 | 45.6 | 20.5
|  $\bar{\nu}_{\mu}$CC$1\pi^{0}$ | 9 | 10.2 | 7.0 | 10.3 | 6.9
$T_{\pi}$ (Shape) |  $\nu_{\mu}$CC1$\pi^{\pm}$ | 7 | 2.5 | 2.5 | 2.3 | 2.4
|  $\nu_{\mu}$CC$N\pi^{\pm}$ | 7 | 31.2 | 28.9 | 29.4 | 27.7
|  $\nu_{\mu}$CC$1\pi^{0}$ | 7 | 30.9 | 27.1 | 29.9 | 32.0
|  $\bar{\nu}_{\mu}$CC$1\pi^{0}$ | 7 | 16.6 | 15.7 | 16.0 | 18.7
$\theta_{\pi}$ (Shape) |  $\nu_{\mu}$CC1$\pi^{\pm}$ | 14 | 13.0 | 13.4 | 12.6 | 12.6
|  $\nu_{\mu}$CC$N\pi^{\pm}$ | 14 | 6.9 | 7.0 | 6.2 | 6.3
|  $\nu_{\mu}$CC$1\pi^{0}$ | 11 | 8.3 | 12.2 | 8.9 | 9.4
|  $\bar{\nu}_{\mu}$CC$1\pi^{0}$ | 11 | 3.4 | 4.4 | 3.5 | 3.7
Total $\chi^{2}$ |  | 148 | 242.3 | 212.2 | 240.7 | 215.7
Table 11: Channel by channel contributions to the $\chi^{2}$ at for the GENIE
tunings with and without the low-$Q^{2}$ correction included.

##  VII Concluding Remarks

We have adjusted the parameters of the GENIE model that are important for pion
production to match MINERvA data in the $\nu_{\mu}$CC1$\pi^{\pm}$,
$\nu_{\mu}$CC$N\pi^{\pm}$, $\nu_{\mu}$CC$1\pi^{0}$ and
$\bar{\nu}_{\mu}$CC$1\pi^{0}$ channels, using the NUISANCE framework. We
incorporate existing results which informs the GENIE model using ANL and BNL
bubble chamber data from scattering off protons and deuterons. Fits of
selected GENIE model parameters were done using the kinematic distributions
$p_{\mu}$, $\theta_{\mu}$, $T_{\pi}$ and $\theta_{\pi}$. Parameter fits were
performed with either the fraction of pions absorbed or the fraction of pions
inelastically scattered in FSI as a floating parameter, with broadly similar
conclusions for the two cases.

The results of the fit (see Table 5) show that the tuning improves the GENIE
pion production model significantly, but tensions remain. The pull on the
ANL/BNL prior demonstrates a tension between MINERvA nuclear target data and
the light-target bubble-chamber data sets used to make the prior, indicating a
deficiency in the GENIE nuclear model which cannot be fixed by modifying the
available reweighting dials. Additionally, fitting to individual MINERvA pion
production channels produces different best-fit parameters, demonstrating that
GENIE cannot describe the different exclusive channels in a consistent manner
with the available dials (shown in Tables 6 and 7). Because the four channels
cover different kinematic regions (see Table 1) and contain different physics
(e.g. different coherent pion production contributions or nonresonant
processes), it is difficult to pinpoint the origin of the discrepancy between
the model and the different MINERvA data sets.

Following experimental hints of discrepancies at low-$Q^{2}$ for a variety of
cross-section measurements on nuclear targets, an additional empirical
low-$Q^{2}$ suppression was introduced and the fits were repeated. Although
the data showed a preference for a strong suppression at low-$Q^{2}$ and the
agreement improved for $\theta_{\mu}$ and $Q^{2}$ distributions, tensions
remain. In particular, fits to individual MINERvA channels still produced
different results, and favor different parameter values for the low-$Q^{2}$
suppression.

The main conclusion of this work is that current neutrino experiments
operating in the few–GeV region should think critically about single pion
production models and uncertainties, as the Monte Carlo models which are
currently widely used in the field are unable to explain multiple data sets,
even when they are from a single experiment.

A key strength of this analysis is its development within the NUISANCE
framework, allowing it to be easily repeated with alternate model assumptions,
neutrino interaction generators, and different data. The developments
presented here will be used in future iterations of this work, as the MINERvA
collaboration works towards a GENIE model that provides a good description of
all their available data, and can be easily applied to other measurements and
experiments.

## Acknowledgements

P.S., L.P., and C.V.C.W. would like to thank the UK Science and Technology
Facilities Council (STFC) for Ph.D. funding support. P.S. acknowledges the
Fermilab Neutrino Physics Center for the scholarship that funded this work and
thanks Fermilab and the MINERvA collaboration for their hospitality during
this work. C.W. acknowledges the support of the Swiss National Science
Foundation and SERI.

This document was prepared by members of the MINERvA Collaboration using the
resources of the Fermi National Accelerator Laboratory (Fermilab), a U.S.
Department of Energy, Office of Science, HEP User Facility. Fermilab is
managed by Fermi Research Alliance, LLC (FRA), acting under Contract No. DE-
AC02-07CH11359. These resources included support for the MINERvA construction
project, and support for construction also was granted by the United States
National Science Foundation under Award No. PHY-0619727 and by the University
of Rochester. Support for participating scientists was provided by NSF and DOE
(USA); by CAPES and CNPq (Brazil); by CoNaCyT (Mexico); by Proyecto Basal FB
0821, CONICYT PIA ACT1413, Fondecyt 3170845 and 11130133 (Chile); by CONCYTEC,
DGI-PUCP, and IDI/IGI-UNI (Peru); and by the Latin American Center for Physics
(CLAF); NCN Opus Grant No. 2016/21/B/ST2/01092 (Poland). We thank the MINOS
Collaboration for use of its near detector data. Finally, we thank the staff
of Fermilab for support of the beam line, the detector, and computing
infrastructure.

## References

  * Mosel and Gallmeister (2017) U. Mosel and K. Gallmeister, Phys. Rev. C 96 (2017).
  * Alvarez-Ruso  _et al._ (2018) L. Alvarez-Ruso _et al._ , Prog. Part. Nucl. Phys. 100 (2018).
  * R. González-Jiménez and Jachowicz (2018) K. N. R. González-Jiménez and N. Jachowicz, Phys. Rev. D 97 (2018).
  * Wilkinson  _et al._ (2016) C. Wilkinson _et al._ , Phys. Rev. D 93 (2016).
  * Mahn  _et al._ (2018) K. Mahn, C. Marshall, and C. Wilkinson, Ann. Rev. Nucl. Part. Sci. 68 (2018).
  * Stowell  _et al._ (2017) P. Stowell _et al._ , J. Instrum. 12 (2017).
  * Andreopoulos  _et al._ (2009) C. Andreopoulos _et al._ , Nucl. Instrum. Methods Phys. Res. 614 (2009).
  * Andreopoulos  _et al._ (2015) C. Andreopoulos _et al._ , arXiv:1510.05494 (2015).
  * Kabirnezhad (2018) M. Kabirnezhad, Phys. Rev. D 97 (2018).
  * Leitner  _et al._ (2009) T. Leitner, O. Buss, L. Alvarez-Ruso, and U. Mosel, Phys. Rev. C 79 (2009).
  * Nieves  _et al._ (2011) J. Nieves, I. R. Simo, and M. J. V. Vacas, Phys. Rev. C 83 (2011).
  * Hernadez and Nieves (2017) E. Hernadez and J. Nieves, Phys. Rev. D 95 (2017).
  * Rodrigues  _et al._ (2016) P. Rodrigues, C. Wilkinson, and K. McFarland, Eur. Phys. J. C 76 (2016).
  * Wilkinson  _et al._ (2014) C. Wilkinson, P. Rodrigues, S. Cartwright, L. Thompson, and K. McFarland, Phys. Rev. D 90 (2014).
  * Adamson  _et al._ (2015) P. Adamson _et al._ (MINOS), Phys. Rev. D 91 (2015).
  * Aguilar-Arevalo  _et al._ (2011) A. A. Aguilar-Arevalo _et al._ (MiniBooNE), Phys. Rev. D 83 (2011).
  * Eberly  _et al._ (2015) B. Eberly _et al._ (MINERvA), Phys. Rev. D 92 (2015).
  * McGivern  _et al._ (2016) C. L. McGivern _et al._ (MINERvA), Phys. Rev. D 94 (2016).
  * Altinok  _et al._ (2017) O. Altinok _et al._ (MINERvA), Phys. Rev. D 96 (2017), with errata of summer 2019.
  * Aliaga  _et al._ (2014) L. Aliaga _et al._ (MINERvA), Nucl. Instrum. Methods Phys. Res., Sect. A 743 (2014).
  * Katori (2008) T. Katori, Ph. D. thesis, University of Indiana (2008).
  * Smith and Moniz (1972) R. A. Smith and E. J. Moniz, Nucl. Phys. B 43 (1972).
  * Bodek and Ritchie (1981) A. Bodek and J. Ritchie, Phys. Rev. D 23 (1981).
  * Nieves  _et al._ (2004) J. Nieves, J. E. Amaro, and M. Valverde, Phys. Rev. C 70 (2004).
  * Gran  _et al._ (2013) R. Gran, J. Nieves, F. Sanchez, and M. J. V. Vacas, Phys. Rev. D 88 (2013).
  * Patrick  _et al._ (2018) C. E. Patrick _et al._ (MINERvA), Phys. Rev. D 97 (2018).
  * Mislivec  _et al._ (2018) A. Mislivec _et al._ (MINERvA), Phys. Rev. D 97 (2018).
  * Aliaga  _et al._ (2016) L. Aliaga _et al._ (MINERvA), Phys. Rev. D 94 (2016).
  * Andreopoulos  _et al._ (2017) C. Andreopoulos _et al._ (GENIE), ”Pre-generated public splines for GENIE 2.12.6”, https://www.hepforge.org/archive/genie/data/2.12.6/ (2017).
  * Rein and Sehgal (1981) D. Rein and L. M. Sehgal, Annals Phys. 133 (1981).
  * Graczyk and Sobczyk (2008) K. M. Graczyk and J. T. Sobczyk, Phys. Rev. D 77 (2008).
  * Bodek and Yang (2003) A. Bodek and U. K. Yang, J. Phys. G 29 (2003).
  * Yang  _et al._ (2009) T. Yang, C. Andreopoulos, H. Gallagher, K. Hoffmann, and P. Kehayias, Eur. Phys. J. C 63 (2009).
  * Koba  _et al._ (1972) Z. Koba, H. B. Nielsen, and P. Olesen, Nucl. Phys. B 40 (1972).
  * T. Sjostrand and Skands (2008) S. M. T. Sjostrand and P. Z. Skands, Comput. Phys. Commun. 178 (2008).
  * Rein and Sehgal (1983) D. Rein and L. M. Sehgal, Nucl. Phys. B 223 (1983).
  * Rein and Sehgal (2007) D. Rein and L. M. Sehgal, Phys. Lett. B 657 (2007).
  * Berger and Sehgal (2009) C. Berger and L. M. Sehgal, Phys. Rev. D 79 (2009).
  * Dytman (2009) S. Dytman, Acta Phys. Polon. B 40 (2009).
  * Arndt  _et al._ (2006) R. A. Arndt, W. J. Briscoe, I. I. Strakovsky, and R. L. Workman, Phys. Rev. C 74 (2006).
  * Buss  _et al._ (2012) O. Buss _et al._ , Phys. Rept. 512 (2012).
  * Radecky and others (ANL) G. M. Radecky and others (ANL), Phys. Rev. D 25 (1982).
  * Kitagaki and others (BNL) T. Kitagaki and others (BNL), Phys. Rev. D 34 (1986).
  * Barish and others (ANL) S. J. Barish and others (ANL), Phys. Rev. D 19 (1979).
  * Kitagaki and others (BNL) T. Kitagaki and others (BNL), Phys. Rev. D 42 (1990).
  * Adler (1968) S. L. Adler, Annals of Physics 50 (1968).
  * Allen  _et al._ (1980) P. Allen _et al._ (Aachen-Bonn-CERN-Munchen-Oxford), Nuclear Physics B 176 (1980).
  * Allasia  _et al._ (1983) D. Allasia _et al._ (Amsterdam-Bologna-Padova-Pisa-Saclay-Torino), Z. Phys. C 20 (1983).
  * Bell.  _et al._ (1978) J. Bell. _et al._ , Phys. Rev. Lett. 41 (1978).
  * James and Roos (1975) F. James and M. Roos, Comput. Phys. Commun. 10 (1975).
  * Berger and Sehgal (2007) C. Berger and L. M. Sehgal, Phys. Rev. D 76 (2007).
  * Singh  _et al._ (1998) S. K. Singh _et al._ , Physics Letters B 416 (1998).
  * Oset and Salcedo (1987) E. Oset and L. L. Salcedo, Nuclear Physics A 468 (1987).
  * Nakamura  _et al._ (2017) S. X. Nakamura _et al._ , Reports on Progress in Physics 80 (2017).
  * Casper (2002) D. Casper, Nucl. Phys. Proc. Suppl. 112 (2002).
  * Przewlocki (2009) P. Przewlocki, Acta Phys. Polon. B 40 (2009).
  * Aquilar-Arevalo  _et al._ (2011) A. A. Aquilar-Arevalo _et al._ (MiniBooNE), Phys. Rev. D 83 (2011).
  * Gallagher (2002) H. Gallagher, Nucl. Phys. Proc. Suppl. 112 (2002).


# The BLUES function method for second-order partial differential equations:
application to a nonlinear telegrapher equation

Jonas Berx<EMAIL_ADDRESS>Joseph O. Indekeu
<EMAIL_ADDRESS>Institute for Theoretical Physics, KU Leuven,
B-3001 Leuven, Belgium

(August 27, 2024)

###### Abstract

An analytic iteration sequence based on the extension of the BLUES (Beyond
Linear Use of Equation Superposition) function method to partial differential
equations (PDEs) with second-order time derivatives is studied. The original
formulation of the BLUES method is modified by introducing a matrix formalism
that takes into account the initial conditions for higher-order time
derivatives. The initial conditions of both the solution and its derivatives
now play the role of a source vector. The method is tested on a nonlinear
telegrapher equation, which can be reduced to a nonlinear wave equation by a
suitable choice of parameters. In addition, a comparison is made with three
other methods: the Adomian decomposition method, the variational iteration
method (with Green function) and the homotopy perturbation method. The matrix
BLUES function method is shown to be a worthwhile alternative for the other
methods.

###### keywords:

BLUES function method , analytic iteration , telegrapher equation

###### MSC:

[2020] 35A35 , 35G20

††journal: Partial Differential Equations in Applied Mathematics

##  1 Introduction

When studying nonlinear differential equations (DEs), it is often necessary to
construct approximate analytical solutions to the equation at hand because an
exact expression is not available. Hence, techniques such as the Adomian
decomposition method (ADM) [1, 2], variational iteration method (VIM) [3],
variational iteration method with Green function (GVIM) [4], homotopy
perturbation method (HPM) [5] or homotopy analysis method (HAM) [6] have been
developed and have been successfully applied to a wide range of problems in
physics and other sciences.

Recently, it has been demonstrated that the theory of Green functions can be
usefully extended to inhomogeneous nonlinear ordinary differential equations
(ODEs), effectively using superposition beyond the linear domain [7, 8, 9].
Consequently, we probed the usefulness of the BLUES function method in the
arena of nonlinear partial differential equations (PDEs) with a first-order
time derivative where the initial condition plays the role of the
inhomogeneous source by multiplication with a Dirac delta point source located
at time $t=0$ [10]. It was shown that the method produces globally convergent
approximants with high accuracy for systems of coupled nonlinear first-order
ODEs when the fixed points of the system are introduced into the linear
operator [11]. However, for second-order time derivatives, the BLUES function
method needs to be modified. In this work we will show that the combination of
the BLUES function method for PDEs with a matrix formalism is required to
generate approximate solutions to PDEs with second- (or higher-) order time
derivatives.

The setup of this paper is as follows. In Section 2 we extend the BLUES
function method to $n$th-order time derivatives and introduce the required
matrix formalism. Next, in Section 3, we apply the matrix BLUES function
method to a nonlinear telegrapher equation and compare the approximants with
those of the VIM, GVIM, ADM and HPM. Finally, in Section 4 we present some
conclusions and thoughts on future research.

##  2 BLUES function method for $n$th-order time derivatives

For nonlinear ODEs with source or sink terms, the BLUES function method has
been studied intensively [7, 8, 9, 11] and anticipated in [12]. For PDEs with
a first-order time derivative, the role of the source was taken up by the
initial condition multiplied by a Dirac point source located at $t=0$ [10].
For $n$th-order time derivatives this approach can be extended to include the
initial conditions of the derivatives of the solution. We propose the
following approach: the $n$th-order in time nonlinear PDE

$\widetilde{\mathcal{N}}^{[n]}_{t,x}\,u=\chi(x,t)$ (1)

with source $\chi(x,t)$ can be decomposed into $n$ first-order coupled PDEs
where the initial conditions for the solution and the derivatives can be
included as sources in their respective constituent equations by a suitable
multiplication with a point source at $t=0$. The BLUES method can subsequently
be applied to the system of first-order in time equations.

For equation (1), this system can be written as a different nonlinear operator
$\mathcal{N}_{t,x}\,$ acting on a vector of solutions
$\boldsymbol{U}=(U_{1},U_{2},...,U_{n})$, i.e.,

$\mathcal{N}_{t,x}\,\boldsymbol{U}(x,t)=\boldsymbol{\psi}(x,t),\;\forall t>0$
(2)

wherein $U_{1}=u$, $U_{2}=\partial u/\partial t$, …,
$U_{n}=\partial^{n}u/\partial t^{n}$. The initial conditions for $t=0$ are
collected in the vector $\boldsymbol{C}$,

$\boldsymbol{U}(x,0)=\boldsymbol{C}(x)\,.$ (3)

Now, $\boldsymbol{\psi}(x,t)=(0,...,0,\chi(x,t))$ is a vector of length $n$
that contains the external source in the last entry and zeroes everywhere
else. We now judiciously decompose the nonlinear operator
$\mathcal{N}_{t,x}\,$ into a linear operator $\mathcal{L}_{t,x}\,$, which
contains time derivatives of at most first order as a consequence of the
decomposition, and a residual operator $\mathcal{R}_{x}\,$, i.e.,
$\mathcal{R}_{x}\,\equiv\mathcal{L}_{t,x}\,-\mathcal{N}_{t,x}\,$, which
contains at least the nonlinear part of $\mathcal{N}_{t,x}\,$ and which has no
explicit time dependence. Thus the action of the linear operator on
$\boldsymbol{U}$ results in the following associated linear coupled system

$\mathcal{L}_{t,x}\,\boldsymbol{U}=\partialderivative{\boldsymbol{U}}{t}-A\boldsymbol{U}=\boldsymbol{\psi},\;\forall
t>0\,,$ (4)

where $A$ is an $n\crossproduct n$ matrix. We now rewrite the system of DEs in
an equivalent form by incorporating the initial condition through
multiplication of $\boldsymbol{C}$ with a Dirac delta source $\delta(t)$ and
including this term on the right-hand-side of the inhomogeneous system, i.e.,

$\mathcal{L}_{t,x}\,\boldsymbol{U}(x,t)=\frac{\partial\boldsymbol{U}}{\partial
t}(x,t)-A\boldsymbol{U}(x,t)=\boldsymbol{\psi}(x,t)\Theta(t)+\boldsymbol{C}(x)\delta(t)\equiv\boldsymbol{\varphi}(x,t),\;\forall
t\geq 0\,,$ (5)

where we have combined the external source $\boldsymbol{\psi}\Theta$ and the
“initial condition source” $\boldsymbol{C}\delta$ into the combined source
$\boldsymbol{\varphi}$.

The solution of this linear system (5) is the following two-variable
convolution

$\begin{split}\boldsymbol{U}(x,t)&=(G\ast\boldsymbol{\varphi})(x,t)\\\
&=\int_{0^{-}}^{t}\int_{\mathbb{R}}G(x-x^{\prime},t-t^{\prime})\left[\boldsymbol{\psi}(x^{\prime},t^{\prime})\Theta(t^{\prime})+\boldsymbol{C}(x^{\prime})\delta(t^{\prime})\right]\mathrm{d}x^{\prime}\mathrm{d}t^{\prime},\;\;\forall
t\geq 0\end{split}$ (6)

where $G(x,t)$ is the Green function matrix for the inhomogeneous linear
system. This $G(x,t)$ satisfies

$\mathcal{L}_{t,x}\,G(x,t)=0,\;\mbox{for}\;t>0,$ (7)

with Dirac-delta initial condition

${\rm lim}_{t\rightarrow 0}\,G(x,t)=\delta(x)\mathbb{1}\,.$ (8)

Adopting the BLUES function strategy, a solution to the nonlinear system (2)
is now proposed in the form of a convolution
$\boldsymbol{U}(x,t)=(B\ast\boldsymbol{\phi})(x,t)$, in which the function
$B(x,t)$, named BLUES function, is taken to be equal to the Green function of
the chosen related linear system, i.e., $B(x,t)\equiv G(x,t)$ and the new
(combined) source $\boldsymbol{\phi}(x,t)$ is to be calculated by systematic
iteration, using the given (combined) source $\boldsymbol{\varphi}(x,t)$. This
procedure starts from the following implicit equation, which makes use of the
action of the residual operator,

$\begin{split}\mathcal{R}_{x}\,(B\ast\boldsymbol{\phi})&=\mathcal{L}_{t,x}\,(B\ast\boldsymbol{\phi})-\mathcal{N}_{t,x}\,(B\ast\boldsymbol{\phi})\\\
&=\boldsymbol{\phi}-\boldsymbol{\varphi}\,.\end{split}$ (9)

To find the solution to the nonlinear system (2), equation (9) can be iterated
to calculate an approximation for $\boldsymbol{\phi}$ in the form of a
sequence in powers of the residual $\mathcal{R}_{x}\,$. By subsequently taking
the convolution product with $B(x,t)$, approximate solutions
$\boldsymbol{U}_{\phi}^{(n)}(x,t)$ to (2) can be found, i.e.,

$\boldsymbol{U}^{(n)}_{\phi}(x,t)=(B\ast\boldsymbol{\phi}^{(n)})(x,t)=\boldsymbol{U}^{(0)}_{\phi}(x,t)+\left(B\ast\mathcal{R}_{x}\,\boldsymbol{U}^{(n-1)}_{\phi}\right)(x,t)\,,$
(10)

where

$\boldsymbol{U}^{(0)}_{\phi}(x,t)=(B\ast\boldsymbol{\phi}^{(0)})(x,t)=(B\ast\boldsymbol{\varphi})(x,t).$
(11)

is the zeroth-order convolution product in which
$\boldsymbol{\phi}^{(0)}\equiv\boldsymbol{\varphi}$. We now apply this to a
second-order-in-time PDE.

##  3 Nonlinear telegrapher equation

From now on we will work with dimensionless functions and variables (which are
obtained after suitable scaling). The linear telegrapher equation with
$\alpha,c\in\mathbb{R}$

$u_{tt}(x,t)+\alpha u_{t}(x,t)-c^{2}u_{xx}(x,t)=0\,,$ (12)

with $u_{tt}\equiv\partial^{2}u/\partial t^{2}$, $u_{t}\equiv\partial
u/\partial t$, etc., is used in a broad spectrum of scientific disciplines,
ranging from electrical transmission in cables, where $u$ is then the voltage
per unit length of cable, to the statistical mechanics of active matter [13,
14] and mathematical biology [15], where $u$ represents a particle or
population density, respectively. Equation (12) interpolates between the wave
equation ($\alpha\rightarrow 0$ with $c$ fixed) and the diffusion equation
($\alpha\rightarrow\infty$ and $c\rightarrow\infty$, with
$c^{2}/\alpha\rightarrow D$ constant) and possesses a finite propagation speed
for disturbances, in contrast to the infinite propagation speed for the
diffusion equation [16], which is a result of the solution being strictly
positive on $\mathbb{R}$ for all $t>0$. If at $t=0$ one assumes e.g., that
$u(x,0)=\delta(x)$ then immediately after $t=0$, the solution jumps to a
nonzero positive value over the entire real line, which is not compact. The
information contained in the initial condition is instantaneously transported
over the entire real line. Consequently, the propagation speed is infinite for
the diffusion equation.

For early times $t\ll T$ where $T=1/\alpha$, the solution of equation (12) is
wave-like, while for later times $t\gg T$, the solution is diffusive.

We now consider the telegrapher equation with a quadratic nonlinearity and
constant forcing, i.e.,

$\begin{split}\widetilde{\mathcal{N}}^{[n]}_{t,x}\,u&=u_{tt}+\alpha
u_{t}-c^{2}u_{xx}+u^{2}=1\,,\\\ u(x,0)&=f(x)\,,\\\
\left.u_{t}(x,t)\right|_{t=0}&=g(x)\,,\end{split}$ (13)

where $f(x)$ and $g(x)$ are the initial conditions and $\chi(x,t)=1$ is the
external source. We will now attempt to generate approximate solutions to (13)
using the formalism developed in section 2.

We can convert this second-order-in-time PDE (13) into two coupled first-
order-in-time PDEs by introducing $v=u_{t}$, i.e.,

$\begin{split}u_{t}&=v\\\ v_{t}&=c^{2}u_{xx}-\alpha v-u^{2}+1\,,\end{split}$
(14)

where $\boldsymbol{\psi}(x,t)=(0,1)^{\intercal}$. We can rewrite the system
(14) with solution vector $\boldsymbol{U}=(u,v)^{\intercal}$ in the form
required by the matrix BLUES function method (2),

$\mathcal{N}_{t,x}\,\boldsymbol{U}(x,t)=\boldsymbol{\varphi}(x,t)\,.$ (15)

The linear operator $\mathcal{L}_{t}\,$ (which is now independent of $x$),
residual $\mathcal{R}_{x}\,$ and source vector $\boldsymbol{\varphi}$ are
given through the following expressions

$\begin{split}\mathcal{L}_{t}\,\boldsymbol{U}&=\begin{pmatrix}u_{t}\\\
v_{t}\end{pmatrix}-A\begin{pmatrix}u\\\ v\end{pmatrix},\\\
\mathcal{R}_{x}\,\boldsymbol{U}&=\begin{pmatrix}0\\\
c^{2}u_{xx}-u^{2}\end{pmatrix},\\\
\boldsymbol{\varphi}&=\delta(t)\begin{pmatrix}f(x)\\\
g(x)\end{pmatrix}+\begin{pmatrix}0\\\ 1\end{pmatrix}\,.\\\ \end{split}$ (16)

The matrix of coefficients $A$ of the linear operator is

$A=\begin{pmatrix}0&1\\\ 0&-\alpha\end{pmatrix}\,,$ (17)

which is independent of $x$, and hence the matrix Green function for the
linear operator is the matrix exponential $G(x,t)=\delta(x)\exp(At)\Theta(t)$,
i.e.,

$G(t)=\delta(x)\begin{pmatrix}1&\frac{1-\mathrm{e}^{-\alpha t}}{\alpha}\\\
0&\mathrm{e}^{-\alpha t}\end{pmatrix}\Theta(t)\,.$ (18)

We can now set up the BLUES iteration procedure. The $n$th approximant for the
solution of (14) is calculated by

$\boldsymbol{U}^{(n)}(x,t)=\boldsymbol{U}^{(0)}(x,t)+\int_{0^{-}}^{t}\int_{\mathbb{R}}G(x-x^{\prime},t-t^{\prime})\mathcal{R}_{x^{\prime}}\,\boldsymbol{U}^{(n-1)}(x^{\prime},t^{\prime})\mathrm{d}x^{\prime}\mathrm{d}t^{\prime}\,,$
(19)

where $\boldsymbol{U}^{(0)}(x,t)$ is the zeroth approximant and is defined as
the convolution product of the Green function and the source
$\boldsymbol{\varphi}(x,t)$, i.e,

$\boldsymbol{U}^{(0)}(x,t)=(G\ast\boldsymbol{\varphi})(x,t)=\left(f(x)+\frac{t}{\alpha}+\frac{1-\alpha
g(x)}{\alpha^{2}}(\mathrm{e}^{-\alpha t}-1),\frac{1}{\alpha}-\frac{1-\alpha
g(x)}{\alpha}\mathrm{e}^{-\alpha t}\right)^{\intercal}\,.$ (20)

For example, for initial conditions $f(x)=1+\sin x$ and $g(x)=0$, the zeroth
approximant is (for $\alpha=1$)

$\boldsymbol{U}^{(0)}(x,t)=\left(\sin
x+t+\mathrm{e}^{-t},1-\mathrm{e}^{-t}\right)^{\intercal}$ (21)

and the first approximant
$\boldsymbol{U}^{(1)}(x,t)=(u^{(1)}(x,t),v^{(1)}(x,t))$ consists of

$\begin{split}u^{(1)}(x,t)&=\sin{x}\left[\mathrm{e}^{-t}(2t-(c^{2}-4))-t^{2}-(c^{2}-2)t+(c^{2}-3)\right]+\frac{1}{2}\cos{2x}\left[\mathrm{e}^{-t}+t-1\right]\\\
&-\frac{1}{6}\left[3\mathrm{e}^{-2t}-\mathrm{e}^{-t}(6t^{2}+12t+9)+(t^{3}-6t^{2}+9t)\right]\,,\\\
v^{(1)}(x,t)&=u^{(1)}_{t}(x,t)\,.\end{split}$ (22)

Note that when a higher-order-in-time PDE is converted into a system of
coupled PDEs, the residual is zero in the individual channels of the “new”
variables, i.e., the derivatives. The matrix formalism therefore decouples
into individual integrations of the elements in the last column of the Green
function matrix with the residual applied to the previous approximant for the
solution $u(x,t)$. Hence, it is only necessary to perform the integration
corresponding to the channel for the solution $u(x,t)$ and not for all the
time derivatives. However, this is not the case when the system is _a priori_
coupled.

In Figs. 1 and 2, we compare the approximate solutions of the matrix BLUES
function method with the VIM, ADM, HPM and GVIM for $\alpha=1$ and $c=1$, with
initial conditions $f(x)=1+\sin x$ and $g(x)=0$. All numerical results are
obtained with Mathematica’s NDSolve function, which uses the numerical method
of lines to discretize the PDE in all but one dimension and then integrates
this semi-discrete problem as a system of ODEs. In Fig. 1, the spatial profile
of the approximants is shown for a fixed time $t=1$. The matrix BLUES function
method and the GVIM, which are almost concident, approximate the numerical
solution quite well. This is confirmed in Fig. 2, where the time evolution of
the approximants for fixed position $x=-\pi/2$ is shown. Note, however, that
all of the implemented methods diverge for $t\rightarrow\infty$. The accuracy
can be increased by considering higher-order iterations. Note that while all
of the methods can quite accurately reproduce the spatial profile at the
global minima of the numerical solution, the performance differs significantly
in the region where the maxima and local minima occur (e.g., around
$x=\pi/2$).

Figure 1: A comparison is made between the numerical solution of equation (13)
(red, full line), the matrix BLUES function method (black, full line), the VIM
(blue, dot-dashed line), the ADM (green, dotted line), the HPM (brown, dot-
dash-dashed line) and the GVIM (orange, dot-dot-dashed line). For all methods
the third iteration or third order result is shown. The inset shows details
around the local minimum at $x=\pi/2$. The time is fixed at $t=1$ and
parameter values are $\alpha=c=1$. Figure 2: A comparison is made between the
numerical solution of equation (13) (red, full line), the matrix BLUES
function method (black, full line), the VIM (blue, dot-dashed line), the ADM
(green, dotted line), the HPM (brown, dot-dash-dashed line) and the GVIM
(orange, dot-dot-dashed line). For all methods the third iteration or third
order result is shown. The position is fixed at $x=-\pi/2$ and parameter
values are $\alpha=c=1$. Figure 3: Spatial plot of the solution of the
nonlinear telegrapher equation (13) for a fixed time $t=1$. A comparison is
made between the numerical solution (red, full line) and the
$n\in\\{0,1,2,3\\}$ matrix BLUES method approximants (black lines). The
parameter values are $\alpha=c=1$.

In Fig. 3, the consecutive matrix BLUES method approximants for
$n\in\\{0,1,2,3\\}$ are compared with the numerical solution at $t=1$. Note
that the local minimum around $x=\pi/2$ as well as the entire solution is well
reproduced by the third approximant.

##  4 Conclusions

In this paper we have extended the BLUES function method that was developed
for nonlinear PDEs with a first-order time derivative to second-order (and
higher) time derivatives by introducing a matrix formalism. We have applied
this extension of the BLUES function method to a second-order-in-time
nonlinear telegrapher equation with a quadratic nonlinearity and a
trigonometric initial condition. From a comparison with the ADM, the VIM, the
HPM and the GVIM, we have shown that the matrix BLUES function method is able
to generate approximants that are a valuable alternative to what can be
obtained through other methods at the same level of iteration or the same
order of approximation. In future research on the BLUES function method, we
envision a possible extension of the method to stochastic (coupled) DEs for
which the noise can play the role of an external source.

## References

  * [1] G. Adomian, Solving Frontier Problems of Physics: The Decomposition Method, 1st Edition, Vol. 60 of Fundamental Theories of Physics, Springer Netherlands, 1994\.  doi:10.1007/978-94-015-8289-6.
  * [2] G. Adomian, A review of the decomposition method and some recent results for nonlinear equations, Math. Comput. Model. 13 (7) (1990) 17 – 43.  doi:https://doi.org/10.1016/0895-7177(90)90125-7.
URL http://www.sciencedirect.com/science/article/pii/0895717790901257

  * [3] J.-H. He, Variational iteration method — some recent results and new interpretations, J. Comput. Appl. Math. 207 (1) (2007) 3 – 17.  doi:https://doi.org/10.1016/j.cam.2006.07.009.
URL http://www.sciencedirect.com/science/article/pii/S0377042706004559

  * [4] S. Khuri, A. Sayfy, Variational iteration method: Green’s functions and fixed point iterations perspective, Appl. Math. Lett. 32 (2014) 28 – 34.  doi:https://doi.org/10.1016/j.aml.2014.01.006.
URL http://www.sciencedirect.com/science/article/pii/S0893965914000433

  * [5] J.-H. He, Homotopy perturbation technique, Comput. Methods Appl. Mech. Eng. 178 (3) (1999) 257–262.  doi:https://doi.org/10.1016/S0045-7825(99)00018-3.
URL https://www.sciencedirect.com/science/article/pii/S0045782599000183

  * [6] S. Liao, Homotopy Analysis Method in Nonlinear Differential Equations, 1st Edition, Springer-Verlag Berlin Heidelberg, 2012.  doi:10.1007/978-3-642-25132-0.
  * [7] J. O. Indekeu, K. K. Müller-Nedebock, BLUES function method in computational physics, J. Phys. A-Math. Theor. 51 (16) (2018) 165201.  doi:10.1088/1751-8121/aab345.
URL https://doi.org/10.1088/1751-8121/aab345

  * [8] J. Berx, J. O. Indekeu, Analytic iteration procedure for solitons and traveling wavefronts with sources, J. Phys. A-Math. Theor. 52 (38) (2019) 38LT01.  doi:10.1088/1751-8121/ab3914.
URL https://doi.org/10.1088%2F1751-8121%2Fab3914

  * [9] J. Berx, J. O. Indekeu, BLUES iteration applied to nonlinear ordinary differential equations for wave propagation and heat transfer, J. Phys. A-Math. Theor. 54 (2) (2020) 025702.  doi:10.1088/1751-8121/abcf57.
URL https://doi.org/10.1088/1751-8121/abcf57

  * [10] J. Berx, J. O. Indekeu, BLUES function method applied to partial differential equations and analytic approximants for interface growth under shear, Phys. Rev. Research 3 (2021) 033113\.  doi:10.1103/PhysRevResearch.3.033113.
URL https://link.aps.org/doi/10.1103/PhysRevResearch.3.033113

  * [11] J. Berx, J. O. Indekeu, Epidemic processes with vaccination and immunity loss studied with the BLUES function method, Physica A 590 (2022) 126724.  doi:https://doi.org/10.1016/j.physa.2021.126724.
URL https://www.sciencedirect.com/science/article/pii/S0378437121009341

  * [12] J. O. Indekeu, R. Smets, Traveling wavefront solutions to nonlinear reaction-diffusion-convection equations, J. Phys. A-Math. Theor. 50 (31) (2017) 315601.  doi:10.1088/1751-8121/aa7a93.
URL https://doi.org/10.1088%2F1751-8121%2Faa7a93

  * [13] T. Demaerel, C. Maes, Active processes in one dimension, Phys. Rev. E 97 (2018) 032604.  doi:10.1103/PhysRevE.97.032604.
URL https://link.aps.org/doi/10.1103/PhysRevE.97.032604

  * [14] G. H. Weiss, Some applications of persistent random walks and the telegrapher’s equation, Physica A: Statistical Mechanics and its Applications 311 (3) (2002) 381–410.  doi:https://doi.org/10.1016/S0378-4371(02)00805-1.
URL https://www.sciencedirect.com/science/article/pii/S0378437102008051

  * [15] W. Alharbi, S. Petrovskii, Critical domain problem for the reaction – telegraph equation model of population dynamics, Mathematics 6 (4) (2018).  doi:10.3390/math6040059.
URL https://www.mdpi.com/2227-7390/6/4/59

  * [16] L. C. Evans, Partial differential equations, 2nd Edition, American Mathematical Society, 2010.