# Comment on “First-Principles Calculation of the Optical Rotatory Power of
Periodic Systems: Modern Theory with Modern Functionals”

Jacques K. Desmarais<EMAIL_ADDRESS>Dipartimento di
Chimica, Università di Torino, via Giuria 5, 10125 Torino, Italy  Université
de Pau et des Pays de l’Adour, CNRS, IPREM, E2S UPPA, Pau, France

###### Abstract

The article [Desmarais et al. Phys Rev. B. 107 224430 (2023)] provided a
generalization of the “modern theory of orbital magnetization” to include non-
local Hamiltonians (e.g. hybrid functionals of the generalized Kohn-Sham
theory) for magnetic response properties. The formulation was employed for the
calculation of the optical rotatory power of periodic systems, where results
indicated inequivalence between sampling of direct and reciprocal spaces for
those calculations far from the complete basis set limit. We show that such
results can be explained by the fact that the reported calculations correspond
to an approximate treatment of the action of the $\boldsymbol{\nabla_{k}}$
operator on Bloch orbitals. In the article, calculations are missing a hidden
“response” contribution to the reciprocal-space derivatives. The missing
response term is shown to (generally) affect the results of calculations of
not only magnetic, but also electric response properties, within the context
of the “modern theory of polarization”. Necessary conditions are provided to
permit avoiding the calculation of the hidden response term.

Desmarais et al. Desmarais et al. (2023) (hereafter referred to as “the
article”) provided a generalization of Ceresoli et al.’s Ceresoli et al.
(2006) “modern theory of orbital magnetization” to include non-local
Hamiltonians and applied the theory to the calculation of the optical rotatory
power (OR) of periodic systems.

The article reported calculations on infinite chains of H2O2, in which visible
differences in the calculated OR are reported for:

  1. i)

the infinite periodic system versus the large finite system

  2. ii)

the 3$\times$ replicated supercell with 3 evenly spaced $\mathbf{k}$ points
versus the 9$\times$ replicated supercell with 1 $\mathbf{k}$ point

In both cases i) and ii), the differences diminish as the calculation
approaches the complete basis set limit. In the article, these differences
(suggestive of a “non-periodic” formulation) were attributed to the gauge-
origin dependence of the first-order magnetic Hamiltonian:Desmarais et al.
(2023)

$H^{(1)}_{\rm
mag}=\frac{1}{2}\left(\mathbf{r}+i\boldsymbol{\nabla_{k}}\right)\wedge\mathbf{p}+{\rm
H.c.}\;.$ (1)

Here we correct this assertion and show that the differences are instead the
result of an approximate treatment therein of the action of the
$\boldsymbol{\nabla_{k}}$ operator. In fact, we show that the same differences
i) and ii) between the infinite periodic vs. large finite system, as well as
between uniform sampling of direct and reciprocal spaces (i.e. uniform
sampling of $\mathbf{k}$-points vs. supercells) are not only obtained for
magnetic properties from Eq. (1), but also electric properties, within a
similar approximation to the action of $\boldsymbol{\nabla_{k}}$. That is, the
same “non-periodic” behaviour is found for calculation of those properties
employing the first-order electric Hamiltonian:Blount (1962); Rérat et al.
(2008); Kirtman et al. (2000); Bishop et al. (2001); Ferrero et al. (2008a, b)

$H^{(1)}_{\rm
ele}=\frac{1}{2}\left(\mathbf{r}+i\boldsymbol{\nabla_{k}}\right)+{\rm
H.c.}\;.$ (2)

which coincides exactly with King-Smith, Vanderbilt and Resta’s “modern theory
of polarization”.King-Smith and Vanderbilt (1993); Resta (1994, 1998);
Springborg and Kirtman (2008) We begin by a review of the state of the art in
application of the $\boldsymbol{\nabla_{k}}$ operator.

The action of $\boldsymbol{\nabla_{k}}$ on Bloch orbitals built from atom-
centered atomic orbitals (AOs) $|\mu^{\mathbf{g}}\rangle$:

$|\psi_{i}\left(\mathbf{k}\right)\rangle=\sum_{\mu}C_{\mu,i}\left(\mathbf{k}\right)|\phi_{\mu}\left(\mathbf{k}\right)\rangle=\sum_{\mu}C_{\mu,i}\left(\mathbf{k}\right)\sum_{\mathbf{g}}e^{i\mathbf{k}\cdot\mathbf{g}}|\mu^{\mathbf{g}}\rangle$
(3)

is trivial to apply on the $|\phi_{\mu}\left(\mathbf{k}\right)\rangle$ part.
The problem of calculating the derivative of the orbital coefficients
$\boldsymbol{\nabla_{k}}C_{\mu,i}\left(\mathbf{k}\right)$ is more subtle.

In general, an expansion:

$\boldsymbol{\nabla_{k}}C_{\mu,i}\left(\mathbf{k}\right)=\sum_{l}^{\rm
all}C_{\mu,l}\left(\mathbf{k}\right)Q_{l,i}\left(\mathbf{k}\right)$ (4)

with (as of yet undetermined) coefficients $Q_{i,l}\left(\mathbf{k}\right)$
provides a solution. To find the coefficients
$Q_{i,l}\left(\mathbf{k}\right)$, the derivative is typically applied to the
Kohn-Sham (KS) single-particle equation
$\hat{F}|\psi_{i}\left(\mathbf{k}\right)\rangle=\epsilon_{i,\mathbf{k}}|\psi_{i}\left(\mathbf{k}\right)\rangle$
yielding:Kirtman et al. (2000); Bishop et al. (2001)

$Q_{i,l}\left(\mathbf{k}\right)=\frac{K_{i,l}\left(\mathbf{k}\right)-\epsilon_{l,\mathbf{k}}R_{i,l}\left(\mathbf{k}\right)}{\epsilon_{l,\mathbf{k}}-\epsilon_{i,\mathbf{k}}}\quad
l\neq i$ (5a) in which $\mathbf{K}$ and $\mathbf{R}$ are the derivatives of
the KS Hamiltonian $\mathbf{F}$ and basis-function overlap $\mathbf{S}$
matrices at fixed orbital coefficients: $\displaystyle
K_{i,l}\left(\mathbf{k}\right)\to
i\sum_{\mathbf{g}}\mathbf{g}e^{i\mathbf{k}\cdot\mathbf{g}}F_{i,l}\left(\mathbf{g}\right)$
(5b) $\displaystyle
R_{i,l}\left(\mathbf{k}\right)=i\sum_{\mathbf{g}}\mathbf{g}e^{i\mathbf{k}\cdot\mathbf{g}}S_{i,l}\left(\mathbf{g}\right)$
(5c)

Here we show that Eq. (5b) is, in fact, an approximation to the full
$K_{i,l}\left(\mathbf{k}\right)$ (i.e. the $\to$ should be replaced by an
approximate equal sign $\approx$), in that “orbital-relaxation” or “response”
contributions to the derivatives are dropped. This approximation explains the
discrepencies i) and ii) in the approximate calculation of responses to
external electromagnetic fields.

To develop the exact treatment, let us now approach the problem of the
analytical calculation of $\mathbf{k}$-space derivatives from the perspective
of linear-response theory. We consider a small displacement $\mathbf{h}$ away
from $\mathbf{k}$ and write the displaced KS single-particle equation, as well
as orthonormality condition:

$\displaystyle\mathbf{F}\left(\mathbf{k}+\mathbf{h}\right)\mathbf{C}\left(\mathbf{k}+\mathbf{h}\right)=\mathbf{S}\left(\mathbf{k}+\mathbf{h}\right)\mathbf{C}\left(\mathbf{k}+\mathbf{h}\right)\boldsymbol{\varepsilon}\left(\mathbf{k}+\mathbf{h}\right)$
(6a)
$\displaystyle\mathbf{C}^{\dagger}\left(\mathbf{k}+\mathbf{h}\right)\mathbf{S}\left(\mathbf{k}+\mathbf{h}\right)\mathbf{C}\left(\mathbf{k}+\mathbf{h}\right)=\boldsymbol{1}$
(6b)
$\displaystyle\mathbf{C}\left(\mathbf{k}+\mathbf{h}\right)\equiv\mathbf{C}^{(0)}\left(\mathbf{k}\right)\mathbf{Q}\left(\mathbf{k}+\mathbf{h}\right)$
(6c)

with (generally) non-canonical displaced Lagrange multipliers
$\boldsymbol{\varepsilon}\left(\mathbf{k}+\mathbf{h}\right)$ and (as of yet)
undetermined coefficients $\mathbf{Q}\left(\mathbf{k}+\mathbf{h}\right)$.
Expanding all quantities in a power series around the point $\mathbf{k}$:

$\displaystyle\mathbf{F}\left(\mathbf{k}+\mathbf{h}\right)$ $\displaystyle=$
$\displaystyle\mathbf{F}^{(0)}\left(\mathbf{k}\right)+\mathbf{h}\
\mathbf{F}^{(1)}\left(\mathbf{k}\right)+\dots$ (7a)
$\displaystyle\mathbf{C}\left(\mathbf{k}+\mathbf{h}\right)$ $\displaystyle=$
$\displaystyle\mathbf{C}^{(0)}\left(\mathbf{k}\right)+\mathbf{h}\
\mathbf{C}^{(1)}\left(\mathbf{k}\right)+\dots$ (7b)
$\displaystyle\mathbf{S}\left(\mathbf{k}+\mathbf{h}\right)$ $\displaystyle=$
$\displaystyle\mathbf{S}^{(0)}\left(\mathbf{k}\right)+\mathbf{h}\
\mathbf{S}^{(1)}\left(\mathbf{k}\right)+\dots$ (7c)
$\displaystyle\boldsymbol{\varepsilon}\left(\mathbf{k}+\mathbf{h}\right)$
$\displaystyle=$
$\displaystyle\boldsymbol{\epsilon}^{(0)}\left(\mathbf{k}\right)+\mathbf{h}\
\boldsymbol{\varepsilon}^{(1)}\left(\mathbf{k}\right)+\dots$ (7d)
$\displaystyle\mathbf{Q}\left(\mathbf{k}+\mathbf{h}\right)$ $\displaystyle=$
$\displaystyle\mathbf{1}+\mathbf{h}\
\mathbf{Q}^{(1)}\left(\mathbf{k}\right)+\dots$ (7e)

