It's just a big trig, sinusoidal, Fourier series thing:
$$\begin{align} y(t) &= \sum_{k=0}^{K} a_k \big( A \cos(\omega t) \big)^k \\ \\ &= \sum_{n=0}^{K} b_n \cos(n \omega t) \\ \end{align}$$
The result I got (I think I fixed my error with the limits in swapping summations) is:
apparently, the DC component ($n=0$) of the Fourier cosine series is:
$$ b_0 = a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \frac{(2k)!}{(k!)^2} \left(\frac{A}{2}\right)^{2k} $$
The coefficient for the even harmonic ($n$ even) terms:
$$ b_{n} = 2 \sum_{k=\frac{n}{2}}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \frac{(2k)!}{\left(k+\frac{n}{2}\right)! \left(k-\frac{n}{2}\right)!} \left(\frac{A}{2}\right)^{2k} $$
And the coefficient for the odd harmonic ($n$ odd) terms:
$$ b_{n} = 2 \sum_{k=\frac{n-1}{2}}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} a_{2k+1} \frac{(2k+1)!}{(k+\frac{n+1}{2})!(k-\frac{n-1}{2})!} \left(\frac{A}{2}\right)^{2k+1} $$
where $\lfloor u \rfloor$ is the floor() function, which is the largest integer no larger than the argument $u$.
$$ \lfloor u \rfloor \le u < \lfloor u \rfloor + 1 \qquad \qquad \lfloor u \rfloor \in \mathbb{Z}, u \in \mathbb{R}$$ or $$ u-1 < \lfloor u \rfloor \le u \qquad \qquad \lfloor u \rfloor \in \mathbb{Z}, u \in \mathbb{R}$$
Since the OP wants a reference, I would say the use the sums found in Prudnikov, Integrals and Series, vol. 1. 4.4.1.16 and 4.4.1.17. Ultimately they are not any more profound than the binomial theorem. In this particular situation they are, with a little algebra to clean them up
$$ \cos^{2n}(x)=2^{1-2n}\sum_{k=1}^{n}\binom{2n}{n-k}\cos{(2\,k\,x)} + 2^{-2n}\binom{2n}{n}$$ $$\cos^{2n+1}(x)=2^{-2n}\sum_{k=0}^{n}\binom{2n+1}{n-k}\cos{((2\,k+1)x)} $$
The OP's procedure basically amounts to inserting these expressions in the original sum and interchanging the summations.