\section{Introduction}

Quantum Chromo-Dynamics (QCD) is the fundamental theory describing the interactions between quarks and gluons (partons), leading to their observed bound states (hadrons) and the strong nuclear force. At high spatial resolution (momentum scale), the QCD coupling constant becomes small (asymptotic freedom~\cite{Gross:1973id,Politzer:1973fx}), and quark and gluon interactions can be calculated perturbatively (pQCD). This leads to the emergence of these partons as effective degrees of freedom in the description of hard processes like deep inelastic scattering where the observed cross section can be described approximately as an incoherent sum of scattering cross sections on individual point-like and structureless partons. On the other hand, at low momenta and long distance scales, the interaction becomes strong and a perturbative treatment is no longer possible. Instead, physical processes can be best described in terms of effective hadronic degrees of freedom, {\it e.g.}, the excitation of resonant hadronic states. By varying the resolution of a probe from short to long distances, physical cross sections displays a transition from the partonic to the hadronic domains. It remains an important question whether there is a region where both pictures apply simultaneously, {\it i.e.}, whether a parton-based description can reproduce the data in the  kinematic region of hadronic resonances, at least on average.  This phenomenon is  known as Quark-Hadron Duality~\cite{Bloom:1970xb, Bloom:1971ye, DeRujula:1976baf, Melnitchouk:2005zr}.  While strong evidence for duality has been found, it is important to fully test the applicability of this concept in the case where spin degrees of freedom are present, and for different hadronic systems. If quark-hadron duality can be firmly established and its applicability quantitatively described, one can use  measurements of hadronic observables to improve contraints on the parton structure of these hadrons. For instance, measurements of nucleon structure functions that are sensitive to high parton momentum fraction $x$ are very difficult at high energies, which limits our knowledge of the very important behavior of  the underlying Parton Distribution Functions (PDFs) as $x \rightarrow 1$.  If the requirement of avoiding the region of nucleon resonances can be relaxed in a controlled manner,  data taken at lower energies could contribute invaluable information on this asymptotic behavior. 

In the present paper, we present new results on tests of duality in proton and neutron spin structure functions. Following this introduction, we introduce the relevant formalism and theoretical concepts, describe the data set we analyzed as well as the phenomenologically extracted spin structure functions from the JAM QCD global analysis to which we compare these data,
and then present results and conclusions.

\section{Theoretical background}
 
 In this paper, we focus on quark-hadron duality in polarized inclusive electron scattering off polarized nucleon targets. In the single photon exchange approximation, an electron with four momentum $l$ scatters with final momentum $l'$ 
 from a nucleon with momentum $p$ by exchanging a space-like virtual photon with momentum $q=l-l'$ (see Fig. \ref{fig:fey1}). 
 
 \begin{figure}[htb]
\centering
\includegraphics[width=6cm]{feyn.pdf}
\caption{Feynman Diagram for inclusive electron scattering off a nucleon target. $W$ is the invariant mass of
the unobserved final state $X$. All other symbols are explained in the text.}
\label{fig:fey1}
\end{figure}

The invariant cross sections can then be written as \cite{Collins:2011zzd}
\begin{align}
    E'\frac{d\sigma}{d^3 l'} = \frac{2\alpha^2}{s Q^4}L_{\mu\nu} W^{\mu\nu}
\end{align}
where $Q^2=-q^2$ is the virtuality of the exchanged photon, $L_{\mu\nu}$ is the leptonic tensor and $W_{\mu\nu}$ is the hadronic tensor. The latter can be written as a linear combination of unpolarized structure functions $F_{2,L}$  and the polarized structure functions $g_{1,2}$. The polarized structure functions can be experimentally accessed 
by measuring cross sections differences of the form 
\begin{align}
    d\sigma^{\downarrow \Uparrow}
   -d\sigma^{\uparrow \Uparrow}
\end{align}
where ${\downarrow \Uparrow}$ and ${\uparrow \Uparrow}$ corresponds to anti-parallel and parallel beam and target spin configurations, respectively.

In the kinematics of moderate $x= Q^2/2P\cdot q$ and $Q^2$  much larger than hadronic mass scales, the $g_1$ structure function can be approximated in collinear factorization schematically  as 
\begin{align}
g_1(x,Q^2)=& \sum_i \int_x^1 \frac{d\xi}{\xi} \Delta f_{i/N}(\xi,Q^2) 
            \Delta H_i\left(\frac{x}{\xi},\alpha_S(Q^2)\right)\notag\\
            &+O\left(\frac{m}{Q}\right)\; .
