How can zero point energy vacuum be Lorentz invariant? What distributions of electromagnetic oscillating fields are Lorentz invariant? How can oscillating electromagnetic fields look the same regardless of inertial frame?
According to Marshall and Boyer a cubic distribution like 
$ \rho(\omega) d\omega = \frac{\hbar\omega^3}{2 \pi^2c^3} d\omega \tag{1}$ 
is Lorentz invariant. How?
The cubic distribution is somewhat explained on WikiPedia page on zero-point energy.
Boyer just assumes that the spectrum is Lorentz invariant. This of course tells nothing. There might not be any such spectrum at all. 
Quotes from Boyer.
From the abstract:

The hypotheses assume (a) the existence, at the absolute zero of temperature, of classical
  homogeneous fluctuating electromagnetic radiation with a Lorentz-invariant spectrum;

Page  1375

Making no initial assumptions about the form of this energy beyond the Lorentz
  invariance of the spectrum...

Page 1377 

We have postulated that the zero-point spectrum is Lorentz-invariant

Marhsall makes similar assumptions.
From the abstract:

The random pure radiation field postulated in an earlier paper is set up in a relativistically invariant manner.

Page 538 

In effect, the requirement that the zero-temperature energy spectrum be left invariant under a Lorentz transformation ...

Page 545

The random field there postulated has here been defined in a relativistically invariant manner.

Where the previous paper is 
MARSHALL, T. W. Proc. Roy. Soc. Ser. A, 276 (1963), 475-491.
 A: 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$.
A: This isn't a complete answer, just some remarks about what should be invariant, if there really is something invariant here.  I use $c = 1$ to simplify things.
I'm expecting that the EM radiation energy-momentum tensor should be of the following shape (incoherent superposition of EM plane waves):
\begin{equation}\tag{1}
T_{ab} \propto \int k_a \, k_b \, f(\omega) \, d^3 k,
\end{equation}
where $k^a$ is a ligth-like 4-vector (so $k_a \, k^a = 0$).  The energy-momentum of EM radiation should have a null trace.  Expression (1) satisfy this constraint, and cannot gives something like a Lorentz invariant enery-momentum $T_{ab}^{\text{vac}} = u \, \eta_{ab}$ (which implies pressure $p = -\, u$).  The function $f(\omega)$ is arbitrary, but is defined in the special frame for which radiation is isotropic.  Here:
\begin{equation}\tag{2}
\omega = u_a^{\text{iso}} \, k^a
\end{equation}
is the frequency in this special frame and could be defined to be a Lorentz invariant (the frequency as measured by the observer of 4-velocity $\smash{u_{\text{iso}}^a}$, for which radiation is isotropic!).  It's easy to verify that pressure of the radiation is $p = \smash{\frac{1}{3}} \, u$, in this frame:
\begin{align}
u \equiv T_{00} &\propto \int \omega^2 \, f(\omega) \, d^3 k = 4 \pi \int \omega^4 \, f(\omega) \, d\omega. \tag{3} \\[12pt]
T_{0 i} &\propto \int \omega^2 \, n_i \, f(\omega) \, d^3 k = 0. \qquad \text{(from isotropy)} \tag{4} \\[12pt]
T_{ij} &\propto \int \omega^2 \, n_i \, n_j \, f(\omega) \, d^3 k = \int \omega^4 \, \big( \frac{4 \pi}{3} \, \delta_{ij} \big) \, f(\omega) \, d\omega \equiv \frac{1}{3} \, u \, \delta_{ij}. \tag{5}
\end{align}
Expression (5): $T_{ij} = \frac{1}{3} \, u \, \delta_{ij}$, implies that pressure is $p = \frac{1}{3} \, u$, as it should for random radiation, and that the full energy-momentum of incoherent EM radiation cannot be Lorentz invariant (i.e $T_{ab}^{\text{rad}} \ne u \, \eta_{ab}$).  In any other inertial frame (frame $\tilde{\mathcal{R}}$, for which $u_{\text{iso}}^a \ne \{ 1, 0, 0, 0 \}$), the radiation has a dipolar pattern since (2) gives
\begin{equation}\tag{6}
f(\omega) = f\big(\gamma \, (1 - v \cos{\vartheta}) \, \tilde{\omega} \, \big),
\end{equation}
which implies an anisotropy of radiation in that frame.  Then $\tilde{T}_{0i} \ne 0$ and there's a directional flux.
However, it doesn't mean that expression (3) cannot have a Lorentz invariant spectrum.  We recover
\begin{equation}\tag{7}
du \propto \rho(\omega) \, d\omega \propto \omega^3 \, d\omega
\end{equation}
in the special case $f(\omega) \propto \omega^{- 1}$ (which I think is weird, since it implies some infrared divergence for $f(\omega)$, and ultraviolet divergence for the total energy density $u$!).  This isn't enough to show that this special spectrum of $u$ is Lorentz invariant.  This is why my answer is currently incomplete.
Notice that, according to expression (1), the product $f(\omega) \, d^3 k$ is a Lorentz scalar, while $f(\omega)$ and $d^3 k$ separately are not Lorentz scalars.  Also, it can be proved that $d^3 k / \omega$ is a Lorentz invariant measure in the momentum space of the radiation.  Then $f(\omega) \, \omega$ should be a Lorentz scalar too.  If $f(\omega) \propto \omega^{-1}$, then obviously $f(\omega) \, \omega = \text{cste}$!  
I feel personally unsure and puzzled by the previous paragraph.
A: Assuming that the vacuum energy consists of radiation that is Doppler shifted when moving relative to leads to the spectral distribution not being Lorentz invariant. The longitudinal Doppler shift
$$ \omega'  {\displaystyle ={\sqrt {\frac {1-\frac{v}{c} }{1+\frac{v}{c} }}}\,\omega} $$
leads the spectrum to become
$$\rho(\omega') d\omega' = {\left(\frac {1-\frac{v}{c}}{1+\frac{v}{c} }\right)}^{\frac{3}{2}} \frac{\hbar\omega^3}{2 \pi^2c^3} d\omega $$
This is not Lorentz invariant. It differs with a factor ${\left(\frac {1-\frac{v}{c}}{1+\frac{v}{c} }\right)}^{\frac{3}{2}}$. Serious consequences are that vacuum would have drag due to differences in radiation pressure and that there is only one inertial system.
This result is commonly avoided by postulating Lorentz invariance but what kind of fields could vacuum consist of if not radiation fields that can be Doppler shifted?
