The goal is to compute $\hat{f}[\varphi]$ for $f=cos^4(x)$ and $\varphi \in \mathscr{S\mathbb{(R}})$ (Schwartz-Space)
My plan was to switch integrals with Fubini (since $\varphi \in \mathscr{S(\mathbb{R})} $ we have $\hat\varphi \in \mathscr{S(\mathbb{R})} $ and hence $\int_\mathbb{R\times\mathbb{R}}|cos^4(x)\varphi(k)e^{-ikx}|dxdk<\infty$ )
$$\hat{f}[\varphi]=f[\hat{\varphi}]=\int_\mathbb{R} cos^4(x) \hat\varphi(x) dx=\int_\mathbb{R}cos^4(x) \int_\mathbb{R}\varphi(k)e^{-ikx}dk \, dx \overset{Fubini}{=} \int_\mathbb{R}\left(\int_{\mathbb{R}}cos^4(x)e^{-ikx}dx \right)\varphi(k)dk= \frac{1}{2^4} \int_\mathbb{R}\left(\int_{\mathbb{R}}(e^{ix}+e^{-ix})^4e^{-ikx}dx \right)\varphi(k)dk=...=\int_\mathbb{R}g(k) \varphi(k)dk=g[\varphi] \Rightarrow \hat f=g$$ but I don't see how to fill the dots since by developing the binomial $(e^{ix}+e^{-ix})^4=e^{4ix}+4e^{2ix}+6+4e^{-2ix}+e^{-4ix}$ I get non-integrable functions.
You're right, that (at least without using lots of imagination) the change of order of integration is not legitimate, and the computation gets stuck. We have to treat $\cos^4 x$ as a tempered distribution, as you anticipate.
As you anticipate, expressing $$\cos^4 x\;=\;\Big({1\over 2}(e^{ix}+e^{-ix})\Big)^4\;=\;{1\over 2^4}(e^{4ix}+4e^{2ix}+6+4e^{-2ix}+e^{-4ix}) $$ in terms of exponentials helps isolate the issues: it suffices to determine the Fourier transform of (the tempered distribution) $e^{icx}$, with real $c$. Of course, trying to do this directly runs into a similar issue as you already found.
One way or another, we might know/anticipate that Fourier transforms of (purely imaginary) exponentials are (constant multiples of) translates of Dirac $\delta$. Thus, to pin down some of the details/constants, we compute the inverse Fourier transform of "$\delta(x-x_o)$", meaning $\delta$ translated by $x_o$, even though pointwise values do not make sense.
A worthwhile point regarding rigor is that, since $\delta(x-x_o)$ is a compactly-supported distribution, its inverse Fourier transform is (up to normalization) obtained just by evaluating it on exponentials: $$ \check{\delta}(x-x_o)(\xi)\;=\;\delta(x-x_o)(e^{ix\xi}) \;=\; e^{ix_o\xi} $$ Maybe you want to divide by $2\pi$ or $\sqrt{2\pi}$, depending on your normalization. By Fourier inversion for tempered distributions, up to the $2\pi$ normalization stuff, $$ \widehat{e^{ix_o x}} \;=\; \delta(x-x_o) $$ So the Fourier transform of $\cos^4 x$ is a linear combination of $\delta(x-4)$, $\delta(x-2)$, $\delta(x)$, $\delta(x+2)$. and $\delta(x+4)$.