We calculate the CMB $\mu$-distortion and the angular power spectrum of its
cross-correlation with the temperature anisotropy in the presence of the
non-Gaussian neutrino isocurvature density (NID) mode. While the pure Gaussian
NID perturbations give merely subdominant contribution to $<\mu>$ and vanishing
$< \mu T>$, the latter quantity can be large enough to be detected in the
future when the NID perturbations $\mathcal S(\mathbf x)$ are proportional to
the square of a Gaussian field $g(\mathbf x)$, i.e. $\mathcal S({\mathbf
x})\propto g^2({\mathbf x})$. In particular, large $< \mu T>$ can be realized
since Gaussian-squared perturbations can yield a relatively large bispectrum,
satisfying the constraints from the power spectrum of CMB anisotropies, i.e.
$\mathcal{P}_\mathcal{SS}(k_0) \sim\mathcal{P}_g^2(k_0)\lesssim10^{-10}$ at
$k_0=0.05$ Mpc$^{-1}$. We also forecast constraints from the CMB temperature
and E-mode polarisation bispectra, and show that
$\mathcal{P}_g(k_0)\lesssim10^{-5}$ would be allowed from Planck data. We find
that $< \mu >$ and $|l(l+1)C^{\mu T}_l|$ can respectively be as large as
$10^{-9}$ and $10^{-14}$ with uncorrelated scale-invariant NID perturbations
for $\mathcal{P}_g(k_0)=10^{-5}$. When the spectrum of the Gaussian field is
blue-tilted (with spectral index $n_g \simeq 1.5$), $< \mu T>$ can be enhanced
by an order of magnitude.