\label{e.factorization}
\end{align}
Here the sum runs over all parton flavors $i$. The term $\Delta H_i$ is the target-independent short-distance partonic coefficient function calculable in pQCD in powers of the strong coupling $\alpha_S$ and is convoluted with  the spin-dependent 
Parton Distribution Function (PDF) $\Delta f$ in the variable $\xi$. The factorization theorem is valid up to corrections of the order $m/Q$ where $m$ is a generic hadronic mass scales. The $\xi$ variable is the light-cone momentum fraction of partons relative to the parent hadron,
{\it  i.e.}, $\xi=k^+/p^+$. At leading order in pQCD, the hard factor $\Delta H_i$ is proportional to $\delta(x-\xi)$; hence the structure function $g_1$ has a leading order sensitivity to PDFs at $\xi=x$. Beyond the leading order however, 
the physical structure function receives PDF contributions in the range $x<\xi<1$ due to the convolution in 
Eq.~\ref{e.factorization}.
The scale dependence on $Q^2$ in $\Delta f$ is governed by the DGLAP evolution equations stemming from the renormalization of parton densities and  are  given as 
\begin{align}
    \frac{d\Delta f_i}{d\ln \mu^2}(\xi,\mu^2) = \sum_j \int_{\xi}^1 \frac{dy}{y} 
        \Delta P_{ij}\left(\frac{\xi}{y},\alpha_S(\mu^2)\right) \Delta f_j(y,\mu^2)
\end{align}
where $\Delta P_{ij}$ are the Altarelli–Parisi space-like splitting functions. Finally we remark that the  structure function $g_2$ has no leading power contribution. 

Since the focus of our study is the behavior of $g_1$ in the large-$x$, moderate $Q^2$ regime, it is important to utilize 
a QCD global analysis frameworks that has a maximal kinematical overlap in $x$ to allow us to study duality with minimal extrapolation. In \cite{Sato:2016tuz}, the Jefferson Lab Angular Momentum Collaboration (JAM)  carried out a comprehensive analysis of the double spin asymmetries in DIS with an extend kinematic coverage in $x$ and $Q^2$. 
To cover lower $Q^2$ in this analysis, it was necessary to go beyond the leading power factorization and include additional power corrections (see last expression in Eq.~\ref{e.factorization}) such as higher twist effects as well as target mass corrections using QCD operator product expansion of DIS in the moment space. 
In the following, we utilize the inferred $g_1$ from the JAM global analysis that has a kinematic convergence up to $x\sim 0.7$ and use DGLAP backward evolution to access the resonance region at high-$x$ and lower $Q^2$. 

\begin{figure}[htb]
\centering
\includegraphics[width=8cm]{example.pdf}
    \caption{
    Schematic dependence of the measured structure function $F_2$ in inelastic electron scattering off the nucleon on the variable $\omega^\prime = W^2/Q^2 +1$, which is close to $1/x$ at large $Q^2$. Panels (a) through (d) are for increasing four-momentum transfer $Q^2$. As can be observed, the resonance excitations of the nucleon are most prominent at low $Q^2$, while at higher $Q^2$ the curve for $F_2$ approaches the scaling limit (dashed line), hence indicating a transition to quark-hadron duality in this observable. Reproduced from the paper by E. D. Bloom and F. J. Gilman~\cite{Bloom:1970xb}, with the permission of AIP Publishing.}
   
\label{fig:BloomGil}
\end{figure}

For moderate final hadronic state masses, $W<2-3$ GeV, 
the cross section typically exhibits multiple resonance
peaks that appear when the target is excited into other baryonic states before later decaying into final state products.
This is illustrated in Fig.~\ref{fig:BloomGil} for the $F_2$ structure function. 
This so-called resonance region can be best described in terms of hadronic degrees of freedom, where the cross section is expressed  in terms of transition strengths to the various nucleon resonances,  together with non-resonant hadron production contributions~\cite{HillerBlin:2019jgp}.

