Nucleon Structure Functions from the NJL-Model Chiral Soliton
I. Takyi H. Weigel







Institute for Theoretical Physics, Physics Department, Stellenbosch University,
Matieland 7602, South Africa


We present numerical simulations for unpolarized and polarized structure functions
in a chiral soliton model. The soliton is constructed self-consistently from
quark fields from which the structure functions are extracted. Central to the 
project is the implementation of regularizing the Dirac sea (or vacuum) contribution 
to structure functions from first principles. We discuss in detail how sum rules
are realized at the level of the quark wave-functions in momentum space. 
The comparison with experimental data is convincing for the polarized structure 
functions but exhibits some discrepancies in the unpolarized case. The vacuum
contribution to the polarized structure functions is particularly small.



Nucleon Structure Functions from the NJL-Model Chiral Soliton
    I. Takyi, H. Weigel
    December 30, 2023
=============================================================





§ INTRODUCTION


Perhaps the most convincing evidence for the quark substructure of baryons emerges
from Deep Inelastic Scattering (DIS). The conjuction of perturbative Quantum Chromo
Dynamics (QCD) and the parton model successfully explains the wealth of DIS data collected 
over the past decades <cit.>. However, these are not fully first principle 
calculations as the hadron wave-functions cannot (yet) be directly computed in QCD. 
Rather, in the spirit of the parton model, quark distributions are parameterized 
and subjected to perturbative QCD analysis <cit.>. On the other hand
there are many phenomenological models based on various aspects of QCD that (attempt to)
describe (static) properties of hadrons with particular focus on the nucleon. Examples
are the non-relativistic quark model <cit.>, the MIT bag model <cit.>
relativistic quark-diquark models <cit.>, chiral soliton models 
<cit.>; just to name a few. In principle any of these approaches 
should also be capable to predict nucleon structure functions. This is particular 
challenging for chiral soliton models that are formulated as bosonized action functionals 
hiding the quark substructure of hadrons. In this respect the Nambu-Jona-Lasino (NJL)
or chiral quark soliton model <cit.> is special: based
on a quark self-interaction the bosonization process can be traced step by step.

In the past nucleon structure functions have indeed been computed from the chiral quark 
soliton model and mainly two approaches were followed. Within the valence quark only 
approximation <cit.> the observation that though the Dirac sea 
(or vacuum) is essential to form the soliton, the vacuum contribution to nucleon properties 
is only moderate thereby justifying neglecting its contribution to the structure functions.
Though this approximation has empirical support, it is formally incomplete and does not follow 
from a systematic expansion scheme. In parallel, studies on the quark distributions were performed 
<cit.>
that identified the model quark degrees of freedom with those using the operator product expansion 
in the analysis of DIS. These studies included vacuum contributions. These are plagued 
by ultra-violet divergences requiring an a posteriori implementation of regularization
by a single Pauli-Villars subtraction imposed onto the distributions.  Unfortunately, this 
is not stringent since there are terms in the action that do not undergo 
regularization, e.g. to maintain the axial anomaly as measured by the decay of 
the neutral pion into two photons. Generally, the model has quadratic divergences (most notably
the gap equation) and a single subtraction may or may 
not be sufficient to remove all divergences. In the present project we therefore 
include the vacuum contributions to the nucleon structure functions (i) without identifying 
the NJL model quarks field with those of QCD, and (ii) by implementing the mandatory regularization 
already at the level of the defining action functional. The first issue is addressed by noting 
that the model emulates the chiral symmetry of QCD and thus produces the same symmetry currents,
in particular the electromagnetic one. To address the second issue we recall that DIS is 
described by the hadronic tensor W_μν which is the Fourier transform of the nucleon 
matrix element of the commutator of two electromagnetic current operators. Though there is 
no direct implementation of this commutator when bosonizing the NJL model, we take advantage
of its relation to the Compton tensor T_μν which itself is computed from the time-ordered 
product of these currents. Time-ordered products are straightforwardly included in a path integral
formalism within which bosonization is conducted. Regularizing this path integral by multiple 
Pauli-Villars subtractions proves most appropriate because it allows to trace the quark that 
carries the large momentum in the Bjorken limit. This formalism was developed already some
time ago <cit.> but its numerical simulation has been long outstanding. It will
be central to the study presented here.

This paper is organized as follows: In Section <ref> we introduce the NJL model with
emphasis on describing the regularization procedure in Minkowski space. The formulation in 
Minkowski space is advantageous to identify the absorptive part of T_μν which leads to 
W_μν. In Section <ref> we review the formalism from Ref. <cit.> 
of how to obtain the hadronic tensor in this model and in particular the role of the Bjorken 
limit. The NJL soliton description is explained in Section <ref>. We discuss the formalism 
to obtain the structure functions via the hadronic tensor from quark spinors that self-consistently
interact with the soliton in Section <ref>. Subsequently (Section <ref>)
we describe the way in which this formalism builds up the sum rules. Numerical results
are presented in Section <ref>. This analysis also includes the perturbative 
evolution to the scale at which experimental data are taken. As stressed above,
the model structure functions are identified from symmetry currents, not by 
equating QCD degrees of freedom. However, the evolution makes this unavoidable for 
the lack of any sensible alternative.  We briefly conclude and summarize in 
Section <ref>. Finally we leave technical details to four Appendixes.

Some preliminary results extracted from this paper have been put forward in
Ref. <cit.>




§ THE MODEL



We formulate the regularized action of the bosonized Nambu-Jona-Lasino (NJL) model
in Minkowski space as the sum of three pieces

    𝒜_ NJL   =𝒜_ R+𝒜_ I
    +1/4G∫ d^4x  tr[S^2+P^2+2m_0S]
    
    𝒜_ R   =-iN_C/2∑_i=0^2 c_i  Tr  log[- 𝐃𝐃_5 +Λ_i^2-iϵ] ,𝒜_ I   =-iN_C/2 Tr  log[-𝐃(𝐃_5)^-1-iϵ]  .

The subscripts R and I for real and imaginary refer to the respective properties 
after Wick-rotation to Euclidean space. In Minkowski space this corresponds to the 
use of two distinct Dirac operators

    𝐃   =-(S+γ_5P) + +γ_5 
    =:𝐃^(π)+ +γ_5 𝐃_5   =--(S+γ_5P) - +γ_5 
    =:𝐃_5^(π)+ +γ_5 .

Here, as in the local part of Eq. (<ref>), S and P are scalar and pseudoscalar 
fields that represent physical particles. Furthermore v_μ and a_μ are source fields. 
When expanding the action with respect to these sources the linear terms couple to the 
vector and axial vector currents. Since 𝒜_ I is (conditionally) finite 
in the ultra-violet, only 𝒜_ R undergoes regularization. Its dynamical 
content is quadratically divergent[The cosmological constant, computed from 
setting all fields to their vacuum expectation values, is quartically divergent and must
be subtracted when contributing.] and we thus require two subtractions. In the 
Pauli-Villars scheme they are implemented as

    c_0=1  ,  Λ_0=0  ,  ∑_i=0^2c_i=0
       and  ∑_i=0^2c_iΛ_i^2=0 .

In practice we will reduce the number of regularization parameters by assuming 
Λ_1  →Λ_2 = Λ which translates into the general prescription
∑_i c_i f(Λ_i^2)= f(0)-f(Λ^2)+Λ^2 f' (Λ^2 ). Nevertheless
we always write the regularization as in Eq. (<ref>).

The variation of the action with respect to the scalar field S yields the gap equation

    1/2G(m-m_0)
    =-4iN_C m∑_i=0^2c_i
    ∫d^4k/(2π)^4[-k^2+m^2+Λ_i^2-iϵ]^-1

that determines the vacuum expectation value of the scalar field ⟨ S⟩ =m.
Replacing S by ⟨ S⟩ in 𝐃 shows that m is the fermion mass and
is thus called the constituent quark mass. The chiral field U defines the non-linear 
representation for the isovector pion field π 

    S+ P=m U =m  exp[g/mπ·τ] .

We obtain the pion propagator by expanding 𝒜 to quadratic order and extract
the pion mass from the pole

    m_π^2=1/2Gm_0/m1/2N_CΠ(m_π^2)

and requiring a unit residuum at the pole fixes the quark-pion coupling

    1/g^2=4N_C d/dq^2[q^2 Π(q^2)]|_q^2=m_π^2  .

The above polarization function is

    Π(q^2)=∫_0^1 dx Π(q^2,x)      with    Π(q^2,x)=-i∑_i=0^2 c_i d^4k/(2π)^4 [-k^2-x(1-x)q^2+m^2+Λ_i^2-iϵ]^-2 .

Finally we get the pion decay constant f_π from the coupling to the axial current.
To this end we expand 𝒜 to linear order in both π and a_μ.
The result is f_π=4N_C m g Π(m_π^2). Using the empirical data m_π=138 MeV
and f_π=93 MeV fixes two of the three (m_0, G and Λ) model constants.
It is customary to use the constituent quark mass m as the single tunable parameter.



§ STRUCTURE FUNCTIONS IN THE NJL MODEL



The hadronic tensor of DIS is the matrix element of the commutator 
of electro-magnetic current operators

    W_μν(q;H) = 1/4π∫ d^4x  e^ q·ξ ⟨ H | [ J_μ(ξ),J_ν^†(0)]
    | H ⟩  ,

where q is the momentum of the virtual photon and H refers to either the pion or 
the nucleon target[The hadronic tensor can be written as the matrix element of the 
commutator for the lowest energy hadron in a given baryon number sector.]. 
This tensor is decomposed into Lorentz structures whose coefficients are form factors that
turn into the structure functions in the so-called Bjorken limit. Labeling the 
target momentum by p this limit is defined as

    Q^2=-q^2→∞     with  
    x=Q^2/2p· q   fixed .

Often x is referred to as the Bjorken variable.

By the optical theorem, W_μν is proportional to the absorptive part of the 
Compton amplitude

    W_μν(q;H)=1/2π  Im  T_μν(q;H)  .

The latter is the  matrix element of a time-ordered
product of the currents

    T_μν(q;H) = ∫ d^4ξ e^ q·ξ ⟨ H| T{J_μ(ξ)
    J_ν^†(0)}| H ⟩ .

With this relation we implement regularization from first principles because the time-ordered
product is immediately obtained from the action by functional differentiation

    T{J_μ(ξ) J_ν(0)}=
    δ^2/δ v^μ(ξ)δ v^ν(0)𝒜_ NJL(v)|_v_μ=0 ,

