GHZ entangled numerical shadow of a matrix $A$ of dimension $d$ is defined as a probability ditribution $P_A(z)$ on the complex plane, supported on the maximally entangled numerical range $W^\mathrm{ent}(A)$. $$ P_A(z) := \int_{\Omega} {\rm d} \mu(\psi) \delta\Bigl( z-\langle \psi|A|\psi\rangle\Bigr), $$ where $\mu(\psi)$ denotes the unique unitarily invariant (Fubini-Study) measure on the set $$ \Omega=\{\ket{\psi} \in \mathbb{C}^d: \ket{\psi} = \frac{1}{\sqrt{2}} \bigotimes_{i=1}^d U_i \left( \ket{0}^{\otimes d} + \ket{1}^{\otimes d} \right), \;\braket{\psi}{\psi}=1\}, $$ where $U_i \in SU(2)$