It is not a priori obvious how this resonant behavior is related to the underlying degrees of freedom of all hadrons, quarks and gluons, and their description in terms of PDFs, perhaps augmented by higher-twist terms in the OPE. This is addressed by the concept of Quark-Hadron Duality that was first introduced in a publication by Bloom and Gilman in 1970~\cite{Bloom:1970xb,Bloom:1971ye}.  They found that the $F_2$ structure function measured in the nucleon resonance region  approaches a smooth ``scaling curve'' as $Q^2$ increases, with the resonant troughs and peaks approximately averaging out to match an extrapolation of the deep inelastic  structure function at high W into the resonance region (see Fig.~\ref{fig:BloomGil}). 

In particular, Bloom and Gilman proposed that integrals over specific ranges in  $\omega^\prime = 1/x + M^2/Q^2$ (or just over $x$) of either the extrapolated DIS fits or the experimental data in the resonance region would give similar results. The case where the limits of integration cover only 100-200 MeV on either side of a single resonance peak is referred to as local duality, as opposed to global duality which covers the entire resonance region from threshold to $W = 2$ GeV, potentially
also including the elastic peak. In either case, the relation can be summarized as
\begin{align}
 \int_{x_1(W_1,Q^2)}^{x_2(W_2,Q^2)}dx\ F_2^{res}(x,Q^2)
=\int_{x_1}^{x_2}dx\ F_2^{DIS}(x,Q^2),
\end{align}
where $F_2^{res}$ is the structure function {\bf measured} in the resonance region, while $F_2^{DIS}$ is 
extrapolated from a QCD global analysis. 
Here,
\begin{align}
x(W,Q^2) = \frac{Q^2}{W^2-M^2+Q^2} .
\label{xofW}
\end{align}


Since the initial discovery by Bloom and Gilman in 1970, considerable progress has been made in the measurement of unpolarized structure functions at low to moderate $Q^2$ and $W$ and their interpretation in terms of quark-hadron duality, notably at the Thomas Jefferson National Accelerator Facility (also known as Jefferson Lab) \cite{Niculescu:2000tk, Niculescu:2000tj,CLAS:2003iiq, Niculescu:2015wka, Christy:2007ve, Tvaskis:2016uxm, JeffersonLabE00-115:2009jll, Malace:2009dg}. 

In addition to this, spin dependent structure functions  in the same kinematic region have also been studied. Experiments at SLAC in the late 70's provided the first resonance region measurements for polarized proton-electron scattering \cite{Baum:1980mh, Baum:1983ha}. These experiments hinted at the applicability of Bloom-Gilman Duality to proton spin structure functions. They were followed in the early 90's by further experiments at SLAC which expanded the $g_1$ and $g_2$ measurements to the neutron as well \cite{Anthony:1993uf} \cite{Abe:1998wq}. In the latter half of the 90's, both DESY (the HERMES Collaboration, \cite{Airapetian:2002rw}) and Jefferson Lab (Hall B, \cite{CLAS:2002knl} and Hall A, \cite{Meziani:2004ne}) contributed to the investigation of duality in spin structure functions with increased kinematics. The new century brought additional high-precision experiments at Jefferson Lab (Halls A \cite{Zheng:2004ce, Zheng:2003un,JeffersonLabE01-012:2008kgk,Parno:2014xzb, Flay:2016wie}, B \cite{CLAS:2006ozz, CLAS:2008xos, CLAS:2006hmk, CLAS:2015otq, Fersch:2017qrq, Fersch:2018ydg} and C \cite{Wesselmann:2006mw}). 
Studying quark-hadron duality in the spin sector is important, since polarization dependent observable can have both positive and negative sign, and hence offer a more stringent test of duality.
In the present paper, we are presenting a new comparison of the most comprehensive data set on spin structure functions in the transition region between hadronic and partonic degrees of freedom, from the EG1b experiment~\cite{Fersch:2017qrq,CLAS:2015otq}, to the recent JAM QCD global analysis~\cite{Sato:2016tuz} at high $x$. 


\section{Input data}

