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.

We performed ab-initio Particle-In-Cell (PIC) simulations of a pulsar magnetosphere with electron-positron plasma produced only in the regions close to the neutron star surface. We study how the magnetosphere transitions from the vacuum to a nearly force-free configuration. We compare the resulting force-free like configuration with ones obtained in a PIC simulation where particles are injected everywhere as well as with macroscopic force-free simulations. We found that although both PIC solutions have similar structure of electromagnetic fields and current density distributions, they have different particle density distribution. In fact in the injection from the surface solution, electrons and positrons counterstream only along parts of the return current regions and most of the particles leave the magnetosphere without returning to the star. We also found that pair production in the outer magnetosphere is not critical for filling the whole magnetosphere with plasma. We study how the current density distribution supporting the global electromagnetic configuration is formed by analyzing particle trajectories. We found that electrons precipitate to the return current layer inside the light cylinder and positrons precipitate to the current sheet outside the light cylinder by crossing magnetic field lines contributing to the charge density distribution required by the global electrodynamics. Moreover, there is a population of electrons trapped in the region close to the Y-point. On the other hand the most energetic positrons are accelerated close to the Y-point. These processes can have observational signatures that, with further modeling efforts, would help to distinguish this particular magnetosphere configuration from others.

Based on the Fermi observational data we reveal meaningful constraints for the dependence of the macroscopic conductivity $(\sigma)$ of dissipative pulsar magnetosphere models on the corresponding spin-down rate, $\dot{\mathcal{E}}$. Our models are refinements of the FIDO (Force-Free Inside, Dissipative Outside) models whose dissipative regions are restricted on the equatorial current-sheet outside the light-cylinder. Taking into account the observed cutoff-energies of all the Fermi-pulsars and assuming that a) the corresponding $\gamma-$ray pulsed emission is due to curvature radiation at the radiation-reaction-limit regime and b) this emission is produced at the equatorial current-sheet near the light-cylinder, we show that the \emph{Fermi}-data provide clear indications about the corresponding accelerating electric-field components. A direct comparison between the \emph{Fermi} cutoff-energies and the model ones reveals that $\sigma$ increases with $\dot{\mathcal{E}}$ for high $\dot{\mathcal{E}}$-values while it saturates for low ones. This comparison indicates also that the corresponding gap-width increases toward low $\dot{\mathcal{E}}$-values. Assuming the Goldreich-Julian flux for the emitting particles we calculate the total $\gamma-$ray luminosity $(L_{\gamma})$. A comparison between the dependence of the Fermi $L_{\gamma}$-values and the model ones on $\dot{\mathcal{E}}$ indicates an increase of the emitting particle multiplicity with $\dot{\mathcal{E}}$. Our modeling guided by the \emph{Fermi}-data alone, enhances our understanding of the physical mechanisms behind the high energy emission in pulsar magnetospheres.

We explore the emission properties of a dissipative pulsar magnetosphere model introduced by Kalapotharakos et al. (2014), comparing its high energy light curves and spectra, due to curvature radiation, with data collected by the Fermi LAT. The magnetosphere structure is assumed to be near the force-free solution. The accelerating electric field, inside the light-cylinder, is assumed to be negligible, while outside the light-cylinder it rescales with a finite conductivity ({\sigma}). In our approach we calculate the corresponding high energy emission by integrating the trajectories of test particles that originate from the stellar surface, taking into account both the accelerating electric field components and the radiation reaction forces. First we explore the parameter space assuming different value sets for the stellar magnetic field, stellar period, and conductivity. We show that the general properties of the model are in a good agreement with observed emission characteristics of young {\gamma}-ray pulsars, including features of the phase resolved spectra. Second we find model parameters that fit each pulsar belonging to a group of eight bright pulsars that have a published phase-resolved spectrum. The {\sigma} values that best describe each of the pulsars in this group show an increase with the spin-down rate $(\dot{E})$ and a decrease with the pulsar age, expected if pair cascades are providing the magnetospheric conductivity. Finally, we explore the limits of our analysis and suggest future directions for improving such models.