Previous analyses of the Verwey transition in magnetite have been placed on a firm footing by a consistent analysis of the statistical properties of an array of quantum states associated with ferrous or ferric ions in octahedrally coordinated interstices of the spinel lattice. This collection is represented by an assembly of bonds and sites, Individual sites can be in one of three configurations: empty, as in the Fe3+ state, trapped, or polaronic if in the Fe2+ state. On neglect of high-energy states involving electron occupation of neighboring lattice sites, one arrives at an analytic equation of state for the order parameter in its dependence on temperature, which can be solved numerically. This equation is sufficiently flexible to handle both the first- and second-order transitions by appropriate changes in parameters, The present theory rationalizes the experimentally observed changes in the order of the Verwey transition that result from alterations in the sample stoichiometry, (C) 1999 Academic Press.