where v_μ is the photon field introduced by minimal coupling. Within the NJL model the 
photon couples to the quarks inside the hadron. As discussed comprehensively in 
Ref. <cit.> the evaluation of T_μν becomes feasible in the Bjorken 
limit. Then the quark propagator that carries the large photon momentum can be 
identified and thus be taken to be that of a free massless fermion. Thus the functional 
derivative from Eq. (<ref>) when applied to the real part simplifies to 
differentiating

    𝒜_Λ, R^(2,v)=
    -iN_C/4∑_i=0^2c_i
     Tr {(-𝐃^(π)𝐃^(π)_5+Λ_i^2)^-1[ Q^2()^-1𝐃^(π)_5
    -𝐃^(π)(()^-1)_5
     Q^2]}.

At this point it is important to explain the crucial role of the subscript `5' 
attached to the second term in square brackets of Eq. (<ref>). For this 
second term we have to recall that the (inverse) derivative operator in
()^-1 is actually
associated with the expansion of 𝐃_5. When comparing this
γ_5–odd operator to the ordinary Dirac operator in Eq. (<ref>)
one observes immediately that 𝐃_5 has a relative sign between the 
derivative operator i and the axial source γ_5. 
Therefore the axial–vector component of
(()^-1)_5
requires a relative sign.
With S_μρνσ=g_μρg_νσ
+g_ρνg_μσ-g_μνg_ρσ, that is

    γ_μγ_ργ_ν
    =S_μρνσγ^σ
    -iϵ_μρνσγ^σγ^5
       while  
    (γ_μγ_ργ_ν)_5
    =S_μρνσγ^σ+
    iϵ_μρνσγ^σγ^5  .

In Ref. <cit.> this modification was formally shown to be
consistent with the affected sum rules. In Section <ref> we see on the level
of the momentum space quark wave-functions that the structure functions
computed on the basis of Eq. (<ref>) indeed fulfill the sum rules.

Similarly, in the Bjorken limit, the imaginary part becomes

    𝒜_Λ, I^(2,v)=
    -iN_C/4 Tr {(-𝐃^(π)𝐃^(π))^-1[ Q^2()^-1𝐃^(π)_5
    +𝐃^(π)(()^-1)_5
     Q^2]}  .

These expression are still quite formal and we will use them to obtain nucleon 
structure functions in Section <ref>. We emphasize that these expressions
are directly deducted from the regularized action in Eq. (<ref>) and that 
no further assumption about the regularization has been made.

In Ref. <cit.> it has been shown that applying this formalism to the
pion relates its structure function to the spectral function from Eq. (<ref>)
as F(x)=5/9 (4N_C g^2) d/dp^2[p^2Π(p^2,x)]|_p^2=m_π^2, 
a result that was previously obtained from the analysis of the handbag diagram
in Refs. <cit.>.



§ NJL MODEL SOLITON



We construct the soliton from static meson configurations by introducing
a Dirac Hamiltonian h via

    𝐃^(π)=β(∂_t-h)    and  𝐃^(π)_5=(-∂_t-h)β  .

Its diagonalization

    hΨ_α = ϵ_αΨ_α  ,

yields eigen-spinors Ψ_α=∑_βV_αβΨ^(0)_α 
(Ψ^(0)_α are free Dirac spinors in a spherical basis, see Appendix 
<ref>) and energy eigenvalues ϵ_α.  
The hedgehog configuration minimizes the action in the unit baryon number sector and 
introduces the chiral angle Θ(r) via

    h=α·p +β  m  U_5(r) 
         where     
    U_5(r)= exp[r̂·τ γ_5 Θ(r)]  .

With the boundary conditions Θ(0)=-π and 
lim_r→∞Θ(r)=0 the diagonalization, Eq. (<ref>) yields
a distinct, strongly bound level, Ψ_ v, referred to as the valence 
quark level <cit.>. Its (explicit) occupation ensures 
unit baryon number. The functional trace in 𝒜_R is computed 
as an integral over the time interval T and a discrete sum over 
the basis levels defined by Eq. (<ref>). In the limit T→∞ 
the vacuum contribution to the static energy is then extracted 
from 𝒜_R→-TE_ vac. Collecting pieces, we obtain
the total energy functional as <cit.>

    E_ tot[Θ]=
    N_C/2[1+ sign(ϵ_ v)]
    ϵ_ v
    -N_C/2∑_i=0^2 c_i ∑_α{√(ϵ_α^2+Λ_i^2)
    -√(ϵ_α^(0)2+Λ_i^2)}
    +m_π^2f_π^2∫ d^3r  [1- cos(Θ)]  .

Here we have also subtracted the vacuum energy associated with the non-dynamical 
meson field configuration Θ≡0 (denoted by the superscript) that is 
often called the cosmological constant contribution. This subtraction will also 
play an important role for the unpolarized isosinglet structure function as it enters
via the momentum sum rule. The soliton profile is then obtained as the profile 
function Θ(r) that minimizes the total energy E_ tot self-consistently 
subject to the above mentioned boundary conditions.

This soliton represents an object which has unit baryon number but neither good quantum 
numbers for spin and flavor (isospin). Such quantum numbers are generated by canonically 
quantizing the time-dependent collective coordinates A(t) which parameterize the 
spin-flavor orientation of the soliton. For a rigidly rotating soliton the Dirac operator 
becomes, after transforming to the flavor rotating frame <cit.>,

    𝐃^(π)=Aβ(∂_t - 
    Ω·τ-h)A^†   and  𝐃^(π)_5=A(-∂_t + 
    Ω·τ -h)β A^†  .

Actual computations involve an expansion with respect to
the angular velocities 

    A^†d/dt A = /2Ω·τ ,

which, according to the quantization rules, are replaced by the spin operator 

    Ω⟶1/α^2 J  .

The constant of proportionality is the moment of inertia 

    α^2=N_C/4[1+ sign(ϵ_ v)]
    ∑_β v|⟨ v|τ_3|β⟩|^2/ϵ_β-ϵ_ v
    +N_C/8∑_αβ∑_i=1^2c_i
    |⟨α|τ_3|β⟩|^2/ϵ_α^2-ϵ_β^2{ϵ_α^2+ϵ_αϵ_β+2Λ_i^2/√(ϵ_α^2+Λ_i^2)
    -ϵ_β^2+ϵ_αϵ_β+2Λ_i^2/√(ϵ_β^2+Λ_i^2)} ,

which is of the order N_C. With Eq. (<ref>) 
the expansion in Ω is thus equivalent 
to the one in 1/N_C. The nucleon wave-function becomes a (Wigner D) 
function of the collective coordinates. A useful relation in 
computing matrix elements of nucleon states is <cit.>

    ⟨ N |1/2 tr(A^†τ_i Aτ_j) |N⟩ =
    -4/3⟨ N | I_i J_j | N⟩  .




§ STRUCTURE FUNCTIONS FROM SOLITON



We first repeat the relation between the structure functions and the 
hadronic tensor of the nucleon. Its symmetric combination
W_μν^S(q)=1/2(W_μν(q) + W_νμ(q)),
is parameterized by two form factors[Factors of the nucleon mass, M,
occur for dimensional reasons.]

    W_μν^S(q)     = M W_1(ν,Q^2) ( -g_μν
    + q_μ q_ν/q^2)+ W_2(ν,Q^2)/M( p_μ-p· q/q^2 q_μ)
    ( p_ν - p · q/q^2 q_ν) .

In the Bjorken limit, Eq. (<ref>), these form factors turn into the 
unpolarized structure functions that we extract by appropriate projections:

    F_1(x) = -1/2g^μν W_μν^S(q)
         and       F_2(x)=-xg^μνW_μν^S(q).

It must be noted that the Callan Gross relation, i.e. F_1(x)=2xF_2(x), is satisfied
in this case by construction. Similarly, the anti-symmetric part is also parameterized by 
two form factors

    W_μν^A(q)     = ϵ_μνλσ q^λ{ M G_1 (ν, Q^2) s^σ + G_2(ν, Q^2)/M( (p · q) s^σ - (q · s) p^σ)  }.

In the Bjorken limit these form factors yield the structure functions

    𝗀_1(x)= M^2ν G_1(ν, Q^2)
         and    𝗀_2(x) = M ν G_2(ν, Q^2) .

The longitudinal, 𝗀_1(x), and transverse,
𝗀_T(x)=𝗀_1(x) + 𝗀_2(x), structure functions are 
extracted from the hadronic tensor using the projection operators

    𝗀_1(x) = /2Mϵ^μνρσq_σ p_ρ/q · s W_μν^A(q)   and  q∥s
    𝗀_T(x) = -/2 Mϵ^μνρσ s_ρ p_σ
    W_μν^A(q)   and  q⊥s.


To obtain the hadronic tensor for the nucleon in the soliton model, the functional traces in 
Eqs. (<ref>) and (<ref>) are computed using the basis defined by the 
self-consistent soliton, Eq. (<ref>). This calculation has been detailed in 
Ref. <cit.> that we adopt directly. We start with the leading 
order in 1/N_C to the vacuum (or sea) contribution to W_μν

    W^ (s)_μν(q)     = - M N_cπ/8∫ dω/2π∑_α∫ d^3 ξ∫dλ/2 π  e^ Mx λ
       ×⟨ N,s| {[ 
    Ψ_α(ξ) 𝒬_A^2 γ_μγ_νΨ_α (ξ + λê_3) 
     e^-λω
    -Ψ_α (ξ) 𝒬_A^2 γ_νγ_μΨ_α (ξ - λê_3)  e^λω] . f_α^+ (ω)|_ p   +[Ψ_α (ξ) 𝒬_A^2 
    (γ_μγ_ν)_5 Ψ_α 
    (ξ - λê_3)  e^-λω 
    - Ψ_α (ξ ) 𝒬_A^2  
    ( γ_νγ_μ)_5 
    Ψ_α (ξ + λê_3)  e^λω] 
    . f_α^- (ω)|_ p}| N,s⟩ .
 
Here n^μ=(1,0,0,1)^μ is the light-cone vector of the 
photon momentum while 𝒬_A=A^†𝒬 A denotes
the flavor rotated quark charge matrix from which we compute nucleon 
matrix elements as in Eq. (<ref>). Furthermore

    f_α^± = ∑_i=0^2 c_i ω±ϵ_α/ω^2 - ϵ_α^2 - Λ_i^2 + ϵ±ω±ϵ_α/ω^2 - ϵ_α^2 +ϵ,
 
are Pauli-Villars regularized spectral functions. The subscript `p' indicates 
their pole contributions that we will explain below.

For the vacuum contribution to the isosinglet unpolarized structure function we then 
obtain

    [F_1^I=0(x)]_s     = 5 π/72 M N_c  ∫dω/2 π∑_α∫dλ/2 π e^ Mx λ( ∑_i=0^2 c_i 
    ω + ϵ_α/ω^2 - ϵ_α^2 - Λ_i^2 
    + ϵ)_ p   ×∫ d^3 ξ {Ψ_α^† (ξ) 
    (1 - α_3) Ψ_α (ξ + λê_3)  e^-ωλ 
    - Ψ_α^† (ξ) (1 -α_3) Ψ_α 
    (ξ - λê_3)  e^ωλ}.
 