For a detailed study of duality, one needs a dense set of data that cover the entire resonance region (conventionally from $W = 1.072$ GeV to 2 GeV) in fine $W$ bins, for a large number of bins in $Q^2$.  The most comprehensive such data set was collected by the ``EG1b'' experiment carried out with CLAS \cite{CLAS:2003umf} at Jefferson Lab during 2000-2001 \cite{CLAS:2006ozz, CLAS:2008xos, CLAS:2006hmk, CLAS:2015otq,Fersch:2017qrq}. The experiment used the polarized electron beam from the  Continuous wave Electron Beam Accelerator Facility (CEBAF) at Jefferson Lab,  with beam energies of 1.6, 2.5, 4.2, and 5.7 GeV. Together with the large acceptance of CLAS, this set of beam energies yielded a large kinematic reach (with partially overlapping regions), covering nearly 2 orders of magnitude in $Q^2$ ($Q^2 = 0.06...5$) and $W$ from threshold to about 3 GeV.  A particular advantage of the wide acceptance of CLAS is that the data could be sorted into a pre-determined grid of $Q^2$ and $W$, with no need to interpolate between different data points. 
 
The polarized nucleon targets were provided in the form of irradiated frozen ammonia and deuterated ammonia for measurements of proton and deuteron asymmetries, respectively. The target was polarized through Dynamic Nuclear Polarization and reached a polarization along the beam
direction of approximately 75\% for the protons and 30\% for the deuterons \cite{Keith:2003ca}.

The measured double-spin asymmetries were converted into spin
structure functions $g_1(W,Q^2)$ using a phenomenological fit to
the world data on polarized and unpolarized structure functions. In the
case of the neutron structure function $g_1^n$, a folding prescription \cite{Kahn:2008nq} was used to relate the measured spin structure function of the deuteron to $g_1^n$ for each kinematic point. This yielded the first data set of un-integrated neutron spin structure functions in the resonance region. Details about the experiment, the data analysis and the complete data sets can be found in \cite{CLAS:2015otq,Fersch:2017qrq}.

Extrapolated pQCD predictions for $g_1^p$ and $g_1^n$, which are compared to the resonance region data in this paper, are taken from 
the JAM15 fits \cite{Sato:2016tuz} of the world data on inclusive spin observables, including the EG1b data {\em outside} the resonance region ({\it i.e.}, for $W > 2$ GeV). The JAM fits used a novel iterative Monte Carlo fitting method that utilizes data resampling techniques and cross-validation for a robust determination of the uncertainty band of the fitted PDFs as well as any observables predicted from the fit. A total of 2515 data points from 35 experiments and 4 facilities (CERN, SLAC, DESY, and JLab) were included in the fit. 


\section{Analysis}


\begin{table}[h!]
\centering
\begin{tabular}{|l| c| c|} 
 \hline
\# & Lower $W$ limit & Upper $W$ limit\\
 \hline
1 & 0.939 & 2\\ 
2 & 1.08 & 2\\ 
3 & 0.939 & 1.38\\ 
4 & 1.08 & 1.38\\ 
5 & 1.38 & 1.58\\
6 & 1.58 & 1.82\\
7 & 1.82 & 2\\
 \hline
\end{tabular}
\caption{$W$ Ranges, in GeV}
\label{table:W}
\end{table}

In this paper, we investigate two different but related tests of duality:
a direct comparison between truncated integrals over measured spin structure functions and extrapolated pQCD fits,
each covering a specific range in final-state mass $W$, and 
a study of the approach to scaling for $g_1$ averaged over a set of narrow ranges in $x$. 

For the first test, we select
seven different ranges of $W$ as shown in table \ref{table:W}. 
The first two of these cover the entire ``canonical'' resonance region,
$W < 2$ GeV, either including (1) or excluding (2) the elastic contribution,
to test ``global'' duality.
The remaining 5 ranges cover specific prominent resonance peaks visible
in inclusive unpolarized cross section data (see Fig.~\ref{fig:BloomGil}):
the Delta resonance ($\Delta(1232)3/2^+$)  (again either combined with the elastic peak (3)
or without it (4)), the region of the N$(1440)1/2^+$, N$(1520)3/2^-$, and 
N$(1535)1/2^-$
resonances (5), the region of the N$(1680)5/2^+$ and nearby resonances (6), and the remaining region
up to $W = 2$ GeV (7) which does not exhibit a strong peak in the inclusive spin-
averaged cross section but is known to contain several Delta-resonances that 
tend to have negative (virtual) photon asymmetries. We test whether local
duality holds in each of these individual resonance regions.

