{'fraction_non_alphanumeric': 0.025487256371814093, 'fraction_numerical': 0.014092953523238382, 'mean_word_length': 4.700986193293885, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Electron positron pair flow and current composition in the pulsar magnetosphere

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.


INTRODUCTION
In the last decade, large steps have been taken in the understanding of the pulsar magnetosphere problem on both theoretical and observational sides. On the observational side the launch of the Fermi γ-ray Space Telescope (Atwood et al. 2009) led to the detection of 205 γ-ray pulsars (updated to 2016 February 22). On the theoretical side, the increased power of contemporary computers allowed verification and exploration of theoretical ideas with computationally expensive techniques. Global magnetospheric solutions have been obtained in ideal force-free electrodynamics (e.g. Contopoulos et al. 1999;Spitkovsky 2006), dissipative electrodynamics (e.g. Li et al. 2012;Kalapotharakos et al. 2012;Kalapotharakos et al. 2014;Kalapotharakos et al. 2017), and then with particle-in-cell simulations (PIC, Philippov & Spitkovsky 2014;Chen & Beloborodov 2014;Cerutti et al. 2015;Philippov et al. 2015a;Belyaev 2015a;Belyaev 2015b;Philippov et al. 2015b;Cerutti et al. 2016a;Cerutti et al. 2016b;Philippov & Spitkovsky 2017;Kalapotharakos et al. 2018). PIC simulations can simulate the pulsar magnetosphere from first principles, in contrast with dissipative electrodynamics simulations. The previous PIC works with pulsar magnetospheres addressed important problems such as the global magnetospheric currents and their composition, the dissipative processes and the electromagnetic emission, the role of pair production and of general relativity. In this paper, we focus on the magnetospheric structure that arises when particles are supplied only at the neutron star surface. We also inject a larger range of pair multiplicity than previous works which makes the simulation more realistic and allows us to explore the dependence of magnetosphere properties on particle injection rate. This paper is structured as follows: in § 2, we present the simulation setup; in § 3, we present how the magnetosphere transitions from vacuum to the force-free limit injecting particles only from the surface, in § 4, we compare the force-free solution obtained injecting particles from the surface to the one obtained injecting particles everywhere. We compare both macroscopic quantities and particle trajectories. In § 5 we discuss our findings; § 6 is the conclusion with a summary of the results and outlooks.

SIMULATION SETUP AND METHODS
The explicit electromagnetic PIC technique (Birdsall & Langdon 1991) works in a cycle, with the electromagnetic fields pushing some particles representing the plasma (usually called macro particles) determining the current that later enter as sources in the Maxwell equations to modify the fields. The simulations presented in this paper are obtained with the electromagnetic and relativistic PIC code C-3PA presented in Kalapotharakos et al. (2018), where simulations with particles injected everywhere in the domain are discussed. The only major difference is that we modified the code in order to save separately the contribution of the two particle species to the currents: where J pos and J ele represent the positron and electron currents, respectively. We will use these quantities mainly in section 4.1.
Here we want to analyze simulations of the pulsar magnetosphere with pair plasma supplied only at the surface of the star. The neutron star has effective radius of 0.36R LC (R LC = c/Ω is the light cylinder, where c is the speed of light and Ω the angular frequency of the neutron star) with conductive boundary conditions implemented below 0.28R LC and a kernel layer between 0.28R LC and 0.36R LC (Kalapotharakos et al. 2018). We inject particles according to the local magnetization in each cell, σ M = B 2 /8πnm e c 2 (B the magnetic field, n the number density, m e the electron mass), rather than at a fixed rate. Regulating the injection based on local magnetization prevents over-injection in the closed field regions where particles are trapped and the their number density increases more quickly than in the open field regions. Moreover, in our simulations we check that σ M > 10 everywhere except in the pulsar current sheet, where the magnetic energy is converted into particle kinetic energy. This is necessary in order to study a well magnetized plasma, such as in a real pulsar magnetosphere. Such an injection with respect to magnetization helps achieve this result. Particles are injected until the local magnetization σ M in each cell is below the threshold: Where r 0 is the stellar radius and r is the radial coordinate in spherical coordinates. Σ 0 is a quantity that is globally known in the simulation and that is assigned an initial value (which is, for example 2400). The particle supply is regulated to achieve a prescribed global injection rate, F. Every 10 time steps, we count all the particles injected at that time step and we compare this rate to F. If this rate is lower/higher than F, we decrease/increase Σ 0 . It takes less than a third of a stellar rotation to achieve a desired F. The unit of measure of F is: where ρ GJ is the Goldreich-Julian charge density in the pulsar polar cap (Goldreich & Julian 1969), where α is the inclination angle between the rotation axis and the magnetic moment, B 0 is the magnetic field at the star pole, q e is the electron charge, and is the area of the polar cap. The factor 2 accounts for the two poles. We introduce also the unit F 0 GJ = F GJ / cos α representing the ρ GJ for an aligned rotator. Pairs are injected with zero velocity. The simulation domains are cubes of side 9.6R LC with the neutron star rotating at the center. A perfectly matched layer (PML) is implemented at the outer boundary of the domain (Berenger 1994(Berenger , 1996Kalapotharakos & Contopoulos 2009). The main limitations of explicit PIC algorithms are the (temporal) resolution of the plasma frequency ω p and the (spatial) resolution of the skin depth λ sd . Not resolving ω p generates a numerical instability, while not resolving λ sd causes numerical heating in the plasma. We use a grid size of d = 0.02R LC and a time step small enough to resolve ω p everywhere in the domain with at least three time steps (dt = 0.003R LC /c). We do not resolve λ sd approximately in a sphere of radius 0.9R LC centered on the star. We observe that the numerical heating gives a small Lorentz factor to the particles. Even though this heating could push the particle energies up to γ ≈ 30, we find that only 1.7% of the particle population reach γ > 20. The simulations in this paper are obtained with a γ max ∼ 500 where Ω is the angular frequency of the neutron star rotation, r 0 the stellar radius, B 0 the magnetic field at the star radius. γ max is the Lorentz factor of an electron accelerated through all the voltage between the center of the polar cap of an aligned rotator and the last open field line (e.g. Rudermann & Sutherland 1975). The value of B 0 and then γ max is necessary to resolve the plasma frequency in the simulation, because a realistic value of B 0 would increase the characteristic charge density at the star surface (Equation 4), therefore, the time and space resolution needed would not be computationally accessible. The particles in our simulations are subjected to an enhanced radiation reaction to ensure a rapid cooling of the perpendicular momentum (Kalapotharakos et al. 2018). Low B 0 and small dt guarantee that the gyro frequency is resolved everywhere in the simulation. If not stated otherwise, when we show 2D slices of our simulations they are in the µ µ µ -Ω plane, where µ µ µ is the magnetic moment and Ω is the angular velocity that lies on the rotational axis. All the visualizations are obtained with VisIt (Childs et al. 2012).

FORMATION OF A FORCE-FREE LIKE MAGNETOSPHERE
The possibility to fill the entire magnetosphere and make it nearly force-free everywhere injecting particles only from the surface has been shown in Cerutti et al. (2016a). However, Cerutti et al. (2016a) focused mostly on the high-energy emission and not on the magnetosphere structure and its dependence on injection rate. As we increase the injection rate at the stellar surface we expect to find many different magnetosphere configurations ranging from charge-separated magnetospheres (Krause-Polstorff & Michel 1985a;Spitkovsky & Arons 2002) to a close to force-free solution (Contopoulos et al. 1999;Spitkovsky 2006).

Electromagnetic energy and Poynting flux
An important question we try to answer is how close the electromagnetic field structure of a magnetosphere is to the one of the force-free limit. The force-free configuration is characterized by the value of energy stored in the electromagnetic fields (see for example Bellan 2006). In each simulation we evaluated the average of the electromagnetic field energy density over the volume of a spherical shell starting outside the boundary layer of the star and extending up to 2.5 R LC . We compared this result with the same quantity evaluated in a macroscopic force-free electrodynamics simulation (Kalapotharakos et al. 2012). We found that our global PIC models presenting both the Poynting flux close to the theoretical value and the electromagnetic field energy density close to the force-free electrodynamics simulations provide the current and charge density distribution and field structure close to the ideal force-free ones. We tested this method to the simulations with injection everywhere and then applied it to the simulations with injection from the surface. The simulations that inject particles over the entire computational domain showed that the average electromagnetic energy density value remains constant over a broad range of injection rates after reaching the force-free value. In Figures 1 and 2 we show the average electromagnetic energy density versus the injection rate for pulsars with α = 45 • . We plot also the average magnetic energy density, however the two quantities show a very similar behavior. In this paper we focus on magnetospheres obtained injecting particles close to the star surface and we tested the validity of the The empirical expression ofĖ for a force-free magnetosphere is (Spitkovsky 2006) In both these cases,Ė is totally carried away by the Poynting flux. In Figure 4 we show the Poynting flux evaluated through spherical shells for pulsars with α = 45 • with increasing injection rate. The Poynting flux decreases along the radial direction because it gets dissipated by J · E in the volume enclosed by the shell. However, we note that the dissipation, especially for the higher injection rates, takes place close to the R LC , and beyond that near the equatorial current sheet. In Figure 4 we see that the maximum dissipation happens for F = 3.5F GJ and does not exceed 15% of the Poynting flux at the surface. The maximum percentage in dissipation happens for F = 0.5F GJ (∼ 20%) but this is because the absolute value of the Poynting flux is lower. The dissipation decreases monotonically going toward the force free solution from F = 3.5F GJ . The dissipation for the nearest to force-free solution is ∼ 6%. These results are consistent with the γ-ray efficiency in Kalapotharakos et al. (2018) for magnetospheres with particles supplied everywhere.

Screening of the accelerating electric fields
Studying how the magnetosphere changes its configuration while the number of particles injected increases is interesting because it shows possible configurations in which a pulsar can operate. However, exact solutions reflecting what happens in a real pulsar can be obtained only having a detailed knowledge of the pair production processes in the magnetosphere, works like Timokhin & Arons (2013) go in this direction. In our study we can observe how the magnetosphere behaves if the particles are injected only at the surface, but without capturing the specific physical mechanism behind the injection. The magnetic field structure, shown in Figure 5, begins to resemble the force-free structure already from F ∼ 5F GJ , while for F = 0.5F GJ it is very close to the vacuum one. The evolution of the current configuration seems to follow the magnetic field structure, in particular once the magnetic field resembles the force-free one, a clear separatrix/Y-point/current sheet configuration is present. As accelerating electric field we consider E 0 (Gruzinov 2008;Li et al. 2012) that is defined as: igure 3. In this figure we show the close to force-free simulations obtained injecting particles from the surface. We show cases that cover the whole range of inclination angles α. The color is the divergence of the electric field, while the field lines are the magnetic field lines projected onto the poloidal plane. The 0 • case is obtained with F = 5FGJ = 5F 0 GJ , the 45 • case with F = 12.5FGJ = 8. , Poynting flux and dissipation for α = 45 • Poy:F = 0.5F GJ Dis:F = 0.5F GJ Poy:F = 3.5F GJ Dis:F = 3.5F GJ Poy:F = 6.5F GJ Dis:F = 6.5F GJ Poy:F = 9.5F GJ Dis:F = 9.5F GJ Poy:F = 12.5F GJ Dis:F = 12.5F GJ force-free Ė with E 0 ≥ 0. The electric field gets gradually screened with higher particle injection, but we can see that the regions that are hardest to screen are the polar cap outflow region and the separatrix region. When the current sheet region is formed, it never gets completely screened.

COMPARISON BETWEEN THE SOLUTIONS APPROACHING THE FORCE-FREE LIMIT
In this section, we present the solution approaching the force-free limit (F = 12.5F GJ ) obtained with particles supplied only at the surface for α = 45 • , as described in section 3.1.

The macroscopic quantities
In Figure 6 we compare the solution obtained with force-free electrodynamics to two PIC simulations approaching the force-free limit: one injecting particles everywhere in the simulation domain and one injecting particles only close to the star surface. The ∇ · E (which represents the charge density), the total current and its sign are very similar in all the solutions. The situation changes when we look in detail at how these structures are sustained. The two solutions are obtained with a different F. This is due to the different way in which the injected particles are distributed in the two injection schemes. For injection everywhere, particles are placed at larger radii, filling a greater volume than the injection from the surface configuration. The more the particles are injected at larger distance, the less they contribute to the force-freeness of the inner magnetosphere because many quickly leave the domain. However, the injection scheme is not crucial to obtain a force-free configuration, as we can see from the very similar values of charge and current densities. We define the multiplicity M as the number of particles present per charge at a given location.
When the particles are injected at the surface, the multiplicity reaches higher values close to the star, but lower values at larger radii. This can be seen in Figure 7. As we described in section 2, our code accounts for the different current contributions of electrons and positrons. We plot their absolute value components in Figure 8.
In the case when particles are injected everywhere, electron currents are present mostly in the negatively charged regions, while positron currents are present mostly in the positively charged regions (for the charge density plot see Figure 6). When we inject particles from the surface, we notice that the electrons (positrons) have an important current contribution even in positively (negatively) charged regions (Figure 8,9). As we saw in Figure 6, the total current densities are very similar, but the difference in electron and positron current densities indicates that when we inject particles only from the surface there are regions where electrons and positrons are streaming in the same direction with electron and positron currents almost cancelling each other. This is possible if these particles are injected in a zone where the accelerating electric field is screened enough to not reverse their initial velocities acquired by the numerical heating. The same qualitative behavior is present in the electron-positron pair cascades at the polar cap. In fact, in polar cap cascades most of the pairs are pro-duced above the pair formation front wth some initial Lorentz factors (Harding & Muslimov 2001;Timokhin & Harding 2015). When particles are injected everywhere in the domain, they are supplied wherever they are needed. Instead when particles are injected only from the surface, they need to arrange themselves in a different way to satisfy the current and charge density requirements of the magnetosphere. In Figure 9 where we plot J · B/B for the electron and positron components, this scenario becomes clear. For the simulation with injection from the surface, we can see that J·B/B indicates counter streaming flows (where the two components have the same color in the same region) only on the negative branch of the separatrix and in a thin layer just above the neutron star surface, where particles are injected. In the same simulation, there is also a clear component of positrons flowing out from the polar cap together with the electrons that is very weak in the simulation with the particles injected everywhere. In the electron component the positive branch of the separatrix that connects the star surface to the Ypoint changes sign: in the surface injection case, on this branch the electrons are flowing outward from the star, while in the injection everywhere case they are flowing in. Therefore, we see that the zones with availability of pairs (where the pair creation happens) greatly influence the underlying currents of the single species and this has important consequences that we will outline in the Discussion (Section 5).

The particle trajectories
We studied the trajectories of the particles in our simulation approaching the force-free limit with injection from the surface (F = 12.5F GJ ). First we describe trajectories followed by the majority of the particles. In general the most energetic particles (the ones that reach γ from ∼ 50 up to ∼ 180 in the 45 • close to force-free simulation) are mainly positrons accelerated along the field lines that constitute the separatrix\Y-point\current sheet complex. The particles gain most of their energy in the proximity of the Y-point. At intermediate energies (γ ∼ 40) we find the electrons flowing from the polar cap. Then at low energies (γ < 30) we find the bulk of the flow with electrons and positrons generally flowing out together. As expected, positrons are dominant in positively charged regions and electrons are dominant in negatively charged regions. In Figure  10, we see some examples of these trajectories. All the trajectories are shown in the corotating frame. Studying PIC simulations of the pulsar magnetosphere is interesting because it can provide solutions to problems that are present in the force-free electrodynamics force-free electrodynamics PIC injection everywhere PIC surface injection Figure 6. The divergence of the electric field, the absolute value of the current density and the projection of the current density on the magnetic field lines for magnetospheres close to the force-free limit obtained with three different kinds of techniques: force-free electrodynamics, PIC with particles injected everywhere in the domain, PIC with particles supplied only at the surface. As it can be seen, these plots are very similar.
limit. One of these problems, is how the current structure of the force-free configuration could be sustained by particles in real pulsars, outside the strict force-free limit. It is reasonable to assume that the field structure of a pulsar magnetosphere is stationary, therefore, the amount of charge in the magnetosphere should remain constant. Because of the charge conservation it follows that the current leaving the star should be balanced by a current entering the star. When Contopoulos et al. (1999) found the first force-free solution for a dipolar magnetic field, the currents were going from the star to infinity and coming from infinity to the star, one through the polar cap flow, the other mainly through the current sheet and separatrix, and a smaller part on a few open magnetic field lines close to the last open magnetic field lines. The surface charge density of the current sheet has some puzzling features. The charge of an aligned force-free magnetosphere at the Y-point should be negative inside the Y-point and positive outside of it (Lyubarskii 1990;Timokhin 2006). However, the current is continuous through the Y-point, but its composition should change to obtain a charge of a different sign. 1976) was thought to provide electrons flowing backward and positrons flowing outward where the charge density changes sign. Studying the particle trajectories, we have found that the pulsar magnetosphere does not need pair production in an outer gap to fill the magnetosphere and maintain the charge density distribution of the separatrix/Y-point/current sheet complex. As we said above, there is an outward flow of electrons from the polar cap. Some of them have γ ∼ 1 and flow very close to the region where the current changes sign. There a low electric field drags part of the low energy electron distribution into the returning current and separatrix region, where the majority of them form the returning current by the electric field that reverse their velocity ( Figure 11). The particles with a higher energy are not affected by this because this accelerating electric field is too small. This phenomenon does not happen at a specific height, but it appears continuous up to the Ypoint. The crossing of field lines happens in a region where λ sd is not resolved, therefore this noise could in principle affect the trajectories. We checked a randomly selected sample of all the electron trajectories starting from the region where the electrons that turn back originate. We found that ∼ 70% of the electrons coming from this region are turned back inside the light cylinder. This behavior is different from the action expected of kicks due to random noise; thus we believe that this phenomenon results from non-fluctuating, low electric fields.
Another phenomenon concerns the outgoing electrons. Some electrons remain stuck at the Y-point and they circle all around the light cylinder, see Figure 12 (a similar behavior was shown also in Cerutti et al. (2016a)). This happens because the solution tends toward the force-free one where the sign of the charge density is negative before the Y-point and positive after it. The resulting electric field accelerates positrons and deflects the electrons. In this motion the electrons get energized, 40 γ 90. Once they are in this regime electrons have two possibilities: either falling back toward the star (and they mix with the electrons of Figure 11) or flying out following other field lines (they do not usually fly far out into the current sheet). When they fall back they lose their energy by radiation reaction (the accelerating fields in that region are not strong enough to sustain the Lorentz factor they had reached). Beyond the 1.5R LC there are very few electrons that turn back and the number of these returning ones decrease drastically with distance, in contrast with what was presented in Cerutti et al. (2015). However, the Cerutti et al. (2015) simulation setup is different from ours; for example, they do not include radiation reaction.
To complete the picture, we must understand the origin of the positrons that support the charge density change of sign through the Y-point. Most of the positrons in the current sheet come along the separatrix; but extra positrons are needed in the current sheet to account for the current of the returning electrons inside the Y-point. They come from the polar cap flow (they are flying out with the electrons), close to the returning current region and the separatrix and they cross field lines outside the light cylinder to enter into the positively charged region and then the current sheet ( Figure 13). We checked that positrons indeed cross magnetic field lines (for the returning electrons it was obvious because of the shapes of their trajectories) looking at the cosine of the angle electrons positrons injection everywhere surface injection Figure 8. Electronic and positronic modulus of the current densities for PIC simulations close to the force-free limit: one is with particles supplied only at the star's surface, the other is with particles injected everywhere. Notice that the scale values are different: in the case of injection from the surface the two components are much more intense.
between the particle momentum and the local electric field outside of the current sheet p · E pE (12) with p particle momentum and E electric field. We use this criterion because in the force-free limit, where particles flow exactly along the magnetic field lines and E⊥B, p has components along B and E × B. E × B keeps the particle on the rotating field line. In this limit, the Expression 12 is always 0. In Figure 13, we see that the expression 12 becomes singificantly > 0. This happens in regions where E · B BE (where E is mostly perpendicular to B). In fact, we identified in our PIC simulations the regions where E has significant components parallel to B (the yellow opaque volume in Figure  13) using E · B BE > 0. 15 (13) and the field line crossing happens outside of this region. The region defined by the Expression 13 traces quite well the regions of reconnecting B. The last three kinds of trajectories that we just described involve the crossing of magnetic field lines. The theoretical gyro radius of these low energy particles is very small while the use of strong radiation reaction forces makes it even smaller. Therefore, the corresponding crossing of the magnetic field lines is not driven by a large gyro radius due to the use of a low magnetic field, but by unscreened electric fields. So far, we showed trajectories for α = 45 • case. However, the most well studied case is the aligned rotator and we looked for the same trajectories in this case as well and we show them in Figure 14. In this case, the electrons circling around the Y-point form a cloud of negative charge that appears as an increase of the nega-electrons positrons injection everywhere surface injection Figure 9. Electronic and positronic projection of the current density on the magnetic field lines for PIC simulations close to the force-free limit: one is with particles supplied only at the star's surface, the other is with particles injected everywhere.
tive charge density where the separatrix touches the Ypoint. This behavior was first noted in force-free electrodynamics simulations by Timokhin (2006). Therefore, we can safely say that this current structure is sustained by particles crossing field lines mainly inside 1.0R LC and for sure inside 2.0R LC . The main mechanisms are all driven by low electric fields operating on the low energy part of the particle distribution. For the nearly orthogonal rotator case, the structure of the current sheet is very different (Kalapotharakos et al. 2012), therefore we decided to reserve this study for the future. Summarizing, we can say that the pulsar magnetosphere structure approaching the force-free solution with particles injected from the surface has these features: 1) Electrons and positrons stream outward together in the polar cap outflow.
2) The electrons that flow back to the star cross field lines, either from the polar cap outflow into the return-ing current region inside the light cylinder or after circling around the Y-point.
3) Positrons flow out on the separatrix and get accelerated close to the Y-point into the current sheet. Some positrons enter in the current sheet beyond the Y-point crossing field lines. In Figure 15 we show a sketch for an aligned rotator to help the reader understand the particle trajectories outlined above. The action of the non ideal electric field is indicated.

