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}$$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.