and taking the derivative of both sides of Eq. (7c) gives:

$\mathbf{S}^{(1)}\left(\mathbf{k}\right)=\left.\frac{\partial\mathbf{S}\left(\mathbf{k}+\mathbf{h}\right)}{\partial\mathbf{h}}\right|_{\mathbf{h=0}}$
(8)

then, inserting Eqs.(7b) and (7e) into Eq. (6c) yields:

$\mathbf{C}^{(1)}\left(\mathbf{k}\right)=\mathbf{C}^{(0)}\left(\mathbf{k}\right)\mathbf{Q}^{(1)}\left(\mathbf{k}\right)$
(9)

At this point, it is convenient to define the matrices:

$\mathbf{K}^{(1)}\left(\mathbf{k}\right)=\left[\mathbf{C}^{(0)}\left(\mathbf{k}\right)\right]^{\dagger}\mathbf{F}^{(1)}\left(\mathbf{k}\right)\mathbf{C}^{(0)}\left(\mathbf{k}\right)$
(10)

and:

$\displaystyle\mathbf{R}^{(1)}\left(\mathbf{k}\right)=\left[\mathbf{C}^{(0)}\left(\mathbf{k}\right)\right]^{\dagger}\mathbf{S}^{(1)}\left(\mathbf{k}\right)\mathbf{C}^{(0)}\left(\mathbf{k}\right)$
(11)

Inserting Eqs. (7) and (9) in Eq. (6) and collecting terms in the first order,
then left-multiplying by
$\left[\mathbf{C}^{(0)}\left(\mathbf{k}\right)\right]^{\dagger}$ and using
Eqs. (10) and (11) leads directly to the first order perturbation equation:

$\displaystyle\mathbf{K}^{(1)}\left(\mathbf{k}\right)+\boldsymbol{\epsilon}^{(0)}\left(\mathbf{k}\right)\mathbf{Q}^{(1)}\left(\mathbf{k}\right)=\mathbf{R}^{(1)}\left(\mathbf{k}\right)\boldsymbol{\epsilon}^{(0)}\left(\mathbf{k}\right)$
(12) $\displaystyle+$
$\displaystyle\mathbf{Q}^{(1)}\left(\mathbf{k}\right)\boldsymbol{\epsilon}^{(0)}\left(\mathbf{k}\right)+\boldsymbol{\varepsilon}^{(1)}\left(\mathbf{k}\right)$

Eq. (12) must be solved self-consistently for
$\mathbf{Q}^{(1)}\left(\mathbf{k}\right)$, under the condition of
orthonormality. The requisite first order orthonormality condition may now be
written by inserting Eqs. (7b), (7c), (9) and (11) into Eq. (6b) to get:

$\left[\mathbf{Q}^{(1)}\left(\mathbf{k}\right)\right]^{\dagger}+\mathbf{Q}^{(1)}\left(\mathbf{k}\right)=-\mathbf{R}^{(1)}\left(\mathbf{k}\right)$
(13)

The standard non-canonical solution procedure Karna and Dupuis (1991); Pople
et al. (1979) provides, by taking advantage of the fact that the occ-virt
blocks of perturbed Lagrange multipliers are vanishing:

$\varepsilon^{(1)}_{ia}\left(\mathbf{k}\right)=\varepsilon^{(1)}_{ai}\left(\mathbf{k}\right)=0\quad
i\in occ,a\in virt$ (14)

the following solution for the occ-virt and virt-occ blocks of
$\mathbf{Q}^{(1)}$:

$Q^{(1)}_{ia}\left(\mathbf{k}\right)=\frac{K^{(1)}_{ia}\left(\mathbf{k}\right)-\epsilon^{(0)}_{a}\left(\mathbf{k}\right)R^{(1)}_{ia}\left(\mathbf{k}\right)}{\epsilon^{(0)}_{a}\left(\mathbf{k}\right)-\epsilon^{(0)}_{i}\left(\mathbf{k}\right)}\quad
i\in occ,a\in virt$ (15a) and for the occ-occ block (by imposing Hermiticity):
$\displaystyle Q^{(1)}_{ij}\left(\mathbf{k}\right)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}R^{(1)}_{ij}\left(\mathbf{k}\right)\quad i,j\in occ$
(15b)

with (here and elsewhere) exactly analogous expressions for the virt-virt
block. Eq. (15b) is consistent with the following non-canonical matrices of
Lagrange multipliers:

$\varepsilon^{(1)}_{ij}\left(\mathbf{k}\right)=K^{(1)}_{ij}\left(\mathbf{k}\right)-\frac{1}{2}\left(\epsilon^{(0)}_{i}+\epsilon^{(0)}_{j}\right)R^{(1)}_{ij}\left(\mathbf{k}\right)\quad
i,j\in occ$ (16)

In Eqs. (15a) and (16) $K^{(1)}_{ll^{\prime}}\left(\mathbf{k}\right)$ is the
first-order perturbed KS Hamiltonian matrix:

$\displaystyle iK^{(1)}_{l,l^{\prime}}\left(\mathbf{k}\right)$
$\displaystyle=$
$\displaystyle-\sum_{\mathbf{g}}\mathbf{g}e^{i\mathbf{k}\cdot\mathbf{g}}F_{l,l^{\prime}}^{(0)}\left(\mathbf{g}\right)+\sum_{\mathbf{g}}e^{i\mathbf{k}\cdot\mathbf{g}}$
(17) $\displaystyle\times$
$\displaystyle\sum_{\mu\nu}C_{\mu,l}^{(0)\ast}\left(\mathbf{k}\right)C_{\nu,l^{\prime}}^{(0)}\left(\mathbf{k}\right)V^{(1)}_{\mu\nu}\left(\mathbf{g}\right)$

with a first term representing the contribution as in Eq. (5b) from standard
approaches, and the second term is an additional “orbital-relaxation” or
“response” correction. Thus, the response term is proportional to the KS
potential $V^{(1)}$ depending on the derivative of the reduced density matrix
coefficients $P^{(1)}_{\mu\nu}$:

$\displaystyle P^{(1)}_{\mu\nu}\left(\mathbf{g}\right)$ $\displaystyle\approx$
$\displaystyle\frac{\partial}{\partial\mathbf{h}}\left.\frac{2}{\Omega}\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{\Omega}\Re\
e^{i\left[\mathbf{k}+\mathbf{h}\right]\cdot\mathbf{g}}P_{\mu\nu}\left(\mathbf{k}+\mathbf{h}\right)\right|_{\mathbf{h=0}}$
$\displaystyle=$ $\displaystyle\frac{2}{\Omega}\Re\
e^{i\mathbf{k}\cdot\mathbf{g}}\sum_{i}^{\rm occ}\sum_{l}^{\rm
all}C_{\mu,l}^{(0)}\left(\mathbf{k}\right)iQ_{l,i}^{(1)}\left(\mathbf{k}\right)\left[C_{\nu,i}^{(0)}\left(\mathbf{k}\right)\right]^{\ast}$
$\displaystyle+$ $\displaystyle
C_{\mu,i}^{(0)}\left(\mathbf{k}\right)i\left[Q_{i,l}^{(1)}\left(\mathbf{k}\right)\right]^{\ast}\left[C_{\nu,l}^{(0)}\left(\mathbf{k}\right)\right]^{\ast}$
(27)

with summation being over quadrature points in the volume $\Omega$ of the
first Brillouin zone (FBZ).

Once $\mathbf{Q}^{(1)}$ has been obtained from a non-canonical solution of Eq.
(12) we need to transform the Bloch orbitals at point $\mathbf{k}+\mathbf{h}$
to canonical ones in order to use them as field-free orbitals for perturbation
by an electric or magnetic field. We can transform the Bloch orbitals to
canonical form by finding the unitary matrix $\mathbf{T}$ that diagonalizes
the occ-occ (or virt-virt) block of the matrix of Lagrange multipliers
$\boldsymbol{\varepsilon}_{OO}$, that is
$\mathbf{T}^{\dagger}_{O}\boldsymbol{\varepsilon}_{OO}\mathbf{T}_{O}=\boldsymbol{\epsilon}_{OO}$.
This means that we need to solve the following eigenvalue equation:

$\boldsymbol{\varepsilon}_{OO}\left(\mathbf{k}+\mathbf{h}\right)\mathbf{T}_{O}\left(\mathbf{k}+\mathbf{h}\right)=\mathbf{T}_{O}\left(\mathbf{k}+\mathbf{h}\right)\boldsymbol{\epsilon}_{O}\left(\mathbf{k}+\mathbf{h}\right)$
(28)

To obtain the orbital energy $\epsilon_{i}$ and corresponding eigenvector
$\mathbf{T}_{i}$ at point $\mathbf{k}+\mathbf{h}$ in reciprocal space.

Once the matrix $\mathbf{T}_{O}$ has been obtain, returning to (28), we find:

$\displaystyle\mathbf{F}\left(\mathbf{k}+\mathbf{h}\right)\mathbf{C}^{\prime}\left(\mathbf{k}+\mathbf{h}\right)$
$\displaystyle=$ $\displaystyle\mathbf{S}\left(\mathbf{k}+\mathbf{h}\right)$
(29) $\displaystyle\times$
$\displaystyle\mathbf{C}^{\prime}\left(\mathbf{k}+\mathbf{h}\right)\boldsymbol{\epsilon}\left(\mathbf{k}+\mathbf{h}\right)$

where:

$\mathbf{C}^{\prime}\left(\mathbf{k}+\mathbf{h}\right)=\mathbf{C}\left(\mathbf{k}+\mathbf{h}\right)\mathbf{T}\left(\mathbf{k}+\mathbf{h}\right)$
(30)

Then, defining:

$\mathbf{T}\left(\mathbf{k}+\mathbf{h}\right)=\mathbf{1}+\mathbf{h}\
\mathbf{T}^{(1)}\left(\mathbf{k}\right)+\dots$ (31)

and proceeding as in Eqs. (7)-(9), we obtain:

$\mathbf{C}^{\prime(1)}\left(\mathbf{k}\right)=\mathbf{C}^{(0)}\left(\mathbf{k}\right)\mathbf{Q}^{\prime(1)}\left(\mathbf{k}\right)$
(32)

with:

