The Madelung potential, as defined in Eq. 9, is a function of the nuclear configuration of the system. Therefore, thermal motions of particles in the system will lead to fluctuations in U i . At equilibrium, these fluctuations can be quantified in terms of a probability distribution function, P (U α ), where α is the species type for particle i. The width of P (U α ) can be specified by the variance, σ 2 α = ⟨(δU i ) 2 ⟩, where δU i = (U i -⟨U α ⟩). For simulation studies, this quantity can be evaluated by studying the statistics of U i for an ensemble of particles of type α sampled from equilibrium simulations. These fluctuations drive charge thermalization and are directly related to the solvent reorganization energy of Marcus theory. The long-range nature of Coulomb interactions presents a challenge when studying the fluctuations of U i . Specifically, due to the r -1 scaling, U i for individual configurations does not converge over a finite length scale. While the strength of contributions to U i decreases as r -1 , the volume element (and thus the number of contributing species) associated with a given distance increases as r 2 . The consequences of this scaling competition can be appreciated by considering a radially-resolved formulation of the Madelung potential for a spherically symmetric particle,