First, let me say that you can show this pretty easily be writing the convolution as a matrix vector product, with the matrix being a circulant one defined by one of the vectors. The result falls out due to the DFT diagnolizing circulant matrices.
Anyway, you can also show this directly substituting the discrete convolution formula, and playing with indices. We need one result before we do this, however. I'm assuming from here out that we're in $\mathbb{R}^N$, and that addition is mod $N$.
Proposition
Let $x^m$ be the vector $x$ shifted by $m$, i.e., $x^m_j = x_{m+j}$. Then $\langle x, y \rangle = \langle x^m, y^m \rangle$.
Note, this says that
$$
\langle x, y \rangle = \sum_{j=0}^{N-1} x_j y_j = \sum_{j=-m}^{N-1-m} x_{j+m} y_{j+m} = \sum_{j=0}^{N-1} x_{j+m} y_{j+m} = \langle x^m, y^m \rangle
$$
proof
Without loss of generality, assume $m \in \{0, 1, \ldots N-1\}$. Let $P^m$ be the identity matrix whose rows are circularly permuted by $m$.
To be precise, rows $m, m+1,\ldots, N-1$ of $P^m$ are the rows $0, 1, \ldots, N-1-m$ rows of the identity, while rows $0,1, \ldots, m-1$ are rows $N-m, \ldots, N-1$ of the identity.
For example,
$$
P^2 = \begin{bmatrix} e_{N-2} \\ e_{N-1} \\ e_{0} \\ e_{1} \\ \vdots \\ e_{N-3} \end{bmatrix}
$$
where $e_j$ is the $j$th row of the identity.
It should be clear that $x^m = P^m x$. Furthermore, $P^m$ is unitary so it preserves inner products, i.e. $\langle x, y \rangle = \langle P^mx, P^my \rangle$.
$$\tag*{$\blacksquare$}$$
This should be intuitively clear. Since we're working with mod $N$, the shift just rearranges the terms in the sum circularly, not changing the overall value of the sum.
Now, back to the main result. Let's call $c_n = (\mathbf{a}*\mathbf{b})_n$.
\begin{align*}
\hat{c}_m &= \sum_{n=0}^{N-1} \exp(-2\pi imn / N) c_n \\
&= \sum_{n=0}^{N-1} \exp(-2\pi imn / N) \sum_{k=0}^{N-1} a_kb_{n-k} \\
&= \sum_{k=0}^{N-1} a_k \sum_{n=0}^{N-1} \exp(-2\pi imn / N) b_{n-k}
\end{align*}
Now, let $n'=n-k$. So,
\begin{align*}
\hat{c}_m &= \sum_{k=0}^{N-1} a_k \sum_{n'=-k}^{N-1-k} \exp(-2\pi im(n'+k) / N) b_{n'} \\
&= \sum_{k=0}^{N-1} a_k \exp(-2\pi imk / N) \sum_{n'=-k}^{N-1-k} \exp(-2\pi imn' / N) b_{n'} \\
&= \left( \sum_{k=0}^{N-1} a_k \exp(-2\pi imk / N) \right) \left(\sum_{n=0}^{N-1} \exp(-2\pi imn / N) b_{n} \right) \\
&= \hat{a}_m \hat{b}_m
\end{align*}
where the second to last line is due to the above proposition.
