(possibly indefinite) integral of product of error functions

As given in this question, the convolution of a uniform and Gaussian distribution gives

$$(*) \quad \quad p(t) = \frac{1}{(b-a)\sqrt{2 \pi \sigma^2}} \int_{a+\mu}^{b+\mu} \exp \left\{-\frac{(\tau -t)^2}{2\sigma^2} \right\} d\tau$$

which gives (I believe)

$$(**) \quad \quad p(t) = \frac{1}{2(b-a)} \left( \text{erf}(\frac{b+\mu-t}{\sigma\sqrt{2}}) - \text{erf}(\frac{a+\mu-t}{\sigma\sqrt{2}}) \right)$$

For simplicity's sake, let's set $\mu=0$. I'm interested in integrating (hopefully indefinite, but possibly definite) the product of two of these distributions, with different initial uniforms given by $\{a,b\}$ and $\{c,d\}$, and different noise levels, given by $\sigma_1$ and $\sigma_2$, over the same domain, $D$. So

$$(***) \int_{t \in D} \frac{1}{4(b-a)(d-c)} \left( \text{erf}(\frac{b-t}{\sigma_1\sqrt{2}}) - \text{erf}(\frac{a-t}{\sigma_1\sqrt{2}}) \right)\left( \text{erf}(\frac{d-t}{\sigma_2\sqrt{2}}) - \text{erf}(\frac{c-t}{\sigma_2\sqrt{2}}) \right)dt$$

is this solvable when $D=\mathbb{R}$, or is there even an exact antiderivative.

It's a little tricky since when multiplying it out and breaking it into parts, each one of the parts is a product of error functions and not integrable.

P.S. The motivation for this is a notion of a "soft overlap" between intervals/products of intervals. To see this, if the normalization of the uniform is taken away so that the uniform PDF becomes the indicator function of an interval, then this exact integral when evaluated with $\sigma_1=\sigma_2=0$ so that the Gaussian becomes a Dirac delta, should give the area of the overlap of the two intervals.

P.P.S. If this interval is indeed intractable, are their other noise models / CDFs besides Gaussian for which a similar integral can be evaluated?

EDIT: Possible Solution:

I have also put up a possible solution up as an answer and will select it unless someone points out an error.

Answer

I think this is a possible solution.

For the case where we want to integrate over all of $\mathbb{R}$ (and probably for other cases, but this is simplest), it looks like we can solve the integral by going back to the initial representation of the convolution of a uniform and a Gaussian, taking the product there, and exchanging orders of integration.

That is, solve

$$\frac{1}{(b-a)(d-c)2\pi\sigma_1\sigma_2} \int_{c}^{d}\int_{a}^{b}\int_{-\infty}^{\infty} \exp \left\{-\frac{(\tau -t)^2}{2\sigma_1^2} \right\}\exp \left\{-\frac{(\rho -t)^2}{2\sigma_2^2} \right\} dt d\tau d\rho$$

The innermost integral is just a cross-correlation (or convolution with one of the means negated) of two Gaussians, with solution

$$\begin{align*} &\frac{1}{2\pi \sigma_1\sigma_2}\int_{-\infty}^{\infty} \exp \left\{-\frac{(\tau -t)^2}{2\sigma_1^2} \right\}\exp \left\{-\frac{(\rho -t)^2}{2\sigma_2^2} \right\} dt\\ &=\frac{1}{\sqrt{2\pi(\sigma_1^2+\sigma_2^2)}}\exp\left(-\frac{(\tau-\rho)^2}{2(\sigma_1^2+\sigma_2^2)} \right) \end{align*}$$

So now we should be able to integrate this Gaussian function twice to get the desired quantity. The antiderivative of $\text{erf}$ is $z\text{erf}(z)+\exp(-z^2)/\sqrt{\pi} + C$, which will be needed on the final integration.

However, this looks like a pain, but Mathematica gives the solution to

$$\int _a^b\int _c^d\exp \left(-\frac{(\tau-\rho)^2}{2 \left(\sigma_1^2+\sigma_2^2\right)}\right) d\rho d\tau$$

as

$$\sqrt{\frac{\pi }{2}} \sqrt{\sigma_1^2+\sigma_2^2} \left(\sqrt{\frac{2}{\pi }} \sqrt{\sigma_1^2+\sigma_2^2} \left(e^{-\frac{(b-c)^2}{2 \left(\sigma_1^2+\sigma_2^2\right)}}-e^{-\frac{(a-c)^2}{2 \left(\sigma_1^2+\sigma_2^2\right)}}\right)+\sqrt{\frac{2}{\pi }} \sqrt{\sigma_1^2+\sigma_2^2} \left(e^{-\frac{(a-d)^2}{2 \left(\sigma_1^2+\sigma_2^2\right)}}-e^{-\frac{(b-d)^2}{2 \left(\sigma_1^2+\sigma_2^2\right)}}\right)+(a-c) \text{erf}\left(\frac{c-a}{\sqrt{2} \sqrt{\sigma_1^2+\sigma_2^2}}\right)+(d-a) \text{erf}\left(\frac{d-a}{\sqrt{2} \sqrt{\sigma_1^2+\sigma_2^2}}\right)-(b-c) \text{erf}\left(\frac{c-b}{\sqrt{2} \sqrt{\sigma_1^2+\sigma_2^2}}\right)+(b-d) \text{erf}\left(\frac{d-b}{\sqrt{2} \sqrt{\sigma_1^2+\sigma_2^2}}\right)\right)$$

so I believe the final answer should be:

$$\frac{1}{2(b-a)(d-c)}\left(\sqrt{\frac{2}{\pi }} \sqrt{\sigma_1^2+\sigma_2^2} \left(e^{-\frac{(b-c)^2}{2 \left(\sigma_1^2+\sigma_2^2\right)}}-e^{-\frac{(a-c)^2}{2 \left(\sigma_1^2+\sigma_2^2\right)}}\right)+\sqrt{\frac{2}{\pi }} \sqrt{\sigma_1^2+\sigma_2^2} \left(e^{-\frac{(a-d)^2}{2 \left(\sigma_1^2+\sigma_2^2\right)}}-e^{-\frac{(b-d)^2}{2 \left(\sigma_1^2+\sigma_2^2\right)}}\right)+(a-c) \text{erf}\left(\frac{c-a}{\sqrt{2} \sqrt{\sigma_1^2+\sigma_2^2}}\right)+(d-a) \text{erf}\left(\frac{d-a}{\sqrt{2} \sqrt{\sigma_1^2+\sigma_2^2}}\right)-(b-c) \text{erf}\left(\frac{c-b}{\sqrt{2} \sqrt{\sigma_1^2+\sigma_2^2}}\right)+(b-d) \text{erf}\left(\frac{d-b}{\sqrt{2} \sqrt{\sigma_1^2+\sigma_2^2}}\right)\right)$$