Consider a sinusoidal input x(t) = x_m \cos(\omega t) u(t) that produces a steady state output y_{ss}(t) and a network that has the transfer function H(s). Then the Laplace transform of the output is: Y(s) = H(s) L \left ( x_m \cos(\omega t) u(t) \right ) = H(s) L \left ( \left ( \frac{x_m}{2} e^{j \omega t} + \frac{x_m}{2} e^{-j \omega t}\right )u(t) \right )
So
Y(s) = \frac{x_m}{2} H(s) \frac{1}{s - j \omega} + \frac{x_m}{2} H(s) \frac{1}{s + j \omega}.
I was taught this in the context of electrical systems. My professor then went on to say that the poles of H(s) influence the transient solution, while the poles of \frac{1}{s \pm j \omega} make up the steady state solution. This feels kind of hand-wavy, but it makes sense when you consider that the transient solution is influenced by the particular system, but the steady state solution (long term) is obviously influenced by the driving force (for instance, consider pushing a spring with a sinusoidally varying force). In this case, if the driving force is a sinusoid, then the laplace transform of the driving force will look like the terms above. So we have Y(s) = \frac{k_1}{s - j \omega} + \frac{k_2}{s + j \omega} + (terms \ involving \ poles \ of \ H(s)) which implies that k_1 = Y(s)(s - j \omega) evaluated at s = j \omega, and by extension, k_2 = k_1^*, or the complex conjugate of k_1. So that means k_1 = \frac{x_m}{2} H(j \omega), k_2 = k_{1}^*.
Finally, you have to remember this little shortcut when taking inverse Laplace transforms: if Y(s) = \frac{k_1}{s + a - i b} + \frac{k_2}{s + a + i b}, then y(t) = 2 |k_1| e^{-a t} \cos(b t + \theta). In this case, a = 0, and \theta is the angle of the complex number k_1.
Putting those things together, when we take the inverse Laplace transform of the equation above, we get y(t) = x_m |H(j \omega)| \cos(\omega t + \theta) + (transient \ terms)
Or, ignoring the transient terms, y_{ss}(t) = x_m |H(j \omega)| \cos(\omega t + \theta) where \theta is the angle of H(j \omega)
