  Three photon decay rate of $J/\psi$ is studied using two $N_f=2$ twisted mass
gauge ensembles with lattice spacings $a\simeq 0.085$ fm (I) and $0.067$
fm(II). Using a new method, only the correlation functions directly related to
the physical decay width are computed with all polarizations of the initial and
final states summed over. Our results for such rare decay on the two ensembles
are: $\mathcal{B}_{I,II}(J/\psi\rightarrow 3\gamma)=(1.614 \pm 0.016 \pm
0.261)\times 10^{-5},(1.809 \pm 0.051 \pm 0.295)\times 10^{-5}$ where the first
errors are statistical and the second are estimates from systematics. We also
propose a method to analyze the Dalitz plot of the corresponding process based
on the lattice data which can provide direct information for the experiments.
