NIntegrate gives large value as result

Question

I have been trying to perform a Numerical Integration using the NIntegrate function. The problem I am facing is that on using the Integrand as a Bessel Function of Imaginary component, the result just fails to converge and the answer just blows up to some large number. None of the methods seem to give the correct result. The sample Integrand which I am trying to Integrate is:

fun[b1_, b2_, x1_, x2_] := 
 BesselK[0, Sqrt[-x1*b1]]*BesselI[0, Sqrt[-x2*b2]]

NIntegrate[
 fun[b1, b2, x1, x2], {b1, 0, Infinity}, {b2, 0, Infinity}, {x1, 0, 
  1}, {x2, 0, 1}, Method -> "AdaptiveMonteCarlo", AccuracyGoal -> 10]

    

During evaluation of In[113]:= 

Out[113]= -8.1181*10^8 - 7.19087*10^9 I    

---

This is the complete integrand I am trying to work on. Please do check for this integrand.

ha[b1_, b2_, x1_, 
      x2_] := ((b1 - b2 + 2)*
         BesselI[0, Sqrt[\[Beta]a[x1, x2, rJ\[Psi], rb, mB]]*b1]*
         BesselK[0, 
          Sqrt[\[Beta]a[x1, x2, rJ\[Psi], rb, mB]]*b2] + \[Theta][
          b2 - b1]
         BesselI[0, Sqrt[\[Beta]a[x1, x2, rJ\[Psi], rb, mB]]*b2]*
         BesselK[0, Sqrt[\[Beta]a[x1, x2, rJ\[Psi], rb, mB]]*b1])*
      BesselK[0, Sqrt[\[Alpha][x1, x2, rJ\[Psi], mB]*b1]]

The functions in the integrand are:

   \[Alpha][x1_, x2_, rJ\[Psi]_, 
      mB_] := -((x1 - x2)*(x1 - x2*rJ\[Psi]^2))*mB^2
    \[Beta]a[x1_, x2_, rJ\[Psi]_, rb_, 
      mB_] := -((1 - x2)*(1 - x2*rJ\[Psi]^2) - rb^2)*mB^2

mB, rJpsi,rb, etc are constant values.

Answer 1

The integral under consideration diverges. The one splits into the product of two double improper integrals. One of these is

Integrate[BesselK[0, Sqrt[-x1*b1]], {b1, 0, Infinity}, {x1, 0, 1}]

which performs

"Integral of -(2/b1)+(2\I\BesselK[1,I\Sqrt[b1]])/Sqrt[b1] does not converge on {0,[Infinity]}"

and

$\int_0^{\infty } \frac{-2+2 i \sqrt{\text{b1}} K_1\left(i \sqrt{\text{b1}}\right)}{\text{b1}} \, d\text{b1}$