$\mathbf{Q}^{\prime(1)}\left(\mathbf{k}\right)=\mathbf{Q}^{(1)}\left(\mathbf{k}\right)+\mathbf{T}^{(1)}\left(\mathbf{k}\right)$
(33)

We can calculate $\mathbf{T}^{(1)}\left(\mathbf{k}\right)$ by solving the
eigenvalue Eq. (28) by non-degenerate Rayleigh-Schrodinger perturbation theory
giving:

$\mathbf{T}^{(1)}_{i}\left(\mathbf{k}\right)=\sum_{j\neq
i}\mathbf{T}^{(0)}_{j}\left(\mathbf{k}\right)\frac{\left[\mathbf{T}^{(0)}_{j}\left(\mathbf{k}\right)\right]^{\dagger}\boldsymbol{\varepsilon}^{(1)}\left(\mathbf{k}\right)\mathbf{T}^{(0)}_{i}\left(\mathbf{k}\right)}{\epsilon^{(0)}_{i}-\epsilon^{(0)}_{j}}$
(34)

We note in passing that degenerate or quasi-degenerate states would require
appropriate modification of our treatment.

Then, inserting Eqs. (14) and (16) in Eq. (34) we obtain:

$T^{(1)}_{ia}\left(\mathbf{k}\right)=0\quad i\in occ,a\in virt$ (35a)
$T^{(1)}_{ji}\left(\mathbf{k}\right)=\frac{\varepsilon^{(1)}_{ji}\left(\mathbf{k}\right)}{\epsilon^{(0)}_{i}\left(\mathbf{k}\right)-\epsilon^{(0)}_{j}\left(\mathbf{k}\right)}\quad
i\in occ,j\in occ$ (35b)

Note that substitution of Eq. (35b) for occ-occ and virt-virt blocks, along
with Eqs. (14)-(16) into Eq. (33) gives $\mathbf{Q}^{\prime(1)}$ that exactly
coincide with the canonical approach of Eq. (5), apart from the correction due
to the response term:

$Q^{\prime(1)}_{ll^{\prime}}\left(\mathbf{k}\right)=\frac{K^{(1)}_{ll^{\prime}}\left(\mathbf{k}\right)-\epsilon^{(0)}_{l^{\prime}}\left(\mathbf{k}\right)R^{(1)}_{ll^{\prime}}\left(\mathbf{k}\right)}{\epsilon^{(0)}_{l^{\prime}}\left(\mathbf{k}\right)-\epsilon^{(0)}_{l}\left(\mathbf{k}\right)}$
(36)

Indeed, Eq. (36) is exactly identical to Eq. (5a), except that
$K^{(1)}_{ll^{\prime}}$ now includes also the second term in Eq. (17) (the
response correction).

At this point, we note three important situations in which the new response
term correction vanishes:

  1. 1.

For a sufficiently fine sampling of the FBZ
$\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{\Omega}\to\int_{\Omega}d\mathbf{k}$
becomes an integral and
$P^{(1)}_{\mu\nu}\left(\mathbf{g}\right)=\frac{\partial}{\partial\mathbf{k}}\left.\mathcal{F}\left[P_{\mu\nu}\left(\mathbf{k}+\mathbf{h}\right)\right]\right|_{\mathbf{h=0}}$
vanishes because the Fourier transform
$\mathcal{F}\left[P_{\mu\nu}\left(\mathbf{k}+\mathbf{h}\right)\right]$ is
independent of $\mathbf{k}$.

  2. 2.

For a “real $\mathbf{k}$ point” (i.e. a point $\mathbf{k}_{\rm real}$ in which
$\sin\left(\mathbf{k}_{\rm real}\cdot\mathbf{g}\right)$ vanishes, that is to
say at the $\Gamma$ point or at the edge of the FBZ), the Hermiticity of
$iQ_{i,l}^{(1)}\left(\mathbf{k}\right)+\frac{i}{2}R_{i,l}^{(1)}\left(\mathbf{k}\right)$
leads to $P^{(1)}_{\mu\nu}\left(\mathbf{k}_{\rm real}\right)=0$

  3. 3.

At the complete basis set limit, the response term vanishes by virtue of the
Hellmann-Feynman theorem

From the above discussion, we expect the response term to be significant for
calculations employing few $\mathbf{k}$ points away from $\Gamma$ and away
from the edge of the FBZ, while employing a small basis set. Below we show
numerical results explicitly demonstrating this behaviour for responses to
electromagnetic fields. In practice, we admit that an actual calculation of
the response term is not computationally practical, because Eq. (Comment on
“First-Principles Calculation of the Optical Rotatory Power of Periodic
Systems: Modern Theory with Modern Functionals”) requires an independent self-
consistent calculation of $\mathbf{P}^{(1)}$ for each $\mathbf{k}$ point. In
fact, the numerical results below suggest that a fine sampling of
$\mathbf{k}$-space (resulting in a vanishing response term, regardless of the
basis set) is probably much cheaper and computationally convenient.

To show the effect of the missing response term for reciprocal space
derivatives, we may choose any response property to external electric and/or
magnetic fields that require the $Q_{i,l}$ matrix. Two such properties are
dynamical hyperpolarizabilities (a third order property to an electric field
perturbation) and OR (mixed second order property to magnetic and electric
field perturbations). Calculations are performed on infinite chains of H2O2 at
a light wavelength of 300 nm, employing the local density approximation and a
developer’s version of the Crystal23 code.Erba et al. (2022) For the
hyperpolarizabilities, we use the unmodified algorithm of Ref. Maschio et al.,
2015 and report second-harmonic generation (SHG) coefficients with three
electric perturbations along the x y and z Cartesian axes. For OR, we report
the trace of the optical rotation tensor, employing our gauge-origin
independent velocity formulation.Desmarais et al. (2023); Rérat and Kirtman
(2021)

Table 1: Second-harmonic generation coefficients in [a.u.] along XYZ (first three rows) per formula unit in infinite chains of H2O2 for different supercell and $\mathbf{k}$-point combinations, as well as basis sets. The last three columns report a similar analysis for the trace of the OR tensor in [${}^{\circ}/mm.$]. cells/$\mathbf{k}$-points | 1/500 | 2/500 | 1/2 | 2/1 | 1/3 | 3/1
---|---|---|---|---|---|---
cc-pvDz | 36.7 | 37.0 | -25.7 | -25.7 | 52.3 | 49.1
cc-pvTz | 710.7 | 625.5 | -235.7 | -235.7 | -1208 | -1212
cc-pvQz | 103.2 | 101.3 | -33.64 | -33.64 | 96.5 | 101.2
cc-pvDz | 47.9 | 47.9 | 19.6 | 19.6 | 58.1 | 47.3
cc-pvTz | 47.8 | 47.9 | 27.9 | 27.9 | 53.4 | 60.8
cc-pvQz | 53.5 | 53.5 | 38.7 | 38.7 | 62.0 | 62.7

Results are presented in Table LABEL:tab:gk for different basis sets and
supercell/$\mathbf{k}$-point combinations. The first three rows of the table
report SHG coefficients and the last three rows are for OR. At the limit of
many $\mathbf{k}$ points (first two columns), the response term is vanishing,
and consistent results are obtained per H2O2 unit with either a single cell or
a 2$\times$ supercell. This is especially true for large basis sets (another
criterion for diminishing the otherwise missing response term). For instance,
with 500 $\mathbf{k}$-points and the quadruple zeta basis set, we obtain SHG
coefficients of 101.3 and 103.2 a.u. (single and doubled cells, respectively),
but a larger disagreement of 710.7 vs. 625.5 a.u. with the triple zeta basis
set. More details are provided at Ref. qd_, . Calculations with one or two
$\mathbf{k}$ points include only “real” points (i.e. points $\mathbf{k}_{\rm
real}$ for which $\sin\left(\mathbf{k}_{\rm real}\cdot\mathbf{g}\right)$ is
vanishing) and the otherwise missing response term is also vanishing. Thus, we
always obtained perfectly consistent results both for SHG and OR between the
2/1 and 1/2 calculations (i.e. 2/1 meaning double cell and 1 $\mathbf{k}$
point, 1/2 meaning single cell and 2 $\mathbf{k}$ points). The 1/3
calculation, on the other hand, includes points away from $\Gamma$, and away
from the edge of the Brillouin zone, and at such points the response term is
not vanishing. In this case, we obtain significant differences between the 1/3
and 3/1 calculations, particularly for small basis sets (e.g. 58.1 vs. 47.3
${}^{\circ}/mm.$ for OR and a double-zeta basis set). The differences are
diminished by employing larger basis sets (e.g. 62.0 vs. 62.7 ${}^{\circ}/mm.$
for OR and a quadruple-zeta basis set). Thus, we may avoid explicit
calculation of the costly (and complicated) response term to reciprocal space
derivatives for electric and magnetic properties in different ways: namely, i)
by employing many $\mathbf{k}$ points and a small cell or ii) by employing
only $\mathbf{k}$ points at $\Gamma$ and the edge of the Brillouin zone and a
large supercell and/or iii) by employing a large basis set.

Acknowledgements I am grateful to Profs. Michel Rérat, Bernard Kirtman and
Michael Springborg for valuable discussions.

