I will construct a response in the only way I see possible, given the concerns of access to the source as mentioned in the comments. I will provide an in detail outline of how photon spectra transform under LTs, and try to directly address the question concerning the Lorentz invariance of the photon spectrum.
In general a photon spectrum ${dN_\gamma \over dE_\gamma d\Omega_\gamma}$ is not a lorentz invariant (as I will explain further below). However, there is an invariant that is close to this, which may be what the authors are getting at.
In particle physics, if we want to know the photon spectrum in the lab frame ${dN_\gamma \over dE_\gamma d\Omega_\gamma}$ , what we do is use the fact that the two observers must agree on the number of photons emitted $dN_\gamma$ in order to transform the (usually easier to find) distribution in the rest frame of, say, a decaying particle $S'$, back to to the lab frame.
So how doest this work? Well, suppose that the photon spectrum is given by function ${dN_\gamma' \over dE'_\gamma d\Omega'_\gamma}$ in the parent particle's rest frame, and suppose that the spectrum in the lab frame is given by ${dN_\gamma \over dE_\gamma d\Omega_\gamma}$. Using the fact that $dN_\gamma$ is an invariant (i.e. $dN_\gamma = dN'_\gamma$), it follows then that
$$dN_\gamma = dN_\gamma' \implies {dN_\gamma \over dE'_\gamma d\Omega'_\gamma} dE_\gamma' d\Omega_\gamma' = {dN_\gamma \over dE_\gamma d\Omega_\gamma} dE_\gamma d\Omega_\gamma$$
Now this could be precisely what the authors meant by Lorentz invariant spectrum. However, this wouldn't be very useful if we didn't know how to get from one distribution to the other. Luckily--given that $dN_\gamma = dN_\gamma'$--we do, as is beautifully outlined here in Jacobians and Relativistic Kinematics. The summary is that the Lorentz transformation of a photon spectrum is just a coordinate transformation subject to the constraint that $dN_\gamma = dN_\gamma'$. Then, we can relate the two distributions via the jacobian that they pick up. That is, (please allow me to drop the subscripts $\gamma$'s) we have that
$$\boxed{ \frac{dN}{dE d\Omega}= \left\lvert \frac{\partial(E', \Omega')}{\partial(E, \Omega)} \right\rvert \frac{dN}{dE' d\Omega'}}$$
This relation is exactly how the two spectra are related, and is of immense usefulness (as it has been to me when calculating photon spectra for dark matter annihilation spectra).
At any rate: Based on what I have just detailed, it seems to me that by "Lorentz invariant photon" spectrum is somewhat of a misnomer, it is not the distribution ${dN_\gamma \over dE_\gamma d\Omega_\gamma}$ that is Lorentz invariant (as it picks up a jacobian), but rather it is the number of photons that is invariant considered in some interval of phase space.
Edit: Using the Lorentz Invariant Phase Space element (LIPS). Note if instead we choose to parametrize our distribution with momentum instead of energy, then all of the previous arguments are valid if we replace $dE_\gamma d\Omega_\gamma \to d^3 k$ (same for the primes). This induces the substitution
$$ {dN_\gamma \over dE_\gamma d\Omega_\gamma} \to {dN_\gamma \over d^3k} $$
and then we may reuse the arguments above. Hence, our photon distributions now are just
\begin{align*}
\begin{cases}
& {dN_\gamma \over d^3 k} \qquad \text{(Lab Frame)}\\
& {dN_\gamma \over d^3 k'} \qquad \text{(Rest Frame of Parent Particle)}
\end{cases}
\end{align*}
It therefore follows that the invariant quantities are now
$$dN_\gamma = dN_\gamma' \implies {dN_\gamma \over d^3 k'} d^3k' = {dN_\gamma \over d^3 k} d^3k$$
where of course the phase space elements $d^3 k \neq d^3 k'$ and hence are not in general Lorentz invariant. Now, this is inconvenient so in field theory what we do is define the Lorentz Invariant Phase Space (LIPS) element by
$$ (d\Pi)_{LIPS} \equiv \prod_{i=1}^n \frac{d^3k_i}{2E_{k_i}}$$
for $n$ final state particles, and $E_{k_i} = \sqrt{\mathbf{k}_i^2 + m_i^2}$. It can be shown, as it is in this answer,that this indeed is a lorentz invariant quantity. Therefore, using this as the phase space measure, we have the relation that
$$dN_\gamma = dN_\gamma' \implies {dN_\gamma \over d^3 k'} \frac{d^3k'}{2E_{k'}} = {dN_\gamma \over d^3 k} \frac{d^3k}{2E_{k}}$$
for a single particle final state (i.e. $n=1$). What this means is that the two observers must agree on the quantities $ \frac{d^3k'}{2E_{k'}}$ and $\frac{d^3k}{2E_{k}}$. However, this does do not imply that need to agree on the quantities ${dN_\gamma \over d^3 k'}$ and ${dN_\gamma \over d^3 k}$. But from particle conservation, they still must agree on the product
$$dN_\gamma =\overbrace{ \underbrace{{dN_\gamma \over d^3 k} }_{\text{Not Lorentz Invariant}}\frac{d^3k}{2E_{k}}}^{\text{Lorentz Invariant}}. $$
This is what is meant by "Lorentz Invariant Distribution". It would seem that if they agree on the phase space element and on the product of the phase space then surely they must agree on ${dN_\gamma \over d^3 k}$, right? Well, no. Two different integrands can yield the the same number.
Edit 2: As @Cham pointed out in the comments, given the above facts, we can then consider your distribution. Namely, that
$$ \hbar \omega \underbrace{\frac{1}{c^3}\omega^2 \ d\omega d\Omega}_{d^3k} = \hbar\omega \ d^3 k $$
So we see that
$$\frac{dN_\gamma}{d^3 k} = \hbar \omega$$
Also, our LIPS element is given by
$$ \frac{d^3k}{2\omega}$$
which means that
$$dN_\gamma =\frac{d^3k}{2\omega} \frac{dN_\gamma}{d\omega d\Omega_\gamma} = \hbar\omega \frac{d^3k}{2\omega} = d^3 k \frac{\hbar}{2}. $$
So using the conservation of particle number we have that
$$ \hbar\omega \frac{d^3k}{2\omega} = \hbar\omega' \frac{d^3k'}{2\omega'} \implies \hbar\ d^3k = \hbar\ d^3k'.$$
So now we can see that when using the Lorentz Invariant Phase Space that the "distribution" $\propto \hbar$ is invariant. I use quotation marks because technically the distribution is w.r.t the phase space element and so really is linear in $\omega$.
