Paleointensities are key to understanding the formation and evolution of Earth and are determined from rocks which record magnetic fields upon cooling; however, experimental protocols for estimating paleointensities frequently fail. The primary reason is that laboratory protocols assume that rocks are dominated by uniformly magnetized, single-domain grains, instead of much more common non-uniformly magnetized grains. Our model for larger grains shows a multiplicity of stable domain states; with preferred states changing as a function of temperature. We show that domain state distribution depends on the thermal history of the sample—in nature and the laboratory. From numerical thermomagnetic modeling, we show that particles with non-uniform domain states will theoretically fail standard experimental paleointensity protocols, preventing us from determining reliable ancient geomagnetic field intensities. We propose that recognizing this type of behavior, and the resulting bias, will yield more reliable paleointensity records, and a better understanding of the Earth.