(He^3) spectral function

The spectral function 𝒫 is parametrized by 3 scalar functions for each isospin state t = 0, 1 of the spectator system as 
(𝒫^t)(E,p,S)=(1/2)((Subscript[((f^t)), 0])+(Subscript[((f^t)), 1]) σ·S+(Subscript[((f^t)), 2]) \(T\&^\)(S)), 
where \(T\&^\)(S)=\(σ·p\&^\) \(S·p\&^\)-(1/3)σ·Sand S is the spin vector and \(p\&^\)=p/|p|.
The proton and neutron spectral functions:
\(TraditionalForm\`𝒫\^p = 𝒫\^\(t = 0\) + 𝒫\^\(t = 1\), \ 
𝒫\^n = 
        2  𝒫\^\(t = 1\)\)

The functions ((Subscript[f, i])^t)with i = 0, 1, 2 and t = 0, 1in terms of two-body and three-body break-up functions.
Isospin 0:
((Subscript[f, i])^0)(E, 
          p) = (Subscript[d, i])(p) δ((Subscript[E, B]) + E - (Subscript[E, d])) + ((Subscript[c, i])^0)(
            (Subscript[E, f]), p), 
where (Subscript[E, B]) = -7.72 MeVis the He-3 binding energy, (Subscript[E, d]) = -2.2246 MeV is the deuteron binding energy, E is the separation energy, and (Subscript[E, f]) = (Subscript[E, B]) + E is the energy of the two-nucleon continuum.
Isospin 1:
((Subscript[f, i])^1)(E, p) = ((Subscript[c, i])^1)((Subscript[E, f]), p). 
In this case the deuteron state does not contribute and the spectral function is determined by the two-nucleon continuum.

The energy integral of the spectral function 𝒫 gives the momentum distribution in the corresponding spin state.

The functions (Subscript[f, i]) in terms of the spectral function:
(Subscript[f, 0])= tr[𝒫]
(Subscript[f, 1])= tr[(Subscript[σ, z]) 𝒫]
(Subscript[f, 2])= \(\(9\/2\) tr[T\&^\ 𝒫]\)

The normalization of the spectral function ([dp]=((dEd)^3)p):
∫[dp]  tr[(𝒫^(p, n))(E,p,S)]=(Z, N), where Z = 2and N = 1for the proton and neutron.
∫[dp]  \(tr[\(σ\_z\) 𝒫\^\(p, n\)]\)=\(⟨σ\_z\^\(p, n\)⟩\), the average proton and neutron polatizations.
∫[dp]  \(tr[
                T\&^\ 𝒫\^\(p, n\)]\)=\(⟨T\&^\^\(p, n\)⟩\), the average proton and neutron tensor polatizations.

Define parameters and load spectral function

<< Graphics`Graphics`
<< Graphics`Graphics3D`
<< NumericalMath`NIntegrateInterpolatingFunct`

Off[NIntegrate::"slwcon"]

hc: 	MeV Fm conversion constant
hc2: 	Gev^2 Fm^2 conversion constant
mN:	Nucleon mass (MeV)
ed:	Deuteron binding energy (MeV)
eB:	Binding energy of He3 (MeV)

hc = 197.327053;
hc2 = (hc/1000. )^2 ;
mN = (938.27231 + 939.56563)/2.  ;
ed = -2.2246 ;
eB = -7.72;

sty = {{}, Hue[1], Hue[0.3]};

intops = {PrecisionGoal -> 4,  MinRecursion -> 3,  
          MaxRecursion -> 33};

SetDirectory["e:/ser/phys/spin/code/SpFn"]

Load SS spectral function

Get["ss_spfn.m"]

? d0 d1 d2 c0t0 c1t0 c2t0 c0t1 c1t1 c2t1

Load KPSV spectral function

Get["kpsv_spfn.m"]

?dkpsv ckpsv pup eup emn pmn

pup

eup

Normalizations and sum rules

?pmin pmax emin1 emax

pmin
pmax
emin1
emax

pmax = 6.5;

Normalization of spin-independent part:

SS:

nd0 = (4*Pi)*
          NIntegrate[p^2*d0[p], {p, pmin, 1. , pmax}, 
            intops // Evaluate];
nc0t0 = (4*Pi)*
          NIntegrate[
            p^2*c0t0[e, p], {p, pmin, 1. , pmax}, {e, emin1, 10. , emax}, 
            intops // Evaluate];
nc0t1 = (4*Pi)*
          NIntegrate[
            p^2*c0t1[e, p], {p, pmin, 1. , pmax}, {e, emin1, 10. , emax}, 
            intops // Evaluate];
n0p = nd0 + nc0t0 + nc0t1;
n0n = 2*nc0t1;

    
Print["No. protons = ", n0p, "\tNo. neutrons = ", 
      n0n, 
"\nnd0 = ", nd0, ", nc0t0 = ",  
      nc0t0, ", nc0t1 = ", nc0t1 ]

KPSV:

ndp = (4*Pi)*
          NIntegrate[p^2*dkpsv[0, p]*2, {p, pmn, 1. , pup}, 
            intops // Evaluate];
ncp = (4*Pi)*
          NIntegrate[
              p^2*ckpsv[0, 1, e, p]*2
, {p, pmn, 1. , 
                pup}, {e, emn, 10. , eup}, intops // Evaluate]/
            hc;
ncn = (4*Pi)*
          NIntegrate[
              p^2*ckpsv[0, 2, e, p]
, {p, pmn, 1. , 
                pup}, {e, emn, 10. , eup}, intops // Evaluate]/
            hc;

    
Print["No. protons = ", ndp + ncp, "\tNo. neutrons = ", 
      ncn, 
"\nndp = ", ndp, ", ncp = ",  
      ncp, ", ncn = ", ncn ]

Normalization of spin part:

SS:

nd1 = (4. *Pi)*
          NIntegrate[p^2*d1[p], {p, pmin, 1. , pmax}, 
            intops // Evaluate];
nc1t0 = (4. *Pi)*
          NIntegrate[
            p^2*c1t0[e, p], {p, pmin, 1. , pmax}, {e, emin1, 10. , emax}, 
            intops // Evaluate];
nc1t1 = (4. *Pi)*
          NIntegrate[
            p^2*c1t1[e, p], {p, pmin, 1. , pmax}, {e, emin1, 10. , emax}, 
            intops // Evaluate];
n1p = nd1 + nc1t0 + nc1t1;
n1n = 2*nc1t1;

    
Print["Proton polarization = ", 
      n1p, "\tNeutron polarization = ", 
      n1n, 
"\nnd1 = ", nd1, ", nc1t0 = ",  
      nc1t0, ", nc1t1 = ", nc1t1 ]

KPSV:

ndp = (4*Pi)*
          NIntegrate[p^2*dkpsv[1, p]*2, {p, pmn, 1. , pup}, 
            intops // Evaluate];
ncp = (4*Pi)*
          NIntegrate[
              p^2*ckpsv[1, 1, e, p]*2
, {p, pmn, 1. , 
                pup}, {e, emn, 10. , eup}, intops // Evaluate]/
            hc;
ncn = (4*Pi)*
          NIntegrate[
              p^2*ckpsv[1, 2, e, p]
, {p, pmn, 1. , 
                pup}, {e, emn, 10. , eup}, intops // Evaluate]/
            hc;

Print["Proton polarization = ", 
      ndp + ncp, "\tNeutron polarization = ", 
      ncn, 
"\nndp = ", ndp, ", ncp = ",  
      ncp, ", ncn = ", ncn ]

Normalization of tensor part:

SS:

nd2 = (4*Pi)*
          NIntegrate[p^2*d2[p], {p, pmin, 1. , pmax}, 
            intops // Evaluate];
nc2t0 = (4*Pi)*
          NIntegrate[
            p^2*c2t0[e, p], {p, pmin, 1. , pmax}, {e, emin1, 10. , emax}, 
            intops // Evaluate];
nc2t1 = (4*Pi)*
          NIntegrate[
            p^2*c2t1[e, p], {p, pmin, 1. , pmax}, {e, emin1, 10. , emax}, 
            intops // Evaluate];
n2p = nd2 + nc2t0 + nc2t1;
n2n = 2*nc2t1;

    
Print["Proton tensor pol. = ", 
      n2p, 
"\tNeutron tensor pol. = ", 
      n2n, 
"\nnd1 = ", nd2, ", nc2t0 = ",  
      nc2t0, ", nc2t1 = ", nc2t1 ]

KPSV:

ndp = (4*Pi)*
          NIntegrate[p^2*dkpsv[2, p]*2, {p, pmn, 1. , pup}, 
            intops // Evaluate];
ncp = (4*Pi)*
          NIntegrate[
              p^2*ckpsv[2, 1, e, p]*2
, {p, pmn, 1. , 
                pup}, {e, emn, 10. , eup}, intops // Evaluate]/
            hc;
ncn = (4*Pi)*
          NIntegrate[
              p^2*ckpsv[2, 2, e, p]
, {p, pmn, 1. , 
                pup}, {e, emn, 10. , eup}, intops // Evaluate]/
            hc;

    
Print["Proton tensor polarization = ", 
      ndp + ncp, 
"\tNeutron tensor polarization = ", 
      ncn, 
"\nndp = ", ndp, ", ncp = ",  
      ncp, ", ncn = ", ncn ]

Average kinetic and separation energies. Koltun sum rule.

SS:

kind0 = (4. *Pi)*
          NIntegrate[p^4*d0[p], {p, pmin, 1. , pmax}, 
              intops // Evaluate]/(2*mN)*
          hc^2;
kinc0t0 = (4. *Pi)*
          NIntegrate[
              p^4*c0t0[e, p]
, {p, pmin, 1. , pmax}, {e, 
                emin1, 10. , emax}, intops // Evaluate]/(2*mN)*
          hc^2;
kinc0t1 = (4. *Pi)*
          NIntegrate[
              p^4*c0t1[e, p]
, {p, pmin, 1. , pmax}, {e, 
                emin1, 10. , emax}, intops // Evaluate]/(2*mN)*
          hc^2;
kp = (kind0 + kinc0t0 + kinc0t1)/
          n0p;
kn = (2*kinc0t1)/n0n;

    
Print["Mean proton KE = ", kp, "\tMean neutron KE = ", 
      kn, 
"\nkind0 = ", kind0, ", kinc0t0 = ",  
      kinc0t0, ", kinc0t1 = ", kinc0t1 ]

sep0 = (ed - eB)*nd0 + (4. *Pi)*
            NIntegrate[
              p^2*(e - eB)*c0t0[e, p]
, {p, pmin, 1. , 
                pmax}, {e, emin1, 10. , emax}, 
              intops // Evaluate];
sep1 = (4. *Pi)*
          NIntegrate[
            p^2*(e - eB)*c0t1[e, p]
, {p, pmin, 1. , 
              pmax}, {e, emin1, 10. , emax}, 
            intops // Evaluate];
sp = (sep0 + sep1)/n0p;
sn = (2. *sep1)/n0n;
Print["Mean proton SE = ", sp, "\tMean neutron SE = ", 
      sn, 
"\nsep0 = ", sep0, ", sep1 = ",  
      sep1]

\(Print["Mean nucleon SE = ", \((meansep = \((sep0 + 3*sep1)\)/3; 
        meansep)\), 
"
Mean nucleon KE = ", \((meankin = \((kind0 + kinc0t0 + 3*kinc0t1)\)/3; 
        meankin)\), 
"
CM energy = ", 
      meankin/2, 
"
Koltum sum rule:
", "⟨K⟩-⟨E⟩-CM = ", 
      meankin - meansep - 
        meankin/2, 
\*""\n2(Subscript[E, B])/nucleon = "", 
      2*eB/3]\)

KPSV:

kdp = (4*Pi)*
          NIntegrate[p^4*dkpsv[0, p]*2, {p, pmn, 1, pup}, 
              intops // Evaluate]/(2*mN)*
          hc^2;
kcp = (4*Pi)*
          NIntegrate[
              p^4*ckpsv[0, 1, e, p]*2
, {p, pmn, 1, 
                pup}, {e, emn, 10. , eup}, intops // Evaluate]/(2*mN)*
          hc^2/hc;
kcn = (4*Pi)*
          NIntegrate[
              p^4*ckpsv[0, 2, e, p]
, {p, pmn, 1, pup}, {e, 
                emn, 10. , eup}, intops // Evaluate]/(2*mN)*
          hc^2/hc;
kp = (kdp + kcp)/2;
kn = kcn;

Print["Mean proton KE = ", kp, "\tMean neutron KE = ", 
      kn, 
"\nkdp = ", kdp, ", kcp = ",  
      kcp, ", kcn = ", kcn ]

sp = (ed - eB)*(4*Pi)*
            NIntegrate[p^2*dkpsv[0, p], {p, pmn, 1, pup}, 
              intops // Evaluate] + (4*Pi)*
            NIntegrate[
                p^2*(e - eB)*ckpsv[0, 1, e, p], {p, pmn, 1, pup}, {e, 
                  emn, 10, eup}, intops // Evaluate]/
              hc;
sn = (4*Pi)*
          NIntegrate[
              p^2*(e - eB)*ckpsv[0, 2, e, p], {p, pmn, 1, pup}, {e, emn, 
                10, eup}, intops // Evaluate]/
            hc;

Print["Mean proton SE = ", sp, "\tMean neutron SE = ", 
      sn]

Momentum distributions

(4*Pi)*
      NIntegrate[p^2*d0[p], {p, pmin, 1. , pmax}, 
        intops // Evaluate]
(4*Pi)*
      NIntegrate[p^2*dkpsv[0, p]*2, {p, pmn, 1. , pup}, 
        intops // Evaluate]

Plot[{d0[p], dkpsv[0, p]*2}
, {p, pmin, 
          1. }
, PlotRange -> All
, 
        PlotStyle -> sty];

Spin-independent distributions

Clear[mom0t0, mom0t1]

d0ss[p_?NumericQ] := 
        d0[p] + NIntegrate[c0t0[e, p], {e, emin1, 10. , emax}, 
            intops // Evaluate];
d1ss[p_?NumericQ] := 
        NIntegrate[c0t1[e, p], {e, emin1, 10. , emax}, 
          intops // Evaluate];
dpss[p_] := 4*Pi*(d0ss[p] + d1ss[p]);
dnss[p_] := 
        8*Pi*d1ss[p];

dpkpsv[p_?NumericQ] := 
        4*Pi*(dkpsv[0, p] + 
              NIntegrate[ckpsv[0, 1, e, p], {e, emn, 10. , eup}, 
                  intops // Evaluate]/hc);
dnkpsv[p_?NumericQ] := 
        4*Pi*NIntegrate[ckpsv[0, 2, e, p], {e, emn, 10. , eup}, 
              intops // Evaluate]/hc;

NIntegrate[p^2*dpss[p], {p, pmin, 1. , pmax}, 
      intops // Evaluate]
NIntegrate[p^2*dpkpsv[p], {p, pmn, 1. , pup}, 
      intops // Evaluate]
NIntegrate[p^2*dnss[p], {p, pmin, 1. , pmax}, 
      intops // Evaluate]
NIntegrate[p^2*dnkpsv[p], {p, pmn, 1. , pup}, 
      intops // Evaluate]

Plot[{dpss[p], dpkpsv[p]*2}
, {p, pmin, 
          1. }
, PlotRange -> All
, 
        PlotStyle -> sty];

Plot[{dnss[p], dnkpsv[p]}
, {p, pmin, 
          1. }
, PlotRange -> All
, 
        PlotStyle -> sty];

Spin distributions

mom1t0[p_?NumericQ] := 
        d1[p] + NIntegrate[c1t0[e, p], {e, emin1, 10. , emax}, 
            intops // Evaluate];
mom1t1[p_?NumericQ] := 
        NIntegrate[c1t1[e, p], {e, emin1, 10. , emax}, 
          intops // Evaluate];
mom1p[p_] := 4*Pi*(mom1t0[p] + mom1t1[p]);
mom1n[p_] := 8*Pi*mom1t1[p];

NIntegrate[p^2*mom1p[p], {p, pmin, 1. , pmax}, 
      intops // Evaluate]
NIntegrate[p^2*mom1n[p], {p, pmin, 1. , pmax}, 
      intops // Evaluate]

Plot[{mom1p[p]/2, 
          mom1n[p], (mom1p[p] + mom1n[p])/3}
, {p, 
          pmin, 1. }
, 
        PlotRange -> All
, 
        PlotStyle -> sty];

Plots

Momentum distributions:

Plot[{d0[p], d1[p], d2[p]}
, {p, pmin, 
        1.5}
, PlotRange -> All 
, 
      PlotStyle -> sty]

e = 1. ;
Plot[{c0t0[e, p], c1t0[e, p], c2t0[e, p]}
, {p, pmin, 
        1.5}
, PlotRange -> All 
, 
      PlotStyle -> sty ]
e =.

e = 10. ;
Plot[{ c0t1[e, p], c1t1[e, p] , c2t1[e, p]}
, {p, 
        pmin, 1.5}
, 
      PlotRange -> All   
, 
      PlotStyle -> 
        sty  
]
e =.

Energy distributions:

p = 0.2;
Plot[{c0t0[e, 
          p], (c1t0[e, p]), (c2t0[e, p])}
, {e, 
        emin1, 15. }
, 
      PlotRange -> All 
, 
      PlotStyle -> sty]
p =.

p = 0.2;
Plot[{c0t1[e, p], c1t1[e, p], c2t1[e, p]} 
, {e, 
        emin1, 5. }
, 
      PlotRange -> All  
, 
      PlotStyle -> sty]
p =.

Energy-momentum distributions:

\(Plot3D[
c0t0[e, p] (*\ \(-c1t1[e, p]\)\ *) , {p, 
        pmin, 1. }, {e, emin1, 10. }, PlotPoints -> 40, 
      AxesLabel -> {"p", "E", "Spectral Fn"}
, 
      PlotRange -> 
        All
 (*\ \(,\)\(Shading -> 
          False\)\ *) 
 (*\(,\)\(ViewPoint -> {1.2, 0.9, 
            1.1}\)*) 
]\)