PULSARS ARE ROTATING MAGNETIZED NEUTRON STARS WHICH PRODUCE BROADBAND PULSED EMISSION. OBSERVATIONS WITH FERMI GAMMA-RAY SPACE TELESCOPE HAVE UNCOVERED MORE THAN 100 GAMMA-RAY EMITTING PULSARS DRAMATICALLY EXPANDING THE NUMBER OF SUCH SOURCES KNOWN AND PRODUCING UNPRECEDENTED DATA ON LIGHT CURVES AND SPECTRA IN THE GAMMA-RAY BAND. THE MAGNETOSPHERIC PROCESSES THAT LEAD TO THIS EMISSION ARE VERY INTERESTING AND CURRENTLY NOT UNDERSTOOD IN DETAIL. THE MAIN DIFFICULTY IS THE ABSENCE OF A SELFCONSISTENT MODEL OF THE PULSAR MAGNETOSPHERE THAT CAN PREDICT THE GLOBAL SHAPE OF THE MAGNETIC FIELD AND THE STATE OF THE PLASMA THROUGHOUT THE MAGNETOSPHERE. THE STANDARD PICTURE ENVISIONS THE MAGNETOSPHERE AS FILLED WITH PLASMA DUE TO SPACE-CHARGE-LIMITED FLOW OF CHARGES FROM THE SURFACE OF THE STAR. THESE CHARGES ARE ACCELERATED BY STRONG ELECTRIC FIELDS INDUCED BY RAPID STELLAR ROTATION. AS CHARGES ARE CURVED WHILE TRAVELING ALONG THE MAGNETIC FIELDS THEY EMIT CURVATURE RADIATION WHICH CAN PAIR-PRODUCE AND FILL THE MAGNETOSPHERE WITH PAIR PLASMA. THE CURRENTS AND CHARGES IN THIS PLASMA AFFECT MAGNETIC FIELD STRUCTURE AND PARTICLE ACCELERATION. UNTIL RECENTLY NO SELF-CONSISTENT SOLUTIONS OF THE PULSAR ELECTRODYNAMICS THAT ACCOUNTED FOR PLASMA SUPPLY AND PRODUCTION EXISTED. RECENTLY THE DEVELOPMENT OF MULTIDIMENSIONAL KINETIC SIMULATIONS HAS BROUGHT THE GOAL OF SELF-CONSISTENTLY CALCULATING PLASMA SUPPLY AND MAGNETOSPHERIC SHAPE OF PULSARS CLOSER TO REALITY. IN THIS PROPOSAL WE PLAN TO PERFORM THREE-DIMENSIONAL KINETIC SIMULATIONS OF PULSAR MAGNETOSPHERES USING PARTICLE-IN-CELL (PIC) METHOD. PIC IS AN ESTABLISHED TECHNIQUE FOR PLASMA SIMULATIONS THAT REPRESENTS THE PLASMA AS CHARGED MACROPARTICLES THAT ARE ADVANCED IN ELECTROMAGNETIC FIELDS. THE FIELDS ARE EVOLVED WITH PARTICLE CURRENTS AS SOURCES COMPLETING THE CYCLE. WE RECENTLY DEMONSTRATED THAT SELF-CONSISTENT PIC SIMULATIONS WITH PARTICLE EXTRACTION FROM THE STAR AND PAIR PRODUCTION IN THE MAGNETOSPHERE LEAD TO SOLUTIONS WHICH RESEMBLE THE FORCE-FREE STATE BUT CARRY INFORMATION ABOUT PARTICLE ACCELERATION AND BEAMING. WE WILL PERFORM A NUMERICAL SURVEY OF PULSAR MAGNETOSPHERES WITH A REALISTIC PAIR PRODUCTION PRESCRIPTION IN ORDER TO ASCERTAIN THE REGIMES OF PULSAR PARAMETER SPACE WHERE MAGNETOSPHERES WITH ABUNDANT PLASMA OCCUR. WE WILL ALSO ADDRESS THE FATE OF PULSARS NEAR THE DEATH LINE AND A SPECIAL CLASS OF WEAK PULSARS WITH LIMITED PAIR FORMATION. WE WILL SIMULATE THE PHYSICS OF RELATIVISTIC RECONNECTION WITH PAIR FORMATION IN ORDER TO PREDICT THE GAMMA-RAY LUMINOSITY AS A FUNCTION OF SPINDOWN ENERGY LOSS RATE FOR A RANGE OF PULSARS. WE WILL ALSO ANALYZE THE PLASMA DISTRIBUTION FUNCTION FOR EVIDENCE OF PLASMA INSTABILITIES THAT CAN LEAD TO COHERENT RADIO EMISSION. FINALLY WE WILL STUDY THE ENERGY AND PARTICLE DISTRIBUTION IN THE PULSAR WIND AND THE WAY THESE CAN TRANSFORM WITH WIND PROPAGATION. WE ALSO WILL SIMULATE THE INTERACTION OF PULSAR WINDS WITH OBSTACLES AND MAGNETOSPHERIC INTERACTIONS DURING NEUTRON STAR MERGERS. 50 YEARS AFTER DISCOVERY OF PULSARS THE STUDY OF PULSAR MAGNETOSPHERES IS NOW UNDERGOING AN EXCITING TRANSFORMATION WHERE A CALCULABLE STANDARD MODEL IS BEING FORMULATED. THE ADVENT OF GLOBAL KINETIC SIMULATIONS THAT WE WILL UTILIZE IN THIS RESEARCH WILL ENABLE US TO EXPLORE COMPLETE AB-INITIO MODELS OF PULSAR MAGNETOSPHERES AND THEIR EMISSION SOLVING ONE OF THE OLDEST UNSOLVED PROBLEMS IN ASTROPHYSICS.