Is every non-negative test function the limit of a sequence of sums of squares of test functions? Let $0\leq f\in\mathscr{D}(\mathbb{R}^n)$. As shown e.g. by J.-M. Bony, F. Broglia, F. Colombini and L. Pernazza, Nonnegative functions as squares or sums of squares, J. Funct. Anal. 232 (2006) 137-147 (see also this MO question and this math.SE question), not every such an $f$ is of the form $f=g^2$ for $g\in\mathscr{D}(\mathbb{R}^n)$ or even a finite sum of such. On the other hand, as mentioned in the paper above, C. Fefferman and D. H. Phong sketched a proof (On positivity of pseudo-differential operators, Proc. Natl. Acad. Sci. U.S.A. 75 (1978) 4673-4674) of the fact that any $0\leq f\in\mathscr{C}^\infty(\mathbb{R}^n)$ can be written as a sum of squares $$f=\sum^k_{j=1}g_j^2$$ with $g_j\in\mathscr{C}^{1,1}(\mathbb{R}^n)$ (i.e. $g_j$ is a differentiable function whose derivatives are locally Lipschitz) for all $1\leq j\leq k$ for some $k\in\mathbb{N}$. This fact was a key ingredient of the proof of the important inequality for scalar pseudodifferential operators with non-negative symbols that bears their name. For a modern, more detailed proof of the above formula, see N. Lerner, Some Facts About the Wick Calculus, in L. Rodino, M.W. Wong (eds.), Pseudodifferential Operators: Quantization and Signals, Lectures given at the C.I.M.E. Summer School held in Cetraro, Italy June 19–24, 2006 (Springer Lecture Notes in Mathematics 1949, 2008), pp. 135-174, particularly Theorem 5.2, pp. 167-172, and the discussion right after Theorem 1.1 in the paper by Bony et alii above on page 139. It immediately follows by multiplication by squares of smooth bump functions that the $g_j$'s may be chosen to be compactly supported if $f\in\mathscr{D}(\mathbb{R}^n)$. Bony et alii showed above that this regularity for $n\geq 4$ is sharp.
All this leads naturally to the following

Question: Is every $0\leq f\in\mathscr{D}(\mathbb{R}^n)$ the limit of a 
  sequence of sums of squares in $\mathscr{D}(\mathbb{R}^n)$ in the latter's topology?

In other words, is the cone of non-negative elements of $\mathscr{D}(\mathbb{R}^n)$ the (sequential) closure of the cone of sums of squares in $\mathscr{D}(\mathbb{R}^n)$?
I am particularly interested in arguments that do not rely on the result by Fefferman and Phong.

EDIT - Follow-up subquestion: (suggested by André Henriques) Is every $0\leq f\in\mathscr{D}(\mathbb{R}^n)$ the limit of an increasing sequence of sums of squares in $\mathscr{D}(\mathbb{R}^n)$ in the latter's topology? Particularly, can $f$ be written as $$f=\sum^\infty_{j=1}g_j^2$$ with $g_j\in\mathscr{D}(\mathbb{R}^n)$ for all $j$ and convergence in $\mathscr{D}(\mathbb{R}^n)$?

 A: Given  $f\in\mathcal{D}(\mathbb{R}^n)$, $f\ge0$ choose a $g\in\mathcal{D}(\mathbb{R}^n)$, $g\ge0$ such that $\mathrm{supp}(f)\subset \{g>0\}$, and let $\epsilon>0$. Then $\sqrt{f+\epsilon^2 g^2}$ is $C^\infty$: by composition, at any point where $g(x)>0$, and because it locally coincides with $\epsilon g$, at any point where $g(x)=0$. 
