[a]Rajan Gupta

# Isovector Axial Charge and Form Factors of Nucleons from Lattice QCD

###### Abstract

I present an overview of the calculations of the isovector axial vector form
factor of the nucleon, $G_{A}(Q^{2})$, using lattice QCD. Based on a
comparison of results from various collaborations, a case is made that lattice
results are now consistent within 10%. A similar level of uncertainty is found
also in the axial charge $g_{A}^{u-d}$, the mean squared axial charge radius,
$\langle r_{A}^{2}\rangle$, the induced pseudoscalar charge $g_{P}^{\ast}$,
and the pion-nucleon coupling $g_{\pi NN}$. These lattice results for
$G_{A}(Q^{2})$ are already compatible with those obtained from the recent
MINER$\nu$A experiment but lie 2-3$\sigma$ higher than the phenomenological
extraction from the old $\nu$-deuterium bubble chamber scattering data for
$Q^{2}>0.3$ GeV2. Fits to our data show that the dipole ansatz does not have
enough parameters to parameterize the form factor over the range $0\leq
Q^{2}\leq 1$ GeV2, whereas even a $z^{2}$ truncation of the $z$-expansion or a
low order Padé are sufficient. Looking ahead, lattice QCD calculations will
provide increasingly precise results over the range $0\leq Q^{2}\lesssim 1$
GeV2, and MINER$\nu$A-like experiments will extend the range to $Q^{2}\sim 2$
GeV2 or higher. To increase precision of lattice data to the percent level,
new developments are needed to address two related issues: the exponentially
falling signal-to-noise ratio in all nucleon correlation functions and
removing excited state contributions. Nevertheless, even with the current
methodology, significant reduction in errors is expected over the next few
years with higher statistics data on more ensembles closer to the physical
point.

##  1 Introduction

The axial charge, $g_{A}^{u-d}$, gives the strength of the coupling of the
weak current to the nucleons. It has been determined very accurately from the
asymmetry parameter $A$ (relative to the plane defined by the directions of
the neutron spin and the emitted electron) in the decay distribution of the
neutron, $n\to p+e^{-}+{\overline{\nu}}_{e}$. The best determination of the
ratio of the axial to the vector charge, $g_{A}/g_{V}$, comes from using (i)
polarized ultracold neutrons (UCN) by the UCNA collaboration, $1.2772(20)$ [1,
2], and (ii) cold neutron beam by PERKEO III, $1.27641(45)(33)$ [3, 4]. Note
that, in the SM, $g_{V}=1$ up to second order corrections in isospin breaking
[5, 6] as a result of the conservation of the vector current.

The axial charge enters in many analyses of nucleon structure, of the Standard
Model (SM), and in probes of beyond-the-SM (BSM) physics [7, 8]. For example,
it enters in the relation between the Cabibbo-Kobayashi-Maskawa (CKM) matrix
element $V_{ud}$ and the neutron lifetime, $\tau_{n}$. High precision
extraction of $V_{ud}$, knowing $\tau_{n}$ and $g_{A}$, is important for the
test of the unitarity of the first row of the CKM matrix [9, 10, 11]. It is
needed in the analysis of neutrinoless double-beta decay [12] and in the rate
of proton-proton fusion [13], the first step in the thermonuclear reaction
chains that power low-mass hydrogen-burning stars like the sun.

The axial-vector form factor (AVFF) gives the dependence of this coupling on
the momentum squared transferred by the weak current to the nucleon. It is an
input in the theoretical calculation of the neutrino-nuclei scattering cross-
section needed for the analysis of neutrino oscillation experiments [14, 15,
16]. The cleanest experimental measurement would be from scattering neutrinos
off liquid hydrogen targets, however, these are not being carried out due to
safety concerns. Extractions from ongoing neutrino scattering experiments
(T2K, NOvA, MINERvA, MicroBooNE, SBN) have uncertainty due to those in the
cross-section and the incoming flux, and from the lack of precise
reconstruction of the final state of the struck nucleus. Conversely,
uncertainty in the AVFF feeds into the uncertainty in the incoming neutrino
energy that is needed to determine the oscillation parameters.

The MINER$\nu$A experiment [17], using a separation based on kinematics, has
recently extracted the axial-vector form factor of the nucleon from the
charged current elastic scattering process
${\overline{\nu}}_{\mu}H\to\mu^{+}n$ in which the free proton in the hydrogen
(H) (part of the hydrocarbon in the target) gets converted into a neutron.
This opens the door to direct measurements of the nucleon axial-vector form
factor without the need for extraction from scattering off nuclei, whose
analysis involves nuclear corrections that have unresolved systematics. On the
theoretical front, lattice QCD provides the best method for first principal
non-perturbative predictions with control over all sources of uncertainty [15,
14].

A recent comparison [18] of results for the AVFF from lattice QCD [19], the
MINER$\nu$A experiment [17], and the phenomenological extraction from
neutrino-deuterium data [20] showed that in the near term the best prospects
for determining the AVFF will be a combination of lattice QCD calculations and
MINER$\nu$A-like experiments. Lattice QCD will provide the best estimates for
$Q^{2}\lesssim 0.5$ GeV2, and be competitive with MINER$\nu$A for $0.5\lesssim
Q^{2}\lesssim 1$ GeV2. For $Q^{2}\gtrsim 1$ GeV2, new ideas are needed for
robust predictions using lattice QCD.

The goal of theory efforts in support of neutrino oscillation experiments is
robust calculations of the cross-section for targets, such as ${}^{12}C$,
${}^{16}O$, and ${}^{40}Ar$, being used in experiments. This involves a four
step process: a precise determination of the AVFF, nuclear models of the
ground state of the nuclei from which the neutrino scatters, the intra-nucleus
evolution of the struck nucleon using many-body theory to include complex
nuclear effects up to $\approx 5$ GeV for the DUNE experiment, and the
evolution of the final state particles to the detectors. The overall program
requires complete implementation of these within Monte Carlo neutrino event
Generators [14, 15, 16] with full uncertainty quantification in each step. The
output of the generators is the essential input required by experimentalists
for determining neutrino oscillation parameters from current and future
experiments.

Here I review the status of lattice QCD calculations of the axial charge,
$g_{A}^{u-d}$ and the AVFF. In addition, note that the flavor diagonal axial
charges $g_{A}^{u,d,s,c,b}$ provide the contribution of each quark flavor to
the spin of the nucleon, whose calculation is computationally more expensive
due to the additional disconnected contributions. Current status of results
for these nucleon charges has been reviewed by the Flavor Lattice Averaging
Group (FLAG) in 2019 and 2021 [21, 22]). Including results post FLAG 2021 [23,
24, 25, 26, 27], the values from the various calculations with 2+1- and
2+1+1-flavors of sea quarks lie in the ranges $1.22\lesssim
g_{A}^{u-d}\lesssim 1.32$, $0.74\lesssim g_{A}^{u}\lesssim 0.89$,
$-0.48\lesssim g_{A}^{d}\lesssim-0.38$, and $-0.06\lesssim
g_{A}^{s}\lesssim-0.025$. There have been no substantial new results for
flavor diagonal charges since the FLAG reports, so I will not discuss them
further in this work.

Based on the results in Refs. [19, 23, 27, 26, 28, 29], I present the case
that lattice results for AVFF over the range $0<Q^{2}\leq 1$ GeV2 are also
available with $\lesssim 10\%$ uncertainty and agree with MINER$\nu$A results
to within a combined sigma as discussed in Ref. [18] but disagree with the
neutrino-deuterium results for $Q^{2}>0.3$ GeV2. At the same time, I also
highlight the need for much higher statistics and better control over excited
state contributions to nucleon correlators in lattice calculations for the
uncertainty to be reduced to the percent level.

The outline of this review is as follows. I will summarize the methodology and
steps in the calculation of the axial and pseudoscalar form factors in Sec. 2.
This includes a discussion of the nucleon 3-point correlation functions
calculated in Sec. 2.1, removing possible excited state contributions (ESC) in
Sec. 2.2, and how form factors are then obtained from them in Sec. 2.3. I then
review the operator constraint imposed on the three form factors, the axial,
$G_{A}(Q^{2})$, the induced pseudoscalar, $\widetilde{G}_{P}(Q^{2})$, and the
pesudoscalar $G_{P}(Q^{2})$ by the axial Ward-Takahashi (also referred to in
literature as the partially conserved axial current (PCAC)) identity in Sec.
2.4, and how it provides a data driven method for validating the enhanced
contributions of multihadron, $N\pi$, excited states. These enhanced excited
state contributions are due to the coupling of the axial and pseudoscalar
currents to a pion, i.e., the pion pole dominance hypothesis. Extrapolation of
the lattice results to the physical point defined by the continuum ($a=0$) and
infinite volume ($M_{\pi}L\to\infty$) limits at physical light quark masses in
the isospin symmetric limit, i.e., $m_{u}=m_{d}$ set using the neutral pion
mass ($M_{\pi^{0}}=135$ MeV) is discussed in Sec. 2.5. A consistency check on
the extraction of the axial charge is discussed in Sec. 2.6. I will then
review the results for the AVFF obtained by the various lattice collaborations
after extrapolation to the physical point in Sec. 3, and the comparison of
lattice QCD result, the recent MINER$\nu$A data, and the phenomenological
extraction from the old neutrino-deuterium scattering data along with my
perspective on future improvements in Sec. 4. The concluding remarks are given
in Sec. 5.

##  2 Calculation of the axial vector form factors using lattice QCD

Figure 1: Quark line diagrams for the gauge invariant time-ordered correlation
functions $C^{\text{2pt}}(\bm{p};\tau)$ and $C_{J}^{\rm 3pt}(\bm{q};t,\tau)$.
The gluon lines are shown only to indicate that all possible gluon exchanges
between quarks are included, i.e., it is a fully non-perturbative calculation.
The electromagnetic and axial form factors are calculated by inserting the
vector, $V_{\mu}$, and axial, $A_{\mu}$, currents, respectively, with momentum
q at times $t$ in between the nucleon source and sink separated by time
$\tau$. The construction of the source for the sequential propagator $P_{\rm
seq}$ is shown schematically in the right panel. The two lower quark lines
within the dotted region are tied at the sink by $\cal N$, leaving the spin
and color indices of the third spinor to serve as the source.

