Improve speed of NIntegrate for functions which are sums of complex exponentials

Question

I need to improve the efficiency for numerically evaluating a matrix element, which uses NIntegrate of the several functions (which are all sums of complex exponentials).

Creating the functions from eigenvectors of some matrix Hamiltonian:

Clear[eEigens, eEnergy, eFunc, hEigens, hEnergy, hFunc]

eEigens[β_, γ_, dtoR_, f_] := 
 eEigens[β, γ, dtoR, f] = Block[{},
   m = 10;
   ham = Table[i^2 KroneckerDelta[i, j], {i, -m, m}, {j, -m, m}]
     + Table[(-β*dtoR*(1 - γ) - I f)*
       KroneckerDelta[i, j + 1], {i, -m, m}, {j, -m, m}]
     + Table[(-β*dtoR*(1 - γ) + I f)*
       KroneckerDelta[i, j - 1], {i, -m, m}, {j, -m, m}]
     + Table[-(1/4) β (1 + γ)*
       KroneckerDelta[i, j + 2], {i, -m, m}, {j, -m, m}]
     + Table[-(1/4) β (1 + γ)*
       KroneckerDelta[i, j - 2], {i, -m, m}, {j, -m, m}];

 eEigs = 
    SortBy[Transpose[Eigensystem[ham, Method -> "Banded"]], First] 
   ] (* Diagonalize hamltonian and store results *)

eEnergy[β_, γ_, dtoR_, f_, n_] := 
  eEnergy[β, γ, dtoR, f, n] = 
  eEigens[β, γ, dtoR, f][[n + 1]][[1]] 

eFunc[β_, γ_, dtoR_, f_, n_] := 
  eFunc[β, γ, dtoR, f, n] = 
   Chop[1./Sqrt[N[2 π]]*
      eEigens[β, γ, dtoR, f][[n + 1]][[2]]].Table[
     Exp[I*j*φ], {j, -m, m}];

hEigens[μ_, β_, γ_, dtoR_, f_] := 
 hEigens[μ, β, γ, dtoR, f] = Block[{},
   m = 10;
   ham = Table[μ*i^2 KroneckerDelta[i, j], {i, -m, m}, {j, -m, 
       m}]
     + Table[(β*dtoR*(1 - γ) - I f)*
       KroneckerDelta[i, j + 1], {i, -m, m}, {j, -m, m}]
     + Table[(β*dtoR*(1 - γ) + I f)*
       KroneckerDelta[i, j - 1], {i, -m, m}, {j, -m, m}]
     + Table[
      1/4 β (1 + γ)*KroneckerDelta[i, j + 2], {i, -m, 
       m}, {j, -m, m}]
     + Table[
      1/4 β (1 + γ)*KroneckerDelta[i, j - 2], {i, -m, 
       m}, {j, -m, m}];
   hEigs = 
    SortBy[Transpose[Eigensystem[ham, Method -> "Banded"]], First] 
   ] (* Diagonalize  hamltonian *)

hEnergy[μ_, β_, γ_, dtoR_, f_, n_] := 
 hEnergy[μ, β, γ, dtoR, f, n] = 
  hEigens[μ, β, γ, dtoR, f][[n + 1]][[
   1]] (* Read stored results to give eigenvalues *)

hFunc[μ_, β_, γ_, dtoR_, f_, n_] := 
  hFunc[μ, β, γ, dtoR, f, n] = 
   Chop[1./Sqrt[N[2 π]]*
      hEigens[μ, β, γ, dtoR, f][[n + 1]][[
       2]]].Table[Exp[I*j*φ], {j, -m, m}];

So eFunc and hFunc are sums of complex exponentials for each eEigs and hEigs respectively, with coefficients given by the corresponding eigenvectors eEigens and hEigens.

I then create the following function, which will then be evaluated 100's of times to populate a matrix.

integral[μ_?NumericQ, β_?NumericQ, γ_?NumericQ, 
  dtoR_?NumericQ, f_?NumericQ, n0_, M0_, n_, M_] := NIntegrate[
    eFunc[β, γ, dtoR, f, n]*
     hFunc[μ, β, γ, dtoR, f, M]*
     Conjugate[eFunc[β, γ, dtoR, f, n0]]*
     Conjugate[hFunc[μ, β, γ, dtoR, f, M0]
      ], {φ, -π, π}, AccuracyGoal -> 10,
    Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule", "SymbolicProcessing" -> 0}
    ] 

This takes about 0.02 seconds on my laptop:

integral[0.4, 2., 1.02, 20., 0., 1, 1, 2, 2] // AbsoluteTiming
(* {0.020465, -0.00181996} *)

Is there a quicker way of doing this?

I then create a final matrix for example:

finalmatrix[μ_, β_, γ_, dtoR_, f_, v0_, x_, max_] :=
   Rationalize[
   Block[{},
    g = Table[(-N[2 π]*v0)/(
       eEnergy[β, γ, dtoR, f, jj] + 
        hEnergy[μ, β, γ, dtoR, f, kk] - x)*
       integral[μ, β, γ, dtoR, f, ii, ff, jj, kk], 
{ii, 0, max}, {ff, 0, max}, {jj, 0, max}, {kk, 0, max}];
    h = ArrayReshape[g, {(max + 1)^2, (max + 1)^2}] - 
      IdentityMatrix[(max + 1)^2]
    ]
   , 0];

which takes far too long for the purposes of subsequent calculations

finalmatrix[0.4, 2., 1.02, 20., 0., -0.2, x, 9]; // AbsoluteTiming
(* {103.234, Null} *)

Any help will be hugely appreciated!