  The particle-in-cell approach has proven effective at modeling neutron star
and black hole magnetospheres from first principles, but global simulations are
plagued with an unrealistically small separation between the scales where
microphysics operates and the system-size scales due to limited numerical
resources. A legitimate concern is whether the scale separation currently
achieved is large enough, such that results can be safely extrapolated to
realistic scales. In this work, our aim is to explore the effect of scaling
physical parameters up, and to check whether salient features uncovered by pure
kinetic models at smaller scales are still valid, with a special emphasis on
particle acceleration and high-energy radiation emitted beyond the light
cylinder. To reach this objective, we develop a new hybrid numerical scheme
coupling the ideal force-free and the particle-in-cell methods, to optimize the
numerical cost of global models. We propose a domain decomposition of the
magnetosphere based on the magnetic field topology using the flux function. The
force-free model is enforced along open field lines while the particle-in-cell
model is restricted to the reconnecting field line region. As a proof of
concept, this new hybrid model is applied to simulate a weak millisecond pulsar
magnetosphere with realistic scales using high-resolution axisymmetric
simulations. Magnetospheric features reported by previous kinetic models are
recovered, and strong synchrotron radiation above 100MeV consistent with the
Fermi-LAT gamma-ray pulsar population is successfully reproduced. This work
further consolidates the shining reconnecting current sheet scenario as the
origin of the gamma-ray emission in pulsars, as well as firmly establishes
pulsar magnetospheres as at least TeV particle accelerators.
