Quark model calculation of $\eta \to l^+ l^-$ to all orders in the bound   state relative momentum
  The electromagnetic box diagram for the leptonic decays of pseudoscalar mesons in the quark model is evaluated to all orders in ${\bf p} / m_q$, where ${\bf p}$ is the relative three-momentum of the quark-antiquark pair and $m_q$ is the quark mass. We compute $B_P \equiv \Gamma(\eta \to l^+ l^-) / \Gamma(\eta \to \gamma\gamma)$ using a popular nonrelativistic (NR) harmonic oscillator wave function, and with a relativistic momentum space wave function that we derive from the MIT bag model. We also compare with a calculation in the limit of extreme NR binding due to Bergstr\"om. Numerical calculations of $B_P$ using these three parameterizations of the wave function agree to within a few percent over a wide kinematical range. We find that the quark model leads in a natural way to a negligible value for the ratio of dispersive to absorptive parts of the electromagnetic amplitude for $\eta \to \mu^+ \mu^-$ (unitary bound). However we find substantial deviations from the unitary bound in other kinematical regions, such as $\eta,\pi^0 \to e^+ e^-$. These quark models yield $B(\eta \to \mu^+\mu^-) \approx 4.3 \times 10^{-6}$, within errors of the recent SATURNE measurement of $5.1 \pm 0.8 \times 10^{-6}$, $B(\eta \to e^+ e^-) \approx 6.3 \times 10^{-9}$, and $B(\pi^0 \to e^+ e^-) \approx 1.0 \times 10^{-7}$. 
B. Margolis, J. Ng, M. Phipps, and H. D. Trottier