\begin{table}[h!]	
\centering
\begin{tabular}{|c|c|c|}						
\hline						
Lower	$Q^2$	&	Upper	$Q^2$	&	Central $Q^2$\\
\hline		
0.92		&	1.10		&	1.00\\
1.10		&	1.31		&	1.20\\
1.31		&	1.56		&	1.43\\
1.56		&	1.87		&	1.71\\
1.87		&	2.23		&	2.04\\
2.23		&	2.66		&	2.43\\
2.66		&	3.17		&	2.91\\
3.17		&	3.79		&	3.47\\
3.79		&	4.52		&	4.14\\
4.52		&	5.40		&	4.94\\
5.40		&	6.45		&	5.90\\
\hline
\end{tabular}	
\caption{Experimental $Q^2$ Ranges, in GeV$^2$}
\label{table:Q}
\end{table}		

For each of these $W$ ranges, our analysis process is the same. Experimental data are first sorted into bins of $Q^2$ with limits shown in Table \ref{table:Q}.
The extrapolated ``pseudo-data'' points provided by the JAM collaboration
have been calculated from their PDF fits at the central $Q^2$ values of each bin.

\begin{figure}[htb!]
\centering
\includegraphics[width=0.4\textwidth]{3Dbins.pdf}
\caption{Representation of the experimental data set used in this analysis. The measured data points are binned
in bins in $Q^2$, as indicated by the shaded area for the example of the bin 1.87 GeV$^2 < Q^2 <$2.23 GeV$^2$;
see also Table~\ref{table:Q}. The truncated integrals are then formed over specific regions in $W$ as spelled out
in Table~\ref{table:W}}.
\label{fig:diagram}
\end{figure}

As a second step, both the measured data and the JAM ``pseudo-data'' for each $Q^2$ bin are
filtered to only include points within one of the given seven $W$ regions. 
This is done by
mapping the edges of each of our 7 regions to the corresponding 
values for $x$, following Eq.~\ref{xofW}.
Both can then be integrated over the corresponding $x$-ranges to yield
the truncated first moments of $g_1$,
\[ \bar{\Gamma}_1(\Delta W, Q^2) = \int_{x_1(W_1,Q^2)}^{x_2(W_2,Q^2)}dx\ g_1(x,Q^2), \]
for each of our 7 regions, as a function
of $Q^2$ bin.

 For the $W$ regions (1) and (3), the contribution from the elastic peak ($W=M$) was added
 ``by hand'' to the measured integrals to extend the
 truncated integrals up to $x_2 =1$, since the EG1b data do not reach
 below a value of $W = 1.072$ GeV.
 This elastic contributions comes in the form
\[g_1^{el} = \frac{1}{2}\frac{G_EG_M+\tau G^2_M}{1+\tau}\delta (x-1) ,\]
where $G_M = F_1 + F_2$ and $G_E = F_1 - \frac{Q^2}{4M^2}F_2$ are the magnetic and electric Sachs form factors \cite{Walker:1993vj}. The corresponding 
integrals for the JAM ``pseudo-data'' are simply integrated up to $x = 1$
but don't contain an elastic contribution (since they are based on DIS pQCD fits).

For the experimental data, the statistical and experimental errors are added in quadrature into the integration and displayed with corresponding error bars. 
For the JAM fits, we add or subtract one ``theoretical
standard deviation'' to
the central values, and the whole range between the corresponding integrals is displayed as band. Results for protons and neutrons are shown together on the same plot. 

For our second investigation, we define a sequence of bins in $x$, each
with a width of $\Delta x = 0.05$.
The measured $g_1$ points within each of these $x$-bins are averaged over each of the $Q^2$ bins and the averages are plotted {\it vs.} the nominal  $Q^2$ values. Again,
the JAM ``pseudo-data'' are treated in the same way and shown as bands
together with the data. 



\section{Results}

\begin{figure}[htb!]
\centering
\includegraphics[width=0.4\textwidth]{combined1.pdf}
\caption{[Color Online] Test of global duality. We show integrals
$\bar{\Gamma}_1(\Delta W, Q^2)$ of the spin structure
function $g_1$
over the entire resonance region, as a function
of $Q^2$.
In the top panel, we
include the elastic
peak: $0.938$ GeV $ < W < 2$ GeV, while the bottom
panel shows only the inelastic part,
$1.072$ GeV $ < W < 2$ GeV.
The top (red) bands 
and data points (circles) are for the proton, and the bottom (blue) bands and data points (triangles) are for the neutron. The data points are shown with statistical
and systematic uncertainties added in quadrature (error bars). The solid bands show the full prediction from the
extrapolated JAM fit, including target mass and higher
twist contributions. The striped band (proton) and the cross-hatched band (neutron) show the results including only the leading twist contribution.}
\label{fig:C1+C2}
\end{figure}