## References

  * Desmarais et al. (2023) J. K. Desmarais, B. Kirtman, and M. Rérat, Phys. Rev. B 107, 224430 (2023).
  * Ceresoli et al. (2006) D. Ceresoli, T. Thonhauser, D. Vanderbilt, and R. Resta, Phys. Rev. B 74, 024408 (2006).
  * Blount (1962) E. Blount, in _Solid state physics_ (1962), vol. 13, pp. 305–373.
  * Rérat et al. (2008) M. Rérat, M. Ferrero, E. Amzallag, I. Baraille, and R. Dovesi, in _J. Phys. Conf. Ser._ (2008), vol. 117, p. 012023.
  * Kirtman et al. (2000) B. Kirtman, F. L. Gu, and D. M. Bishop, J. Chem. Phys. 113, 1294 (2000).
  * Bishop et al. (2001) D. M. Bishop, F. L. Gu, and B. Kirtman, J. Chem. Phys. 114, 7633 (2001).
  * Ferrero et al. (2008a) M. Ferrero, M. Rérat, R. Orlando, and R. Dovesi, J. Comp. Chem. 29, 1450 (2008a).
  * Ferrero et al. (2008b) M. Ferrero, M. Rérat, R. Orlando, and R. Dovesi, J. Chem. Phys. 128, 014110 (2008b).
  * King-Smith and Vanderbilt (1993) R. King-Smith and D. Vanderbilt, Phys. Rev. B 47, 1651 (1993).
  * Resta (1994) R. Resta, Rev. Mod. Phys. 66, 899 (1994).
  * Resta (1998) R. Resta, Phys. Rev. Lett. 80, 1800 (1998).
  * Springborg and Kirtman (2008) M. Springborg and B. Kirtman, Phys. Rev. B 77, 045102 (2008).
  * Karna and Dupuis (1991) S. Karna and M. Dupuis, J. Comput. Chem. 12, 487 (1991).
  * Pople et al. (1979) J. Pople, R. Krishnan, H. Schlegel, and J. S. Binkley, Int. J. Quant. Chem. 16, 225 (1979).
  * Erba et al. (2022) A. Erba, J. K. Desmarais, S. Casassa, B. Civalleri, L. Donà, I. J. Bush, B. Searle, L. Maschio, L. Edith-Daga, A. Cossard, et al., J. Chem. Theor. Comput. (2022).
  * Maschio et al. (2015) L. Maschio, M. Rérat, B. Kirtman, and R. Dovesi, J. Chem. Phys. 143 (2015).
  * Rérat and Kirtman (2021) M. Rérat and B. Kirtman, J. Chem. Theor. Comput. 17, 4063 (2021).
  * (18) The agreement with the single and doubled cells and many k points, is always much closer for OR than for SHG (e.g. 47.8 vs. 47.9 with triple-zeta and 53.5 vs. 53.5 with quadruple-zeta basis sets for OR compared to 710.7 vs. 625.5 and 101.3 vs. 103.2 for SHG). This is, in part, because of our numerically more robust treatment of quasi-degeneracies in the calculation of $Q_{i,l}$ coefficients for OR, employing a third-order self-consistent Brillouin-Wigner perturbation theory approach, which will be presented elsewhere. On the other hand, for SHG, we simply use the unmodified algorithm of Maschio et al., in which quasi-degenerate $Q_{i,l}$ coefficients are set to zero.


#  The Superscaling Variable and Neutrino Energy Reconstruction,
From Theoretical Predictions to Experimental Limitations

D. Douqa  S. Bordoni  L. Giannessi  Université de Genève, Section de Physique,
DPNC, 1205 Genève, Switzerland  A. Nikolakopoulos  Theoretical Physics
Department, Fermilab, Batavia, IL 60510, USA  F. Sánchez  C. Schloesser
Université de Genève, Section de Physique, DPNC, 1205 Genève, Switzerland

###### Abstract

We introduce the novel approach of using the superscaling variable as an
observable and an analysis tool in the context of charged current neutrino-
nucleus interactions. We study the relation between the superscaling variable
and the removal energy, in addition to other fundamental parameters of the
neutrino-nucleus interaction models. In the second half of the paper, we
discuss the experimental viability of this measurement following a study of
neutrino energy and missing momentum reconstruction. We show that the
superscaling variable is measurable in neutrino interaction experiments
provided that the proton is detected in the final state. We discuss the
resolution of this measurement, and the limitation imposed by the proton’s
detection threshold.

##  I Introduction

The studies of neutrino-nucleus interactions are entering a new stage,
motivated by long-baseline experimental programs, in which the statistical
uncertainties will diminish and thus nuclear effects – contributing to the
systematic error – have to be kept well under control [1]. The incomplete
theoretical knowledge of neutrino-nucleus interactions influences various
stages of experimental analyses. For instance, for the future Hyper-Kamiokande
water Cherenkov detector [2], the method for reconstructing the neutrino
energy will be mainly based on the kinematics of the outgoing muon, which is
the only particle observed, assuming that the reaction-mechanism is two-body
quasi-elastic (QE) scattering on a bound nucleon.

In the context of CCQE interactions, the superscaling approach has been
advocated to enhance some of the most advanced models of neutrino-nucleus
interactions [3]. Superscaling was introduced initially to describe electron-
nucleon scattering where it is a very well established feature observed in the
data [4, 5]. The experimental determination of the superscaling behaviour of
neutrino interactions will be fundamental to the development of the
aforementioned models.

There are many uncertainties in the modelling of neutrino-nucleus interactions
[1], among them the description of the initial state nucleus, dynamical
phenomena such as Random Phase Approximations (RPA) or 2p2h interactions,
Pauli blocking and nuclear re-scattering. The removal energy —energy absorbed
by the unobserved nuclear system during the interaction— has emerged among
them as one of the largest systematic errors on T2K neutrino oscillations [6].
A direct measurement of the average removal energy will be very beneficial to
reduce systematic errors and to add credibility to neutrino-nucleus
interaction models. However, in flux-averaged data, the removal energy is not
directly observable as the neutrino energy is unknown on an event-to-event
basis. We present the reconstructed superscaling variable as a potential
observable to address this uncertainty experimentally.

We introduce theoretical predictions and discuss the relation between the
superscaling variable and the removal energy, employing a range of neutrino-
nucleus interaction models. These models span from the conventional
Relativistic Fermi Gas (RFG), for which the superscaling variable was
originally formulated, to more sophisticated Mean Field approaches. We perform
a study under the conditions of the T2K experiment [7], using the so-called
topological event selection. Topologies are based on the counting of particles
emerging from the nucleus after nuclear re-scattering (NrS). In this study we
concentrate on the CC0$\pi$1p topology, with a muon, a proton and no mesons in
the final state. We show that a measurement of the reconstructed superscaling
variable provides sensitivity to the removal energy in a $1$-nucleon knockout
event (CC1p1h).

##  II Event Kinematics and Neutrino Energy Reconstruction

In this section we discuss the kinematics of the CC1p1h reaction, and the
neutrino energy reconstruction in a $1\mu 1p$ event. We use several nuclear
models for this study, and expose the differences between their treatments.

We consider the semi-inclusive reaction, in which a proton and a muon are
detected in coincidence

$\nu_{\mu}(k_{\nu})+A(k_{A})\rightarrow\mu(k_{\mu})+p(k_{p})+B(k_{B})$ (1)

where $A$ is the initial nucleus, and $B$ is the residual hadronic system,
which remains undetected. Particle four-vectors are labelled between
parentheses, where
$k_{i}=(E_{i}=\sqrt{M_{i}^{2}+\lvert\vec{k}_{i}\rvert^{2}},\,\vec{k}_{i})$.

The energy and momentum transfer to the nucleus are $\omega=E_{\nu}-E_{\mu}$,
and $\vec{q}=\vec{k}_{\nu}-\vec{k}_{\mu}$, respectively.

Momentum and energy conservations are expressed as

$\displaystyle\vec{q}$ $\displaystyle=$ $\displaystyle\vec{k}_{p}+\vec{k}_{B}$
(2) $\displaystyle\omega-E_{p}$ $\displaystyle=$ $\displaystyle
M_{B}-M_{A}+\xi^{*}+T_{B}$ (3)

where $\xi^{*}$ is the excitation energy of the final nucleus, $M_{B}$ is its
ground state mass and $T_{B}$ is its kinetic energy. A further simplification
can be performed using the definition of the experimental nucleon separation
energy $S_{N}=M_{B}-M_{A}+M_{N}$, giving

$\omega-E_{p}+M_{N}=S_{N}+\xi^{*}+T_{B}$ (4)

The missing momentum is defined as $\vec{p}_{m}\equiv-\vec{k}_{B}$, while the
missing energy is defined as

$E_{m}=\omega-T_{p}-T_{B}=S_{N}+\xi^{*}+(M_{p}-M_{N})$ (5)

For what follows, it is useful to define

$\displaystyle\tilde{E}_{m}$ $\displaystyle=$ $\displaystyle
E_{m}-(M_{p}-M_{N})$ (6) $\displaystyle=$ $\displaystyle\omega-E_{p}+M_{N}$
(7) $\displaystyle\equiv$ $\displaystyle S_{N}+\xi^{*}$ (8)

where we neglect the kinetic energy of the recoiling nucleus as it becomes
small when the mass of the residual system is large with respect to its
momentum. The missing energy can be associated with the removal energy $S$,
the energy needed to remove a nucleon from the nucleus, and we will treat them
as equivalent ($S\equiv\tilde{E}_{m}$).

###  II.1  Monte Carlo Models

NEUT [8] provides a very good environment to predict experimental observables.
It includes a complete simulation of resonance production, shallow and deep
inelastic scattering, meson exchange currents and charged and neutral currents
both for neutrinos and anti-neutrinos, and for several nuclei. In addition, it
includes an intra-nuclear cascade model, developed in [9] and recently tuned
to available scattering data [10].

In order to study the model-dependence of our procedure, we use the following
interaction models implemented in NEUT for QE scattering.

  1. 1.

The Relativistic Fermi Gas (RFG) model, which uses the Smith-Moniz
parametrization [11].

  2. 2.

The Local Fermi Gas (LFG) model of Refs. [12, 13]. It includes corrections for
long-range correlations, final-state interactions, and the Coulomb potential
experienced by the charged lepton. The final-state nucleon kinematics are
generated from the LFG momentum distribution assuming quasi-free scattering,
as described and applied to neutrino scattering data in Ref. [14].

  3. 3.

The Spectral Function (SF) model. This is a calculation in the Plane Wave
Impulse Approximation (PWIA) that uses the realistic Rome spectral function
[15, 16]. The implementation is detailed in Ref. [17].

In addition to the NEUT native models, we include the Energy-Dependent
Relativistic Mean Field (EDRMF) [18] and Relativistic Plane Wave Impulse
Approximation (RPWIA) [19] models where NEUT’s intra-nuclear cascade was used
to simulate nuclear re-scattering and final state interactions (FSI) following
the same approach presented in [20].

Within the Fermi gas models, interaction with an initial on-shell nucleon with
momentum $\vec{k}_{N}$ is considered. In this case $\vec{p}_{m}=\vec{k}_{N}$
and

$\omega-E_{p}+M_{N}=T_{F}-T_{N}+S_{N}$ (9)

where $T_{F}$ is the kinetic energy of a nucleon at the Fermi level. One thus
has

$\tilde{E}_{m}=T_{F}-T_{N}+S_{N}$ (10)

where for the RFG model, $T_{F}$ is a constant.

