The rare decay J/ψ→3γ, analog to Ortho-positronium decaying to 3γ in quantum electrodynamics, can provide a high precision test for the non-perturbative quantum chromodynamics. Such a decay process was first observed by CLEO collaboration in 2008 and then by BESIII in 2013. However, the relevant theoretical researches are very limited due to the dominant non-perturbative effects. We propose to use lattice QCD to study this problem. To this end, a new method has been proposed, that only the correlation functions directly related to the physical decay width are computed with all polarizations of the initial and final states summed over, to avoid the complicated decomposition for the matrix element. Using this new method, we present the first lattice result for this rare decay. Such a new method has also been applied for the decay ηc→2γ, and we obtain a lattice result that is consistent with the experimental one within two standard deviations for the first time. In the work of three-photon decay, we also put forward a scheme to analyze the Dalitz plot of the corresponding process based on the lattice data which can provide direct information for the relevant experiments.