Use the moebius inversion formula to obtain $a_n = \sum_{d|n} \mu (d) 2^{\frac{n}{d}}$. For an arbitrary prime $p$, let $p^k$ be the highest power of the prime $p$ dividing $n=p^k m$. Then $\sum_{d|n} \mu (d) 2^{\frac{n}{d}} = \sum_{p\nmid d|n} \mu (d) 2^{\frac{n}{d}} + \sum_{p|d|n} \mu (d) 2^{\frac{n}{d}} = \sum_{p\nmid d|n} \mu (d)( 2^{\frac{mp^k}{d}} - 2^{\frac{mp^{k-1}}{d}})$. But $2^{\frac{mp^k}{d}} - 2^{\frac{mp^{k-1}}{d}}= 2^{\frac{mp^{k-1}}{d}} \cdot \left( (2^m)^{\frac{p^{k-1}(p-1)}{d}} -1\right)$ is divisible by $p^k$ by the theorem of Euler-Fermat (in case of $p=2$, we still have $p^k=2^k|2^{2^{k-1}} |2^{\frac{mp^{k-1}}{d}}$), so the whole sum is divisible by $p^k$. But since $p$ was arbtrary, the sum is divisible by $n$.

Alternatively, let $b_n$ be the number of aperiodic sequences $\left(x_1,x_2,\ldots,x_n\right) \in \{1,2\}^n$ of length $n$ (i.e., there exists no $d \mid n$ such that $x_k = x_l$ whenever $k \equiv l \pmod{d}$). It is obvious that $n \mid b_n$ for every $n$ since the rotation \[\left(x_1,x_2,\ldots,x_{n-1},x_n\right) \mapsto \left(x_n,x_1,x_2,\ldots,x_{n-1}\right)\] preserves aperiodicity. The rest is showing that $a_n=b_n$ by observing that $\left\{b_n\right\}$ obeys the same identity that defines $\left\{a_n\right\}$. Indeed, the integer $2$ in the definition may be replaced by any other integer and the result still holds.