Experimental evidences have indicated that cholesterol may adapt highly regular lateral distributions (i.e., superlattices) in a phospholipid bilayer. We investigated the formations of superlattices at cholesterol mole fraction of 0.154, 0.25, 0.40, and 0.5 using Monte Carlo simulation. We found that in general, conventional pairwise-additive interactions cannot produce superlattices. Instead, a multibody (nonpairwise) interaction is required. Cholesterol superlattice formation reveals that although the overall interaction between cholesterol and phospholipids is favorable, it contains two large opposing components: an interaction favoring cholesterol-phospholipid mixing and an unfavorable acyl chain multibody interaction that increases nonlinearly with the number of cholesterol contacts. The magnitudes of interactions are in the order of kT. The physical origins of these interactions can be explained by our umbrella model. They most likely come from the requirement for polar phospholipid headgroups to cover the nonpolar cholesterol to avoid the exposure of cholesterol to water and from the sharp decreasing of acyl chain conformation entropy due to cholesterol contact. This study together with our previous work demonstrate that the driving force of cholesterol-phospholipid mixing is a hydrophobic interaction, and multibody interactions dominate others over a wide range of cholesterol concentration.