We present 294 pulsars found in GeV data from the Large Area Telescope LAT on the Fermi Gammaray Space Telescope. Another 33 millisecond pulsars MSPs discovered in deep radio searches of LAT sources will likely reveal pulsations once phaseconnected rotation ephemerides are achieved. A further dozen optical andor Xray binary systems colocated with LAT sources also likely harbor gammaray MSPs. This catalog thus reports roughly 340 gammaray pulsars and candidates 10 of all known pulsars compared to 11 known before Fermi. Half of the gammaray pulsars are young. Of these the half that are undetected in radio have a broader Galactic latitude distribution than the young radioloud pulsars. The others are MSPs with six undetected in radio. Overall 236 are bright enough above 50 MeV to fit the pulse profile the energy spectrum or both. For the common twopeaked profiles the gammaray peak closest to the magnetic pole crossing generally has a softer spectrum. The spectral energy distributions tend to narrow as the spindown power dotE decreases to its observed minimum near 10SUP33SUP erg sSUP1SUP approaching the shape for synchrotron radiation from monoenergetic electrons. We calculate gammaray luminosities when distances are available. Our allsky gammaray sensitivity map is useful for population syntheses. The electronic catalog version provides gammaray pulsar ephemerides properties and fit results to guide and be compared with modeling results.
['Binary pulsars', 'Millisecond pulsars', 'Pulsar timing method', 'Rotation powered pulsars', 'Pulsars', 'Gamma-rays', 'Celestial objects catalogs', 'Neutron stars', 'Radio pulsars', '153', '1062', '1305', '1408', '1306', '637', '212', '1108', '1353', 'Astrophysics - High Energy Astrophysical Phenomena']
The Third Fermi Large Area Telescope Catalog of GammaRay Pulsars
{'The Third Fermi Large Area Telescope Catalog of Gamma-ray Pulsars': "D. A. Smith, 1, 2 S. Abdollahi, 3 M. Ajello, 4 M. Bailes, 5 L. Baldini, 6 J. Ballet, 7 M. G. Baring, 8 C. Bassa, 9 J. Becerra Gonzalez, 10 R. Bellazzini, 11 A. Berretta, 12 B. Bhattacharyya, 13 E. Bissaldi, 14, 15 R. Bonino, 16, 17 E. Bottacini, 18, 19 J. Bregeon, 20 P. Bruel, 21 M. Burgay, 22 T. H. Burnett, 23 R. A. Cameron, 19 F. Camilo, 24 R. Caputo, 25 P. A. Caraveo, 26 E. Cavazzuti, 27 G. Chiaro, 26 S. Ciprini, 28, 29 C. J. Clark, 30, 31 I. Cognard, 32, 33 A. Corongiu, 22 P. Cristarella Orestano, 12, 34 M. Crnogorcevic, 35, 25 A. Cuoco, 16, 17 S. Cutini, 34 F. D'Ammando, 36 A. de Angelis, 37 M. E. DeCesar, 38 S. De Gaetano, 15, 14 R. de Menezes, 16, 39 J. Deneva, 38 F. de Palma, 40, 41 N. Di Lalla, 19 F. Dirirsa, 42 L. Di Venere, 15 A. Dom'ınguez, 43 D. Dumora, 1 S. J. Fegan, 21 E. C. Ferrara, 35, 44, 25 A. Fiori, 6 H. Fleischhack, 45, 25, 44 C. Flynn, 5, 46 A. Franckowiak, 47 P. C. C. Freire, 48 Y. Fukazawa, 49 P. Fusco, 14, 15 G. Galanti, 26 V. Gammaldi, 50, 51 F. Gargano, 15 D. Gasparrini, 28, 29 F. Giacchino, 28, 29 N. Giglietto, 14, 15 F. Giordano, 14, 15 M. Giroletti, 36 D. Green, 52 I. A. Grenier, 53 L. Guillemot, 32, 33 S. Guiriec, 54, 25 M. Gustafsson, 55 A. K. Harding, 56 E. Hays, 25 J.W. Hewitt, 57 D. Horan, 21 X. Hou, 58, 59 F. Jankowski, 32 R. P. Johnson, 60 T. J. Johnson, 38 S. Johnston, 61 J. Kataoka, 62 M. J. Keith, 63 M. Kerr, 64 M. Kramer, 48, 63, 65 M. Kuss, 11 L. Latronico, 16 S.-H. Lee, 66 D. Li, 67, 68 J. Li, 69, 70 B. Limyansky, 60 F. Longo, 71, 72 F. Loparco, 14, 15 L. Lorusso, 14, 15 M. N. Lovellette, 73 M. Lower, 61 P. Lubrano, 34 A. G. Lyne, 63 Y. Maan, 13 S. Maldera, 16 R. N. Manchester, 61 A. Manfreda, 6 M. Marelli, 26 G. Mart'ı-Devesa, 74 M. N. Mazziotta, 15 J. E. McEnery, 25, 35 I.Mereu, 34, 12 P. F. Michelson, 19 M. Mickaliger, 63 W. Mitthumsiri, 75 T. Mizuno, 76 A. A. Moiseev, 44, 35 M. E. Monzani, 19, 77 A. Morselli, 28 M. Negro, 78, 25 R. Nemmen, 39 L. Nieder, 30, 31 E. Nuss, 79 N. Omodei, 19 M. Orienti, 36 E. Orlando, 80, 19 J. F. Ormes, 81 M. Palatiello, 72, 71, 82, 37 D. Paneque, 52 G. Panzarini, 14, 15 A. Parthasarathy, 48 M. Persic, 72, 83 M. Pesce-Rollins, 11 R. Pillera, 14, 15 H. Poon, 49 T. A. Porter, 19 A. Possenti, 22 G. Principe, 71, 72, 36 S. Rain'o, 14, 15 R. Rando, 18, 84, 85 S. M. Ransom, 86 P. S. Ray, 64 M. Razzano, 6 S. Razzaque, 87, 54 A. Reimer, 74 O. Reimer, 74 N. Renault-Tinacci, 88 R. W. Romani, 19 M. S'anchez-Conde, 51, 50 P. M. Saz Parkinson, 60 L. Scotton, 79 D. Serini, 15 C. Sgr'o, 11 R. Shannon, 46 V. Sharma, 44 Z. Shen, 89 E. J. Siskind, 90 G. Spandre, 11 P. Spinelli, 14, 15 B. W. Stappers, 63 T. E. Stephens, 91, 25", 'Fermi-LAT Collaboration': '0  \nFigure B-1. Same as Figure 9, for PSR J0007+7307.  \n<!-- image -->  \nFigure B-2. Same as Figure 9, for PSR J0633+1746 (the Geminga pulsar).  \n<!-- image -->  \n0  \n0  \nFigure B-3. Same as Figure 9, for PSR J1311 -3430. The ⋆ symbol indicates that we lacked information needed for absolute phase alignment of the radio and gamma-ray profiles (see Section 5.1).  \n<!-- image -->  \nFigure B-4. Same as Figure 9, for PSR J2229+6114. The ⋆ symbol indicates that we lacked information needed for absolute phase alignment of the radio and gamma-ray profiles (see Section 5.1).  \n<!-- image -->  \n0  \nFigure B-5. Same as Figure 12, for PSR J0007+7303. Here, the P3 (bridge) component was too weak to fit in the different energy intervals, hence the points are not shown (see Section 5.3).  \n<!-- image -->  \nFigure B-6. Same as Figure 12, for PSR J0633+1746 (the Geminga pulsar).  \n<!-- image -->  \nFigure B-7. Same as Figure 12, for PSR J1311 -3430.  \n<!-- image -->  \nFigure B-8. Same as Figure 12, for PSR J2229+6114.  \n<!-- image -->  \nFigure B-9. Typical gamma-ray spectral energy distributions (SEDs). The plots were made using the python script sed plotter.py and data files in the Online Material. The SEDs for all the gamma-ray pulsars detected as point sources (TS > 25) are in the Online Material as well. Similar to those provided with 4FGL-DR3, they use the refits done for some pulsars, described in Section 6. When b fr exists, both b 23 and b fr fits are shown, with the b fr uncertainty envelope.  \n<!-- image -->', 'ABSTRACT': 'We present 294 pulsars found in GeV data from the Large Area Telescope (LAT) on the Fermi Gamma-ray Space Telescope . Another 33 millisecond pulsars (MSPs) discovered in deep radio searches of LAT sources will likely reveal pulsations once phase-connected rotation ephemerides are achieved. A further dozen optical and/or X-ray binary systems co-located with LAT sources also likely harbor gamma-ray MSPs. This catalog thus reports roughly 340 gamma-ray pulsars and candidates, 10% of all known pulsars, compared to ≤ 11 known before Fermi . Half of the gamma-ray pulsars are young. Of these, the half that are undetected in radio have a broader Galactic latitude distribution than the young radio-loud pulsars. The others are MSPs, with 6 undetected in radio. Overall, ≥ 236 are bright enough above 50 MeV to fit the pulse profile, the energy spectrum, or both. For the common two-peaked profiles, the gamma-ray peak closest to the magnetic pole crossing generally has a softer spectrum. The spectral energy distributions tend to narrow as the spindown power ˙ E decreases to its observed minimum near 10 33 erg s -1 , approaching the shape for synchrotron radiation from monoenergetic electrons. We calculate gamma-ray luminosities when distances are available. Our all-sky gamma-ray sensitivity map is useful for population syntheses. The electronic catalog version provides gamma-ray pulsar ephemerides, properties and fit results to guide and be compared with modeling results.  \nKeywords: catalogs - gamma rays: observations - pulsars: general - stars: neutron', '1. INTRODUCTION': "Fewer than a dozen gamma-ray pulsars were known when the Fermi Gamma-ray Space Telescope was launched on 2008 June 11, and the extent and diversity of the population and its role in Galactic dynamics were subject to debate (Thompson 2008). Fermi 's primary instrument, the Large Area Telescope (LAT, Atwood et al. 2009), quickly established that the gamma-ray population is large and varied and is the dominant GeV gamma-ray source class in the Milky Way (Abdo et al. 2010a, The First Fermi LAT Catalog of Gamma-ray Pulsars, hereafter 1PC). The 46 pulsars in 1PC (6 months of data) grew to 132 in 2PC (3 years of data), the Second Fermi -LAT gamma-ray pulsar catalog (Abdo et al. 2013). This third gamma-ray pulsar catalog (based on 12 years of data) characterizes 294 confirmed gamma-ray pulsars, and tabulates 33 MSPs for which gamma-ray pulsations have not yet been seen but  \nlikely will be once accurate rotation ephemerides are established. We further tabulate LAT sources likely to reveal new 'spider' MSPs, and LAT sources co-located with known pulsars, some of which may ultimately reveal gamma-ray pulsations. Roughly 340 gamma-ray pulsars and candidates are thus now known, or about 10% of the > 3400 currently known pulsars (see Table 1).  \nThese results build on much previous work. GeV pulsations from the Crab were glimpsed using a balloon-borne instrument at the start of the 1970s (McBreen et al. 1973), followed by the pulsed detection of Vela by the SAS-2 satellite (Thompson et al. 1975). The COS-B satellite improved the measurements (Swanenburg et al. 1981). In the 1990s EGRET on the Compton Gamma-Ray Observatory (CGRO) saw six GeV pulsars (Thompson et al. 1999) while a rare pulsar that is brighter below 100 MeV than above, PSR B1509 -58 was detected with COMPTEL on CGRO (Kuiper et al. 1999; Abdo et al. 2010b). EGRET data revealed three other strong candidates: PSRs J0659+1414 (Ramanamurthy et al. 1996) and J1048 -5832 (Kaspi et al. 2000), and the millisecond pulsar PSR J0218+4232 (Kuiper et al. 2000). AGILE discovered gamma-ray pulsations from PSR J2021+3651 before Fermi 's launch (Halpern et al. 2008). All 11 of these pulsars were quickly confirmed using LAT data, and are noted in Figure 1. Geminga, seen with EGRET, was undetected at radio wavelengths (Bignami & Caraveo 1996; Abdo et al. 2010c) and has turned out to be the prototype of about half of the young gamma-ray pulsars. The 294 pulsars reported here are more numerous than the 271 sources, all object classes combined, in the 3 rd EGRET source catalog (Hartman et al. 1999).  \nFigure 1 shows that the discovery rate since launch is steady. Table 1 breaks the numbers down by category. Long-term radio observations by the 'Pulsar Timing Consortium' (Smith et al. 2008) enabled about half the discoveries and, importantly, also allowed for an unbiased sample of pulsars not seen in gamma rays (Smith et al. 2019). The 'Pulsar Search Consortium' (Ray et al. 2012), later joined by FAST (Li et al. 2018a; Wang et al. 2021) and the TRAPUM 1 project on MeerKAT (Clark et al. 2023a), discovered large numbers of radio millisecond pulsars (MSPs) at the positions of unidentified gamma-ray sources, leading to the subsequent detection of gamma-ray pulsations. Gamma-ray blind searches of unidentified sources revealed radio-quiet pulsars that make up a quarter of the current sample (see e.g. Clark et al. 2017; Wu et al. 2018).  \nThe discovery rate is sustained by innovations in how we detect pulsations. Everimproving blind search algorithms (Pletsch & Clark 2014; Nieder et al. 2020a) allowed the first discovery of a radio-quiet MSP (Clark et al. 2018). The increasingly sophisticated use of an optical companion's orbit reduces the parameter space searched using Einstein@home 2 to discover gamma-ray MSPs in binary systems with perturbed orbits (Nieder et al. 2022). Another example is a method that allows photon weighting  \n(see Section 2) even if the unpulsed gamma-ray source is undetected (Bruel 2019; Smith et al. 2019). Underlying all analysis efforts are the improved sensitivity and low energy reach afforded by the Pass8 reconstruction method (Atwood et al. 2013; Bruel et al. 2018).  \nAs a result, not just the numbers but the variety of gamma-ray pulsars continues to grow. The minimum spindown power ˙ E is now 20 × lower than the pre-launch expectation of 10 34 erg s -1 (Smith et al. 2008) 3 . The fastest known field MSP, PSR J0952 -0607, was found at the location of a gamma-ray source and subsequently timed with LAT data (Bassa et al. 2017b; Nieder et al. 2019). A third globular cluster, NGC 6652, was found to have gamma-ray emission dominated by a single MSP, PSR J1835 -3259B (Gautam et al. 2022b; Zhang et al. 2022), with a fourth recently reported, PSR J1717+4308A in M92 (Zhang et al. 2023). The gamma-ray flux and pulse profile of PSR J2021+4026 in the γ Cygni supernova remnant (SNR) changed during mode transitions in 2011 and 2018 (Allafort et al. 2013; Razzano et al. 2023). The LAT sees more than forty 'spider' MSPs and a dozen candidates (see Section 7), compact binary systems where the pulsar wind ablates its companion. Spiders fall into two categories: 'black widows' have companion masses 0 . 01 M ⊙ < M c < 0 . 05 M ⊙ and orbital periods P B < 10 h whereas 'redbacks' have M c > 0 . 2 M ⊙ and P B < 1 d. Gamma-ray timing of 35 stable MSPs for over 12 years usefully constrains the intensity of gravitational waves from super massive black hole binaries in the hearts of distant galaxies. The upper limit should become a measurement in the coming years (Ajello et al. 2022). Following the methods initially applied to PSR J1555 -2908 (Nieder et al. 2022), we may be poised to detect planets in multi-year orbits in several compact binary MSP systems.  \nThis heterogeneous population can be classified by comparing the spin period ( P ) and the period derivative ( ˙ P ), shown in Figure 2. Throughout this paper, we call pulsars in the main population 'young' to distinguish them from the much older MSPs, thought to be spun up to rapid periods via accretion from a companion (Alpar et al. 1982), although e.g. the accretion-induced collapse of white dwarfs might also create MSPs (Gautam et al. 2022a).  \nAll known gamma-ray pulsars are rotation-powered (RPPs): LAT has not yet detected accretion-powered pulsars nor the magnetars that populate the upper-right portion of the P ˙ P plane, for which the dominant energy source is magnetic field decay (Parent et al. 2011). An interesting exception is a LAT detection of a few photons for a few minutes from an extragalactic magnetar giant flare (Ajello et al. 2021a). The locations of all 294 gamma-ray pulsars on the sky are shown in Figure 3. The P ˙ P diagram shows diagonal lines of constant ˙ E , τ c , and B S derived from the timing information as follows. For an orthogonal rotator, the magnetic field on the neutron star surface at the magnetic equator (the rotation pole)  \nis B S = (1 . 5 I 0 c 3 P ˙ P ) 1 / 2 / (2 πR 3 NS ) ≃ 3 . 2 × 10 19 G √ P ˙ P . The 'characteristic age' τ c = P/ (2 ˙ P ) assumes that magnetic dipole braking is the only energy-loss mechanism, that the magnetic moment and inclination do not change, and that the initial spin period was much less than the current period. τ c thus approximates true age well for some young pulsars, and poorly for MSPs. We set the neutron star radius to R NS = 10 km, and c is the speed of light in a vacuum.  \nThe 4 th Fermi -LAT source catalog (Abdollahi et al. 2020), and specifically Data Release 3 (DR3, Abdollahi et al. 2022, hereafter 4FGL) 4 characterizes 6658 point and extended sources using 12 years of LAT data. Half of the sources are various blazar classes of active galactic nuclei, but a third remain unassociated with objects known at other wavelengths. Radio and gamma-ray pulsation searches at the positions of unidentified sources have yielded fully half of the gamma-ray pulsars. The discoveries continue, as detailed in Sections 3 and 7. The 4FGL spectral, flux, and variability measurements are used throughout this work for the pulsar searches and characterization.  \nWe provide the pulsar catalog as in FITS 5 and spreadsheet file formats, along with other information, available in the journal electronic supplement as well as on the Fermi Science Support Center (FSSC) servers 6 . Appendix D describes the online material.  \n- 4 https://fermi.gsfc.nasa.gov/ssc/data/access/lat/12yr catalog  \nFigure 1. Cumulative number of known gamma-ray pulsars, beginning with the launch of Fermi . The crosses show the numbers included in the first (1PC) and second (2PC) catalogs of LAT pulsars and their publication dates. Some key discoveries are highlighted. See also Table 1.  \n<!-- image -->", '2. OBSERVATIONS': "Atwood et al. (2009) describe the Fermi LAT, and Abdo et al. (2009a), Ackermann et al. (2012a), and Ajello et al. (2021b) report on-orbit performance. Fermi carries another instrument, the Gamma-ray Burst Monitor (GBM; Meegan et al. 2009), which was not used to prepare this catalog.  \nThe LAT is a pair-production telescope composed of a 4 × 4 grid of towers. Each tower consists of a stack of tungsten foil converters interleaved with silicon-strip particle tracking detectors, mated with a hodoscopic cesium-iodide calorimeter. A segmented plastic scintillator anti-coincidence detector covers the grid to help discriminate charged particle backgrounds from gamma-ray photons. The LAT field of view is ∼ 2.4 sr. For most of the Fermi mission, the primary operational mode has been a sky survey where the satellite rocks between a pointing above the orbital plane and one below the plane after each orbit. In this mode, the entire sky is imaged every two orbits ( ∼ 3 hours) and any given point on the sky is observed ∼ 1 / 6 th of the time. For one year beginning 2013 December, the normal survey mode was changed to favor exposure to the Galactic center region 7 . Survey mode was again modified after a  \nFigure 2. Pulsar spindown rate, ˙ P , versus the rotation period P . Green dots indicate young, radio-loud (RL) gamma-ray pulsars and blue squares show 'radio-quiet' (RQ) pulsars, defined as S 1400 < 30 µ Jy, where S 1400 is the radio flux density at 1400 MHz. Red triangles are millisecond gamma-ray pulsars. Black dots indicate pulsars phase-folded in gamma rays without significant pulsations. Phase-folding was not done for pulsars shown by gray dots. Orange triangles are radio MSPs discovered at the positions of previously unassociated LAT sources, hidden by red triangles when gamma pulsations were subsequently found. The rest are listed in Table 6, and plotted at ˙ P = 5 × 10 -22 when ˙ P is unavailable. The solid black diagonal is the radio deathline of Eq. 4 of Zhang et al. (2000). Shklovskii corrections to ˙ P have been applied only to gamma-ray MSPs with measured proper motion (see Section 4.3).  \n<!-- image -->  \nsolar panel rotation drive stopped moving on 2018 March 16, detailed in Section 2.1 of Ajello et al. (2021b).  \nThe LAT is sensitive to gamma rays with energies E from 20 MeV to over 300 GeV, with an on-axis effective area of ∼ 8000 cm 2 above 1 GeV. Multiple Coulomb scattering of the electron-positron pairs created by converted gamma rays degrades the perphoton angular resolution, with the average 68% containment radius varying as E -0 . 8 from 5 · at 100 MeV to 0 . · 1 at 10 GeV 8 .  \nThis energy-dependent point spread function (PSF), the backgrounds from the complex diffuse emission and nearby sources, and the source spectrum are encapsulated in the photon weights ( w i , Kerr 2011), which give the probability that a photon", '3 rd Fermi LAT Pulsar Catalog': 'Table 7 (continued)  \nNote -The best known distances of the 150 young pulsars detected by the Fermi LAT. For the DM method, the reference is for the Dispersion Measurement and the distance was obtained with the YMW16 model for the Galactic electron density (Yao et al. 2017). DMM is the distance at which the YMW16 DM reaches 95% of its maximum value for that line of sight. Method P is parallax, K is kinetic (HI absorption), KPLK combines the two and corrects for Lutz-Kelker bias (Verbiest et al. 2012), X indicates an X-ray measurement of the hydrogen column density, and O englobes optical and other methods.  \n<!-- image -->  \nT able 8 (c ontinue d)  \n<!-- image -->  \nT able 8 (c ontinue d)  \n<!-- image -->', '3. DISCOVERY AND TIMING': "Inclusion in the main catalog requires a statistically significant pulsation detection in the Fermi LAT gamma-ray data. The following subsections describe the paths to detection, as well as a brief description of gamma-ray timing that often strengthens the initial signal.  \nIn all cases, detecting and characterizing pulsations requires a rotation ephemeris (or a 'timing model') to convert photon arrival times t i to neutron star rotational phases ϕ i . This 'folding', which 'stacks' photons at the same rotational phase, allows a pulse to rise above the background since ≪ 1 photons are collected per pulse even from the brightest pulsars (Kerr 2022). We detect pulsations with the H-test (de Jager et al. 1989; de Jager & Busching 2010), a statistical test for discarding the null hypothesis that a set of photon phases is uniformly distributed. For N γ gamma rays, the m = 20 harmonic weighted version of the H-test statistic (Kerr 2011) is  \nH ≡ max( Z 2 mw -4 × ( m -1) , 1 ≤ m ≤ 20) , (1)  \nFigure 3. Pulsar sky map in Galactic coordinates (Hammer projection). Symbols as in Figure 2.  \n<!-- image -->  \nwith  \nZ 2 w ≡ 2 ∑ i w 2 i m =20 ∑ k =1 ( α 2 kw + β 2 kw ) , (2)  \nand α kw and β kw the empirical trigonometric coefficients α kw ≡ ∑ N γ i =1 w i sin(2 πkϕ i ) and β kw ≡ ∑ N γ i =1 w i cos(2 πkϕ i ). The w subscripts indicate that this is a photonweighted version of the test, and the w i are the photon weights, evaluated via spectral analysis (see also 2PC, Abdo et al. 2013). For m = 20, which we adopt universally, the cumulative distribution function for H in the asymptotic limit is P ( H ≥ x ) = exp( -0 . 398405 x ) (Kerr 2011), and 3 σ , 4 σ , and 5 σ thresholds correspond to H = 14, H = 24, and H = 36. The H-test is unbinned, well-suited to the extremely sparse gamma-ray pulsar data: the LAT often detects only one photon in tens of thousands (millions, for MSPs) of pulsar rotations. Bruel (2019) gives corrections for small photon counts. Pulsars with narrow, sharp peaks are easier to detect than pulsars with broad peaks (see Figure 5 of Hou et al. 2014). All pulsars in our sample are detected with m = 8. Using m = 20 incurs little computational cost and insures sensitivity to putative exotic profile shapes. Kerr (2011) also showed that large m does not cause false positive detections.  \nFigure 4 highlights some aspects of gamma-ray phase folding. Figure 9 and Appendix B show other example profiles. The top-most frame shows a weighted phase histogram, duplicated over a second rotation. Section 5 describes the profile fit overlaid in blue. The phase-aligned 1 . 4 GHz radio pulse overlaid in red comes from the radio timing observations used to create the rotation ephemeris, in this specific case by Parthasarathy et al. (2019). The horizontal dashed line shows the gamma-ray  \nbackground level, estimated from the photon weights as 1 N bin ( ∑ i w i -∑ i w 2 i ), the sum of the expected contribution to the weights of background photons not associated with the pulsar. A phase histogram baseline exceeding the background level may indicate the presence of unpulsed magnetospheric emission. Section 5.1 gives details.  \nThe next frame below it shows the phase drifting after the last radio time of arrival used to model the neutron star rotation, indicated by the green horizontal dashed line. The start of this timing model's 'validity' range, before Fermi 's launch, is not shown. Pulsars with irregular spindown, as for this young high-˙ E pulsar, require extra parameters to model the rotation, and accuracy of the extrapolation past the validity range rapidly degrades. Stable pulsars can be modeled with few parameters, often accurately predicting the neutron star rotation for years before and/or after validity. The right-hand frame shows the weighted H-test increasing as data accumulated over the years, a nearly straight line for most pulsars. Changes in slope can result from phase drifts, as in this case, or from increased background due to e.g. a nearby flaring blazar (see Section 6.6), or from changes in the LAT's exposure to the pulsar's sky position. Exposure per unit time increased for this pulsar during 2014 (mid-year was MJD 56810) when LAT pointed more frequently towards the Galactic center, visible in the Time versus Phase plot, without however affecting the H-test growth. Smith et al. (2019) show simulations for which Poisson fluctuations in the very low rate of photon arrivals for pulsars near detection threshold significantly perturb the H-test time evolution. Slope variations due to pulsar flux changes nearly never happen: rare exceptions are the young pulsar PSR J2021+4026 (Allafort et al. 2013; Razzano et al. 2023) and transitional MSPs like PSR J1023+0038 (see Stappers et al. 2014, and Appendix A) and J1227 -4853 (Roy et al. 2015; Johnson et al. 2015), which lack detectable gamma-ray pulsation during their accretion states.  \nIf different analyses were tried, the number of trials must be accounted for in the pulsation significance calculation, requiring a higher H-test value to claim a detection. Before the advent of photon weighting, we sometimes varied the minimum energy and/or the angular extent of our dataset, to explore the pulsar's spectral hardness and local background level. Using a LAT source's measured spectrum to calculate weights allows one single trial for this parameter space.  \nWhen the source spectrum is too faint to measure, we use 'simple weights' with a parameter µ w of unknown optimal value, and thus, possible additional trials. The µ w parameter is the logarithm of the energy at which the log-Gaussian weighting function peaks. Since Smith et al. (2019) we have been using six trials for simple weight gamma-ray pulsation searches using radio or X-ray ephemerides, and an Htest > 25 ( p = 4 . 7 × 10 -5 , or > 4 . 1 σ ) detection threshold. We fold the entire data set 3 times, using µ w = (3 . 2 , 3 . 6 , 4 . 0), and 3 more times restricting the data to the ephemeris validity period. Although these 6 trials are not independent, a conservative estimate for the chance probability for a false positive detection with a H-test > 25 threshold is 6 × 4 . 7 × 10 -5 , or 0.3 for a sample of 1000 pulsars. In 2PC we required  \nime (MJD)  \nT  \nFigure 4. Top: gamma-ray phase histogram for PSR J1648-4611 discovered by Kramer et al. (2003), overlaid with the 1 . 4 GHz profile (red) obtained during Parkes radio telescope timing (Parthasarathy et al. 2019). The blue curve shows a fit to the histogram, and the horizontal dashed line is an estimate of the background level (see Section 5.1). Bottom-left: Phaseogram over the course of the mission. Dots indicate photons, with the grayscale set according to the photon weight. Right: H-test significance accumulated over the course of the mission. The green horizontal dashed line shows when the ephemeris validity finishes (see Section 3.1).  \n<!-- image -->  \nH-test > 36 ( > 5 σ ). Smith et al. (2019) discovered 16 new gamma-ray pulsars, in part because of this refined, relaxed threshold. Several more pulsars in this catalog were later found in the same way.  \nThis prescription may miss some pulsars if the ephemeris extrapolates well enough that the accumulation of additional data beyond the validity range yields H> 25, but poorly enough that the full data yield H< 25. Such cases also arise from positive fluctuations, so we avoid further trials factors by limiting our search to the two combinations just described.  \nAs further explained in Section 3.4, when the ephemeris is established independently of the LAT, using radio or X-ray observations (Section 3.1), the pulsar can be up to 20 times fainter in gamma rays than if an unidentified LAT point source (Section 3.2) guides a 'blind' gamma-ray search. Once pulsations are found during some epoch of LAT data, gamma-ray timing can extend and improve the ephemeris (Section 3.5).  \nUsing the above detection criteria, we report 294 gamma-ray pulsars in the breakdown tabulated in Table 1 with their distribution on the sky shown in Figure 3. The 150 young pulsars are named in Table 2, and the 144 MSPs are in Table 3. The new pulsars in 3PC were mostly more 'difficult' to discover than those reported in the earlier catalogs. For pulsars found with just a few foldings using a radio ephemeris, discussed in Section 3.1, difficulty stems from their extreme faintness ( < 1 photon per month): an analysis using many trials (for example, from exploring different data intervals) yields significance too low to distinguish from statistical fluctuations. Furthermore, we continue to discover radio MSPs at the positions of unidentified FGL sources, described in Section 3.3. These are difficult both in radio, because of eclipses, scintillation, and intrinsic faintness, and in gamma rays if their signal-to-noise ratio is low and the ephemeris validity is brief, that is, if the radio observations cover only a few years of the LAT mission. The blind searches, discussed in Section 3.4, are now finding MSPs, including those with short, varying orbital periods, feared as perhaps impossible before launch (Ransom 2007). These searches require many orders of magnitude more computing power than for young pulsars. The catalog includes 44 'spider' systems, and in Section 7.1 we tabulate several candidates co-located with gamma-ray sources, for which pulsations will likely be seen in the coming years. Timing models for 'noisy' pulsars require many parameters which presents difficulty in maintaining a coherent timing solution, e.g. with frequent radio monitoring. Our methodical pursuit of ever more difficult gamma-ray pulsars in the LAT's unbiased all-sky data set means that our notion of gamma-ray pulsars is greatly enriched compared to 1PC.  \nTable 2 . Some parameters of young LAT-detected pulsars  \nTable 2 continued on next page  \nFermi-LAT Collaboration  \nTable 2 (continued)  \nTable 2 continued on next page", 'rd LAT Pulsar Catalog': '3 Fermi Table 10 (continued)  \nNote -Same as Table 9, for the gamma-ray MSPs.', '3.1. Using Known Rotation Ephemerides': "The most straightforward way to search for gamma-ray pulsations is to fold the LAT data using an existing rotation ephemeris that is valid for the full 12 or more years on orbit. To that end, radio and X-ray astronomers shared > 1400 rotation ephemerides with the LAT team, ∼ 40% of the over 3400 known rotation-powered pulsars (mostly from the ATNF Pulsar Catalog v1.69 9 , Manchester et al. 2005, as well as some unpublished pulsars, see Table 1). Table 4 includes radio telescopes that contributed timing models used in the LAT pulsar searches. These timing models yielded the discovery of ∼ 110 gamma-ray pulsars. A special effort was made to search for LAT pulsations from the energetic ( ˙ E > 10 34 erg s -1 ) subset of the pulsar population, representing about 10% of all pulsars. We have folded > 90% of the high-˙ E pulsars (see Figure 5), revealing gamma-ray pulsations in about two thirds of them. Importantly, the ephemerides provided by the radio community sample the entire P -˙ P plane. The > 1200 pulsars that were phase-folded without seeing gamma-ray pulsations are shown as black dots in Figure 2. Folding so many pulsars gives a largely unbiased view of known pulsars which may emit gamma rays not visible from earth, discussed by Johnston et al. (2020).  \nThe availability of so many ephemerides is possible because of the long-term timing campaigns of several radio telescopes, summarized in Smith et al. (2019). Astronomers originally organized support for the LAT as a 'Pulsar Timing Consortium' (Smith et al. 2008) and have continued as the Fermi mission proceeds. We aim for P/ 50 ephemeris accuracy (0 . 02 in phase), which we have generally achieved: possible gamma-ray pulse width smearing at this level is unlikely to impede pulsation discovery, given that the narrowest known pulse is P/ 33 (0 . 03 in phase) for PSR J1959+2048 (Guillemot et al. 2012a). Gamma-ray pulse widths and background rates are such that phase histograms generally require far fewer than 100 bins (Figure 4 uses 50 bins). Of the ∼ 110 gamma-ray pulsars found by phase-folding using radio ephemerides, 20 are too faint for the 4FGL-DR3 catalog: they rise above the background only when the photons accumulate in a restricted phase range, and thus probe a population beyond that accessible by the source catalog methods. We do not include additional spectral analysis of these sources, so they are absent from the analysis of Section 6.  \nFully describing the rotation of young pulsars, which have rotational irregularities known as timing noise, over a timespan of many years often requires many noisesmoothing parameters, which may be covariant with parameters that have a physical interpretation. That is to say, 'whitened' solutions used for our analysis and distributed with this paper accurately calculate rotational phase, but should be used with care when studying braking indices, proper motions, or other pulsar properties.  \nFor a few pulsars, X-ray timing was either combined with radio timing to construct an ephemeris (e.g. PSR J2022+3842, Limyansky 2022), or was used exclusively (e.g. PSR J0540 -6919 in the Large Magellanic Cloud, Marshall et al. 2016). Section 7.1 further addresses searches for GeV pulsations in X-ray bright pulsars.  \nFinally, the timing models enable phase alignment between observations at different wavelengths from different observatories. The absolute phase reference is given by the Tempo2 parameters for the arrival time (TZRMJD) and location (TZRSITE) of a reference pulse of a particular frequency (TZRFREQ). The dispersion measure (DM) determines the frequency-dependent delays in the interstellar medium and is necessary to align radio and high-energy data. Carefully-aligned pulse profiles provide information about the relative geometry of the different emission regions.", '3.2. Pulsar Search Targets': "About half of the gamma-ray pulsars were discovered in searches around LAT sources with pulsar-like properties. In addition, many LAT sources are co-located with previously known radio and/or X-ray pulsars (see Table 5). Section 7 discusses the prospects for finding gamma-ray pulsations from the latter, as well as from 'spider'-like binary systems found in LAT unidentified sources, likely to yield still more MSP discoveries. This section describes how we select and improve candidate  \ntargets for the deep radio searches and the gamma-ray blind searches described in Sections 3.3 and 3.4.", '3.2.1. Pulsar-like LAT Sources': "The LAT source catalog lists 'associations' and 'identifications' for two-thirds of the sources. 'Association' means that the probability that an object known at some other wavelength is responsible for the LAT emission is estimated to be > 80%, where the prior for the Bayesian probability is the sky distribution for each population class (Abdollahi et al. 2020). 'Identification' requires either detection of gammaray pulsations (in the case of pulsars) or an additional match with observations at another wavelength. Examples are simultaneous blazar flares or supernova remnant morphology. Pulsars are the largest class of identified sources, assigned 'PSR' or 'MSP' classes in the 4FGL catalog (capital letters indicate identification). A spatial match between a known pulsar and a catalog source without gamma-ray pulsations yields association classes 'psr' and 'msp' (lower case indicates association without confirmed identification). The 4FGL catalog table has a second association column for low probability ( < 20%) associations. Table 5 is complementary to the association process, listing simple co-locations without regard for non-pulsar populations.  \nAdditional observables are used to rank candidates in the long list of unidentified sources. Pulsar spectral energy distributions (SEDs), νF ν = E 2 dN dE , have sharp cutoffs in the GeV range (Section 6.1) that distinguish them from other categories of gamma-ray sources. The 4FGL catalog includes a variability index for each source, which can help distinguish between pulsars, whose fluxes are stable (see Section 6.6), and the 10 × more common blazars, which can exhibit gamma-ray flares. Early successes in selecting pulsar candidates from among the unidentified LAT catalog sources exploited SED curvature versus flux variability correlations (Ackermann et al. 2012b). Subsequent works used machine learning methods to refine the selections (Mirabal et al. 2012; Lee et al. 2012; Saz Parkinson et al. 2016; Wu et al. 2018; Luo et al. 2020; Finke et al. 2021), or visual inspection and ranking of the LAT source spectra (Camilo et al. 2015). Regardless of the ranking scheme, these lists of pulsar-like unassociated LAT sources have provided a large number of candidate pulsar positions that have been targeted by radio, X-ray, and gamma-ray searches. As one recent example, PSR J1653 -0158 (Nieder et al. 2020b) was recently found in the brightest remaining pulsar candidate of Saz Parkinson et al. (2016).  \nIn 4FGL-DR3, about 80% of the unidentified sources have error ellipses with semimajor axes between 0 . 05 and 0 . 15 degrees. Improved localization can be valuable for the high-frequency radio observations needed to pierce the thick electron column density (i.e. high DM values) in, for example, the direction of the Galactic center. Improved localization also reduces the computing cost of blind period searches, and it can allow a match with optical or X-ray sources necessary to motivate time-consuming observations and analyses (see Section 3.2.2). Radio telescope arrays can coherently combine data from multiple antennas to form tied array beams which are typically  \nmuch smaller than the primary beam. These also benefit from improved localization, since computational resources generally limit the number of synthesized beams that can be formed and searched.  \nFigure 5. Left-hand plots: the lines show the fractions of known field pulsars that have been gamma-folded, versus spindown power ˙ E (top) and, when a distance estimate exists, the heuristic gamma-ray flux G h defined in Section 7.1 (bottom). The points show the fractions of folded pulsars for which pulsations have been detected. The vertical dotted line approximates the minimum detectable integral energy flux G 100 . Right-hand plots: Numbers of known and gamma-ray detected field pulsars versus ˙ E . (A quarter of ATNF psrcat pulsars have unmeasured ˙ E .) An earlier version of this figure appears in Laffon et al. (2015).  \n<!-- image -->  \n∝", '3.2.2. Optical, X-ray, and Radio Studies to Refine Targets': "Fewer than 1% of known neutron stars pulse in the optical band, so optical observations of gamma-ray sources might seem unproductive. However, relative to pulsars discovered in radio surveys, the LAT pulsar sample includes many more members of the subclass of interacting binaries-black widows, redbacks, transitional MSPs, and related systems. The companions in these systems can have relatively bright optical magnitudes and (importantly) often exhibit modulation at the orbital period  \nin brightness, color and radial velocity. Consequently, searches for periodic optical sources within LAT source error regions have identified many candidate gamma-ray pulsars. Some of these have been followed up with radio detections, while others have enabled direct gamma-ray periodicity detections by using the precise position and orbital parameters to reduce the vast parameter space to be searched, as detailed in Section 3.4. Others remain strong candidates that are likely powering the gamma-ray sources, but have yet to be confirmed as pulsars (see Section 7). Some examples of the power of optical studies of LAT source regions include:  \n- · Optical and X-ray studies strongly indicated that 1FGL J2339.7 -0531 hosts a black-widow MSP (Romani & Shaw 2011; Kong et al. 2012). After LAT blind searches found no pulsar, Ray et al. (2020) discovered radio pulsations. PSR J2339 -0533 turned out to be a redback with a companion mass outside the parameter space covered by the LAT searches. The resulting ephemeris yielded gamma-ray pulsations as well (Pletsch & Clark 2015).\n- · Optical studies by Romani (2012) led to the first discovery of a millisecond pulsar, PSR J1311 -3430, in a blind search of LAT data (Pletsch et al. 2012a). Detection of radio pulsations followed shortly (Ray et al. 2013). Two further binary MSPs have been discovered in this way: the black-widow PSR J1653 -0158 (Nieder et al. 2020a) from an optical candidate discovered by Romani et al. (2014) and Kong et al. (2014); and the redback PSR J2039 -5617 (Clark et al. 2021) from an optical candidate discovered by Romani (2015) and Salvetti et al. (2015).\n- · Four candidate gamma-ray redbacks identified in optical and X-ray searches were recently identified as radio MSPs by TRAPUM (Halpern et al. 2017; Li et al. 2018b; Swihart et al. 2020; Au et al. 2023), listed in Table 6.\n- · Optical studies of an unidentified LAT source by Strader et al. (2015) revealed an accretion disk in a binary system with a 5 . 4 d period. Camilo et al. (2016) subsequently discovered PSR J1417-4402 using the Parkes radio telescope. First thought to be a redback, its orbit is wide and the wind interactions differ from other spiders. Swihart et al. (2018) argued that it is a possible progenitor of normal field MSP binaries and dubbed it the 'huntsman' as the archetype of an emerging class. Unfortunately, it is difficult to time and a phase-connected ephemeris which would yield gamma-ray pulsations is not anticipated soon.  \nNeutron stars can also power X-ray emission either directly from their surfaces or magnetospheres, or via shock emission when the pulsar wind interacts with a companion or the interstellar medium (Gentile et al. 2014; Marelli et al. 2015; Coti Zelati et al. 2020). The non-thermal shock spectra, sometimes modulated at the orbital period, provide strong signatures. Consequently, X-ray surveys also provide valuable information as demonstrated by varied studies of the error ellipses of unidentified  \nLAT sources. For example, the Swift XRT systematically studied a large number of LAT unassociated source regions 10 . This dataset has been mined using machine learning techniques to identify promising pulsar candidates for radio follow up (Kaur et al. 2019; Kerby et al. 2021).  \nPulsars are steep-spectrum radio point sources (Jankowski et al. 2018). This was used to identify the source (4C 21.53) that was later found to be the first millisecond pulsar (PSR B1937+21; Backer et al. 1982). The availability of radio surveys covering nearly the whole sky was recently harnessed to look for pulsar candidates associated with LAT source regions (Frail et al. 2016, 2018; Bruzewski et al. 2021; Ray et al. 2022). These candidate lists have high precision positions from radio interferometric imaging and well-measured fluxes, enabling both gamma-ray blind searches with reduced parameter spaces and radio pulsation searches that discovered several new MSPs.", '3.3. Targeted Radio Searches': "Shortly after LAT operations began, the Fermi Pulsar Search Consortium (PSC) was organized, bringing together the LAT team and many international collaborators with pulsar searching expertise and access to the world's largest radio telescopes (Ray et al. 2012). The aim was to coordinate radio pulsation searches of LAT sources, and to see which pulsars found in gamma-ray blind searches are truly radio quiet. The PSC has been exceptionally successful at efficiently discovering millisecond pulsars in searches of LAT sources. Beginning in 2010, the TRAPUM collaboration formed to exploit MeerKAT (Stappers & Kramer 2016). It is not formally part of the PSC but includes searches of LAT sources. The FAST radio telescope in China also joined the searches under a separate arrangement. Radio astronomers are provided with early access to LAT source lists as well as the parameters of new pulsars discovered in gamma-ray blind searches (see Section 3.4).  \nLAT gamma-ray source localizations are much better than was possible with EGRET, making it often possible to search an entire 95% confidence error region with a single radio telescope pointing. Of course, the primary beam of a radio telescope depends on the observation frequency and the dish diameter. Thus radio astronomers could choose telescope/frequency combinations optimized for the error regions being searched. Initially, the telescopes were all single-dish (Green Bank Telescope (GBT), Parkes, Effelsberg, Nan¸cay Radio Telescope (NRT), Arecibo Telescope, and the Lovell Telescope at Jodrell Bank), but array telescopes were soon added, the first being the Giant Metrewave Radio Telescope (GMRT), followed by LOFAR and MeerKAT (see Table 4).  \nArray-based telescopes offer increasingly powerful and flexible methods for pulsar searches. Incoherently combining the beams allows for very efficient but less sensitive  \nTable 4. Radio telescopes searching for new pulsars in unidentified LAT sources, and timing discoveries.  \nsearches of large error boxes. Some modern arrays have backends powerful enough to search hundreds of coherent (tied-array) beams arranged to cover a large primary beam, enabling extremely efficient full sensitivity searches of large error boxes (Clark et al. 2023a). Gated imaging (Roy & Bhattacharyya 2013) or tied-array beamforming (Stappers & Kramer 2016; Sanidas et al. 2019) can localize newly-discovered pulsars extremely well, without the need for long-term timing campaigns.  \nA key feature of radio pulsar searches, compared to LAT blind searches, is that it is computationally feasible to maintain sensitivity to short-period binary pulsars even when none of the system parameters are known (as described in Section 3.4, binary MSPs can be found in gamma-ray blind searches when approximate orbital parameters are known from optical studies). Because the observation durations are generally short compared to the orbital period, this can be done through Fourier search techniques that assume a constant acceleration (Ransom et al. 2002), or a constant jerk (Andersen & Ransom 2018). Recently, an analysis of PSC search data  \nusing the jerk search method resulted in the discovery of 11 new MSPs, at least one of which would likely have been missed in a standard acceleration search (Tabassum et al. 2021).  \nAt least 110 new pulsars have been discovered in the radio searches targeting LAT unassociated sources. All except 8 are MSPs - MSPs dominate the high Galactic latitudes where radio surveys have been less thorough. At present, only one of the slow pulsars has been confirmed as a gamma-ray pulsar (PSR J2030+3641; Camilo et al. 2012). Given the many observations, some chance coincidences are expected, especially for telescopes with wide fields of view (particularly arrays of smaller telescopes combined incoherently, like the GMRT, e.g. Bhattacharyya et al. 2022). Table 6 lists 39 MSPs found in the targeted radio searches for which gamma-ray pulsations have not yet been detected. An X indicates 6 where it has been concluded that the MSP is not associated with the target source. Most of the others were discovered after 2016, so the lack of gamma-ray pulsations from these sources is likely due to the lag in obtaining a rotation ephemeris with which to fold LAT data.  \nDeveloping such an ephemeris is an involved process, as the timing model must extrapolate well over the long span of LAT data required for the accumulated signal to become significant. Especially if the position is not well constrained by radio imaging or an X-ray or optical counterpart, then the radio campaign requires a dense set of observations to measure the spin period well and longer intervals (0.5-1.0 yr) to measure the position and ˙ P . Pulsars with orbital period variability may require  \nyears of monitoring. Such follow-up campaigns can require multiple proposal cycles. In most cases, after getting a good timing solution, LAT pulsations have been detected from the MSPs. If pulsations are detected in a fraction of the LAT data, the solution can often be extended to cover the full mission using LAT timing methods (see Section 3.5). We expect that a large fraction of the recent MSPs in Table 6 will be identified as gamma-ray pulsars in the future. Some of the pulsars in Table 5 may also someday reveal gamma-ray pulsations.", '3.4. Blind Periodicity Searches': "In addition to the detection of gamma-ray pulsations from known radio or X-ray pulsars via folding, the LAT data also allow for the discovery of entirely new pulsars through direct searches for their gamma-ray pulsations, without prior knowledge of their spin periods or spin-down rates. These searches complement and contrast the efforts described above in Sections 3.1 and 3.3, providing access to a population of otherwise undetectable pulsars. However, these searches are less sensitive than following-up a known radio or X-ray pulsar, and are extremely computationally demanding: searching for isolated young pulsars/MSPs in a single Fermi -LAT source typically takes a few days/weeks, respectively, using tens of thousands of CPU cores, while a search for a single binary MSP can run for months on thousands of GPUs.  \nOf the pulsars in this catalog, 81 were first discovered in the LAT data, compared to 36 in 2PC. The majority are young pulsars, located near the Galactic plane. Notably, this includes 10 LAT-discovered millisecond pulsars of which three are in binary systems, whereas 2PC contained just one such object. That the discovery rate from direct gamma-ray pulsation searches has not fallen off is perhaps surprising when one considers how the computational cost and sensitivity of these searches increases with the data duration.  \nThe computational cost depends on the dimensionality and volume of the parameter space which must be searched, and the required density of search points covering this volume. At minimum, these searches are 2-dimensional, covering the range of spin frequencies (or periods) and spin-down rates spanned by the known pulsar population. A sub-arcsecond precision position-much better than can be provided by the LATis required to correct for the Earth's orbital Romer delay. Unless one is available from multi-wavelength observations, an additional 2-dimensional grid of sky locations covering the localization region must be searched. The required positional precision in each dimension is proportional to the spin frequency, f . This would make searches for MSPs thousands of times more costly than searching for young pulsars, were this not partially mitigated by the far smaller range of spin-down rate ˙ f over which MSPs are found.  \nWhen searching for pulsars in circular binaries, at least three further parameters are required (the orbital period, phase and radius), with two additional parameters required for eccentric systems. The search grid density in all binary parameters depends on f , again making rapidly spinning millisecond pulsars substantially more costly to find. The binary search parameter space is in fact so large that entirely uninformed searches are as yet infeasible, with searches only becoming possible when a candidate binary system is identified and its orbital parameters constrained via multi-wavelength observations.  \nThe required search-grid density in most dimensions depends on the duration T obs of the data over which the signal is integrated (Pletsch & Clark 2014). In f , and ˙ f , the number of grid points scales linearly and quadratically, respectively, with T obs . For T obs < 1 yr, the number of sky positions which must be searched scales quadratically with T obs , but this saturates at longer timescales. For binary pulsar searches, the required density of orbital period search points also scales linearly with T obs (Nieder et al. 2020a). The result is that the computing cost scales with at least T 3 obs for isolated pulsars and T 4 obs for binaries. This rapid scaling leads to fully-coherent searches of data lasting more than a few months being entirely infeasible.  \nA compromise is achieved using 'semi-coherent' search methods, where the signal powers from many shorter windows (of duration T coh , typically a few weeks long) are combined incoherently (Atwood et al. 2006; Pletsch et al. 2012b). For a semi-coherent search, the numbers of grid points in the frequency and spin-down dimensions are each reduced by a factor of T obs /T coh , while the number of sky positions that must  \nbe searched is reduced by a factor of (1 yr) 2 /T 2 coh (assuming T coh < 1 yr < T obs ). This semi-coherent search is performed as an initial scan over the entire parameter space, with more sensitive fully-coherent stages being used to 'follow-up' a smaller number of semi-coherent candidates to increase their significance.  \nThe compromise is that the sensitivity of a semi-coherent search is far worse than that of a fully-coherent search (e.g. for pulsations from a known pulsar). When folding with a known ephemeris, the signal-to-noise ratio is proportional to T 1 / 2 obs . In a semi-coherent search, this scaling becomes ( T coh T obs ) 1 / 4 . Combining this with the computing cost scaling, we see that the recovered signal-to-noise ratio (a proxy for sensitivity) increases very slowly as the computing cost, C ( T obs , T coh ), increases. For fixed T obs > 1 yr, S / N ∝ C ( T coh ) a , where a = 1 / 8 if positional information is available, or a = 1 / 16 otherwise. For fixed T coh < 1 yr, S / N ∝ C ( T obs ) b where b = 1 / 4 for isolated pulsars, and b = 1 / 8 for binary pulsars.  \nTwo additional effects further limit this sensitivity. First, the vast parameter space which must be searched over introduces a large 'trials' factor, meaning that pulsed signals need to be detected with far higher test statistic values in order to be statistically significant. Where pulsars can be found in single-trial ephemeris folding searches at H = 25, the associated trials factor of a semi-coherent search raises this significance threshold. To be statistically significant, these searches require the coherent power at the fundamental frequency to be P 1 ≳ 120. For typical gamma-ray pulse profile shapes, this corresponds to H ≳ 200. T coh can be chosen such that signals above this level will be detectable in the initial semi-coherent stage.  \nSecond, the accumulation of a pulsed signal can be disrupted by more complex timing behavior such as timing noise, glitching, or orbital period variations in the case of redback binary systems. However, statistical noise still accumulates over the full data duration. For timing behavior that results in a signal being 'well behaved' only over a timescale τ , this means that the signal-to-noise is further reduced by a factor of √ τ/T obs . This latter effect does not occur in isolated millisecond pulsars, whose spindown rates are extremely stable over time, whereas energetic young pulsars, which tend to exhibit the most timing noise and glitches, will be the worst affected. Clark et al. (2017) predicted that the detection rates of isolated MSPs would likely increase relative to those of young pulsars as a result.  \nIn summary, as the LAT data duration increases, these searches become more expensive, and the minimum detectable flux density decreases more slowly than that of an ephemeris-folding search. The ratio between these sensitivities depends on the background photon flux, and hence on the location of the targeted source in the sky, but varies between a factor of ∼ 3-4 for the ongoing Einstein@Home pulsar surveys (with T obs = 10yr, T coh = 2 22 s). The difference in sensitivities between these two methods is emphasized by the difference in the gamma-ray fluxes of the radio-loud and radio-quiet populations shown in Figure 6.  \nTo tackle these growing challenges, search efforts require increasing computing resources and methodological developments to increase efficiency and sensitivity. Fortunately, the available computing power does increase over time as technology progresses, and as search efforts have shifted to larger computing systems. The first successful searches in the early Fermi mission either searched the center of the localization region for each source, or targeted the location of a plausible X-ray counterpart (Abdo et al. 2009b; Saz Parkinson et al. 2010). After the first two years of the Fermi mission, these discoveries became less frequent. Newer searches (Pletsch et al. 2012b,c) mitigated dependence on event selection criteria and source localization by weighting events and shifting to large computing clusters to enable searching over a grid of positions. These searches were later migrated to the Einstein@Home distributed volunteer computing system, allowing for longer T coh and hence more sensitivity (Pletsch et al. 2013; Clark et al. 2015, 2016, 2017, 2018).  \nAdditional sensitivity can be gained not just through computing power, but also through careful study. While semi-coherent search methods were developed prior to Fermi 's launch (Atwood et al. 2006), improvements to their sensitivity and efficiency have been made throughout the LAT mission (Pletsch & Clark 2014; Nieder et al. 2020a). Improvements to the LAT event reconstruction and background rejection increase the recoverable signal strengths (Atwood et al. 2013; Bruel et al. 2018), as do higher-level processing developments such as photon weighting (Kerr 2011; Bruel 2019). An unbinned likelihood method improved LAT's sensitivity at all time scales, useful for a broad range of astrophysical targets in addition to pulsars (Kerr 2019).  \nExtensive multi-wavelength searches for possible counterparts at which to target pulsation searches continue to lead to new discoveries. To give some recent examples: targeted searches of X-ray sources discovered energetic gamma-ray pulsars within two supernova remnants (PSRs J1111 -6039 and J1714 -3830), and searches targeting steep-spectrum radio sources associated with LAT sources discovered two isolated MSPs (Frail et al. 2018).  \nSearches for periodic optical and X-ray counterparts to pulsar-like LAT sources have identified several promising black-widow and redback candidates, detailed in Section 7. The precise orbital constraints from long-term optical observations are key inputs for binary searches, although they remain the most computationally demanding projects. Nevertheless, the development of efficient search grid designs (Nieder et al. 2020a) and a GPU-accelerated 11 search algorithm deployed on Einstein@Home has resulted in two binary MSP discoveries (Nieder et al. 2020b; Clark et al. 2021), in addition to PSR J1311 -3430 (Pletsch et al. 2012b), which was included in 2PC.  \nThe complementary nature of the blind-search and ephemeris-folding populations is highlighted by the number of radio detections of blind-search pulsars. In 2PC, there were radio detections of just four pulsars discovered in gamma-ray searches.  \nThis fraction remains small in this catalog, only 11 radio detections for 81 gammaray blind-search pulsars, despite deep radio follow-ups. Of those detected, most have radio flux densities lying below typical radio survey detection thresholds, which are above the 30 µ Jy limit we use to define 'radio quiet', shown in Figures 6 and 7. This is perhaps unsurprising: most pulsars discovered in gamma-ray searches are young pulsars located in the Galactic plane, which has been extensively surveyed by radio telescopes, and therefore we would expect most radio-loud pulsars in this region to be known. Most gamma-ray pulsars found off the Galactic plane are millisecond pulsars, of which most are found in binary systems, and therefore inaccessible to gammaray pulsation searches unless orbital constraints from an optical counterpart can be derived. A notable addition since 2PC are the 6 MSPs that remain undetected in radio observations, including the black-widow PSR J1653-0158 (Nieder et al. 2020b). In such spider binaries, the dense ionized wind can obscure the pulsed radio signal for much, or at times all, of the orbit, hindering radio detection.", '3.5. Timing Pulsars using LAT data': "The timing validity range of a recently discovered radio or X-ray pulsar will likely not cover the entire LAT data set, limiting the statistics available for profile characterization, and limiting the discovery potential for glitches, profile changes, spin-ups (Marshall et al. 2015) or other rare behavior. Often, pulsar timing using LAT data allows an accurate phase-connected ephemeris covering the entire mission duration.  \nLAT pulsar timing has significantly improved, with new techniques to handle timing noise and faint pulsars developed by Ray et al. (2011) and Kerr et al. (2015). They use pulse time of arrival (ToA) estimation and time-domain models of timing noise. For faint pulsars, ToA estimates for short data segments are not possible, and we adopted the unbinned methods used in Ajello et al. (2022). They optimize the timing model by maximizing the likelihood, using PINT (Luo et al. 2021) to compute the rotational phase for each photon. Parameter uncertainties are derived from the likelihood curvature near the optimum. Particularly useful for MSPs, we have been able to estimate some precise proper motions directly from LAT data, necessary for the Doppler corrections presented in Section 4.3.  \nWhen both radio and LAT timing of comparable quality are available, we use the radio ephemeris. LAT timing was used for 142 of the pulsars in this catalog. Twothirds of those were made with the ToA method. The fraction of timing solutions using the unbinned method is increasing, as analysts adopt the tool, as more difficult pulsars are added to the sample, and as the datasets lengthen. The ephemerides we used to fold the 294 pulsars are included in the supplementary on-line data, at https:// fermi.gsfc.nasa.gov/ssc/data/access/lat/3rd PSR catalog/ (password protected into August 2023), and at https://fermi.gsfc.nasa.gov/ssc/data/access/lat/ephems/. See Appendix D.  \nTable 5 . Galactic pulsars colocated with LAT sources, i.e. within f = ∆ θ/r 95 < 1 . 2 of a 4FGL-DR3 source, and having efficiency ϵ = L γ / ˙ E < 10, assuming the pulsar distance (73 matches were overluminous). For association classes 'psr' or 'msp' we list the pulsar regardless of the efficiency. These pulsars show no gamma-ray pulsations, although many have not been gamma-ray folded, indicated by no ˙ E value or by * in the PSR column. ∆ θ is the angular separation between the pulsar and the source, and r 95 is the 95% confidence level semi-major axis of the LAT error ellipse. E peak is the peak of the 4FGL spectral energy distribution when the PLEC fit was reported in 4FGL-DR3. The last column lists the 4FGL-DR3 association for the LAT source, 'sic' means that it matches the pulsar name. V indicates Variability Index > 24 . 7 . Pulsars with a 'g' suffix come from the FAST Galactic plane pulsar survey (GPPS) (Han et al. 2021).  \nrd LAT Pulsar Catalog  \n3 Fermi  \nTable 5 (continued)", '4. THE GAMMA-RAY PULSARS': "The search methods yielded gamma-ray pulsations for 294 pulsars. Table 1 breaks down the numbers by different types. Here we compile their radio fluxes, distances, and, for the MSPs, their proper motions, used to Doppler-correct their spindown properties (Shklovskii 1970). Figure 3 shows the sky distribution of the Fermi pulsars. The LAT's uniform sky coverage is important for detecting MSPs, as these relatively nearby sources are approximately uniformly distributed in Galactic latitude.", '4.1. Radio Intensities': "We tabulate measured or derived 1400 MHz flux density values ( S 1400 ) in Tables 2 and 3. ATNF psrcat provides flux densities S ν from among 27 different radio frequencies ν ; thus as a starting point, we used psrqpy (Pitkin 2018) to extract all S ν values from psrcat , as well as the radio spectral index, SPINDX . When available, we use S 1400 . Otherwise, we calculate S 1400 = S ν ( 1400 ν ) α with α = SPINDX when present, and α = -1 . 7 otherwise. We use the S ν value corresponding to the frequency ν closest to 1400 MHz. We used α = -1 . 7 in 2PC, a compromise among published values: -1 . 6 from Lorimer et al. (1995), -1 . 8 ± 0 . 2 from Maron et al. (2000), and -1 . 60 ± 0 . 03 from Jankowski et al. (2018). For MSPs, Spiewak et al. (2022) found -1 . 92 ± 0 . 06 with a standard deviation of 0 . 6, and Dai et al. (2015) found -1 . 81 ± 0 . 01. The small formal uncertainties are potentially misleading: for example, the value from Jankowski et al. (2018) is for a sample of 400 pulsars with a log-normal distribution of standard deviation 0 . 54, the distribution being sensitive to e.g. the frequency range used, yielding  \nTable 6. 39 MSPs discovered in radio searches of LAT unidentified sources, yet to show gamma-ray pulsations. Six pulsars with X in the 4FGL column are unrelated to the gamma-ray source that guided the radio search. Localization of the two pulsars with no 4FGL name is as yet inadequate for source matching. Distances use DM and the YMW16 model (Yao et al. 2017). Three unpublished pulsars were found using methods as in Camilo et al. (2015).  \na Search targeted uncatalogued LAT source 605N-0275 (T. Burnett, priv. comm.)  \nNote - The Discovery column takes ' Fermi ' entries from http://astro.phys.wvu.edu/GalacticMSPs/ GalacticMSPs.txt. See also http://www.trapum.org/discoveries/  \nFigure 6. The radio flux density at 1400 MHz versus the integral energy flux above 100 MeV. Gamma-ray and radio fluxes are essentially uncorrelated. Dark borders on the triangle symbols indicate MSPs discovered in deep radio searches of previously unidentified LAT sources by the PSC (see Section 3.3). The horizontal line at 30 µ Jy shows our convention to distinguish radio 'loud' from radio 'quiet' pulsars.  \n<!-- image -->  \nthe spread in average indices. Using multiple S ν values to fit a power law gives S 1400 values generally compatible with those found above, and possibly more accurate in some cases. However, discrepant S ν values lead to unrealistic results in some cases, hence we use the simple prescription described above.  \nThe psrcat S ν values have accumulated over decades, using a variety of instruments and analysis conventions. Frequency extrapolation introduces additional biases. Fortunately, since v1.68 psrcat includes S 1400 values from two recent studies. Spiewak et al. (2022) used MeerKAT to obtain S 1400 for 189 MSPs 12 , and Parent et al. (2022) analyzed Arecibo data from the PALFA survey for 93 MSPs. We also replace the Hessels et al. (2011) values with the updated results reported by Bangale et al. (2023). For a small number of pulsars, we obtained as-yet-unpublished S ν measurements, noted in the Tables.  \nIn addition to pulsars which are well characterized from radio surveys, those discovered in blind searches are typically searched deeply for radio pulsations (Camilo et al. 2009; Abdo et al. 2010d; Saz Parkinson et al. 2010; Ray et al. 2011, 2012; Pletsch  \net al. 2012b; Nieder et al. 2020b). The 81 blind discoveries were followed by 11 radio detections, 6 of which qualify as 'radio-loud', with S 1400 > 30 µ Jy (see below). Radio measurements are not yet available for recently discovered LAT pulsars.  \nFigure 6 shows that radio flux densities are essentially uncorrelated with the integrated gamma-ray energy fluxes, G 100 . Population studies aim to constrain a powerlaw model of radio luminosity in which L r ∝ ˙ E s . Johnston & Karastergiou (2017) favor s = 1 / 4 but recognize that a range of s values match observations, because of the large dispersion. Posselt et al. (2023) find an even smaller index. Since L γ ∝ √ ˙ E (see Figure 23), flux correlation between the two wavelengths is naturally expected. However, the fraction of neutron star spindown power carried off by radio waves is thousands of times lower, or more, than that of gamma rays, and radio emission is so weakly connected to the braking mechanism that the correlation is invisible.  \nFigure 7 shows how the radio flux densities of the gamma-ray pulsars compare with the overall pulsar population, as well as their lack of distance dependence. As in 2PC, we define a pulsar as 'radio-loud' if S 1400 > 30 µ Jy, and 'radio-quiet' if the measured flux density is lower. The horizontal lines in Figures 6 and 7 show the threshold. Also as in 2PC, the diagonal line in Figure 7 shows an alternate threshold at pseudo-luminosity 100 µ Jy-kpc 2 .  \nWhether a radio or a gamma-ray beam (or both) sweeps through our line of sight from Earth depends on the beam shapes and orientations. Sieber et al. (1975) already showed that radio beam size often increases with decreasing frequency. Grießmeier et al. (2021) searched 27 radio-quiet northern gamma-ray pulsars at 150 MHz but obtained no new detections. This result is compatible with the current understanding of these pulsars' orientations based on their emission geometries, which are derived from LAT light curves and plotted in Johnston et al. (2020); Grießmeier et al. (2021), and which are consistent with the radio beam not crossing the Earth's line of sight.", '4.2. Distances': "Converting integrated energy fluxes G 100 to emitted luminosities L γ = 4 πd 2 f Ω G 100 requires the source distances d . (We set the beaming fraction f Ω = 1, see Section 6.7.) Knowing pulsar distances also allows mapping neutron star distributions relative to the Galaxy's spiral arms, as in Figure 8, and evaluating their scale height above the plane. The various methods to estimate pulsar distances offer widely differing accuracies. Tables 7 and 8 list the distances that we use, the estimation methods, and the references.  \nDispersion measure (DM) is the column density of free electrons along the path from Earth to the pulsar, in units of pc cm -3 . The electrons delay the radio pulse arrival by ∆ t = DM( pν 2 ) -1 where ν is the observation frequency in MHz and p = 2 . 410 × 10 -4 MHz -2 pc cm -3 s -1 . Given a model for the electron density n e in the various structures of our Galaxy, integrating DM = ∫ d 0 n e dl along the line of sight dl yields the distance d for which DM matches the radio measurement. We use mainly the YMW16  \nmodel (Yao et al. 2017), but often consult the NE2001 model when studying particular lines of sight (Cordes & Lazio 2002). Significant discrepancies between model and true distances occur along some lines of sight, more than three times the probable true distance in the worst cases. Examples are PSR J0614 -3329 (Bassa et al. 2016a); PSR J2021+3651 (Abdo et al. 2009c); PSR J0248+6021 (Theureau et al. 2011); and PSR J0908 -4913 (Johnston & Lower 2021). Figure 4 of Hou et al. (2014) illustrates one way to identify pulsar lines of sight that have DM model peculiarities.  \nRadio spectra can show Doppler-shifted neutral hydrogen ( Hi ) absorption or emission lines due to clouds along the line of sight. A Galactic rotation model converts the Doppler-inferred radial velocity to a cloud distance estimate, and the pulsar distance is constrained if it appears to be in or between clouds. This 'kinematic' method is tagged with a K in the tables. The X tag indicates cases where X-ray absorption below 1 keV yields a column density that can be compared with one obtained from 21 cm radio surveys.  \nTo obtain the distances we first query the ATNF psrcat DIST parameter, which is most often DIST DM, the dispersion measure distance using YMW16. We adopt a 40% uncertainty on DIST DM as per Figure 12 and Table 5 of Yao et al. (2017). When psrcat has a DIST A value, obtained by another method, we use that instead  \nFigure 7. Radio flux density at 1400 MHz versus pulsar distance. The horizontal line at 30 µ Jy shows our convention to distinguish radio 'loud' from radio 'quiet' pulsars. The diagonal line shows a threshold in pseudo-luminosity of 100 µ Jy - kpc 2 . The pulsars at lower-right are assigned distance limits along the Milky Way's rim in Figure 8. DM for the pulsars on the line at 25 kpc saturate the YMW16 model.  \n<!-- image -->  \nFigure 8. Gamma-ray pulsar positions projected onto the Milky Way model of Reid et al. (2009). The MSPs that appear coincident with the Galactic center, PSRs J1823 -3021A and J1835 -3259B in the globular clusters NGC 6624 and 6652, lie well above the Galactic plane. Distance uncertainties are not shown, but can be large for the more distant objects. The squares with arrows indicate the lines of sight toward pulsars for which no distance estimates exist, placed at the distances where 95% of the electron column density has been integrated in the YMW16 model. Other symbols are as in Figure 2.  \n<!-- image -->  \nand propagate the method and reference to our tables. The exception is that we exclude the distances deduced by Wang (2011) using the LAT G 100 values and the correlation between L γ and ˙ E shown in Section 6.7, since they would bias our attempts to better determine that correlation. If the psrcat parallax parameter PX is present, DIST A = 1/PX and we use the reference for PX. Parallax is measured only for relatively nearby pulsars, using X-ray or optical images, radio interferometric imaging, or accurate timing 13 . When converting parallax measurements to distances, the LutzKelker (LK) effect is an overestimate of parallax values, and hence an underestimate of distances, due to the larger volume of space traced by smaller parallax values (Lutz & Kelker 1973). We use the LK-corrected distance estimates determined by Verbiest et al. (2012). Good estimates of the distance to an optical companion in a binary system are more frequent than in 2PC because we have more spider MSPs than before, and because of the wealth of information provided by Gaia (Gaia Collab.  \net al. 2018). For two of our MSPs, PSRs J1227 -4853 and J1431 -4715, Jennings et al. (2018) find distances using Gaia companion parallaxes. They use a spatial distribution prior specific to Galactic pulsars, i.e., not isotropic, which means that the LK-like correction brings some pulsars closer to Earth when they are towards the anti-center, as opposed to farther for standard LK. For these two MSPs, we simply use 1/PX, the same within uncertainties as the Jennings et al. (2018) distances.  \nFor several pulsars, measured DM exceeds DM max , the maximum modeled by Yao et al. (2017) for that line of sight, by between 1 and 50%. The YMW16 model therefore returns the distance to the Galaxy's edge. Modeled DM saturates much closer to Earth, increasing slowly thereafter. For example, towards PSR J1327 -0755, 95% of DM max occurs at 3 . 7 kpc, compared to 25 kpc returned by the code, and in better agreement with the 1 . 7 kpc NE2001 distance. Therefore, for pulsars with DM > DM max , we use the YMW16 distance for 95% of DM max . Other pulsars cited in this catalog for which this occurs are PSR J1102+02 (Table 6), and PSRs J1302 -03, J2051+4434g and J2052+4421G (Table 5).  \nFor 63 pulsars we have no DM measurement or other distance estimator. We again calculate the YMW16 distance for 95% of DM max , but use it as an upper limit. The method is 'DMM' in Tables 7 and 8, and these radio-quiet pulsars are illustrated in Figure 8.  \nBoth Espinoza et al. (2013) and 2PC discussed the problem of PSR J0610 -2100 having gamma-ray efficency η = L γ / ˙ E > 1, using the DM distances of 3 . 3 and 4 kpc. van der Wateren et al. (2022) used LOFAR timing and optical modelling for this system to find d = 2 . 2 ± 0 . 7 kpc. Espinoza et al. (2013) point out that the line of sight to the pulsar grazes the electron overdensities bordering the local bubble, in both electron models. A slight extension of the modeled bubble yields the increase in modeled DM needed to reconcile the distances. The closer distance for PSR J0610 -2100 leads to η ≈ 50%. It also reduces the pulsar's transverse velocity to v T ≈ 220 km s -1 , so that this MSP is no longer among the fastest known.  \nTable 7 . Distance estimates for young LAT-detected pulsars.  \nTable 7 continued on next page  \n3 rd Fermi LAT Pulsar Catalog  \nTable 7 (continued)", 'T able 8 c ontinue d on next p age': "T able 8 (c ontinue d)  \n<!-- image -->  \nT able 8 (c ontinue d)  \n<!-- image -->  \nT able 8 (c ontinue d)  \n<!-- image -->  \nNote -Columns 2 and 3 giv e the distances for the 144 MSPs detected with the F ermi LA T, and the metho d used to find them: P -parallax; LK -corrected for the Lutz-Kelk er effect; DM -from disp ersion measure using the YMW16 mo del ( Y ao et al. 2017 ) ; O -optical and other metho ds. F or DM, the references in Column 4 giv e the DM measuremen t. DMM is a limi t on the pulsar distance, corresp onding to 95% of the maxim um YMW16 DM for that line of sigh t. Prop er motion measuremen ts µ are listed in C olumn 5, tak en f rom the references in Column 6, often the 'Same' as in Column 4. 'LA T' refers to gamma-ra y timing, describ ed in Section 3.5 . 'NR T', 'JBO', 'PKS', and 'GBT ' indicate that the measuremen ts w ere m a de for this w ork using data from the Nan¸ ca y , Lo v ell, P a rk es and Green Bank radio telescop es. Columns 7 and 8 list the con tributions of the Shklo vskii effect ˙ P shk and of acceleration due to the Galactic r otation, ˙ P gal to the in trinsic spindo wn rate ˙ P in t . The last colum n s list the in trinsic spi n do wn p o w er ˙ E in t = ˙ E (1 -ξ ). Se e text for ho w cases with a missing or insignifican t µ measuremen t, or a distance lim it, are handled.", '4.3. Doppler corrections': "Many pulsar characteristics depend on the intrinsic spin period P int and spindown rate ˙ P int . The Doppler shift of the observed period is P = (1 + v R /c ) P int , where v R is the pulsar's radial velocity along the unit vector n from the solar system. The Doppler correction to ˙ P is obtained by differentiating the equation and separating the effects of the system's proper motion (Shklovskii 1970) from the acceleration due to Galactic rotation (the treatment by Damour & Taylor (1991) is particularly clear). The result is  \n˙ P int = ˙ P -˙ P shk -˙ P gal (3)  \nwith  \nand  \n˙ P shk = 1 c µ 2 dP = k ( µ masyr -1 ) 2 ( d kpc ) ( P s ) (4)  \n˙ P gal = 1 c n · ( a 1 -a 0 ) P (5)  \nwhere k = 2 . 43 × 10 -21 for pulsar distance d and proper motion transverse to the line of sight µ . In 2PC we obtained the accelerations a 1 of the pulsar and a 0 of the Sun using the Galactic potential model of Carlberg & Innanen (1987) and Kuijken & Gilmore (1989) . Here, we instead use GalPot 14 with the PJM17 best.Tpot parameter set (McMillan 2017).  \nThe corrections are negligible for the young gamma-ray pulsars, which all have ˙ P > 10 -17 . However, for MSPs ˙ P int can differ noticeably from observed ˙ P . Expressing the correction magnitude as  \nξ = ( ˙ P shk + ˙ P gal ) / ˙ P = ( ˙ P -˙ P int ) / ˙ P (6)  \nleads to ˙ E int = ˙ E (1 -ξ ). For µ ≲ 10 mas yr -1 and d ≲ 1 kpc, ˙ P shk and ˙ P gal have comparable magnitudes and opposite signs. For larger µ 2 d , ˙ P shk dominates.  \nTable 8 lists d, µ , their uncertainties δd, δµ , and the resulting corrections. When we have no µ measurement and only a limit for d , we leave the last columns empty. Parentheses around the µ ± δµ value indicate that the measurement is insignificant, δµ/µ > 3, and we express the remaining columns as limits calculated using µ max < µ +2 δµ . In parentheses in the ˙ P gal column are the minimum and maximum values for that direction for distances from 500 pc to the limit. Given a distance but no proper motion value, we calculate ˙ P gal and limits on ξ , using ˙ P shk ≥ 0. There are 8 such MSPs, of which PSR J0605+3757 is listed first. In the same spirit, e.g. PSR J1335 -5656 has measured proper motion but unknown distance and we combine the limit for ˙ P shk with the extremum from the ˙ P gal range.  \nTable 8 lists a single MSP with well-measured µ but (1 -ξ ) < 0, implying implausible spin-up, ˙ E int < 0. It is PSR J1024 -0719, with a binary companion in an atypical  \nmulti-century orbit that leads to mis-measured ˙ P int (Bassa et al. 2016b; Kaplan et al. 2016; Guillemot et al. 2016). Six other MSPs have ˙ E int ≲ 2 × 10 33 erg s -1 , some with nominally small uncertainties. The lowest is ˙ E int = 9 +21 -8 × 10 32 erg s -1 for PSR J1946+3417. Improved precision and a larger sample will clarify the gamma-ray 'deathline' for MSPs, below which emission in the LAT range ceases. A useful tool to understand the uncertainty on ˙ E int are plots of µ vs d showing the constraints suggested by the neutron star's transverse velocity v T = µd , minimum ˙ E int , and so on. A powerful constraint is to require ˙ E int > L γ (Section 6.7 addresses luminosity L γ ). Examples of µ vs d plots, or similar, are Figure 11 in 2PC (which includes PSR J1024 -0719), Figure 1 in Guillemot et al. (2016), and Figure 4 in Clark et al. (2018).", '5.1. Gamma-ray and Radio Light Curves': "Figure 9 shows an example of gamma-ray light curves in different energy bands for the MSP PSR J0030+0451. Figure 4 and Appendix B show further examples, overlaid with radio profiles when available. The supplementary online material includes all the pulse profiles constructed for this work.  \nAfter preparation according to Section 2, the photon arrival times were converted to rotational phase using timing models obtained from pulsar timing observations conducted in radio or X-rays, or by direct timing of the LAT photons, as described in Section 3. In a number of cases, the timing models' validity interval only partially overlapped with the time range of the LAT dataset. To avoid contaminating off-pulse phase regions by extrapolating inaccurate models, we restricted the LAT datasets to the validity intervals in the cases of (1) timing solutions containing high-order ( ≥ 3) frequency derivatives suggestive of strong timing noise; (2) 'IFUNC' or 'WAVE' terms (also used to mitigate timing noise); or (3) 'ORBIFUNC' or high-order ( ≥ 2) orbital frequency derivatives for pulsars in binary systems, or black widow pulsars that are often subject to timing instabilities. For other objects, generally stable MSPs, all LAT photons were used and we checked that the pulse phases do not drift outside the ephemeris validity intervals (see Figure 4). Additionally, we conservatively excluded 10 days of data before and after pulsar glitch epochs, to avoid contamination of offpulse phase intervals caused by imperfect glitch modeling. At the top of all the plots like Figure 9, t int gives the number of years of LAT data retained.  \nPulse profiles as in Figures 4 and 9 were generated by computing histograms of the rotational phases ϕ j , weighted by the photon weights w j . The weighted phase histograms repeat over two pulsar rotations to clarify structures near phase 0. The uncertainty for the i th histogram bin containing N i photons was calculated as σ 2 i = ∑ N i j =1 w 2 j + ( max j ∈ [1 ,N i ] w j ) 2 . The first term is the sample variance, while the second term adds in quadrature the maximum possible variance for a bin with one photon. This is an ad hoc correction for bins with few counts, where by chance the observed sample variance may be far below the typical value: this prescription provides  \nmore realistic error bars than that used in 2PC. The number of bins in the weighted histograms was chosen based on the H -test TS above 0.1 GeV: 25 bins for pulsars with H < 100, 50 bins when 100 ≤ H < 1000, 100 bins when 1000 ≤ H < 10000, and so on.  \nDashed horizontal lines in the figures show the background contribution from neighboring point sources and from diffuse emission. As in 2PC Section 5.1, it was estimated as b = ( ∑ N j =1 w j -∑ N j =1 w 2 j ) /N b in where N is the number of sources, excluding the pulsar, and N b in is the number of histogram bins. The background uncertainty is dominated by systematic errors in the diffuse background normalization, estimated by changing the background normalization by ± 6%, also shown as dashed horizontal lines. A pulsar wind nebula (PWN) may surround the pulsar, but with a power-law spectrum extending to high energies. If omitted from the source model, it appears as an excess above background, especially at higher energies. Section 6.2 lists the pulsars for which we added a PWN component to our source model, in addition to the Crab and Vela. 2PC Section 7 describes a search for off-pulse emission, which measured off-pulse spectra for pulsars with an excess above the background level.  \nFor gamma-ray pulsars with no radio emission, we rotated the pulse profiles such that the phase of the first emission component (as determined from the light curve fit above 0.1 GeV, described in the next Section) is at phase 0.1 (see Figure B-1 in the Appendix for an example). When radio profiles are available, a radio fiducial phase at infinite frequency was determined (see below) and placed at phase 0. The gamma-ray light curves were rotated accordingly, to preserve the relative radio/gamma-ray phase offset, including DM delays (see Section 4.2). We determined the uncertainty on the radio/gamma-ray phase offset stemming from the DM uncertainty neglecting DM time-derivative terms, and show it in the figures with an horizontal error bar at topright. Figure B-3 in the Appendix shows a case where a large DM uncertainty leads to a large relative phase uncertainty. For some pulsars the DM and its uncertainty were not determined during the radio timing analysis; for these objects we used external DMinformation, generally from the ATNF pulsar catalog or separate analyses. Those cases are highlighted with a marker underneath the phasing error bar, as illustrated in Figure B-4. In a very few other cases, the radio data analysis did not yield the information necessary to determine the relative radio/gamma-ray phase offset. For these light curves we indicated the arbitrary relative phase offset with a ⋆ symbol.  \nIdentifying the 'fiducial phase' at which the magnetic pole of the neutron star crosses the line of sight is of interest. For some pulsars, this can be inferred from radio polarimetry via the rotating vector model (Radhakrishnan & Cooke 1969), but we do not consider such models here. As noted in 2PC, for pulsars with a single, symmetric radio peak, the peak maximum can be chosen as an estimator for the fiducial phase. However, many radio profiles in our sample do not consist simply of single, symmetric peaks, so that the choice of the fiducial phase is not always trivial. We adopted the following 'Gaussian smoothing' approach which reduces complex  \nemission structures by replacing them with the best-fitting Gaussian equivalent. Each radio profile was inspected visually and one peak or peak complex was chosen as the phase reference. This reference peak or set of peaks was then fitted with a single asymmetric Gaussian curve, or with two asymmetric Gaussian curves, depending on the morphology. Single peaks or 'patchy' profiles with more than two components were fitted with single asymmetric Gaussian functions. Peaks with two components at their edges were fitted with two asymmetric Gaussians. In the first case, the fiducial phase was simply defined as the best-fit peak phase, and in the second case the fiducial phase was defined as the average of the two best-fit peak phase values.  \nWhile the above procedure for defining the radio fiducial phase works well in most cases, a number of radio profiles in the sample are too complex for the procedure to produce an estimate. We inspected the radio profiles and classified them into four categories:  \n- 1. radio profiles for which one can readily define one peak as the main peak, which has a simple, single or double structure. The procedure for defining the radio fiducial phase should have produced robust results.\n- 2. radio profiles for which one can readily define one peak as the main peak, whose structure is complex (i.e., patchy profiles). The procedure is believed to be robust.\n- 3. one can define a radio peak as the main one, but cannot readily determine whether it is isolated, or part of a broader structure with other peaks. The radio fiducial phase estimate may not be robust for these objects.\n- 4. the profile morphology is too complex for the fiducial phase to be unambiguously chosen, based on the profile alone.  \nThe categories, or 'radio classes', assigned to the different pulsars are listed in Tables 9 and 10. The fiducial phases were finally placed at phase 0 in the radio and gamma-ray pulse profiles, and the gamma-ray phases were rotated accordingly.", '5.2. Gamma-ray Light Curve Fitting Above 0.1 GeV': 'Fits of the LAT weighted histograms above 0.1 GeV are displayed in the Figures as magenta lines (see e.g. Figure 9). We followed a similar approach as in 2PC for fitting the LAT pulse profiles, by representing the gamma-ray light curves as wrapped probability density functions (pdfs) of ϕ ∈ [0 , 1), and constructing the pdfs for the individual pulsars as linear combinations of N peaks unimodal distributions:  \nf ( ϕ ) = N peaks ∑ i =1 n i g i ( ϕ ) +   1 -N peaks ∑ i =1 n i   . (7)  \nIn the above expression, n i terms represent normalization factors, and the g i terms represent individually normalized distributions and the term in parentheses is a uni-', '5.3. Energy-resolved Light Curve Fitting': 'The best-fit profiles above 0.1 GeV were then used to fit gamma-ray pulse profiles in narrower energy ranges. To this end, we divided the LAT datasets into logarithmically-spaced energy intervals, and fitted the pulsar pulse profiles in each energy interval independently. The number of energy intervals per decade was determined by the H -test TS above 0.1 GeV: one per decade for pulsars with H -test TS values between 100 and 1000, two per decade when 1000 ≤ H ≤ 10000, and so on. Profile fits are again done only if H -test > 100 in an energy band: 167 pulsars have at least two such intervals. Profile parameters were allowed to vary within the same ranges as for the fits over the entire energy range. By default, the number of components for a given pulsar was the same as for the fit over the entire energy range. However, profile fits in the narrow energy ranges often require fewer components than for the entire range. We fitted the profiles with N peaks > 2 components in the individual energy ranges with all possible component combinations, and selected the combination favored by the BIC in most energy intervals. For most objects the number of peaks remained the same, for the others the third component (typically, bridge emission) was discarded. PSR J0007+7303, shown in Appendix B, is an example. The peak shapes, Gaussian or Lorentzian, are the same in the sub-bands as in the energy-integrated profile.', '6. PHASE-AVERAGED ENERGY SPECTRA': "Two significant improvements enhance our pulsar spectral results relative to 2PC. The first is simply that we have more pulsars, and more photons per pulsar. The second is a change in our expression for the exponentially cutoff power law spectrum, a characteristic of pulsars in the GeV domain. We replace the classic sub-exponentially cutoff power-law model ('PLEC') with the new 'PLEC4' expression introduced in 4FGL-DR3. Degeneracies between parameters in PLEC4 are reduced and the likelihood fits converge better, yielding smaller uncertainties. However, some subtleties of the new expression merit the detailed discussion below.  \nIn this Section we first review the motivations behind the exponentially cutoff power law spectrum, highlighting the differences between PLEC used in 2PC and PLEC4 used here. We then describe our results, using the 4FGL-DR3 phase-integrated spectral analysis available for 255 pulsars, and explain some refinements of the spectral analyses that help characterize the pulsars' emission properties, finally yielding 255 pulsars with spectral analyses in this catalog. Tables 13 and 14 list parameter values, also available in the electronic version of the catalog. Section 6.5 highlights some features of individual pulsar results. Section 6.6 addresses pulsar flux stability, and in Section 6.7 we combine integrated energy fluxes and distances to obtain gamma-ray luminosities and efficiencies.", '6.1. Pulsar Spectral Shapes': "As we will show, the general GeV pulsar spectrum is peaked around E p = 1 . 5 GeV (see also Figure B-9 in Appendix B). It is natural to relate this to the synchrotron spectrum radiated by an ultrarelativistic primary particle (Pacholczyk 1970; Dermer & Menon 2009)  \nE 2 dN/dE ≡ S ( E,E c ) ∝ E 4 3 exp[ -E/E c ] ≡ E 2 -Γ sr exp[ -E/E c ] . (13)  \nHere, the asymptotic spectral index is Γ sr = 2 / 3 and the cutoff energy E c depends on the specific radiation mechanism, e.g. the field line curvature radius ρ with E c ∝ γ 3 e /ρ for curvature radiation within the light cylinder (Romani 1996). Farther out, for example in current sheets, the curvature of the particle trajectory may be quite different. Other emission mechanisms, such as inverse Compton, can also produce peaked spectra. In general, pulsars accelerate particles to a range of energies and the observed spectrum is a broader superposition of such basic SED shapes. A main goal of gamma-ray spectroscopy is to invert this superposition and reveal properties of the underlying particle energy distribution. In a curvature radiation scenario, e.g., the observed E p will roughly track the bulk maximum values of E c , while the width of the peak and the low-energy spectrum ( E ≪ E p ) encode information about the range of E c . A nearly monoenergetic distribution of E c will produce a narrow SED peak and a low-energy S ( E ) ∝ E 4 / 3 , while a broad range of E c will produce a broad gamma-ray peak and Γ 0 of ∼ 1-2 (Romani 1996).  \nMost previous analyses of LAT pulsar spectra have used the PLEC model:  \ndN dE ≡ N ( E ) = N 0 ( E E 0 ) -Γ 0 exp [ -( E E c ) b ] , (14)  \nwhich is closely related to a monoenergetic curvature emission or synchrotron radiation spectrum, but with the cutoff strength modulated by the parameter b . Blends of particle energies and averaging over pulse phases both tend to broaden the peak, yielding 'sub-exponential' shapes with b < 1 (Abdo et al. 2010h). One innovation since 2PC, adopted beginning with the 8-year 4FGL catalog and justified therein, is to fix b = 2 / 3 instead of b = 1 for pulsars with lower TS, as this choice better matches the observed spectra.  \nA good spectral model should have the flexibility to capture the features evoked above, but in most cases measuring them all is not possible. The LAT sensitivity peaks at ∼ 1 GeV and the SED slope at 100 MeV, 2 -Γ 100 , can only be measured for the brightest pulsars. Furthermore, it is only in the cases of sharply-peaked spectra which quickly transition to a low-energy power law that Γ 100 may be a good proxy for the asymptotic low-energy slope Γ 0 . Above 1GeV, pulsar spectra cut off quickly, making it difficult to measure the cutoff shape and strength. We show below that for most pulsars, the only reliable features which can be measured are the position and width of the SED peak. Fortunately, it is these features which strongly correlate with underlying pulsar properties like ˙ E , so they are likely to be the most important in capturing information about the pulsar emission mechanism.  \nThe PLEC formula suffers from large covariance between the model parameters (Abdollahi et al. 2022), magnifying the effect of statistical uncertainty. To reduce these correlations, the spectral analysis of 4FGL-DR3 adopted the PLEC4 model ( PLSuperExpCutoff4 in the Fermitools 15 , specifically designed to minimize the covariance between the model parameters at a reference energy, E 0 :  \ndN dE = N 0 ( E E 0 ) -Γ+ d b exp [ d b 2 ( 1 -( E E 0 ) b )] (15)  \n≡ N 0 x -Γ+ d b exp [ d/b 2 ( 1 -x b )] , x = E E 0 . (16)  \nFigure 15. Realizations of the PLEC4 model for three values of the peak curvature, d p = 0 . 2 , 0 . 6 , 1 . 3, and the asymmetry parameter b = 0 . 3 , 0 . 6 , 1 . 3. We set E p = E 0 = 1GeV, which uniquely determines Γ and Γ 0 . The SEDs are normalized to 1 at E p but offset by 50% for clarity. Each curve is labeled with the tuple ( d , b , Γ 0 ). The shaded regions show the difference between the SED and the asymptotic power law, normalized to 0.1 GeV.  \n<!-- image -->  \nFor PLEC4, N 0 is the flux density at E 0 . The PLEC4 spectral shape is exactly the same as for the PLEC model, as can be seen by identifying  \nPLEC ⇐⇒ PLEC 4  \nE c = E 0 ( b 2 d ) 1 b (17)  \nΓ 0 = Γ -d b (18)  \nN 0 = N 0 exp ( d b 2 ) . (19)  \nIn the fits presented below, the reference energies E 0 are confined to a relatively narrow energy range with a typical value of 1.6 GeV.  \nThe PLEC4 model parameters Γ and d are related to the logarithmic derivatives of the model. The gamma-ray convention defines the photon spectral index as  \nΓ E = -∂ log N ∂ log E = Γ -d b ( 1 -x b ) = Γ 0 + d b x b , (20)  \nfrom which it clear that Γ is the photon index at the reference energy ( x = 1) and Γ 0 is the asymptotic spectral index ( x → 0). The spectral curvature can be defined as -∂ 2 log N ∂ log E 2 = dx b , so d is the SED curvature at E 0 .  \nFigure 15 shows some example realizations of the PLEC4 model and in particular highlights the role of d in producing narrow (large d ) or broad (small d ) peaks. b controls the overall shape of the spectrum, with larger values of b producing greater asymmetry. For b ≪ 1 the spectrum becomes more symmetric, with both a slower cutoff and a slower convergence to a power law at low energies. b > 1 produces asymmetric spectra with sharp, 'super-exponential' cutoffs and power-law behavior at low energy. When b < 0, the cutoff and power-law regimes are reflected. For many of these examples, the transition to power-law behavior has not yet happened at 100 MeV. In this example, the local spectral index at 100 MeV agrees with the asymptotic index to within 0 . 1 for b = 0 . 6, but for b = 0 . 3 it differs by 0.3 ( d = 0 . 2), 0.9 ( d = 0 . 6), and 1.9 ( d = 1 . 3).  \nThe energy at which the SED peaks (if the peak exists, is nonzero, and is finite) is  \nx p = E p E 0 = [ 1 + b d (2 -Γ) ] 1 b = [ b d (2 -Γ 0 ) ] 1 b . (21)  \nNote that b < 0 can yield a valid E p if Γ 0 ≥ 2. E p is a key physical property of interest, and we can obtain a second key property, the width of the SED at its peak, by evaluating the curvature at x p :  \ndx b p ≡ d p = d + b (2 -Γ) = b (2 -Γ 0 ) . (22)  \nFor reference, the maximum expected spectral curvature, obtained for monoenergetic synchrotron radiation, is d p = 4 / 3.  \nWe expect these quantities associated with the SED peak to be well-determined and relatively independent because E p is generally close to E 0 . However, Equation 22 indicates that we cannot measure the peak curvature and the asymptotic index independently when b is fixed, and since we expect d p to be the better-determined quantity, we cannot use the PLEC4 model to extrapolate the asymptotic spectral index, Γ 0 . Γ 100 could be used as a proxy for Γ 0 (when E p ≫ 100 MeV and the peak is not too broad), but it too suffers from bias: for a typical E 0 = 1 . 6 GeV and with b = 2 / 3,  \nΓ 100 = Γ 0 +0 . 23 d = Γ -1 . 26 d. (23)  \nThis formulation all but maximizes the correlation between the model parameters Γ and d , undoing the benefits of PLEC4, and like Γ 0 , Γ 100 will potentially be biased if b = 2 / 3 does not describe the true SED.  \nThus, to use the PLEC4 model to estimate reliable physical properties, we need to answer two primary questions:  \n- 1. For bright pulsars where b can be measured, does the PLEC4 model allow robust estimation of other physical parameters of interest, like Γ 100 ? That is, can PLEC4 reproduce the true shapes of pulsar SEDs?\n- 2. For fainter pulsars, for which b must be fixed to a canonical value, does the reduced parameter space introduce a bias on physical parameters?", '6.2. Spectral Model Fits': "We now address these questions using phase-averaged spectral fits. All spectral results in this section rely on the 4FGL-DR3 spectral fitting methodology. Some results are incorporated directly from that catalog, and we have extended these results using the DR3 model as a baseline. In particular, we have  \n- · Performed spectral fits in which b was free to vary for 119 pulsars with TS > 1000 (c.f. 28 pulsars with TS > 10,000 in 4FGL-DR3).\n- · Fit the 28 b -free pulsars of 4FGL with a b = 2 / 3 model for comparison.\n- · Re-fit PSRs J0922+0638, J1757 -2421, J1913+1011, and J2205+6012 with a PLEC4 model, as opposed to the power law used in 4FGL-DR3.\n- · Added a power-law component to the spectral models for pulsars with TeV counterparts (PSRs J0007+7303, J0205+6449, J1016 -5857, J1119 -6127, J1714 -3830, J1833 -1034, J1907+0602, and J2032+4127) to encapsulate potential emission from a pulsar wind nebula and re-optimized the PLEC4 model. Crab, Vela and a half-dozen other pulsars were already colocated with TeV power-law sources in 4FGL.\n- · Directly measured the low-energy spectral indices with a power-law fit of the 50-300 MeV band, for the 42 pulsars with TS > 25 in that band and TS > 1000 overall.  \nOf the 255 pulsars, 251 produce a spectral fit with a well-defined E p .  \nHere we quickly note a possible point of confusion. The reference energy E 0 in Eq. 15 is used in fitting models and evaluating covariance between parameters. From these results a 'pivot energy' can be determined at which the uncertainty on the flux density normalization N 0 is minimized. The results presented here require the covariance matrices from the likelihood fits and so use the reference energy, and thus parameter values will generally differ from those in the FITS files accompanying 4FGL-DR3, which have been scaled to the pivot energy. bf A script provided with the Online Material shows how to correctly unpack the parameters and covariance from the electronic catalog.", '6.2.1. Distribution of b': 'First, we consider the distribution of best-fit b values, plotted against ˙ E in Figure 16. It is clear that b does not take on a constant value, so we used maximum likelihood  \nFigure 16. The best-fit values of b , with 1 σ uncertainties, shown as a function of ˙ E . Maximum likelihood linear models are shown for each class and for the combined population. The centroid indicates the model fit, while the distance between the dashed lines indicates the intrinsic scatter ( ± σ i , see main text). The dark shaded region indicates the same quantities for the full population with measured values of b .  \n<!-- image -->  \nto fit a linear model b = b 0 + m b log 10 ( ˙ E/ 10 36 erg s -1 ) to each class and to the full population. We assumed that each measured value follows a Gaussian distribution with the measurement uncertainty, and we also included an intrinsic scatter, σ i , which is added in quadrature to the measurement error.  \nWe excluded the Crab pulsar from the fits. In Figures 17 through 21 its wellmeasured parameters generally deviate from those of most other pulsars. With its high ˙ E , including it changes the fit results substantially, such that they no longer describe most other gamma-ray pulsars well. If the Crab is atypical only in having both high ˙ E and good statistics, then perhaps our simple linear fit breaks down for high ˙ E . But the Crab is also atypical in being a rare pulsar with much stronger flux density in the MeV range than in the GeV range, as we further discuss in Section 8.1.  \nThe results, depicted in Figure 16, are (1) there is strong evidence (log likelihood improves by > 40) of intrinsic scatter, i.e. σ i > 0, and (2) there is no evidence (log likelihood improves only by a few) for a different linear relation amongst the classes. Therefore, we report a single model for all pulsars, finding b 0 = 0 . 42, m B = -0 . 15, and σ i = 0 . 15.  \nThese values are intended as a guide only: the simple model does not fully capture the observed variations, and the fitting approximates the measurement uncertainty on b as a Gaussian rather than using the full likelihood surface. Nonetheless, it is clear that many pulsars fit with a b = 2 / 3 model will suffer appreciable biases. In particular, estimates of Γ 0 and Γ 100 for MSPs (with ˙ E ∼ 10 34 erg s -1 ) will tend to be too soft, and those of RL pulsars (with ˙ E > 10 36 erg s -1 ) too hard. We analyze these biases more quantitatively below. With reference to Figure 15, we also see that the distributions of b indicate that MSP spectra hew more closely to the power-law with exponential cutoff model, while young pulsars, particularly high-˙ E RL pulsars, have substantial curvature over the full LAT band and more gradual high-energy cutoffs. This is suggestive of MSPs possessing an intrinsically narrower radiating particle distribution than do their RL and RQ counterparts; how this might connect to the smaller magnetospheres of MSPs remains to be understood.  \nIn Appendix C, we study the change in model parameters between the 119 pulsars which have b -free and b = 2 / 3 fits available. In particular, we find that d and Γ agree well between the models, while Γ 100 is irrecoverably biased. The physical properties d p and E p are robustly measured in both b -free and b = 2 / 3 fits and we conclude that we can make inferences on their distribution using the entire population of pulsars. On the other hand, the low-energy spectrum is unreliable, and values inferred from the b = 2 / 3 model could be very far from the true values.  \nFigure 17 establishes the converse, that the b -free PLEC4 model yields a good estimate of Γ 100 , at least for the small sample for which model comparison is possible. The figure shows the indices resulting from the 119 low-energy fits, of which 42 have TS > 25 in the 50-300 MeV band. The geometric mean energy of this band is 122 MeV, but the improvement of LAT sensitivity with energy, along with the use of likelihood weights, means the effective energy is greater. Nevertheless, the agreement is good for all 3 pulsar classes 16 , and consequently we conclude that the b -free Γ 100 measurements are reliable. If additionally E p ≫ 100 MeV and d p is not too small, Γ 100 may approximate Γ 0 well. We discuss this further below.', '6.3. Physical Properties': "We can thus make inferences about three physical properties: d p and E p for the full population of 251 pulsars, and Γ 100 for the 119 b -free pulsars. Unless otherwise noted, results from b -free fits are preferred and presented whenever available.  \nFirst, we consider the distribution of spectral curvature, d p , or equivalently, (inverse) SED peak widths, with results shown in Figure 18. d p depends strongly on ˙ E , varying in width by a factor of about 3 over the main population. The narrowest peaks-primarily of MSPs, but with members of each class-saturate at about 4 / 3, corresponding to a monoenergetic synchrotron spectrum. Only PSR J1827 -1446, a  \nFigure 17. The comparison between spectral indices measured directly with a low-energy fit (y-axis) and the predicted local index at 100 MeV from b -free models.  \n<!-- image -->  \nFigure 18. Dependence of the spectral peak curvature d p on ˙ E , showing a well-defined trend. The horizontal dashed line at 4/3 indicates the expected peak curvature for a monoenergetic curvature radiation spectrum. The lefthand panel shows results for all 251 pulsars (i.e. mixes b -free and b = 2 / 3 fits), while the righthand panel shows only the 116 pulsars with a b -free fit (and a well-defined peak energy). The solid line indicates a best-fit relation to the merged sample (lefthand panel).  \n<!-- image -->  \nFigure 19. Dependence of the spectral peak energy on ˙ E, with low-˙ E pulsars having a peak energy roughly twice that of high-˙ E pulsars. The samples and best-fit model are as in the left panel of Figure 18.  \n<!-- image -->  \nradio-quiet pulsar, lies out of family with d p = 2 . 6 ± 0 . 5 and one of the largest E p values (see below). To quantify the relation, we fit a simple power-law model to the full population (left panel, excluding Crab), finding d p = 0 . 46( ˙ E/ 10 36 erg s -1 ) -0 . 12 with an additional scatter of 17% added in quadrature.  \nNext, we consider the distribution of peak energies, E p , shown in Figure 19. As with d p , there is a trend with ˙ E , which we characterized as E p = 1 . 1 GeV( ˙ E/ 10 36 erg s -1 ) -0 . 09 with an additional scatter of 30%. This ˙ E evolution leads to a distinct difference between the pulsar classes, with median values of 1.7 GeV (MSP), 1.2 GeV (RQ), and 0.9 GeV (RL). Other population variables, such as the magnetic field at the surface or light cylinder, do not reduce the larger intrinsic scatter.  \nd p and E p are also correlated (Figure 20), but this relation appears to be entirely captured by ˙ E evolution. Using the single-variable models determined above removes the correlation to the level of the intrinsic scatter in E p . In general, we arrive at a picture of larger peak widths at lower peak energies, a paradigm clearly illustrated by Crab-like pulsars with very broad peaks with, presumably, substantial contributions from secondary particles. The large widths are such that some pulsars, like the Crab  \nFigure 20. Dependence of the spectral curvature d p on the peak energy, E p . Panel selections are as in Figure 18.  \n<!-- image -->  \nand Vela, emit into the TeV range (HESS Collab. et al. 2023), implying very high particle energies, which emission models need to accommodate.  \nFinally, we consider Γ 100 , a proxy for the asymptotic spectral index. As an illustration of the possibility of making an incorrect inference about the physical properties of the pulsar population, in the first panel of Figure 21 we show the results of 251 b = 2 / 3 fits, revealing an apparently significant trend with hardening spectral indices for increasing ˙ E . However, this result is largely due to the bias relating Γ 100 to the spectral curvature and thus reproduces the trend shown in Figure 18. In the right panel, we show the b -free fits, which we demonstrated provide relatively unbiased estimators for Γ 100 . The trend is reduced or absent: the bulk of the pulsars with ˙ E > 10 35 erg s -1 have a spectral index near 1 . 2, and somewhat lower ( ≈ 1) for lower ˙ E . Below this (possible) threshold, harder spectra become apparent, with values approaching the monoenergetic limit of 2/3. These pulsars have the highest values of E p , and so for them Γ 100 may be closer to Γ 0 .  \nBoth observed SED shapes and pulse profile shapes should depend on the magnetic field strength and configuration, as well as on 'geometry', meaning for example the inclination α of the magnetic axis relative to the rotation axis; the inclination ζ of the rotation axis relative to the line of sight; and the radius of the light cylinder R LC = cP/ 2 π . We thus searched for correlations between measured profile and spectral parameters. Figure 22 suggests that pulsars with classic Vela- or Crab-like gamma-ray pulse profiles, with ∆ ≲ 0 . 5 rotations, have typical SED peak energies. The highest E p values seem to occur mainly for pulsars where the second gammaray peak follows the first more closely. This is perhaps related to the evidence that the gamma pulse farther from the magnetic pole generally has the harder spectrum (Figure 13).  \nFigure 21. Dependence of the low-energy spectral index Γ 100 on ˙ E . The left panel shows fits in which every pulsar has b = 2 / 3, revealing a strong but erroneous trend. The right panel shows only the 119 unbiased b -free fits. The horizontal dashed line indicates the monoenergetic limit, 2/3.  \n<!-- image -->  \n<!-- image -->  \nFigure 22. Left: The SED peak energy E p versus the phase separation ∆ between the two principal gamma-ray peaks. The Crab and PSR J1730-3350 are off-scale with E p ≈ 50 MeV. Right: E p distributions for pulsars with closely- (top) and widely- (bottom) spaced peaks. Pulsars with closely-spaced peaks perhaps have the highest E p values.  \n<!-- image -->  \nWe summarize our main results here:  \n- · The broadband shape, governed by the parameter b , evolves with ˙ E , yielding slowly-decaying spectra for high-˙ E pulsars and sharply-cutoff spectra for low-˙ E pulsars.\n- · Pulsar spectral energy densities (SEDs; E 2 dN/dE ) are most effectively characterized by the curvature at the SED peak, d p , which is inversely proportional to the SED peak width. d p anti-correlates with ˙ E with low (17%) scatter: see Figure 18.\n- · The energy at which the SED peaks, E p , also varies inversely with ˙ E with somewhat larger (30%) scatter: see Figure 19.  \n/negationslash  \n- · Consequently, the spectra of MSPs are more sharply peaked than those of young pulsars, and reach the highest peak energies.\n- · The asymptotic ( E → 0) spectral index, Γ 0 lies out of the LAT band. The spectral index near 100 MeV, Γ 100 , is close to Γ 0 when E p ≫ 100 MeV, and it can be estimated for 116 of the brightest pulsars (Figure 21). These indices appear to take on an almost universal value of about 1.2, though low-˙ E pulsars, particularly MSPs, may have harder spectra that approach the 'monoenergetic limit' of 2/3.\n- · Thus, for our LAT pulsars, nearly all of the rotational energy that is converted to electromagnetic radiation is concentrated in a single decade of the spectrum.\n- · The SED peak energy E p seems to be highest for pulsars with pulse phase separations ∆ < 0 . 4 (Figure 22). This, coupled with the observation that the second pulse tends to have a harder spectrum than the first (Figure 13), is new empirical evidence that the same underlying magnetosphere configurations affect both profile and spectral shapes.", '6.5. Notes on Individual Pulsars': 'Twenty gamma-ray pulsars in Tables 13 and 14 are flagged with † symbols, meaning that no 4FGL-DR3 source with Test Statistic > 25 is colocated with the pulsar. The pulsed signature is detectable in these cases, in spite of a low gamma-ray flux or large local background.  \nGuillemot et al. (2013) provide spectral parameters for the double pulsar PSR J0737 -3039A via a dedicated analysis. PSR J1513 -5908 (B1509 -58) is brightest in the MeV domain, with a spectral cutoff below 100 MeV, measured by Kuiper et al. (2018). The colocated extended source 4FGL J1514.2 -5909e corresponds to the pulsar wind nebula MSH 15 -52 and its spectral parameters do not apply to the pulsar. Similarly, PSR J1023+0038 is a transitional pulsar, and the 4FGL-DR3 spectral results are dominated by the data obtained after it transitioned to its accreting state, where the (unpulsed) gamma-ray flux is greatly enhanced.  \nPSR J1748-2815 was co-located with a TS ∼ 300 source in the LAT 8-year catalog. Prior to 4FGL-DR1 the LAT team adopted a new model for the diffuse Galactic gamma-ray emission that is brighter around the Galactic center, and the source disappeared, although its pulsations remain significant with the model weights used here.  \nFor PSR J1531 -5610, 2PC mistakenly reported a luminosity L γ over ten times smaller than the value in Table 13. An outlier in 2PC, with corrected L γ its efficiency ϵ = L γ / ˙ E becomes quite typical.  \nTable 13 . Spectral fitting results for young LAT-detected pulsars.  \nTable 13 continued on next page', 'Table 13 (continued)': 'Table 13 continued on next page  \nTable 13 continued on next page  \nTable 13 (continued)  \nTable 14 . Spectral fitting results for LAT-detected millisecond pulsars.', 'Table 14 (continued)': 'Table 14 continued on next page  \nTable 14 continued on next page', '6.6. Searches for Pulsar Flux Variability': "In 4FGL-DR3, fluxes for sources with Variability Index < 24 . 7 are considered stable. Three pulsars exceed this threshold: PSR J2021+4026, cited in Section 3 with Variability Index = 285, and the two transitional MSPs J1023+0038 and J1227 -4853 with Variability Index = 58 and 1470, respectively, which exhibit large flux changes when switching between accretion- and rotation-powered states, but appear to be stable emitters like all other gamma-ray MSPs when in the rotationpowered state. In contrast, ˙ E of PSR J0540 -6910 in the Large Magellanic Cloud increased by a third, with no corresponding change in the gamma-ray flux (Marshall et al. 2015), contrary to what the correlation between L γ and ˙ E seen in Figure 23 suggests ( Variability Index = 11 . 1). A search for variability in our pulsar sample yielded no new cases.  \nThe linear increase in pulsed significance (weighted Htest) with time of PSR J2256 -1024 changes abruptly near MJD 56600. It has Variability Index = 18 . 3. A dedicated analysis confirms that the pulsar flux is constant, but that nearby 4FGL J2250.0+3825, associated with a blazar, has a very bright flare peaking at the same time. The softer significance slope is thus due to increased background rather than decreased signal. Exposure changes due to pointing changes also cause Htest vs MJD slope changes.  \nPSR J2116+1345 is an MSP discovered using the Arecibo radio telescope in a search of a pulsar-like LAT source. Gamma-ray pulsations were found using Arecibo and GMRT timing after we froze the list of pulsars reported in this catalog. Its corresponding source, 4FGL J2117.0+1344, has Variability Index = 16 . 2. Here too, pulsed significance is not linear with time. Analysis showed that the decreased rate  \nof significance growth between MJD 56900 and 57900 is due to the flaring activity of 4FGL J2035.4+1056 (PKS 2032+107) which is 10.5 ° away from the pulsar.  \nClark et al. (2023b) discovered dips in the gamma-ray flux of 'spider' MSPs at orbital phases when the companion star lies between Earth and the pulsar. The 4FGLDR3 fluxes are integrated in both rotational and orbital phase, these eclipses hence do not affect Variability Index . Orbital modulations (see for example An et al. 2020, or the references cited in Table 15) similarly do not influence Variability Index .", '6.7. Luminosity': "Most of the low-frequency Poynting flux generated by the dipole radiation of the rotating neutron star converts to a wind of electrons and ions towards and/or beyond the light cylinder, carrying away angular momentum and energy, braking the rotation. About 10% of the known pulsar population-those in this catalog-convert much of their spindown power into intense gamma-ray beams. Pulsar energetics contribute a piece to the broader puzzle of Galactic ecology, as recently illustrated by the PeV energy detection by LHAASO of the Crab nebula (Cao et al. 2021) and by the Fermi observation of GeV-band halos around 3 bright HAWC sources (Di Mauro et al. 2021), and gamma-ray pulsars allow the most direct look at the acceleration mechanisms.  \nIntegrating the LAT spectral energy flux E dN dE above 100 MeV yields G 100 , which is tabulated in 4FGL. The range of G 100 values for the gamma-ray pulsars is apparent in Figure 6. We can then calculate the luminosity  \nL γ = 4 πd 2 f Ω G 100 , (24)  \nusing the distance d to the pulsar, and the beaming fraction f Ω .  \nFigure 23 shows that the luminosity (Eq. 24) is highly correlated with the available braking power. However, at a given ˙ E there are two decades of dispersion in the observed efficiency η = L γ / ˙ E , shown in Figure 24. Distance errors cause some of the dispersion. The total radiated power likely depends on the inclination α between the neutron star's magnetic and rotation axes, adding more pulsar-to-pulsar variation. The factor f Ω in Eq. (24) is the ratio of total beam power averaged over the whole sky to the power in the beam slice illuminating the Earth averaged in phase. It is defined in 2PC Eq. 16, and by Romani & Watters (2010), and depends on both α and on the inclination ζ between the rotation axis and the line of sight. Different models predict different pulse profile shapes and thus f Ω , and in any case observational difficulties cause large uncertainties in the α, ζ model inputs. Lacking f Ω ( α, ζ ), our choice of f Ω = 1 likely contributes to the observed η dispersion. A key test of a successful theory of pulsar radiation is its predictions of L γ for given pulsars. Kalapotharakos et al. (2023) generally find f Ω < 1 for the LAT pulsar sample, which would make our L γ values over-estimates.  \nThe upper diagonal line in Figure 23 shows 100% efficiency, L γ = ˙ E , for converting spindown power into gamma rays. Progress since 2PC in determining pulsar  \ndistances and MSP proper motions has eliminated pulsars more than 1 σ above the line. Furthermore, R NS = 10 km is an outdated convention, and ˙ E ∝ M NS R 2 NS is ≈ 2 times larger than what we use, reducing all η values. The neutron star mass values M NS in our sample probably have variance < 20% (see for example Figure 5 in Shamohammadi et al. 2023), contributing imperceptibly to the luminosity dispersion.  \nThe lower diagonal shows the 'heuristic' gamma-ray pulsar luminosity,  \nL h γ = ϵ γ ˙ E ≈ ˙ E √ 10 33 erg s -1 ˙ E = 10 33 erg s -1 √ ˙ E 10 33 erg s -1 , (25)  \nmotivated by the idea that above some minimum open field-line voltage V ≃ √ 10 -5 ˙ E volts, electron-positron cascades occur in an acceleration region with a self-governing size that produces a gamma-ray efficiency ϵ γ ∝ 1 /V (Arons 1996). Recent particle-incell simulations of pulsar magnetospheres also support the idea of a critical pair production rate below which pulsar magnetospheres become 'dead' (e.g. Kalapotharakos et al. 2018; Philippov & Spitkovsky 2018; Chen et al. 2020). The scaling of L γ with ˙ E depends on whether gamma rays are radiated primarily via synchrotron or curvature processes (Kalapotharakos et al. 2019). We discuss our results in the context of these simulations further in Section 8.", '7. THE PULSARS NOT SEEN': "Here we address known pulsars not seen with the LAT, in spite of characteristics shared with gamma-ray emitters. PSR J1846 -0258 (Section 7.2) is a unique case of a pulsar seen in LAT data only below 100 MeV (Kuiper et al. 2018), not using the methods of this catalog. Finally, we present the LAT's sensitivity for future pulsar discoveries.", '7.1. Radio and X-ray pulsars': "Figure 5 shows that about 5% and 35% of young pulsars and MSPs, respectively, are seen with the LAT in the ˙ E ≈ 10 33 erg s -1 decade, climbing to over 60% and 80% at high ˙ E . Removing pulsars discovered using LAT sources as seeds, the low ˙ E fractions decrease to about 3% and 25%. 2PC Table 13 listed 28 pulsars with ˙ E > 10 36 erg s -1 that were not seen at the time, of the 64 then known. Since 2PC, half of these 28 have been detected and are in this catalog, while the total number increased to 71. The majority of those not yet seen are young X-ray pulsars, discussed below.  \nWith a distance estimate, we calculate the pulsar's 'heuristic' flux G h = L h γ / 4 πd 2 , which converts to the units of the spectral measurements using G 100 /G h = 3 . 33 × 10 -27 = 10 -26 . 48 . Figure 5 shows the fraction of detected pulsars rising steadily with G h above the LAT detection threshold near G 100 ≈ 10 -12 erg cm -2 s -1 (detailed in Section 7.3). Here, removing pulsars discovered in searches of unidentified LAT sources moves the MSP points closer to those of the young pulsars. Smith et al. (2019) discuss whether the non-detection of lower ˙ E pulsars is due to a deathline,  \nFigure 23. Gamma-ray luminosity L γ = 4 πf Ω d 2 G 100 in the 0.1 to 100 GeV energy band versus spindown power ˙ E . The vertical error bars from the statistical uncertainty on the energy flux G 100 are colored in the electronic journal version. The vertical error bars due to the distance uncertainties are black-dashed, and generally larger. Doppler corrections (Section 4.3) have been applied to MSPs with known proper motions, leading to visible horizontal error bars in many cases, and when the corrections yield negative ˙ E , we place the limit at ˙ E > 2 × 10 32 erg s -1 for clarity. Luminosity upper limits use distance upper limits, while ˙ E limits result from insignificant proper motion measurements or undetermined distances. We use a light color shade to plot these limits in order to de-emphasize them, and we do not include them in the tallies in the legend. The upper diagonal line indicates 100% conversion of spindown power into gamma-ray flux: pulsars above this line, may have smaller distance d , and/or the assumed beam correction f Ω ≡ 1 is wrong. The lower diagonal line indicates the heuristic luminosity L h γ ∝ √ ˙ E , to guide the eye. For the Crab and PSR J0540 -6919 in the Large Magellanic Cloud, at far right, the upper points include the X-ray energy flux.  \n<!-- image -->  \n˙  \nmeaning that the emission mechanisms cease below a minimum braking power, or comes from the small number of pulsars close to Earth and the low detection fraction for low G h .  \nFor the non-detections at high ˙ E , for the young pulsars, Johnston et al. (2020) find that beam geometry and pulsar inclination dominate (the beam misses the Earth). Rookyard et al. (2017) observe that broad radio pulses are a good proxy for the inclinations unfavorable to gamma detection. An example is the young pulsar PSR J1302-6350 (B1259-63), with large √ ˙ E/d 2 and co-located with a LAT source (see  \nFigure 24. Gamma-ray efficiency η = L γ / ˙ E versus spindown power ˙ E . The error bars are as in Figure 23. The dashed diagonal shows η ∝ ˙ E -1 / 2 , corresponding to the heuristic luminosity L h γ , to guide the eye. The correlation is more striking than in 2PC.  \n<!-- image -->  \nTable 5), yet showing no gamma-ray pulsations (Abdo et al. 2011; Smith et al. 2019). GeV flares occur during periastron every 3.4 years, when the neutron star moves out of the stellar disk of the companion star (Chang et al. 2018). Here, the absence of gamma-ray pulsations is likely due to the beam orientation. Some gamma-ray pulsars are detected at much lower efficiency than expected from the general trend (e.g. PSRs J0659+1414, J0940 -5428 and J1740+1000 near ˙ E = 4 , 23 , and 190 × 10 34 erg s -1 , respectively, in Figure 24). If these objects are typical gamma-ray pulsars, then Earth lies well outside of the standard strong beam, and a weak fringe of this beam or an alternate emission component may be responsible for such 'sub-luminous' gamma-ray pulsars (Romani et al. 2011).  \nDifferent factors contribute to the larger fraction of MSP gamma-ray detections. For gamma-rays emitted in the inner magnetosphere, the small light cylinders broaden the beams. However, current sheet emission models predict similar beams for MSPs and young pulsars. MSPs have lower ˙ E , and thus lower L γ than young pulsars (see Figure 23) but they are also generally closer to Earth. Figure 19 shows that E p is slightly higher for MSPs, facilitating detection, as was highlighted by Kalapotharakos et al. (2017) using 2PC data. The details of which MSPs are not seen by the LAT remain to be clarified.  \nAnother category of as-yet unseen pulsars is those for which the radio rotation ephemerides do not suffice for multi-year gamma folding. An example is the 33 radio MSPs found in deep radio searches targeting unidentified LAT sources that have not yet been confirmed as gamma-ray pulsars (Tables 1 and 6). It also occurs for pulsars found in a survey, but never included in a long-term timing campaign: psrcat lists 278 pulsars with 30 < P < 800 ms but without ˙ P . Nearly half were recently discovered by the GPPS 17 and CRAFT 18 surveys, both on the FAST telescope (Han et al. 2021), with improved measurements ongoing. Some of the others have been added to the timing program for Fermi at the Parkes radio telescope, when their declination and radio brightness allow. Of those found to have large ˙ E , a fraction may be detected in LAT data in the coming years.  \nX-ray timing with RXTE enabled the gamma-ray discovery of PSR J0540 -6919 in the Large Magellanic Cloud (LMC) (Ackermann et al. 2015). Subsequent Swift timing revealed a rare 30% increase in ˙ E and a near-zero braking index (Marshall et al. 2016). XMM-Newton timing assisted the LAT pulsation searches for PSR J2022+3842 (Limyansky 2022).  \nSeveral of the high-˙ E pulsars unseen in LAT data mentioned above are radio-quiet X-ray pulsars. Timing noise for such young pulsars is generally high, and while Xray observations made over a few years with XMM or Chandra yield astrometry and the spindown rate, they cannot give the phase-connected ephemeris necessary for multi-year gamma-folding. Thus, for many years, only pulsars which have low timing noise and/or are bright enough for routine observations with RXTE or Swift could be gamma-ray phase-folded with the same level of confidence as we have for our thousand radio ephemerides.  \nNICER has provided new capabilities, allowing precision X-ray timing for more pulsars than previously possible. An example is the study of PSR J1813 -1749, the third pulsar in 2PC Table 13, by Ho et al. (2020). The pulsar has no 4FGL counterpart and gamma folding yields a null result. Ho et al. (2022) present NICER timing for 5 other high-˙ E X-ray pulsars, unseen with the LAT.  \nX-ray and optical studies of unidentified 4FGL sources have led to the pulsar candidates listed in Table 15. Many of these have led to MSP discoveries in radio pulsation searches using the methods described in Section 3.3 or gamma-ray MSP discoveries in restricted blind searches (Section 3.4), once the orbital parameters are better constrained by further multi-wavelength studies. Transitional MSPs, of which PSR J1023+0038 is the archetype (see Stappers et al. 2014, and Appendix A) , may 'turn on' and show pulsations at some future date. In a similar vein, Hare et al. (2019) identified a point source in the heart of an X-ray pulsar wind nebula that merits pulsation searches.  \nTable 15. Candidate 'spider' and transitional MSPs (tMSPs) discovered in deep searches of previously unidentified LAT catalog sources. BW means black widow, RB means redback, and DStMSP means that a 'disk state' is presumed to have obscured pulsed detections to date. 'PL' in place of the SED peak energy (in GeV) indicates that the best spectral fit uses a power law instead of log parabola. Two tMSPs seen to pulse in gamma rays are PSRs J1023+0038 and J1227 -4853. All LAT names are 4FGL except for 2FGL J0846.0+2820, detected only in 2009.", '7.2. PSR J1846-0258': "PSR J1846-0258 is a young ( τ ≈ 723 yr), high magnetic field (5 × 10 13 G) pulsar located in the Kes 75 SNR (Gotthelf et al. 2000). Although radio quiet (Archibald et al. 2008), it has been timed regularly in X-rays by RXTE, Swift , INTEGRAL, and NICER. PSR J1846-0258 behaves largely as a rotation-powered pulsar, but has twice exhibited magnetar-like outbursts (Gavriil et al. 2008; Blumer et al. 2021). These outbursts have garnered interest in this source as a 'magnetar-pulsar transitional object', possibly able to shed light on the evolutionary relationship between rotation-powered pulsars and magnetars. Additionally, it is one of 18 'soft' gammaray pulsars (Kuiper & Hermsen 2015), with strong emission in the MeV range. The X-ray spectrum suggested that the LAT would be able to detect the pulsed highenergy tail, aiding research into the soft pulsar emission mechanism (Torres 2018).  \nPSR J1846-0258 is the only pulsar with a published gamma-ray detection (pulsed 4 . 2 σ , via unweighted Z 2 m =1 , see Eq. 2 in Sec. 3) that is not detected using the methods of this catalog (Kuiper et al. 2018). It has no 4FGL-DR3 counterpart source, and thus no L γ measurement. Following Kuiper et al. (2018)'s prescription, we confirmed this detection. It increases to 4 . 3 σ using Pass 8 data. As discussed in Section 3, timing solutions may produce significance vs time curves which peak outside the ephemeris validity range. This is indeed the case with PSR J1846-0258, with a maximum HTest ( Z 2 m =1 ) of 4 . 4 σ (5 . 0 σ ) occurring ∼ 484 days after the end of the published 3295 day timing solution (including a 153 day loss of coherence). We also used NICER to construct an independent ∼ 770 day timing solution outside the time range covered  \nby Kuiper et al. (2018), which we applied to similarly filtered LAT data. This did not yield a LAT detection, which given the 4 × shorter span of the NICER data is consistent with expectations from the previous detection (Limyansky 2022).", '7.3. Flux Upper Limits and Sensitivity': 'Discovery of most of the pulsars in this list started with an initial phase-integrated source detection of the type applied to all Fermi -LAT catalog entries. We summarize the 4FGL procedure: for an all-sky model defined by pointlike (Kerr 2010) which includes all current sources, a new trial source is introduced at the position of each of HEALPix nside=256 ( ∼ 1/4 deg) pixels. Its flux is optimized using each of five fixed trial spectral shapes. Three are power laws, but two are peaked, at 1 and 2 GeV. The one with the highest TS, if greater than 10, is used as a seed to add to the model. The model is re-optimized with now free spectral parameters, including the new source or sources. All sources are also relocalized. All that survive with TS > 10 are used for additional iterations. Finally the list is used as seeds for the full catalog gtlike analysis. After this stage only sources with TS > 25 are included in the final catalog. Similarly, 2PC started from the 3FGL catalog, which however used only a single power-law trial spectrum.  \nConstruction of an estimated sensitivity threshold map, shown in Figure 25, followed a similar procedure. For each HEALPix pixel location a trial source was introduced, now with a pulsar-like spectrum, and the flux likelihood function determined. Instead of optimizing it as before, the value corresponding to its 95% cumulative probability was recorded for that point. A more precise procedure would have instead inserted simulated trial data with varying intensity, applying the full subsequent analysis procedure, to determine the threshold curve.  \nWe assume that our procedure yields the correct position dependence. Converting the sensitivity estimates to a 50% completeness threshold requires an ad hoc scale adjustment, which we obtain by comparing the distribution of sensitivity values to the detected pulsar fluxes, shown in Figure 26 as a function of latitude. This detection threshold analysis considers only phase-integrated source detections.  \nAs discussed in 2PC and by Hou et al. (2014), and further demonstrated by the 15 pulsars in this catalog that were detected via their pulsations but have TS below the 4FGL threshold, narrow gamma-ray pulses can facilitate source discovery. Such pulsars appear in Figure 26 near the shaded area.  \nFinally, a requirement for inclusion in the catalog is that no new candidate source be within 0 · . 5 of an existing one. Therefore, in the map in Figure 25, we set such pixels around each 4FGL source to nan . The all-sky sensitivity threshold map is provided in FITS format with the online material.  \nFigure 25. Hammer-Aitoff projection of the estimated LAT 12-year sky-survey energy flux pulsar detection sensitivity.  \n<!-- image -->  \nFigure 26. Integral energy flux from 0.1 to 100 GeV, G 100 , versus Galactic latitude b (scaled as b 0 . 65 for clarity). The shaded area shows the 5% to 95% percentile range of the 12-year pulsar detection sensitivity in the sky directions at the corresponding latitude. The points correspond to the gamma-ray pulsars in this catalog.  \n<!-- image -->  \nFigure 27. One realization of the cutoff energy, ϵ 37 , hi , takes on a nearly constant value, shown here in the estimates for the b -free sample.  \n<!-- image -->', '8.1. Phase-averaged Spectral Trends': "In Section 6, we established that the energy where the SED culminates, E p , and the SED's curvature at that energy, d p , are robust quantities that can be reliably estimated when using the general b -free SED model or the restricted b = 2 / 3 case. E p , d p , and b all evolve with ˙ E , but the scatter of b and E p is much larger compared to d p . For 90% of the pulsar sample, 0.5 GeV <E p < 2.5 GeV, with young, high-˙ E pulsars at the bottom of the range and older, unrecycled pulsars and MSPs at the top.  \nA consequence of the opposed ˙ E evolution of E p and d p is a nearly-universal cutoff energy, ϵ 37 , hi , at which the SED has fallen 1 /e ≈ 0 . 37 from the peak (Figure 27). Another way of looking at this universality is to consider individual spectral shapes. Figure 28 shows the best-fit SED shapes for the b -free pulsars, encoded according to ˙ E . The region between ϵ 37 , hi and ϵ 37 , lo (the mirror 1/ e point below E p ) is highlighted, allowing a visual estimate of the effective maximum energy and width. The upper peak bounds cluster in the 3-10 GeV range, independently of ˙ E . Thus, we have the simple descriptive result: almost all γ -ray pulsars emit their power within a two-decade (0.1-10 GeV) envelope, with the Crab-like, flat-spectrum, high-˙ E pulsars nearly filling the envelope, while sharply-peaked MSPs concentrate their power into  \nFigure 28. The spectral energy density, S ( E ), for 116 pulsars with an estimate of b . All pulsars have a valid E p , marked with a circle ( b > 0) or a cross (six pulsars with b < 0). The spectra are normalized such that S ( E p ) = ˙ E . Curve colors in ˙ E decades guide the eye. Curves are dashed below ϵ 37 , lo and end at ϵ 37 , lo .  \n<!-- image -->  \nonly the higher decade. In contrast with this absolute energy scale, the spectral curvature as estimated by d p varies strongly with ˙ E . Selection bias, if any, contributes weakly to this result, as we discuss below.  \nIn Section 6, we also showed that the low-energy spectral index Γ 100 could not be measured reliably for b = 2 / 3 models, and of the population of 116 pulsars for which it can be measured, the values tend to cluster around 1.2 with a scatter of a few tenths and only a very modest evolution with ˙ E . This result is clearly visible in Figure 28, in which it can further be seen that many of these SEDs have not transitioned to a power law at 100 MeV, making these measurements effectively upper limits on the asymptotic spectral index. The low-energy spectral shape is thus another frequent property of gamma-ray pulsars, taking on values which are well in excess of the monoenergetic synchrotron radiation spectral index (Γ 0 = Γ sr = 2 / 3) for all but few of the lowest-˙ E MSPs.  \nAs stated, the curvature evolves strongly with ˙ E , and depends on the spectral steepness beyond E p and ϵ 37 , hi . Indeed, Lyutikov et al. (2012) and Richards & Lyutikov (2018) argue that the shape of the high-energy spectral tail is particularly rich in information about the highest energies reached by the particle acceleration mechanisms,  \nand thus, whether and when the radiation-reaction limit is reached. The tail's shape can reveal, for example, the presence of an inverse Compton 'rebound', and also informs about the best candidates for > 10 GeV emission to be observed by Cherenkov telescopes. Barnard et al. (2022) calculate the very high energy Vela spectrum radiated by primary particles in the pulsar magnetosphere and current sheet. Measuring the high-energy spectral shape with Fermi alone is, however, challenging: the photon rates decrease rapidly, and possible contributions from a surrounding PWN and/or SNR must be handled carefully. We envisage future studies dedicated specifically to these measurements.  \nThe low- and high-energy spectral shapes, and thus the curvature, depend on the underlying acceleration of leptons; on variations in the local conditions that govern the balance between acceleration and radiation, observational convolutions of which depend on pulse phase; and on the extent and re-acceleration of pair cascades. Disentangling these effects may not be possible with only GeV data, but the robust evolution of spectral shape demonstrated here is an important input in constraining models of pulsar emission. Future successful theories should also reproduce the approximately constant spectral indices at 100 MeV, the weak decline of E p with pulsar ˙ E , curvatures which track ˙ E , and the overall shape which evolves from flat, Crab-like spectra to MSP-like spectra with well-defined power-law indices and sharp peaks. These measured trends also provide new ingredients for population syntheses: by adopting a particular spectral shape for a given synthesized ˙ E , more reliable estimates of the LAT sensitivity and selection effects can be obtained.  \nThe physical implication of these results may be somewhat diluted since the spectra are phase-averaged, with multiple distinct acceleration regions (whose observational sampling depends on the pulse phase) contributing to and broadening the measured spectral features. On the other hand, we can point to the relatively tight clustering of d p about the ˙ E scaling relation, which suggests that these accelerators must have broadly similar properties. Indeed, while we established a general trend of harder trailing peaks in Section 5.3, these differences do not seem to be large enough to add appreciable scatter to the observed d p values. Furthermore, at least some MSPs approach the shape of monoenergetic curvature radiation, further limiting the possible variations in spectral shape with phase. On the other hand, the preference for b < 1 for high-˙ E pulsars may indicate a richer dependence on phase.  \nKalapotharakos et al. (2019) examine the mechanisms and configurations at play in the current sheet just outside the light cylinder, and conclude that curvature radiation dominates over synchrotron radiation yielding a 'Fundamental Plane' (FP) such that L γ ∝ ϵ 4 / 3 c ut B 1 / 6 S ˙ E 5 / 12 . The spectral cutoff energy ϵ c ut tracks the pulsar's highest photon energies, which reflect the highest attained particle energies. Kalapotharakos et al. (2022) developed a proxy for ϵ c ut called ϵ c 1 , obtained by imposing a PLEC shape with b = 1 on the 4FGL spectra and then finding the energy ϵ c 1 > E p where the SED falls to 10% of its maximum. This is similar to our ϵ 37 , h i shown in Figure 27, and  \nwe include PLEC E10 b23 and PLEC E10 bfr in the 3PC catalog FITS file, with the important difference that we use the b = 2 / 3 and b -free PLEC4 fits that give the best LAT photon data likelihoods. Kalapotharakos et al. (2022, 2023) thus improved the FP data fits as compared to their earlier work, finding good agreement with their FP exponent predictions. Further improvement could include selecting pulsars with the most reliable distances, and reviewing the choice of the ϵ c ut estimator in light of what we established here regarding measurement bias, and improving the uncertainties used in the FP fits. Applying f Ω corrections to L γ for pulsars with reliable estimates of the inclination angles α, ζ ; including Shklovskii corrections for the MSPs ; and perhaps correcting ˙ E for pulsars with neutron star mass measurements may also improve FP tests.", '8.2. Selection Bias': "Are our results representative of all gamma-ray pulsars? Here we address whether observation bias might skew our E p distribution, and prospects for discovering pulsars with lower E p . In 4FGL-DR3, over 2000 sources have E p < 500 MeV. They were fit with the LogParabola function rather than PLEC4, but that affects the E p values only marginally. The LAT detection and analysis chain is clearly sensitive to SEDs with maxima below what we see in our sample.  \nHalf of our gamma-ray pulsars were discovered using radio rotation ephemerides of pulsars discovered independently of the LAT, code 'R' 19 , generally in radio surveys covering large sky areas, with selection biases a priori unrelated to gamma-ray properties. Conceivably, neutron star orientations such that both radio and gamma-ray beams sweep the Earth could correlate with gamma-ray spectral shape. The 'X' pulsars in our sample mitigate such a putative radio-selection bias somewhat.  \nOn the other hand, the average E p of 'U' pulsars is higher than for 'R' pulsars, with 'P' MSPs having higher E p than the 'G' pulsars. MSPs generally have lower ˙ E than young radio-quiet pulsars, so this reflects the E p vs ˙ E trend we report. Our target selection procedures may have tended to favor larger E p , to be addressed in future work. Because the E p -˙ E anti-correlation is fairly modest with substantial scatter, its slope does not change substantially when excluding LAT-induced discoveries.  \nWe thus confirm that our highly diverse pulsar sample is striking in its spectral homogeneity: almost all gamma-ray pulsars are 'GeV pulsars', with 0 . 5 < E p < 2 . 5 GeV and following the same ˙ E trends.  \nExceptions exist. Upon its EGRET discovery, Ramanamurthy et al. (1996) highlighted the low gamma-ray efficiency of PSR J0659+1414. Romani et al. (2011) cast it as the archetype of a hypothesized category of pulsars sub-luminous in gamma rays. It stands out here as having very low E p and d p (Figure 20) and the lowest ϵ 37 , hi (the isolated dot at lower left in Figure 27).  \nAnother clear exception is the Crab: although its measured E p = 90 ± 80 MeV and spectral shape might be construed as a high-˙ E continuation of the trend in Figure 28, its true SED peak occurs in the soft gamma-ray band (Kuiper & Hermsen 2015). A further two pulsars, PSRs B1509 -58 and J1846 -0258, have spectral shapes similar to GeV pulsars but with values of E p well outside the continuum of GeV pulsars. Both pulsars have ˙ E ≈ 10 37 erg s -1 and little or no LAT signal above 100 MeV. Kuiper & Hermsen (2015) highlight them as archetypes of a pulsar population with SED peaks in the MeV range, and show 18 high-˙ E pulsars with SEDs rising in the keV range. Half of the 18 are in this catalog, as are 2 of 5 similar pulsars they discuss. PSR B1509 -58 has the best-defined MeV peak. They argue that some of the others may not be detected by the LAT because their SEDs also peak below the LAT's range, as for PSR J1846 -0258, at the very limit of our instrument's sensitivity (see Section 7.2). However, the SEDs shown by Kuiper & Hermsen (2015) leave open an alternate explanation: PSRs J0205+6449 and J2229+6114 indeed have keV slopes similar to PSR J1846 -0258, yet are bright up to a GeV. The as-yet undetected pulsars may simply be too faint, located as they are in very high background regions. Phase-connected rotation ephemerides for a few of the radio-quiet, noisy pulsars are unavailable. Others may be GeV pulsars, but with gamma-ray beams that simply miss the Earth. Harding & Kalapotharakos (2017) further discuss the MeV pulsars.  \nTo see whether the paucity of E p < 500 MeV pulsars reflects magnetospheric emission processes or stems from observational limits, we start with a tally of LAT candidate sources. There are 238 unidentified non-variable LAT sources with E p < 500 MeV, undetected beyond 10 GeV, and significantly curved (well fit with a LogParabola function). Of these, 122 have localisation ellipse semi-major axes with 95% C.L. < 10 arcmin, suitable for radio searches (see Table 4). Restricting the search to the 87 sources with | b | < 10 · would favor young, energetic pulsars. Blind searches for radio-quiet pulsars require a large gamma-ray count rate, even more so if timing noise, glitch rates, or background are high. Unfortunately, only 4 (25) of the 238 targets have TS > 200 (100). The 4 sources are 4FGL J0340.4+5302, J0426.5+5434, J1749.8-0303 and J1805.7+3401. All were searched several times for radio pulsations but none have been reported to date.  \nMany high-˙ E pulsars have already been found in radio and X-ray searches, effective at finding such sources because high-˙ E pulsars tend to have broad radio beams, and because searches can target the supernova remnants near to which such young sources likely remain. But high-˙ E pulsars are young, and therefore rare, limited by the Galactic supernova rate. What is rare is far, and thus faint. The subset of high-˙ E pulsars for which only the spectral tail reaches the LAT's sensitivity range is thus even rarer. Kuiper et al. (2018)'s detection of PSR J1846 -0258 importantly demonstrates the LAT's ability to see such objects. Learning how many more E p < 500 MeV pulsars exist will likely have to wait for a sensitive, large field-of-view MeV instrument such as AMEGO-X (Caputo et al. 2022) and ASTROMeV (De Angelis et al. 2021). It would  \nallow the missing SED segments for the known and candidate lowE p pulsars to be filled in, and to search for pulsations from new ones. But again, the ultimate low E p population is likely to be small compared to that for GeV pulsars, thus we propose that the trends shown here do capture the spectral character of most gamma-ray pulsars.", '8.3. Light Curve Trends': "The 237 energy-integrated profile fits in this catalog generally confirm the broad categories of light curve shapes known since EGRET, consolidated in 1PC and 2PC: a majority of gamma-ray pulsars have two principal peaks separated by ∆ ≈ 0 . 4 ± 0 . 15 rotations (see Figure 10). For radio-detected pulsars, the first peak generally trails the radio pulse by δ ≈ 0 . 2 ± 0 . 2 in phase. Gamma-ray emission from near or beyond the light cylinder can produce the observed anti-correlation between ∆ and δ , as discussed in 2PC and explicitly reproduced by the models of Kalapotharakos et al. (2014, 2023). About 40 pulsars have a single broad gamma-ray peak. Very few pulsars deviate from these trends: only two young pulsars in Table 9 are flagged as 'complex', appearing to have a third peak half-way between the 'standard' two peaks (one is EGRET pulsar B1055-52). The number increases to 7 for the MSPs in Table 10. Overall, the MSP and young pulsar profiles are more similar than they differ.  \nOn the other hand, 167 pulsars have enough accumulated photons to allow profile fits in multiple energy sub-bands, fits not done for 2PC. For the first time, we quantify how profiles evolve with energy for a large sample. Figures 13 and 14 highlight the 28 brightest profiles, fit in at least 6 energy sub-bands. Reminiscent of the spectral results, profile morphology and evolution are remarkably homogeneous given the large diversity of pulsar types and environments. Here follow two examples.  \nIn the second panel from the top 20 , the energy-resolved fit summaries show the phase evolution for P1 and P2, and P3 for the few cases where it is a strong component. Changes in the peak separation ∆ are nearly always small, ≲ 0 . 1 at most. For < 5 pulsars, such as PSR J2229+6114 shown in Figure B-4, ∆ appears to change with energy. However, inspection of the profiles in Figure B-4 shows that individual profile components vary to accommodate slight pulse shape changes, but ∆ in fact varies little.  \nThe summary plots also illustrate the other example. Whereas P3 for the Vela pulsar increases by 0 . 1 in phase over the LAT energy range, visibly detaching from P1 and moving across the bridge region, this turns out to be the only known pulsar where a gamma-ray peak phase varies substantially with E γ . Three pulsars, PSRs J0633+1736 (Geminga), J1747-2958, and J1836+5925 ('the next Geminga'), show P3 evolution in the summary plots. Again, inspection of the profiles shows it comes from weak broadening of P2 rather than a peak drift as seen for Vela.  \nSuch weak dependence of pulse phase on energy does not necessarily imply that the region illuminating an observer is uniform. Relativistic aberrations and time delays create caustics such that photons emitted at different distances from the neutron star remain in narrow phase ranges. Thus, whether different electron energies occur at different magnetospheric depths or not, we could see little light curve evolution with photon energy. Nevertheless, since the spectral energy distributions suggest populations of gamma-ray emitting electrons with a narrow range of maximum energies for the lowest ˙ E pulsars, broadening gradually for systems with more available power, the light curves also encourage vision of emission zones along caustic viewing lines that change little with electron energy. Figure 22 attempts to relate spectral and profile properties: it suggests that when the line of sight cuts across the emission region such that the viewer sees two closely-spaced peaks (small ∆), then the line of sight sometimes samples a region with the highest electron energies (large E p ).  \nIndeed, how the line of sight intersects the emission region is key. Beginning with Watters et al. (2009), and continuing as mastery of magnetospheric processes evolves (e.g. Venter et al. 2009; Contopoulos & Kalapotharakos 2010; Bai & Spitkovsky 2010; Arka & Dubus 2013; Kalapotharakos et al. 2014; Harding 2016; Cerutti et al. 2016; Kalapotharakos et al. 2018; Philippov & Spitkovsky 2018; P'etri 2019; Cerutti et al. 2020; Kalapotharakos et al. 2023), modelers have provided 'atlases' allowing comparison between observed and calculated pulse profiles. Sky maps show beam intensity versus neutron star latitude and longitude, that is, as a function of the inclination ζ between the rotation axis and the line of sight, and the rotational phase ϕ , with other parameters fixed. A cut across a sky map yields the profile sampled by a given line of sight. The dominant 'other' parameter is the magnetic inclination α , but e.g. period or magnetic field strength have more or less importance for different models. Some theories favor emission from plasmoids in current sheets beyond the light cylinder, where magnetic reconnection accelerates electrons to high energies to then radiate synchrotron gamma rays (P'etri & Lyubarsky 2007; Cerutti et al. 2020). Other models favor electric field acceleration near the current sheet and curvature radiation. In either case the picture is quite different from the various 'gaps' that previously dominated the field (Philippov & Kramer 2022), and the resulting sky maps differ in consequence. The particle-in-cell (PIC) calculations of the relativistic magnetohydrodynamic (MHD) processes consume substantial computer resources. Realistic atlases are currently out of reach, with progress in kinetic plasma simulations expected in the coming years.  \nFigure 29 histograms the gamma-ray duty cycle D , the total phase for which the fitted profile is > 25% of the range from its baseline to the maximum. (The online catalog also includes D with 10% and 50% thresholds, and psrcat includes radio pulse width.) Hou et al. (2014) highlighted some pulsars with large duty cycle, and discussed how detection sensitivity depends on D . With two exceptions, the young  \nradio-loud pulsars have D < 0 . 4. In contrast, 3 times more young radio-quiet pulsars have 0 . 3 < D < 0 . 5 as loud ones.  \nRadio detection occurs when | ζ -α | ≲ ρ , that is, along a broad diagonal swath of the plane of neutron star inclinations ( α, ζ ) (see Figure 3 of Johnston et al. 2020). The opening angle of the radio cone is ρ = 3 √ πh/ 2 Pc , for emission height h . Gamma-ray detection favors the α > 30 · and ζ > 30 · part of the plane. Finally, long duty cycles arise when the line of sight slices a skymap in a ζ range with emission extended in ϕ , but these 'flat' beam shapes favor small α (see e.g. Figure 15 of Kalapotharakos et al. 2018). The tallies in Figure 29 are in good qualitative agreement with this simple geometric beam model. Thus, study of the high duty cycle pulsars may help constrain ( α, ζ ) and/or discriminate between models. The small light cylinders characteristic of MSPs make for broader beams, the likely reason for the MSP D tail in Figure 29.  \nIf D could constrain the 'beaming fraction' f Ω used to calculate luminosity L γ (Eq. 24), the dispersion in the L γ versus ˙ E could be reduced. Large D pulsars, with their distinctive skymaps, could be particularly useful. But they are rare: the α range producing the right beam shape is limited, and the line of sight must fall into a narrow ζ range. PSR J1857+0143, the young radio-loud outlier in Figure 29, may be an example. Serylak et al. (2021) show its polarization position angle versus phase at 1 . 4 GHz, but unfortunately scattering excludes a rotating vector model (RVM) estimate of ( α, ζ ) at this frequency while its narrow radio pulse width poorly constrains the RVM at higher frequencies. Table 6 points out unidentified LAT sources co-located with known pulsars for which we have not seen gamma-ray pulsations: some may be D ≈ 1 pulsars and could benefit from a ( α, ζ ) study.  \nThe differences in spectral hardness (Figure 13) between the peaks nearer (P1) and farther (P2) from the magnetic pole, from the observer's point of view, is another powerful discriminant, mentioned above. Harding et al. (2021) and Barnard et al. (2022) address this specifically. Early LAT studies measured spectra in as many phase bins as the statistics would allow (Abdo et al. 2010c,g,h; DeCesar 2013). The observed phase evolution of the spectral parameters is smooth, and the range of parameter values is moderate (the cutoff energy can vary by a factor of 2). A simple analysis, measuring on-peak spectral shape in the peak phase intervals we provide, is easier to perform than a full phase-resolved spectral analysis that samples light curves at many points. Furthermore, a larger sample of pulsars is bright enough for on-pulse measurements than for full phase-resolved analysis, yet a simple analysis is well-suited to the current model accuracy. Such future work is encouraged.", '8.4. Luminosity and Pulsar Populations': "Our large gamma-ray pulsar sample, along with the known radio pulsars, is a starting point for a census of the total Milky Way population. Characterizing the population sheds light on many open questions. These include the supernova rate and neutron star evolution. The large number of 'spider' binary MSPs caused a renais-  \nFigure 29. Duty cycle D of the gamma-ray pulse profiles (fraction of a rotation that the fitted profile exceeds 25% of the maximum amplitude, measured from the background level). More young radio quiet pulsars have large duty cycles than do the young radio loud pulsars.  \n<!-- image -->  \nsance in the study of binary evolution and pulsar recycling (Strader et al. 2019; van der Merwe et al. 2020; Swihart et al. 2022b). The number and luminosity of high-energy pulsars helps determine the relative contributions of pulsars and putative dark matter annihilation to the observed cosmic positron spectrum (Venter et al. 2015; Cholis & Krommydas 2022) and, via pulsar wind nebulae, to the leptonic cosmic ray spectrum generally. Applying emission models to a synthetic population and then simply tallying the observed numbers of, for example, radio-loud versus radio-quiet pulsars tests predictions of beam geometry (Johnston et al. 2020).  \nThe strong dependence of L γ on ˙ E (Figure 23) and the LAT's sensitivity (Figure 25) are key inputs. Key tests of a synthesis's fidelity are the pulsar latitude distribution (Figure 3) and the cumulatative logN-logS distribution of our pulsar sample, shown in Figure 30. The slopes in the logN-logS plot are related to the spatial distributions of the sources: a flux-limited sample of a population of sources having a narrow luminosity distribution has slope -3/2 if the population is isotropic, and -1 if the distribution is disk-like. Indeed, in Figure 30 the MSPs have a steeper slope, being more isotropic, while the young pulsars have a slope close to -1. The turnover at low fluxes is the detection threshold, lower for objects away from the intense background  \nFigure 30. The cumulative number of pulsars detected with > 100 MeV flux over G 100 , 'logN-logS'. The lines show slopes of -1 and -3 / 2.  \n<!-- image -->  \nin the plane. Distance uncertainties do not affect the logN-logS plot. However, an extrapolation from the observed to the unresolved population needs to take selection biases into account. One example is that most of our young radio-loud pulsars were discovered in radio surveys near the plane - results sensitive to details of the highlatitude population may suffer from bias. Another is that target sources for gammaray blind searches for radio-quiet pulsars have a higher threshold than shown in Figure 25, explained in Section 3.4.  \nThe scatter around the L γ versus ˙ E correlation (or, equivalently, of efficiency versus ˙ E shown in Figure 24) is a combination of physics and uncertainties. The 'physics' bias includes beaming, since f Ω = 1 in our L γ computation. We use a canonical moment of inertia I = 10 45 g cm 2 for all pulsars, but the true average is probably twice that, since neutron star radii are closer to 14 km than to 10 km, shifting the xaxis of all ˙ E plots to the right. Furthermore, especially MSPs have a range of masses, of order only 60%, neglected in our single I value. Luminosity may also depend on parameters other than ˙ E , such as the inclination angles, which would contribute to the scatter. The correlation is stronger at higher spindown powers, and is less pronounced for the MSPs. It remains to be understood whether this is due to high energy processes 'shrinking' as ˙ E tends toward the 'deathline' near 10 33 erg s -1 , or is because e.g. of the compactness of the MSP magnetosphere.  \nHighly polarized radio beams generally occur only for pulsars with ˙ E > 10 34 erg s -1 (Weltevrede & Johnston 2008; Johnston & Kerr 2018). Smith et al. (2019) suggested it is an unlikely coincidence that both radio polarization and gamma-ray emission  \nrequire the same minimum spin parameters for independent reasons - that these two phenomena share the same deathline may mean that the same electron populations produce both signals.", '8.5. Towards 4PC and 400 gamma-ray pulsars': 'A small number ( ∼ 20) of our gamma-ray pulsars are below the 4FGL catalog detection threshold. Thus we expect that the bulk of future pulsed discoveries will correspond to the thousand or more currently unidentified LAT sources. Requiring non-variability and a pulsar-like spectrum reduces this to a few hundred. It thus appears unlikely that the LAT will ever more than double the size of the current sample.  \nHowever, another hundred gamma-ray pulsars is possible, and dozens are likely. Machine-learning methods applied to the 4FGL sources typically find ∼ 70 pulsar-like unassociated objects with | b | > 10 · (Mirabal & Bonaca 2021). Future high-sensitivity radio surveys with SKA are predicted to find 27 , 000 young pulsars and 3 , 000 MSPs (Keane et al. 2015). Radio timing solutions for those co-located with LAT source would allow gamma-ray folding, at least for the years following the radio discovery. A radio position, period, and orbital period for the binaries would improve the sensitivity and decrease the computing cost for gamma-ray blind searches. Without waiting for SKA, TRAPUM at MeerKAT and FAST are already finding new pulsars. Sustained timing of the new discoveries will require thoughtful management of available resources. The insights we have achieved in this work into the shapes of pulsar spectra will help quantify how many of the unidentified LAT sources are likely to be pulsars without pulsed detections.  \nWe intend to maintain the online catalog (see the footnote at the end of Section 1). We will add new gamma-ray pulsars as they accumulate. If new analyses are performed, e.g. on-peak spectral analysis, we will add those results.  \nPulsar observations enable investigations over a spectacular range of topics, from fundamental physics to variations of the interstellar medium, including several aspects of stellar evolution, and more. Our gamma-ray pulsar sample is so large that even subsamples selected for specific studies have more members than did the entire collection acquired by Compton Gamma-ray Observatory , over twenty years ago. Discoveries enabled by exploiting the sample will expand for many years to come, at the same time that Fermi LAT keeps finding more, and better characterizes those already known.', 'A. RADIO TIMING FOR THIS WORK': "PSR J0653+4706: PSR J0653+4706 was discovered with LOFAR (van Haarlem et al. 2013) at 135 MHz in a targeted pulsar survey towards unidentified Fermi LAT gamma-ray sources (Pleunis et al. 2017; Bassa et al. 2017b). The survey used semicoherent dedispersion to minimize the effects of dispersive smearing at these low frequencies while retaining sensitivity to MSPs (Bassa et al. 2017a). The pulsar has a 4.76 ms spin period and a dispersion measure of DM = 25 . 74 ± 0 . 14 pc cm -3 , placing it at a distance of ∼ 900 pc with both the NE2001 and the YMW16 models. It is in a 5.84 day binary system with a M c ≳ 0 . 21 M ⊙ binary companion (Bassa et al. 2018). The timing solution used in the catalog is obtained from LOFAR timing observations covering a 2.8 yr timespan. Subsequently, a 5 . 2 year solution using Jodrell Bank Observatory L-band observations confirmed the model.  \nThe LOFAR data indicate a flux density of 5 ± 2 mJy at 150 MHz. An 1800 s Green Bank Telescope observation at 820 MHz, with 200 MHz bandwidth, yielded a 37 . 2 σ detection. For 2 polarizations, gain G = 2 K Jy -1 , system temperature T sys = 29 K, and a 12% duty cycle seen in both the 150 and 820 MHz profiles, we obtain a flux density of 260 µ Jy, calculated using the PRESTO s um profiles.py routine.  \nAnalysis of 101 positive detections at 1400 MHz with the Lovell telescope at Jodrell Bank Observatory yields flux densities of 89 ± 11 and 10 ± 3 µ Jy for the main peak and interpulse, respectively. The uncertainties are statistical, and include variations due to interstellar scintillation. The total is S 1400 = 99 ± 23 µ Jy, having added the 20% systematic calibration uncertainty in quadrature. The spectral index obtained from these three measurements agrees well with the average value of 1 . 7 of the pulsar population.  \nThe JBO ROACH timing data show no significant pulsar proper motion. We found PMDEC = -1 . 5(2) and PMRA = 0 . 5(1 . 8). The Doppler correction is hence 0 . 98 < (1 -ξ ) < 1, with (1 -ξ ) = ˙ E int / ˙ E as in Section 4.3.  \nPSR J1023+0038: PSR J1023+0038 was the first transitional millisecond pulsar to be discovered (Archibald et al. 2009), and it underwent a transition from a rotationally-powered MSP state to an accretion-powered state in 2013 (Stappers et al. 2014). A faint source coincident with the pulsar was detected by the LAT in the pre-transition data (Nolan et al. 2012), and the flux from this source increased five-fold when accretion began (Stappers et al. 2014).  \nUsing the radio timing ephemeris of Archibald et al. (2013), we folded the pretransition data (from the beginning of the LAT data to 2013 June), discovering significant pulsations, with H = 62. This ephemeris contains a single orbital period derivative to account for long-term trends in the orbital period, but does not account for shorter-timescale orbital phase variations seen in radio timing observations. We  \nrefined this solution using an unbinned timing analysis on the LAT data, also using only one orbital period derivative, resulting in an increased H = 92.  \nWe also folded the post-transition data using a timing model built from X-ray timing by Jaodand et al. (2016) and X-ray and optical timing from Illiano et al. (2023). We fit a polynomial function to the orbital phase deviations measured by Illiano et al. (2023), and folded the LAT data using this orbital model and the spindown model of Jaodand et al. (2016). No pulsations were detected. We also searched for pulsations with slightly different spin periods and spin-down rates, and with slight shifts in orbital phase, using a sliding window approach, but again did not detect any significant pulsations.  \nThis non-detection implies that the flux increase due to accretion is not an enhancement of the pulsed flux (seen in other non-accreting redbacks, e.g. An et al. 2020; Clark et al. 2021) , but is instead an additional unpulsed component. If the five-fold increase in gamma-ray flux is indeed due to unpulsed emission from the accretion process, and the gamma-ray pulsations remain unaffected by the accretion, then the pulsed fraction would be reduced to 20% of the pre-transition level, and this would leave the pulsations undetectable to our analysis. Like PSR J1227 -4853 (Johnson et al. 2015), it therefore remains unclear whether or not this system emits gamma-ray pulsations in the accretion-powered state.  \nPSR J1827 -0849: Applying the methods of Clark et al. (2017) and Section 3.4 to LAT data around the position of a steep-spectrum radio source coincident with unidentified LAT source 3FGL J1827.6-0846, found by Frail et al. (2016), we discovered gamma-ray pulsations with spin period 2 . 24 ms in 2016. The radio positional uncertainties of ∼ 2 '' have greatly reduced the computing cost and allowed us to increase our sliding coherence window length to 2 23 s (or ∼ 97 days) and still perform the search on the ATLAS cluster 21 in about 1 day per source.  \nThe associated radio source has flux density 229 . 9 ± 14 . 8 mJy in the GMRT 150 MHz All-Sky Survey, and 777 ± 126 mJy in VLSS r at 74 MHz. NVSS has no detection, giving S 1400 < 0 . 8 mJy per beam. Table 8 lists an insignificant proper motion estimate from LAT timing as per Section 3.5. C.J. Clark and co-workers are studying other LAT unidentified sources co-located with steep-spectrum radio sources and will detail their methods and results in future work.  \nPSR J1833 -3840: This 1 . 87 ms pulsar was discovered in 2015 by F. Camilo, in observations of unidentified LAT sources using the Parkes radio telescope, using the methods reported by Camilo et al. (2015). Later, Frail et al. (2018) independently highlighted it in a search for steep radio spectrum pulsar candidates in LAT unidentified sources using the GMRT 150 MHz all-sky survey (TGSS ADR1). Fewer than 10% of known pulsars have such high spindown power, ˙ E = 1 . 1 × 10 35 erg s -1 , in the Milky Way or in globular clusters. The orbital period is 0 . 9 days. In 2020,  \nthen-undergraduate Camryn Phillips phase-connected GBT timing observations during PINT code development she was conducting. We phase-folded the gamma-ray data using her ephemeris and obtained the pulsed detection reported here.  \nPSR J1852-1310: The P = 4 . 31 ms isolated pulsar PSR J1852 -1310 was discovered in a 1900 s pointing on MJD 57793 towards a LAT unidentified source with a pulsar-like spectrum, at 800 MHz using the GUPPI backend on the Green Bank radio telescope. It was initially named PSR J1852 -13. Routine timing observations at 1400 MHz with the Lovell telescope at Jodrell Bank Observatory began on MJD 57818. As of MJD 59927, 124 arrival times have been acquired using ROACH, allowing a phase-connected rotation ephemeris, provided with the online material. The spindown power is ˙ E = 4 . 95 × 10 33 erg s -1 . For DM = 44 . 95 ± 0 . 01 pc cm -3 in the direction of this pulsar, the YMW16 model gives distance 1 . 27 ± 0 . 5 kpc. The integral energy flux of co-located 4FGL J1852.2 -1309 thus gives L γ = (5 . 8 ± 1 . 0) × 10 32 erg s -1 , for a gamma-ray efficiency of 11%.  \nGBT timing observations at 1400 MHz between MJD 58698 and 59050 yielded 49 times of arrival. Combining these with the JBO data slightly improves the MSP's proper motion measurement, and we obtain PMDEC = -6 . 6 ± 1 . 4 and PMRA = -1 . 84 ± 0 . 25. The Doppler correction is hence (1 -ξ ) = 0 . 91, as defined in Section 4.3, such that ˙ E i nt = (4 . 6 ± 0 . 2) × 10 33 erg s -1 .  \nGamma-ray phase folding 14 . 5 years of data using the JBO timing model revealed pulsations with test statistic H = 122 using the simple weights parameter µ w = 4 . 3. Model weights increase this to H = 133. H increases steadily for the 8 . 5 years before the start of the timing model, indicating that the MSPs spindown rate is stable, well-described using only spin frequency and its first derivative. The radio source is bright (250 σ in the discovery observation) and the profile has sharp features, hence this MSP may rate inclusion in gravitational wave searches.  \nPSR J1857+0943: This 5 . 36 ms pulsar has long been known for its stable rotation rate (Segelstein et al. 1986). Gamma-ray phase folding using the rotation ephemeris provided by Reardon et al. (2021) yielded H-test > 30, using simple weights. Located in a high background region, ( l, b ) = (42 . 2 , 3 . 12) degrees, the gamma-ray profile is ill-defined but appears to be dominated by a single narrow peak. The timing parallax distance from Table 8 and the integral energy flux of the co-located LAT source 4FGL J1857.2+0941 give L γ = (1 . 27 ± 0 . 31 0 . 6 -0 . 9 ) × 10 33 erg s -1 , where the first uncertainty comes from the flux measurement and the second from the distance, about 30% of the spindown power.  \nPSR J1901 -0125: The same procedure as described for PSR J1827 -0849, above, led to the 2016 discovery of 2 . 79 ms gamma-ray pulsations. F. Camilo then re-analyzed GBT radio data from previous observations of the LAT source position using the LAT rotation ephemeris and discovered radio pulsations and measured the DM. Table 3 lists S 1400 = 2 . 50 mJy, extrapolated from the Frail et al. (2016) GMRT flux density of 362 ± 15 mJy at 150 MHz and radio spectral index -2 . 23. Table 8 lists  \na 2 . 4 σ proper motion measurement using the methods of Section 3.5. Clark and coworkers are studying other LAT unidentified sources co-located with steep-spectrum radio sources and will detail their methods and results in future work.  \nJ2256 -1024: Crowter et al. (2020) found µ = 9 ± 3 mas yr -1 using GBT timing data. NRT has also monitored this pulsar, but with a substantial gap between the end of the GBT data and the beginning of the NRT observations. This gap and limited frequency coverage cause a strong degeneracy between the spindown parameters, time-varying dispersion measure, and orbital frequency evolution. We therefore perform a hybrid analysis, using the methods of Section 3.5, in which we use LAT data to determine the 'long-term' parameters (F0, F1, FB0, FB1, and FB2). We then fix these parameters and fit for the astrometry using the radio data, finding a proper motion of PMRA = 2 . 80 ± 0 . 05, PMDEC = -5 . 30 ± 0 . 11 from which µ = 6 . 0 ± 0 . 1 mas yr -1 . A systematic error comparable to or larger than the statistical uncertainty is likely, due to the inhomogeneous data used in the analysis.", 'B. SAMPLE LIGHT CURVES, PROFILE EVOLUTIONS, AND SPECTRA.': 'Figures B-1 to B-4 show examples of LAT pulsar light curves in different energy bands, with the radio pulse profiles overlaid when available. See Section 5 for details about the construction of these light curves. All the other pulse profiles are available in the supplementary online material, detailed in Appendix D.  \nFigures B-5 to B-8 show examples of energy-resolved light curve fitting analysis summaries, as described in Section 5. The plots for the remaining pulsars are in the supplementary online material.  \nFigure B-9 shows sample spectral energy distributions (SEDs) provided as part of 4FGL-DR3 22 . The supplementary online material for this pulsar catalog includes a compilation of the SEDs of the gamma-ray pulsars, in PDF format. It also includes the spectral parameters with their covariance matrices and the data points necessary to reproduce the SEDs, and a python script that performs the task.', 'C. ADDITIONAL NOTES ON SPECTRAL ANALYSIS': "We study here the distribution of the model parameters and derived physical quantities as they vary among the 116 pulsars which have both b -free and b = 2 / 3 fits available. These results are presented in Figure C-1. In brief, they are  \n- · Γ, the spectral index at the reference energy E 0 , is stable, varying by < 10% over the effective range in b .\n- · The curvature at the reference energy, d , varies more than Γ, by up to ∼ 50%. If considering b < 1, the total range is only about ± 20%.\n- · Γ 100 varies almost linearly with b before saturating at b = 1 . 2. The shift is up to 1 unit harder for pulsars with b near 0, and 0 . 5 units softer for (mostly) MSPs with b ≥ 1 . 2. Thus, clearly, freeing b profoundly alters the inferred low-energy spectral index.\n- · The curvature at the peak, d p , is stable if the intrinsic value of b > 0 . 5. For pulsars with spectra more consistent with small b , the b -free pulsars become sharper, with a difference in log-width ranging up to 2 × . It is not (yet) clear if this difference primarily 'fixes' the true width of the peak or delivers a better low-energy index.\n- · The peak energy E p varies by < 10-15% in the 'core' b range and < 25% overall. The largest deviations are for the MSPs with fit values b > 1 which move to higher E p .  \nIn summary, the physical properties measured at the peak, d p and E p , are robustly measured in both b -free and b = 2 / 3 fits and so we can make inferences on their distribution using the entire pulsar sample. On the other hand, the low-energy spectrum is unreliable, and values inferred from the b = 2 / 3 model could be very far from the true values.  \n/negationslash  \n/negationslash  \n/negationslash  \n/negationslash  \n/negationslash  \nFigure C-1. Comparisons of base and derived parameters between the b = 2 3 and b -free populations. See main text for further discussion.  \n<!-- image -->  \nb", 'D. DESCRIPTION OF THE ONLINE CATALOG FILES': "Online material is available through this Journal as well as at https://fermi.gsfc. nasa.gov/ssc/data/access/lat/3rd PSR catalog/. Here we describe the files posted as of publication. Additional material may be added in the future. The 'xxx' in the file names is a multi-digit timestamp of the file creation date. Contact the authors to request missing information.  \n3PC Catalog+SEDs xxx.fits, 3PC Catalog xxx.xls : The pulsar catalog 'flat file', in FITS and spreadsheet formats. Each has 294 rows, one for each gamma-ray pulsar. The 112 columns contain basic pulsar properties (name, position, spin parameters, distance, radio flux density, et cetera) as well as the spectral and profile parameters described in this article. The FITS version includes a subtable ('HDU') with the detailed spectral fit parameters and covariance matrix needed to reproduce the Spectral Energy Distributions (SEDs), provided in 3PC SEDs.tgz (individual plots) and 3PC AllSEDs.pdf (all plots concatenated into a single file). 3PC SEDPlotter.tgz contains a python script that reads 3PC Catalog+SEDs xxx.fits to generate the SED plots. Most Figures and Tables in this article were created from the FITS and spreadsheet files: the python script Example 3PC CatalogReader.py is a simple example. README 3PC Catalog+SEDs xxx.txt provides more detail .  \n3PC ProfileData xxx.tgz : Contains two files for each pulsar, one in FITS format and the other in simple text (ascii), containing the weighted gamma-ray phase histograms and, when they exist, the fit to the histogram and/or the radio profile. The python script Example 3PC PlotProfileData.py illustrates how to read the files to reproduce the plots provided in 3PC ProfilePlotsPDF xxx.tgz and 3PC ProfilePlotsPNG xxx.tgz . The plots of gamma-ray peak evolution with energy, as in Figures 12 and B-5 to B-8, are provided in 3PC LC Eres pdf xxx.tgz, 3PC LC Eres png xxx.tgz . README ProfileData.txt provides details.  \n3PC TimingModels xxx.tgz : A compressed archive of the 294 rotation ephemerides ('.par files') used to calculate the gamma-ray rotational phases using software such as Tempo2 or PINT. README 3PC TimingModels.txt provides details.  \nPhased, weighted FT1 FITS files: available at https://heasarc.gsfc.nasa.gov/ FTP/fermi/data/lat/catalogs/3PC/photon/. Two subdirectories called 3deg 50MeV and 15deg 20MeV each contain one file for each gamma-ray pulsar, in the standard Fermi LAT data format, with additional columns PHASE and MODEL WEIGHT. The 277 files in 3deg 50MeV include photons with energies > 50 MeV, within 3 · of the pulsar position, and were used to create the gamma-ray phase histograms. They occupy roughly 30 Gb of disk space, with a large variation in single file size depending on the background intensity (generally decreasing away from the Galactic plane) and pulsar intensity. The files in 15deg 20MeV are on average 15 × larger, including photons with energies > 20 MeV within 15 · of the pulsar position. They are suitable for spectral analyses (e.g. using the gtlike Fermi tool) or to explore the  \nlowest LAT photon energies. README 3PC PhotonData .txt describes the FT1 files in more detail.  \n3PC sensitivity xxx.fits : The all-sky Fermi LAT pulsar sensitivity, in FITS format, useful for e.g. population syntheses. Example 3PC SensitivityMap xxx.py is a python script that reads the FITS file and generates a plot as in Figure 25. README 3PC SensitivityMap xxx.txt provides details.  \nFinally, we provide HTML files for each pulsar, to allow a visual summary of the information for each gamma-ray pulsar, named 3PC/J0007+7303 LAT.html and so on.  \nAcknowledgments: The Fermi -LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat 'a l'Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucl'eaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden.  \nAdditional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d' ' Etudes Spatiales in France. This work performed in part under DOE Contract DE-AC0276SF00515.  \nThe Parkes radio telescope is part of the Australia Telescope which is funded by the Commonwealth Government for operation as a National Facility managed by CSIRO. The Green Bank Telescope is operated by the National Radio Astronomy Observatory, a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Arecibo Observatory is part of the National Astronomy and Ionosphere Center (NAIC), a national research center operated by Cornell University under a cooperative agreement with the National Science Foundation. The Nan¸cay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS). The Lovell Telescope is owned and operated by the University of Manchester as part of the Jodrell Bank Centre for Astrophysics with support from the Science and Technology Facilities Council of the United Kingdom. The International LOFAR Telescope and the Westerbork Synthesis Radio Telescope are operated by ASTRON, the Netherlands Institute for Radio Astronomy. This work made extensive use of the ATNF pulsar catalog (Manchester et al. 2005). IG, AKH, and DAS thank the ISSI (Bern, Switzerland) for financial support of Team 459 meetings that helped to improve the present work.  \nFacilities:  \nFermi(LAT)  \nSoftware: Tempo2 (Hobbs et al. 2006), PINT (Luo et al. 2021), astropy (Astropy Collaboration et al. 2018), numpy (Harris et al. 2020)"}