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