Discretization and Dirac Delta function

Question

Consider the following integral: \begin{align} I= \int_{\mathbb{R}^9}\frac{d^3\mathbf{x_1}\,d^3\mathbf{x_2}\,d^3\mathbf{x_3}}{|\mathbf{x}_1|\sqrt{|\mathbf{x}_2||\mathbf{x}_3|}}f(\mathbf{x_2})f(\mathbf{x_3})\delta(\Omega+|\mathbf{x}_1|-|\mathbf{x}_2|)\delta(\Omega+|\mathbf{x}_1|-|\mathbf{x}_3|) \end{align} where $\mathbf{x_j} = (x_j,y_j,z_j)$. This is a nine-dimensional integral. I would typically choose $f(\mathbf{x})$ to be some Gaussian of the form \begin{align} f(\mathbf{x}-\mathbf{x}_0) = \exp -\sum_{j=1}^3\frac{-(x_j-x_{0,j})^2}{2\sigma_j} \end{align} where $\mathbf{x}_{0}$ is some constant vector. Now, this integral is most likely not doable analytically, so I would need NIntegrate. NIntegrate does not understand DiracDelta[x] so the way I would try to do this is using $\text{ImplicitRegion}$$\left[\Omega+|\mathbf{x}_1|-|\mathbf{x}_2| = 0\, \wedge \, \Omega+|\mathbf{x}_1|-|\mathbf{x}_3| = 0\right]$.

Problem: it turns out that this does not work (i.e. Mathematica does not compute this), and the reason (as far as I understand) is this: at the backend, Mathematica needs to perform DiscretizeRegion which is only supported for three-dimensional domain (the error messages show this).

Question: is there a way out of this problem, i.e. perform this numerically without facing problem with DiscretizeRegion?

This is mainly frustrating because if $\sigma_j=\sigma$ for $j=1,2,3$, I could in fact do this integral by hand but for my problem at hand I need at least $\sigma_3\neq \sigma_1,\sigma_2$ which definitely makes it not analytically solvable.

Remark: I tried breaking the problem into a few steps, namely first doing the $\mathbf{x_2}$ integral (integral over $\mathbf{x}_3$ is identical) over region $R$ and then hoping that integrating over $\mathbf{x}_1$ will work. This means

R[y1_, y2_, y3_, s1_, s2_, s3_, gap_] := 
DiscretizeRegion[ImplicitRegion[{Sqrt[x1^2 + x2^2 + x3^2] - Sqrt[y1^2 + y2^2 + y3^2] == gap},
 {{x1, -5 s1, 5 s1}, {x2, -5 s2, 5 s2}, {x3, -5 s3, 5 s3}}]]

where I chose $s_j=\sigma_j$ just to make the integral domain bounded (ideally I want $x_j$ to be integrated over the real line since the original integral is over $\mathbb{R}^9$ instead of $\pm 5s_j$). Then

G[y1_, y2_, y3_, gap_, x0_, y0_, z0_, s1_, s2_, s3_] := NIntegrate[f[x1, x2, x3, x0, y0, z0, s1, s2, s3]/Sqrt[Sqrt[x1^2 + x2^2 + x3^2]], 
Element[{x1, x2, x3}, R[y1, y2, y3, s1, s2, s3, gap]]]

where $G$ essentially is the result after integrating over $\mathbf{x}_2$. The original integral is obtained by integrating $G(\mathbf{x_1})^2/|\mathbf{x_1}|$. Now, $G(\mathbf{x}_1)$ actually gives a definite result for fixed $\mathbf{x}_1:=(y_1,y_2,y_3)$ but I cannot integrate over $y_j$.

Answer 1

In spherical coordinates,

x1 = r1 {Sin[θ1] Cos[φ1], Sin[θ1] Sin[φ1], Cos[θ1]};
x2 = r2 {Sin[θ2] Cos[φ2], Sin[θ2] Sin[φ2], Cos[θ2]};
x3 = r3 {Sin[θ3] Cos[φ3], Sin[θ3] Sin[φ3], Cos[θ3]};

and with the definition (modified a bit from yours because I think you had some typos)

f[x_, y_, z_] = Exp[-((x - x0)^2/(2 σx^2))]*
                Exp[-((y - y0)^2/(2 σy^2))]*
                Exp[-((z - z0)^2/(2 σz^2))];

your integral becomes

Integrate[((f @@ x2) (f @@ x3))/(r1 Sqrt[r2 r3]) *
  r1 Sin[θ1] * r2 Sin[θ2] * r3 Sin[θ3] /. {r2 -> r1 + Ω, r3 -> r1 + Ω},
  {r1, Max[0,-Ω], ∞}, {θ1, 0, π}, {φ1, 0, 2 π},
  {θ2, 0, π}, {φ2, 0, 2 π},
  {θ3, 0, π}, {φ3, 0, 2 π}]

In this way you don't have any Dirac distributions left and you have a regular integral. Further, as $r_2=r_3\ge0$, the denominator is cancelled by the Jacobian and the whole integral becomes non-singular.

To continue, here's the explicit numerical integral, where $\theta_1$ and $\phi_1$ have been integrated exactly (thanks @yarchik!), and we are left with five integration variables, four of which are finite:

4π*NIntegrate[
  f[(r1+Ω)Sin[θ2]Cos[φ2], (r1+Ω)Sin[θ2]Sin[φ2], (r1+Ω)Cos[θ2]]*
  f[(r1+Ω)Sin[θ3]Cos[φ3], (r1+Ω)Sin[θ3]Sin[φ3], (r1+Ω)Cos[θ3]]*
  (r1+Ω)*Sin[θ2]*Sin[θ3],
  {r1, Max[0, -Ω], ∞}, {θ2, 0, π}, {φ2, 0, 2π}, {θ3, 0, π}, {φ3, 0, 2π}]

In this integral there is no cross-talk between $\vec{x}_2$ and $\vec{x}_3$, so we can decompose it further, noticing that the integrals over $\vec{x}_2$ and $\vec{x}_3$ give the same result $F(r_1+\Omega)$:

F[r_?NumericQ] := NIntegrate[f[r*Sin[θ]Cos[φ], r*Sin[θ]Sin[φ], r*Cos[θ]]*Sin[θ], {θ,0,π}, {φ,0,2π}]

4π*NIntegrate[(r1+Ω)*F[r1+Ω]^2, {r1, Max[0, -Ω], ∞}]

In this form, the 9-dimensional integral has been reduced to three dimensions.