# Numerical Integral and Partial Differentiation Of a Complicated Expression

I need help plotting Figures 1 and 2 solving EQ 36 and 37 from https://arxiv.org/pdf/2310.12445. here I am trying to find \Gamma,\phi and their Q values for different a_B values. But the interpolation and derivative are somehow going wrong!

Figure 1(a) matches somehow using this code

n0 = 10^20; la = 45*10^-9; a1a2 = 2.9*10^-9; mb = 14.45*10^-26; ma =
3.82*10^-26; mab =
ma*mb/(ma + mb); p = (2 n0)/\[Pi]^2 (4 \[Pi]^2 (a1a2)^2)/
mab^2;(*Subscript[a, B]=5.3*10^-9*)
(*Ek=k^2/(2mb);gb=(4 \[Pi] aB)/mb;\[Omega]k=Sqrt[Ek^2+2 \
n0*gb*Ek]=Sqrt[k^4/(4 mb^2)+2*n0*(4 \[Pi] aB)/mbk^2/(2 mb)];*)
\[CapitalGamma][t_?NumericQ, aB_] :=
NIntegrate[
p*k^2*Exp[-0.5*la^2*k^2]*k^2/(
2 mb) (1 -
Cos[Sqrt[k^4/(4 mb^2) + 2*n0*(4 \[Pi] aB)/mb k^2/(2 mb)]
t])/(Sqrt[
k^4/(4 mb^2) + 2*n0*(4 \[Pi] aB)/mb k^2/(2 mb)])^3, {k,
0, \[Infinity]}, MaxRecursion -> 44, PrecisionGoal -> 6] // Quiet;
data = ParallelTable[{{t, \[CapitalGamma][t,
0.5*5.3*10^-9]}, {t, \[CapitalGamma][t,
5.3*10^-9]}, {t, \[CapitalGamma][t, 2*5.3*10^-9]}}, {t,
0*10^-3, 2*10^-3, 0.1*10^-3}] // Quiet;
ListLinePlot[{data[[All, 1]], data[[All, 2]], data[[All, 3]]},
PlotRange -> {{0, 1*10^-3}, {0, 0.2}}, InterpolationOrder -> 3,
PlotLegends -> {"0.5\!$$\*SubscriptBox[\(a$$, $$B$$]\)",
"\!$$\*SubscriptBox[\(a$$, $$B$$]\)",
"2\!$$\*SubscriptBox[\(a$$, $$B$$]\)"}, Joined -> True,
PlotStyle -> {{Blue, Thickness[0.005]}, {Orange,
Thickness[0.005]}, {Green, Thickness[0.005]}}, Frame -> True,
FrameLabel -> {{Style["\[CapitalGamma]",
FontFamily -> "Times New Roman", FontSlant -> Italic,
FontWeight -> Bold, FontSize -> 10],
None}, {Style["t", FontFamily -> "Times New Roman",
FontSlant -> Italic, FontWeight -> Bold, FontSize -> 10], None}}]


But after interpolation and derivative, to find Q, does not match

