These files are scripts supporting all results reported in McLaskey et al., (2026). We found: Earthquakes occur in populations with few large events and many small ones, yet many earthquake rupture models employ a smooth stress field and other conditions that prohibit the co-occurrence of large earthquakes and smaller foreshocks and aftershocks. We describe populations of earthquakes that propagate and arrest within a stochastic stress field DeltaTau_pot(x) using a 1D fracture mechanics framework. DeltaTau_pot(x) is characterized by its mean (meanDeltaTau), standard deviation (stdevDeltaTau), and scaling exponent (mDeltaTau), which describes how DeltaTau_pot(x) changes as a function of wavelength. DeltaTau_pot(x) is related to the strain energy that fuels earthquakes, and it embodies the heterogeneous stresses that develop in fault zones with multi-scale geometrical complexity. Realistic populations of earthquakes, with more small ones than large ones, are produced by highly variable fields with stdevDeltaTau (8-80 MPa) that far exceeds meanDeltaTau, resulting in significant sections with highly negative DeltaTau_pot(x). Earthquake populations produced by such variable stress fields exhibit b-values that decrease with increasing meanDeltaTau and stress drops that are primarily correlated with stdevDeltaTau, are independent of magnitude, and may be limited by the strength of rocks at seismogenic depths. Our preferred models suggest 0 ≤ mDeltaTau ≤ 0.25 (mDeltaTau = 1.5 is self-similar, 1.0 is Brownian, 0 is white noise), consistent with expectations based on multiscale roughness measured on exhumed faults. This suggests that DeltaTau_pot(x) must be highly variable and highly negative, even at short wavelengths, a property that strongly influences the earthquake energy budget but is absent from state-of-the-art dynamic rupture models and unresolvable with kinematic finite-fault inversions.