convolution
Let $f,g$ be two continuous functions with compact support. Show that if $f$ and $g$ are not identically $0$, then neither is $f\ast g$.
This statement seems rather elementary, and I would prefer if the proof was also so, i.e avoiding reference to Titchsmarsh convolution theorem, and if possible not using Fourier transforms.
Perhaps this might soothe some of your discomfort.
Smoothing Action
There are many ways that convolution is useful in mathematics. First of all, as you have noted, $$ \mathrm{D}^\alpha\left(f\ast g\right)=\left(\mathrm{D}^\alpha f\right)\ast g\tag{1} $$ This is simply repeated changes of the order of integration and differentiation: $$ \frac{\mathrm{d}}{\mathrm{d}x}\int f(x-t)\,g(t)\,\mathrm{d}t=\int f'(x-t)\,g(t)\,\mathrm{d}t\tag{2} $$ This step can be justified in different ways depending on the context. For instance, if the limit which defines the derivative of $f$, $$ f'(x)=\lim_{h\to0}\frac{f(x+h)-f(x)}{h}\tag{3} $$ converges uniformly, then $(2)$ is valid for all $g\in L^1$.
Convolution combines the smoothness of two functions. That is, if both $f$ and $g$, and their first derivatives are in $L^1$, then the second derivative of their convolution is in $L^1$. This is because $f\ast g = g\ast f$, and so we can use $(2)$ twice to get $$ \frac{\mathrm{d}^2}{\mathrm{d}x^2}(f\ast g)=f'\ast g'\tag{4} $$
Fourier Analysis
Convolution plays an important role in Fourier Analysis. The key formulas demonstrate the duality between convolution and multiplication: $$ \mathscr{F}(f\ast g)=\mathscr{F}(f)\mathscr{F}(g)\quad\text{and}\quad\mathscr{F}(fg)=\mathscr{F}(f)\ast\mathscr{F}(g)\tag{5} $$ There also exists a duality between decay at $\infty$ and smoothness. Essentially, one derivative of smoothness of $f$ corresponds to one factor of $1/x$ in the decay of $\mathscr{F}(f)$, and vice versa.
The product of decaying functions decays even faster; e.g. $x^{-n}x^{-m}=x^{-(n+m)}$. The duality demonstrated in $(5)$ then says that the convolution of smooth functions is even smoother.
The Riemann-Lebesgue Lemma says that for $f\in L^1$, $$ \lim_{|x|\to\infty}\mathscr{F}(f)(x)=0\tag{6} $$ However, this is simply decay with no quantification. About all that can be said about $f,g\in L^1$ is that $f\ast g\in L^1$. However, if $f,g\in L^2$, then $f\ast g$ is continuous.
Summing Dice
Perhaps one of the earliest uses of convolution was in probability. If $f_n(k)$ is the number of ways to roll a $k$ on $n$ six-sided dice, then $$ f_n(k)=\sum_jf_{n-1}(k-j)f_1(j)\tag{7} $$ That is, for each way to achieve $k$ on $n$ dice, we must have $k-j$ on $n-1$ dice and $j$ on the remaining die. Equation $(7)$ represents discrete convolution.
The distribution function for the roll of a single six-sided die is evenly distributed among $6$ possibilities. This has discontinuities at $1$ and $6$ ($n=1$). The distribution function for the sum of two six-sided dice is the convolution of two of the one die distributions. This is continuous, but not smooth ($n=2$). The distribution function for the sum of three six-sided dice is the convolution of the one die and two dice distributions. This is smooth ($n=3$). For each die we add, we convolve one more of the one die distributions and the function gets smoother.
$\hspace{8mm}$
As $n\to\infty$, the distribution approaches a scaled version of the normal distribution: $\frac1{\sqrt{2\pi}}e^{-x^2/2}$.
Not in general. Let Consider the following function, where $n = 1,2,...$:
$$f(x) = \begin{cases} 0 & \text{ if } x < 0 \\ \frac{1}{n^3}& \text{ if }\frac{n^2-n}{2} \leq x < \frac{n^2+n}{2}\end{cases}$$
This function is $L^1$ with total integral $\frac{\pi^2}{6}$. However, it is constant and non-zero on arbitrarily long intervals. Thus for any $\phi$ as specified, there is sufficiently large $N$, such that for all $n >N$, there is some $x$ in the interval $\frac{n^2-n}{2} < x<\frac{n^2+n}{n}$ with $(f * \phi)(x) = \frac{1}{n^3}>0$.
Best Answer
A basic property of the Fourier transform is that it is not possible for $f$ and $\hat{f}$ to be both compactor supported (if $f\not \equiv 0$). To see this, note that
$$ \hat{f}(\xi) = \int_{[-R,R]^n} f(x) \cdot e^{2\pi i \langle x,\xi \rangle} \, dx $$
defines a holomorphic function on $\Bbb{C}^n$, as can be seen by differentiation under the integral sign (here $\mathrm{supp} (f) \subset [-R,R]^n$ for suitable $R>0$), so that the identity theorem for holomorphic functions implies that the support of $\hat{f}$ can not be compact. Even more, $\hat{f}$ can not vanish on any nondegenerate cube.
The convolution theorem implies
$$ \widehat{f\ast g} = \hat{f} \cdot \hat{g}, $$
where none of the two factors on the right can vanish on a nondegenerate cube. By continuity, $\hat{f} \cdot \hat{g}$ can not vanish on any nondegenerate cube.
Hence, $\widehat{f\ast g} \not\equiv 0$, so that $f\ast g \not\equiv 0$.
I am aware that this solution does not satisfy your requirement of avoiding the Fourier transform, but maybe it is better than having no proof at all.