The quark line diagrams for the 2-point, $C^{\text{2pt}}$ and the 3-point
$C_{J}^{\rm 3pt}(\bm{q};t,\tau)$ (with the insertion of the axial, $A_{\mu}$,
(or vector $V_{\mu}$) and pseudoscalar, $P$, currents) correlators are shown
in Fig. 1. The methodology for calculating these is the same in all ongoing
calculations. For $C_{J}^{\rm 3pt}(\bm{q};t,\tau)$, two kinds of quark
propagators are calculated by solving, using Krylov solvers such as conjugate
gradient and accelerated using multigrid [30], the linear equation $DP=\eta$
where $D$ is the Dirac matrix on the lattice and $\eta$ is a source vector.
The first, $P$, is constructed using a delta function or smeared $\eta$ and
shown moving forward from the source location, i.e., from the left blob in the
right panel in Fig. 1. The second is a sequential propagator, $P_{\rm seq}$,
shown moving backwards (red line from the right blob representing the sink)
from $\eta=nucleon$ source with definite momentum ${\textbf{p}}_{f}$. This
nucleon source, shown schematically by the part of the diagram lying inside
the dotted area, is constructed by contracting together two original ($P$)
propagators. The insertion of the current with 3-momentum q between the source
and sink nucleons then reduces to that between the original propagator and the
sequential propagator, again shown schematically by the top line. By momentum
conservation, the source nucleon is then projected to momentum
${\textbf{p}}_{i}={\textbf{p}}_{f}-{\textbf{q}}$. The Euclidean 4-momentum
transfer squared is given by $Q^{2}={\textbf{q}}^{2}-(E_{N}-M_{N})^{2}$.

In current calculations (the standard method) the nucleon interpolating
operator, ${\cal N}$, used is

$\displaystyle{\cal N}(x)=$
$\displaystyle\epsilon^{abc}\left[q_{1}^{aT}(x)C\gamma_{5}\frac{1\pm\gamma_{4}}{2}q_{2}^{b}(x)\right]q_{1}^{c}(x)\,,$
(1)

where $C=\gamma_{4}\gamma_{2}$ and the optional factor $1\pm\gamma_{4}$
projects on to positive parity nucleon states propagating forward/backward in
time for zero momentum correlators. Developing a variational basis of
interpolating operators that includes all states making significant
contributions, including $N\pi$ states, i.e., the holy grail of taming ESC, is
still work under progress [31, 32].

A short description of the six steps in the calculation of the AVFF that are
common to all fermion discretization schemes and independent of the selection
of input simulation parameters is given next in Secs. 2.1–2.6.

###  2.1 Correlation functions $C^{{\rm 2pt}}$ and $C_{J}^{\rm
3pt}({\textbf{q}};t,\tau)$

Two kinds of smeared sources, $\eta$, have been used to generate the original
and sequential quark propagators in most lattice calculations: (i) Wuppertal
[33] and (ii) exponential source [25]. The quark propagators so obtained are
stitched together to construct the gauge invariant time-ordered correlation
functions $C^{\text{2pt}}(\bm{p};\tau)$ and $C_{J}^{\rm 3pt}(\bm{q};t,\tau)$
shown in Fig. 1, whose spectral decompositions are

$\displaystyle C^{\text{2pt}}(\bm{p};\tau)\equiv$
$\displaystyle\langle\Omega\lvert{\cal T}({\cal N}(\tau)\overline{\cal
N}(0))\rvert\Omega\rangle=\sum_{i=0}\lvert
A_{i}^{\prime}\rvert^{2}e^{-E_{i}\tau}\,,$ (2) $\displaystyle C_{J}^{\rm
3pt}(\bm{q};t,\tau)\equiv$ $\displaystyle\langle\Omega\lvert{\cal T}({\cal
N}(\tau)J_{\Gamma}(t)\overline{\cal
N}(0))\rvert\Omega\rangle=\sum_{i,j=0}{A_{i}^{\prime}}{A_{j}^{\ast}}\langle
j\lvert J_{\Gamma}\rvert i^{\prime}\rangle e^{-E_{i}t-M_{j}(\tau-t)}\,,$ (3)

where $J_{\Gamma}=A_{\mu}={\overline{\psi}}\gamma_{\mu}\gamma_{5}\psi$ or
$J_{\Gamma}=P={\overline{\psi}}\gamma_{5}\psi$ is the quark bilinear current
inserted at time $t$ with momentum $\bm{q}$, and $\lvert\Omega\rangle$ is the
vacuum state. In the $C_{J}^{\rm 3pt}(\bm{q};t,\tau)$, using the quantum
mechanical right to left time ordering, the nucleon in the final state
$\langle j\rvert$ is, by construction, projected to zero momentum, i.e.,
$p_{j}=(M,\bm{0})$, whereas the initial state $\lvert i^{\prime}\rangle$ is
projected onto definite momentum $p_{i}=(E,\bm{p})$ with $\bm{p}=-\bm{q}$ by
momentum conservation. The prime in $\lvert i^{\prime}\rangle$ indicates that
this state can have non-zero momentum. Consequently, the states on the two
sides of the inserted operator $J$ are different for all ${\textbf{q}}\neq 0$.
The goal is to extract the ground-state matrix elements (GSME), $\langle
0\lvert J\rvert 0^{\prime}\rangle$, from fits to Eq. (3).

A major challenge in the analysis of all nucleon correlators is the
exponential decay of the signal-to-noise ratio, i.e., as
$e^{-(M_{N}-1.5M_{\pi})\tau}$ with the source-sink separation $\tau$ [34, 35].
With ($O(10^{5})$ measurements), a good signal in
$C^{\text{2pt}}(\bm{p};\tau)$ and $C_{J}^{\rm 3pt}(\bm{q};t,\tau)$ extends to
$\approx 2$ and $\approx 1.5$ fm, respectively.

At these $\tau$, the residual contribution of many theoretically allowed
radial and multihadron excited states are observed to be significant. These
states arise because the standard nucleon interpolating operator ${\cal N}$,
defined in Eq. (1) and used to construct the correlation functions given in
Eqs. (2) and (3), couples to a nucleon and all its excitations with positive
parity including multihadron states, the lowest of which are $N({\bf
p})\pi(-{\bf p})$ with $|{\bf p}_{\rm lowest}|=2\pi/La$ and $N({\bf
0})\pi({\bf 0})\pi({\bf 0})$. Since it is not known, a priori, which excited
states contribute significantly to a given $C_{J}^{\rm 3pt}(\bm{q};t,\tau)$,
the first goal is to develop methods to identify these and remove their
contributions. Operationally this boils down to knowing/determining the
energies $E_{i}$ to put in fits to data using the theoretically rigorous (for
unitary actions) spectral decomposition given in Eq. (3). Note that the
$A_{i}$ are not needed as they come in combinations
${A_{i}^{\prime}}{A_{j}^{\ast}}\langle j\lvert J_{\Gamma}\rvert
i^{\prime}\rangle$, which are fit parameters and never used.

Figure 2: Data for the ratio $C_{J}^{\rm 3pt}(\bm{q},t,\tau)/\sqrt{C^{\rm
2pt}(\bm{q},t)C^{\rm 2pt}(0,\tau-t)}$ that, in the limits $(\tau-t)\to\infty$
and $t\to\infty$, should be independent of $\tau$ and $t$, i.e., lie on a
horizontal line in the center about $t=\tau/2$ with value that is the GSME.
Current data show large ESC, and the grey band is the estimate of GSME given
by the fit to Eq. (3). In each row, the data in each pair of panels are the
same but the fit on the left is without the $N\pi$ state and on the right is
with. The top row (panels 1 and 2) show the data and fit to $J=A_{4}$ with
${\bm{n}}=(0,0,1)$. These two panels illustrate (i) the improvement
($\chi^{2}/dof$) in the fit to $J=A_{4}$ data with the inclusion of the $N\pi$
state and (ii) a very large ESC indicated by the large slope slowly rotating
counterclockwise to the expected horizontal band. The right two panels show
the data and fits to $J=P$ with ${\bm{n}}=(1,0,0)$ that should give $G_{P}$.
These two panels illustrate that the difference in $G_{P}$ with and without
including the $N\pi$ state is about $50\%$ (enhanced ESC) and the
$\chi^{2}/dof$ is better with the $N\pi$ state. Panels in rows two and three
show the data and fit to $J=A_{3}$ with ${\bm{n}}=(1,0,0)$, $(2,0,0),\
(2,1,0),\ (3,1,0)$ that should give $G_{A}$ for $\tau\to\infty$. Each pair of
panels illustrates that the difference in $G_{A}$ with and without including
the $N\pi$ state is a few percent and the $\chi^{2}/dof$ of the two fits is
comparable. Note the change in the behavior: the ${\bm{n}}=(1,0,0)$ data
converge from below, while for (1,1,1) and higher momenta, the data are
rotating clockwise to the expected horizontal line. Also, the fits become less
robust with increasing $\bm{n}$. See Ref. [27] for details on these data and
the fits.

Figure 3: (Left) Results for the energy gaps, labeled $\Delta E_{1}^{A4}$ and
$\Delta M_{1}^{A4}$, for the first excited state extracted from fits to the
$C_{A_{4}}$ correlator. These mass gaps are compared with the first excited
state energy $\Delta E_{1}^{\rm 2pt}$ from four-state fits to the nucleon two-
point correlator. Note that the difference between them (black circles versus
blue triangles), and consequently the difference between the form factors
extracted increases as $M_{\pi}\to 135$ MeV and ${\bm{n^{2}}}\to 0$
(equivalently $Q^{2}\to 0$). (Right) The 1-loop correction to the 3-point
function in $\chi$PT.

Figure 4: The quark line diagrams illustrating the contribution of $N\pi$
states. (Left) the current $A_{\mu}$ annihilates the pion produced by the
source. (Right) The states involved in the transitions: the ground state
process is shown above the quark line diagram and those involving an excited
state on one side of the current insertion are shown below.

###  2.2 Extracting the ground state matrix elements: exposing and
incorporating $N\pi$ states

The most direct way to extract $\langle 0\lvert J\rvert 0^{\prime}\rangle$ is
to make fits to Eq. (3) keeping as many intermediate states as allowed by
data’s precision and demonstrate convergence. The problem is that even
unconstrained 2-state fits are numerically ill-posed. The next option is to
take the $E_{i}$ from $C^{\text{2pt}}(\bm{p};\tau)$ as ${\cal N}$ creates the
same set of states in $C^{\text{2pt}}$ and $C^{\text{3pt}}_{J}$ and input
these in fits to $C_{J}^{\rm 3pt}(\bm{q};t,\tau)$, either by doing
simultaneous fits or via priors within say a bootstrap process to correctly
propagate the errors. Of these, the ground state ${A_{0}^{\prime}}$,
${A_{0}}$, and $E_{0}$ are well-determined from fits to the 2-point function.
Similarly, one would expect $E_{1}$ can also be taken from $C^{\text{2pt}}$.
This was the strategy used until 2017 when it was shown in Ref. [36] that the
resulting form factors do not satisfy the constraint imposed on them by PCAC.
Deviations from PCAC due to discretization effects of about $\approx 5\%$ were
expected, however, almost a factor of two was found on the physical pion mass
ensembles.

The reason was provided by Bär [37, 38] using $\chi$PT: enhanced contributions
to ME from multihadron, $N\pi\,,\ldots$, excited states that have much smaller
mass gaps than of radial excitations, the lowest being $\approx 1230$ versus
N(1440). These states were not evident in fits to
$C^{\text{2pt}}(\bm{p};\tau)$ as they have small amplitudes. A different
approach to analysis that includes the $N\pi$ states was needed.

