Long story short, I only get $C(n,f,\varepsilon) \, a^{-n (k-\varepsilon n-\frac{1}{2})}$ as a bound for $C^k$-functions for any $\varepsilon > 0$.
I get a nice bound for every $f$ with absolutely summable Fourier coefficients (the set of such functions is called the Wiener-Algebra and it is a super-set of $C^1$). Unfortunately the bounds are in terms of the Fourier coefficients and not in terms of a $C^k$-norm, though there are connections, as I mention at the end.
I am assuming $a \in \mathbb{N}$. You will have to add the bound for $|I(f,a,n) - I(f,\lceil a \rceil,n)|$, where
\begin{align*}
& I(f,a,n) = \int_0^1 f(x)f(ax)\cdots f(a^nx) \, dx
\end{align*}
yourself.
For $f \in C^k$, we can of course use the Fourier series
\begin{align*}
& f(x) = \sum\limits_{r \in \mathbb{Z}} \hat{f}_r e^{i2\pi rx}
\end{align*}
with the Fourier coefficients
\begin{align*}
& \hat{f}_r = \int_0^{1} f(x) e^{-i2\pi r x} \, dx \ .
\end{align*}
We can simply calculate
\begin{align*}
& I(f,a,n) = \sum\limits_{r_0,...,r_n \in \mathbb{Z}} \bigg( \prod\limits_{j=0}^n \hat{f}_{r_j} \bigg) \underbrace{\int_0^1 e^{i2\pi x (r_0+ar_1+...+a^nr_n)} \, dx}_{= \mathbb{1}_{0 = r_0+ar_1+...+a^nr_n}} \ .
\end{align*}
The assumption $0 = \int_0^1 f(x) \, dx$ is equivalent to $\widehat{f}_0=0$ and one can inductively show that $0 = r_0+ar_1+...+a^nr_n$ implies that (either $r_j=0$ for some $j$, which then does not contribute to the sum, or) at least one $t \in \{0,...,n\}$ must satisfy $|r_t| \geq a$. It follows that
\begin{align*}
& |I(f,a,n)| \leq \sum\limits_{\substack{r_0,...,r_n \in \mathbb{Z} \\ 0 = r_0+ar_1+...+a^nr_n}} \prod\limits_{j=0}^n |\hat{f}_{r_j}|\\
& \leq (n+1) \sum\limits_{\substack{r_0,...,r_n \in \mathbb{Z} \\ |r_0| \geq a}} \prod\limits_{j=0}^n |\hat{f}_{r_j}|\\
& \leq (n+1) \bigg( \sum\limits_{r \in \mathbb{Z}} |\widehat{f}_r| \bigg)^n \sum\limits_{|r_0| \geq a} |\widehat{f}_{r_0}| \ .
\end{align*}
Since the Fourier coefficients $\widehat{(f')}_r$ of $f'$ are $ir \hat{f}_r$, one can inductively show $\sum\limits_{r \in \mathbb{Z}} |\widehat{f}_r \, r^k|^2 < \infty$ for all $f \in C^k$, which also implies $|\widehat{f}_r| \, |r|^{k+\frac{1}{2}} \xrightarrow{|r| \rightarrow \infty} 0$.
Edit: For $C^k$-functions one can even get something like the $O(\frac{1}{a^n})$-bound you were looking for, though the involved constant is dependent on $f$ and $n$. We make use of the fact that for all $f \in C^k$ there is a $c=c_{f}>0$ such that $|\hat{f}_r| \leq \frac{c}{|r|^{k+\frac{1}{2}}}$ for all $r \neq 0$.
We prove first that the equality
\begin{align*}
& 0 = r_0+ar_1+...+a^nr_n \tag{1}
\end{align*}
with all $r_j$ being different from zero must imply $|r_0| \cdots |r_n| > \frac{a^n}{n^n}$. As $r_n \neq 0$, we have $a^n \leq |r_0+ar_1+...+a^{n-1}r_{n-1}| \leq |r_0|+...+a^{n-1}|r_{n-1}|$. There must be some $j<n$ such that $|r_j| \geq a^{n-j} \frac{a}{n}$ or else we for $a \geq 3$ would have $|r_0|+...+a^{n-1}|r_{n-1}| < a^{n}$, which would be a contradiction. Set $n^{(2)} = j-1$ then either $r_0+ar_1+...+a^{n^{(2)}}r_{n^{(2)}}$ is zero, in which case we proceed as above to find some $j^{(2)}<n^{(2)}$ with $|r_{j^{(2)}}| \geq a^{n^{(2)}-j^{(2)}} \frac{a}{n}$, or $r_0+ar_1+...+a^{n^{(2)}}r_{n^{(2)}}$ is not zero, in which case (1) implies $a^{n^{(2)}+1} \leq |r_0+ar_1+...+a^{n^{(2)}}r_{n^{(2)}}|$ and we can proceed as above to find a $j^{(2)}\leq n^{(2)}$ with $|r_{j^{(2)}}| \geq a^{n^{(2)}-j^{(2)}} \frac{a}{n}$. We can repeat this induction until $j^{(K)}$ is $0$ and with the lower bounds for $|r_{j^{(\bullet)}}|$ arrive at $|r_0| \cdots |r_n| > \frac{a^n}{n^n}$.
We now sharpen our earlier bound on $I(f,a,n)$ into
\begin{align*}
& |I(f,a,n)| \leq \sum\limits_{\substack{r_0,...,r_n \in \mathbb{Z} \\ 0 = r_0+ar_1+...+a^nr_n}} \prod\limits_{j=0}^n |\hat{f}_{r_j}|\\
& \leq \sum\limits_{\substack{r_0,...,r_n \in \mathbb{Z}\setminus\{0\} \\ |r_0| \cdots |r_n| > \frac{a^n}{n^n}}} \prod\limits_{j=0}^n \frac{c}{|r|^{k+\frac{1}{2}}}\\
& = c^n \sum\limits_{p = \frac{a^n}{n^n}+1}^\infty \frac{1}{p^{k+\frac{1}{2}}} \sum\limits_{\substack{r_0,...,r_n \in \mathbb{Z}\setminus\{0\} \\ |r_0| \cdots |r_n| = p}} 1 \ .
\end{align*}
Since $l_p := \lceil\log_2(p)\rceil \leq C_\varepsilon p^\varepsilon$ is an upper bound for the number of prime factors the number $p$ can have, and there are ${l_p+n \choose n} \leq \frac{(l_p+n)^n}{n!} \leq C_{\varepsilon,n} p^{\varepsilon n}$ ways to distribute these primes among $n+1$ factors, we can further bound
\begin{align*}
& |I(f,a,n)| \leq c^n \sum\limits_{p = \frac{a^n}{n^n}+1}^\infty \frac{1}{p^{k+\frac{1}{2}}} 2^{n+1} {l_p+n \choose n}\\
& \leq 2(2c)^n C_{\varepsilon,n} \sum\limits_{p = \frac{a^n}{n^n}+1}^\infty \frac{p^{\varepsilon n}}{p^{k+\frac{1}{2}}} \leq \frac{2(2cC_\varepsilon)^n}{n!} \int_{\frac{a^n}{n^n}}^\infty x^{\varepsilon n-k-\frac{1}{2}} \, dx\\
& = 2(2c)^n C_{\varepsilon,n} \Big[ \frac{1}{\varepsilon n-k+\frac{1}{2}} x^{\varepsilon n-k+\frac{1}{2}} \Big]_{\frac{a^n}{n^n}}^\infty\\
& = 2(2c)^n C_{\varepsilon,n} \frac{1}{k-\varepsilon n-\frac{1}{2}} \Big( \frac{a^n}{n^n} \Big)^{\varepsilon n-k+\frac{1}{2}}\\
& \leq C(n,f,\varepsilon) a^{-n (k-\varepsilon n-\frac{1}{2})} \ .
\end{align*}
