It indeed lies in the Schwartz space, which is the space $$
\mathcal{S} := \left\{f: \mathbb{R} \to \mathbb{R} \,:\,
\forall \alpha,\beta \in \mathbb{N}\,
\sup_{x\in\mathbb{R}} \left|x^\alpha f^{(\beta)}(x)\right| < \infty \right\}
$$
where $f^{(n)}$ denotes the $n$-th derivative of $f$. This definition looks a bit intimidating at first, but it's easy to understand. A function $f$ lies in $\mathcal{S}$ if it decays rapidly, i.e. faster than any polynomial grows (thus $x^nf(x)$ is bounded for every n), and if the same is true for every derivative of $f$.
This space is closed under fourier transform, i.e. if $f \in S$, then so is the fourier transform $\mathcal{F}(f)$. The same is true from the reverse fourier transform. The space of compactly supported $C^\infty$ functions is quite obviously a subspace of $\mathcal{S}$, hence it's image under $\mathcal{F}$ is thus also a subspace of $\mathcal{S}$.
The problem with your reasoing is that you're thinking about pointwise limits. Under those, $\mathcal{S}$ isn't complete, so knowing something about the limit doesn't tell you much about $\mathcal{S}$. There are topologies on $\mathcal{S}$ which make it a complete space, and under those you indeed cannot approximate the step function by functions in $\mathcal{S}$. The "trick" of those topologies is basically to force convergence of all the derivatives, which e.g. prevents you from approximating any discontinous functions.
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}.
$$
Best Answer
Since the inverse Fourier transform operator differs from the forward one merely by a complex conjugation in the integrand, ...
implies