$\newcommand{\bra}[1]{\langle #1 |}$ $\newcommand{\ket}[1]{| #1 \rangle}$ $\newcommand{\braket}[2]{\langle #1 | #2 \rangle}$ $\newcommand{\bbraket}[3]{\langle #1 | #2 | #3 \rangle}$ Although the question asks specifically about a harmonic oscillator, we can understand the meaning of the spectral density by considering a somewhat more general problem.
Consider a quantum system $S$ coupled to an environment $E$.
Focus attention on a single pair of states of $S$ which we denote $\ket{0}$ and $\ket{1}$.
Denote the energy difference between these levels as $E_{10} = E_1 - E_0$ and define the transition frequency by
$$\omega_{10} \equiv E_{10}/\hbar \, .$$
With this notation, the Hamiltonian of $S$ is just
$$H_S = - \frac{\hbar \omega_{10}}{2} \sigma_z \, .$$
Suppose $E$ couples to the $S$ via an interaction Hamiltonian given by
$$ V = g (F \otimes \sigma_x)$$
where here $F$ is an operator acting on $E$ and $\sigma_x$ is the Pauli $x$ operator on the two-level subspace of $S$ consisting of $\ket{0}$ and $\ket{1}$.
We use perturbation theory to compute the affect of this interaction on $S$.
There are several ways to think about perturbation theory.
I prefer to use the rotating frame in which the system and environment self Hamiltonians go to zero and the interaction Hamiltonian becomes
$$V(t) = g ( F(t) \otimes (\cos(\omega_{10} t )\sigma_x + \sin(\omega_{10} t) \sigma_y ) )$$
In order to compute the evolution of $S$ we use first order perturbation theory.
We start with Schrodinger's equation
$$ i \hbar (d/dt)\ket{\Psi(t)} = V(t) \ket{\Psi(t)}$$
which is formally solved by
$$ \ket{\Psi(t)} = \ket{\Psi(0)} - \frac{i}{\hbar} \int_0^t dt' \, V(t') \ket{\Psi(t')} \, . $$
Iterating this equation once gives
$$ \ket{\Psi(t)} = \ket{\Psi(0)} - \frac{i}{\hbar} \int_0^t dt' \, V(t') \ket{\Psi(0)} \, . $$
Suppose $S$ starts in the lower state $\ket{0}$ and $E$ starts in some particular state $\ket{n}$.
In other words, $\ket{\Psi(0)} = \ket{n} \otimes \ket{0}$.
Let us compute the amplitude for $S$ to make a transition to $\ket{1}$:
\begin{align}
\braket{1}{\Psi(t)} =
&= - \frac{i}{\hbar} \int_0^t \bbraket{1}{\tilde{V}(t')}{\Psi(0)} \, dt' \\
&= \frac{-ig}{\hbar} \int_0^t F(t')\ket{n} \otimes \bbraket{1}{e^{i \omega_{10} t'} \sigma_+ + e^{-i \omega_{10} t'} \sigma_-}{0} \, dt' \\
&= \frac{-ig}{\hbar} \int_0^t F(t')\ket{n} e^{i \omega_{10} t'} \, dt' \, . \qquad (1)
\end{align}
Note that the expression on the left hand side is a state vector for the environment.
This is because we have only taken the inner product with respect to a final state of $S$, but have not specified a final state of $E$.
Note also that our initial state $\ket{\Psi(0)}$ had the environment in a specific state $\ket{n}$.
Let's consider a more interesting case where the environment is initially in a mixed state $\rho_E = \sum_n \rho_{nn} \ket{n}\bra{n}$.
To accommodate this into our analysis we form the density matrix for the environment state given in equation (1) and sum it over the initial states of the environment:
\begin{equation}
\rho(t) = \left(\frac{g}{\hbar}\right)^2 \sum_n \rho_{nn} \int_0^t \int_0^t \, dt' \, dt'' \, F(t')\ket{n}\bra{n}F(t'') e^{i \omega_{10} (t'-t'')} \, .
\end{equation}
Now, to finally get a probability $p_1(t)$ for the system to be in $\ket{1}$ averaged over all final states of the environment, we trace over the environmental states:
\begin{align}
p_1(t)
&= \left(\frac{g}{\hbar}\right)^2 \sum_{nm} \int_0^t \int_0^t \, dt' \, dt'' \, \rho_{nn} \bbraket{m}{F(t')}{n} \bbraket{n}{F(t'')}{m} e^{i \omega_{10} (t' - t'')} \\
&= \left(\frac{g}{\hbar}\right)^2 \sum_{nm} \int_0^t \int_0^t \, dt' \, dt'' \, \rho_{nn} \bbraket{n}{F(t'')}{m} \bbraket{m}{F(t')}{n} e^{i \omega_{10} (t' - t'')} \\
&= \left(\frac{g}{\hbar}\right)^2 \sum_n \int_0^t \int_0^t \, dt' \, dt'' \, \rho_{nn} \bbraket{n}{F(t'') F(t')}{n} e^{i \omega_{10} (t' - t'')} \\
&= \left(\frac{g}{\hbar}\right)^2 \int_0^t \int_0^t \, dt' \, dt'' \, \langle F(t'') F(t') \rangle e^{i \omega_{10} (t' - t'')} \, .
\end{align}
Changing variables to $\tau \equiv t'' - t'$ and assuming the average of $F$ is stationary gives
\begin{equation}
p_1(t)
= \left(\frac{g}{\hbar}\right)^2 \int_0^t \int_{-t'}^{t-t'} \, dt' \, d\tau \, \langle F(\tau) F(0) \rangle e^{-i \omega_{10} \tau} \, .
\end{equation}
If the correlation function of $F$ is very short then we can extend the limits of integration of $\tau$ to positive and negative infinity.
Doing this and defining the spectral density as
\begin{equation}
S_{FF}(\Omega) \equiv \int_{-\infty}^\infty \, dt \, \langle F(t) F(0) \rangle e^{i \Omega t}
\end{equation}
we find
\begin{equation}
p_1(t) = \left( \frac{g}{\hbar} \right)^2 \int_0^t S_{FF}(-\omega_{10})
\end{equation}
which we can interpret as an upward transition rate
\begin{equation}
\Gamma_{\uparrow} = \frac{dp_1}{dt} = \left(\frac{g}{\hbar}\right)^2 S_{FF}(-\Omega) \, .
\end{equation}
Similar arguments can be used to show
\begin{equation}
\Gamma_\downarrow = \left(\frac{g}{\hbar}\right)^2 S_{FF}(\Omega) \, .
\end{equation}
Therefore, we have shown that when a system is connected to a quantum environment via an operator $F$, the quantum spectral density $S_{FF}$ characterizes the upward and downward transition rates induced on the system by the environment.
We note that positive frequencies in the spectral density give downward transitions in the system (absorption by the environment) while negative frequencies give upward transitions (emission by the environment).
Bonus:
For those readers interested in quantum electronics, the quantum spectral density of current for an arbitrary linear element with admittance $Y(\omega)$ which is in thermal equilibrium is
$$S_{II}(\omega) = \hbar \omega \Re Y(\omega) \left[ \coth \left( \frac{\beta \hbar \omega}{2}\right) +1 \right]$$
where $\beta \equiv 1 / k_b T$.