DISCUSSION
In this section we discuss the results we obtained in the context of pulsar magnetosphere theory. First of all, we should note again the limitations of PIC codes in capturing the physical quantities of the pulsar magnetosphere. The open field voltage from Equation 6 for a pulsar with B 0 ∼ 10 12 G and P ∼ 0.1s would be γ max ∼ 10 9 , while we and all the previous studies with PIC simulations (e.g Figure 10. Most common particle trajectories in the simulation approaching force-free with particles supplied close to the surface. The color on the trajectories represents the Lorentz factor (γ). A) is a "not so highly" accelerated positron. B) is a positron flowing out from the polar cap flow at low energy. C) is an electron flowing out from the polar cap flow at low energy. D) is a "highly accelerated" positron. E) is an intermediate energy electron flowing out from the polar cap. In the picture on the bottom, we have a volume rendering of E0 (Equation 9, 10) that identifies the current sheet: we can see that the difference in acceleration between A) and D) is due to the strength of E0 on the trajectory. This non uniformity in E0 is found only through PIC simulations and it can be useful to model the γ-ray emission. All the trajectories are in the corotating frame. Figure 11. Electron trajectories falling back on the star from different heights. The Lorentz factor γ is the color on the trajectory. We added a red arrow representing the magnetic moment, because the zoomed region could be difficult to identify. Note that the color scale has a different range with respect to all the others shown in this paper for trajectories. This is because of the low energy of these particles. The trajectories are in the corotating frame. 1000. This is necessary because we need to resolve ω p everywhere in our system, as we explained in Section 2. Therefore, the particle energy distribution is squeezed into a narrow range and it cannot be simply linearly stretched or shifted to higher energies 1 . Some of the previously referenced works on PIC pulsar magnetospheres how their results depend on the number of particles injected. We showed that many properties and accelerating gaps are connected to the number of particles injected into the magnetosphere, and for high injection rate the magnetosphere reached the force-free limit for the whole range of inclination angles (0 • , 45 • , 85 • ). In our case, we are confident in claiming that the current composition we discussed is characteristic of a magnetosphere approaching the force-free limit with particles supplied over the whole stellar surface. In Section 3.2, we show that a maximum in the dissipation occurs at an intermediate injection rate between the charge separated solutions and the force-free case (Figure 4). Using dissipative models, Kalapotharakos et al. (2012) also found a maximum in dissipation at an intermediate conductivity. We can qualitatively associate the increasing conductivity of these models with the increasing injection rate in our simulations. Gruzinov (2013); Contopoulos (2016b) have recently proposed weak pulsars, that are magnetosphere configurations that present a larger dissipation than the force-free magnetosphere. These solutions are expected from a particle supply only at the neutron star surface, therefore they should be comparable to our simulations. We identified the solution with the highest dissipation at an intermediate F.
Considering also the previous results, we suggest that a weak pulsar magnetosphere originates for these intermediate F. A similar behavior is reported in Cerutti et al. (2015) for an aligned rotator and in Kalapotharakos et al. (2018) for simulations that inject particles over the entire computational domain. From the study of the macroscopic quantities and confirmed by the study of particle trajectories, we found that if particles are injected at the surface there are only a few regions of counter streaming particles. This is important because it was not clear if the currents in the magnetosphere were built of counter streaming species or not. For example, the photon-photon pair production in the current sheet (Lyubarskii 1996) would be inhibited if the electrons and positrons flow out in the same direction, as is the case in our simulation of plasma injection from the surface. Photon-photon pair production in the current sheet is implemented with simple prescriptions in other works (e.g Chen & Beloborodov 2014;Philippov & Spitkovsky 2017). We think that local simulations of this phenomenon in the pulsar current sheet and at the Y-point are needed to address this issue more Figure 13. Positron trajectories flowing from the polar cap into the current sheet. The star with the magnetic field lines is plotted on the background. The yellow opaque surface is the pulsar current sheet. We used an opaque profile instead of a volume rendering profile (as in Figure 10) in order to facilitate seeing that the trajectories shown are initially outside the current sheet. We selected the current sheet according to Equation 13. On the left panel the color on the trajectories is the normalized projection of p on E, as defined in Equation 12. We can see that the value is larger than 0.5 in many parts of the trajectories. On the right panel the color is the divergence of the electric field along the trajectories. We can see that these positrons are transitioning from a negatively charged region, to a positively charged region. The trajectories are in the corotating frame.
carefully. However, there are also other mechanisms that could trigger the pair production, like considering other sources of photons, as was discussed for the outer gap (Chiang & Romani 1994). Another phenomenon impacted by this effect is the hypothesis of generation of the radio emission through the two stream instability (e.g. Usov 2002 for a review). The two stream instability could still occur in the returning current region on the polar cap rim (note that for 45 • we found only one of the two branches to have counter streaming particles, see 4.1) and below the pair formation front (e.g. Harding & Muslimov 1998) that is not resolved by this simulation. When we look at the energetics of the most common particle trajectories (section 4.2) we see that the highest energy particles gain most of their energy close to the Y-point and they are outflowing positrons. Outflowing, energetic particles in the current sheet can produce light curves and spectra (Brambilla et al. 2015;Kalapotharakos et al. 2014Kalapotharakos et al. , 2017) that match well with those of the Fermi pulsars (Abdo et al. 2013). In these works, the particle acceleration was allowed for outflowing particles only after the R LC where the Y-point is located and the particles were injected only at the surface. In Figure 5, we see that for magnetospheres that are far from the force-free limit, there can be some acceleration below the Y-point along the separatrix and above the polar cap. This probably indicates that the young γray pulsars we selected for Brambilla et al. (2015) have a magnetosphere close to the force-free limit. However, some millisecond γ-ray pulsars could have emission coming from these lower altitude gaps (Johnson et al. 2014), and their spectra would not suffer any magnetic pair attenuation because of the lower magnetic field. Other energetic particles are the electrons that are circling on the Y-point but they have lower energy than these outgoing positrons. These electrons, and the other particles that we showed crossing field lines, naturally gain pitch angles, thus breaking the ideal force-free limit. This makes them natural candidates for the non-thermal synchrotron emission observed at MeV energies and in the hard x-rays (e.g. Kuiper & Hermsen 2015) and it would explain the misalignment with the GeV emission that is observed in certain cases (e.g. Marelli et al. 2014). Obviously the electrons that are circling around the Ypoint are more promising candidates, but it is difficult to give final answers when nine orders of magnitude are squeezed into three. It would be extremely interesting to see at which energies this crossing of field lines happens for real pulsars with γ max ∼ 10 9 , that would result in high energy particles with γ ∼ 10 7 because of the radiation reaction. However, these kinds of works are Figure 14. The particle trajectories that we presented above, but for the aligned rotator. For clarity, we removed the azimuthal component of their trajectory. The magnitude of the azimuthal components were similar to the trajectories shown in the Figures  11, 12, 13, for the 45 • case. The red dot indicates where the particle is injected. All these trajectories had large azimuthal components. The color is a label to help distinguish their intricate trajectories. A) is one of the electrons that starts from the polar cap outflow and get turned back into the separatrix and the returning current. B) is one of the electrons that reaches the Y-point, circles for a while and then flies out following another field line. C) is one of the electrons that reaches the Y-point, circles with huge azimuthal components, and then falls back to the star losing energy. D) is one of the positrons that starts in the polar cap flow, close to the separatrix, crosses field lines and then enters the current sheet.
helpful because they indicate a direction to follow and new hypotheses to be tested that did not emerge previously. Future missions looking at pulsars in the MeV band (e.g. De Angelis et al. 2017;Moiseev et al. 2015;Hunter et al. 2014) could help unveil the mystery and constrain the models. Looking at the particle trajectories showed also how a current structure close to the one of the force-free configuration can be sustained injecting particles only from the surface of the neutron star. Probably magnetospheres with pair production at the Y-point and/or in the current sheet or in other locations would settle on a different configuration, closer to the simulations where particles are injected everywhere. These scenarios will produce different signatures in the heating of the polar cap in addition to the heating generated by the pair production below the pair formation front (e.g. Harding & Muslimov 2001). These signatures could be potentially observed and discriminated with NICER (Gendreau et al. 2012;Özel et al. 2016).

