Numerical integration converging too slowly with Tanhyberbolic function

Question

I am trying to make the plots using the second equation from the  appendix B.5 from the paper (https://arxiv.org/pdf/cond-mat/0008249.pdf) which is basically integration of polarization function.

My code is given as below:

    In[150]:= d = 1;
    t = 10^-5 d;
    avn = 0.875;
    q0 = 0.5;
    v = 0.5 d;
    ufc = 0.1 d;
    \[Epsilon]f0 = -0.736;
    \[Epsilon]f = 0.0985;
    \[Mu] = 0.0962;
    a = 0.0846;
    avnf = 0.5 - a^2;
    avnc = avn - avnf;
    \[Nu] = 10^-7;
    m = (3 \[Pi]^2/Sqrt[2])^(2/3) 1/(2 d);
    kval[\[Epsilon]_] := Sqrt[2 m (\[Epsilon] + d)];
    ds[\[Epsilon]_] := 3/(4 Sqrt[2] d) Sqrt[(\[Epsilon] + d)/d];
    
    bepsk[\[Epsilon]_] := \[Epsilon] + ufc avnf;
    bepskq [\[Epsilon]_, q_, \[Theta]_] := 
      bepsk[\[Epsilon]] + kval[\[Epsilon]] q Cos[\[Theta]]/m + q^2/(2 m);
    rk[\[Epsilon]_] := 
      Sqrt[(\[Epsilon]f - bepsk[\[Epsilon]])^2 + (2 v a)^2];
    rkq[\[Epsilon]_, q_, \[Theta]_] := 
      Sqrt[(\[Epsilon]f - bepskq[\[Epsilon], q, \[Theta]])^2 + (2 v a)^2];
    akq[\[Epsilon]_, q_, \[Theta]_] := v a/rkq[\[Epsilon], q, \[Theta]];
    bkq[\[Epsilon]_, q_, \[Theta]_] := -v a/rkq[\[Epsilon], q, \[Theta]];
    ak[\[Epsilon]_] := v a/rk[\[Epsilon]];
    bk[\[Epsilon]_] := -v a/rk[\[Epsilon]];
    uk2 [\[Epsilon]_] := 
      0.5 (1 - (\[Epsilon]f - bepsk[\[Epsilon]])/rk[\[Epsilon]]);
    ukq2 [\[Epsilon]_, q_, \[Theta]_] := 
      0.5 (1 - (\[Epsilon]f - bepskq [\[Epsilon], q, \[Theta]])/
         rkq[\[Epsilon], q, \[Theta]]);
    vk2 [\[Epsilon]_] := 
      0.5 (1 + (\[Epsilon]f - bepsk[\[Epsilon]])/rk[\[Epsilon]]);
    vkq2 [\[Epsilon]_, q_, \[Theta]_] := 
      0.5 (1 + (\[Epsilon]f - bepskq [\[Epsilon], q, \[Theta]])/
         rkq[\[Epsilon], q, \[Theta]]);
    
    In[178]:= 
    E1k[\[Epsilon]_] := 
      0.5 (\[Epsilon]f + bepsk[\[Epsilon]] + rk[\[Epsilon]]);
    E1kq[\[Epsilon]_, q_, \[Theta]_] := 
      0.5 (\[Epsilon]f + bepskq[\[Epsilon], q, \[Theta]] + 
         rkq[\[Epsilon], q, \[Theta]]);
    E2k[\[Epsilon]_] := 
      0.5 (\[Epsilon]f + bepsk[\[Epsilon]] - rk[\[Epsilon]]);
    E2kq[\[Epsilon]_, q_, \[Theta]_] := 
      0.5 (\[Epsilon]f + bepskq[\[Epsilon], q, \[Theta]] - 
         rkq[\[Epsilon], q, \[Theta]]);
    nf1kq[\[Epsilon]_, q_, \[Theta]_] := 
      0.5 (1 - Tanh[((E1kq[\[Epsilon], q, \[Theta]] - \[Mu])/(2 t))]);
    nf2kq[\[Epsilon]_, q_, \[Theta]_] := 
      0.5 (1 - Tanh[((E2kq[\[Epsilon], q, \[Theta]] - \[Mu])/(2 t))]);
    nf1k[\[Epsilon]_] := 0.5 (1 - Tanh[((E1k[\[Epsilon]] - \[Mu])/(2 t))]);
    nf2k[\[Epsilon]_] := 
      0.50 (1 - Tanh[((E2k[\[Epsilon]] - \[Mu])/(2 t))]);
    
    In[186]:= 
    part1[\[Epsilon]_, q_, \[Theta]_] := (
      vkq2[\[Epsilon], 
        q, \[Theta]] vk2[\[Epsilon]] (nf1kq[\[Epsilon], q, \[Theta]] - 
         nf1k[\[Epsilon]]))/(
      E1k[\[Epsilon]] - E1kq[\[Epsilon], q, \[Theta]] + I \[Nu]);
    part2[\[Epsilon]_, q_, \[Theta]_] := (
     vkq2[\[Epsilon], 
       q, \[Theta]] uk2[\[Epsilon]] (nf1kq[\[Epsilon], q, \[Theta]] - 
        nf2k[\[Epsilon]]))/(
     E2k[\[Epsilon]] - E1kq[\[Epsilon], q, \[Theta]] + I \[Nu]); 
    part3[\[Epsilon]_, q_, \[Theta]_] := (
     ukq2[\[Epsilon], 
       q, \[Theta]] vk2[\[Epsilon]] (nf2kq[\[Epsilon], q, \[Theta]] - 
        nf1k[\[Epsilon]]))/(
     E1k[\[Epsilon]] - E2kq[\[Epsilon], q, \[Theta]] + I \[Nu]); 
    part4[\[Epsilon]_, q_, \[Theta]_] := (
     ukq2[\[Epsilon], 
       q, \[Theta]] uk2[\[Epsilon]] (nf2kq[\[Epsilon], q, \[Theta]] - 
        nf2k[\[Epsilon]]))/(
     E2k[\[Epsilon]] - E2kq[\[Epsilon], q, \[Theta]] + I \[Nu]);
    
    In[188]:= 
    Mff[\[Epsilon]_, q_, \[Theta]_] := 
      ds[\[Epsilon]] (part1[\[Epsilon], q, \[Theta]] + 
         part2[\[Epsilon], q, \[Theta]] + part3[\[Epsilon], q, \[Theta]] +
          part4[\[Epsilon], q, \[Theta]]);
    Pff[q_] := 
      NIntegrate[ -Mff[\[Epsilon], q, \[Theta]] Sin[\[Theta]], {\[Theta], 
        0, \[Pi]}, {\[Epsilon], -d, d}, MaxRecursion -> 60, 
       Method -> {GlobalAdaptive, MaxErrorIncreases -> 8000}, 
       WorkingPrecision -> 10];
    
    In[190]:= 
    qa = 0;
    qb = 2;
    Ns = 10;
    h = (qb - qa)/(Ns - 1);
    res = {};
    For[i = 0, i < Ns, i++,
     q = qa + h *i;
      AppendTo[res, q];
     res]
    
    data =  {};
    For[i = 1, i <= Ns, i++,
       AppendTo[data , {res[[i]] // N, Re[Pff[res[[i]]]]}]];
    
    
    
    In[198]:= ListLinePlot[
     Table[{data[[i, 1]], data[[i, 2]]}, {i, 1, Length[data]}], 
     Frame -> True, PlotStyle -> Blue, PlotRange -> Full, 
     FrameLabel -> {"q", "Polre"}]

)

I also followed https://mathematica.stackexchange.com/questions/85694/numerical-integration-converging-too-slowly but with no success.

My function involves Tan hyberbolic function and is the reason it is converging slowly.

Runnnig the code gave me an error

During evaluation of In[190]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. During evaluation of In[190]:= General::stop: Further output of NIntegrate::slwcon will be suppressed during this calculation.

Could someone help me out?