The state space form of the harmonic oscillator consists of the following system of differential equations $$\begin{array}{rcl} \dot{q_1}(t) &=& - \omega q_2(t) \\ \dot{q_2}(t) &=& \omega q_1(t) \\ \end{array}$$ This can be solved by Laplace-transforming the differential equations to a set of equations $$\begin{array}{rcl} s Q_1(s) - q_1(0) &=& - \omega Q_2(s) \\ s Q_2(s) - q_2(0) &=& \omega Q_1(s) \\ \end{array}$$ Solving this for $Q_1(s), Q_2(s)$ gives a solution in the Laplace domain, which can be converted back to the time domain. The solutions looks like $$Q_i(s) = \frac{n_i(s)}{s^2 + \omega^2}$$ where the numerators $n_i(q)$ are first order polynomials in terms of the initial conditions $q_i(0)$. The interesting part is the denominator $s^2 + \omega^2 = (s-j\omega)(s+j\omega)$ which determines the \emph{shape} of the result. Splitting the $Q_i(s)$ in partial fractions then gives the sum of first order rational functions proportional to $(s-j\omega)^{-1}$ and $(s+j\omega)^{-1}$, which are Laplace transforms of $e^{j\omega t}$ and $e^{-j\omega t}$ respectively. The denominator can also be obtained as $\det(Is-A)$ where $A$ is the system matrix $$A = \matrix{0 & -\omega \\ \omega & 0}.$$