On the other hand, in the LFG model, $T_{F}$ depends on the local nuclear
density, and an average over the density is performed. The value of $T_{F}$ is
hence smeared out, losing the direct relation between missing energy and
missing momentum of the RFG.

The RPWIA and EDRMF approaches [20, 19, 21] use a RMF model to describe the
initial state [22, 23]. In this case, a nucleon is knocked out from a nuclear
shell, which has a fixed missing energy. For carbon, considered in this work,
these are the $s_{1/2}$ and $p_{3/2}$ shells, with occupancies of two and four
nucleons respectively.

The SF calculation uses a more realistic energy and momentum density, by
implementing the experimental observation that the energy of the nuclear
shells is smeared out, and that the shell-model states are depopulated [24,
25, 26]. The strength missing from the mean-field appears at high missing
energy and momentum and is modelled in the local-density approximation [15,
16].

These approaches represent an increasingly realistic description of the
nuclear spectral function. And the EDRMF results, that use the same RMF
initial state as the RPWIA, are used to gauge the effect of nucleon final
state interactions. Consistent comparisons in the PWIA using these different
spectral functions, and detailed discussion on the kinematics, are presented
in Ref. [27].

###  II.2 Reconstruction of Full Kinematics

The full kinematics of the CC1p1h event, where the target nucleon is a neutron
($\nu_{\mu}+n\rightarrow\mu+p$), can be reconstructed using two different
approaches. When the missing energy and momentum are fully correlated, as in
the RFG model, we can perform reconstruction by energy-momentum conservation
(REP). In this case, Eq. 9 can be solved for the longitudinal missing
momentum, giving

$p_{m}^{L}=\frac{(p_{m}^{T})^{2}+M_{N}^{2}-(E_{\mu}+E_{p}-k_{\mu}^{L}-k_{p}^{L}+S_{\rm
REP})^{2}}{2(E_{\mu}+E_{p}-k_{\mu}^{L}-k_{p}^{L}+S_{\rm REP})},$ (11)

where the superscript $L$ denotes the components in the direction of the
neutrino beam, and

$p_{m}^{T}=(\vec{k}_{p}+\vec{k}_{\mu})^{T},$ (12)

is the missing momentum transverse to the beam. The neutrino energy can then
be obtained, using the fact that $E_{\nu}=|\vec{k}_{\nu}|$, from

$E_{\nu}=k_{\mu}^{L}+k_{p}^{L}-p_{m}^{L}$ (13)

Here, $S_{\rm REP}$ is a constant representing the average $\overline{\omega-
E_{p}+E_{N}}$ (or $\overline{T_{F}+S_{N}}$), which in the RFG model can be
identified as the fixed removal energy.

The distribution of $S_{\rm REP}$ is shown in the top panel of Fig. 1. It is
clear that while for the RFG model this leads to an exact reconstruction, this
is not the case for the other more realistic models that we consider.

Figure 1: The true distributions of $\omega-E_{p}+E_{N}$ (REP) and $\omega-
E_{p}+M_{N}$ (RE) for the five models for true CC1p1h events.

A more realistic approach, is to consider $\tilde{E}_{m}\approx constant$,
hence performing reconstruction solely by energy conservation (RE)

$E_{\nu}=E_{\mu}+E_{p}+S_{\rm RE}-M_{N}$ (14)

with $S_{\rm RE}$ being a constant, which corresponds to the average missing
(removal) energy $\overline{\omega-E_{p}+M_{N}}$ (or
$\overline{S_{N}+\xi^{*}}$).

This assumption considers $\tilde{E}_{m}$ to be independent from the missing
momentum, which is natural for the RPWIA, EDRMF and SF approaches in the low
missing (removal) energy region, where the bound-states do not satisfy a
dispersion relation based on free nucleon kinematics. The only dependence of
$\tilde{E}_{m}$ on the missing momentum comes from the recoil energy of the
residual nucleus. This is seen clearly in the bottom panel of Fig. 1, where
the tails of the peaks in the EDRMF and RPWIA calculations are due to the
recoil energy. In the LFG model, the average $\tilde{E}_{m}$ corresponds to
$\overline{T_{F}-T_{N}+S_{N}}$ from Eq. 10. Since the value of $T_{N}$ in the
LFG model peaks around the local value of $T_{F}$, there is a smaller
dispersion of the average $\overline{T_{F}-T_{N}+S_{N}}$ compared to
$\overline{T_{F}+S_{N}}$, which can be observed in Fig. 1.

This implies that the RE method would perform better in the reconstruction of
the neutrino energy for the SF, LFG, EDRMF and RPWIA models compared to the
REP method, the latter being only suitable for the RFG.

The average removal energy is computed from the averages of the distributions
in Fig. 1. Taking into consideration the optimal reconstruction method for
each model, the SF, LFG, EDRMF and RPWIA models have $S_{\rm RE}$ values of
26.45, 28.39, 29.02, and 29.20 MeV, respectively. On the other hand, the RFG
model has an $S_{\rm REP}$ value of 24.73 MeV.

To perform the neutrino energy reconstruction, we need to make an assumption
for a constant average removal energy. Based on the distributions in Fig. 1,
we set $S_{\rm RE}=S_{\rm REP}=28$ MeV for all models. We will depart from
this assumption later in the paper to explore the impact of changing the
average removal energy on the superscaling variable distribution.

The results for the reconstructed neutrino energy, using the RE and REP
methods are shown in Fig. 2. For the RE method (bottom panel), the
distributions are mainly contained between $-4\%$ and $2\%$ for the LFG and SF
models, while the EDRMF and RPWIA have a distinct concentration of strength at
$-3\%$ due to the $s$-shell contribution. The REP method, shown in the top
panel, exhibits a much broader but very similar distribution for all these
models. However, the RFG is the outlier, with a significantly better
resolution in REP, as expected.

Figure 2: The reconstructed neutrino energy resolution via REP (top) and RE
(bottom) for true CC1p1h events.

Figure 3: The reconstructed missing momentum resolution via REP (top) and RE
(bottom).

The reconstructed missing momentum is shown in Fig. 3, again the RFG is the
outlier, while for the other models we obtain fairly symmetric distributions
for RE, with a full width at half maximum of around $5\%$.

It is clear that the RE method is more robust and provides better results for
most models. Therefore, we will be using the RE method for all models, keeping
in mind that the RFG model is a special case due to its simplicity.

Lastly, we note that these assumption will become less reliable when the
signal has a large contribution from processes other than one-nucleon
knockout. However, as we will show in section V, the bulk of these
contributions can be removed by a kinematic cut on the reconstructed missing
momentum.

##  III  The Superscaling variable

When an interaction between a particle and a many-body system involves energy
and momentum transfers only to individual constituents of the complex system,
the inclusive cross section can be approximated as a single-nucleon cross
section times a specific function of the energy and momentum transfers
$f(\omega,|\vec{q}|)$. Scaling, of the first kind, occurs when that function
becomes independent of both $\omega$ and $\vec{q}$ explicitly. The scaling
function depends on the kinematics only through a single quantity
$\psi(\omega,|\vec{q}|)$ [28, 29].

The superscaling variable $\psi^{\prime}$ was first introduced by Donnelly and
Sick [4] as an evolution of the scaling variable $\psi$ used by Alberico et
al. [30] and within the framework of the RFG model. It is defined as:

$\psi^{\prime}(\omega,\vec{q})=\frac{1}{\sqrt{\sqrt{1+\eta_{F}^{2}}-1}}\frac{\lambda-\tau}{\sqrt{(1+\lambda)\tau+\kappa\sqrt{\tau(1+\tau)}}}$
(15)

with

$\displaystyle\eta_{F}$ $\displaystyle=$ $\displaystyle\frac{k_{F}}{M_{N}}$
(16) $\displaystyle\kappa$ $\displaystyle=$
$\displaystyle\frac{|\vec{q}\,|}{2M_{N}}$ (17) $\displaystyle\lambda$
$\displaystyle=$ $\displaystyle\frac{\omega-E_{\rm\rm shift}}{2M_{N}}$ (18)
$\displaystyle\tau$ $\displaystyle=$ $\displaystyle\kappa^{2}-\lambda^{2}$
(19)

where, $M_{N}$ is the neutron mass, $k_{F}$ is the Fermi momentum (fixed to
228 MeV for Carbon [31]) and $E_{\rm shift}$ is a shift energy used to make
the quasi-elastic peak coincide with $\psi^{\prime}=0$. Presumably, $E_{\rm\rm
shift}$ includes information on the separation energy, the average removal
energy of nucleons in the nucleus as well as aspects of final-state
interactions like RPA, which can influence the removal energy distribution by
altering the cross-section strength as a function of the momentum transfer in
the reaction. In the following, we use the reported value for Carbon of 20 MeV
[31].

The discussion in the following sections relies on the analysis of the special
case when $\psi^{\prime}=0$. This condition is achieved when:

$\displaystyle\lambda-\tau=\lambda-(\kappa^{2}-\lambda^{2})=0$ (20)

The $\lambda$ parameter has a direct dependency with $E_{\rm shift}$ and the
removal energy $S$—assumed to be equivalent to the missing energy
$\tilde{E}_{m}$ in section II—through the energy transfer $\omega$. From Eq.
7:

$\displaystyle\omega=E_{\nu}-E_{\mu}=E_{p}+S-M_{N}$ (21)

By solving Eq. 20 and considering the positive solution as the physically
meaningful one, we arrive at the following relationship

$\displaystyle S-E_{\rm shift}=\sqrt{M_{N}^{2}+|\vec{q}|^{2}}-E_{p}$ (22)

For vanishing missing momentum, Eq. 22 becomes

$\begin{split}S-E_{\rm shift}&=\sqrt{M_{N}^{2}+|\vec{k}_{p}|^{2}}-E_{p}\\\
&\approx\frac{1}{2}\frac{M_{N}^{2}-M_{p}^{2}}{E_{p}}\end{split}$ (23)

where it is clear that $E_{\rm shift}$ approaches the average value of $S$
(previously denoted as $S_{\rm RE}$) for $|\vec{p}_{m}|=0$, with a maximum
difference of 1.29 MeV, attained when the proton is at rest.

On the other hand, for $|\vec{p}_{m}|>0$, one obtains

