NIntegrate gives large value as result

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.

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