In this Section, we  present the results of
our two tests of duality. We begin by comparing the
integral of $g_1$ over the entire resonance region
(test of ``global duality''), see Fig.~\ref{fig:C1+C2}.
 In the top panel of Fig.~\ref{fig:C1+C2}, we integrate
the data over the entire region $0.938$ GeV $ < W < 2$ GeV,
including the elastic peak contribution. Correspondingly,
the JAM PDF fit is integrated up to $x=1$. 
For the proton
(upper bands and data points), a clear (and non-trivial)
agreement between data and PDF prediction is observed, starting around
$Q^2 = 1.4$ GeV$^2$. Initially, the data agree best with the
JAM prediction that includes higher twist contributions and target
mass corrections (solid band), but at higher $Q^2$ they
tend to be a bit lower and closer to the striped
band (which includes only leading twist).
For the neutron (lower bands and data points), the predictions from
the PDF fits as well as the data are mostly consistent
with zero, and there is little difference between the
extrapolation including higher twist and target mass corrections (solid band) and the one including only leading twist 
(cross-hatched band). There may be a slight tendency for
the data to fall below the PDF bands at high $Q^2$,
which would agree with the observation that the $d-$quark
polarization appears to remain negative up to the highest
$x$ measured~\cite{Zheng:2004ce}; however, the experimental uncertainties
are too large to make a definite statement.

\begin{figure}[htb!]
\centering
\includegraphics[width=0.4\textwidth]{combined2.pdf}
\caption{[Color Online] Test of local duality in
the region of the $\Delta$ resonance. Again,
in the top panel we
include the elastic
peak: $0.938$ GeV $ < W < 1.38$ GeV, while the bottom
panel shows only the inelastic part 
$1.072$ GeV $ < W < 1.38$ GeV.
All symbols are the same as in  Fig.~\ref{fig:C1+C2}}
\label{fig:C3+C4}
\end{figure}


We repeat our test of ``global duality'' by also comparing
the integrals excluding the elastic contribution (and, 
correspondingly, the region $x \rightarrow 1$), see
lower panel of
Fig.~\ref{fig:C1+C2}. While the prediction based on the
extrapolated PDF fit changes only slightly (as expected,
since $g_1 \rightarrow 0$ as $x \rightarrow 1$), we
see that the data for the proton, excluding the elastic contribution, fall
short of that expectation up to rather high $Q^2 > 3.5$ 
GeV$^2$. (The agreement for the neutron is only slightly worse).
This finding indicates that the elastic contribution must
be included to get a fairly rapid convergence (low threshold)
for global duality. This is due to the need to compensate
for the rather striking breakdown of ``local duality'' in the region of the 
Delta-resonance.

This breakdown can be seen in Fig.~\ref{fig:C3+C4}. Here, we integrate only over the peak
of the lowest-lying Delta resonance, either with the
addition (top panel of Fig.~\ref{fig:C3+C4}) of the elastic peak, or
without it (bottom panel of Fig.~\ref{fig:C3+C4}). In the latter case, both
the proton and the neutron data are either negative or
close to zero, while the PDF extrapolation for both
is positive (significantly so for the proton). Again,
agreement between proton data and PDF extrapolation only
begins around $Q^2 > 3$ GeV$^2$. This is due to the well-known
fact that the excitation of the Delta resonance is dominated
by a $M1$ transition, for which the final state helicity
$3/2$ has a stronger coupling than the final state helicity
$1/2$, leading to a negative (virtual)
photon asymmetry $A_1$ and, in consequence, a negative
value for $g_1$. Interestingly, this stark deviation
from the PDF extrapolation is more than compensated by adding
in the elastic peak (top panel of Fig.~\ref{fig:C3+C4}) , which has $A_1 = 1$ by definition. 
In that case, the proton data overshoot the pQCD prediction, while the neutron data are largely
in agreement.
\begin{figure}[htb!]
\centering
\includegraphics[width=0.4\textwidth]{C5.pdf}
\caption{[Color Online] Local duality test in
the region of the $N(1440)1/2^+$, $N(1520)3/2^-$, and 
$N(1535)1/2^-$ resonances. All symbols as in Fig.~\ref{fig:C1+C2}.}
\label{fig:P5}
\end{figure}

