Edit: I realized that, in your question, the Fourier series decomposition of the temperature should go as (please confirm the braces after $n$)
\begin{equation}
T(t) = T_0 + \sum_{n=1}^\infty T_n e^{in (\omega t - \phi)}
\end{equation}
Original answer: The only equation that contains any physical link between $T(t)$ and $S(t)$ is the one you have described:
\begin{equation}
\sigma T^4(t) = (1-A)S(t).
\end{equation}
Taking the time derivative and then dividing both sides by $\sigma T^4(t)$, you get
\begin{align}
4\frac{\sigma T^3(t)}{\sigma T^4(t)} \frac{dT(t)}{dt} &= \frac{(1-A)}{\sigma T^4(t)} \frac{dS(t)}{dt}\\
& \\
\Rightarrow S(t) \frac{d T(t)}{dt} &= \frac{1}{4} T(t) \frac{dS(t)}{dt},
\end{align}
where in the second equation, I have replaced $\sigma T^4(t)$ using the above equation and simplified further.
Here, you can substitute the respective Fourier expansions and obtain
\begin{align}
S_0\left ( 1 + \sum_{m=1}^\infty \frac{S_m}{S_0} e^{im\omega t} \right )
\left( \sum_{n=1}^\infty T_n (i n \omega) e^{in(\omega t - \phi)} \right ) =
\frac{1}{4}T_0
\left( 1 + \sum_{n=1}^\infty \frac{T_n}{T_0} e^{i n(\omega t - \phi)}\right) \\
\times \left( \sum_{m=1}^\infty S_m (i m\omega)e^{i m \omega t}\right )
\end{align}
\begin{align}
\Rightarrow \left ( 1 + \sum_{m=1}^\infty \frac{S_m}{S_0} e^{im\omega t} \right )
\left( \sum_{n=1}^\infty \frac{T_n}{T_0} (i n \omega) e^{in(\omega t - \phi)} \right ) = \frac{1}{4}
\left( 1 + \sum_{n=1}^\infty \frac{T_n}{T_0} e^{i n(\omega t - \phi)}\right) \\
\times \left( \sum_{m=1}^\infty \frac{S_m}{S_0} (i m \omega) e^{i m \omega t}\right )
\end{align}
which, upon re-ordering the terms and comparing the coefficients for each unique Fourier mode, should give you the required result.
EDIT: The above expression can be re-written as
\begin{equation}
\sum_{m,n=0}^\infty \frac{S_m}{S_0} \frac{T_n}{T_0} (i n \omega) e^{i(m+n)\omega t -in\phi} = \frac{1}{4}
\sum_{m,n=0}^\infty \frac{S_m}{S_0} \frac{T_n}{T_0} (i m \omega) e^{i(m+n)\omega t -in\phi}
\end{equation}
On the RHS, exchange $m$ with $n$ (they are merely dummy indices) and bring the expression to the LHS, you get
\begin{equation}
\sum_{m,n=0}^\infty \left ( \frac{S_m}{S_0} \frac{T_n}{T_0} e^{-in\phi} - \frac{1}{4} \frac{S_n}{S_0}\frac{T_m}{T_0} e^{-im\phi} \right ) (i n \omega) e^{i(m+n)\omega t} = 0
\end{equation}
You can further simplify by the following set of transformations
\begin{align}
k &= m + n\\
m &= k-n
\end{align}
and by requiring that the coefficient of $e^{ik\omega t}$ is ZERO $\forall k$, you can obtain the required relation. In fact, it will become
\begin{equation}
\sum_{k=0}^\infty \left ( \sum_{n = 0}^k i n \omega \left [
\frac{S_{k-n}}{S_0} \frac{T_n}{T_0} e^{-in\phi} - \frac{1}{4} \frac{S_n}{S_0}\frac{T_{k-n}}{T_0} e^{-i(k-n)\phi}
\right ]\right) e^{i k \omega t} = 0
\end{equation}
\begin{equation}
\Rightarrow \sum_{n = 0}^k i n \omega \left [
\frac{S_{k-n}}{S_0} \frac{T_n}{T_0} e^{-in\phi} - \frac{1}{4} \frac{S_n}{S_0}\frac{T_{k-n}}{T_0} e^{-i(k-n)\phi}
\right ] = 0
\end{equation}
for every value of $k$. Beyond this point, it would be more convenient to put in the values and build a recursion relation. For instance, $k = 1$ gives
$$
\frac{T_1}{T_0} = \frac{1}{4} \frac{S_1}{S_0} e^{i \phi}
$$
Take just the real part (or rather let the Fourier series index run from $-\infty$ to $+\infty$) and you get your expression (in terms of the $cos \phi$ which is equivalent to $1/\sqrt{\dots}$). Using the relation for $n=1$, proceed to $k = 2$ and so on. You'll start seeing a pattern in the recursive relations.
