  We study the scaling asymptotics of the eigenspace projection kernels
$\Pi_{\hbar, E}(x,y)$ of the isotropic Harmonic Oscillator $- \hbar ^2 \Delta +
|x|^2$ of eigenvalue $E = \hbar(N + \frac{d}{2})$ in the semi-classical limit
$\hbar \to 0$. The principal result is an explicit formula for the scaling
asymptotics of $\Pi_{\hbar, E}(x,y)$ for $x,y$ in a $\hbar^{2/3}$ neighborhood
of the caustic $\mathcal C_E$ as $\hbar \to 0.$ The scaling asymptotics are
applied to the distribution of nodal sets of Gaussian random eigenfunctions
around the caustic as $\hbar \to 0$. In previous work we proved that the
density of zeros of Gaussian random eigenfunctions of $\hat{H}_{\hbar}$ have
different orders in the Planck constant $\hbar$ in the allowed and forbidden
regions: In the allowed region the density is of order $\hbar^{-1}$ while it is
$\hbar^{-1/2}$ in the forbidden region. Our main result on nodal sets is that
the density of zeros is of order $\hbar^{-\frac{2}{3}}$ in an
$\hbar^{\frac{2}{3}}$-tube around the caustic. This tube radius is the
`critical radius'. For annuli of larger inner and outer radii $\hbar^{\alpha}$
with $0< \alpha < \frac{2}{3}$ we obtain density results which interpolate
between this critical radius result and our prior ones in the allowed and
forbidden region. We also show that the Hausdorff $(d-2)$-dimensional measure
of the intersection of the nodal set with the caustic is of order $\hbar^{-
\frac{2}{3}}$.
