We study the luminosity gap, Δm₁₂, between the first- and second-ranked galaxies in a sample of 59 massive ( 10¹⁵ M ) galaxy clusters, using data from the Hale Telescope, the Hubble Space Telescope, Chandra and Spitzer. We find that the Δm₁₂ distribution, p(Δm₁₂), is a declining function of Δm₁₂ to which we fitted a straight line: p(Δm₁₂) -(0.13 ± 0.02)Δm₁₂. The fraction of clusters with 'large' luminosity gaps is p(Δm₁₂≥ 1) = 0.37 ± 0.08, which represents a 3σ excess over that obtained from Monte Carlo simulations of a Schechter function that matches the mean cluster galaxy luminosity function. We also identify four clusters with 'extreme' luminosity gaps, Δm₁₂≥ 2, giving a fraction of More generally, large luminosity gap clusters are relatively homogeneous, with elliptical/discy brightest cluster galaxies (BCGs), cuspy gas density profiles (i.e. strong cool cores), high concentrations and low substructure fractions. In contrast, small luminosity gap clusters are heterogeneous, spanning the full range of boxy/elliptical/discy BCG morphologies, the full range of cool core strengths and dark matter concentrations, and have large substructure fractions. Taken together, these results imply that the amplitude of the luminosity gap is a function of both the formation epoch and the recent infall history of the cluster. 'BCG dominance' is therefore a phase that a cluster may evolve through and is not an evolutionary 'cul-de-sac'. We also compare our results with semi-analytic model predictions based on the Millennium Simulation. None of the models is able to reproduce all of the observational results on Δm₁₂, underlining the inability of the current generation of models to match the empirical properties of BCGs. We identify the strength of active galactic nucleus feedback and the efficiency with which cluster galaxies are replenished after they merge with the BCG in each model as possible causes of these discrepancies.