$\displaystyle S-E_{\rm
shift}\approx\frac{p_{m}^{2}-2\vec{k}_{p}\cdot\vec{p}_{m}}{2E_{p}}$ (24)

where the equation remains more complex with dependencies on the missing and
final nucleon momenta through the residual dependency on $|\vec{q}|$.

As such, $E_{\rm shift}$ should be a good estimator, within a few MeV, of the
average removal energy, at low missing momenta.

Note that in the analysis of inclusive electron scattering data [31], the
value of $E_{\rm shift}$ is determined by requiring that $\psi^{\prime}=0$ at
the experimentally observed quasi-elastic peak. In the RFG, the peak of the
cross section is indeed obtained for $\psi^{\prime}=0$ when $E_{\rm shift}$ is
equal to the separation energy used in the model. The presence of final-state
interactions, RPA effects, and Pauli blocking among others, can break this
simple relation. In any case, the peak position of $\psi^{\prime}$ is still
strongly correlated with the average removal energy. In the following section,
we study in detail the distribution of the superscaling variable, and its
dependence on kinematic variables.

##  IV True Monte Carlo Predictions

In this section, we perform a $\psi^{\prime}$ calculation and characterization
using samples of true CC1p1h events, without including nuclear re-scattering
(NrS).

Figure 4 shows the distribution of the superscaling variable from the event
kinematics of the five models, where $\psi^{\prime}$ is calculated using the
true $\omega$ and $\vec{q}$ and with a fixed $E_{\rm shift}$ of 20 MeV. The
estimated peak position of $\psi^{\prime}$ for each model is obtained by
taking the mean of a Gaussian fit to the core of the $\psi^{\prime}$
distribution (maximum-1$\sigma$, maximum+1$\sigma$).

Figure 4: The true superscaling variable for true CC1p1h events, for the five
models. The legend shows the peak positions for each model.

Since the definition of $\psi^{\prime}$ (Eq. 15) is derived from the RFG
model, the distribution of $\psi^{\prime}$ in the RFG follows the expected
behaviour with a peak at zero contained between $\pm 1$. The deviation from
the assumptions of RFG, particularly when it comes to the momentum of the
target nucleon and removal energy distribution, creates the positive tails
observed in the other models.

Further details can be inferred from the relation between $\psi^{\prime}$ and
the true missing momentum, equivalent to the momentum of the target neutron in
the case of LFG, as shown in Fig. 5. According to Eq. 23, in the absence of
Fermi momentum ($p_{m}=0$ MeV/c) the condition for $\psi^{\prime}=0$ entails
that the difference $S-E_{\rm shift}$ can be up to 1.29 MeV. When using the
true $\psi^{\prime}$, we do not make any assumptions on the average removal
energy, but we fix $E_{\rm shift}$ to 20 MeV. This gives a narrow
$\psi^{\prime}$ distribution around zero, as seen in Fig. 6, with a width that
is correlated to the width of the removal energy distribution (Fig. 1, bottom
plot). The small displacement of the peak position from zero for low $p_{m}$
can be attributed to the difference $S-E_{\rm shift}$ being larger than 1.29
MeV when using $E_{\rm shift}=20$ MeV while the average removal energy for LFG
is 28.39 MeV.

Figure 5: Distribution of events as a function of the superscaling variable
and the missing momentum. The width of $\psi^{\prime}$ increases for larger
values of $p_{m}$. Figure 6: Distribution of the superscaling variable in
different slices of missing momentum.

When taking slices of higher $p_{m}$, as shown in Fig. 6, the distribution of
$\psi^{\prime}$ flattens and does not exhibit a clear peak anymore. This
occurs because in the LFG the separation energy is not a single fixed number
for each $p_{m}$, and its dependence on $p_{m}$ is different than what is
found in the RFG.

Figure 7 shows the event distribution of the LFG model in terms of the true
$\psi^{\prime}$ and the removal energy defined in Eq. 7 and calculated using
true event kinematics. The plot shows that $\psi^{\prime}$ goes to zero when
the removal energy is around 20 MeV (equal to $E_{\rm shift}$). For large
values of $\omega-E_{p}+M_{N}$, there is a bias towards high values of
$\psi^{\prime}$ which produces an asymmetric distribution, as seen in Fig. 4,
and shifts the peak position towards larger values with a maximum shift of
0.25. The shift in peak position, as well as the increasing asymmetry of
$\psi^{\prime}$ for increasing slices of removal energies can seen in Fig. 8.
For events in these regions, a larger $E_{\rm shift}$ than the one used would
be required to shift the QE peak to $\psi^{\prime}=0$.

The widening of the $\psi^{\prime}$ distribution with increasing removal
energy is correlated to the missing momentum through the model-dependent
relation between $S$ and $p_{m}$. In the SF model, all values of $p_{m}$ are
accessible in each slice of $S$, and therefore this widening effect is not
present.

Figure 7: The true superscaling variable as a function of the removal energy
defined as $\omega-E_{p}+M_{N}$ for the LFG model. The points represent the
peak position of $\psi^{\prime}$ for each bin of the removal energy. Figure 8:
The true $\psi^{\prime}$ distribution for different slices of the removal
energy for the LFG model. Figure 9: The peak positions of the true
$\psi^{\prime}$ distribution as a function of 100 MeV/c slices of the true
momentum transfer $|\vec{q}|$ for the five models.

Another source of difference between the models is the relation between the
distribution of $\psi^{\prime}$ and the momentum transfer, $\vec{q}$. In Fig.
9 we take 100 MeV/c slices of the true $|\vec{q}|$ and plot the
$\psi^{\prime}$ peak position as a function of the mid point for each slice.
In the SF, LFG, RFG and RPWIA models, the peak position of the distribution of
$\psi^{\prime}$ is seen to have a similar dependence on $|\vec{q}|$, namely a
positive shift in the peak position at low-$|\vec{q}|$. The EDRMF model is the
outlier here, the peak position shows little dependence on the momentum
transfer. This is to be attributed to the final-state potential included in
this calculation. Indeed, a similar difference in the peak position of the
scaling function in the RPWIA and EDRMF approaches as function of $|\vec{q}|$
can be seen for electron-scattering cross sections in Ref. [32].

Finally, we performed a search for the $E_{\rm shift}$ value that gives a
$\psi^{\prime}$ distribution peaking at zero for $E_{\rm shift}$ between 20
and 50 MeV. Figure 10 shows that the relation between the peak position and
$E_{\rm shift}$ is approximately linear, with a nearly universal slope for the
models, but differing intercepts.

Figure 10: The peak positions of each $\psi^{\prime}$ distribution as a
function of the $E_{\rm shift}$ value used in its calculation for the five
models.

The $E_{\rm shift}$ values that give a $\psi^{\prime}$ distribution peaking at
zero for each model are listed in Tab. 1. As expected, the RFG model arrives
at $\psi^{\prime}=0$ at an $E_{\rm shift}$ value closest to the 20 MeV used in
the calculation [31], while the other models have different results based on
their individual removal energy distributions. Comparing these values with the
average removal energy (the means in Fig. 1), we observe that they are within
a few MeVs, except for the RFG model where $\overline{\omega-E_{p}+E_{N}}$ is
more consistent with the model’s assumptions, and the EDRMF model where the
final-state potential is known to shift the position of the quasi-elastic peak
to smaller energy transfer, in line with electron scattering data [18]. This
means that the distribution peaks at lower $\psi^{\prime}$, as seen in Fig. 4.

Table 1: The values of $E_{\rm shift}$ that give a $\psi^{\prime}_{\rm true}$ distribution peaking at zero for each model, compared to the average removal energy defined as $\overline{\omega-E_{p}+M_{N}}$. Model |  $E_{\rm shift}|_{\psi^{\prime}=0}$ [MeV] |  $\overline{\omega-E_{p}+M_{N}}$ [MeV]
---|---|---
SF | 26.19 | 26.45
LFG | 24.07 | 28.39
RFG | 20.62 | 9.91
EDRMF | 16.16 | 29.02
RPWIA | 23.95 | 29.20

This measurement of $E_{\rm shift}$ relies on using true Monte Carlo (MC)
information to obtain the true $\psi^{\prime}$ distribution. This is not
achievable experimentally since energy and momentum transfer are not
observable in a neutrino experiment. However, by taking a fixed average value
of the removal energy ($S$), we can reconstruct the full kinematics of a
1$\mu$1p event, and obtain an approximate superscaling variable. We will later
perform another measurement of $E_{\rm shift}$ using the reconstructed
$\psi^{\prime}$ and under more experimental-like conditions.

##  V  Experimental-like Event Selection

For a more realistic sample of events from an experimental point of view, we
use the selection criteria for the CC0$\pi 1p$ topology as done for the T2K
experiment [33]. This selection requires the detection of a muon, and a single
proton in the final state, with $|\vec{k_{p}}|>0.45$ GeV and
$\cos{\theta_{p}}>0.4$ to reflect the proton detection threshold and detector
acceptance. Events with multiple protons, or pions, exiting the nucleus are
rejected.

This selection is applied to events following nuclear re-scattering, where we
use the intra-nuclear cascade model implemented in NEUT for all five models.

For the NEUT models (SF, LFG and RFG), we distinguish five different samples
that pass our selection according to their production mechanism:

  1. 1.

CC1p1h, charged current one particle one hole, when the proton emerges from
the nucleus with no NrS.

  2. 2.

CC1p1h+NrS, charged current one particle one hole, when the proton undergoes
NrS.

  3. 3.

CC2p2h, charged current two particles two holes where the second proton is
undetectable.

  4. 4.

CCRes, charged current resonance, where the pion is absorbed by the nucleus.

  5. 5.

CCDis, charged current deep inelastic scattering, where all pions are absorbed
by the nucleus.

In the case of EDRMF and RPWIA, only CC1p1h and CC1p1h+NrS are included as
they are not native NEUT models.

The modeling of the selection and the prediction from the LFG model in NEUT
are in very good agreement with T2K and MINERvA experimental results as shown
in [14].

Additionally, we also introduce a cut on the reconstructed missing momentum as
proposed in [17] and the first experimental applications described in [34].
When including background events that pass the CC0$\pi 1p$ selection, the
reconstructed missing momentum shows a very distinctive shape with CC1p1h
events concentrated below the Fermi level (at around 300 MeV/c), see Fig. 11.
This variable, thus, provides a very efficient cut to reduce the contamination
of NrS, CC2p2h and CCRes and increase the purity of CC1p1h events in the
selected sample.

