Nuclear convolution

Clear["Global`*"]
<< Graphics`

      
<< Graphics`MultipleListPlot`
<< Graphics`Legend`

Off[NIntegrate::"slwcon"];
  
Off[General::spell1];

Initialization parameters

\(\(\( (*\(PhysPath = "~/phys/";\)*) \)\(
\)\(\(PhysPath = "e:/ser/phys/";\)
 (*
PDFPath = 
        PhysPath <> "pdf/a02m/"; 

      PionPath = 
        PhysPath <> "pdf/pion/"; 
DeuteronPath = 
        PhysPath <> "deuteron/"; 
WFPath = 
        DeuteronPath <> "wf/"; 

      SpecFnPath = PhysPath <> "SpecFn/"; 

      ConvolutionPath = 
        PhysPath <> "convolution/";
*) 

    \(packagepath = PhysPath <> "spin/code/";\)\)\)\)

Options for numerical integration and graphics drawing

intopts = {MaxRecursion -> 30, 
          PrecisionGoal -> 4};

xmin = 1.*^-6;
xmax = 0.98;

colsty0 = {{}, RGBColor[1, 0, 0]
, 
          RGBColor[0, 1, 0], RGBColor[0, 0, 1]
, 
          RGBColor[1, 0.7, 0], RGBColor[1, 0, 1]
, 
          RGBColor[0, 1, 1], 
          RGBColor[0, 1, 0.5]};
colsty1 = {Dashing[{Dash}]
, {Dashing[{Dash}], 
            RGBColor[1, 0, 0]}
, {Dashing[{Dash}], 
            RGBColor[0, 1, 0]}
, {Dashing[{Dash}], 
            RGBColor[0, 0, 1]}
, {Dashing[{Dash}], 
            RGBColor[1, 1, 0]}};

symsty = {PlotSymbol[Box,  
            Filled  ->  False], 
PlotSymbol[
            Box], 
MakeSymbol[
            Disk[{0, 0}, 
              Offset[{2.5, 2.5}] ]], 
MakeSymbol[
            Circle[{0, 0}, 
              Offset[{2.5, 2.5}] ]], 
PlotSymbol[
            Triangle], 
PlotSymbol[
            Diamond], 
PlotSymbol[Star]};

hc: 	MeV Fm conversion constant
hc2: 	(Gev^2) (Fm^2) conversion constant
mN:	Nucleon mass in (Fm^(-1))

hc = 197.327053;
hc2 = (hc/1000. )^2 ;
mN = (938.27231 + 939.56563)/2. /
          hc ;
mN2 = mN^2 ;

Nucleon structure functions

Test PDFs

xuv[x_, ___] := 2.0*x^0.5*(1 - x)^3;
xdv[x_, ___] := x^0.5*(1 - x)^4;
xqval[x_, ___] := xuv[x] + xdv[x];
xqsea[x_, ___] := (24. /5. )*0.25*
          x^(-0.2)*(1 - x)^5;
f2[x_, ___] := 
        5. /18. *(xqval[x] + xqsea[x]);
f2p[x_, ___] := 
        4. /9.  *xuv[x] + 1. /9.  *xdv[x] + 
          4. /3.  * 0.25*
            x^(-0.2)*(1 - x)^5;
f2n[x_, ___] := 
        4. /9.  *xdv[x] + 1. /9.  *xuv[x] + 
          4. /3.  * 0.25*
            x^(-0.2)*(1 - x)^5;
xf1[x_, ___] := f2[x] /2. ;
xf3[x_, ___] := x^0.5*(1 - x)^4  ;

qsq = 10. ;
\(LogLinearPlot[{xuv[x, qsq], xdv[x, qsq], 
        f2[x, qsq]}
, {x, 1.0*^-4, 1}
, 
      PlotRange -> All
, 
      PlotStyle -> {Hue[1], Hue[0.3], {}}
, 
      PlotLabel -> \*""Nucleon PDFs at \!\(Q\^2\) = "" <> 
          ToString[qsq]
\ ]\)
qsq =. 
x =.

Target mass corrections

Clear[TMC2, TMC1, TMC3, TMC2e, TMC1e, TMC3e, x, q2, p2, y]

Georgi-Politzer TMC expressions in terms of Nachtmann variable:

tmcint[fLT_, x_?NumberQ, q2_, p2_] := 
        NIntegrate[fLT[y, q2, p2]/y^2,  {y, x, 1}, 
          PrecisionGoal -> 4];

TMC2[f0LT_, f1LT_, x_?NumberQ, q2_?NumberQ, 
          p2_?NumberQ]  := (tau = 1 + 4 p2 x^2/q2;  
          xi = 2 x/(1 + Sqrt[tau]); 
x^2*
              f0LT[xi, q2, p2]/(xi^2 tau^1.5)  +  
            6*p2*x^3/(q2*tau^2)* 
              tmcint[f1LT, xi, q2, 
                p2]
);
TMC1[xf1LT_, f2LT_, x_?NumberQ, q2_?NumberQ, 
          p2_?NumberQ]  :=  (tau = 1 + 4 p2 x^2/q2;  
          xi = 2 x/(1 + Sqrt[tau]); 
x^2*
              xf1LT[xi, q2, p2]/(xi^2* tau^0.5)  +   p2*
              x^3/(q2*tau)*
              tmcint[f2LT, xi, q2, 
                p2] 
);
TMC3[xf3LT_, f2LT_, x_?NumberQ, q2_?NumberQ, 
          p2_?NumberQ]  :=  (tau = 1 + 4 p2 x^2/q2;  
          xi = 2 x/(1 + Sqrt[tau]); 
x^2*
              xf3LT[xi, q2, p2]/(xi^2* tau)  +  
            2*p2*x^3/(q2*tau^1.5)* 
              tmcint[f2LT, xi, q2, p2]
);

Expansion of GP TMC expressions in powers of 1/(Q^2).
We keep only  1/(Q^2) term

TMC2e[f0LT_, f1LT_, x_, q2_, 
          p2_]  := (
(1 - 4*x^2*p2/q2)*
              f0LT[x, q2, p2]  +  
            x^3*p2/q2* (
6*tmcint[f1LT, x, q2, p2] - 
                  D[f0LT[x, q2, p2], 
                    x]
)
);
TMC3e[f0LT_, f1LT_, x_, q2_, 
          p2_]  := (
(1 - 2*x^2*p2/q2)*
              f0LT[x, q2, p2]  +  
            x^3*p2/q2* (
2*tmcint[f0LT, x, q2, p2] - 
                  D[f0LT[x, q2, p2], 
                    x]
)
);
TMC1e[f0LT_, f1LT_, x_, q2_, 
          p2_ ]  := (
f0LT[x, q2, p2]  +  
            x^3*p2/q2* (
2*tmcint[f1LT, x, q2, p2] - 
                  D[f0LT[x, q2, p2], 
                    x]
)
);

Target mass corrected structure functions:

f2TMC[x_, q2_, p2_]  :=  
        TMC2[f2, f2, x, q2, p2];

f2onTMC[x_, q_, p_]  :=  
        TMC2[f2, f2, x, q, mN2];

Deuteron wave function and form-factors

Get[ packagepath  <> "DeuteronWF.m" ]

? setdeuteron psi ed

setdeuteron[1];

Deuteron binding energy eD (in (fm^(-1))) and momentum distribution (Subscript[n, d])(k,  cosθ)=|(Subscript[Ψ, d])(k)\( | \^2\)

eD = ed ;

psi0[k_] := psi[0, 0, k];
psi2[k_] := psi[0, 2, k];

ndParis[k_, y_]  := psi0[k]^2 + psi2[k]^2 ;

Mean momentum squared and separation energy in the deuteron:

ksqD = 
        NIntegrate[ 
          k^4*ndParis[k, y], {k, 0, 
            7. }, {PrecisionGoal  ->  5}];
ekinD = ksqD/(2  mN);
epsD = eD - ekinD;
virD = 2*(epsD - ekinD)/mN;

Deuteron elastic form-factor SD

(*
zmn = 0;
zmx = 1;
SD[k_?NumericQ, wf2_]:=Block[{r = Log[1/z]},
NIntegrate[wf2[r]*Sin[k*r]/(k*r*z),  {z, zmn, zmx}, 
                MaxRecursion -> 30, PrecisionGoal -> 4 ]];
*)

Mean nucleon light-cone momentum in the deuteron:

1 + (epsD + 2/3*ekinD)/mN

virD

Off-shell corrections

Generic off-shell correction.
Note that dimensional quantities in off-shell functions are in Fm units.

off[sf_, δf_, x_, q2_, p2_]  := 
        sf[x, q2*
              hc2]*(1  +  (((p2 - mN2))/mN2)*δf[
                  x] );

Phenomenological off-shell function

x1 = 0.050  ; x2 = 0.448 ; 
    cN = 8.10;

deltaf[x_, ___] := 
        cN*(x - x1) (x - x2) (1 + x2 - x);

Off-shell structure function

f2off[x_, q2_, p2_] := off[f2, deltaf, x, q2, p2];

Combined TMC and off-shell corrections

f2offTMC[x_, q2_, p2_]  :=  
        TMC2[f2off, f2off, x, q2, p2];

Tests

Get[ packagepath <> "NuclearConvolution.m"]

?GenNConv2

?NConv2

Nuclear structure functions:

xgrid = Range[0.2, 0.9, 0.05];

qsq = 10. ;
q2 = qsq/hc2;

Table of SF values on the xgrid:

f2tbl = 
        Table[{
    #
, 
                f2[#, qsq]
, 
                f2onTMC[#, q2, mN2]
, 
                GenNConv2[2, ndParis, eD, mN, #, qsq, f2]
, 
                GenNConv2[2, ndParis, eD, mN, #, qsq, 
                  f2off]
, 
                GenNConv2[2, ndParis, eD, mN, #, qsq, 
                  f2offTMC]
} &@ xgrid[[i]], {i, 1, 
            Length[xgrid]}];

The same as a function of (W^2) = (M^2) + (Q^2)(1/x - 1)

tblw = 
        Table[{( w2 = mN2*hc2 + qsq*(1/#1 - 1); 
                  w2), #2, #3, #4 , #5, #6 } & @@ f2tbl[[i]], {i, 
            1, Length[f2tbl]}];

Ratios:

r200 = 
        Table[{#1, #4/#2} & @@ f2tbl[[i]], {i, 1, 
            Length[f2tbl]}];
r2off = 
        Table[{#1, #5/#2} & @@ f2tbl[[i]], {i, 1, 
            Length[f2tbl]}];
r2tot = 
        Table[{#1, #6/#3} & @@ f2tbl[[i]], {i, 1, 
            Length[f2tbl]}];
rw = 
        Table[{#1, #6/#3} & @@ tblw[[i]], {i, 1, 
            Length[tblw]}];

Plot the EMC ratios:

\(\(grx = 
        MultipleListPlot[r200, \ r2off, \ r2tot
, 
          PlotJoined -> \ True
, 
          SymbolShape -> None
, 
          PlotStyle\  -> \ 
            colsty0
 (*\(,\)\(PlotRange -> {{0.05, 
                  0.9}, {0.95, 1.15}}\)*) 
, 
          AxesLabel -> {x, ""}
, 
          PlotLabel -> \*""((Subscript[F, 2])^D)/((Subscript[F, 2])^N) at (Q^2)= "" <> ToString[
                qsq] <> "
Paris potential"
];\)\)

Plots of the EMC effect in D as a function of (W^2)

\(\(grw = 
        MultipleListPlot[rw
, 
          PlotJoined -> \ True
, 
          SymbolShape -> None
, 
          PlotStyle\  -> \ 
            colsty0
 (*\(,\)\(PlotRange -> {{0.05, 
                  0.9}, {0.95, 
                  1.15}}\)*) 
 (*\(,\)\(Frame -> 
              True\)*) 
 (*\(,\)\(FrameLabel -> {\*""(W^2)"", ""}\)\ *) 
, 
          AxesLabel -> {W\^2, ""}
, 
          PlotLabel -> \*""((Subscript[F, 2])^D)/((Subscript[F, 2])^N) at (Q^2)= "" <> ToString[
                qsq] <> "
Paris potential"
];\)\)

Clear[xgrid, qsq, q2];