Estimate of $\int_y^{+∞}\frac{[2-({\rm erf}(x+a)+{\rm erf}(x-a))]^2}{e^{-(x+a)^2}+e^{-(x-a)^2}}\,\mathrm dx$ $\DeclareMathOperator{\erf}{erf}\DeclareMathOperator{\erfc}{erfc}$I am trying to find either a closed form expression or an $y$-dependent analytical estimate for
$$\int_y^{+\infty}\frac{[2-(\erf(x+a)+\erf(x-a))]^2}{e^{-(x+a)^2}+e^{-(x-a)^2}}\,\mathrm dx.$$
This appeared in my research. It is related to the diffusion equation under a specific potential. All I could come up with is to define
$$f(x)=\erfc(x+a)+\erfc(x−a)$$
and then I see that this integral is
$$-\frac{2}{\sqrt{\pi}}\int_{y}^{+\infty}\frac{f(x)^{2}}{\frac{\mathrm df(x)}{\mathrm dx}}\,\mathrm dx.$$
Any suggestions are welcome!
 A: $\DeclareMathOperator{\erf}{erf}\def\d{\mathrm{d}}\def\e{\mathrm{e}}\def\paren#1{\left(#1\right)}$Assume $a > 0$. For any $ν \in \mathbb{R}$, by l'Hospital's rule,\begin{gather*}
\lim_{x → +∞} \frac{\displaystyle\int_x^{+∞} t^ν \exp(-t^2) \,\d t}{x^{ν - 1} \exp(-x^2)} = \lim_{x → +∞} \frac{-x^ν \exp(-x^2)}{-(2x^ν - (ν - 1) x^{ν - 2}) \exp(-x^2)}\\
= \lim_{x → +∞} \frac{x^2}{2x^2 - (ν - 1)} = \frac{1}{2},
\end{gather*}
thus when $x → +∞$,\begin{gather*}
1 - \erf(x) = \frac{2}{\sqrt{π}} \int_x^{+∞} \exp(-t^2) \,\d t = \frac{1}{\sqrt{π} x} \exp(-x^2) + o\paren{ \frac{1}{x} \exp(-x^2) }, \tag{1}\\
\int_x^{+∞} \frac{1}{t^2} \exp(-t^2) \,\d t \sim \frac{1}{2x^3} \exp(-x^2). \tag{2}
\end{gather*}
Note that $\exp(-(x + a)^2) = o(\exp(-(x - a)^2))$ ($x → +∞$). Therefore when $x → +∞$,\begin{align*}
&\mathrel{\phantom=} 2 - \erf(x + a) - \erf(x - a)\\
&= (1 - \erf(x - a)) + (1 - \erf(x + a))\\
&\stackrel{(1)}{=} \frac{\exp(-(x - a)^2)}{\sqrt{π} (x - a)} + o\paren{ \frac{\exp(-(x - a)^2)}{x - a} } + \frac{\exp(-(x + a)^2)}{\sqrt{π} (x + a)} + o\paren{ \frac{\exp(-(x + a)^2)}{x + a} }\\
&= \frac{\exp(-(x - a)^2)}{\sqrt{π} (x - a)} + o\paren{ \frac{\exp(-(x - a)^2)}{x - a} } \sim \frac{\exp(-(x - a)^2)}{\sqrt{π} (x - a)},
\end{align*}
which implies that$$
\frac{(2 - \erf(x + a) - \erf(x - a))^2}{\exp(-(x + a)^2) + \exp(-(x - a)^2)} \sim \frac{\paren{ \dfrac{\exp(-(x - a)^2)}{\sqrt{π} (x - a)} }^2}{\exp(-(x - a)^2)} = \frac{\exp(-(x - a)^2)}{π(x - a)^2}.
$$
Therefore when $y → +∞$,\begin{gather*}
\int_y^{+∞} \frac{(2 - \erf(x + a) - \erf(x - a))^2}{\exp(-(x + a)^2) + \exp(-(x - a)^2)} \,\d x\\
\sim \int_y^{+∞} \frac{\exp(-(x - a)^2)}{π(x - a)^2} \,\d x = \frac{1}{π} \int_{y - a}^{+∞} \frac{1}{t^2} \exp(-t^2) \,\d t \stackrel{(2)}{\sim} \frac{\exp(-(y - a)^2)}{2π(y - a)^3}.
\end{gather*}