Table 2 shows the detailed sample composition for the five models before and
after introducing the cut at $p_{m}^{\rm RE}=300$ MeV/c.

Figure 11: A stacked histogram of the reconstructed missing momentum via RE
for the LFG model. Table 2: The fraction of each interaction type in the
original sample of events (top row) compared to the sample after applying the
cut on the missing momentum (bottom row) for each model.

Model | CC1p1h | CC1p1h+NrS | CC2p2h | CCRes | CCDis
---|---|---|---|---|---
SF | 0.575 | 0.143 | 0.148 | 0.127 | 0.007
0.833 | 0.076 | 0.032 | 0.057 | 0.001
LFG | 0.560 | 0.171 | 0.097 | 0.125 | 0.007
0.839 | 0.085 | 0.023 | 0.052 | 0.001
RFG | 0.636 | 0.164 | 0.085 | 0.109 | 0.006
0.860 | 0.075 | 0.019 | 0.044 | 0.001
EDRMF | 0.764 | 0.236 | – | – | –
0.897 | 0.103 | – | – | –
RPWIA | 0.770 | 0.230 | – | – | –
0.901 | 0.099 | – | – | –

##  VI Experimental Observables and the Reconstructed Superscaling Variable

In this section we take a practical point of view, and discuss a measurement
of the reconstructed superscaling variable, based on the energy reconstruction
in 1$\mu$1p events presented in section II.2. We make use only of observables
accessible in neutrino detectors, and include re-scattering in the intra-
nuclear cascade model as discussed in the previous section.

###  VI.1 The Superscaling Variable

The reconstructed superscaling variable ($\psi^{\prime}_{\rm RE}$) can be
calculated from Eq. 15 following the RE energy reconstruction to obtain the
energy and momentum transfer. Figure 12 shows the $\psi^{\prime}_{\rm RE}$
distribution for the LFG model as a function of the true value for all events
passing the CC$0\pi 1p$ selection. We observe strong correlation between the
true and reconstructed values corresponding to CC1p1h events (distinguished in
red). Under the CC1p1h line, the correlation disappears for background events
(in blue) as the reconstruction does not account for their energy losses
resulting in the wrong energy-momentum transfer estimation and leading to an
underestimated $\psi^{\prime}_{\rm RE}$ calculation. This can shift the peak
position of $\psi^{\prime}_{\rm RE}$ and produce negative tails. However,
imposing the cut on the reconstructed missing momentum can be used to limit
these events. The distribution of $\psi^{\prime}_{\rm RE}$ before and after
the cut for the LFG model is shown in Fig. 13, where the purity of CC1p1h (no
NrS) events has increased from 0.56 to 0.84. We will apply this cut to the
selected sample of events in what follows.

Figure 12: Reconstructed vs true $\psi^{\prime}$ for the LFG model, where
event kinematics were reconstructed via RE and no cuts were applied to
minimize the background. True CC1p1h events (red) fall on the diagonal, while
all other events (blue) have lower reconstructed $\psi^{\prime}$ compared to
the true value. Figure 13: The reconstructed superscaling variable for the LFG
model. Applying the cut $p_{m}^{\rm RE}<300$ MeV/c eliminates a good fraction
of the background reducing the positive tails. Figure 14:  The reconstructed
energy transfer $\omega$ versus the reconstructed momentum transfer
$|\vec{q}|$. Overlaid are the contour lines representing
iso-$\psi^{\prime}_{\rm RE}$ values. Figure 15: The reconstructed superscaling
variable for the five models. The cut on $p_{m}^{\rm RE}$ was applied to
minimize the background.

The values of the reconstructed energy and momentum transfers used in the
calculation of $\psi^{\prime}_{\rm RE}$ are shown in Fig. 14. The cutoff at
low $\omega$ is a consequence of the requirement of a 450 MeV/c proton in the
final state, which corresponds to a kinetic energy of about 140 MeV. The non-
linearity of the iso-$\psi^{\prime}_{\rm RE}$ lines means that negative values
of the superscaling variable are more affected by this limitation on the
proton momentum than positive ones, which introduces another source of
asymmetric in the distribution of $\psi^{\prime}_{\rm RE}$. Another
consequence of the acceptance cut on the proton momentum is that it eliminates
a region of phase space that particularly highlights differences between the
models in their $\vec{q}$ dependence (low $|\vec{q}|$ in Fig. 9) as well as in
their treatment of Pauli blocking. This is a delicate experimental issue that
will be partially circumvented by the new generation of near detectors [35]
that reduce the proton detectability threshold.

Using the RE method, and the aforementioned selection and $p_{m}^{\rm RE}$ cut
on all models, Fig. 15 shows the $\psi^{\prime}_{\rm RE}$ distributions for
the five considered models. Compared to the true distributions and peak
positions (Fig. 4), all reconstructed distributions exhibit the positive
shifts discussed above. The average shift in the peak position for all five
models is 0.068, with maximum shifts observed in the EDRMF and RFG models.

In the case of the RFG, the large difference between true and reconstructed
$\psi^{\prime}$ is mainly due to the energy reconstruction. On the other hand,
the final state potential of the EDRMF model—previously shown to cause
$\psi^{\prime}$ to peak at lower values in section IV—has a diminishing effect
when the proton detection threshold is applied and the low $|\vec{q}|$ region
is eliminated. This causes the peak position of the EDRMF model to align with
that of the RPWIA.

In the true Monte Carlo predictions we observed that the distribution of
$\psi^{\prime}$ loses its distinguished peak and becomes more flat when taking
slices of increasing missing momenta. This effect can also be observed in the
reconstructed case, see Fig. 16, and would affect the shape of the
$\psi^{\prime}_{\rm RE}$ distribution for experimental conditions where only a
certain part of the kinematical phase space is accessible. A good
reconstruction of the missing momentum allows for identifying the regions with
low $p_{m}$ where the $\psi^{\prime}$ distribution is more narrow and the peak
position can be determined with less bias.

Figure 16: The reconstructed superscaling variable for different slices of the
reconstructed missing momentum. $\psi^{\prime}_{\rm RE}$ becomes more flat for
increasing $p_{m}^{\rm RE}$.

###  VI.2 The Removal and Shift Energies

Figure 17: The peak position of the reconstructed $\psi^{\prime}$ as a
function of $E_{\rm shift}$. The average removal energy is fixed to 28 MeV in
the energy reconstruction for all models.

Using $\psi^{\prime}_{\rm RE}$, we can perform another search for the $E_{\rm
shift}$ values required to shift the peak position to zero for each model.
Figure 17 shows the peak positions of $\psi^{\prime}_{\rm RE}$ for each value
of $E_{\rm shift}$. Compared to Fig. 10, the slopes for all five models are
now closer in value. This is a consequence of the detectability threshold of
the proton momentum, which eliminates the low $|\vec{q}|$ region, minimizing
the differences between models and leading to more similar slope values.

To distinguish between the effect of the energy reconstruction and that of the
sample selection and detector acceptance cuts, we perform the same exercise
using the reconstructed energy and momentum transfer but on true CC1p1h
events. Table 3 shows the $E_{\rm shift}$ values required to shift the
$\psi^{\prime}_{\rm RE}$ peak to coincide with zero for the experimental-like
event sample (left) compared to the true CC1p1h event sample (center). On the
right, we also include the results obtained in Tab. 1 for the true
$\psi^{\prime}$.

Table 3: The values of $E_{\rm shift}$ that give a $\psi^{\prime}_{\rm RE}$
peaking at zero for a realistic sample of CC0$\pi$1p events (with the
$p_{m}^{\rm RE}$ and proton acceptance cuts) and true CC1p1h events. Compared
to the $E_{\rm shift}$ that give a $\psi^{\prime}_{\rm true}$ peaking at zero.

Model |  $E_{\rm shift}|_{\psi^{\prime}_{\rm RE}=0}$ [MeV] |  $E_{\rm shift}|_{\psi^{\prime}_{\rm true}=0}$ [MeV]
---|---|---
| CC0$\pi$1p* | CC1p1h | CC1p1h
SF | 36.64 | 25.28 | 26.19
LFG | 31.94 | 25.59 | 24.07
RFG | 38.34 | 25.37 | 20.62
EDRMF | 34.93 | 16.24 | 16.16
RPWIA | 34.23 | 21.88 | 23.95
*CC0$\pi$1p with proton acceptance cuts and $p_{m}^{\rm RE}<$300 MeV/c.

We observe that larger $E_{\rm shift}$ values are required for
$\psi^{\prime}_{\rm RE}$ in the experimental-like event sample. This is to
compensate for the effect observed in Fig. 14 where the proton acceptance cut
eliminates more events on the negative side of $\psi^{\prime}$ compared to the
positive, as well as the inclusion of NrS and CC0$\pi$1p background events
that occupy the positive side of $\psi^{\prime}$.

On the other hand, for true CC1p1h events without any cuts, we obtain similar
values of $E_{\rm shift}$ for the true and reconstructed $\psi^{\prime}$. This
means that the energy reconstruction and our assumption on the average removal
energy has a small effect on the $\psi^{\prime}_{\rm RE}$ peak position and
the required $E_{\rm shift}$, and the error is dominated by the proton’s
momentum detection threshold and background contamination.

##  VII  Conclusions

We presented the superscaling variable as an observable in charged current
neutrino-nucleus interactions, comparing the shape and characteristics of the
$\psi^{\prime}$ distribution for five neutrino-nucleus interaction models.
Using true Monte Carlo information, we studied the superscaling variable
dependencies with fundamental parameters of the models such as the removal
energy, the target nucleon momentum and the momentum transfer of the
interaction.

We have shown that $\psi^{\prime}$, and particularly the peak position of the
distribution, provides information on the removal energy of nucleons. In the
peak region, where $\psi^{\prime}=0$, the average removal energy is related to
the shift energy, $E_{\rm shift}$, required to shift the quasi-elastic peak to
zero, at low missing momentum

To comment on the determination of $\psi^{\prime}$ in neutrino interaction
experiments where the neutrino energy is not directly measurable, we
introduced two neutrino energy reconstruction methods built on different
assumptions for the average removal energy. We concluded that reconstruction
by energy-momentum conservation (REP), where the missing energy and momentum
are fully correlated, performs better for the RFG model due to its simplicity,
while reconstruction by energy conservation (RE), where the missing (removal)
energy is assumed to be constant, has the advantage in all other models.

