Let $Q = [0,1]^2$. For sake of notation, let $$ f^{(i,j)}(x,\xi) = \frac{\partial^{i+j}}{\partial x^i \partial \xi^j}f(x,\xi). $$ Fix some non-negative integer $k$. Moreover let $f\in C^k(Q)$ if $$ \|f\|_{C^k} = \sum_{i+j \leq k}\sup_{(x,\xi)\in Q}|f^{(i,j)}(x,\xi)| < \infty. $$ Additionally, let $m = (m_1,m_2) \in \mathbf{Z}^2$. Define $ E_m = E_m(x,\xi) = e^{2\pi i m_1 x}e^{-2\pi i m_2 \xi}. $
$\textbf{Question.}$ Is $\{E_m\}_{m\in \mathbf{Z}^2}$ dense in $C^k(Q)$ for $k\geq 0$.
$\textbf{Ideas.}$ The question can be rephrased as follows. For any $f\in C^k(Q)$ and $\epsilon > 0$ does there exist some double trigonometric polynomial $$ p_k(x,\xi) = \sum_{m\in F}c_m E_m,\quad F \subseteq \mathbf{Z},\; |F| <\infty. $$ such that $$ \|f-p_k\|_{C^k} < \epsilon. $$
I am confident that for $k=0$ this result is already been established elsewhere. I've tried a similar technique in Are Trigonometric Functions Dense in $C^k(S^1)?$. However, the multi-index derivatives make it complicated to carry over. I've also thought about some potential "convolution" argument, but I haven't given it much time. Any advice would be appreciated or a counterexample if this does not hold.