2026ApJ...999..204H

X-ray pulse-profile modeling of millisecond pulsars offers a direct route to measuring neutron star masses and radii, thereby constraining the dense-matter equation of state. However, standard analyses typically rely on ad hoc hotspot parameterizations rather than self-consistent physical models. While connecting surface heating directly to the magnetospheric geometry provides a more natural physical pathway, computing global magnetospheric solutions is too computationally expensive to perform on-the-fly during parameter inference. In this work, we bridge this gap by deriving fully analytic, first-principles expressions for surface return currents in mixed dipole-quadrupole magnetospheres. Working within force-free electrodynamics, we generalize the field-aligned current invariant Λ, the crucial scalar that maps the far-zone magnetic structure to the near-zone heating rate, from the standard dipole approximation to arbitrary quadrupolar configurations. We demonstrate that even when the quadrupole component is subdominant in the far zone (the mixing regime), using a dipole-based heating prescription fails to capture the significant enhancement or suppression of the return-current density on the polar cap. Our consistent quadrupole-aware framework reveals that these multipolar currents redistribute the surface heating, leading to systematic discrepancies in predicted pulse profiles that are amplified by atmosphere beaming and can reach ∼30% near pulse peaks. These results provide a rigorous analytic foundation for mapping global magnetic geometry to surface heating in multipolar magnetospheres, enabling physically consistent inference beyond the idealized dipole approximation.