It turned out that 2-state fits to $C^{\rm 3pt}_{A_{4}}$ exposed these states
and provided a data-driven method [39]. These fits confirmed that the lowest
of the tower of $N({\bf p})\pi(-{\bf p})$ states makes a very significant
contribution. By itself, $C^{\rm 3pt}_{A_{4}}$ is dominated by excited states
and fits to it using the $E_{i}$ from $C^{\text{2pt}}$ gave very poor
$\chi^{2}/dof$ . Making fits leaving $E_{1}$ a free parameter dramatically
improved the $\chi^{2}/dof$ (compare the left two panels in the top row of
Fig. 2), and the resulting output values of $E_{1}$ on the ${\bf p}=0$ (${\bf
p}$) side of operator insertion were roughly consistent with $N({\bf
p}=1)\pi({\bf p}=-1)$ ($N(0)\pi({\bf p})$) as shown in Fig. 3 (left),
reproduced from Ref. [39]. An illustration of the current understanding of the
process giving GSME and of those involving the lowest excited states
contributing is shown in Fig. 4 (right). The caption of Fig. 2 points out some
of the features of ESC observed in current data and the efficacy of fits to
$C_{J}^{\rm 3pt}$ with and without including the lowest $N\pi$ state to remove
the ESC.

The fact that there is enhancement of ME in the axial channel has been
understood for over 60 years as the “pion pole dominance” (PPD) hypothesis. On
the lattice, the creation of a $N\pi$ state by $\cal N$ is suppressed by $V$,
the 3-d volume, compared to just the nucleon as each state has a normalization
factor of $1/V$ for a point (local) source $\cal N$. The axial current can,
however, couple to this pion, and because the pion is light, this coupling can
occur anywhere in the time slice at which the current is inserted with
momentum $\bm{q}$ (see Fig. 4 (left)). This gives a factor of $V$ enhancement,
approximately cancelling the normalization factor $1/V$ [37, 38]. This
enhanced contribution to the ME when the pion comes on-shell is an artifact
that has to be removed. Note that since energy is not conserved on the
lattice, both the neutron and the pion can come on shell, however, since
momentum is conserved, possible excited states must have the same total
momentum as the created neutron state. PPD tells us that the axial current
with momentum $\bm{q}$ can be viewed as the insertion of a pion with $\bm{q}$,
and this has a large coupling to the nucleon. These processes are illustrated
in Fig. 4.

Having identified large contributions from the $N({\bf p}=1)\pi({\bf p}=-1)$
state, certainly in the extraction of $\widetilde{G}_{P}$ and $G_{P}$, the
question is—do we need to include other multihadron and radial excited state
contributions if we want results with percent level precision? What about in
$G_{A}$? Note that, in addition to the enhanced contribution shown in Fig. 4,
$\chi$PT also indicates that the 1-loop contributions due to the diagram shown
in Fig. 3 (right) (again a $N\pi$ contribution) could be $O(5\%)$ in all the
five $C_{J}^{\rm 3pt}$. Thus the $N\pi$ state could also be significant for
extracting $G_{A}$ and the charge $g_{A}$ (the GSME $\langle 0\lvert
A_{3}\rvert 0^{\prime}\rangle$) from $C_{A_{3}}$ at the percent precision
desired. Based on these arguments, it is clear that one needs at least 3-state
fits to the five $C_{J}^{\rm 3pt}$ in Eq. (3)—the ground state, the $N\pi$
state and the third that effectively accounts for all other excited state
contributions.

In my evaluation, details of the fits made to remove ESC are the most
significant differences between the calculations performed by the different
collaborations. With the current methodology, higher statistics data are badly
needed to improve these fits, i.e., include more states in the fits, and
gradually reduce the dependence on exactly how the analyses are done.

A very important point to remember in such analyses is that fit parameters in
a truncated ansatz (the $A_{i}$ and $E_{i}$ in say a 3-state fit in our case)
try to incorporate the effects of all contributions. Thus the connection
between parameters coming out of fits to a truncated Eq. (3) and physical
states made in Figs. 3 and 4 is very approximate at best. Thus, when I write
“$N\pi$ state” I really only mean an $E_{1}$ close to that of a non-
interacting $N(1)\pi(-1)$ state and/or an $N(0)\pi(0)\pi(0)$ state that is
essentially degenerate on our lattices.

###  2.3 Extracting the form factors

Once the GSME, $\langle 0\lvert J\rvert 0^{\prime}\rangle$, have been
extracted, their Lorentz covariant decomposition into the axial $G_{A}$,
induced pseudoscalar $\widetilde{G}_{P}$, and pseudoscalar $G_{P}$ form
factors is

$\displaystyle\langle N(\bm{p}_{f})\lvert A_{\mu}(\bm{q})\rvert
N(\bm{p}_{i})\rangle$
$\displaystyle={\overline{u}}(\bm{p}_{f})\left[G_{A}(Q^{2})\gamma_{\mu}\gamma_{5}+q_{\mu}\gamma_{5}\frac{\widetilde{G}_{P}(Q^{2})}{2M}\right]u(\bm{p}_{i})\,,$
(4) $\displaystyle\langle N(\bm{p}_{f})\lvert P(\bm{q})\rvert
N(\bm{p}_{i})\rangle$
$\displaystyle={\overline{u}}(\bm{p}_{f})\left[G_{P}(Q^{2})\gamma_{5}\right]u(\bm{p}_{i})\,,$
(5)

where $u(\bm{p})$ is the nucleon spinor with momentum $\bm{p}$,
$q=p_{f}-p_{i}$ is the momentum transferred by the current,
$Q^{2}=-q^{2}=\bm{q}^{2}-(E({\bm{p}}_{f})-E({\bm{p}}_{i}))^{2}$ is the space-
like four momentum squared transferred. On the lattice, the discrete momenta
are ${\bm{p}}=2\pi{\bm{n}}/La=2\pi(n_{x},n_{y},n_{z})/La$ with
$n_{i}\in\pm\\{0\ldots L\\}$. The spinor normalization used is

$\displaystyle\sum_{s}u(\bm{p},s)\overline{u}(\bm{p},s)$
$\displaystyle=\frac{E(\bm{p})\gamma_{4}-i\gamma\cdot\bm{p}+M}{2E(\bm{p})}\,.$
(6)

It is important to note that the excited states have to be removed from the
correlation functions,$C$, which have the spectral decomposition given in Eq.
(3), and not from the form factors, i.e., after the decompositions. Eqs. (4)
and (5) are only valid for GSME, i.e., $\lvert N(\bm{p}_{i})\rangle$ and
$\langle N(\bm{p}_{f})\rvert$ are assumed to be ground states of the nucleon.
If there are residual ESC, then additional “transition” form factors have to
be included in the rhs of Eqs. (4) and (5).

Assuming GSME have been extracted, and choosing the nucleon spin projection to
be in the “3” direction, the explicit form of the decompositions in Eqs. (4)
and (5) become

$\displaystyle C_{A_{\\{1,2\\}}}(\bm{q})\to$ $\displaystyle
K^{-1}\left[-q_{\\{1,2\\}}q_{3}\frac{\widetilde{G}_{P}}{2M}\right]\,,$ (7)
$\displaystyle C_{A_{3}}(\bm{q})\to$ $\displaystyle
K^{-1}\left[-q_{3}^{2}\frac{\widetilde{G}_{P}}{2M}+(M+E)G_{A}\right]\,,$ (8)
$\displaystyle C_{A_{4}}(\bm{q})\to$ $\displaystyle
K^{-1}q_{3}\left[(M-E)\frac{\widetilde{G}_{P}}{2M}+G_{A}\right]\,,$ (9)
$\displaystyle C_{P}(\bm{q})\to$ $\displaystyle K^{-1}q_{3}G_{P}\,,$ (10)

where the kinematic factor $K\equiv\sqrt{2E(E+M)}$. In each case, data with
all equivalent momenta that have the same ${\bm{q}}^{2}$ are averaged to
improve the statistical signal. These correlation functions are complex
valued, and the signal, for the CP symmetric theory, is in
$\mathop{\mathrm{Im}}C_{A_{i}}$, $\mathop{\mathrm{Re}}C_{A_{4}}$, and
$\mathop{\mathrm{Re}}C_{P}$.

It is clear that $G_{P}$ is determined uniquely from $C_{P}$ (Eq. (10)), and
for certain momenta $\widetilde{G}_{P}$ from $C_{A_{\\{1,2\\}}}$ using Eq.
(7). The $C_{A_{3}}(\bm{q})$ and $C_{A_{4}}(\bm{q})$ give linear combinations
of $G_{A}$ and $\widetilde{G}_{P}$, and Eq. (8) gives only $G_{A}$ when
$q_{3}=0$.

###  2.4 Satisfying PCAC

The non-singlet PCAC relation between bare axial, $A_{\mu}(x)$, and
pseudoscalar, $P(x)$, currents is:

$\partial_{\mu}A_{\mu}=2\widehat{m}P\,,$ (11)

where the quark mass parameter $\widehat{m}\equiv Z_{m}m_{ud}Z_{P}Z_{A}^{-1}$
includes all the renormalization factors, and $m_{ud}=(m_{u}+m_{d})/2=m_{l}$
is the light quark mass in the isospin symmetric limit. Using the
decomposition in Eqs. (4) and (5) of GSME, the PCAC relation requires that the
three form factors $G_{A}$, $\widetilde{G}_{P}$, and $G_{P}$ must satisfy, up
to discretization errors, the relation

$2M_{N}G_{A}(Q^{2})-\frac{Q^{2}}{2M_{N}}\widetilde{G}_{P}(Q^{2})=2\widehat{m}G_{P}(Q^{2})\,,$
(12)

on each ensemble. All pre Ref. [36] calculations did not check this relation
and missed observing that the data showed large deviations. Calculations
subsequent to Ref. [39] that include the lowest mass gap state, $N({\bf
p}=1)\pi({\bf p}=-1)$, in the analysis, obtain form factors that satisfy PCAC
to within $\approx 10\%$ already at lattice spacing of $a\approx 0.9$ fm. (The
ETMC result is an exception as explained in Ref. [23]). An illustration of the
size of the deviation from unity of
$R_{1}+R_{2}\equiv\frac{2\widehat{m}G_{P}(Q^{2})}{M_{N}G_{A}(Q^{2})}+\frac{Q^{2}\widetilde{G}_{P}(Q^{2})}{4M_{N}^{2}G)_{A}(Q^{2})}$,
without and with the lowest $N\pi$ state included is shown in Fig. 5 taken
from Ref. [19].

Figure 5: (Left) Results for $R_{1}+R_{2}$ on 10 ensembles from fits to
$C_{J}^{\text{3pt}}$ without including the $N\pi$ state, i.e., the spectrum
taken from fits to $C^{\rm 2pt}$. (Right) Including the $N\pi$ state. For PCAC
to be satisfied, $R_{1}+R_{2}$ should be unity up to discretization errors.
The dotted lines show the 5% deviation band.

