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 \langle \mu \rangle and vanishing \langle \mu T \rangle , the latter quantity can be large enough to be detected in the future when the NID perturbations \mathcal { S } ( \bm { x } ) are proportional to the square of a Gaussian field g ( \bm { x } ) , i.e . \mathcal { S } ( { \bm { x } } ) \propto g ^ { 2 } ( { \bm { x } } ) . In particular , large \langle \mu T \rangle 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 } ) \raisebox { -2.15 % pt } { $\ > \stackrel { \textstyle < } { \sim } \ > $ } 10 ^ { -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 } ) \raisebox { -2.15 pt } { $\ > \stackrel { \textstyle < } { \sim } \ > $ } 10 % ^ { -5 } would be allowed from Planck data . We find that \langle \mu \rangle 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 ) , \langle \mu T \rangle can be enhanced by an order of magnitude .