Applying the RE method on a more realistic event sample that passes the
CC0$\pi$1p selection, and utilizing a cut on the missing momentum to increase
the purity of CC1p1h events, we were able to obtain the reconstructed
superscaling variable with an average displacement of 0.068 for the peak
position with respect to the true distributions among the five models.

Finally, we have shown that while the energy reconstruction is sufficient to
measure $\psi^{\prime}$ experimentally, the limitations of current detector
technologies in measuring low momentum protons impose a bias on the
$\psi^{\prime}$ distribution requiring a larger $E_{\rm shift}$ to make the
peak occur at $\psi^{\prime}=0$. Even though it is difficult to estimate the
true average removal energy, the procedure outlined here—due to the simplicity
of the energy reconstruction method and robustness of the results—allows to
make a novel model-independent measurement which may be used to constrain or
distinguish interaction models. Additionally, future and upgraded detectors
are expected to have lower detection thresholds that allow for more accurate
measurements [35].

###### Acknowledgements.

This work was supported by the Swiss National Foundation Grant No.
200021_85012 and Fermi Research Alliance, LLC under Contract No. DE-
AC02-07CH11359 with the U.S. Department of Energy.

## References

  * Alvarez-Ruso  _et al._ [2018] L. Alvarez-Ruso _et al._ , NuSTEC White Paper: Status and challenges of neutrino–nucleus scattering, Prog. Part. Nucl. Phys. 100, 1 (2018).
  * Abe  _et al._ [2018a] K. Abe _et al._ (Hyper-Kamiokande), Hyper-Kamiokande Design Report, (2018a), arXiv:1805.04163 [physics.ins-det] .
  * Gonzaléz-Jiménez  _et al._ [2014] R. Gonzaléz-Jiménez, G. Megias, M. Barbaro, J. Caballero, and T. Donnelly, Extensions of Superscaling from Relativistic Mean Field Theory: the SuSAv2 Model, Phys. Rev. C 90, 035501 (2014).
  * Donnelly and Sick [1999] T. Donnelly and I. Sick, Superscaling of inclusive electron scattering from nuclei, Phys. Rev. C 60, 065502 (1999).
  * Barbaro  _et al._ [1998] M. Barbaro, R. Cenni, A. De Pace, T. Donnelly, and A. Molinari, Relativistic y - scaling and the Coulomb sum rule in nuclei, Nucl. Phys. A 643, 137 (1998).
  * Abe  _et al._ [2020] K. Abe _et al._ (T2K), Constraint on the matter–antimatter symmetry-violating phase in neutrino oscillations, Nature 580, 339 (2020).
  * Abe  _et al._ [2011] K. Abe _et al._ , The t2k experiment, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 659, 106–135 (2011).
  * Hayato and Pickering [2021] Y. Hayato and L. Pickering, The NEUT neutrino interaction simulation program library, The European Physical Journal Special Topics 230, 4469 (2021).
  * Salcedo  _et al._ [1988] L. L. Salcedo, E. Oset, M. J. Vicente-Vacas, and C. Garcia-Recio, Computer Simulation of Inclusive Pion Nuclear Reactions, Nucl. Phys. A484, 557 (1988).
  * Pinzon Guerra  _et al._ [2019] E. S. Pinzon Guerra _et al._ , Using world charged $\pi^{\pm}-$nucleus scattering data to constrain an intranuclear cascade model, Phys. Rev. D99, 052007 (2019).
  * Smith and Moniz [1972] R. Smith and E. Moniz, Neutrino reactions on nuclear targets, Nuclear Physics B 43, 605 (1972).
  * Nieves  _et al._ [2004] J. Nieves, J. E. Amaro, and M. Valverde, Inclusive quasi-elastic neutrino reactions, Phys. Rev. C 70, 055503 (2004), [Erratum: Phys.Rev.C 72, 019902 (2005)].
  * Gran  _et al._ [2013] R. Gran, J. Nieves, F. Sanchez, and M. J. Vicente Vacas, Neutrino-nucleus quasi-elastic and 2p2h interactions up to 10 GeV, Phys. Rev. D88, 113007 (2013).
  * Bourguille  _et al._ [2021] B. Bourguille, J. Nieves, and F. Sánchez, Inclusive and exclusive neutrino-nucleus cross sections and the reconstruction of the interaction kinematics, Journal of High Energy Physics 2021, 10.1007/jhep04(2021)004 (2021).
  * Benhar  _et al._ [1994] O. Benhar, A. Fabrocini, S. Fantoni, and I. Sick, Spectral function of finite nuclei and scattering of gev electrons, Nuclear Physics A 579, 493 (1994).
  * Benhar  _et al._ [2005] O. Benhar, N. Farina, H. Nakamura, M. Sakuda, and R. Seki, Electron- and neutrino-nucleus scattering in the impulse approximation regime, Phys. Rev. D 72, 053005 (2005).
  * Furmanski and Sobczyk [2017] A. P. Furmanski and J. T. Sobczyk, Neutrino energy reconstruction from one muon and one proton events, Phys. Rev. C 95, 065501 (2017).
  * González-Jiménez  _et al._ [2019] R. González-Jiménez, A. Nikolakopoulos, N. Jachowicz, and J. M. Udías, Nuclear effects in electron-nucleus and neutrino-nucleus scattering within a relativistic quantum mechanical framework, Physical Review C 100, 10.1103/physrevc.100.045501 (2019).
  * González-Jiménez  _et al._ [2022] R. González-Jiménez, M. B. Barbaro, J. A. Caballero, T. W. Donnelly, N. Jachowicz, G. D. Megias, K. Niewczas, A. Nikolakopoulos, J. W. Van Orden, and J. M. Udías, Neutrino energy reconstruction from semi-inclusive samples, Phys. Rev. C 105, 025502 (2022).
  * Nikolakopoulos  _et al._ [2022] A. Nikolakopoulos, R. González-Jiménez, N. Jachowicz, K. Niewczas, F. Sánchez, and J. M. Udías, Benchmarking intranuclear cascade models for neutrino scattering with relativistic optical potentials, Phys. Rev. C 105, 054603 (2022).
  * Franco-Patino  _et al._ [2022] J. M. Franco-Patino, R. González-Jiménez, S. Dolan, M. B. Barbaro, J. A. Caballero, G. D. Megias, and J. M. Udias, Final state interactions in semi-inclusive neutrino-nucleus scattering: Applications to the T2K and MINER$\nu$A experiments, Phys. Rev. D 106, 113005 (2022).
  * Sharma  _et al._ [1993] M. Sharma, M. Nagarajan, and P. Ring, Rho meson coupling in the relativistic mean field theory and description of exotic nuclei, Physics Letters B 312, 377 (1993).
  * Horowitz and Serot [1981] C. Horowitz and B. D. Serot, Self-consistent hartree description of finite nuclei in a relativistic quantum field theory, Nuclear Physics A 368, 503 (1981).
  * Benhar  _et al._ [1993] O. Benhar, V. R. Pandharipande, and S. C. Pieper, Electron-scattering studies of correlations in nuclei, Rev. Mod. Phys. 65, 817 (1993).
  * Kelly [2005] J. J. Kelly, Relativistic distorted-wave impulse approximation analysis of ${}^{12}\mathrm{C}(e,{e}^{{}^{\prime}}p)$ for ${Q}^{2}<2(\mathrm{GeV}/c){}^{2}$, Phys. Rev. C 71, 064610 (2005).
  * Dutta  _et al._ [2003] D. Dutta _et al._ , Quasielastic $(e,{e}^{{}^{\prime}}p)$ reaction on ${}^{12}\mathrm{C},^{56}\mathrm{Fe}$, and ${}^{197}\mathrm{Au}$, Phys. Rev. C 68, 064603 (2003).
  * Van Orden and Donnelly [2019] J. W. Van Orden and T. W. Donnelly, Nuclear theory and event generators for charge-changing neutrino reactions, Phys. Rev. C 100, 044620 (2019).
  * Day  _et al._ [1990] D. B. Day, J. S. McCarthy, T. W. Donnelly, and I. Sick, Scaling in inclusive electron-nucleus scattering, Annual Review of Nuclear and Particle Science 40, 357 (1990).
  * Amaro  _et al._ [2020] J. E. Amaro, M. B. Barbaro, J. A. Caballero, R. González-Jiménez, G. D. Megias, and I. R. Simo, Electron- versus neutrino-nucleus scattering, Journal of Physics G: Nuclear and Particle Physics 47, 124001 (2020).
  * Alberico  _et al._ [1988] W. M. Alberico, A. Molinari, T. W. Donnelly, E. L. Kronenberg, and J. W. Van Orden, Scaling in electron scattering from a relativistic fermi gas, Phys. Rev. C 38, 1801 (1988).
  * Maieron  _et al._ [2002] C. Maieron, T. W. Donnelly, and I. Sick, Extended superscaling of electron scattering from nuclei, Phys. Rev. C 65, 025502 (2002).
  * González-Jiménez  _et al._ [2020] R. González-Jiménez, M. B. Barbaro, J. A. Caballero, T. W. Donnelly, N. Jachowicz, G. D. Megias, K. Niewczas, A. Nikolakopoulos, and J. M. Udías, Constraints in modeling the quasielastic response in inclusive lepton-nucleus scattering, Phys. Rev. C 101, 015503 (2020).
  * Abe  _et al._ [2018b] K. Abe _et al._ (T2K), Characterization of nuclear effects in muon-neutrino scattering on hydrocarbon with a measurement of final-state kinematics and correlations in charged-current pionless interactions at T2K, Phys. Rev. D 98, 032003 (2018b).
  * Abratenko  _et al._ [2023] P. Abratenko _et al._ (MicroBooNE), Measurement of nuclear effects in neutrino-argon interactions using generalized kinematic imbalance variables with the MicroBooNE detector, (2023), arXiv:2310.06082 [nucl-ex] .
  * Dolan  _et al._ [2022] S. Dolan _et al._ , Sensitivity of the upgraded T2K Near Detector to constrain neutrino and antineutrino interactions with no mesons in the final state by exploiting nucleon-lepton correlations, Phys. Rev. D 105, 032010 (2022).