To summarize, satisfying the PCAC relation in Eq. (12) provides a strong and
necessary constraint on the extraction of the three axial form factors.
$\chi$PT analysis by Bär [37, 38] and data driven validation in Ref. [39, 29,
31] show that the lowest, $N({\bf p}=1)\pi({\bf p}=-1)$ and $N({\bf
p}=0)\pi({\bf p}=1)$, states makes a large contribution on the two sides,
respectively, and need to be included in the analysis. For percent level
precision, the next question is—what other states need to be included? Current
analyses include up to three states, where the third state, if its parameters
are left free, effectively tries to account for all residual ESC. Such fits
have been implemented in different ways. For example, in Ref. [29], the $N\pi$
state is hardwired and the third state is taken to be the lowest excited state
in fits to $C^{\rm 2pt}$. In Refs. [39, 27, 19], a simultaneous fit to all
five $J=A_{\mu}$ and $P$ correlators is made wherein the $A_{4}$ correlator
fixes $E_{1}$ to approximately the non-interacting energy of the $N\pi$ state.
Over time, with much higher statistics data, results from different
collaborations using different methods should converge as more more excited
states are included.

Figure 6: (Left) Comparison of the nucleon axial-vector form factor
$G_{A}\left(Q^{2}\right)$ as a function of $Q^{2}$, the momentum transfer
squared, obtained by the PNDME 23 [19] shown by the turquoise band; RQCD 19
[29] (light faun band); ETMC 21 [23] (faun band); NME 22 [27] (light brown
band); and Mainz 22 [26] (brown band). The $\nu$D band is the fit to the old
neutrino-deuterium data taken from Ref. [20].

###  2.5 Extrapolating lattice AVFF to the physical point for use in
phenomenology

The next step, once ESC have been removed and form factors have been extracted
from GSME on each ensemble, is to extrapolate these data to the physical point
and provide a parameterized form for $G_{A}$ and $\widetilde{G}_{P}$ that can
be used in phenomenology. The challenge is that the discrete set of
$Q^{2}_{i}$ values at which data are obtained are different on each ensemble.

One, simple to implement, way consists of the following three steps:

  1. 1.

Parameterize the $G_{A}(Q^{2}_{i})$ data on each ensemble. Depending on the
number of $Q^{2}_{i}$ values, it could be a suitably truncated $z$-expansion
or a Padé. The expected $1/Q^{4}$ asymptotic behavior can be built in by using
sum rules in the $z$-expansion [40] or through an $\\{n,n+2\\}$ Padé in
$Q^{2}/M_{N}^{2}$. The Mainz collaboration [26] combines the removal of ESC at
various values of $Q^{2}_{i}$ and the $Q^{2}$ parameterization on a given
ensemble to include correlations.

  2. 2.

Pick $n$ values of momenta, $Q_{k}^{2}$, over a range, say $0\leq Q^{2}\leq 1$
GeV2. Extrapolate the data at each of these values of $Q^{2}_{k}$ using a
simultaneous fit in $\\{m_{q},a,M_{\pi}L$} to the physical point. A typical
ansatz used for such chiral-continuum-finite-volume (CCFV) extrapolations is

$\displaystyle g(M_{\pi},a,M_{\pi}L)$
$\displaystyle=g_{0}+c_{1}a+c_{2}M_{\pi}^{2}+\frac{c_{3}M_{\pi}^{2}\exp{(-M_{\pi}L)}}{\sqrt{M_{\pi}L}}\,,$
(13)

where I have kept only the lowest order corrections in each of the
$\\{m_{q},a,M_{\pi}L\\}$ variables and assumed that discretization errors
start at $O(a)$.

  3. 3.

Having obtained the form factor in the continuum limit at the $n$ points,
$Q^{2}_{k}$, carry out the final parameterization, using again a truncated
$z$-expansion or a Padé.

This 3-step process can be done within a single bootstrap procedure to
propagate errors as has been done in Ref. [27, 19] to produce the NME and
PNDME results shown in Fig. 6. Or these steps can be combined, especially if
there are correlations between them. For example, to account for correlations
between the coefficients of the CCFV fits at different values of $Q^{2}_{k}$
in step 2.

The plots in Fig. 6 provide two comparisons. In panels on the left, physical
point results from the RQCD [29, 24], ETMC [23], NME [27], and Mainz [26]
collaborations are compared against those from PNDME [19]. On the right they
are overlaid and compared to the phenomenological extraction from the old
neutrino-deuterium bubble chamber data [20]. The PNDME, RQCD and NME data
mostly overlap, whereas the ETMC and Mainz data overlap and fall off slightly
slower for $Q^{2}\gtrsim 0.3$ GeV2. On the other hand, the neutrino-deuterium
($\nu$D) data [20] falls off much faster for $Q^{2}\gtrsim 0.2$ GeV2. Overall,
as shown in the right plot, the five lattice QCD estimates are consistent
within $1\sigma$ and lie about $2\sigma$ above the $\nu$D band for
$Q^{2}\gtrsim 0.3$ GeV2.

There also are results from CalLAT [41, 16], PACS [25, 42], and LHP+RBC+UKQCD
[43] collaborations, which have not been included in the comparison because
they have not been extrapolated to the physical point. The Fermilab
collaboration [32] has embarked on the much harder problem of calculating
transition matrix elements as well, e.g., $N\to\Delta$ or $N\to N\pi$.

From the analysis of the NME and PNDME data, my understanding is that the
differences in exactly how the ESC are handled by the various collaborations
and the consequent uncertainty in the final results should be considered work
in progress. The uncertainty from the differences in the overall procedure for
parameterization and CCFV extrapolation is, I believe, smaller because the
data do not show large dependence on any of the three parameters
$\\{m_{q},a,M_{\pi}L\\}$, especially for $a\lesssim 0.1$ fm and
$M_{\pi}L\gtrsim 4$, as illustrated in Fig. 7 [27, 19]. Hopefully, the next
generation calculations will shed light on, and possibly resolve, the various
differences.

Figure 7: The data for renormalized $G_{A}(Q^{2})$ on 13 HISQ ensembles show
small variation in $\\{a,M_{\pi},M_{\pi}L\\}$. The $a06m135$ data are
statistics limited. Figure reproduced from Ref. [19].

Other findings in Refs. [27, 19] are (i) the dipole ansatz
$G_{A}(Q^{2})=\frac{g_{A}}{(1+cQ^{2}/M_{N}^{2})^{2}}$ gives poor fits (very
low $p$ values) to data on many ensembles. Our conclusion, therefore, is that
the lattice data already show that the dipole ansatz does not have enough
parameters to capture the $Q^{2}$ behavior over the range $0\leq Q^{2}\leq 1$
GeV2. (ii) The PPD relation between $G_{A}$ and $\widetilde{G}_{P}$ works very
well.

Figure 8: (Left) Comparison of the parameterized nucleon axial-vector form
factor $G_{A}(Q^{2})$ versus $Q^{2}$ up to $2$ GeV2 obtained from (i) fit to
the deuterium bubble-chamber data [20] shown by blue solid lines with error
band; (ii) fit to recent MINERvA antineutrino-hydrogen data [17], shown by
black dashed lines and turquoise error band; and (iii) lattice QCD result
obtained by the PNDME Collaboration [19] shown by red solid lines without a
band. (Right) A comparison of antineutrino-nucleon charged current elastic
differential cross section using AFF from (i) lattice QCD by the PNDME
collaboration [19] (red bands) and (ii) the deuterium bubble-chamber data [20]
(black bands) with the MINERvA antineutrino-hydrogen data [17] (black
circles). These figures are taken from Ref. [18].

###  2.6 Consistency check in the extraction of the axial charge $g_{A}^{u-d}$

There are two ways in which one can extract the axial charge $g_{A}^{u-d}$.
The first is from the forward matrix element using $C_{A_{3}}$ in Eq. (8) with
${\bm{q}}=0$ and the second is by extrapolating the form factor
$G_{A}(Q^{2}\neq 0)$ to $Q^{2}=0$. I am considering them as separate because
the extraction from the forward matrix element is computationally clean:
$C_{A_{3}}({\bm{q}}=0)$ has the smallest errors and verification of the
symmetry of the data about $\tau/2$ is a good test. The errors grow with
${\bm{q}}$ as shown in Fig. 2. On the other hand, $G_{A}(Q^{2})$ is
constrained by being part of the PCAC relation, Eq. (12), that has to be
satisfied. The two results must agree after CCFV extrapolation. Based on the
data in Ref. [19], I conclude

  * •

The difference between $g_{A}^{u-d}$ extracted without and with including
$N\pi$ states is $O(\approx 6\%)$, i.e., $1.218(39)\to 1.294(48)$ on including
one (the lowest) $N\pi$ state in the analysis. Note that the errors in each
result are $O(\approx 3\%)$.

  * •

The difference between $g_{A}^{u-d}\equiv G_{A}(Q^{2}\to 0)$ extracted by
extrapolating $G_{A}(Q^{2})$ data obtained without and with including the
lowest $N\pi$ state is also $O(\approx 6\%)$, i.e., $1.213(39)\to 1.289(56)$.
Again, the errors in each are $O(\approx 4\%)$.

Thus, for each of the two cases, without and with $N\pi$ states, we get
consistent estimates for $g_{A}^{u-d}$ from the two methods, however results
including the $N\pi$ state are about $6\%$ larger. This difference is
consistent with the expected $\sim 5\%$ 1-loop correction to the charge in
$\chi$PT, however, it is roughly one combined $\sigma$. It therefore needs
validation. My pick for the final result is the analysis including the $N\pi$
state since it gives form factors that satisfy the PCAC relation. Higher
precision data are needed to further clarify the other significant ESCs and
how to include them.

Collaboration | $g_{A}^{u-d}$ |  $\langle r_{A}^{2}\rangle$ fm2 | $g_{P}^{\ast}$ | $g_{\pi NN}$
---|---|---|---|---
PNDME 23 | 1.292(53)(24) | 0.439(56)(34) | 9.03(47)(42) | 14.14(81)(85)
RQCD 19/23 | $1.284^{28}_{27}$ | 0.449(88) | 8.68(45) | 12.93(80)
ETMC 23 | 1.283(22) | 0.339(46)(6) | 8.99(39)(49) | 13.25(67)(69)
PACS 23 | 1.264(14)(1) | 0.316(67) |  |
Mainz 22 | 1.225(39)(25) | 0.370(63)(16) |  |
NME 21 | 1.32(6)(5) | 0.428(53)(30) | 7.9(7)(9) | 12.4(1.2)
CalLat 18 | 1.271(10)(7) |  |  |
PNDME 18 | 1.218(25)(30) |  |  |
Mainz 19 | 1.242(25)(${}^{0}_{-30})$ |  |  |
$\chi$QCD 18 | 1.254(16)(30) |  |  |
Table 1: Comparison of $g_{A}$, $\langle r_{A}^{2}\rangle$, $g_{P}^{\ast}$ and
$g_{\pi NN}$ from recent calculations: PNDME 23 [19], RQCD 19/23 [29, 24],
ETMC 23 [28], PACS 23 [42, 25], Mainz 22 [26], and NME 21 [27]. Lattice
results for all charges are usually quoted in the $\overline{\rm MS}$ scheme
at scale 2 GeV, however, note that the renormalization of $g_{A}^{u-d}$ is
trivial [21]. Earlier results for $g_{A}^{u-d}$ with 2+1+1-flavor simulations
by CalLat 18 [44] and PNDME 18 [45], and with 2+1-flavor simulations by Mainz
19 [46] and $\chi$QCD 18 [47] that entered in the averages compiled in the
FLAG report 2021 [21] are given below the dividing line.

