What would be the most accurate way to decompose $\mathrm{erf}(x-y)$ into pairs of separable functions? Or in another word, we would like to find $f_i(x)$ and $g_i(y)$ such that with $$ \mathrm{erf}(x-y) \approx \sum_{i=1}^N f_i(x)g_i(y), $$ the approximation error decreases if we increase $N$.
What I have done
I have considered using Taylor series of $\mathrm{erf}$ around $0$, which is $$ \mathrm{erf}(x - y) = \frac{2}{\sqrt{\pi}}\sum_{n=0}^\infty \frac{(-1)^n (x-y)^{2n+1}}{n!(2n+1)}, $$ where $(x - y)^n$ can be expanded and be decomposed into $\sum_i f_i(x) g_i(y)$. However, this naive Taylor series approach requires a lot of decomposed terms to get a reasonable accuracy, especially for large $|x - y|$. Also, for each Taylor series term (i.e., for each value of $n$), it requires $n + 1$ decomposition terms (e.g., $(x - y)^2$ can be decomposed into 3 terms: $x^2 - 2xy + y^2$). Therefore, including 10 Taylor series terms, for example, will actually produce 110 decomposition terms.
I'm looking for some decomposition that can converge relatively quickly within a few terms. For example, $\exp(-(x-y)^2 / 2)$, can be decomposed into: $$ \exp(-(x-y)^2 / 2) = \sum_{n=0}^\infty \frac{1}{n!} \left(x^n e^{-x^2/2}\right) \left(y^n e^{-y^2/2}\right). $$ The expression above can be obtained by expanding the square inside exp, and taking the Taylor series for $e^{xy}$. In this exponential case, 10 terms is enough to get error less than $10^{-3}$ for reasonable values of $x$ and $y$. However, I can't find the similar decomposition terms for $\mathrm{erf}$. Any helps are appreciated!