First-Principles Polar-Cap Currents in Multipolar Pulsar Magnetospheres

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 \emph{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 $\Lambda$, 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 sub-dominant 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 $\sim 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.