Dear Mathematica community, in a study about a Bose fluid, the vacuum density of energy (using coordinates to reduce the problem in 1D) is:

emod[k_, \[Mu]_, g0_, g1_, \[Lambda]_] := 
 k^2 (k^2/2 (k^2/2 + 
      2 \[Mu] (g0 + (4 \[Pi] g1)/(k^2 + \[Lambda]^2))/(
       g0 + (4 \[Pi] g1)/\[Lambda]^2)))^(1/2)

and after a renormalization procedure, in order to absorb the UV divergence (k->infinity), the density becomes:

emodrin1[k_, \[Mu]_, g0_, g1_, \[Lambda]_] := 
 k^2 (k^2/2 (k^2/2 + 
       2 \[Mu] (g0 + (4 \[Pi] g1)/(k^2 + \[Lambda]^2))/(
        g0 + (4 \[Pi] g1)/\[Lambda]^2)))^(1/2) - 
  1/2 k^4 - (
   g0 \[Mu])/(g0 + (4 g1 \[Pi])/\[Lambda]^2)  k^2 - ((
    4 g1 \[Pi] \[Mu])/(g0 + (4 g1 \[Pi])/\[Lambda]^2) - (
    g0^2 \[Mu]^2)/(g0 + (4 g1 \[Pi])/\[Lambda]^2)^2)

My aim is to numerically integrate:

NIntegrate[emodrin1[k, 1, 3, -1, 1], {k, 0, \[Infinity]}]

the parameter are chosen in order to don't have problem with the square root ( any other good choice could be fine), but my version (Mathematica 10.0 Student Editon) is not able to solve the calculation, I also tried changing the integration method, the WorkingPrecision and the Maxrecursion but without any results. So is there a smarter way to solve the integral? In addition I tried to plot:

Plot[emodrin1[k, 1, 3, -1, 1], {k, 0, 1000}, 
PlotRange -> {-0.0003, 0.0003}]

And the result is strange because for high k there are a lot of oscillations that present a fractal behavior, if I zoom in a region of high k the oscillations have the same aspect. In addition the amplitude of oscillations is increasing with k. Because of the renormalization procedure ( I am quite confident it is right) I expect that the function is going to 0 for k to infinity enough fast to be integrated. May you help me with this problems? are they correlated? Thanks you for the attention and the answers.