We develop a new method that allows us to deal with two-photon intermediate states in a lattice QCD calculation. We apply this method to perform a first-principles calculation of the $\pi^0 \rightarrow e^+ e^-$ decay amplitude. Both the real and imaginary parts of amplitude are calculated. The imaginary part is compared with the prediction of optical theorem to demonstrate the effectiveness of this method. Our result for the real part of decay amplitude is $19.68(52)(1.10) \ \text{eV}$, where the first error is statistical and the second is systematic.