  We present 3D global kinetic pulsar magnetosphere models, where the charged
particle trajectories and the corresponding electromagnetic fields are treated
self-consistently. For our study, we have developed a cartesian 3D relativistic
particle-in-cell code that incorporates the radiation reaction forces. We
describe our code and discuss the related technical issues, treatments, and
assumptions. Injecting particles up to large distances in the magnetosphere, we
apply arbitrarily low to high particle injection rates and get an entire
spectrum of solutions from close to the Vacuum-Retarded-Dipole to close to the
Force-Free solution, respectively. For high particle injection rates (close to
FF solutions) significant accelerating electric field components are confined
only near the equatorial current sheet outside the light-cylinder. A judicious
interpretation of our models allows the calculation of the particle emission
and consequently the derivation of the corresponding realistic high-energy
sky-maps and spectra. Using model parameters that cover the entire range of
spin-down powers of Fermi young and millisecond pulsars, we compare the
corresponding model $\gamma$-ray light-curves, cutoff energies, and total
$\gamma$-ray luminosities with those observed by Fermi to discover a dependence
of the particle injection-rate, $\mathcal{F}$, on the spin-down power,
$\dot{\mathcal{E}}$, indicating an increase of $\mathcal{F}$ with
$\dot{\mathcal{E}}$. Our models guided by Fermi observations provide
field-structures and particle distributions that are not only consistent with
each other but also able to reproduce a broad range of the observed
$\gamma$-ray phenomenology of both young and millisecond pulsars.