##  3 Comparison of charges obtained by various lattice collaborations

The results for the axial charge, $g_{A}^{u-d}$, the charge radius squared,
$\langle r_{A}^{2}\rangle$, the induced pseudoscalar coupling $g_{P}^{\ast}$,
and the pion-nucleon coupling $g_{\pi NN}$ , extracted from $G_{A}$ and
$\widetilde{G}_{P}$ by various collaborations using the relations

$\displaystyle G_{A}(Q^{2})$ $\displaystyle=g_{A}(1-\frac{\langle
r_{A}^{2}\rangle}{6}Q^{2}+\cdots)\,,$ (14) $\displaystyle g_{P}^{\ast}$
$\displaystyle\equiv\frac{m_{\mu}}{2M_{N}}\widetilde{G}_{P}(Q^{\ast 2})\,,$
(15) $\displaystyle g_{\pi NN}$ $\displaystyle\equiv\;\lim_{Q^{2}\to-
M_{\pi}^{2}}\frac{M_{\pi}^{2}+Q^{2}}{4M_{N}F_{\pi}}{\widetilde{G}}_{P}(Q^{2})\,.$
(16)

are summarized in Table 1. Here $m_{\mu}$ is the muon mass and $Q^{\ast
2}=0.88m_{\mu}^{2}$ is the energy scale of muon capture, and $F_{\pi}=92.9$
MeV is the pion decay constant.

The results show about $10\%$ variation in $g_{A}^{u-d}$, $g_{P}^{\ast}$, and
$g_{\pi NN}$ and about $25\%$ in $\langle r_{A}^{2}\rangle$. Part of this is
likely due to different methodologies used in the analysis, in particular how
and if the lowest $N\pi$ state is included in the analysis. These results will
improve steadily over time with higher statistics data.

##  4 Comparison of differential cross-section using lattice AVFF with
MINER$\nu$A data

A comparison of the antineutrino-nucleon charged-current elastic cross
sections calculated using predictions of AVFF from lattice (PNDME 23 [19]) and
neutrino-deuterium analysis [20] with MINERvA measurement [17, 48] is
presented in Fig. 8 reproduced from Ref. [18]. A $\chi^{2}$ test was performed
to determine the significance of the differences between the three. No
significant difference was found between MINERvA-lattice QCD (PNDME) and
between MINERvA-deuterium results. A $\approx 2.5\sigma$ tension was, however,
found between the PNDME and the deuterium results. Based on data shown in Fig.
6, the deviation in the deuterium-Mainz and deuterium-ETMC will be even
larger. To assess the scope for future progress, three regions of $Q^{2}$ with
different prospects for the extraction of AVFF from lattice QCD and MINERvA-
like experiments were identified.

For $Q^{2}\lesssim 0.2\leavevmode\nobreak\ \mathrm{GeV}^{2}$, LQCD predictions
and fits to the deuterium bubble-chamber data are in good agreement. In this
region, the experimental errors in the measurement on hydrogen by MINERvA are
large, whereas the errors in the parameterization of the deuterium bubble
chamber data are smaller. The $\nu$D result has often been used as a
benchmark, however, note that there is unresolved uncertainty in the deuterium
data as discussed in Ref. [20]. Also, no new deuterium data are expected in
the near-term, so I do not comment on its future prospects. Lattice QCD data
are competitive and will improve steadily. This region will be well-
characterized by the axial charge, the axial charge radius, and well-
parameterized by a low-order $z$-expansion or a Padé.

For $0.2\leavevmode\nobreak\ \mathrm{GeV}^{2}\lesssim Q^{2}\lesssim
1\leavevmode\nobreak\ \mathrm{GeV}^{2}$, the AVFF from PNDME has the smallest
errors and the predicted differential cross section lies above the hydrogen
and $\nu$D values, i.e., the same ordering as for the AVFF shown in Fig. 8
(left). Future improvements in both the hydrogen data and lattice calculations
will provide robust cross-checks in this region.

The region $Q^{2}\lesssim 0.5\leavevmode\nobreak\ \mathrm{GeV}^{2}$ is where
lattice QCD data, even with current methodology, will improve rapidly as more
simulations are done closer to $M_{\pi}=135$ MeV, $a\to 0$ and on larger
volumes because, in these limits and for given statistics, the value of
$Q^{2}|_{\rm max}$ with a good signal (usually taken to be some fixed lattice
momentum $\textbf{n}^{2}$) decreases.

For $Q^{2}\gtrsim 1\leavevmode\nobreak\ \mathrm{GeV}^{2}$, current LQCD data
have larger statistical errors and systematic uncertainties—discretization and
residual excited state contributions. With the current methodology, the
lattice AVFF comes mostly from simulations with $M_{\pi}\gtrsim 300$ MeV on
$a<0.1$ fm ensembles [19]. If the dependence on $\\{a,M_{\pi}\\}$ is mild, as
has been observed so far, then these data are useful by themselves. With
higher statistics and improved actions, one can push the lattice momenta
${\bm{n}}^{2}|_{\rm max}$ higher and perhaps reliably reach $Q^{2}\sim 2\
\mathrm{GeV}^{2}$. Nevertheless, new methods are definitely needed to get data
at $Q^{2}\gtrsim 2\ \mathrm{GeV}^{2}$ from simulations with physical pion
masses, $M_{\pi}\approx 135$ MeV and $a<0.1$ fm. Similarly, improvements in
MINERvA and follow on experiments are needed to cover the full range of
$Q^{2}$ relevant for DUNE.

##  5 Concluding remarks

Extensive calculations of the AVFF are being carried out by at least the
following nine lattice QCD collaborations: PNDME [19], RQCD [29, 24], ETMC
[23], NME [27], Mainz [26], CalLAT [41, 16], PACS [25, 42], LHP+RBC+UKQCD
[43], and Fermilab [32]. As shown in Sections 3 and 2.5, we now have results
to within 10% precision. The major uncertainty comes from resolving and
removing excited-state contributions.

The good news is that the methodology for the calculation of the correlation
functions, $C^{\rm 2pt}$ and $C_{J}^{\rm 3pt}$, is robust. The bad news is
that the exponentially falling signal to noise ratio in them means that ESC
are large at source-sink separations possible in today’s calculations. Second,
it is also clear that multihadron, $N\pi$, excited states give large
contributions which must be removed. Unfortunately, it is not yet known how
many of these states need to be included in the analysis for percent level
precision. The operator constraint that the form factors satisfy the PCAC
relation in Eq. (12) provides a valuable check, so it must be carried out in
all calculations. The third challenge is getting data at large $Q^{2}$ because
the discretization and statistical errors grow with $Q^{2}$ on a given
ensemble. Also, the $Q^{2}|_{\rm max}$ (the largest lattice momenta
${2\pi{\bm{n}}}/{La}$ with a good signal to noise ratio for fixed statistics)
decreases as simulations are done closer to the physical point. Thus, to get
data for $Q^{2}\gtrsim 1$ GeV2 on physical pion mass ensembles will
need/benefit from new methods and very high statistics.

As first step towards percent level precision, my estimate is that a factor of
ten increase in statistics will reduce the statistical errors to a level that
will provide much more clarity in removing ESC. Similary new developments,
including variational methods [31] with multihadron states and momentum
smearing [49], will improve the calculations and extend the range of $Q^{2}$.
I anticipate continued improvements in both, statistics and methods, will
provide LQCD predictions of AVFF for nucleons in the range $Q^{2}\lesssim
2\leavevmode\nobreak\ \mathrm{GeV}^{2}$ (hopefully higher) with percent level
precision by about 2030, in concert with DUNE producing data.

## Acknowledgments

Many thanks to my collaborators, Yong-Chull Jang, Sungwoo Park, Oleksandr
Tomalak, Tanmoy Bhattacharya, Vincenzo Cirigliano, Huey-Wen Lin, Emanuele
Mereghetti, Santanu Mondal, and Boram Yoon for making the PNDME and NME
results happen. Also to all who presented relevant results at Lattice 2023 and
to André Walker-Laud for valuable comments. Our calculations used the Chroma
software suite [50] and resources at (i) the National Energy Research
Scientific Computing Center, a DOE Office of Science User Facility supported
by the Office of Science of the U.S. Department of Energy under Contract No.
DE-AC02-05CH11231; (ii) the Oak Ridge Leadership Computing Facility, which is
a DOE Office of Science User Facility supported under Contract DE-
AC05-00OR22725, through awards under the ALCC program project LGT107 and
INCITE award HEP133; (iii) the USQCD collaboration, which is funded by the
Office of Science of the U.S. Department of Energy; and (iv) Institutional
Computing at Los Alamos National Laboratory. R. Gupta was partly supported by
the U.S. Department of Energy, Office of Science, Office of High Energy
Physics under Contract No. DE-AC52-06NA25396 and by the LANL LDRD program.