A: It is possible I have misunderstood the topology on $\mathcal{D}(\mathbb{R}^n)$; I assumed it means that all the derivatives must converge. With that understanding, let $a \in (0,3]$.  I claim that  $x^4 y^2+y^4 z^2+z^4 x^2 - a x^2 y^2 z^2$ cannot be written as a sum $\sum f_i(x,y,z)^2$ of $C^3$ functions. Note that this function is nonnegative by the arithmetic-geometric mean inequality. 
If we had $x^4 y^2+y^4 z^2+z^4 x^2 - a x^2 y^2 z^2 = \sum f_i(x,y,z)^2$ then the degree $k$ part of the Taylor series of the right hand side would converge to the degree $k$ part of the Taylor series of the left hand side. From this we see that the $f_i$ vanish up to third order and, writing $c_i$ for the cubic part of the Taylor series of $f_i$, we have $x^4 y^2+y^4 z^2+z^4 x^2 - a x^2 y^2 z^2 = \sum c_i^2(x,y,z)$. This is the standard example of a nonnegative polynomial which cannot be written as a sum of squares, so now we just hav to show that the same proof shows that it cannot be written as an infinite sum of squares.
There are no $x^6$, $y^6$ or $z^6$ terms on the right hand side, so the $c_i$ have no $x^3$, $y^3$ or $z^3$ terms. Similarly, looking at the coefficients of $x^2 y^4$, $y^2 z^4$ and $z^2 z^4$, we see that there are no $x y^2$, $y z^2$ or $z x^2$ terms. So $c_i(x,y,z) = p_i x^2 y + q_i y^2 z + r_i z^2 x + s_i x y z$. But then $-a = \sum s_i^2$, a contradiction.
A: The last counterexample  posted here led me to think a bit on the theme "What can be done then?". I came out with the following simple crank machine, that crumbles a positive $C^1_c$ function $f$ into a $C^1$ convergent series of squares, and that could also work in other ordered algebras (in fact I suspect it may be a known procedure). We may assume with no loss of generality $0\le f\le1/2$.
Consider the sequence of functions:
$$\begin{cases} g_1=f\\
g_{n+1}=g_n-g_n^2
\end{cases}$$
It is immediate by the definition that $f=\sum_{k=1}^n g_k^2 +g_{n+1}$, so that we have $f=\sum_{k=1}^\infty g_k^2$ in various senses, provided that the remainder  $g_{n}$ converges to zero in the same sense.
Uniform convergence of $g_n$: We have $0\le g_n\le1$ and
${1\over g_{n+1}}={1\over g_{n}}+{1\over 1-g_{n}}\ge {1\over g_{n}}+1,$ hence ${1\over g_{n}}\ge n$ and $g_n\le {1\over {n}}$.
Uniform convergence of $\partial_i g_n$: We have $\partial_i g_{n+1}=(1-2g_n)\partial_ig_n$. Since $0\le 1-2g_n\le1-g_n={g_{n+1}\over g_n}\le 1$ we also have $|\partial_i g_{n+1}|\le {g_{n+1}\over g_n}|\partial_ig_n|\le |\partial_ig_n|$. Iterating the latter inequalities we get
$$|\partial_i g_{n}|\le {g_{n}\over f}|\partial_if|\le|\partial_i f|$$
which immediately implies the uniform convergence of $\partial_i g_{n}$ to $0$ (indeed $\partial_i g_{n}$ converges uniformly to zero on any set $\{f\ge\epsilon\}$ and it is dominated by $|\partial_i f|$).
$$*$$
Rmk if for instance $f(x)=x^2+o(x^2)$ for $x\to0$ then also $g_n(x)=x^2+o(x^2)$ so that $g_n''(0)$ can't converge to zero, as illustrated in the last post.
On the other hand, by David Speyer's example, we know that, in general, for a smooth $f$, no series can't be convergent in $C^3$.
A: This is a comment on Pietro Majer's answer, with graphics. The first figure shows $x^2$ (red) and the first five functions of Pietro's method approaching it (blue). The second figure shows the double derivatives thereof. As you can see, the functions are converging but their second derivatives are not.