\begin{figure}[htb!]
\centering
\includegraphics[width=0.4\textwidth]{C6.pdf}
\caption{[Color Online] Local duality test in
the region of the $N(1680)5/2^+$ resonance. All symbols as in Fig.~\ref{fig:C1+C2}.}
\label{fig:P6}
\end{figure}


\begin{figure}[htb!]
\centering
\includegraphics[width=0.4\textwidth]{C7.pdf}
\caption{[Color Online] Local duality test for
the remainder of the resonance region, 1.82 GeV $< W <$ 2 GeV.
All symbols as in Fig.~\ref{fig:C1+C2}.}
\label{fig:P7}
\end{figure}

The next two proton resonance regions (Figs.~\ref{fig:P5}
and \ref{fig:P6}) show remarkably good agreement between
the data and the extrapolated PDF bands
(in particular the extrapolations including
higher twist), indicating that
``local duality'' works well for these resonances.
The remaining region, up to $W = 2$ GeV, shows again
a deviation of the data which tend to lie below the
PDF fits. Once again, this is consistent with the 
assumption that this region has a strong contribution from
various Delta-resonances, where the helicity-3/2 contribution
dominates at small $Q^2$. It is remarkable, coming back
to Fig.~\ref{fig:C1+C2}, how the negative deviations in the
lowest and highest $W$ regions (both
populated by Delta resonances) are compensated by the
inclusion of the elastic peak to get a rather rapid
approach to ``global duality''.

\begin{figure}[htb!]
\centering
\includegraphics[width=0.4\textwidth]{35-40.pdf}
\caption{[Color Online]
Approach towards the scaling limit for the structure function
$g_1(x,Q^2)$, averaged over the $x$-bin of $0.35 < x < 0.4$,
as a function of $Q^2$ bin.
Top: Data and JAM bands are shown multiplied with the average
$x = 0.375$; symbols are as in Figs.~\ref{fig:C1+C2}-\ref{fig:P7}.
The vertical dashed line indicates the limit $W=2$ GeV of the resonance
region, which lies to the left. 
Bottom: kinematic location of all data points from EG1b in
the $x$ {\it vs.} $Q^2$ plane. The grey band indicates the
interval in $x$ over which the data are averaged, and the red 
vertical lines indicate the nominal central values of each 
$Q^2$ bin.}
\label{fig:35}
\end{figure}

\begin{figure}[htb!]
\centering
\includegraphics[width=0.4\textwidth]{70-75.pdf}
\caption{[Color Online] Same as Fig.~\ref{fig:35} except for the
bin $0.7 < x < 0.75$.}
\label{fig:70}
\end{figure}

In our second analysis, we are integrating the EG1b data and JAM
PDF fits over fixed intervals in $x$ for each of our $Q^2$ bins,
to study the approach towards scaling for the structure function
$g_1(x,Q^2)$. The integrals are divided by the bin width $\Delta x = 0.05$
to obtain the average $g_1$ and then multiplied by the bin centroid in $x$ for better visibility - see Figs.~\ref{fig:35}-\ref{fig:70}. 
In contrast to the previous analysis, we include in these figures {\em all} data from
EG1b, from both the resonance and the DIS region, with the
boundary between the two indicated by the vertical dashed line at 
$Q^2=(W^2- m^2)/(1/x-1)$, with $W = 2$ GeV. The bottom half
of each figure shows all of the EG1b data points, where the points
included in each integral are those lying within the grey bands.

For the lower $x$ bin (Fig.~\ref{fig:35}), both the proton data
and the neutron data approach a smooth, largely flat behavior
at increasing $Q^2$, indicating the onset of scaling roughly 
around $Q^2 = 2$ GeV$^2$ which corresponds to $W > 2$ GeV. As in the earlier
duality studies, the neutron data have larger uncertainties and
are fully consistent with the expectations from the JAM PDF
results, while the convergence of the proton data appears to be a
bit slower, and they are closer to the bands including higher twist. In both cases there is no strong resonance
structure visible beyond $W =2$ GeV, while the proton data clearly
exhibit some non-trivial structure below $W =2$ GeV.
Since in this case, we are not integrating over any particular
structure in the resonance region but over a fixed $x$ bin, 
different resonances contributed differently to each data point,
spoiling the ``local duality'' agreement observed above for
some of the prominent resonances.