## References

  * [1] UCNA Collaboration collaboration, M. Mendenhall et al., _Precision measurement of the neutron $\beta$-decay asymmetry_, _Phys.Rev._ C87 (2013) 032501 [1210.7048].
  * [2] UCNA collaboration, M. A. P. Brown et al., _New result for the neutron $\beta$-asymmetry parameter $A_{0}$ from UCNA_, _Phys. Rev._ C97 (2018) 035505 [1712.00884].
  * [3] B. Märkisch et al., _Measurement of the Weak Axial-Vector Coupling Constant in the Decay of Free Neutrons Using a Pulsed Cold Neutron Beam_ , _Phys. Rev. Lett._ 122 (2019) 242501 [1812.04666].
  * [4] D. Mund, B. Maerkisch, M. Deissenroth, J. Krempel, M. Schumann, H. Abele et al., _Determination of the Weak Axial Vector Coupling from a Measurement of the Beta-Asymmetry Parameter A in Neutron Beta Decay_ , _Phys. Rev. Lett._ 110 (2013) 172502 [1204.0013].
  * [5] M. Ademollo and R. Gatto, _Nonrenormalization Theorem for the Strangeness Violating Vector Currents_ , _Phys.Rev.Lett._ 13 (1964) 264.
  * [6] J. F. Donoghue and D. Wyler, _Isospin breaking and the precise determination of $V_{ud}$_, _Phys.Lett._ B241 (1990) 243.
  * [7] T. Bhattacharya, V. Cirigliano, S. D. Cohen, A. Filipuzzi, M. Gonzalez-Alonso, M. L. Graesser et al., _Probing Novel Scalar and Tensor Interactions from (Ultra)Cold Neutrons to the LHC_ , _Phys. Rev. D_ 85 (2012) 054512 [1110.6448].
  * [8] A. N. Ivanov, M. Pitschmann, N. I. Troitskaya and Y. A. Berdnikov, _Bound-state $\beta$$-$decay of the neutron re-examined_, _Phys. Rev. C_ 89 (2014) 055502 [1401.7809].
  * [9] A. Czarnecki, W. J. Marciano and A. Sirlin, _Neutron lifetime and axial coupling connection_ , _Phys. Rev. Lett._ 120 (2018) 202002.
  * [10] A. Czarnecki, W. J. Marciano and A. Sirlin, _Radiative Corrections to Neutron and Nuclear Beta Decays Revisited_ , _Phys. Rev._ D100 (2019) 073008 [1907.06737].
  * [11] A. Czarnecki, W. J. Marciano and A. Sirlin, _Pion beta decay and Cabibbo-Kobayashi-Maskawa unitarity_ , _Phys. Rev._ D101 (2020) 091301 [1911.04685].
  * [12] M. Horoi and A. Neacsu, _Shell model study of using an effective field theory for disentangling several contributions to neutrinoless double- $\beta$ decay_, _Phys. Rev. C_ 98 (2018) 035502 [1801.04496].
  * [13] B. W. Carroll and D. A. Ostlie, _An Introduction to Modern Astrophysics_. 2nd (international) ed., 2007.
  * [14] L. Alvarez Ruso et al., _Theoretical tools for neutrino scattering: interplay between lattice QCD, EFTs, nuclear physics, phenomenology, and neutrino event generators_ , 2203.09030.
  * [15] USQCD collaboration, A. S. Kronfeld, D. G. Richards, W. Detmold, R. Gupta, H.-W. Lin, K.-F. Liu et al., _Lattice QCD and Neutrino-Nucleus Scattering_ , _Eur. Phys. J._ A55 (2019) 196 [1904.09931].
  * [16] A. S. Meyer, A. Walker-Loud and C. Wilkinson, _Status of Lattice QCD Determination of Nucleon Form Factors and their Relevance for the Few-GeV Neutrino Program_ , _Ann. Rev. Nucl. Part. Sci._ 72 (2022) 205 [2201.01839].
  * [17] MINERvA collaboration, T. Cai et al., _Measurement of the axial vector form factor from antineutrino–proton scattering_ , _Nature_ 614 (2023) 48.
  * [18] O. Tomalak, R. Gupta and T. Bhattacharya, _Confronting the axial-vector form factor from lattice QCD with MINERvA antineutrino-proton data_ , _Phys. Rev. D_ 108 (2023) 074514 [2307.14920].
  * [19] Precision Neutron Decay Matrix Elements (PNDME) collaboration, Y.-C. Jang, R. Gupta, T. Bhattacharya, B. Yoon and H.-W. Lin, _Nucleon isovector axial form factors_ , _Phys. Rev. D_ 109 (2024) 014503 [2305.11330].
  * [20] A. S. Meyer, M. Betancourt, R. Gran and R. J. Hill, _Deuterium target data for precision neutrino-nucleus cross sections_ , _Phys. Rev. D_ 93 (2016) 113015 [1603.03048].
  * [21] Flavour Lattice Averaging Group (FLAG) collaboration, Y. Aoki et al., _FLAG Review 2021_ , _Eur. Phys. J. C_ 82 (2022) 869 [2111.09849].
  * [22] Flavour Lattice Averaging Group collaboration, S. Aoki et al., _FLAG Review 2019: Flavour Lattice Averaging Group (FLAG)_ , _Eur. Phys. J. C_ 80 (2020) 113 [1902.08191].
  * [23] C. Alexandrou, S. Bacchio, M. Constantinou, J. Finkenrath, R. Frezzotti, B. Kostrzewa et al., _Nucleon axial and pseudoscalar form factors using twisted-mass fermion ensembles at the physical point_ , 2309.05774.
  * [24] G. S. Bali, S. Collins, S. Heybrock, M. Löffler, R. Rödl, W. Söldner et al., _Octet baryon isovector charges from $N_{f}=2+1$ lattice QCD_, 2305.04717.
  * [25] PACS collaboration, R. Tsuji, N. Tsukamoto, Y. Aoki, K.-I. Ishikawa, Y. Kuramashi, S. Sasaki et al., _Nucleon isovector couplings in Nf=2+1 lattice QCD at the physical point_ , _Phys. Rev. D_ 106 (2022) 094505 [2207.11914].
  * [26] D. Djukanovic, G. von Hippel, J. Koponen, H. B. Meyer, K. Ottnad, T. Schulz et al., _Isovector axial form factor of the nucleon from lattice QCD_ , _Phys. Rev. D_ 106 (2022) 074503 [2207.03440].
  * [27] Nucleon Matrix Elements (NME) collaboration, S. Park, R. Gupta, B. Yoon, S. Mondal, T. Bhattacharya, Y.-C. Jang et al., _Precision nucleon charges and form factors using (2+1)-flavor lattice QCD_ , _Phys. Rev. D_ 105 (2022) 054505 [2103.05599].
  * [28] C. Alexandrou et al., _Nucleon axial and pseudoscalar form factors from lattice QCD at the physical point_ , _Phys. Rev. D_ 103 (2021) 034509 [2011.13342].
  * [29] RQCD collaboration, G. S. Bali, L. Barca, S. Collins, M. Gruber, M. Löffler, A. Schäfer et al., _Nucleon axial structure from lattice QCD_ , _JHEP_ 05 (2020) 126 [1911.13150].
  * [30] R. Babich, J. Brannick, R. C. Brower, M. A. Clark, T. A. Manteuffel, S. F. McCormick et al., _Adaptive multigrid algorithm for the lattice Wilson-Dirac operator_ , _Phys. Rev. Lett._ 105 (2010) 201602 [1005.3043].
  * [31] L. Barca, G. Bali and S. Collins, _Toward N to N $\pi$ matrix elements from lattice QCD_, _Phys. Rev. D_ 107 (2023) L051505 [2211.12278].
  * [32] A. V. Grebe and M. Wagman, _Nucleon-Pion Spectroscopy from Sparsened Correlators_ , 2312.00321.
  * [33] S. Gusken, U. Low, K. H. Mutter, R. Sommer, A. Patel and K. Schilling, _Nonsinglet Axial Vector Couplings of the Baryon Octet in Lattice QCD_ , _Phys. Lett. B_ 227 (1989) 266.
  * [34] G. Parisi, _The Strategy for Computing the Hadronic Mass Spectrum_ , _Phys. Rept._ 103 (1984) 203.
  * [35] G. P. Lepage, _The Analysis of Algorithms for Lattice Field Theory_ , in _Boulder ASI 1989:97-120_ , pp. 97–120, 1989, http://alice.cern.ch/format/showfull?sysnb=0117836.
  * [36] R. Gupta, Y.-C. Jang, H.-W. Lin, B. Yoon and T. Bhattacharya, _Axial Vector Form Factors of the Nucleon from Lattice QCD_ , _Phys. Rev._ D96 (2017) 114503 [1705.06834].
  * [37] O. Bär, _Nucleon-pion-state contamination in lattice calculations of the axial form factors of the nucleon_ , in _36th International Symposium on Lattice Field Theory (Lattice 2018) East Lansing, MI, United States, July 22-28, 2018_ , 2018, 1808.08738.
  * [38] O. Bär, _$N\pi$ -state contamination in lattice calculations of the nucleon axial form factors_, _Phys. Rev. D_ 99 (2019) 054506 [1812.09191].
  * [39] Y.-C. Jang, R. Gupta, B. Yoon and T. Bhattacharya, _Axial Vector Form Factors from Lattice QCD that Satisfy the PCAC Relation_ , _Phys. Rev. Lett._ 124 (2020) 072002 [1905.06470].
  * [40] G. Lee, J. R. Arrington and R. J. Hill, _Extraction of the proton radius from electron-proton scattering data_ , _Phys. Rev. D_ 92 (2015) 013013 [1505.01489].
  * [41] A. S. Meyer et al., _Nucleon Axial Form Factor from Domain Wall on HISQ_ , _PoS_ LATTICE2021 (2022) 081 [2111.06333].
  * [42] PACS collaboration, R. Tsuji, Y. Aoki, K.-I. Ishikawa, Y. Kuramashi, S. Sasaki, K. Sato et al., _Nucleon form factors in $N_{f}=2+1$ lattice QCD at the physical point : finite lattice spacing effect on the root-mean-square radii_, 2311.10345.
  * [43] LHP, RBC,, UKQCD collaboration, S. Ohta, _Nucleon isovector form factors from domain-wall lattice QCD at the physical mass_ , in _40th International Symposium on Lattice Field Theory_ , 11, 2023, 2311.05894.
  * [44] CalLat collaboration, C. C. Chang et al., _A per-cent-level determination of the nucleon axial coupling from quantum chromodynamics_ , _Nature_ (2018)  [1805.12130].
  * [45] R. Gupta, Y.-C. Jang, B. Yoon, H.-W. Lin, V. Cirigliano and T. Bhattacharya, _Isovector Charges of the Nucleon from 2+1+1-flavor Lattice QCD_ , _Phys. Rev._ D98 (2018) 034503 [1806.09006].
  * [46] T. Harris, G. von Hippel, P. Junnarkar, H. B. Meyer, K. Ottnad, J. Wilhelm et al., _Nucleon isovector charges and twist-2 matrix elements with $N_{f}=2+1$ dynamical Wilson quarks_, _Phys. Rev. D_ 100 (2019) 034513 [1905.01291].
  * [47] [$\chi$QCD 18] J. Liang, Y.-B. Yang, T. Draper, M. Gong and K.-F. Liu, _Quark spins and Anomalous Ward Identity_ , _Phys. Rev._ D98 (2018) 074505 [1806.08366].
  * [48] F. Irani, M. Goharipour, H. Hashamipour and K. Azizi, _New insight on the nucleon structure from recent MINERvA measurement of the antineutrino-proton scattering cross-section_ , 6, 2023.
  * [49] G. S. Bali, B. Lang, B. U. Musch and A. Schäfer, _Novel quark smearing for hadrons with high momenta in lattice QCD_ , _Phys. Rev. D_ 93 (2016) 094515 [1602.05525].
  * [50] SciDAC Collaboration, LHPC Collaboration, UKQCD Collaboration collaboration, R. G. Edwards and B. Joo, _The Chroma software system for lattice QCD_ , _Nucl.Phys.Proc.Suppl._ 140 (2005) 832 [hep-lat/0409003].


# Multiband Transport in Bilayer Graphene at High Carrier Densities

Dmitri K. Efetov  Patrick Maher  Simas Glinskis  Philip Kim  Department of
Physics, Columbia University New York, NY 10027

###### Abstract

