The $3$rd equation line, where you introduced the variable $e$, is not necessary and it could be removed without affecting the argument. Also note the term $\zeta_n$ in your question is not a primitive $n$th root of unity, but just a general $n$th root of unity.
As you correctly state we need to show the double summation in the $2$nd equation line equals $\sum_{j=1}^n\cos \frac{2j\pi}{n}$, and that to do that it suffices to show that the two sets $S_1$ and $S_2$ are equal. The $S_1 \subseteq S_2$ part is straightforward, however the other direction is not so obvious, but it can be arrived at by making use of the totient function $\phi$.
Instead of showing $S_2 \subseteq S_1$ to complete the proof of $S_1 = S_2$, since we are dealing with finite sets we can instead show $|S_1| = |S_2|$, ie. $|S_1| = n$, since we already know $S_1 \subseteq S_2$.
So, using a slightly different definition from yours, for a fixed $n \geq 1$ let :
\begin{eqnarray*}
S_1 & = & \{\textstyle{\frac{k}{d}} : d \mid n, k \in [1, d], (k, d) = 1 \}, \\
S_2 & = & \{ \textstyle{\frac{j}{n}} : j \in [1, n] \}
\end{eqnarray*}
We can write $S_1$ as a union of the subsets $S_1^{(d)}$ that come from particular divisors $d$ of $n$ :
$$
S_1 = \bigcup_{d \mid n} S_1^{(d)}
$$
where, with $n \geq 1$ fixed and $d$ a divisor of $n$ :
$$
S_1^{(d)} = \{ \textstyle{\frac{k}{d}} : k \in [1, d], (k, d) = 1 \}.
$$
Then we have $|S_1^{(d)}| = \phi(d) \; \forall d$, and so if the $S_1^{(d)}$'s are mutually disjoint, forming a partition of $S_1$, then using the divisor sum of the totient (eg here or [1, p26, Theorem 2.2]) :
$$
|S_1| = \sum_{d \mid n} |S_1^{(d)}| = \sum_{d \mid n} \phi(d) = n = |S_2|
$$
as required. Thus it suffices to show the $S_1^{(d)}$'s are mutually disjoint. Consider an $x \in S_1^{(d_1)} \cap S_1^{(d_2)}$. Then there exists $k_1 \in [1, d_1]$ with $(k_1, d_1) = 1$, and $k_2 \in [1, d_2]$ with $(k_2, d_2) = 1$ such that :
$$
x = \frac{k_1}{d_1} = \frac{k_2}{d_2}.
$$
But these are fractions in their simplest form so their numerators and denominators must be the same, due to the uniqueness of the simplest form of a fraction, and hence in particular $d_1 = d_2$. Thus the $S_1^{(d)}$'s are mutually disjoint. Alternatively we could say $k_1 d_2 = k_2 d_1$ and then in terms of the 'multisets' which are the prime decompositions, the primes of $d_2$ must all be contained within $d_1$, as $d_2$ shares no primes with $k_2$, so that '$d_2 \subseteq d_1$', and likewise as $d_1$ shares no primes with $k_1$, so '$d_1 \subseteq d_2$', and thus $d_1 = d_2$ - which is how the uniqueness of the simplest form is proved.
The same method of proving $S_1 = S_2$ can be used to solve a more general problem (see [1, Problem 2.14, p48]), namely if :
\begin{eqnarray*}
F(n) & = & \sum_{k = 1}^n f(k/n) \\
\mbox{and} \hspace{2em} F^*(n) & = & \sum_{ \substack{k = 1 \\ (k, n) = 1} }^n f(k/n)
\end{eqnarray*}
then $F$ and $F^*$ are related by $F^* = \mu * F$.
From this we can arrive at the following relation for the sum of the primitive $n$th roots of unity (from which the relation you require can be deduced by taking the real part) :
$$
\sum_{ \substack{k = 1 \\ (k, n) = 1} }^n e^{2 \pi i k / n} = \mu(n), \hspace{3em} \forall n \geq 1
$$
by putting $f(x) = e^{2 \pi i x}$, so that the LHS is $F^*(n)$ and :
$$
F(n) = \sum_{k = 1}^n e^{2 \pi i k / n} \hspace{1em} = \mbox{sum of $n^{\mathrm{th}}$ roots of unity.}
$$
But for $n > 1$, the sum of the $n$th roots of unity $\{ \alpha_r \}$ is minus the coefficient of $z^{n - 1}$ in the complex polynomial $z^n - 1 = (z - \alpha_1) \cdots (z - \alpha_n)$, thus :
$$
F(n) = \left\{
\begin{array}{ll}
1, & \mbox{when $n = 1$}, \\
0, & \mbox{when $n > 1$}
\end{array}
\right.
$$
which equals the identity $\delta$ in the Dirichlet ring, and hence $F^* = \mu$, as required.
References
[1] Tom M. Apostol (1976), Introduction to Analytic Number Theory, Springer-Verlag.
