There's a nice proof of this in Hörmander's book "The Analysis of Linear Partial Differential Operators I" (Theorem 7.1.14). The main tool needed is a distributional version of Fubini's theorem (Theorem 5.1.1 in Hörmander's book). I'll summarize the main points:
Preliminary step 1: Fubini's theorem
If $u$ and $v$ are distributions on ${\mathbb R}^m$ and ${\mathbb R}^n$ respectively, then there is a unique product distribution $u \otimes v$ on ${\mathbb R}^{m+n}$ characterized by the condition
$$
(u \otimes v)(\varphi \otimes \psi) = u(\varphi) v(\psi)
$$
for all $\varphi \in {\mathcal D}({\mathbb R}^m)$ and $\psi \in {\mathcal D}({\mathbb R}^n)$, where we write $(\varphi \otimes \psi)(x,y) := \varphi(x) \psi(y)$. Moreover, $u \otimes v$ can be evaluated on an arbitrary test function $\varphi \in {\mathcal D}({\mathbb R}^{m+n})$ by
$$
(u \otimes v)(\varphi) = u(x \mapsto v(\varphi(x,\cdot))) = v(y \mapsto u(\varphi(\cdot,y))).
$$
(Note that there is a slightly nontrivial exercise required to show that both expressions on the right make sense, e.g. that $x \mapsto v(\varphi(x,\cdot))$ defines a smooth compactly supported function on ${\mathbb R}^m$. This rests mainly on the fact that by uniform continuity, $x \mapsto \varphi(x,\cdot)$ is a continuous map to the space of test functions, and other arguments of this sort.) If you write down this formula in the case where $u$ and $v$ are given by locally integrable functions, you'll find that it follows easily from the classical Fubini's theorem.
If you know that the product distribution is unique, then the formula follows by verifying directly that both expressions on the right hand side define distributions that satisfy the defining property of a product distribution. Uniqueness can be proved via mollification: if $u \otimes v$ were not unique, then there would exist a nontrivial distribution $w$ on ${\mathbb R}^{m+n}$ such that $w(\varphi \otimes \psi) = 0$ for all $\varphi \in {\mathcal D}({\mathbb R}^m)$ and $\psi \in {\mathcal D}({\mathbb R}^n)$. Choose approximate identities, i.e. sequences of smooth functions $\rho_j : {\mathbb R}^m \to [0,\infty)$ and $\sigma_j : {\mathbb R}^n \to [0,\infty)$ with shrinking compact supports near $\{0\}$ that converge in the space of distributions to $\delta$-functions. Then the classical Fubini theorem implies that the sequence $\rho_j \otimes \sigma_j : {\mathbb R}^{m+n} \to [0,\infty)$ also defines an approximate identity in the same sense, and it follows that the sequence of smooth functions $(\rho_j \otimes \sigma_j) * w$ converges to $w$ in the space of distributions. But those functions are all $0$ due to the defining property of $w$, thus $w=0$.
Preliminary step 2: polynomial growth
Before we can view the function $g(\xi) := u(\chi(x) e^{-i x \xi})$ as a plausible candidate for the Fourier transform of $\chi u$, we need to know that it behaves reasonably enough at infinity to define a tempered distribution. As I indicated in my comments on the previous answer, $g$ is definitely not a Schwartz function in general, but one can show that it has polynomial growth. Perhaps the quickest way is to use standard properties of the Fourier transform and rewrite $g$ as
$$
g(\xi) = \left( ( {\mathcal F}\chi)^- * {\mathcal F}^*u\right)(-\xi),
$$
where I'm using the notation $f^-(x) := f(-x)$. As the convolution of a Schwartz function with a tempered distribution, it follows from standard results about the convolution that this function has polynomial growth.
The main argument
As stated in the question, we need to prove that the relation
$$
\int_{{\mathbb R}^n} u(\chi(x) e^{-i x \xi}) \phi(\xi) \, d\xi = u\left( \int_{{\mathbb R}^n} \chi(x) e^{- i x \xi} \phi(\xi)\, d\xi \right)
$$
holds for every $u \in {\mathcal D}'({\mathbb R}^n)$, $\chi \in {\mathcal D}({\mathbb R}^n)$ and $\phi \in {\mathcal S}({\mathbb R}^n)$. By step 2, we already know that both sides give well-defined tempered distributions when regarded as functionals of $\phi$, so by density, it will suffice to assume $\phi \in {\mathcal D}({\mathbb R}^n)$. The key observation now is that by the theorem in step 1, both sides can be identified with
$(u \otimes 1)(f)$, where $1 \in {\mathcal D}'({\mathbb R}^n)$ is the distribution $1(\varphi) := \int_{{\mathbb R}^n} \varphi(x)\, dx$ and $f \in {\mathcal D}({\mathbb R}^{m+n})$ is given by
$$
f(x,\xi) := \chi(x) \phi(\xi) e^{-i x \xi}.
$$