We report a multiband transport study of bilayer graphene at high carrier
densities. Employing a poly(ethylene)oxide-CsClO4 solid polymer electrolyte
gate we demonstrate the filling of the high energy subbands in bilayer
graphene samples at carrier densities $|n|\geq 2.4\times 10^{13}$ cm-2. We
observe a sudden increase of resistance and the onset of a second family of
Shubnikov de Haas (SdH) oscillations as these high energy subbands are
populated. From simultaneous Hall and magnetoresistance measurements together
with SdH oscillations in the multiband conduction regime, we deduce the
carrier densities and mobilities for the higher energy bands separately and
find the mobilities to be at least a factor of two higher than those in the
low energy bands.

###### pacs:

73.63.b, 73.22.f, 73.23.b

Multiband transport is common for many complex metals where different types of
carriers on different pieces of the Fermi Surface (FS) carry electrical
currents. Conduction in this regime is controlled by the properties of the
individual subbands, each of which can have distinct mobilities, band masses,
and carrier densities. Other changes to the single-band conduction model
include inter-band scattering processes and mutual electrostatic screening of
carriers in different subbands, which alters the effective strength of the
Coulomb potential and hence adjusts the strength of electron-electron and
electron-charged impurity interactions.

To understand electronic conduction in this regime, it is desirable to study
the properties of the individual bands separately and compare these to the
properties in the multiband regime. This was achieved in 2-dimensional
electron gases (2DEGs) formed in GaAs quantum wells GaAs82 , where the
subbands can be continuously populated and depopulated by inducing parallel
magnetic fields. In these 2DEGs, an increased overall scattering rate due to
interband scattering was observed upon the single- to multiband transition,
GaAs90 ; GaAs97 , along with changes in the effective Coulomb potential which
led to the observation of new filling factors in the fractional quantum Hall
effect Shayegan10 .

Figure 1: (a) The tight-binding band structure of bilayer graphene for
interlayer asymmetries $\Delta=0$ eV (gray) and $\Delta=0.6$ eV (black). (b)
Schematic view of the double gated device, consisting of the SiO2/Si back gate
and the electrolytic top gate. Debye layers of Cs+ or ClO${}_{4}^{-}$ ions are
formed $d\sim 1$ nm above the bilayer and the gate electrode, respectively.
(c) Longitudinal resistivity and Hall resistance of the bilayer graphene
device at $T=2$ K as a function of back gate voltage $V_{bg}$ for 3 different
fixed electrolyte gate voltages $V_{eg}=$ -1.7, -0.4, and 1 V from left to
right, corresponding to predoping levels of $n_{H}=$ (-2.9, 0, 2.9)$\times
10^{13}$ cm-2. Inset shows an optical microscope image of a typical Hall Bar
device (the scale bar corresponds to 5 $\mu$m). Figure 2: Landau fan diagram
of the differential longitudinal resistivity
$\mathrm{d}\rho_{xx}/\mathrm{d}n_{H}$ for 3 different density ranges at $T=2$
K as a function of the Hall density $n_{H}$ and the magnetic field. (center)
The SdH oscillations in the LES converge at the CNP and flatten out at higher
$n_{H}$ due to decreasing LL separation. (left and right) For
$|n_{H}|>2.6\times 10^{13}$ cm-2 additional SdH oscillations appear,
originating at the resistivity spikes (red regions) that mark the onset of the
HES.

Bilayer graphene (BLG)Geim07 ; Geim08 ; McCann06 ; Rotenberg06 ; Novoselov07 ,
with its multiband structure and strong electrostatic tunability, offers a
unique model system to investigate multiple band transport phenomena. BLG’s
four-atom unit cell yields a band structure described by a pair of low energy
subbands (LESs) touching at the charge neutrality point (CNP) and a pair of
high energy subbands (HESs) whose onset is $\sim\pm$0.4 eV away from the CNP
(Fig.1(a)). Specifically, the tight binding model yields the energy dispersion
McCann06 :

$\epsilon_{1,2}^{\pm}(k)=\pm\sqrt{\frac{\gamma_{1}^{2}}{2}+\frac{\Delta^{2}}{4}+v_{F}^{2}{k}^{2}\pm\sqrt{\frac{\gamma_{1}^{4}}{4}+v_{F}^{2}{k}^{2}(\gamma_{1}^{2}+\Delta^{2})}},$
(1)

where the upper and lower index indicates the conduction (+) and valence (-);
and LES (1) and HES (2), $k$ is the wave vector measured from the Brillouin
zone corner, $v_{F}\approx$106 m/s is the Fermi velocity in single layer
graphene, $\gamma_{1}\approx 0.4$ eV is the interlayer binding energy, and
$\Delta$ is the interlayer potential asymmetry. Interestingly, since a
perpendicular electric field $E$ across the sample gives rise to an interlayer
potential difference $\Delta$, it opens up a gap in the spectrum of the LES
Rotenberg06 ; band gap ; Fai09  and is furthermore predicted to adjust the
onset energy of the HES. Whereas the LESs have been widely studied, the HESs,
with their expected onset density of $n\gtrsim 2.4\times 10^{13}$ cm-2
Rotenberg06 , have thus far not been accessed in transport experiments. This
can mainly be attributed to the carrier density limitations set by the
dielectric breakdown of the conventional SiO2/Si back gates, which do not
permit the tuning of carrier densities above $n\approx 0.7\times 10^{13}$ cm-2
($\epsilon_{F}\approx$ 0.2 eV).

In this letter, we report multiband transport in bilayer graphene. Using an
electrolytic gate, we were able to populate the HES of bilayer graphene,
allowing for both the LES and HES to be occupied simultaneously. The onset of
these subbands is marked by an abrupt increase of the sample resistivity, most
likely due to the opening of an interband scattering channel, along with the
appearance of a new family of Shubnikov-de Haas (SdH) oscillations associated
with the HES. A detailed analysis of the magneto- and Hall resistivities in
combination with the HES SdH oscillations in this regime enables us to
estimate the carrier mobilities in each subband separately, where we observe a
two-fold enhanced mobility of the HES carriers as compared to the LES carriers
at the same band densities.

Bilayer graphene devices were fabricated by mechanical exfoliation of Kish
graphite onto 300 nm thick SiO2 substrates, which are backed by degenerately
doped Si to form a back gate. The samples were etched into a Hall bar shape
with a typical channel size of $\sim$ 5 $\mu$m and then contacted with Cr/Au
(0.5/30 nm) electrodes through beam lithography (Fig. 1(c) inset). In order to
access the HES we utilized a recently developed solid polymer electrolyte
gating techniqueelectrolyte ; Iwasa08 ; Fai09 ; me10 , which was recently
shown to induce carrier densities beyond values of $n>10^{14}$ cm-2 me10  in
single layer graphene. The working principle of the solid polymer electrolyte
gate is shown in Fig.1(b). Cs+ and ClO${}_{4}^{-}$ ions are mobile in the
solid matrix formed by the polymer poly(ethylene)oxide (PEO). Upon applying a
gate voltage $V_{eg}$ to the electrolyte gate electrode, the ions form a thin
Debye layer a distance $d\sim 1$ nm away from the graphene surface. The
proximity of these layers to the graphene surface results in huge capacitances
per unit area $C_{eg}$, enabling extremely high carrier densities in the
samples. While CsClO4 has almost the same properties as the typically used
LiClO4 salt, we find a reduced sample degradation upon application of the
electrolyte on top of the sample, resulting in considerably higher sample
mobilities.

One major drawback of the electrolyte gate for low temperature studies is that
it cannot be tuned below $T<250$ K, where the ions start to freeze out in the
polymer and become immobile (though leaving the Debye layers on the bilayer
surface intact) Iwasa08 ; me10 . A detailed study of the density dependent
transport properties at low temperatures can therefore be quite challenging.
In order to overcome this issue, we employ the electrolyte gate just to
coarsely tune the density to high values ($|n|<$1014 cm-2) at $T=300$ K,
followed by an immediate cooldown to $T=2$ K. We then use the standard SiO2/Si
back gate to map out the detailed density dependence of the longitudinal sheet
resistivity $\rho_{xx}$ and the Hall resistance $R_{H}$, from which we extract
the total carrier density of the sample $n_{H}=B/eR_{H}$, with $B$ the
magnetic field and $e$ the electron charge. Here we find the back gate
capacitance to be $C_{bg}=141$ aF/$\mu$m2, almost unaltered by the presence of
the Debye layers on top of the sample.

In this experiment, we have measured $\rho_{xx}$ and $R_{H}$ of more than 3
BLG devices as a function of the back gate voltage $V_{bg}$ at various fixed
$V_{eg}$ corresponding to the wide density range of $n_{H}\sim\pm 8\times
10^{13}$ cm-2. Fig.1(c) shows $\rho_{xx}$ and $R_{H}$ for a representative
device for 3 selected cool-downs at $V_{eg}=$ -1.7, -0.4, 1 V from left to
right, corresponding to a pre-doping level of $n_{H}=$(-2.9, 0, 2.9)$\times
10^{13}$ cm-2. For low doping levels ($V_{eg}=$-0.4 V, Fig.1(c) middle) we
observe the expected Dirac Peak in $\rho_{xx}$ and the ambipolar transition of
$R_{H}$ as $V_{bg}$ sweeps through the CNP. Away from the CNP, $\rho_{xx}$ and
$R_{H}$ decrease as $|n_{H}|$ increases, as was observed before in BLG samples
Novoselov07 . For the strongly pre-doped gate sweeps however (Fig.1(c) left
and right), we observe a rather unexpected non-monotonic feature in the sample
resistivity. Instead of a monotonic decrease of $\rho_{xx}$ with increasing
$|n_{H}|$, it exhibits an abrupt increase by $\sim 10\%$ symmetrically at both
electron and hole sides at $n^{*}\sim|n_{H}|=2.6\times 10^{13}$ cm-2, a
carrier density which is consistent with theoretical expectations for the
onset density of the HES McCann06 ; Rotenberg06 . A similarly increasing
resistivity at the opening of a new subband was previously observed in 2D
electron systems formed in wide GaAs quantum wells GaAs82 ; GaAs90 ; GaAs97 ,
where it was attributed to a decreased overall scattering time $\tau$ due to
the opening of an additional inter-band scattering channel as the new subbands
are populated. Such an inter-band scattering mechanism between the LES and the
HES is also expected to give rise to a resistivity increase upon filling of
the HES in BLG samples, making it a likely candidate for the origin of the
observed resistivity increase. However, considering the strong differences
between the 2D electron gases in GaAs quantum wells and in BLG, including
vastly different densities of states, mobilities, and electron energies, more
theoretical work needs to be done to conclusively determine the origin of this
resistivity increase.

