I have a Gaussian kernel that I wish to evaluate $$\int_{0}^{\infty} \int_{0}^{2\pi} \exp{[-r^2 - r_\alpha^2 +2 r r_\alpha \cos{(\theta-\theta_\alpha)}]} [r -r_\alpha \cos{(\theta-\theta_\alpha)}]f(r,\theta) ~ \mathrm{d}\theta ~\mathrm{d}r $$
We are working in polar coordinates in 2D, hence, $r_i=(r\cos\theta,r\sin\theta)$. The functions within the kernel $f(r,\theta)$ actually come from a tensor which looks like $[r_ir_jr_k....]/r^2$, starting from the second tensor, $f(r,\theta)=[\cos^2\theta,\cos\theta\sin\theta,\sin^2\theta]$. Then they could be $f(r,\theta)=[r\cos^3\theta,r\cos^2\theta\sin\theta,r\sin^2\theta\sin\theta,r\sin^3\theta]$ for the next higher order terms and so on.
For the first few terms $f(r,\theta)=[\cos^2\theta,\cos\theta\sin\theta,\sin^2\theta]$, I am able to get Mathematica to give me an output. For e.g., $$\int_{0}^{\infty} \int_{0}^{2\pi} \exp{[-r^2 - r_\alpha^2 +2 r r_\alpha \cos{(\theta-\theta_\alpha)}]} [r -r_\alpha \cos{(\theta-\theta_\alpha)}][\cos^2\theta,\cos\theta\sin\theta,\sin^2\theta] ~ \mathrm{d}\theta ~\mathrm{d}r = \left[\frac{\pi}{2},0,\frac{\pi}{2}\right]. $$ However, it struggles for the next higher order terms like $(r\cos^3 \theta)$ etc. Hence, I am wondering if I can see how Mathematica is doing this integral and can do it on paper itself. Also, if someone knows how to do this by hand, I will be grateful to know how this may be done. Thanks for any help, much appreciated.
Mathematica Code: Integrate[(Exp[-r^2 + 2*r*a*Cos[\[Theta] - b]])*(r - a*Cos[\[Theta] - b])*(Cos[\[Theta]])^2, {r, 0, \[Infinity]} , {\[Theta], , 2*Pi}, Assumptions -> a \[Element] PositiveReals && b \[Element] Reals]
(Trick changing $\theta = \theta + b$ makes Mathematica run faster)
Integrate[(Exp[-r^2 + 2*r*a*Cos[\[Theta]]])*(r-a*Cos[\[Theta]])*(Cos[\[Theta] + b])^2, {r, 0, \[Infinity]} , {\[Theta], b, b+2*Pi}, Assumptions -> a \[Element] PositiveReals && b \[Element] Reals]