  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.