n0 = 10^20; la = 45*10^-9; a1a2 = 2.9*10^-9; mb = 14.45*10^-26; ma =
3.82*10^-26; mab =
ma*mb/(ma + mb); p = (2 n0)/\[Pi]^2 (4 \[Pi]^2 (a1a2)^2)/mab^2;
\[CapitalGamma][t_?NumericQ, aB_] :=
NIntegrate[
p*k^2*Exp[-0.5*la^2*k^2]*k^2/(
2 mb) (1 -
Cos[Sqrt[k^4/(4 mb^2) + 2*n0*(4 \[Pi] aB)/mb k^2/(2 mb)]
t])/(Sqrt[
k^4/(4 mb^2) + 2*n0*(4 \[Pi] aB)/mb k^2/(2 mb)])^3, {k,
0, \[Infinity]}, MaxRecursion -> 44, PrecisionGoal -> 6] // Quiet;
data = Flatten[
ParallelTable[{t, aB, \[CapitalGamma][t, aB]}, {t, 0.1*10^-3,
2*10^-3, 0.1*10^-3}, {aB, 0.5*5.3*10^-9, 2*5.3*10^-9,
0.1*5.3*10^-9}], 1] // Quiet;
gt = Interpolation[data, InterpolationOrder -> 3] // Quiet;
qnd = ParallelTable[{t,
aB, (aB*Derivative[0, 1][gt][t, aB])^2/(Exp[
2*\[CapitalGamma][t, aB]] - 1)}, {t, 0.1*10^-3, 2*10^-3,
0.1*10^-3}, {aB, 0.5*5.3*10^-9, 2*5.3*10^-9, 0.1*5.3*10^-9}] //
Quiet;
set1 = {#[[1]], #[[3]]} & /@ qnd[[All, 1]];(*Taking data {t,0.5aB}*)
set2 = {#[[1]], #[[3]]} & /@ qnd[[All, 6]];(*Taking data {t,aB}*)
set3 = {#[[1]], #[[3]]} & /@ qnd[[All, 16]];(*Taking data {t,2aB}*)
ListLinePlot[{set1, set2, set3},
PlotRange -> {{0.1*10^-3, 2*10^-3}, {0, 0.1}},
InterpolationOrder -> 3,
PlotLegends -> {"0.5\!$$\*SubscriptBox[\(a$$, $$B$$]\)",
"\!$$\*SubscriptBox[\(a$$, $$B$$]\)",
"2\!$$\*SubscriptBox[\(a$$, $$B$$]\)"}, Frame -> True,
ImageSize -> Medium,
FrameLabel -> {{Style[
"\!$$\*SubscriptBox[SuperscriptBox[\(Q$$, $$\[DoubleVerticalBar]\$$], SubscriptBox[$$a$$, $$B$$]]\)", FontFamily -> "Times New Roman",
FontSlant -> Italic, FontWeight -> Bold, FontSize -> 10],
None}, {Style["t", FontFamily -> "Times New Roman",
FontSlant -> Italic, FontWeight -> Bold, FontSize -> 10], None}}]


On the other hand for \phi integral does not match with paper

It is coming huge for me in 10^34 order . Here is the code for \phi

n0 = 10^20; la = 45*10^-9; a1a2 = 2.9*10^-9; mb = 14.45*10^-26; ma =
3.82*10^-26; mab =
ma*mb/(ma + mb); p = (2 n0)/\[Pi]^2 (4 \[Pi]^2 (a1a2)^2)/mab^2;
\[Phi][t_?NumericQ, aB_] :=
NIntegrate[
p*k^2*Exp[-0.5*la^2*k^2]*k^2/(
2 mb) (Sin[
Sqrt[k^4/(4 mb^2) + 2*n0*(4 \[Pi] aB)/mb k^2/(2 mb)] t] -
t Sqrt[k^4/(4 mb^2) + 2*n0*(4 \[Pi] aB)/mb k^2/(2 mb)])/(Sqrt[
k^4/(4 mb^2) + 2*n0*(4 \[Pi] aB)/mb k^2/(2 mb)])^3, {k,
0, \[Infinity]}, MaxRecursion -> 44, PrecisionGoal -> 6] // Quiet;
data = ParallelTable[{{t, \[Phi][t, 0.5*5.3*10^-9]}, {t, \[Phi][t,
5.3*10^-9]}, {t, \[Phi][t, 2*5.3*10^-9]}}, {t, 0*10^-3,
2*10^-3, 0.1*10^-3}] // Quiet;

ListLinePlot[{data[[All, 1]], data[[All, 2]], data[[All, 3]]},
PlotRange -> All, InterpolationOrder -> 3,
PlotLegends -> {"0.5\!$$\*SubscriptBox[\(a$$, $$B$$]\)",
"\!$$\*SubscriptBox[\(a$$, $$B$$]\)",
"2\!$$\*SubscriptBox[\(a$$, $$B$$]\)"}, Joined -> True,
PlotStyle -> {{Blue, Thickness[0.005]}, {Orange,
Thickness[0.005]}, {Green, Thickness[0.005]}}, Frame -> True,
FrameLabel -> {{Style["\[CapitalPhi]",
FontFamily -> "Times New Roman", FontSlant -> Italic,
FontWeight -> Bold, FontSize -> 10],
None}, {Style["t", FontFamily -> "Times New Roman",
FontSlant -> Italic, FontWeight -> Bold, FontSize -> 10], None}}]


Q for \phi is uncalculatable

n0 = 10^20; la = 45*10^-9; a1a2 = 2.9*10^-9; mb = 14.45*10^-26; ma =
3.82*10^-26; mab =
ma*mb/(ma + mb); p = (2 n0)/\[Pi]^2 (4 \[Pi]^2 (a1a2)^2)/mab^2;
\[Phi][t_?NumericQ, aB_] :=
NIntegrate[
p*k^2*Exp[-0.5*la^2*k^2]*k^2/(
2 mb) (Sin[
Sqrt[k^4/(4 mb^2) + 2*n0*(4 \[Pi] aB)/mb k^2/(2 mb)] t] -
t Sqrt[k^4/(4 mb^2) + 2*n0*(4 \[Pi] aB)/mb k^2/(2 mb)])/(Sqrt[
k^4/(4 mb^2) + 2*n0*(4 \[Pi] aB)/mb k^2/(2 mb)])^3, {k,
0, \[Infinity]}, MaxRecursion -> 44, PrecisionGoal -> 6] // Quiet;
data = Flatten[
ParallelTable[{t, aB, \[Phi][t, aB]}, {t, 0.1*10^-3, 2*10^-3,
0.3*10^-3}, {aB, 0.5*5.3*10^-9, 2*5.3*10^-9, 0.1*5.3*10^-9}],
1] // Quiet;
gt = Interpolation[data, InterpolationOrder -> 3] // Quiet;
qnd = ParallelTable[{t, (5.3*10^-9*
Derivative[0, 1][gt][t, 5.3*10^-9])^2*
Exp[-2*\[Phi][t, 5.3*10^-9]]}, {t, 0.1*10^-3, 2*10^-3,
0.3*10^-3}] // Quiet;
ListLinePlot[qnd, PlotRange -> {{0, 2*10^-3}, {0, 2}},
InterpolationOrder -> 3]

• Why not to use nondimensional form of this equations? Commented Jul 9 at 4:04
• Hello@AlexTrounev! It's been a year. I hope you are well. Actually I tried and posted somehow the dimensionless analogous integration here: mathematica.stackexchange.com/questions/304814/…. But the answare does not match either with the paper. Commented Jul 9 at 5:08
• Here@AlexTrounev from the first program after integration the produced graph "\Gamma" seems close to paper. But in the second program for Q, when after interpolation and derivative I am getting an oscillatory graph rather than a smooth one. Commented Jul 9 at 5:16
• We need here some analytical job to compute expressions around t=0. Commented 2 days ago