Deuteron wave function.

The deuteron wave function is parameterized in a standard way as follows: 
                                          (Subscript[Ψ, m])(r) = (1/(Sqrt[4  π]))(((u(r))/r) + ((w(r))/r)(((Subscript[S, 12])((OverHat[r])))/((Sqrt[8]))))(Subscript[χ, 1  m]), 
where u(r) and w(r) are the S- and D-wave radial w.f. respectively, (Subscript[χ, 1  m]) is the spin 1 w.f. with  (Subscript[S, z]) = m,  \(\(S\_12\)(r\&^)\) = 3(Subscript[σ, 1])·\(r\&^\)(Subscript[σ, 2])·\(r\&^\)-(Subscript[σ, 1])·(Subscript[σ, 2])  with (OverHat[r])=(Overscript[r, →])/r is the tensor operator  constructed from the Pauli matrices and which acts as a projector onto the l = 2 state. The wave function normalization: \(TraditionalForm\`∫\ \(\(∣\)\(Ψ\)\( ∣ \^2\)\(\(d\^3\) r\)\) = \(∫\_0\%∞ 
              dr\ \((u\^2 + w\^2)\)\  = 1\)\).
In the momentum space the wave function reads as follows:
                                         (Subscript[Ψ, m])(p) = \(\(\@\(2  π\^2\)\) \((\(ψ\_0\)(
                  p)\  + \ \(\(ψ\_2\)(
                    p)\) \(\(S\_12\)(p\&^)\)\/\@8)\) χ\_\(1  m\)\), 
Normalization:  \(TraditionalForm\`∫\ \(\(∣\)\(Ψ\)\( ∣ \^2\)\(d\^3 
                  p\/\((2  π)\)\^3\)\) = \(∫\_0\%∞ 
              dp\ \(\(p\^2\)(ψ\_0\^2 + ψ\_2\^2)\) = 1\)\).

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

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

hc is Mev fm  conversion constant,
mN is the isoscalar nucleon mass in fm^-1

hc = 197.327053; 
mN =  (938.27231 + 939.56563)/2. /hc;

sty = {
{Dashing[{}], 
            Thickness[0.005]}
, {Dashing[{}], 
            Thickness[0.005], 
            RGBColor[1, 0, 0]} 
, {Dashing[{}], 
            Thickness[0.005], 
            RGBColor[0, 1, 0]}
};

\(\(\(\ \)\( (*\ \(SetDirectory["e:/user/ser/phys/deuteron/WF"];\)\ *) \)\(
\)\( (*\ 
      set\ path\ to\ package\ *) \)\(
\)\(\(dpath = "e:/ser/phys/deuteron/WF/";\)

    Get[\ dpath <> "DeuteronWF.m"\ ]\)\)\)

? setdeuteron psi ed

setdeuteron[1];

u[r_] = psi[1, 0, r];
w[r_] = psi[1, 2, r];
psi0[k_] := psi[0, 0, k];
psi2[k_] := psi[0, 2, k];

The wave function squared:
|Ψ(r)\( | \^2\)=(1/(4  π (r^2)))[
            (u^2) + (w^2) + (3  (cos^2) θ - 1) (u w/(Sqrt[2]) - 
                    (w^2)/4)]

psisq[r_, 
          y_] := (u[r]^2 + 
                w[r]^2 + 
(3*y^2 - 
                      1)*(u[r]*w[r]/Sqrt[2] - w[r]^2/4))/
            r^2/(4*Pi);

Tests of the w.f.

Check normalization

Print["Integration in coordinate space:"]
NIntegrate[(u[r]^2 + w[r]^2), {r, 0, 
        Infinity}, {PrecisionGoal  ->  4}]

    
Print["Integration in momentum space:"]
NIntegrate[ 
      q^2 (psi0[q]^2 + psi2[q]^2), {q, 0, 
        Infinity}, {PrecisionGoal  ->  4}]

NIntegrate[
        psisq[r, y]*r^2, {r, 0, Infinity}, {y, -1, 
          1}, {PrecisionGoal -> 3}]*2*Pi

Check the asymptotics at r -> 0

\(\(\(\(\(Series[
      u[r], {r, 0, 3}]\)\(\ \)\( (*\(//\)\(Chop\)*) \)\(
\)
    \)

    Series[w[r], {r, 0, 3}]\)\(\ \)\( (*\(//\)\(Chop\)*) \)\)\)

Static parameters of the deuteron

D-state probability, (Subscript[P, D]) = ∫dr (w^2)(
                r) = ∫dq (q^2) ((Subscript[ψ, 2])^2)(q).

\(\(\(Print[\*""\!\(P\_D\) = ""]

    NIntegrate[\ 
      q^2*psi2[q]^2, {q, 0, 
        Infinity}, {PrecisionGoal\  -> \ 4}]\)\(\ \)\)\)

SD-interference, (Subscript[P, SD]) = ∫dq (q^2) (Subscript[ψ, 0])(
              q) (Subscript[ψ, 2])(q).

\(\(\(Print[\*""\!\(P\_SD\) = ""]

    NIntegrate[\ 
      q^2*psi0[q]*\ psi2[q], {q, 0, 
        Infinity}, {PrecisionGoal\  -> \ 4}]\)\(\ \)\)\)

Average polarization ⟨(Subscript[σ, z])⟩=(1/2)⟨(Subscript[(Ψ^†)((Subscript[σ, 1])+(Subscript[σ, 2])), z])Ψ⟩ 
and average tensor polarization ⟨T⟩=(1/2)⟨(Ψ^†)(p)(((Subscript[σ, 1])·p+(Subscript[σ, 2])·p)(Subscript[p, z])/(p^2)-(1/3)(Subscript[((Subscript[σ, 1])+(Subscript[σ, 2])), z]))Ψ(p)⟩
⟨(Subscript[σ, z])⟩ = 
        1 - (3/2) (Subscript[P, D])
\(TraditionalForm\`⟨T⟩ = \(1\/3\) \((P\_D + \(\@2\) P\_SD)\)\)

\(\(Print[\*""⟨(Subscript[σ, z])⟩ = "", 

        1 - 3. /2. *
            NIntegrate[\ 
              q^2*psi2[q]^2, {q, 0, 
                Infinity}, {PrecisionGoal\  -> \ 
                  4}]\ 
];\)\)

\(\(Print["⟨T⟩ = ", 
\((NIntegrate[\ 
                q^2\ *psi2[q]^2, {q, 0, 
                  Infinity}, {PrecisionGoal\  -> \ 4}]\  + 
              Sqrt[2]*NIntegrate[\ 
                  q^2*psi0[q]*psi2[q], {q, 0, 
                    Infinity}, {PrecisionGoal\  -> \ 4}])\)/
          3. 
];\)\)

Quaqdrupole moment, (Subscript[Q, d]) = ∫dr (r^2) w((Sqrt[8]) u - w)/
              20.

\(\(Print[\*""\!\(Q\_d\) = "", 
	\ NIntegrate[\ 
            r^2\ w[r]\ \((Sqrt[8]\ u[r] - w[r])\), \ {r, 0, 
              Infinity}, {PrecisionGoal\  -> \ 4}]
/
          20];\)\)

Average number density in the deuteron:
(Subscript[ρ, d])=\(TraditionalForm\`2 \(∫\(d\^3\) r\)\  | Ψ(
          r)\( | \^4\)\)

rho1 = 
      2*NIntegrate[
          psisq[r, y]^2*r^2, {r, 0, Infinity}, {y, -1, 
            1}, {PrecisionGoal -> 3}]*2*Pi

The root-mean-square radius of the deuteron.

\(\(Print[\*""⟨(r^2)⟩^(1/2)  = "", 

        rms = Sqrt[\ 
              NIntegrate[
                r^2*\((u[r]^2 + w[r]^2)\), {r, 0, 
                  Infinity}, {PrecisionGoal\  -> \ 4}]\ ]/
            2\ 
];\)\)

Average density in terms of r.m.s.:  ρ = 6/(4  π(R^3)), where (R^2) = 
        5 ⟨(r^2)⟩/3

Rem = (5. *rms^2/3. )^0.5; 
rho2 = 6. /(4*Pi*Rem^3)
rho2/rho1

Another density parameter: \(TraditionalForm\`2 | Ψ(0)\( | \^2\)\)

Series[2*psisq[r, 0], {r, 0, 2}] // Chop

Average momentum squared ⟨(q^2)⟩and kinetic energy ⟨(q^2)/2  
          M⟩ for different waves

\(\(\(Print["Paris:"]

    Print[\*""⟨(q^2)⟩_S = "", 
\((q2s = 
          NIntegrate[\ 
            q^4*psi0[q]^2, {q, 0, Infinity}, {PrecisionGoal\  -> \ 4}]; 
        q2s)\), \*"" (fm^(-2)) ""]

    Print[\*""⟨(q^2)⟩_D = "", 
\((q2d = 
          NIntegrate[\ 
            q^4*psi2[q]^2, {q, 0, Infinity}, {PrecisionGoal\  -> \ 4}]; 
        q2d)\)\ , \*"" (fm^(-2)) ""]

    Print[\*""⟨(q^2)/2(Subscript[m, N])⟩_S = "", q2s/\((2\ mN)\)\ hc, "  MeV"]

    Print[\*""⟨(q^2)/2(Subscript[m, N])⟩_D = "", 
      q2d/\((2\ mN)\)\ hc, "  MeV"\ ]\ \ \ 

    Print["⟨T⟩ = ", \((q2d + 
              q2s)\)/\((2\ mN)\)\ hc, "  MeV"\ ]\)\(\ \ \ \)\)\)

Plots

\(\(Plot[{u[r], w[r]}, {r, 0.001, 
          3. }
 (*\(,\)\(Frame\  -> 
            True\)\ *) 
, 
        PlotLabel -> "S & D waves in coordinate space"
, PlotRange\  -> \ All\ 
\ , 
        PlotStyle\  -> \ sty\ 
];\)\)

\(\(Plot[{psi0[q], psi2[q]*20}
, {q, 0, 
          3. }
 (*\(,\)\(Frame\  -> 
            True\)\ *) 
, 
        PlotLabel -> "S & D w.f. in momentum space"
, 
        PlotRange\  -> \ All\ 
, 
        PlotStyle\  -> \ sty\ 
];\)\)