The electronic structure of the LES and HES can be further investigated by
studying the effect of the magnetic field $B$ on the longitudinal resistivity
$\rho_{xx}(B)$ in the various density ranges. Fig.2 shows the Landau fan
diagram of the differential sheet resistivity
$\mathrm{d}\rho_{xx}/\mathrm{d}n_{H}$ as a function of $B$ and $n_{H}$. Close
to the CNP (Fig.2 center) the SdH oscillations in the two LES are quite
pronounced, but with increasing density their amplitude quickly decays as the
energy separation of the Landau Levels (LL) decreases. Above the onset of the
HES (Fig.2 left and right), marked by the “spikes” of increased resistivity
(here the red regions) however, we observe another set of SdH oscillations
which form LL fans converging into the onset point of the HES.

In order to analyze the SdH oscillations, we now plot the $\rho_{xx}(B)$
traces for various fixed $n_{H}$ as a function of the inverse magnetic field
$B^{-1}$. Fig. 3(a) displays three exemplary traces at different $n_{H}$ above
the onset density of the HES. All traces show periodic oscillations in
$B^{-1}$ allowing us to obtain the SdH density,
$n_{SdH}=\frac{4e}{h}\Delta(B^{-1})$, assuming that each LL is both spin and
valley degenerate. Whereas for all $|n_{H}|<n^{*}$ we find that the obtained
$n_{SdH}\approx n_{H}$, indicating that the SdH oscillations are solely from a
single band (i.e., the LES), for $|n_{H}|>n^{*}$ the obtained $n_{SdH}$ values
are much smaller than the simultaneously measured $n_{H}$ values. This
behavior can be well explained by assuming that these SdH oscillations reflect
only the small fraction of charge carriers lying in the HES. For
$|n_{H}|>n^{*}$ we hence are able to extract the occupation densities of the
LES ($n_{LES}$) and HES ($n_{HES}$) from $n_{LES}=n_{H}-n_{SdH}$ and
$n_{HES}=n_{SdH}$. Fig. 3(b) shows the $|n_{LES}|$ and $|n_{HES}|$ in this
regime as a function of the total carrier density $|n_{H}|$. For each fixed
$V_{eg}$, the obtained $|n_{LES}|$ and $|n_{HES}|$ increase as $|n_{H}|$
increases (adjusted by $V_{bg}$), for both electrons and holes. Interestingly,
we notice that the $|n_{LES}(n_{H})|$ are slightly larger for larger
$|V_{eg}|$ while the trend is opposite for the HES, i.e. $|n_{HES}(n_{H})|$
are smaller for larger $|V_{eg}|$, even though their $n_{H}$ values are in
similar ranges. These general trends can be explained by an increase of the
interlayer potential difference $\Delta$ for increased values of $|V_{eg}|$,
which are predicted by the tight-binding model in Eq.1 to result in an
increase of the onset density (energy) of the HES.

Figure 3: (a) Exemplary traces of the longitudinal resistivity as a function
of inverse magnetic field at fixed values of $n_{H}$ beyond the onset of the
HES. (b) Carrier densities inferred from the SdH oscillations vs. the overall
Hall densities $n_{H}$, from 4 cooldowns at different set electrolyte gate
voltages $V_{eg}=$ -2 V (red), -1.7 V (orange), 1 V (purple), 1.4 V (blue).
(bottom) $|n_{HES}|$ vs. $n_{H}$, fitted with theoretical expectations for the
HES. (top) $|n_{LES}|$ vs. $n_{H}$, fitted with theoretical expectations for
the LES. Line traces correspond to theoretical fits for different values of
$\Delta=$ 0.31 eV (red), 0.17 eV (orange), 0.13 eV (purple), 0.26 eV (blue).

While a precise quantitative determination of the expected shift in the onset
density of the HES as a function of $V_{eg}$ and $V_{bg}$ requires a self-
consistent calculation of $\Delta(V_{eg},V_{bg})$ and would go beyond the
scope of this paper, we can still qualitatively test the above prediction.
This is possible since $\Delta$ is mostly controlled by $V_{eg}$, which has a
much stronger coupling to the BLG sample than the $V_{bg}$, thus allowing us
to approximately treat $\Delta$ as a constant for fixed $V_{eg}$. Since the
experimental traces displayed in Fig.3(b) correspond to different values of
$V_{eg}$ but the same ranges of $V_{bg}$, $\Delta$ is different for each trace
and can be extracted from the theoretical fits from Eq. 1, with $\Delta$ as
the only fitting parameter. Indeed for all 4 traces we find good agreement
with the theoretical fits; we clearly observe an enhanced onset density
(energy) for the traces with larger set potential differences across the
sample, which is in good qualitative agreement with theoretical predictions.

We now turn our attention to the transport properties of BLG in the limit of
$n_{H}>n^{*}$. The filling of these sub-bands creates a parallel transport
channel in addition to the one in the LES, thus defining the transport
properties in this regime by two types of carriers with distinct mobilities
$\mu_{1,2}$, effective masses $m^{*}_{1,2}$ and subband densities $n_{1,2}$
(here the index corresponds to the LES(1) and HES(2)) Fuhrer08 ; vanHouten88 .
In sharp contrast to a single band Drude model, where $\rho_{xx}(B)$ does not
depend on the $B$ field, in a two-carrier Drude theory it is expected to
become strongly modified, resulting in a pronounced $B$ field dependence
Ashcroft76 :

$\rho_{xx}(B)=\frac{n_{1}\mu_{1}+n_{2}\mu_{2}+(n_{1}\mu_{1}\mu_{2}^{2}+n_{2}\mu_{2}\mu_{1}^{2})B^{2}}{e((n_{1}\mu_{1}+n_{2}\mu_{2})^{2}+\mu_{1}^{2}\mu_{2}^{2}(n_{1}+n_{2})^{2}B^{2})},$
(2)

Fig.4(a) shows magnetoresistance traces for different fixed Hall densities
$n_{H}$. Close to the CNP, where only the LES are populated (Fig.4(a) black
trace), the $\rho_{xx}(B)$ traces are nearly flat as expected from the one-
fluid Drude theory. When the density is increased and the HES starts to fill
up, however, we observe a smooth transition to an approximately parabolic $B$
field dependence, resulting in a strong increase of $\rho_{xx}$ of up to
25$\%$ from 0 T to 8 T. Using the previously extracted carrier densities in
the two bands $n_{1,2}$ we can now fit the $\rho_{xx}(B)$ traces with the two-
carrier Drude model in Eq.2, with the mobilities of the two subbands
$\mu_{1,2}$ as the only fitting parameters. As shown in Fig.4(b) the
experimental finding are in excellent agreement with the theory, allowing us
to deduce the values of $\mu_{1,2}$ with good accuracy. Moreover, the ability
to extract the mobilities of the HES allows us now to characterize the HES in
more detail.

In Fig.4(c) we plot the extracted mobilities of the HES $\mu_{2}$ against the
carrier density in the HES $n_{2}$ and compare it to the mobilities $\mu_{1}$
of the LES at a similar range of subband densities in the LES $n_{1}$. We find
that the mobilities in the HES are at least a factor of two higher than those
in the LES. Considering that the effective carrier masses are similar for the
LES and the HES, this feature of the HES may be due to the enhanced screening
of charged impurity scatterers at higher carrier densities, effectively
reducing the scattering rate of the HES carriers on these scatterers. A more
detailed theoretical study is required, however, to undertake a quantitative
analysis of this problem.

Figure 4: (a) Longitudinal resistivity $\rho_{xx}(B)$ as a function of
magnetic field for $n_{H}=$ (0.84, 2.88, 3.00, 3.31, 3.83, 4.21, 6.96)$\times
10^{13}$ cm-2, from top to bottom. The $\rho_{xx}(B)$ traces undergo a smooth
transition from a nearly $B$ independent behavior when the HES is empty (below
$n_{H}<n^{*}\sim 2.6\times 10^{13}$ cm-2), to a strong, non-trivial $B$
dependence when the HES is occupied. (b) An exemplary $\rho_{xx}(B)$ trace at
$n_{H}=$6.96$\times 10^{13}$ cm-2 and $n_{SdH}=$1.42$\times 10^{13}$ cm-2 with
accompanying fit (dashed line) from Eq.2, using the mobilities in the LES
$\mu_{1}=541$ Vs/cm2 and the HES $\mu_{2}=2428$ Vs/cm2 as fitting parameters.
(c) The mobilities $\mu_{1,2}(n)$ as extracted from $\rho_{xx}(B)$ traces at
various fixed $n_{H}$ as a function of the density in the individual subbands.

In conclusion, using a polymer electrolyte gate we have achieved two-band
conduction in bilayer graphene. We have found that the filling of these bands
above a Hall density of $|n_{H}|>2.4\times 10^{13}$ cm-2 is marked by an
increase of the sample resistivity by $\sim 10\%$ along with the onset of SdH
oscillations. From simultaneous Hall and magnetoresistivity measurements, as
well as the analysis of the SdH oscillations in the two carrier conduction
regime, we have characterized the distinct carrier densities and mobilities of
the individual subbands, where we have found a strongly enhanced carrier
mobility in the HES of bilayer graphene.

The authors thank I.L. Aleiner, E. Hwang and K.F. Mak for helpful discussion.
This work is supported by the AFOSR MURI, FENA, and DARPA CERA. Sample
preparation was supported by the DOE (DE-FG02-05ER46215).

## References

  * (1) H. L. Stormer et. al., Sol. St. Com., 41, 707-709 (1982).
  * (2) G. R. Facer et al., Phys. Rev. B 56, 10036-10039 (1997).
  * (3) D. R. Leadley et al., Semi. Sci. Tech. 5, 1081-1087 (1990).
  * (4) J. Shabani et al., arXiv:1004.0979v2 (2010).
  * (5) A. K. Geim and K. S. Novoselov, Nat Mater 6, 183 (2007).
  * (6) A. K. Geim and P. Kim, Scientific American 298, 68 (2008).
  * (7) E. McCann et. al., Phys. Rev. B 74, 161403 (2006).
  * (8) T. Ohta, et.al., Science 313, 951 (2006).
  * (9) K. S. Novoselov, et. al., Nature Physics 2, 177 - 180 (2006) (2007).
  * (10) E. V. Castro et al., Phys. Rev. Lett. 99, 216802 (2007); J. B. Oostinga et al., Nature Mater. 7, 151 (2008); Y. Zhang et al., Nature (London) 459, 820 (2009);
  * (11) Matthew J. Panzer et.al., Advanced Materials 20, 3177 - 3180 (2008); A. Das et.al., Nature Nanotechnology 3, 210 - 215 (2008); J. Yan et. al., Phys. Rev. B 80, 241417 (2009);
  * (12) K. Ueno et.al., Nature Materials 7, 855 - 858 (2008);
  * (13) Kin Fai Mak et.al., Phys. Rev. Lett. 102, 256405 (2009);
  * (14) Dmitri K. Efetov and Philip Kim, Phys. Rev. Lett. 105, 256805 (2010).
  * (15) Sungjae Cho and Michael S. Fuhrer, Phys. Rev. B 77, 084102 (2008).
  * (16) H. van Houten et.al., Phys. Rev. B 37, 2756-2758 (1988).
  * (17) Neil W. Ashcroft and N. David Mermin, Thomson Learning Inc. (1976).