Here the pole contributions is

    ( ∑_i=0^2c_i 1/ω^2-ϵ_α^2
    -Λ_i^2+ϵ)_ p    =  ∑_i=0^2 c_i
    -π/ω_α[ δ( ω + ω_α) 
    + δ( ω -ω_α) ] ,
 
where (eventually we take the single cut-off limit as described after 
Eq. (<ref>) and thus omit the label i on ω_α)

    ω_α = √(ϵ_α^2 + Λ_i^2) .
 
We recall that the single cut-off approach requires a derivative with 
respect to that cut-off. Of course, this also affects the implicit dependence 
of ω_α on that cut-off.

We introduce the Fourier transform of the quark wave-function as

    Ψ_α (p) = ∫d^3 r/4 π Ψ_α(r) 
     e^r·p.
 
Implementing a full Fourier transform differs from the approaches of Refs. 
<cit.> who used the expansion from diagonalizing 
the Dirac Hamiltonian, Eq. (<ref>). This resulted in discontinuities of the 
numerically computed quark distributions and required a smoothening procedure.

Performing the frequency (ω) and lambda[Technically it is 
advisable to average the photon direction rather than fixing it along 
the z-axis <cit.>.] `λ' integrals gives the 
vacuum contribution of the flavor-singlet unpolarized structure function in the 
nucleon rest frame (RF)

    [ F_1^I=0 (x)]_ s^∓    = 5MN_c/144∑_α∑_i=0^2 c_i  ∫_| M  x_α^±|^∞ p  d p ∫ d Ω_p {±Ψ_α^† (p) 
    Ψ_α (p) - ϵ_α/ω_αM x_α^±/pΨ_α^† (p) 
    p̂·αΨ_α (p) },
 
where   

    M x_α^±  = M x ±ω_α.
 
In the above [F_1^I=0(x)]_ s^± refers to the positive (negative) 
frequency components that are typically referred to as quark and antiquark distributions. 
In our calculation they arise from the two poles of the δ-function in 
Eq. (<ref>). Then the vacuum part of the isoscalar, unpolarized structure function 
becomes

    [F_1^I=0(x)]_ s = [F_1^I=0(x)]_ s^- + [F_1^I=0(x)]_ s^+.

As a matter of fact, this is still not the full result. Substituting free spionors 
(not interaction with the soliton) produces a non-zero result.  This non-zero result 
must also be subtracted. In the discussion of the sum rules below we will see that 
this is nothing but the ϵ^(0)_α type subtraction performed in 
Eq. (<ref>) and may be considered a cosmological constant type contribution.

The valence quark contribution is obtained by replacing the quark levels in 
(<ref>) by the cranked valence level 

    Ψ_ v^( rot) (r,t)     =  e^-ϵ_ vt A(t) 
    {Ψ_ v(r) + 1/2∑_α≠ vΨ_α (r) ⟨α | τ·Ω 
    |  v ⟩/ϵ_ v-ϵ_α} 
    = e^-ϵ_ vt A(t) ψ_ v(r) ,

In the above ψ_ v (r) is the spatial part of the valence quark 
wave-function with the rotational correction included and ϵ_ v is 
the energy eigenvalue of the valence quark level. Noting that the valence quark 
wave-function has positive parity and the pole contribution . f^±|_ pole = -4 πδ( ω∓ϵ_ v) 
gives the valence quark contribution 

    [ F_1^I=0(x) ]_ v^∓    = -5 M N_c/144[1+ sign(ϵ_ v)]
     ∫_M | x_ v^±|^∞ p  d p ∫ d Ω_p 
    {±Ψ_ v^† (p) Ψ_ v 
    (p) - M x_ v^±/pΨ_ v^† 
    (p) p̂·αΨ_ v (p)},

where
M x_ v^± = M x ±ϵ_ v.
Again, we have separated positive and negative frequency components.

The quark spinors Ψ_ν(r) separate into radial and angular pieces <cit.>. 
At the end, the structure functions, as in Eq. (<ref>) are computed as integrals 
over the (Bessel-)Fourier transforms of the radial functions in the quark spinors. In 
Appendix <ref> we list examples explicitly.

In quite an analogous manner, the isovector components of the polarized structure functions 
are extracted from the anti-symmetric combination W_μν^A(q).  Explicitly we find
the vacuum contribution to the longitudinal polarized structure function to be

    [𝗀_1^I=1(x)]_ s^∓   = 
    - M N_c/72 I_3∑_α∑_i=0^2 c_i{∓∫_| M x_α^±|^∞ 
     d p M x_α^±∫ d Ω_pΨ_α^† (p) p̂·τγ_5Ψ_α (p)      - ϵ_α/ω_α∫_| M x_α^±|^∞ d p p^2[ A_±∫ d Ω_pΨ_α^† (p) τ·σΨ_α (p)   
    +  B_±∫ d Ω_pΨ_α^† (p) 
    p̂·τp̂·σΨ_α 
    (p) ] },

and the isovector transverse polarized structure function as 

    [𝗀_T^I=1(x)]_ s^∓   =M N_c/144 I_3∑_α∑_i=0^2 c_i{ϵ_α/ω_α∫_| M x_α^±|^∞ d p p  ∫ d Ω_pΨ_α^† (p) τ·σΨ_α (p)     -ϵ_α/ω_α∫_| M x_α^±|^∞ 
    d p p^2[A_±∫ d Ω_pΨ_α^† (p) 
    τ·σΨ_α (p)+ B_±∫ d Ω_pΨ_α^† (p) p̂·τp̂·σΨ_α (p)]}.

In these formulas we have introduced the abbreviations, see also Eq. (<ref>)

    A_±    = 1/2p( 1 - (Mx_α^±)^2/p^2),
        B_±  = 1/2p( 3(Mx_α^±)^2/p^2 -1 ).

The total (vacuum contribution to the) polarized structure functions is the sum of the 
positive (+) and negative (-) frequency components. Again, some details in terms of 
the Fourier transformed radial functions are given in Appendix <ref>. For 
completeness we also list the formulas for the valence quark contribution <cit.>. 
The contribution to the longitudinal polarized structure function is obtained as

    [𝗀_1^I=1(x)]_ v^∓    = M N_c/72[1+ sign(ϵ_ v)]I_3{∓∫_| M x^±|^∞ d p M x_ v^±∫ d Ω_pΨ_ v^† (p)p̂·τγ_5Ψ_ v (p)     -∫_| M x_ v^±|^∞ d p p^2[ A_±∫ d Ω_pΨ_ v^† 
    (p) τ·σΨ_ v (p)  
    + B_±∫ d Ω_pΨ_ v^† (p) p̂·τp̂·σΨ_ v (p) ] },

and that for the transverse polarized structure function as 

    [𝗀_T^I=1(x)]_ v^∓    = M N_c/144[1+ sign(ϵ_ v)] I_3{∫_| M x_ v^±|^∞ d p p  ∫ d Ω_pΨ_ v^† (p) τ·σΨ_ v (p)      -∫_| M x_ v^±|^∞ d p p^2[ A_±∫ d Ω_pΨ_ v^† 
    (p) τ·σΨ_ v (p)  
    + B_±∫ d Ω_pΨ_ v^† (p) 
    p̂·τp̂·σΨ_ v (p) ] }.

The isovector unpolarized and isoscalar polarized structure functions are subleading 
in the 1/N_C counting. They are also more complicated to compute as they are quartic
in the spinors and involve double sums over the basis states defined Eq. (<ref>).
We refrain from presenting those lengthy expressions here and rather refer the interested 
reader to the Appendixes of Ref. <cit.>.



§ FORMAL DISCUSSION OF SUM RULES



In this section we discuss how the sum rules for the unpolarized and polarized structure 
functions work out when written explicitly in terms of the momentum space eigenspinors 
Ψ_α. In this context it is important to note that we compute 
the structure functions for a localized configuration in its rest frame. Then the 
Bjorken variable has support on the half axis from zero to infinity. Lorentz covariance
is regained by transforming to the infinite momentum frame, cf. Section <ref>.

Sum rules relate integrated structure functions to static observables. In soliton models 
the latter are directly expressed in terms of the eigenspinors, Eq. (<ref>) in 
coordinates space. Typically the sum rules can then be expressed as level-by-level identities. 
The only exception is the momentum (or energy) sum rule. For it to be obeyed it is compulsory 
that the soliton is an extremum of the energy functional, Eq. (<ref>).



 §.§ Momentum sum rule

For the momentum sum rule we require that 36/5∫ dx x F_1(x) produces
the quark contribution to the classical energy, i.e. all but the last integral 
in Eq. (<ref>). First we consider the scalar terms, 
±Ψ_α^† (p) Ψ_α(p),
from the vacuum contribution, Eq. (<ref>)[We adopt the 
notation
⟨α|Ô|α⟩_a 
= ∫_a^∞ d p p ∫ d Ω_p 
Ψ^†_α(p)ÔΨ_α (p).] 

    [ℳ_G^0]_ s    
    = M N_c/4∑_α∑_i=0^2 c_i ∫_0^∞ d x x 
    {⟨α|α⟩_| M x^+| 
    -⟨α|α⟩_| M x^-|}   =M N_c/4∑_α∑_i=0^2 c_i 
    {∫_ω_α/M ^∞d y (y -ω_α/M ) 
    ⟨α|α⟩_ M y  
    -∫_-ω_α/M ^∞ dy (y + ω_α/M ) 
    ⟨α| α⟩_| M y|},     =-N_c/2 ∑_α∑_i=0^2 c_i ω_α∫_0^∞ d y ⟨α| α⟩_ M y 
    =  N_c/2 ∑_α∑_i=0^2 c_i ω_α∫_0^∞ d y  
    y∂/∂ y∫_My^∞ d p p  ∫ d ΩΨ^†_α(p)Ψ_α(p),     =-M N_c/2∑_α∑_i=0^2 c_i ω_α∫_0^∞
    d y  y [ p ∫ dΩ_p 
    Ψ_α^†(p)
    Ψ_α(p)]_p=My =-N_c/2M∑_α∑_i=0^2 
    c_i √(ϵ_α^2+Λ_i^2),

which is 1/M times the vacuum contribution to the classical 
energy. This contribution also includes subtraction of the trivial vacuum 
energy, when there is no soliton. Hence the isoscalar unpolarized structure 
function necessitates the analog subtraction, as indicated earlier.
For the valence contribution the momentum sum rule gives

    [ℳ_G^0]_ v = 36/5∫_0^∞ dx  
    x[F_1^I=0(x)]_ v 
    =N_C/2M[1+ sign(ϵ_ v)] ϵ_ v .

Similarly, integrating the term with the operator
p̂·α gives 

    [ ℳ_G^1]_ s 
    =-M^2 N_c/2∑_α∑_i=0^2 c_i ϵ_α/ω_α∫_0^∞ d y y^2 ⟨α|1/pp̂·α|α⟩_ M y
    =-N_c/6M∑_α∑_i=0^2c_iϵ_α/ω_α⟨α|α·p|α⟩ .
 
Next we use the Dirac Hamiltonian, Eq. (<ref>) to write

    α·p=[r·p,h]
    -mβ[r·p,U_5(r)]

so that 
⟨α|α·p|α⟩=
 mβ⟨α|[r·p,U_5(r)]|α⟩.
Since r·p is the dilatation operator this matrix element 
measures the change of the single particle energy when scaling the soliton extension 
by an infinitesimal amount. Furthermore ϵ_α/ω_α=
∂/∂ϵ_α√(ϵ^2_α+Λ^2_i) so that 
[ℳ_G^1]_ s is the change of the vacuum energy 
when the soliton extension deviates slightly from its stationary point. 
Similarly, the valence quark adds
[ℳ_G^1]_ v=N_C/12M[1+ sign(ϵ_ v)]
⟨ v|[r·p,U_5(r)]| v⟩
to the sum rule. Then
[ℳ_G^1]_ s 
+ [ℳ_G^1]_ v is the coefficient of (λ-1) 
term in the expansion 

    E[U(λx)]=E_0 +(λ-1)E_1 +⋯(λ-1)^l E_l + ⋯

of the classical energy. Since U(x) is a stationary point, E_1=0 
thus verifying the momentum sum rule <cit.>. Obviously the 
momentum sum rule is not saturated level by level; rather it requires summing 
all contributions to this isoscalar unpolarized structure function. Hence this 
sum rule will be a very sensitive test of the numerical simulation.



 §.§ Bjorken sum rule

Here we verify the Bjorken sum rule in our model, which relates the 
isovector polarized structure function 𝗀_1^I=1
to the axial charge <cit.>. First, we show that the term 
in Eq. (<ref>) with the operator 
p̂·τγ_5 integrates to zero

    ∑_α∑_i=0^2 c_i ∫_0^∞d x {M x^- ⟨α|1/pp̂·τγ_5 |α⟩_| M x^- |-M x^+ ⟨α|1/pp̂·τγ_5 |α⟩_| M x^+ |},    =∑_α∑_i=0^2 c_i{∫_-ω_0/M^∞d y  My ⟨α|1/pp̂·τγ_5 |α⟩_| My| 
    -∫_ω_0/M ^∞d y  My⟨α|1/pp̂·τγ_5 |α⟩_My},   =∑_α∑_i=0^2 c_i 
    ∫_-ω_0/M^ω_0/M  d y  My⟨α|1/pp̂·τγ_5 |α⟩_| M y |=0 .
 
There are two contributions without γ_5. The first one contributes

    M N_c I_3/72∑_α∑_i=0^2 c_i 
    ϵ_α/ω_α∫_0^∞ d x [
    ⟨α|pA_+τ·σ|α⟩_| M x^+ |
    +⟨α|pA_-τ·σ|α⟩_| M x^- |]   
    =M N_c I_3/144∑_α∑_i=0^2c_i ϵ_α/ω_α[∫_ω_0/M^∞d y⟨α|τ·σ(1-(My)^2/p^2)|α⟩_ M y
    +∫_-ω_0/M ^∞ dy ⟨α|
    τ·σ(1-(My)^2/p^2)
    |α⟩_| M y |]   =N_c I_3/108∑_i=0^2 c_i ϵ_α/ω_α⟨α|τ·σ|α⟩ .

The term with p̂·τp̂·σ disappears because

    ∫_0^∞ d y∫_My^∞ p^2 d p ∫ dΩ_p 
    (1/p-3(My)^2/p^3)
    Ψ_α^†(p)p̂·τp̂·σΨ_α (p)   =∫_0^∞ dy∫_My^∞ dp ∫ dΩ_p
    ∂/∂ y(py-M^2y^3/p)
    Ψ_α^†(p)p̂·τp̂·σΨ_α (p)   =∫_0^∞ d y (M^2y^3/p-py) 
    [∫ dΩ_p Ψ_α^†(p)p̂·τp̂·σΨ_α(p)]_p=My=0 .

Hence the Bjorken sum rule for the vacuum contribution of the longitudinal polarized 
structure function becomes

    ∫ d x [𝗀_1^p(x) -𝗀_1^n(x)]_ s 
    =N_c /108∑_i=0^2 c_i ϵ_α/ω_0⟨α|τ·σ|α⟩ 
    =1/6[N_C /18∑_α∑_i=0^2 c_iϵ_α/√(ϵ_α^2 + Λ_i^2)⟨α |γ_3γ_5τ_3 | α⟩] .

The object in square brackets is the vacuum contribution to the axial 
charge <cit.>. Similar 
calculations from the valence contribution give

    ∫ d x [𝗀_1^p(x) -𝗀_1^n(x)]_ v 
    = -N_c/54[1+ sign(ϵ_ v)]
    ⟨ v|τ·σ| v⟩= 
    1/6[-N_C/9[1+ sign(ϵ_ v)]⟨ v| 
    γ_3γ_5τ_3|  v⟩]

with the object in square brackets being the valence quark contribution to the axial
charge. This indeed verifies the Bjorken sum rule for the total axial charge.

In an analog, yet much more tedious, calculation the sum rules for the subleading
contributions in the 1/N_C expansion are also verified via level by level 
identities. Details may be found in Ref. <cit.>. We would like to mention 
however, that the Adler sum rule <cit.>, which concerns a structure 
function from neutrino interactions and thus the exchange of a W gauge boson (not 
considered here), measures the isospin of the nucleon. In that case the sum rule is 
not level by level; rather summing this integrated structure function over all levels 
reproduces the moment of inertia, Eq. (<ref>) <cit.>.



§ NUMERICAL RESULTS



In this Section we present our numerical results for the structure functions. These 
results are obtained in a number of subsequent steps. First we construct the coordinate 
space eigenspinors of the self-consistent chiral soliton as described in Section <ref> 
for the parameters listed at the end of Section <ref>. In the second step we 
evaluate the Fourier transform according to Eq. (<ref>). Details of this 
transformation are provided in Appendix <ref>. Essentially the spinors in
momentum space are combinations of spherical harmonic functions of the unit momentum
vector and momentum space radial functions that are Bessel transforms of the 
radial functions in the coordinate space spinors from Section <ref>.
In momentum space the spherical harmonic functions combine to the conserved 
grand spin just as do those in coordinate space. Hence we formally obtain matrix 
elements of operators as, for example α·p̂, in the
very same way as the matrix elements of α·r̂ in 
coordinate space. In the third step the momentum space radial functions are numerically 
integrated to produce the structure functions in the nucleon rest frame. In the next
step they are transformed to the infinite momentum frame <cit.> 
and subsequently the standard (perturbative QCD) evolution to the scale of the 
experimental data is performed to allow for a sensible comparison. We note that
this evolution brings into the game a new model parameter, the scale at which
the evolution commences. We take a single scale for all structure functions.

We test the outcome of our numerical simulations via the sum rules, that is, we compare 
the integrated functions with associated local quantity obtained from the coordinate 
space spinors. To gain acceptable agreement a very fine (equi-distant) grid for 
the radial variable in momentum space is required. Typically we take several 
thousand points on an interval between zero and ten times the physical cut-off, 
Λ. Needless to say that this consumes a large amount of CPU time and 
obtaining (in particular the subleading 1/N_C contributions to) the 
structure functions takes days or weeks on an ordinary desktop PC. Still, there 
are minor numerical inaccuracies as reflected by small oscillations of the structure 
functions around a central value at larger x , cf. figures below.
Working in momentum space, rather than using the expansion coefficients 
V_αβ introduced after Eq. (<ref>) has, however, the advantage that no 
smearing <cit.> procedure is required. 



 §.§ Rest frame results



In Figures <ref> and  <ref> we show typical results for the isoscalar and isovector 
components of the unpolarized structure functions respectively. In this case they have
been obtained using the constituent quark mass of m=400  MeV. We separately show 
the contributions of the discrete valence level, those of the vacuum contributions and their 
sums (labeled as total). 


For the vacuum contribution we find the unexpected result that it dominates the 
valence counterpart. Mainly this originates from the (additional) subtraction
of the non-soliton piece mentioned after Eq. (<ref>). Without that 
subtraction we would not get a finite result, of course. Neither would the
momentum sum rule be fulfilled. However, this piece does not connect to the 
soliton rest frame and it is not clear at all whether or not transformation
of the Bjorken variable should be performed before taking the difference 
between the soliton and non-soliton isoscalar unpolarized structure functions.
Therefore we do not attach much relevance to this large vacuum contribution.
As expected, the vacuum contribution is sub-dominant for the isovector unpolarized 
structure function.


In Figure <ref> we present the unpolarized structure function that enters the 
Gottfried sum rule, i.e. 
F_2^p(x)-F_2^n(x) =2x [F_1^p(x)-F_1^n(x)] as the 
Callan-Gross relation holds in the soliton rest frame. The vacuum contribution turns 
slightly negative at large x which persists when adding the dominating valence piece
to form the total contribution of this structure function. In table <ref> we 
compare our results for the Gottfried sum rule, 
𝒮_G= ∫_0^∞d x/x (F_2^p-F_2^n), for various 
constituent quark masses to the experimental data from the NM Collaboration <cit.>.
Under this integral, the vacuum part is even less significant as its positive and negative
parts compensate. In total, the agreement for the Gottfried sum rule is surprisingly good since 
usually chiral soliton models reproduce empirical data with 30% accuracy <cit.>.


In Figures <ref> and <ref> we show the isoscalar and isovector contributions to 
the longitudinal polarized structure functions 𝗀^I=0,1_1, respectively. In both 
pictures we display the valence and vacuum contributions as well as their sums. Also 
the positive and negative frequency components of the valence and vacuum parts are shown. 
Obviously these structure functions are indeed dominated by their valence contributions and we thus 
a posteriori verify the valence only approximation adopted in Ref. <cit.>.
The only exception is the isoscalar structure function 𝗀_2 for which the valence 
contribution is small by its own due to large cancellations in 
𝗀_2=𝗀_T-𝗀_1. We also recognize some minor oscillations in the 
vacuum contributions at larger x. These occur as remnants of numerical inaccuracies.
We also remark that the present valence quark results do not exactly match those from 
Ref. <cit.> in which a soliton profile from the proper-time regularization
scheme was employed.



We have computed the axial vector and singlet charges on one hand side via the 
respective sum rules, i.e. by integrating the structure functions g^I=1,0_1(x)
and on the other hand via the coordinate space matrix elements of τ_3γ_3γ_5 
and γ_3γ_5 as e.g. in Eqs. (<ref>) and (<ref>).
The comparison in table <ref> serves as test for the numerical accuracy 
which works perfectly in the vector case while some minor discrepancies are observed for 
the axial singlet charge. This is understood as the latter is actually quartic in the
quark wave-functions (two of which are Fourier transformed) and it is also a double
sum over those wave-functions. So even tiny numerical errors  are amplified. 

As is typical for chiral soliton models, the axial vector charge falls short off the 
measured datum by about 30-40% <cit.>. It has been argued that this could be 
remedied by 1/N_C corrections arising from a particular handling of the collective 
coordinate quantization <cit.>. However, these corrections do 
not emerge in the current approach and also lead to inconsistencies with PCAC <cit.>.
On the other hand, the predicted axial singlet charge, which is linked to the proton spin 
problem <cit.>, is well within the errors of the empirical value.









 §.§ Projection and evolution



The soliton picture for baryons employs a localized field configuration which generally
breaks translational invariance. This causes the structure functions not to vanish when
x>1 as would be demanded kinematically. This effect is obvious in the above figures. 
We note that it is not limited to soliton models but is observed, i.e. in the bag
model as well <cit.>. Using light cone coordinates in the bag model in one 
space dimension a mapping of the structure functions from the localized field configuration 
was constructed that annihilated the structure functions for x>1 <cit.>. 
Guided by that construction a Lorentz boost was applied transforming the rest frame structure 
functions to the infinite momentum frame <cit.>

    f_IMF(x) = Θ(1-x)/1-x f_RF(-ln (1-x) ) ,
 
where f_RF refers to any of the structure functions computed in Section <ref>.
In what follows we will omit the label IMF for the boosted structure functions.

Even though we have adopted the high energy Bjorken limit in our kinematical analysis
of the Compton tensor, it must be emphasized that the NJL model is (at best) an approximation 
to QCD at the low mass scale, μ^2 which is thus a hidden parameter in the approach.
To compare with experimental data that are taken at higher energy scales, Q_ exp^2 we adopt 
Altarelli-Parisi (or DGLAP) equations <cit.>, 
to evolve the model structure functions accordingly. To be precise, we integrate

    f(x,t+δ t)=f(x,t)+δ td f(x,t)/d t .

with t = ln( Q^2/Λ^2_QCD) from Q^2=μ^2 to
Q^2=Q_ exp^2. The structure functions from Eq. (<ref>) are 
the initial values and we tune μ^2 for best fit at Q_ exp^2.

Since the isoscalar structure functions are associated with gluon type quantum numbers
they mix under the evolution. We take this into account under the assumption that
the gluon distributions vanish at μ^2. We consider the leading order of the perturbative 
expansion sufficient to estimate the quality of our results. Then the evolution equations 
have the following structure 

    d f^(I=1)(x,t)/d t   =g_QCD(t)/2 πC_R(F)∫_x^1d y/y
    P_qq(y) f^(I=1)(x/y,t),
    d f^(I=0)(x,t)/d t   =g_QCD(t)/2 πC_R(F)∫_x^1d y/y{P_qq(y) f^(I=0)(x/y,t)+P_qg(y) g(x/y,t)
    }, 
    d g(x,t)/d t   =g_QCD(t)/2 πC_R(F)∫_x^1d y/y{P_gg(y) g(x/y,t)+ P_gq(y)f^(I=0)(x/y,t) 
    } ,

where 
C_R(F)= (N_f^2-1)/2N_f is the color
factor for N_f flavors. Furthermore g_QCD(t) = 4π/β_0 t with
β_0=11/3N_C- 2/3N_f is the leading order 
perturbative running coupling constant. Explicit expressions for the splitting functions 
P_qq,…,P_gg, taken from Ref. <cit.> are listed in 
Appendix <ref> for completeness. From the evolved isoscalar and isovector 
components we finally obtain the proton and neutron structure functions as sum
and difference

    f^(p,n)(x,Q^2)=1/2[f^I=1(x,Q^2) ± f^I=0(x,Q^2)] .


We note that applying 
the perturbative QCD scheme to the model structure functions requires the identification
of model and QCD degrees of freedom even though there is no definite reason for doing
so other than the lack of any sensible alternative.

The second polarized structure function 𝗀_2 contains subleading, twist three, 
elements that undergo a different evolution procedure that is also described in 
Appendix <ref>.



 §.§ Comparison with experiment


As in previous calculations within the valence only approximation
<cit.> we take μ^2=0.4 GeV^2 
as initial value in the evolution differential equations. Smaller values contradict 
the perturbative nature of the evolution procedure as g_QCD(t)/2 π becomes 
sizable in view of Λ_QCD^2=0.2  GeV^2. On the other hand, significantly
larger μ^2 values worsen the agreement with experimental data.


In the left panel of Figure <ref> we show the numerical result for the polarized 
proton structure function 𝗀_1 obtained from the evolution equation at 
Q^2=3 GeV^2. We compare our results to experimental results from the E143 
Collaboration <cit.>. At small x the model results are 
somewhat larger than the data, but definitely the gross features are predominantly
reproduced.

For the neutron data are available in terms of the helium structure 
function <cit.>[In Ref. <cit.> direct neutron 
data are only given as the ratio 𝗀_1^n(x)/F_1(x).]

    𝗀_1^^3 He(x)≈ P_n 𝗀_1^n(x)
    +P_p 𝗀_1^p(x) -0.014[𝗀_1^p(x)-4𝗀_1^n(x)] ,

with P_n≈ 0.86 and P_p≈-0.028. From the right panel in Figure <ref>
we see that our model results reproduce the main features of the data: small positive 
values at large x turning negative at moderate x, though the minimum is more pronounced
by the model.



Next we discuss the results for the structure function 𝗀_2(x). As discussed
in Appendix <ref> the twist-2 and -3 pieces must be disentangled within the 
evolution whose result is shown in Figure <ref>. The effect of evolution is small 
for the twist-2 component but essential for the twist-3 counterpart.

When the end point of evolution is reached, the two components are combined to 
𝗀_2^p(x,Q^2). We display the model prediction in Figure <ref> and 
see that the data are well produced. This shows that the higher twist contributions 
cannot be neglected. This suggests that the higher twist contributions cannot be 
neglected.

Recently data were reported for the neutron twist-3 moment

    d_2^n(Q^2)=3∫_0^1 dx  x^2 𝗀^n_2(x,Q^2)

at two different transferred momenta: d_2^n(3.21 GeV^2)=-0.00421±0.00114 and
d_2^n(4.32 GeV^2)=-0.00035±0.00104 <cit.>, where we added the listed
errors in quadrature. For m=400 MeV the model calculation yields -0.00426 and
-0.00409, repsectively. While the lower Q^2 result nicely matches the observed value,
the higher one differs by about three standard deviations. This discrepancy as a function
of Q^2 indicates that the large N_C approximation to evolve 𝗀_2 (cf.
Appendix <ref>) requires improvement.


Finally, in Figure <ref> we display the unpolarized structure function that enters the Gottfried sum rule, 
i.e. F_2^p(x)-F_2^n(x) using the evolution equation. Though the negative 
contribution from the Dirac vacuum (cf. Figure <ref>) around x=1 is tiny in 
the rest frame, it gets amplified when transforming to the infinite momentum frame by the 
factor 1/(1-x) in Eq. (<ref>) thereby worsening the 
agreement with the experimental data collected by the NMC <cit.>.




§ CONCLUSION



We have presented the numerical simulation of nucleon structure functions within the 
NJL chiral soliton model. Central to this analysis has been the consistent implementation
of the regularized vacuum contributions that arise from all quark spinors being distorted
by the chiral soliton. Generally speaking, vacuum contributions should not be omitted
in any quark model as no expansion scheme suppresses them. This is even more the case
for the NJL model because the vacuum part significantly contributes to forming the soliton.

In our analysis we have only identified the symmetry currents of the model with those from QCD, 
not the quark distributions that are bilinear operators which are bilocal in the quark fields. 
Also, it is important to enforce the regularization on the action functional so that the 
regularization prescription for a given structure function is an unambiguous result. This 
increases the predictive power compared to previous similar studies that ”advocated” an 
ad hoc regularization of quark distributions 
<cit.>[Schwinger's proper time regularization
scheme is popular in the context of the NJL chiral soliton <cit.>.
Ref. <cit.> explicitly states that its application to the quark 
distributions is not yet known.].  A first principle regularization is particularly 
important when the sum rule for the structure function does not relate to coordinate space 
matrix elements of the quark fields. As an example we have seen that the isovector 
unpolarized structure functions are not subject to regularization (the explicit, 
lengthy formulas can be obtained from Ref. <cit.>). The prediction for 
the corresponding, so-called Gottfried, sum rule decently matches the empirical value. 
This is a very favorable case for the valence only approximation as the vacuum 
contribution essentially integrates to zero due to an unexpected negative contribution at 
large x. For the isoscalar part we recognized that the subtraction of the zero-soliton 
vacuum contribution has a sizable effect at small x. The emergence of this 
contribution is somewhat surprising as it suggests that the zero-soliton vacuum has 
structure. Yet it is required for convergence as well as fulfilling the momentum sum 
rule. We emphasize that we observe acceptable agreement for polarized proton structure 
functions between our model results and the experimental data. For the polarized structure 
functions our numerically expensive computation indeed showed that the vacuum contribution
is sub-dominant, except maybe for the isoscalar part of 𝗀_2 where the
valence part is tiny by itself. Nevertheless these results support the 
valence only approximation to a large extend.

In both, the unpolarized and polarized cases, the comparison with 
experiment required two additional operations on the model structure functions.
As the soliton is a localized field configuration, translational invariance is lost
and the rest frame structure functions must be Lorentz transformed to the infinite
momentum frame. In turn the results from that transformation are subject to the
perturbative QCD evolution scheme. This brings into the game the hidden parameter
at which to commence the evolution. We took that to be the same as in the 
valence only approximation.

Even though we have separated positive and negative frequency contributions to 
the structure functions we stop short of identifying them as (anti-)quark distributions
that parameterize semi-hard processes <cit.>, like e.g. 
Drell-Yan <cit.>. The reason being that we avoid to identify model 
and QCD quark degrees of freedom at the model scale.
There are many other nucleon matrix elements of bilocal, bilinear quark operators 
for which experimental results or lattice data are available. Examples are chiral odd 
distributions <cit.> or 
quasi-distributions <cit.>. 
It is challenging to see whether quark distributions, or at least some 
of them, can also be formulated and computed with a first principle 
regularization scheme in the NJL chiral soliton model.

This project is is supported in part by the National Research Foundation of South Africa (NRF)
by grant 109497. I. T. gratefully acknowledges a bursary from the Stellenbosch Institute
for Advanced Studies (STIAS).





§ SOLITON MATRIX ELEMENTS



The Dirac Hamiltonian h of the hedgehog field configuration (<ref>) commutes 
with the grand spin operator 

    G=J+ τ/2 = L 
    + σ/2 + τ/2,

which is the operator sum of the total spin J and the isospin τ/2. 
The total spin is the operator sum of the orbital angular momentum L and the 
intrinsic spin σ/2. Since h preserves G, the eigenfunctions of 
the Dirac Hamiltonian are also eigenfunctions of G. The quantum numbers of 
G are G^2=G(G+1) and G_3=M. The respective  eigenfunctions are 
tensor spherical harmonics associated with the grand spin

    [𝒴_LJGM(r̂) ]_si =  ∑_m, s_3, i_3, 
    J_3 C^GM_JJ_3, 1/2 i_3 C^JJ_3_Lm, 1/2 s_3Y_Lm(r̂) χ_s(s_3) χ_i(i_3)

where C^GM_JJ_3, 1/2 i_3 and C^JJ_3_Lm, 1/2 s_3 
are SU(2) Clebsch-Gordon coefficients that describe the coupling of  χ_s and 
χ_i, which are two components spinors and isospinors, respectively, and the spherical 
harmonic functions Y_Lm. 

For a prescribed profile function Θ (r) the numerical diagonalization 
of the Dirac Hamiltonian (<ref>) produces the radial functions g_α^(G,±;i) 
and f_α^(G,±;i)(i=1,2) that feature in eight component spinors <cit.>

    Ψ_α^(G,+)(r) = [   g_α^(G,+;1)(r) 𝒴_GG+1/2GM (r̂); f_α^(G,+;1)(r) 𝒴_G+1G+1/2GM (r̂) ]
    +
    [    g_α^(G,+;2)(r) 𝒴_GG-1/2GM (r̂); -f_α^(G,+;2)(r) 𝒴_G-1G-1/2GM (r̂) ]
    Ψ_α^(G,-)(r) = [ g_α^(G,-;1)(r) 𝒴_G+1G+1/2GM (r̂);  -f_α^(G,-;1)(r) 𝒴_GG+1/2GM (r̂) ]
    +
    [ g_α^(G,-;2)(r) 𝒴_G-1G-1/2GM (r̂);   f_α^(G,-;2)(r) 𝒴_GG-1/2GM (r̂) ].

Here, the second superscript (±) denotes the intrinsic parity defined by the parity 
eigenvalue as (-1)^G× (±1). The radial functions are written as linear combinations 
of spherical Bessel functions that build the free spinors Ψ_α^(0). The order of 
these Bessel functions matches the angular momentum label (first subscript) of the multiplying 
𝒴. The linear combination goes over momenta discretized by pertinent boundary 
conditions at a radius significantly larger than the extension of the profile function 
Θ (r). In Ref. <cit.> the condition that the radial function multiplying 
the 𝒴 with equal orbital angular momentum and grand spin indexes vanished at that 
large distance. In contrast, we impose that condition on the radial function of the upper 
component.  This avoids spurious contributions to the moment of inertia <cit.>.
 
Writing 

    e^p·r= 4π∑_Lm ()^L j_L(pr) 
    Y_Lm^* (r̂) Y_Lm(p̂).

we find the Fourier transform, Eq. (<ref>) of the spinors 

    Ψ_α^(G,+)(p) = ()^G+1[   g_α^(G,+;1)(p) 𝒴_GG+1/2GM (p̂); f_α^(G,+;1)(p) 𝒴_G+1G+1/2GM (p̂) ]
    +
    ()^G+1[   g_α^(G,+;2)(p) 𝒴_GG-1/2GM (p̂); f_α^(G,+;2)(p) 𝒴_G-1G-1/2GM (p̂) ]Ψ_α^(G,-)(p) = -()^G[ g_α^(G,-;1)(p) 𝒴_G+1G+1/2GM (p̂);   f_α^(G,-;1)(p) 𝒴_GG+1/2GM (p̂) ]
    +
    ()^G[   g_α^(G,-;2)(p) 𝒴_G-1G-1/2GM (p̂); ()^Gf_α^(G,-;2)(p) 𝒴_GG-1/2GM (p̂) ] .

The radial functions in momentum space are the Fourier-Bessel transforms

    ϕ_α (p) = ∫_0^∞ dr  r^2 j_L_α(pr) ϕ_α(r) ,

where L_α is the angular momentum associated with the coordinate space
radial wave-function ϕ_α(r). Note that the grand spin spherical harmonic
functions in momentum space are constructed precisely as in coordinate space, just
that the argument is the momentum space solid angle. Note that the intrinsic parity is 
also conserved quantum number. 

The valence quark carries G=0, then only the components 
with J=+1/2 are allowed for the eigenspinor   

    Ψ_α^0,+ (r)=Ψ_ v(r) = [ g_ v(r) 𝒴_0,1/2,0,0 (r̂); f_ v(r) 𝒴_1,1/2,0,0 (r̂) ]

here g_ v(r) =g_α^(0,+;1)(r) etc, are the particular eigenwave-functions. 
The cranking correction associated with the first order rotation (<ref>) 
dwells in the channel with G=1 and negative intrinsic parity 

    Ψ_α^(1,-)(r) = [  g_α^(1)(r) 𝒴_2,3/2,1,M (r̂); -f_α^(1)(r) 𝒴_1,3/2,1,M (r̂) ]
    +
    [ g_α^(2)(r) 𝒴_0,1/2,1,M (r̂); f_α^(2)(r) 𝒴_1,1/2,1,M (r̂) ],

for convenience we have written g_α^(1,-;1)(r) as g_α^(1)(r) etc. 
Taking the Fourier transform of equation (<ref>) gives

    ψ_ v (p)= Ψ_ v(p)+ ∑_α⟨ H_α⟩Ψ_α(p),

where

    Ψ_ v (p) = [ g_ v(p) 𝒴_0,1/2,0,0 (p̂); f_ v(p) 𝒴_1,1/2,0,0 (p̂) ]

and 

    Ψ_α (p) = -[  g_α^(1)(p) 𝒴_2,3/2,1,M (p̂) -
    g_α^(2)(p) 𝒴_0,1/2,1,M (p̂); f_α^(1) (p) 𝒴_1,3/2,1,M (p̂) -
    f_α^(2)(p) 𝒴_1,1/2,1,M (p̂) ].

The “matrix element” ⟨ H_α⟩ arises from perturbatively treating
the collective rotation

    ⟨ H_α⟩ = 1/2⟨α 
    | τ·Ω|  v⟩/ϵ_ v - ϵ_α.




§ UNPOLARIZED STRUCTURE FUNCTIONS AT LEADING ORDER



The level sums (over α) as e.g. in Eq. (<ref>) concern the 
label of the radial function, grand spin (G) and its projection (M). As we average 
of the direction of the virtual photon <cit.>, the matrix elements are 
degenerate in M. This produces the extra factor 2G+1 that we make explicit.


It is then straightforward to compute the matrix elements that appear 
in (<ref>): 

    ∫ d Ω_pΨ_α^† (p) Ψ_α 
    (p)   and
       ∫ d Ω_pΨ_α^† (p) p̂·αΨ_α (p) =  ∫ d Ω_pΨ_α^† (p) p̂·σγ_5Ψ_α (p) .

The positive intrinsic parity of the matrix element of (<ref>) 
is obtained as 

    ∫ d Ω_pΨ_α^† (G,+) (p) 
    Ψ_α^(G,+) (p)  = (2G + 1) 
    ( g_α^(G,+;1) (p)^2 + f_α^(G,+;1) 
    (p)^2 + g_α^(G,+;2)(p)^2 + f_α^(G,+;2) 
    (p)^2),
 
and for the negative intrinsic parity as 

    ∫ d Ω_pΨ_α^†(G,-) (p) 
    Ψ_α^(G,-) (p)  = (2G+1) 
    ( g_α^(G,-;1) (p)^2 + f_α^(G,-;1) 
    (p)^2 + g_α^(G,-;2)(p)^2 + f_α^(G,-;2) 
    (p)^2).

Here the overall factor (2G+1) arises from summing the grand spin projection 
contained in ∑_α. 

From Table (<ref>) the positive intrinsic parity of the matrix element 
of (<ref>) is obtained as 

    ∫ d Ω_p Ψ_α^† (G,+) (p) p̂·σγ_5Ψ_α^(G,+) (p)   =  -2(2G+1) 
    ( g_α^(G,+;1)(p) f_α^(G,+;1)(p)  
    + g_α^(G,+;2)(p) f_α^(G,+;2)(p) ),

and that for the negative intrinsic parity as  

    ∫ d Ω_p Ψ_α^† (G,-) (p) p̂·σγ_5Ψ_α^(G,-) (p)   = -2(2G+1) 
    ( g_α^(G,-;1)(p) f_α^(G,-;1)(p)  
    +g_α^(G,-;2)(p) f_α^(G,-;2)(p) ).

The matrix element from the valence contribution (<ref>) is easily 
obtained, using the definition of the decomposition of the valence wave function 
(<ref>). They are given as 

    ∫ d Ω_pΨ_ v^† (p) 
    Ψ_ v(p)  = g_ v(p)^2 
    + f_ v(p)^2  and
       ∫ d Ω_pΨ_ v^† (p) p̂·σγ_5Ψ_ v(p)  = 
    -2 g_ v(p) f_ v(p)

at leading order 1/N_C.



§ POLARIZED STRUCTURE FUNCTIONS AT LEADING ORDER



Here we list the matrix elements that appear in the vacuum contribution of the 
polarized structure functions, Eqs. (<ref>) and 
(<ref>).
The matrix element to be considered are 

    ∫ d Ω_pΨ_α^† (p) p̂·τγ_5Ψ_α (p) , 
    
       ∫ d Ω_pΨ_α^† (p) p̂·τp̂·σΨ_α (p) 
      and
       ∫ d Ω_pΨ_α^† (p) τ·σΨ_α (p).


The matrix element (<ref>) is computed from the matrix elements from 
Table <ref>: the positive intrinsic parity is obtained as

    ∫ d Ω_p Ψ_α^† (G,+) (p) p̂·τγ_5Ψ_α^(G,+) (p) 
        = 2(g_α^(G,+;1) (p) f_α^(G,+;1) (p) -  
    g_α^(G,+;2)(p) f_α^(G,+;2)(p) )     - 4 √(G(G+1))( g_α^(G,+;1)(p) f_α^(G,+;2)(p)  
    +  g_α^(G,+;2)(p)f_α^(G,+;1) (p) ) ,

and that for negative intrinsic parity as

    ∫ d Ω_p Ψ_α^† (G,-) (p) p̂·τγ_5Ψ_α^(G,-) (p) 
        = 2( g_α^(G,-;1) (p) f_α^(G,-;1) (p) - 
     g_α^(G,-;2)(p) f_α^(G,-;2)(p) )     + 4 √(G(G+1))( g_α^(G,-;1)(p) f_α^(G,-;2)(p)  
    +  g_α^(G,-;2)(p)f_α^(G,-;1) (p) ) .

Also the matrix element (<ref>) is computed from the matrix elements 
from Table <ref>:

the positive intrinsic parity contribution is given as

    ∫ d Ω_p Ψ_α^† (G,+) (p) p̂·τp̂·σΨ_α^(G,+) (p) 
        = (- g_α^(G,+;1) (p)^2 - f_α^(G,+;1) (p)^2  
    +g_α^(G,+;2) (p)^2 + f_α^(G,+;2) (p)^2)     + 4√(G(G+1))( g_α^(G,+;1)(p) g_α^(G,+;2)(p)  
    + f_α^(G,+;1)(p) f_α^(G,+;2) (p) ),

and that for the negative intrinsic parity as

    ∫ d Ω_p Ψ_α^† (G,-) (p) p̂·τp̂·σΨ_α^(G,-) (p) 
        = ( -g_α^(G,-;1) (p)^2 - f_α^(G,-;1) (p)^2  
    + g_α^(G,-;2) (p)^2 + f_α^(G,-;2) (p)^2)      - 4√(G(G+1))( g_α^(G,-;1)(p) 
    g_α^(G,-;2)(p)  + f_α^(G,-;1)(p) 
    f_α^(G,-;2) (p) ).

Furthermore the matrix element (<ref>) is computed from the matrix 
elements from Table <ref>:

the positive intrinsic parity contribution becomes

    ∫ d Ω_p Ψ_α^† (G,+) (p) 
    τ·σΨ_α^(G,+) (p) 
        = (2G+1) ( f_α^(G,+;1) (p)^2 
    + f_α^(G,+;2)(p)^2) - (2G+3) 
    g_α^(G,+;1)(p)^2    + (2G-1) g_α^(G,+;2) (p)^2  + 8 √(G(G+1))g_α^(G,+;1)(p) g_α^(G,+;2)(p) ,

and for negative intrinsic parity becomes

    ∫ d Ω_p Ψ_α^† (G,-) (p) 
    τ·σΨ_α^(G,-) (p) 
        = (2G+1) ( g_α^(G,-;1) (p)^2 + 
    g_α^(G,-;2)(p)^2) - (2G+3) 
    f_α^(G,-;1)(p)^2    - (2G-1) f_α^(G,-;2) (p)^2  - 8 √(G(G+1))f_α^(G,-;1)(p) f_α^(G,-;2) (p) .




§ SPLITTING FUNCTIONS



Here we list the  splitting functions used in Eqs. (<ref>), 
(<ref>) and (<ref>). They are 
different for the isovector, isosinglet and gluon contributions and are 
given as <cit.> 

    P_qq(z)   =1+z^2/(1-z)_+ ,
    P_qg(z)   =1/2C_R(F)[z^2 + (1-z)^2], 
    P_gq(z)   =1+(1-z)^2/z, 
    P_gg(z)   =2[1-z/z+z/(1-z)_++z(1-z)].

They determine the probability for the parton m to emit a parton n such that the 
momentum of the parton m is reduced by the fraction z. The regularized function 
(1-z)^-1_+ is defined under the integral by <cit.>

    ∫_0^1 d z f(z)/(1-z)_+=∫_0^1 d z f(z) -f(1)/1-z.

In the above C_R(F)= (N_f^2-1)/2N_f is the color 
factor for N_f flavors. Also the running coupling constant in the leading order is given 
by g_QCD(t) = 4π/β_0 t with 
β_0=11/3N_c- 2/3N_f being the coefficient of the leading 
term of the QCD β-function. Using the “+ ” prescription, the evolution equations 
for the isovector, isosinglet and gluon contributions become <cit.>

    d f^(I=1)(x,t)/d t   =2C_R(F)/9t{∫_x^1d y/y(1+y^2/1-y)[1/y f^(I=1)(x/y,t) - f^(I=1)(x,t) ]     + [ x + x^2/2 + 2ln(1-x) ] 
    f^(I=1)(x,t) },d f^(I=0)(x,t)/d t   =2C_R(F)/9t{∫_x^1d y/y((1+y^2/1-y)[1/y 
    f^(I=0)(x/y,t) - f^(I=0)(x,t) ]    + 3/4(y^2+(1-y^2)) g(x,t)/y) 
    + [ x + x^2/2 + 2ln(1-x) ] 
    f^I=0(x,t) }, d g(x,t)/d t   =2C_R(F)/9t{∫_x^1d y/y((1+(1-y)^2/y)f^(I=0)(x/y,t)    +9/2(1-y/y + y^2(1-y))
     g(x/y,t)/y + 9/2g(x/y,t)-g(x,t)/1-y)  
    +[3/2 + 9/2ln(1-x) ] g(x,t) }.

Now, since our NJL model calculations do not account for any gluon content in the 
nucleon, we assume in our numerical calculations that, at the initial boundary scale 
μ^2, the gluon content, f(x,t_0)=0 for both the polarized and unpolarized 
structure functions.

Unlike the polarized spin structure function 𝗀_1(x) of the nucleon and the unpolarized 
structure functions, the nucleon's second polarized spin structure function 𝗀_2 involves 
contributions from quark-gluon iterations and quark masses <cit.>. 
According to the standard operator product expansion analysis, these contributions come from the 
twist-3 local operators. It, however, also receives contribution from twist-2 local operators under 
the impulse approximation. Thus, the structure function 𝗀_2 can be written as the sum 
of the twist pieces 

    𝗀_2 (x,Q^2)=𝗀_2^WW(x,Q^2)+𝗀_2 (x,Q^2),
   
where the twist-2 piece is given as <cit.>

    𝗀_2^WW (x,Q^2)=-𝗀_1(x,Q^2) + ∫_0^11/y𝗀_1 (y,Q^2) d y ,

while that of the twist-3 piece is

    𝗀_2 (x,Q^2)=𝗀_1(x,Q^2)+𝗀_2(x,Q^2) 
    -∫_0^11/y𝗀_1 (y,Q^2) d y .

The twist-2 part undergoes the ordinary evolution as in Eq. (<ref>)
and the twist-3 piece is first parameterized by its moments

    M_j(μ^2)=∫ dx  x^j-1𝗀_2 (x,μ^2)

that scale as <cit.>

    M_j(Q^2)=[ln(μ^2)/ln(Q^2)]^γ_j-1/β_0     with    γ_j-1=2N_c[ψ(j)+1/2j+γ_E-1/4] .

Here ψ(j) is the logarithmic derivative of the Γ-function.
Then 𝗀_2 (x,Q^2) is obtained by expressing it in terms
of the evolved moments, i.e. by inverting Eq. (<ref>).

*0.5cm

61
natexlab#1#1bibnamefont#1#1bibfnamefont#1#1citenamefont#1#1url<#>1urlprefixURL 


[Muta(1987)]Mu87
authorT. Muta,
  titleFoundations of Quantum Chromodynamics
  (publisherWorld Scientific, Singapore, year1987).

[Roberts(1990)]Ro90
authorR. G. Roberts,
  titleThe Structure of the Proton
  (publisherCambridge Monographs on Mathematical Physics,
  year1990).

[Ynduráin(1993)]Yn93
authorF. J. Ynduráin,
  titleThe Theory of Quark and Gluon Interactions
  (publisherSpringer-Verlag, Berlin - Heidelberg - New York,
  year1993).

[Pumplin et al.(2002)Pumplin, Stump,
  Huston, Lai, Nadolsky, and Tung]Pumplin:2002vw
authorJ. Pumplin,
  authorD. R. Stump,
  authorJ. Huston,
  authorH. L. Lai,
  authorP. M. Nadolsky,
  and authorW. K. Tung,
  journalJHEP volume07, pages012
  (year2002), hep-ph/0201195.

[Kokkedee(1979)]Ko79
authorJ. J. J. Kokkedee,
  titleThe Quark Model (publisherBenjamin, New
  York, year1979).

[Thomas(1984)]Thomas:1982kv
authorA. W. Thomas,
  journalAdv. Nucl. Phys. volume13,
  pages1 (year1984).

[Maris and Roberts(2003)]Maris:2003vk
authorP. Maris and
  authorC. D. Roberts,
  journalInt. J. Mod. Phys. volumeE12,
  pages297 (year2003), nucl-th/0301049.

[Weigel(2008)]Weigel:2008zz
authorH. Weigel,
  journalLect. Notes Phys. volume743,
  pages1 (year2008).

[Alkofer et al.(1996)Alkofer, Reinhardt,
  and Weigel]Alkofer:1994ph
authorR. Alkofer,
  authorH. Reinhardt,
  and authorH. Weigel,
  journalPhys. Rept. volume265,
  pages139 (year1996), hep-ph/9501213.

[Christov et al.(1996)Christov, Blotz,
  Kim, Pobylitsa, Watabe, Meissner, Ruiz Arriola, and Goeke]Christov:1995vm
authorC. Christov,
  authorA. Blotz,
  authorH.-C. Kim,
  authorP. Pobylitsa,
  authorT. Watabe,
  authorT. Meissner,
  authorE. Ruiz Arriola,
  and authorK. Goeke,
  journalProg. Part. Nucl. Phys. volume37,
  pages91 (year1996), hep-ph/9604441.

[Weigel et al.(1997a)Weigel,
  Gamberg, and Reinhardt]Weigel:1996kw
authorH. Weigel,
  authorL. P. Gamberg,
  and
  authorH. Reinhardt,
  journalPhys. Lett. volumeB399,
  pages287 (year1997a),
  hep-ph/9604295.

[Weigel et al.(1997b)Weigel,
  Gamberg, and Reinhardt]Weigel:1996jh
authorH. Weigel,
  authorL. P. Gamberg,
  and
  authorH. Reinhardt,
  journalPhys. Rev. volumeD55,
  pages6910 (year1997b),
  hep-ph/9609226.

[Diakonov et al.(1996)Diakonov, Petrov,
  Pobylitsa, Polyakov, and Weiss]Diakonov:1996sr
authorD. Diakonov,
  authorV. Petrov,
  authorP. Pobylitsa,
  authorM. V. Polyakov,
  and authorC. Weiss,
  journalNucl. Phys. volumeB480,
  pages341 (year1996), hep-ph/9606314.

[Diakonov et al.(1997)Diakonov, Petrov,
  Pobylitsa, Polyakov, and Weiss]Diakonov:1997vc
authorD. Diakonov,
  authorV. Yu. Petrov,
  authorP. V. Pobylitsa,
  authorM. V. Polyakov,
  and authorC. Weiss,
  journalPhys. Rev. volumeD56,
  pages4069 (year1997), hep-ph/9703420.

[Pobylitsa et al.(1999)Pobylitsa,
  Polyakov, Goeke, Watabe, and Weiss]Pobylitsa:1998tk
authorP. V. Pobylitsa,
  authorM. V. Polyakov,
  authorK. Goeke,
  authorT. Watabe, and
  authorC. Weiss,
  journalPhys. Rev. volumeD59,
  pages034024 (year1999), hep-ph/9804436.

[Wakamatsu and Kubota(1998)]Wakamatsu:1997en
authorM. Wakamatsu and
  authorT. Kubota,
  journalPhys. Rev. volumeD57,
  pages5755 (year1998), hep-ph/9707500.

[Wakamatsu and Kubota(1999)]Wakamatsu:1998rx
authorM. Wakamatsu and
  authorT. Kubota,
  journalPhys. Rev. volumeD60,
  pages034020 (year1999), hep-ph/9809443.

[Wakamatsu(2003a)]Wakamatsu:2003wg
authorM. Wakamatsu,
  journalPhys. Rev. volumeD67,
  pages034005 (year2003a).

[Wakamatsu(2003b)]Wakamatsu:2002kq
authorM. Wakamatsu,
  journalPhys. Rev. volumeD67,
  pages034006 (year2003b),
  hep-ph/0212356.

[Weigel et al.(1999)Weigel, Ruiz Arriola,
  and Gamberg]Weigel:1999pc
authorH. Weigel,
  authorE. Ruiz Arriola,
  and authorL. P.
  Gamberg, journalNucl. Phys.
  volumeB560, pages383 (year1999),
  hep-ph/9905329.

[Takyi and Weigel(2019)]TakyiQNP18
authorI. Takyi and
  authorH. Weigel
  (year2019), noteContribution to the Proc. of Int.
  Conf. Quarks and Nuclear Physics, Tsukuba, Japan, Nov. 2018,
  1903.10748.

[Frederico and Miller(1994)]Frederico:1994dx
authorT. Frederico and
  authorG. A. Miller,
  journalPhys. Rev. volumeD50,
  pages210 (year1994).

[Davidson and Ruiz Arriola(1995)]Davidson:1994uv
authorR. M. Davidson and
  authorE. Ruiz Arriola,
  journalPhys. Lett. volumeB348,
  pages163 (year1995).

[Reinhardt(1989)]Reinhardt:1989st
authorH. Reinhardt,
  journalNucl. Phys. volumeA503,
  pages825 (year1989).

[Adkins et al.(1983)Adkins, Nappi, and
  Witten]Adkins:1983ya
authorG. S. Adkins,
  authorC. R. Nappi, and
  authorE. Witten,
  journalNucl. Phys. volumeB228,
  pages552 (year1983).

[Kahana and Ripka(1984)]Kahana:1984be
authorS. Kahana and
  authorG. Ripka,
  journalNucl. Phys. volumeA429,
  pages462 (year1984).

[Takyi(2019)]TakyiThesis
authorI. Takyi, Ph.D. thesis,
  schoolStellenbosch University (year2019),
  noteto be published on SUNScholar: scholar.sun.ac.za.

[Bjorken(1970)]Bjorken:1969mm
authorJ. D. Bjorken,
  journalPhys. Rev. volumeD1,
  pages1376 (year1970).

[Adler(1966)]Adler:1965ty
authorS. L. Adler,
  journalPhys. Rev. volume143,
  pages1144 (year1966).

[Weigel et al.(1996)Weigel, Gamberg, and
  Reinhardt]Weigel:1996ef
authorH. Weigel,
  authorL. P. Gamberg,
  and
  authorH. Reinhardt,
  journalMod. Phys. Lett. volumeA11,
  pages3021 (year1996).

[Gamberg et al.(1998)Gamberg, Reinhardt,
  and Weigel]Gamberg:1997qk
authorL. P. Gamberg,
  authorH. Reinhardt,
  and authorH. Weigel,
  journalInt. J. Mod. Phys. volumeA13,
  pages5519 (year1998), hep-ph/9707352.

[Arneodo et al.(1994)]Arneodo:1994sh
authorM. Arneodo et al.
  (collaborationNew Muon), journalPhys. Rev.
  volumeD50, pagesR1 (year1994).

[Wakamatsu and Watabe(1993)]Wakamatsu:1993nq
authorM. Wakamatsu and
  authorT. Watabe,
  journalPhys. Lett. volumeB312,
  pages184 (year1993).

[Christov et al.(1994)Christov, Blotz,
  Goeke, Pobylitsa, Petrov, Wakamatsu, and Watabe]Christov:1993ny
authorC. V. Christov,
  authorA. Blotz,
  authorK. Goeke,
  authorP. Pobylitsa,
  authorV. Petrov,
  authorM. Wakamatsu,
  and authorT. Watabe,
  journalPhys. Lett. volumeB325,
  pages467 (year1994), hep-ph/9312279.

[Alkofer and Weigel(1993)]Alkofer:1993pv
authorR. Alkofer and
  authorH. Weigel,
  journalPhys. Lett. volumeB319,
  pages1 (year1993), hep-ph/9308327.

[Aidala et al.(2013)Aidala, Bass, Hasch,
  and Mallot]Aidala:2012mv
authorC. A. Aidala,
  authorS. D. Bass,
  authorD. Hasch, and
  authorG. K. Mallot,
  journalRev. Mod. Phys. volume85,
  pages655 (year2013), 1209.2803.

[Barnett et al.(1996)]Barnett:1996hr
authorR. M. Barnett
  et al. (collaborationParticle Data Group),
  journalPhys. Rev. volumeD54,
  pages1 (year1996).

[Alexakhin et al.(2007)]Alexakhin:2006oza
authorV. Yu. Alexakhin
  et al. (collaborationCOMPASS),
  journalPhys. Lett. volumeB647,
  pages8 (year2007), hep-ex/0609038.

[Jaffe(1975)]Jaffe:1974nj
authorR. L. Jaffe,
  journalPhys. Rev. volumeD11,
  pages1953 (year1975).

[Jaffe(1981)]Jaffe:1980qx
authorR. L. Jaffe,
  journalAnnals Phys. volume132,
  pages32 (year1981).

[Gribov and Lipatov(1972)]Gribov:1972ri
authorV. N. Gribov and
  authorL. N. Lipatov,
  journalSov. J. Nucl. Phys. volume15,
  pages438 (year1972), note[Yad.
  Fiz.15,781(1972)].

[Altarelli and Parisi(1977)]Altarelli:1977zs
authorG. Altarelli and
  authorG. Parisi,
  journalNucl. Phys. volumeB126,
  pages298 (year1977).

[Dokshitzer(1977)]Dokshitzer:1977sg
authorY. L. Dokshitzer,
  journalSov. Phys. JETP volume46,
  pages641 (year1977), note[Zh. Eksp. Teor.
  Fiz.73,1216(1977)].

[Peskin and Schroeder(1995)]peskin1995introduction
authorM. E. Peskin and
  authorD. V. Schroeder,
  titleAn Introduction to Quantum Field Theory
  (publisherWestview Press, year1995).

[Weigel(2000)]Weigel:1998yy
authorH. Weigel,
  journalNucl. Phys. volumeA670,
  pages92 (year2000), hep-ph/9902390.

[Abe et al.(1995)]Abe:1994cp
authorK. Abe et al.
  (collaborationE143), journalPhys. Rev. Lett.
  volume74, pages346 (year1995).

[Abe et al.(1998)]Abe:1998wq
authorK. Abe et al.
  (collaborationE143), journalPhys. Rev.
  volumeD58, pages112003
  (year1998), hep-ph/9802357.

[Flay et al.(2016)]Flay:2016wie
authorD. Flay et al.
  (collaborationJefferson Lab Hall A), journalPhys.
  Rev. volumeD94, pages052003
  (year2016), 1603.03612.

[Abe et al.(1996)]Abe:1995dc
authorK. Abe et al.
  (collaborationE143), journalPhys. Rev. Lett.
  volume76, pages587 (year1996),
  hep-ex/9511013.

[Gribov et al.(1983)Gribov, Levin, and
  Ryskin]Gribov:1984tu
authorL. V. Gribov,
  authorE. M. Levin, and
  authorM. G. Ryskin,
  journalPhys. Rept. volume100,
  pages1 (year1983).

[Kenyon(1982)]Kenyon:1982tg
authorI. R. Kenyon,
  journalRept. Prog. Phys. volume45,
  pages1261 (year1982).

[Jaffe and Ji(1991a)]Jaffe:1991kp
authorR. L. Jaffe and
  authorX.-D. Ji,
  journalPhys. Rev. Lett. volume67,
  pages552 (year1991a).

[Jaffe and Ji(1992)]Jaffe:1991ra
authorR. L. Jaffe and
  authorX.-D. Ji,
  journalNucl. Phys. volumeB375,
  pages527 (year1992).

[Ji(2013)]Ji:2013dva
authorX. Ji, journalPhys.
  Rev. Lett. volume110, pages262002
  (year2013), 1305.1539.

[Alexandrou et al.(2017)Alexandrou,
  Cichy, Constantinou, Hadjiyiannakou, Jansen, Steffens, and
  Wiese]Alexandrou:2016jqi
authorC. Alexandrou,
  authorK. Cichy,
  authorM. Constantinou,
  authorK. Hadjiyiannakou,
  authorK. Jansen,
  authorF. Steffens, and
  authorC. Wiese,
  journalPhys. Rev. volumeD96,
  pages014513 (year2017), 1610.03689.

[Broniowski and
  Ruiz Arriola(2017)]Broniowski:2017wbr
authorW. Broniowski and
  authorE. Ruiz Arriola,
  journalPhys. Lett. volumeB773,
  pages385 (year2017), 1707.09588.

[Altarelli et al.(1994)Altarelli, Nason,
  and Ridolfi]Altarelli:1993np
authorG. Altarelli,
  authorP. Nason, and
  authorG. Ridolfi,
  journalPhys. Lett. volumeB320,
  pages152 (year1994), note[Erratum: Phys.
  Lett. B325, 538 (1994)], hep-ph/9311255.

[Jaffe(1990)]Jaffe:1989xx
authorR. L. Jaffe,
  journalComments Nucl. Part. Phys. volume19,
  pages239 (year1990).

[Ji and Chou(1990)]Ji:1990br
authorX.-D. Ji and
  authorC.-h. Chou,
  journalPhys. Rev. volumeD42,
  pages3637 (year1990).

[Jaffe and Ji(1991b)]Jaffe:1990qh
authorR. L. Jaffe and
  authorX.-D. Ji,
  journalPhys. Rev. volumeD43,
  pages724 (year1991b).

[Wandzura and Wilczek(1977)]Wandzura:1977qf
authorS. Wandzura and
  authorF. Wilczek,
  journalPhys. Lett. volume72B,
  pages195 (year1977).