In contrast to the above, the data for the $x$-range 
$0.7 < x < 0.75$ are all inside the resonance region, since
the kinematic reach of EG1b (with a maximum beam energy of
5.7 GeV) was not sufficient to reach the limit $W = 2$ GeV
for this rather high $x$. Both the proton and the neutron
data show prominent resonance structure at low $Q^2$, while
approaching a more smooth behavior (and the JAM PDF predictions)
at moderately higher $Q^2$, above  $Q^2 \approx 3$ GeV$^2$. 
This corresponds to an average $W$ of only 1.4 GeV, still well
into the resonance region. Thus it appears as if the approach
to scaling may set on early at larger $x$ values, which would be very beneficial for the goal of extracting the behavior of
spin structure functions at large $x$, a topic of continuing high 
interest~\cite{Liu:2019vsn}. 
For tables for all $x$ ranges and for all data plotted in Figs.~\ref{fig:C1+C2}-\ref{fig:P7}
see the Supplemental Material at [URL will be inserted by publisher]. 

\section{Conclusion}

In this paper, we present the most detailed study of quark-hadron duality in
the spin structure function $g_1$ to date, for both the proton and the neutron. We study several different formulations of duality,
and find that duality seems to hold much better (at smaller momentum
transfer) in some cases than in others. In particular, we conclude
the following:
\begin{itemize}
    \item When forming integrals over kinematic regions 
    corresponding to specific resonance peaks, we observe good
    agreement between the measured data and the extrapolation from
    pQCD fits whenever
     several resonances with different
    spins contribute, {\it i.e.} the for the regions $W = 1.38$ GeV...1.58 GeV 
    (including the N$(1440)1/2^+$, N$(1520)3/2^-$, and 
N$(1535)1/2^-$ resonances) and $W = 1.58$ GeV...1.82 GeV (several
higher-lying resonances). In contrast, in the region dominated by the ground-state Delta-resonance ($W =$ 1.08 GeV...1.38 GeV) and the region $W = 1.82$ GeV...2 GeV with several Delta resonances, we
observe a much slower approach of the measured integrals towards
the extrapolated PDF fits with $Q^2$. This is likely due to the fact
that, at least at moderately low $Q^2$, for the excitation of Delta resonances
the transition to final-state helicity 3/2 dominates.
\item If we integrate over the entire resonance region up to $W=2$,
including the elastic peak at $W = 0.938$ GeV, a rather rapid convergence
towards the extrapolated PDF fits is observed: Global duality seems to work
in spin structure functions.
\item If instead we integrate over a fixed bin in $x$, with different
resonances contributing at different $Q^2$, we find that for lower
values of $x$, the transition with $Q^2$ towards a smooth scaling curve occurs 
only if the value of $Q^2$ is high enough so that $W > 2$ GeV. Conversely,
for the highest $x$ values, we observe that the approach towards a smooth
scaling curve (and the extrapolated PDF fits) occurs even below $W = 2$ GeV,
albeit at a higher $Q^2$. This may be due to the fact that at higher $Q^2$,
resonant and non-resonant contributions with different asymmetries average
out, leading to a ``precocious'' approach to scaling (or a different
form of local duality). This observation supports the idea that, for high
enough $x$ and $Q^2$, even data in the resonance region may be used
to constrain (polarized) parton distribution functions. 
Being able to include data in the resonance
region and {\it a fortiori} at moderate $W^2 < 10$ GeV$^2$ - a 
limit often imposed on PDF fits - will help with the goal
to pin down more precisely the quark polarization of both
types of valence quarks in the limit $x \rightarrow 1$, which
is still an open question at this time.
\end{itemize}

\acknowledgments{}
We would like to thank our collaborators on CLAS experiment EG1b 
and the members of the JAM collaboration for their help with the
data presented. We are particularly indebted to W. Melnitchouk for his valuable comments and to
J. Ethier, who provided us with the extrapolated JAM results.
This work was supported by the Department of Energy under Contract
DE-FG02-96ER40960 (ODU). The work of N.S. was supported by the DOE, Office of Science, Office of Nuclear Physics in the Early Career Program.


