Consider the following function dsPrimakovdtnucleus:

Amo = 96;
sval = 0.9;
fmtoGeVm1 = 5;
R1 = Sqrt[(1.23*Amo^(1/3) - 0.6)^2 + 7/3 Pi^2*0.52^2 - 5*sval^2];
FHelm[Z_, q2_] = 3*Z*BesselJ[1, Sqrt[q2]*R1*fmtoGeVm1]/(Sqrt[q2]*R1*fmtoGeVm1)*Exp[-((q2*sval^2*fmtoGeVm1^2)/2)];
aEMval = 1/137;
(*Differential cross section*)
sinv[Eg_, mt_] = 2*mt*Eg + mt^2;
dsPrimakovdtnucleus[ma_, t_, Eg_] = (2*aEMval*FHelm[Zmo, -t]^2*mZ^4/(t^2*(mZ^2 - s)^2 (t - 4*mZ^2)^2) (ma^2*t (mZ^2 + s) -ma^4 mZ^2 - t*((mZ^2 - s)^2 + s*t)) /. {s -> sinv[Eg, mZ]} // Simplify) /. {mZ -> mMo};

Here, t is a variable, while ma, Eg are positive parameters, with Eg > Egmin, where

Egmin[ma_] = (ma^2 + 2*mMo*ma)/(2*mMo);

I need to integrate dsPrimakovdtnucleus over t which ranges within tmin and tmax:

tmin[ma_, Eg_] = 
       mt]) - ((sinv[Eg, mt] - mt^2)/(2*Sqrt[sinv[Eg, mt]]) + 
      Sqrt[(sinv[Eg, mt] + ma^2 - mt^2)^2/(4*sinv[Eg, mt]) - ma^2])^2/.{mt->mMo};
tmaxTemp[ma_, Eg_, mt_] = 
       mt]) - ((sinv[Eg, mt] - mt^2)/(2*Sqrt[sinv[Eg, mt]]) - 
      Sqrt[(sinv[Eg, mt] + ma^2 - mt^2)^2/(4*sinv[Eg, mt]) - ma^2])^2;
tmaxLimit[ma_, Eg_] = 
  Assuming[Eg > 0 && mt > 0, 
   Simplify[Normal@Series[tmaxTemp[ma, Eg, mt], {ma, 0, 4}]]];
tmax[ma_, Eg_] = 
  If[ma/Eg > 0.01, Evaluate[tmaxTemp[ma, Eg, mt]/.{mt->mMo}], 
   Evaluate[tmaxLimit[ma, Eg]]];

Note that I have expanded tmaxTemp in the limit ma/Eg << 1.

The integral is

sPrimakovNucleus[ma_, Eg_] := 
 If[Eg > Egmin[ma], 
   dsPrimakovdtnucleus[ma, t, Eg], {t, 
    tmin[ma, Eg], tmax[ma, Eg]}], 0]

It behaves well until the ratio ma/Eg becomes small. Say, sPrimakovNucleus[1, 10] returns 0.73 without any issues, while sPrimakovNucleus[1, 50] returns

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in t near {t} = {-0.00010003}. NIntegrate obtained 0.129474 and 0.122071 for the integral and error estimates.

The reason is that the integrand includes a pole, and the upper integration limit becomes very close to it, being regularized by the limit tmaxLimit. If it is very small (i.e., ma/Eg << 1), then, probably, Mathematica fails.

Could you please tell me how to avoid this problem?

Edit: added missing definitions.

Update: the problem was that once ma/Eg becomes much less than one, only a very narrow integration domain contributes to the total value. So I may manually reduce tmax[ma, Eg] to say 0.1 tmax[ma,Eg]. It would be better however to automatize this procedure...

  • $\begingroup$ Do we need definitions for mMo and Zmo to be able to test this to see if a fix works? $\endgroup$
    – Bill
    Commented Feb 2, 2023 at 2:04
  • $\begingroup$ @Bill : oh, I am totally sorry. Just added them. $\endgroup$ Commented Feb 2, 2023 at 7:06
  • 1
    $\begingroup$ Change all decimal numbers to exact rationals to allow precision greater than MachinePrecision. Then add two options to NIntegrate, WorkingPrecision->64 and MaxRecursion->32 That allowed me to evaluate sPrimakovNucleus[1, 50] without warning or error messages. Experiment with those two option values to see what you really need. $\endgroup$
    – Bill
    Commented Feb 3, 2023 at 7:05

1 Answer 1


Your integrand has an exponential-decay factor, which causes "a very narrow integration domain [to contribute] to the total value" when ma/Eg is small. A simple change of variables $t=-e^u$ fixes it.

sPrimakovNucleus[ma_, Eg_] := 
 If[Eg > Egmin[ma], 
   -Exp[u] dsPrimakovdtnucleus[ma, -Exp[u], Eg],
   {u, Log[-tmin[ma, Eg]], Log[-tmax[ma, Eg]]}

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.