$C_0^\infty(X) \otimes C_0^\infty(Y)$ is dense in $C_0^\infty(X\times Y)$ My professor told me that one can approximate a function $\varphi$ in $C_0^\infty(X \times Y)$ via
$$\varphi(x,y)=\lim_{n \to \infty} \sum_i \phi_i^n(x) \psi_i^n(y) $$
where $\phi_i^n$ is sequence in $C_0^\infty(X)$ and $\psi_i^n$ a sequence in $C_0^\infty(Y)$. You can imagine $X,Y$ sufficiently nice, in my case they where open and bounded subsets of $\mathbb{R}^n$, $n \in \mathbb{N}$.
A property like this is essential if one wants to prove a parabolic PDE in distributional sense and in one step separates the time and space test functions. 
I asked for a proof of this statement and he told me that one can find proofs of this statement in books that are older than him and right now, he does not have any reference. Furthermore, I looked up some PDE books like Evans, Renardy, etc, but unfortunately I couldn't find such a statement. 
So I would like to ask you to kindly tell me a source so I can look up the proof.
Thank you so much!
 A: I don't have a reference, but IIRC the following proof sketch can be used:


*

*every function is approximately a sum of bump functions (with small domains)

*bump functions of the form $f(x) g(y)$ are sufficient


If you're careful, you can probably adapt the sketch to show there exists a sequence of functions with
$$ \varphi(x, y) = \sum_{i=0}^{\infty} \phi_i(x) \psi_i(y) $$
The idea being once you pick the terms up to $\phi_k$ and $\psi_k$, you dedicate the rest of the sequence to approximating $\varphi(x,y) - \sum_{i=0}^k \phi_i(x) \psi_i(y)$.
Maybe this form is even easier to work with; e.g. show that there exist finitely many functions so that $\| \varphi(x,y) - \sum_{i=0}^k \phi_i(x) \psi_i(y) \| < \frac{1}{2} \| \varphi(x,y) \|$ by your favorite norm, and recurse.

Maybe there's a snazzy Fourier transform-like approach (wavelet transform?) that lets you just write down such a sum. e.g.


*

*Do a transform in one variable and see the coefficients as functions in the other

*Do a two-dimensional transform, to express a general function as a linear combination of basis functions of shape $f(x) g(y)$.


Alas, my functional analysis isn't up to doing such a thing rigorously, so I can only suggest it.
A: Source: Theorem 39.2 (page 409) in Topological Vector Spaces, Distributions and Kernels by F. Treves.

Remark. Corollay 2, used in the proof, says that "every function $f\in C^m_c(\mathbb R^n)$ is the limit, in $C^m(\mathbb R^n)$, of a
sequence of polynomials".
