Do flat functions with flat Fourier transform exist? Call a smooth function flat (in zero) if $f^{(k)}(0) = 0$ for all $k\ge 1$. Does there exist a flat Schwartz function $f\in\mathcal S(\mathbb R)$ such that its Fourier transform $\hat f = \mathcal F(f)$ is flat, too?
The construction of such a function has proven difficult to me. The simplest example of flat Schwartz functions are Bump functions which are constant around zero. However since, by the Paley-Wiener theorem, the Fourier transform of a bump function is real analytic, it cannot be flat as the Taylor series in zero is a constant.
Also note that the important identity
$$ \int_{-\infty}^{+\infty} f(x) x^k dx = \mu_k \hat f^{(k)}(0) $$
implies that all but the zeroth moment of the fourier transform of a flat function vanish. Typically the graph of a function whose moments are almost all zero looks similar to a sinc function: oscillating indefinitely and decaying as $|x|\to\infty$.
As a second remark, if such a function exists, it can w.l.o.g. be chosen as an Eigenfunction of the Fourier transform since for arbitrary $f$
$$ g =  f + \mathcal F(f) + \mathcal F^2(f) + \mathcal F^3(f) $$
is an Eigenfunction of the Fourier transform, and the property of the derivatives being $0$ in the origin carries over to $g$ (since $\mathcal F^2(f)(t) = f(-t)$). This leads to a slight reformulation
Question: Does there exist an (even) Schwartz function $f\in\mathcal S(\mathbb R)$ such that 


*

*$f$ is a density function: $\int_{-\infty}^{+\infty} f(x) dx = 1$

*All higher moments of $f$ vanish: $\int_{-\infty}^{+\infty} f(x)x^k dx = 0$ for all $k\ge 1$

*$f$ is an Eigenfunction of the Fourier transform: $\mathcal F(f) = f$
Any ideas, remarks, or references to relevant literature are greatly appreciated!
EDIT: 7th JUL
I am still interested in this question and have started a bounty. My progress towards finding such a function has been minimal. Although it is easy to find Schwartz functions that are flat - take e.g. $e^{-x^2-1/x^2}$ - it seems very hard to impossible to control both the moments in time and frequency space simultaneously. It doesn't help that concrete examples like the above have seemingly no closed form Fourier transform. Moreover forcing $f(0)=1$ makes thing even harder. The simplest example of such a  flat schwartz function I was able to find is $\exp(-\exp(x^2-1/x^2))$.
On the other hand there might be an avenue towards disproving the existence of such function. As pointed out above a flat function cannot be real analytic, but there are theorems connecting analyticity of the Fourier transform to properties like its decay rate. However the one I am familiar with (Schwartz-Payley-Wiener) deals with functions of compact support only. There seem to be generalizations with some interesting discussion found here https://mathoverflow.net/questions/23679/fourier-transform-of-analytic-functions
 A: Not sure why you want the $0$th moment to be $1$, this construction makes $f$ Schwartz such that all the moments of $f,\hat{f}$ vanish :
Take $h \in C^\infty_c([1/4,3/4])$ and let 
$$c_n = \int_0^1 h(x) e^{-2i \pi n x}dx, \qquad h(x) = \sum_{n=-\infty}^\infty c_n e^{2i \pi nx}, \qquad \forall k, \ \ \sum_{n=-\infty}^\infty c_n n^k = 0$$
Then take another $g \in C^\infty_c([1/4,3/4])$ and let 
$$f(x) = \sum_{n=-\infty}^\infty c_n g(x-n) \quad \in S(\mathbb{R})$$
With $\displaystyle d_m = \int_{-\infty}^\infty g(x) x^mdx$ we obtain
$$\int_{-\infty}^\infty f(x) x^k dx= \sum_{n=-\infty}^\infty  c_n \int_{-\infty}^\infty g(x) (x+n)^k dx = \sum_{m=0}^k {k \choose m}d_m \sum_{n=-\infty}^\infty  c_n n^{k-m} = 0$$
And the same holds for $\hat{f}$ since $f$ vanishes on $[-1/4,1/4]$.
A: I think a simple inductive construction will work, though the details are somewhat tedious. My function $f(x)=\sum g_j(x-a_j)$ will consist of (widely) separated bumps located at $0<a_1<a_2<\ldots$.
I must show that these can be chosen such that $\int x^n f(x)\, dx=0$ and $f\in\mathcal S$. I start out with an odd $g_1$, so $\int f=0$, which is good, but unfortunately $\int xf\not=0$, which isn't. To address this issue, I use $g_2$, which will be chosen orthogonal to $1$, and such that $\int x(g_1+g_2)=0$ also. In the next step, I use $g_3$, which will be chosen orthogonal to $1,x$, and such that $\int x^2(g_1+g_2+g_3)=0$ also.
The key observation is that the orthogonality that I impose to preserve earlier achievements is translation invariant in the sense that $g(x)$ is orthogonal to $1,x,\ldots , x^N$ precisely if $g(x-a)$ is. This means that by taking $a$ large, I can make $g_n$ as small as I want to in a given step, and this allows me to keep $f\in\mathcal S$ overall. (This requirement is really the list of conditions $|f^{(k)}(x^2+1)^m|\le 1$, $x\ge x_{N,k}$, and I satisfy these by paying attention to the first $N$ of these at a given step, and $N$ is increased gingerly.)
Remark: I'm no expert on these things, but I think it's known that there are measures (and functions, I would think) $f$ such that both $f$ and $\widehat{f}$ vanish on an interval, but I don't know if there are such examples in $\mathcal S$.
A: Here is some numerical evidence to the affirmative.
Let $g(x)=\exp(-x^2/2)$ be the Gaussian function such that $\mathcal Fg=g$. Let us define our function of interest as a linear combination of copies of $g$ stretched by powers of $2$ (factors in blue are mathematically superfluous but improve the numerical conditioning of the system we will form):
$$f_n(x) = \sum_{i=-n}^n \color{blue}{2^{-i^2}}a_ig(2^ix).$$
By choosing $a_{-i}=2^{-i}a_i$ for negative $i$, we satisfy $\mathcal Ff_n=f_n$. We are left with $n+1$ degrees of freedom $a_0,\dots,a_n$, so we may impose $n+1$ equations:
$$\begin{align}
\int_{\infty}^\infty f_n(x)\,\mathrm dx &= 1, \\
\color{blue}{\frac12}\,f''(0) &= 0, \\
&\vdots \\
\color{blue}{\frac{2^{-n^2}}{(2n - 1)!!}}\,f^{(2n)}(0) &= 0.
\end{align}$$
(Note that $f_n$ is an even function so all odd derivatives are automatically zero.) These equations form a linear system in the coefficients $a_i$, in fact a nearly tridiagonal one with diagonal entries converging to $\pm1$, which can be easily solved numerically. Let's look at the resulting functions:

Sure seems to be converging to something.
Here's a table of the values of $a_i$ for different $n$:
$$
\begin{array}{rcccccc}
 n & a_0 & a_1 & a_2 & a_3 & a_4 & a_5 \\
 0 & 0.398942 &   &   &   &   &   \\
 1 & 0.626909 & -0.303956 &   &   &   &   \\
 10 & 0.729176 & -0.465597 & 0.246237 & -0.124822 & 0.0626214 & -0.0313298 \\
 100 & 0.729177 & -0.465598 & 0.246238 & -0.124825 & 0.0626265 & -0.0313400 \\
\end{array}
$$
