We analyze the electromagnetic amplitude for the leptonic decays of pseudoscalar mesons in the quark model, with particular emphasis on η → l + l − (l = e, µ). We evaluate the electromagnetic box diagram for a quark-antiquark pair with an arbitrary distribution of relative three-momentum p: the amplitude is obtained to all orders in p/m q , where m q is the quark mass. We compute B P ≡ Γ(η → l + l − )/Γ(η → γγ) using a harmonic oscillator wave function that is widely used in nonrelativistic (NR) quark model calculations, and with a relativistic momentum space wave function that we derive from the MIT bag model. We also compare with a quark model calculation in the limit of extreme NR binding due to Bergström. Numerical calculations of B P using these three parameterizations of the wave function agree to within a few percent over a wide kinematical range. Our results show 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 η → µ + µ − (unitary bound). However we find substantial deviations from the unitary bound in other kinematical regions, such as η, π 0 → e + e − . Using the experimental branching ratio for η → γγ as input, these quark models yield B(η → µ + µ − ) ≈ 4.3 × 10 −6 , within errors of the recent SAT-URNE measurement of 5.1 ± 0.8 × 10 −6 , and B(η → e + e − ) ≈ 6.3 × 10 −9 . While an application of constituent quark models to the pion should be viewed with particular caution, the branching ratio B(π 0 → e + e − ) ≈ 1.0×10 −7 is independent of the details of the above quark model wave functions to within a few percent.