Fourier transform of $\operatorname{erfc}^3\left|x\right|$

I have not yet managed to obtain the final answer and I am not even sure this is possible. Will keep trying. Below is some work which reduces the Fourier transform to an integral involving only one error function instead of three - this is the simplest expression I've found so far. In case someone would like to continue in this direction, the results were checked numerically.

Let us compute the integral $$I(\omega)=\int_0^{\infty}e^{i\omega x}\operatorname{erfc}^3x\,dx.\tag{1}$$ Obviously, the Fourier transform we are looking for is $\mathcal{F}(\omega)=I(\omega)+I(-\omega)$.

Integrating once by parts, we can write \begin{align}I(\omega)=&-\frac{1}{i\omega}-\frac{3e^{-\omega^2/4}}{i\omega}\left(-\frac{2}{\sqrt{\pi}}\right)\int_0^{\infty}e^{-(x-i\omega/2)^2}\operatorname{erfc}^2x\,dx=\\ =&-\frac{1}{i\omega}-\frac{3e^{-\omega^2/4}}{i\omega}\int_0^{\infty} \operatorname{erfc}^2x\,d(\operatorname{erfc}(x-i\omega/2))=\\ =&-\frac{1}{i\omega}+\frac{3e^{-\omega^2/4}}{i\omega}\operatorname{erfc}(-i\omega/2)-\frac{12e^{-\omega^2/4}}{i\omega\sqrt{\pi}}\int_0^{\infty}\operatorname{erfc}(x-i\omega/2)\operatorname{erfc} x\,e^{-x^2}dx.\end{align} Consider the following one-parameter deformation of the remaining integral: $$J(s,\omega)=\int_0^{\infty}\operatorname{erfc}(x-i\omega/2)\operatorname{erfc} (sx)\,e^{-x^2}dx.$$ Its partial derivative w.r.t. $s$ can be computed in explicit form: \begin{align} \frac{\partial J}{\partial s}=&-\frac{2}{\sqrt{\pi}}\int_0^{\infty}\operatorname{erfc}(x-i\omega/2)\,e^{-(1+s^2)x^2}xdx=\\ =&\frac{1}{\sqrt{\pi}\left(1+s^2\right)}\int_0^{\infty}\operatorname{erfc}(x-i\omega/2)\,d\left(e^{-(1+s^2)x^2}\right)=\\=&\frac{1}{\sqrt{\pi}\left(1+s^2\right)}\left[-\operatorname{erfc}\left(-\frac{i\omega}{2}\right)+\frac{2}{\sqrt{\pi}}\int_0^{\infty}e^{-(1+s^2)x^2 -\left(x-\frac{i\omega}{2}\right)^2}dx\right]=\\ =&\frac{1}{\sqrt{\pi}\left(1+s^2\right)}\left[-\operatorname{erfc}\left(-\frac{i\omega}{2}\right)+\frac{\exp\left\{\frac{1+s^2}{2+s^2}\frac{\omega^2}{4}\right\}\operatorname{erfc}\left(-\frac{i\omega}{2\sqrt{2+s^2}}\right)}{\sqrt{2+s^2}}\right] \end{align} Integrating this back, we find the quantity we are really looking for: \begin{align} J(1,\omega)=&-\int_1^{\infty}\frac{\partial J}{\partial s}ds=\\ =&\frac{\sqrt{\pi}}{4}\operatorname{erfc}\left(-\frac{i\omega}{2}\right) -\frac{1}{\sqrt{\pi}}\int_1^{\infty}\frac{\exp\left\{\frac{1+s^2}{2+s^2}\frac{\omega^2}{4}\right\}\operatorname{erfc}\left(-\frac{i\omega}{2\sqrt{2+s^2}}\right)}{(1+s^2)\sqrt{2+s^2}}ds \end{align} This gives \begin{align} I(\omega)=&-\frac{1}{i\omega}+\frac{12}{i\omega\pi}\int_1^{\infty}\frac{\exp\left\{-\frac{\omega^2}{4(2+s^2)}\right\}\operatorname{erfc}\left(-\frac{i\omega}{2\sqrt{2+s^2}}\right)}{(1+s^2)\sqrt{2+s^2}}ds=\\ =&-\frac{1}{i\omega}+\frac{48}{i\omega\pi}\int_0^{\frac{\sqrt{3}}{6}} \frac{e^{-\omega^2u^2}\operatorname{erfc}(-i\omega u) \, udu}{(1-4u^2)\sqrt{1-8u^2}}. \end{align} Therefore the Fourier transform can be written as $$\boxed{\displaystyle\mathcal{F}(\omega)=\frac{96}{i\omega\pi}\int_0^{\frac{\sqrt{3}}{6}} \frac{e^{-\omega^2u^2}\operatorname{erf}(i\omega u) \, udu}{(1-4u^2)\sqrt{1-8u^2}}}\tag{2}$$

Further thoughts:

Obviously, using (2), one can compute an arbitrary number of terms in the Taylor expansion of $\mathcal{F}(\omega)$. Here are several first terms: \begin{align} \mathcal{F}(\omega)=\frac{6}{\pi^{3/2}}\left[(\pi-2\alpha)-\frac{2+6\pi-15\alpha}{36}\omega^2+\frac{74+144\pi-387\alpha}{8640}\omega^4+\ldots\right], \end{align} where $\alpha=\sqrt{2}\arcsin\sqrt{\frac23}$. This structure of coefficients persists in higher order terms: they are linear combinations of $1$, $\pi$ and $\alpha$ with rational coefficients. If one could identify the corresponding rational sequences, it would probably become possible to sum up the series (maybe to some hypergeometric expressions).

Actually, if we write $\mathcal{F}(\omega)=\sum_{k=0}^{\infty}\mathcal{F}_k\omega^{2k}$, it is possible to derive the following exact expression for the Taylor coefficients: $$\mathcal{F}_k=\frac{1}{(2k+1)!}\frac{6}{\pi^{3/2}}\left[\frac{\partial^k}{\partial p^k}\left(\frac{\pi}{p}-\frac{4\arctan\sqrt{1+p}}{p\sqrt{1+p}}\right)\right]_{p=1}\tag{3}$$ From this we find that the ''$\pi$-contribution'' (corresponding to the derivatives of $\frac{\pi}{p}$ in the last formula) in the above reasoning sums up to $$\mathcal{C}_{\pi}(\omega)= \frac{6e^{-\omega^2/4}}{i\omega}\operatorname{erf}\left(\frac{i\omega}{2}\right).$$

The ''$\alpha$-contribution'' corresponds to the derivatives of the second term in (3) which do not act on $\arctan\sqrt{1+p}$. Here again it is possible (the calculation is somewhat involved but relatively straightforward, so I will only present it if somebody will be explicitly interested) to sum up the corresponding series to $$\mathcal{C}_{\alpha}(\omega)=-\frac{24\arctan \sqrt2}{i\pi\omega}e^{-\omega^2/4}\operatorname{erf}\left(\frac{i\omega}{2\sqrt{2}}\right)$$

There remains ''$1$-contribution'' to determine (working on it). Hence so far the statement is that $$\boxed{\displaystyle\mathcal{F}(\omega)=\mathcal{C}_{\pi}(\omega)+\mathcal{C}_{\alpha}(\omega)+\frac{\omega^2}{3\pi^{3/2}}\sum_{k=0}^{\infty}(-1)^{k+1}r_k\omega^{2k}}$$ where $r_k$ are positive $rational$ numbers (the normalization is chosen so that $r_0=1$).