We study the scaling asymptotics of the eigenspace projection kernels $Π_{\hbar, E}(x,y)$ of the isotropic Harmonic Oscillator $- \hbar ^2 Δ+ |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 $Π_{\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^α$ with $0< α< \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}}$.