CONCLUSIONS
We presented PIC simulations of the pulsar magnetosphere injecting particles only at the surface of the neutron star. We outlined the regime of our magnetosphere simulations, and we showed some properties of the solutions that are in between charge-separated solutions, and the force-free limit. Then we concentrated on a solution approaching the electromagnetic field structure of the force-free solution. We showed the different macroscopic quantities and compared them to a solution approaching the force-free limit where particles are injected everywhere in the domain. Following the differences, we studied the particle trajectories. The main findings are the behavior of the main flow of electrons and positrons in the magnetosphere and the complex mechanisms that sustain the current configuration. We discussed both their theoretical and observational implications, underlining how comparing the same quantities in solutions obtained with different particle supply could be crucial for a theory-observation comparison. Future work could try to reproduce the injection in pulsar magnetospheres that are far from the force-free limit. In order to do this more realistically a more self-consistent injection mechanisms should be implemented, as for example injection dependent on the microphysics of the polar cap pair cascades. As we already mentioned in the text, more detailed studies of the pulsar current sheet (DeVore et al. 2015) and Y-point would be helpful as well. Figure 15. In the sketch above we illustrate for the aligned rotators case the electron and positron trajectories summarized at the end of Section 4.2. The action of the low electric fields acting on low energy particles is highlighted by the green quasi transparent shapes. The stronger electric fields are highlighted by the yellow quasi transparent shapes. A) an electron trajectory that flows out from the polar cap. B) one of the positron trajectories that flow together with the electrons out of the polar cap. C) trajectory of a positron accelerated at high energy by the electric field at the Y-point and in the reconnecting current sheet. D) a trajectory of an electron that arrives close to the Y point and is bounced back by the same electric field that accelerates the positrons at very high energy. This bouncing makes the electron circle all around the RLC and gain energy. Then the electron can either fall back to the star D1), or fly away D2) depending on where this circling around the RLC takes it. E) a trajectory of a low energy electron that crosses magnetic field lines towards the return current region. Eventually, this electron is deflected by a low electric field and returns to the star. Electrons approximately on the same trajectory but with a higher energy are not deflected. F) a trajectory of a positron that flows out of the polar cap at low energy as B), but outside the RLC it crosses the field lines and enters in the current sheet far from the Y-point. The figure is a modified version of a figure from Timokhin (2006).