Exact solution of the integral

Question

Can the following integral be solved exactly?

ClearAll["Global`*"]

Psi[r_, z_, b_, a_] = Exp[-b*z^2]*Exp[-a*r^2];

VB[r_, z_] = 1/(2*Sqrt[r^2 + z^2])*Exp[-3/10*Sqrt[r^2 + z^2]];

PB[b1_, a1_, b2_, a2_] = 
 Integrate[
  Psi[r, z, b2, a2]*VB[r, z]*Psi[r, z, b1, a1]*r, {r, 
   0, \[Infinity]}, {z, -Infinity, Infinity}, 
  Assumptions -> {a1 > 0, b1 > 0, a2 > 0, b2 > 0}]

Answer 1

The integration over r can be done, but not the integration over z. You can see this by splitting the integration

ClearAll["Global`*"]
Psi[r_, z_, b_, a_] = Exp[-b*z^2]*Exp[-a*r^2];
VB[r_, z_] = 1/(2*Sqrt[r^2 + z^2])*Exp[-3/10*Sqrt[r^2 + z^2]];
integrand = Psi[r, z, b2, a2]*VB[r, z]*Psi[r, z, b1, a1]*r

And now

int1 = Integrate[integrand, {r, 0, ∞}, 
  GenerateConditions -> False, 
  Assumptions -> {a1 > 0, b1 > 0, a2 > 0, b2 > 0}]

However, it can't do the above over z, either definite or indefinite

 Integrate[int1,z,Assumptions->{a1>0,b1>0,a2>0,b2>0},GenerateConditions->False]

Integrate[int1,{z,-Infinity,Infinity},Assumptions->{a1>0,b1>0,a2>0,b2>0},GenerateConditions->False]

Same result. Also Rubi can't.

Switching the order over the integration does not help. i.e integration over z first produces no result. So the above is the best you can get for now for analytical result, at least using direct Mathematica Integrate command.

ps. use exact numbers with Integrate. Changed your .3 to 3/10

V 13.2.1

Answer 2

If the special case of $a_1+a_2=b_1+b_2=\text{some positive constant }n$, then

integrand = (Psi[r, z, b2, a2]*VB[r, z]*Psi[r, z, b1, a1]*r // FullSimplify) /. a1 + a2 -> n /. b1 + b2 -> n
2 Integrate[integrand, {r, 0, ∞}, {z, 0, ∞}, Assumptions -> n > 0]

results in

(* 1/(2 n) - (3 E^((9/400)/n) Sqrt[π] Erfc[3/(20 Sqrt[n])])/(40 n^(3/2)) *)

There are other special cases such as $a_1+a_2=0$ and $b_1+b_2=n$ with $n>0$:

2 Integrate[(E^(-n z^2 - 3/10 Sqrt[r^2 + z^2]) r)/(2 Sqrt[r^2 + z^2]), {r, 0, ∞}, {z, 0, ∞}, 
  Assumptions -> n > 0]

(* (5 E^((9/400)/n) Sqrt[π] Erfc[3/(20 Sqrt[n])])/(3 Sqrt[n]) *)

Answer 3

The basic integral of Exp and Erf cannot be easily obtained . Its result contains the little known OwenT function (look up the Mathematica documentation for this) :

Integrate[
Erf[a*x + b]/E^(p^2*x^2), {x,0,Infinity}] = 
(Sqrt[Pi]/p)*((1/2)*Erf[(b*p)/Sqrt[a^2 + p^2]] + 
2*OwenT[(Sqrt[2]*b*p)/Sqrt[a^2 + p^2], a/p]).

From here you get contracting the a's and b's together with using elementary algebra to the result of Nasser' s integral with Exp and Erfc as :

(E^(9/400/a1)*Pi*Erfc[(3*Sqrt[-a1 + b1])/(20*Sqrt[a1]*Sqrt[b1])])/(4*
Sqrt[a1]*Sqrt[-a1 + b1])

Numeric check for b1 > a1 :

a1 = 0.3; b1 = 0.8; {NIntegrate[(Exp[9/(400*a1) + (a1 - b1)*z^2]*
Sqrt[Pi]*Erfc[(3 + 20*a1*z)/(20*Sqrt[a1])])/(4*
Sqrt[a1]), {z, -Infinity, Infinity}], 
(E^(9/400/a1)*Pi*Erfc[(3*Sqrt[-a1 + b1])/(20*Sqrt[a1]*Sqrt[b1])])/(4*
Sqrt[a1]*Sqrt[-a1 + b1])}
 (*  {1.66006, 1.66006}  *)

a nice simple form.

Edit:

an easier way to calculate

int = Integrate[E^(-a z^2)  Erfc[b + c z], {z, -Infinity, Infinity}]

(which resembles the OP's integral after integrating over r) is by differentiation under the integral sign :

Derivative with respect to b :

D[E^(-a z^2)  Erfc[b + c z], b]
(* -((2*E^((-a)*z^2 - (b + c*z)^2))/Sqrt[Pi]) *)

Now integrate over z:

Integrate[-((2*E^((-a)*z^2 - (b + c*z)^2))/Sqrt[Pi]),
{z, -Infinity, Infinity}, 
GenerateConditions -> False]
(* -(2/(E^((a*b^2)/(a + c^2))*Sqrt[a + c^2])) *)

Find the Antiderivative with respect to b to 'cancel' the derivative :

Integrate[-(2/(E^((a*b^2)/(a + c^2))*Sqrt[a + c^2])), b]
(* -((Sqrt[Pi]*Erf[(Sqrt[a]*b)/Sqrt[a + c^2]])/Sqrt[a]) *)

To get the missing constant of integration compare both expressions at b = 0 :

Integrate[Erfc[c*z]/E^(a*z^2), {z, -Infinity, Infinity}, 
GenerateConditions -> False]
(* Sqrt[Pi]/Sqrt[a] *)

and the final result is

int = (Sqrt[Pi]/Sqrt[a])*Erfc[(Sqrt[a]*b)/Sqrt[a + c^2]]

Check by plot:

a = 0.7; c = 0.3; Plot[{NIntegrate[
Erfc[b + c*z]/E^(a*z^2), {z, -Infinity, Infinity}], 
  (Sqrt[Pi]*Erfc[(Sqrt[a]*b)/Sqrt[a + c^2]])/Sqrt[a]}, {b, 0, 1}, 
PlotStyle -> {Blue, Dashed}]

Answer 4

This looks suspiciously like a time-independent quantum mechanics problem, solving for the expectation value of VB. Forgive me, if I am mistaken. If it is, in fact, a quantum mechanics problem, you are approaching the problem incorrectly, and the second psi should be the complex conjugate transpose of the first psi. Further, you are only integrating over two dimensions, which would yield a nonsensical result. Perhaps converting from cylindrical to spherical coordinates would help.