The Fundamental Plane Relation for Gamma-Ray Pulsars Implied by 4FGL

We explore the validity of the recently reported fundamental plane (FP) relation of gamma-ray pulsars using 190 pulsars included in the latest 4FGL-DR3 catalog. This sample number is more than twice as large as that of the original study. The FP relation incorporates 4 parameters, i.e., the spin-down power, $\dot{\mathcal{E}}$, the surface magnetic field, $B_{\star}$, the total gamma-ray luminosity, $L_{\gamma}$, and a spectral cutoff energy, $\epsilon_{\rm cut}$. The derivation of $\epsilon_{\rm cut}$ is the most intriguing one because $\epsilon_{\rm cut}$ depends on the proper interpretation of the available phase-averaged spectra. We construct synthetic phase-averaged spectra, guided by the few existing phase-resolved ones, to find that the best fit cutoff energy, $\epsilon_{\rm c1}$, corresponding to a purely exponential cutoff (plus a power law) spectral form, is the parameter that optimally probes the maximum cutoff energy of the emission that originates from the core of the dissipative region, i.e., the equatorial current sheet. Computing this parameter for the 190 4FGL pulsars, we find that the resulting FP relation, i.e. the gamma-ray luminosity in terms of the other observables, reads $L_{\gamma}=10^{14.3\pm 1.3}(\epsilon_{\rm c1}/{\rm MeV})^{1.39\pm0.17}(B_{\star}/{\rm G})^{0.12\pm 0.03}(\dot{\mathcal{E}}/{\rm erg\;s^{-1}})^{0.39\pm 0.05}{\rm ~erg\;s^{-1}}$; this is in good agreement with both the empirical relation reported by Kalapotharakos et al. (2019) and the theoretically predicted relation for curvature radiation. Finally, we revisit the radiation reaction limited condition, to find it is a sufficient but not necessary condition for the theoretical derivation of the FP relation. However, the assumption of the radiation reaction limited acceleration reveals the underlying accelerating electric field component and its scaling with $\dot{\mathcal{E}}$.


INTRODUCTION
The Fermi Gamma-Ray Space Telescope has provided an unprecedented level of information, deepening our understanding of the γ-ray pulsars considerably. The observational data compiled in the Second Fermi Pulsar Catalog (2PC, Abdo et al. 2013) showed various correlations especially between spectral properties and observed or derived parameters, such as cutoff energy, cut , with light cylinder (LC) magnetic field, B LC and spectral photon index, Γ, with spin-down power,Ė, that actually pointed out the operational regime of particle ac-§2 below). Remarkably, the observed FP is close to the theoretical relation L γ ∝ 4/3 cut B 1/6 Ė 5/12 (2) corresponding to curvature radiation (CR) emission in the radiation reaction limited (RRL) regime that occurs at the equatorial current sheet in the vicinity of the LC. The FP was independently confirmed by Ploeg et al. (2020) who studied MPs, exploring the possibility that unresolved MPs are responsible for the observed GeV photon Galactic center excess (Abazajian 2011; also see, Brandt & Kocsis 2015;Bartels et al. 2016;Hooper & Mohlabeng 2016;Ackermann et al. 2017;Gonthier et al. 2018;Ploeg 2021;Berteaud et al. 2021;Gautam et al. 2021).
(1), relying on the 2PC data. More specifically, we used 88 pulsars with reliable (i.e., published) spectral values (i.e., L γ , cut ) out of the 114 pulsars that were totally compiled in 2PC while currently, more than 270 γ-ray pulsars have been detected 1 .
On the one hand, the exponents in the expressions in Eqs. (1)-(2) indicate greater sensitivity of L γ on cut than B andĖ. On the other hand, the B andĖ values of Fermi pulsars cover a range spanning several orders of magnitude while the corresponding cut values cover a limited range of less than an order of magnitude. Thus, the determination of cut becomes vital for a reliable derivation of the FP. The 2PC cut values, employed in K19, had been derived adopting a power law with a pure exponential cutoff spectral function model. Nonetheless, the Fermi fits indicate that for the brightest pulsars, a sub-exponential function model can statistically describe better the phase-averaged spectra.
In this paper, we first explore the behavior of these different function models, uncovering the parameter values that best correlate with the cutoff energies corresponding to the highest-energy pulsar emission. Then, we update the FP using the spectral analysis data from the 12-year incremental version of the fourth Fermi-LAT catalog of point γ-ray sources (4FGL, Abdollahi et al. 2020; 4FGL-DR3, for Data Release 3, Fermi-LAT collaboration et al. 2022) combined with data from the Australia Telescope National Facility (ATNF) Pulsar Catalog (Manchester et al. 2005) 2 that expand by more than two fold (compared to our previous exploration) the sample of pulsars. Finally, we revisit and further generalize the necessary and sufficient underlying theoretical conditions put forth in K19.
The structure of the paper is as follows. In Section 2, we discuss the behavior of the spectral fitting function models and the related parameters. In Section 3, we present the expanded sample of Fermi pulsars taking into account the 4FGL-DR3 and ATNF data and update the FP relation of γ-ray pulsars. In section 4, we review the RRL requirement and discuss the possible effects on the FP. Finally, in the last section, we discuss our results and summarize our conclusions.
2. FINDING THE CUTOFF ENERGY Ideally, a fitting function model should describe the spectra and provide meaningful parameter values that correspond to the underlying physics. A model that describes the phase-averaged energy-weighted differential spectral energy density (SED) adequately for the vast majority of Fermi pulsars reads where is the photon energy, dN is the number of photons in energy bin ( , + d ), and Γ is the photon index. The parameter b describes the decay of the emitted power at large energies. The pure exponential cutoff corresponds to b = 1 while the brightest Fermi pulsars (e.g., Vela, Crab, Geminga) are well-described by a subexponential cutoff (i.e., b < 1). The apex energy, A , that signifies the photon energy of the SED maximum occurs at Thus, for fixed Γ < 2 and b values, A traces cut (being always proportional to it). However, the Fermi data have shown that Γ positively correlates with the spindown powerĖ (2PC). This implies that following the behavior of cut with A no longer reveals the corresponding cut behavior. Therefore, cut can remain constant withĖ and still A can decrease withĖ (as long as Γ increases withĖ). Moreover, for b ∼ 1, cut signifies the energy at which the rapid decline of the SED begins. However, when b starts deviating from 1, the corresponding cut value no longer discloses that characteristic photon energy scale. In Figure 1, we plot normalized -to the maximum value-mock SEDs color-coded according to b values that range from 0.5 to 1 and for Γ = 2/3 (top panel) and Γ = 1.2 (bottom panel). As b becomes gradually smaller than 1, the cut values, which in Fig. 1 are fixed at 1 (in arbitrary units), do not correspond to the apex energy and the energy level where the rapid decline of the curve occurs. The smaller the Γ value, the more prominent this effect is (see top panel). Note that in 4FGL, 4FGL-DR3 even though function models with spectral shapes equivalent to Eq. (3) are used, no characteristic energy level is mentioned, similarly to cut , except the pivot energy in which the observed spectrum has the smallest error (see equation 4 in 4FGL and Eq. 5 below). The two panels correspond to the indicated Γ, i.e., photon index values. As b decreases cut stops probing the apex energy becoming irrelevant to the energy level where the rapid decrease of the SED begins. The smaller the Γ is the more prominent this effect is. Harding et al. (2018), Torres (2018), Torres et al. (2019), and Harding et al. (2021) have shown that the low-energy part of the Fermi spectrum is mainly synchrocurvature radiation in the synchrotron radiation regime from lower-energy particles while the higherenergy part of the spectrum is synchrocurvature in the CR regime of the high-energy particles. The theoretical scaling of Eq. (2) assumes emission due to CR and thus, the adopted cut in Eqs. (1), (2) should probe the characteristic photon energy scale of the high-energy part of the spectrum.
In that respect, the Γ values and their dependence onĖ reflect the lower-energy particles' properties (i.e., number and energy distribution). On the other hand, the LAT Vela paper, (Abdo et al. 2010a), showed that the observed sub-exponential cutoff of the phaseaveraged spectra possibly reflects the superposition of many purely exponential cutoff spectra corresponding to different phases and magnetosphere regions.
In the phase-averaged spectra, the different phases involve photons from different magnetospheric locations corresponding, in general, to different electric fields, particles energies, and radii of curvature, which introduces ambiguity regarding the best choice for cut that should Figure 2. The Fermi γ-ray light curves together with the exponential cut energies of the corresponding phase-resolved spectra for the Vela, Crab, and Geminga pulsars as indicated in the figure. The black solid and magenta dashed lines denote the photon and energy flux, respectively (left-hand side ordinate). The red lines denote the exponential cut energies (right ordinate). The data have been adapted from DeCesar (2013).
be considered for the γ-ray luminosity function, i.e., FP, relation.
Phase-resolved spectral analysis for a number of bright pulsars (e.g., Vela, Crab, Geminga) have indicated that the cut values corresponding to the phases around the peaks of the γ-ray light curves tend to be higher (DeCesar 2013, see Fig. 2) 3 . The peak phases probe the core-emitting magnetosphere regions considered to be Figure 3. A phase-averaged synthetic spectrum (black line) and the corresponding fit lines considering b-free (green dashed line) and b = 1 (orange dotted line). The red, blue, and yellow lines denote the corresponding phase-resolved contributing spectra of the primary, secondary, and inter/off peak emission components. The big dots and the horizontal lines near the bottom abscissa indicate the corresponding characteristic cutoff energies and ranges (see text for more details). For clarity only 1/4 of the contributing spectra has been plotted.
at the so-called Y-point and in the equatorial current sheet emanating from it (Bai & Spitkovsky 2010;Contopoulos & Kalapotharakos 2010;Kalapotharakos et al. 2014;Cerutti et al. 2016;Philippov & Spitkovsky 2018;Kalapotharakos et al. 2018). In realistic oblique rotators, the axisymmetry breaks and therefore, the Y-point and equatorial current sheet regions are not equally efficient in producing γ-ray emission (Kalapotharakos et al. 2014;Cerutti et al. 2016;Philippov & Spitkovsky 2018;Kalapotharakos et al. 2018). Thus, for any specific object different observers trace, in principle, emission from different magnetosphere regions. Actually, in K19, we claimed that the FP theory implies that many different (though parallel) FPs exist corresponding to these different magnetosphere regions, which, in general, incorporate different particle energies and radii of curvature. These planes contribute, partially to the observed scatter around the FP.
Nonetheless, in typical double-peaked pulse profiles, the cut values corresponding to the different peaks are not always indistinguishable. In phase-averaged spectra, the cut determination corresponding to any of the peaks is nontrivial. Below, we propose treatments of the phase-averaged spectrum that provide energy values that capture the highest exponential cutoff energies found around the Fermi γ-ray light-curve peaks, which supposedly trace the most efficient γ-ray emitting regions for the observer angle values corresponding to the Earth's line of sight.
2.1. Purely exponential cutoff energy We explore here the characteristic energy levels that best probe the maximum cut value, at the peaks of the γ-ray light curves. Unfortunately, even though there are at least a few tens of Fermi pulsars that that are bright enough to allow phase-resolve spectroscopy, currently, there are only half a dozen of bright pulsars with published phase-resolved spectra (Abdo et al. 2010b,c;DeCesar 2013). Thus, in this study, we explore how various characteristic energies correlate with the maximum exponential cutoff energies in the γ-ray light-curve peaks using synthetic spectra.
To explore such behavior, we synthesize mock phaseresolved spectra making considerations that are guided by the phase-resolved spectroscopy results for bright pulsars (DeCesar 2013, see also Figure 2) which then determine the phase-averaged ones.
As mentioned above, the sub-exponential cutoffs can be the result of the superposition of many purely exponential spectra of varying purely exponential cutoff energies, c1 4 , and Γ values that correspond to different magnetosphere regions and particle populations. Thus, one of our simplified assumptions is that the phaseaveraged spectra are the superpositions of many (∼ 100) purely exponential phase-resolved spectra, i.e., Eq. (3) with b = 1, that correspond to different phases. The SED apex energy f The energy where the SED power is f times lower than at A The production of the simulated spectra does not depend so much on the details of the underlying shape of the γ-ray light curves and it only requires a proper (i.e., realistic) consideration of the various spectral components. Thus, we consider that there are three major spectral components, which are those corresponding to two main peaks and one corresponding to the inter-peak and off-peak emission.
We generate 10 4 phase-averaged synthetic spectra in which the percentage of the power 5 corresponding to the inter/off-peak emission (with respect to the total spectrum power) is randomly selected from a uniform distribution that ranges from 10 to 40%. The percentages of the emitting power corresponding to the primary and secondary peaks for each synthetic spectrum are also randomly selected taking into account the constraints: (i) the primary peak carries at least half of the remaining power and (ii) the secondary peak carries a power that is never less than 20% of the total power.
We also assign characteristic (purely exponential) cutoff values, p , s , and i corresponding to the primary, secondary, and inter/off peak emission components, respectively. The maximum characteristic cutoff value among the two peaks X = max( p , s ) can be assigned to any of the two peaks (i.e., primary and secondary) and is selected randomly from a uniform distribution that ranges from 1 to 6 GeV. Note that the primary and secondary peaks are defined based on the corresponding power carried by each peak, i.e., the integrated energy under each peak over a period. We emphasize that this is not necessarily related to the peak amplitudes and the cutoff energies of the corresponding spectra. Observationally, for a majority of Fermi pulsars, the amplitude ratio P2/P1, where P1, P2 are the amplitude of the first and second peak (ordered by the phase lag from the ra-dio peak), increases with the photon energy (see Abdo et al. 2010c). This implies that the highest cutoff spectral energy appears in the second peak (see Figure 2), which however does not guarantee that this peak, i.e., the second one, carries more energy. Thus, in the Crab pulsar, the first peak is the one that carries more energy than the second one despite the fact that the second peak has a higher spectral cutoff energy.
After the X value is selected and assigned to either p or s , the characteristic value corresponding to the other peak is also randomly selected to be between 1/3 and 2/3 of X , which ensures cases where p and s considerably differ and cases where they are close to each other. The phase-resolved spectra indicate the the cutoff energy values corresponding to the inter/off peak emission are either lower than the cutoff values of both peaks or in between these values. Thus, in our model, i is randomly selected to lie in the interval [0.5 min( p , s ), 4/5 X ].
Each of the major emission components runs over a phase interval consisting of spectra with a variety of cutoff energies, i.e., c1 values, (e.g., Fig. 2). Thus, for each emission component, we consider the contribution of, in general, several purely exponential spectra with c1 values that are uniformly distributed over an -interval around the corresponding characteristic values (i.e., p , s , i ). The length of these intervals are randomly selected to be between 1/3 and 2/3 of the corresponding characteristic values, emulating lower and greater variety of spectra within each emission component. We note that the Γ values of the contributing purely exponential spectra are also randomly selected from the interval [0.9, 1.75], which is the interval the vast majority of Fermi pulsar lie. We note that the adopted broad range of Γ values is motivated by the broad range of Γ values seen in the phase-resolved spectra (DeCesar 2013). In any case, it is noted that the results presented below regarding especially the conclusions of the cutoff energy parameters do not depend on the adopted Γ val- Note-The boldface in p, s denotes the highest cutoff energy, i.e., X. ues 6 . We also highlight that monoenergetic e ± produce spectra (from either curvature or synchrotron radiation) with a photon index Γ = 2/3, which is quite different from what is observed in the available phase-resolved spectra (DeCesar 2013). The observed Γ values in the phase-resolved spectra imply the existence of particle populations with a range of smaller energies, i.e., lowenergy pairs, than the high-energy particles that are responsible for the high-energy part of the spectrum and the observed cutoff.
We note that for each component, the number of the purely exponential cutoff spectra, contributing to the synthetic spectrum, depends on the corresponding power percentage. So, for instance for an emission component whose emission power contributes 32% to the total power, we consider 32 contributing (purely exponential cutoff) spectra of equal integrated energy. Superimposing all these spectra, corresponding to the 3 emission components, we construct the corresponding phase-averaged one. Then, we numerically identify the corresponding A value and fit the spectra using the Figure 6. Probability density function histograms for the indicated energy ratios are shown. All the histograms are normalized with respect to the corresponding mean ratio value, which is indicated in each panel. The orange histogram refer to the entire population of the synthetic spectra while the blue ones refer to those with b ≤ 0.75. The gray areas indicate overlapping. function model in Eq. (3) considering both the b free 7 and b = 1 cases. The fitting results also depend on the data extent in energy. Hence, in order to mimic the dif-7 The cutoff energy value in Eq. (3) corresponding to b-free fit will be denoted by cb in the rest of the paper (see also Table 1).

Figure 7.
The 10 value is the energy > A where the SED power is 1/10 of the corresponding maximum power at A.
ferent data extents, all the spectra range from 0.1 GeV up to an energy that is selected randomly and is uniformly distributed from 8 to 30 GeV. In our study, many cutoff energy values, denoted by different symbols, are considered. Thus, in Table 1, we present for clarity all the cutoff energy symbols used in this paper and their description.
In Table 2, we present the phase-resolved spectrum characteristic values for the Vela, Crab, and Geminga pulsars from data adapted from DeCesar (2013). More specifically, we present the powers of the various emission components (i.e., primary, secondary, and inter/off peak) as percentages of the corresponding total powers and the corresponding p , s and i values as these are derived by the phase-resolved spectroscopy (see Figure 2). The properties of these three bright pulsars are commensurate with the choices made for generating synthetic spectra, discussed above.
In Figure 3, we demonstrate the case of one synthetic SED spectrum. The synthetic spectrum denoted by the black line is the result of the addition of 100 phaseresolved spectra. The contributing phase-resolved spectra corresponding to the primary peak, secondary peak, and the inter/off peak emission components are denoted by the red, blue, and yellow color lines, respectively. We note that, for clarity, we plot only one fourth of the contributing spectra, i.e., 25 phase-resolved spectra. The p , s , and i and the corresponding ranges of the c1 contributing to each emission component are also indicated in the figure. The dashed green and dotted orange lines denote the b free and b = 1 fit lines, respectively. As denoted in the figure, the c1 value is superior to cb and A values in probing the X value. In Figure 4, we plot the probability density function histogram of the best-fit values of b for the b-free models for all the 10 4 synthetic spectra. The vast majority of the synthetic spectra have a sub-exponential cutoff (i.e., b < 1) while the distribution peaks at b ≈ 0.85. We note that this maximum is located at higher b values than those that are usually required to fit the observed phaseaveraged spectra of Fermi pulsars (see below the related discussion in Section 3 and Table 3). This implies that the adopted considerations, which have been broadly deduced from a very small sample of available phaseresolved spectra may not accurately describe the underlying statistics of the Fermi pulsar population. Nonetheless, the resulting b values extend down to b ≈ 0.4 with ≈ 25% of all the synthetic spectra having b < 0.75. Thus, below, we examine not only the collective behavior of the entire number of synthetic spectra but also separately the collective behavior of those with b 0.75.
In Figure 5, we plot characteristic synthetic phaseaveraged spectra with X = 3 GeV and for the indicated b values. The black lines correspond to the synthetic spectra, while the dashed green and the dotted orange lines correspond to the b-free and b-fixed (i.e., b = 1) function models, respectively. Similarly to Figure 3, the c1 is superior in probing the corresponding X value. A characteristic SED energy value would probe the behavior of the cutoff energy value that enters the FP relation (see Eq. 2) if the two energy values are proportional to each other. Thus, the substitution of the cutoff energy value in Eqs. (1), (2) by another characteristic SED energy value that is proportional to it would reproduce the same exponent possibly affecting only the normalization factor. In order to see how the various characteristic energy values behave, we calculate for all the synthetic spectra the ratios of these energy values over X . In Figure 6, we plot the probability density histograms of the indicated ratios for the synthetic spectra. We note that in order to make the deviations from the corresponding mean values directly comparable to each other, we have normalized all ratio values to the corre- Figure 9. The 10 values are plotted as a function of the corresponding c1 for all the YPs (black points) and MPs (gray points) for which A exists (see Table 3). The energies sponding mean values, which are also indicated in the panels. The light orange histograms denote the probability densities for all the synthetic spectra while the blue histograms denote only the spectra with b ≤ 0.75. The c1 values are not only these that probe better the actual X but also those that have the smallest dispersion around the mean ratio value. On the other hand, the cb values are those that present the highest spread, having an almost uniform distribution, probing poorly the corresponding X values. The A values seem to behave fairly even though they are inferior to c1 in probing X . More specifically, the mean A / X ratio deviates considerably from 1 while the ratio values present a wider spread around the corresponding mean value. Moreover, A would behave even worse if we had allowed higher Γ values because as Γ goes close to 2 the corresponding A value deviates fast toward very low energy values. On the other hand, the c1 value is not affected by the Γ value. Finally, as we mention in the beginning of this section and as we demonstrate in the next section, the dependence of the spectral photon index onĖ affects A introducing biases making A practically unsuitable for probing the behavior of the cutoff spectral energy in the FP relation.
Another way to probe the high-energy part of the phase-averaged spectrum is the identification of the energy level ( f > A ) where the power of the spectrum reaches a certain fraction, 1/f , of the SED power corresponding to A . We have seen that a fraction of onetenth is adequate because, on the one hand, the corre- Figure 10. The c1 values are plotted as a function of the correspondingĖ for all the YPs (black points) and MPs (gray points). The energy c1 shows an increasing trend withĖ for the lowĖ values, which, however, seems to saturate toward the highĖ values.
sponding 10 is not considerably affected by A and on the other hand, it is not much higher than A and therefore, the flux is still high enough for detection in LAT. In Figure 7, we demonstrate the location of 10 with respect to A in an SED spectrum. In Figure 8, similarly to Figure 6, we plot the probability density histograms of the normalized 10 / X ratios for the synthetic spectra. We see that the 10 / X ratio behaves similarly to the c1 / X ratio and therefore, its use is expected to be able to reproduce the exponent value of the cutoff energy that enters the FP relation (see Eqs. 1 and 2).

FUNDAMENTAL PLANE: FROM 4FGL
TOWARD 3PC In this section, we update the FP of the observed γ-ray pulsars taking into account the data analysis of 4FGL, 4FGL-DR3, data from the ATNF Pulsar Catalog as well as the analysis discussed in the previous section.
The calculation of the best-fit parameters for the FP requires the determination of 4 parameter values (i.e., E, B , cut , L γ ). The 4FGL-DR3 catalog provides the observed energy flux, F E , above 100 MeV, which is easily converted to L γ (i.e., L γ = F E 4πD 2 ) by taking into account the corresponding object distance D, which is taken from the ATNF catalog 8 . This conversion as-sumes, similarly to 2PC, that the beaming factor, f b , is fixed at 1, which implies isotropic emission of a certain level. Furthermore, the ATNF catalog provides thė E value assuming that the stellar moment of inertia is I = 10 45 g cm 2 . We note that for a number of objects, the ATNF catalog provides the ShklovskiiṖ corrections. For these objects, we calculate and use the corresponding correctedĖ values. Finally, the ATNF catalog does not indicateĖ values for 5 4FGL-DR3 objects and therefore, we got theseĖ values from the 2PC catalog.
TheĖ values also lead to the polar stellar surface magnetic field values, B = Ė c 3 P 4 /4π 4 r 6 (1 + sin 2 α), assuming the force-free spin-down power (Spitkovsky 2006;Kalapotharakos & Contopoulos 2009;Pétri 2012;Tchekhovskoy et al. 2013) for a magnetic inclination angle, α = 45 • , and a stellar radius r = 10 6 cm. The 4FGL-DR3 catalog also provides the best-fit parameters of the individual phase-averaged spectra using the function model which for d = 0 takes the form of a pure power law while for d = 0 describes a family of spectral shapes equivalent to that described by the model of Eq. (3), i.e., power law plus sub-exponential cutoff. Comparing Eqs. (3) and (5), it becomes evident that while the corresponding apex energy of the SED reads which is defined if d = 0 and 2b + d − b Γ S > 0. Even though Eqs. (3) and (5) can describe, for d = 0, the same spectral shape families, the form of Eq. (5) and the introduction of parameter d provide more reliable fits in the sense that the best-fit parameters do not correlate strongly especially for the objects with low quality data, i.e., less bright objects (4FGL-DR3).
Bright pulsars indicate a sub-exponential cutoff (i.e., b < 1, usually b 0.75). The function model in Eq. (5) allowed fits with free b even for fainter objects (up to 28). However, freeing b for the rest of the objects was not statistically beneficial and therefore, these objects were fitted with b fixed with a value b = 2/3.
In our study, we use the best-fit parameters provided in the 4FGL-DR3 catalog to reconstruct the individual phase-averaged spectra from 0.1 GeV up to 15 GeV. Then, we fit the reconstructed spectra using the b = 1 model function that provides the c1 values that according to our analysis better probe the high-energy part of the spectra and the corresponding X .
We note that the 4FGL-DR3 catalog contains 255 pulsars. However, the ATNF catalog provides the distance for only 190 of them. Therefore, our analysis incorporates only these 190 pulsars for which the calculation of L γ is feasible. For most pulsars (113 out of 190), the derivation of the distance implements the dispersion measure (DM) using the YMW16 model (Yao et al. 2017). Independent distance estimations are provided for the rest of the pulsars, which are either the only available (for 22 pulsars) or more reliable than the corresponding DM-derived distances (for 55 pulsars). The 22 pulsars without a recorded distance are radio-quiet pulsars since these do not have a DM. The radio-quiet and radio-loud γ-ray pulsars are indeed expected to trace different combinations of α and ζ, i.e., observer angle. Hence, it becomes evident that the sample of 190 pulsars may trace the (α, ζ) values corresponding to the radioloud parameter domain, i.e., small |α − ζ| values, more efficiently because this domain guarantees the derivation of the distance and, therefore, of L γ . The f b value depends, in principle, on ζ, which implies that tracing mainly the ζ values corresponding to those of the radioloud pulsars could bias the L γ values if the f b values of the radio-loud pulsars are systematically different than the average f b values. In any case, the actual f b values remain unknown. In this study, as mentioned above, we assume f b = 1, which may affect only the normalization factor, since it is a population average and not the exponent values of the FP relation.
In Table 3 in the Appendix, we present a set of parameter values, includingĖ, B , c1 , L γ , for the sample of 190 pulsars included in the 4FGL-DR3, which have recorded distances in the ATNF catalog.
The luminosity function, i.e., the FP relation, for the extended sample of the 190 pulsars reads L γ = 10 14.3±1.3 1.39±0.17 c1 B 0.12±0.03 Ė 0.39±0.05 (9) where c1 is measured in MeV, B is measured in G, andĖ and L γ are measured in erg s −1 . The expression in Eq. (9) is consistent with both the FP relation that had taken into account only the 2PC data (see Eq. 1) and the theoretically predicted one (see Eq. 2).
We also derived the FP relation L γ = 10 18.1±2.6 1.17±0.76 c1 B 0.13±0.09 Ė 0.30±0.14 corresponding only to the 22 radio-quiet pulsars, which, despite the large uncertainties that are mainly due to the small number statistics, is compatible with the Eq. 9 implying that the FP relation in Eq. 9 is not drastically affected by this selection bias.
As mentioned above, the 10 values also probe the high-energy parts of the spectra. Thus, the FP relation or luminosity function that incorporates the correspond-  (10) which shows that the exponents are consistent with those in all the previous FP versions supporting the FP underlying theoretical considerations. We also note that the uncertainties in the above equations, as well as those in Eq. 1, are purely statistical, and do not explicitly incorporate systematic uncertainties associated with distances or other variables; however, these systematics may implicitly affect the quoted statistical uncertainty in the luminosity function exponents).
The exponent consistency between the FP relations in Eqs. (9) and (10) implies that the ratio c1 / 10 is practically constant (i.e., proportional relationship between c1 and 10 ). Indeed, as shown in Figure 9, where we plot 10 as a function of c1 , 10 is closely proportional to c1 . We also note that 10 and c1 are not directly related through the same function in the sense that the 10 value is derived from best fit corresponding to the model function in Eq. (5) while the c1 value is the bestfit parameter of the function model in Eq. (3) for b = 1.
In Figure 10, we plot the c1 values as a function of the correspondingĖ values, which shows a weakly increasing trend of c1 withĖ.
In § 2, we discussed that, in general, A is expected to probe relatively well X even though its performance is inferior to those of c1 and 10 . As also mentioned in § 1 and § 2, the greatest issue with A is that it depends, in principle, on the the spectral photon index, Γ (see Eqs. 3 and 6). However, as demonstrated in Figure 11, where we plot Γ as a function ofĖ, the observed Γ values show an increase withĖ (see also 2PC), which affects the behavior of A . The increase of Γ pushes A down (see Eq. 4). We note that the Γ values have a large scatter despite the increasing trend withĖ. Moreover, the apparent scatter is further enhanced by the presence of the high significance sources, i.e., those with b free (denoted by the big dots in Figure 11), which are mainly located below the main diagonal.
In Figure 12, which shows A as a function ofĖ, we see a behavior that actually reflects the behavior of Γ withĖ. More specifically, the A values for the highĖ values, where Γ is higher, are considerably lower than the corresponding c1 values (compared to the situation at lowĖ values). As a consequence A seem stabilized or even decreasing withĖ.
The FP denotes a 3-dimensional (3D) plane embedded in the 4-dimensional (4D) logarithmic space of (Ė, B , c1 , L γ ) and therefore, it is difficult to be visualized. In order to provide a visualization of the FP, we present the FP and the corresponding data in a compacted 3D space. This space involves the variables (X , Y, Z) = (Ė 5/12 B 1/6 , 4/3 c1 , L γ ). In these variables the theoretical FP relation (Eq. 2) reads Z ∝ X Y. Figure 11. The spectral photon index Γ is plotted as a function ofĖ for all the YPs (black points) and MPs (gray points). The Γ values have been calculated using Eq. 6, while the corresponding errors have been calculated considering the error propagation and the values and errors of ΓS, d, b, provided in 4FGL-DR3. The big dots denote the high-significance sources, i.e., b-free sources. The spectral index shows an increasing trend withĖ despite the large spread. Moreover, the MPs seem to have higher Γ values compared to those of YPs, for the sameĖ values.
In Figure 13, we plot the 190 pulsars points in the compacted (X , Y, Z) space seen from two different points of view. The transparent reddish plane corresponds to the theoretical FP (i.e., Z ∝ X Y) while the transparent yellowish plane corresponds to the fitting of the observed points. An edge-on view of the fitting yellowish plane is shown in the right-hand side panel. We note that the scattering around the FP has a standard deviation ∼ 0.33 dex, similar to the one found in K19. Nonetheless, the errors of the fitting parameters are reduced compared to those corresponding to the fitting of the significant smaller sample of objects compiled in K19.

THE FUNDAMENTAL PLANE REVISITED
In K19, one of the main assumptions that led to the theoretical FP relation (i.e., Eq. 2) was that the CR is in the RRL regime. Nonetheless, a more careful investigation shows that the RRL regime requirement is sufficient but not necessary. The RRL regime balances the rate of energy gain to the rate of radiation energy losses. Setting the region of interest near the LC, the Figure 12. The apex energy, A, of the SED spectra is plotted as a function ofĖ for all the YPs (black points) and MPs (gray points). The behavior of A is affected by the behavior of Γ. The plotted errors have been taken from 4FGL-DR3. The increasing Γ toward highĖ affects the behavior of A (compared to that of c1), which especially foṙ E 10 34 erg s −1 show a decreasing trend withĖ. We note that the points with the largest errors that appear mostly for the low A values correspond to pulsars with Γ values close to 2.

RRL regime simply means that
where E BLC is the accelerating electric field in B LC units while γ L and R C denote the Lorentz factor, i.e., particle energy, and radius of curvature, respectively. In K19, we claimed that the γ-ray luminosity of an emitting particle scales as E BLC B LC and therefore as γ 4 L R −2 C , which is true in the RRL regime. However, the emitting power scales always as γ 4 L R −2 C independently of whether this power matches the accelerating one or not. It is the γ 4 L R −2 C scaling that leads to the FP relation (K19), which makes the RRL regime requirement unnecessary. Nonetheless, the RRL regime assumption allows the determination of the accelerating electric field component (i.e., E BLC ), connecting it monotonically to the observed γ-ray efficiency.
In Kalapotharakos et al. (2017, see also K19), assuming the RRL regime and emission near the LC, we derived the accelerating electric field component and its behavior withĖ that revealed a decreasing behavior of E BLC withĖ at highĖ values and a saturation trend toward lowĖ values indicating a fixed dissipation power toward lowĖ and a decreasing dissipation power toward Figure 13. The FP plotted in the reduced dimensinality space (X , Y, Z) = (Ė 5/12 B 1/6 , 4/3 c1 , Lγ). The data points correspond to the YPs (black color) and MPs (gray color), while the reddish plane denotes the theoretical FP (i.e., Z ∝ X Y; see Eq. 2) proposed by K19 and the yellowish plane denotes the the best fit of the observed Fermi pulsars. In the right-hand panel the yellowish plane is shown edge on.
highĖ. Figure 14, which presents E BLC as a function ofĖ is an update of the bottom panel of figure 2 in Kalapotharakos et al. (2017) that had used the sample and data analysis of 2PC. We briefly remind the reader that the expression of the cutoff energy, i.e., c1 = 3c γ 3 L /2R C provides an estimation of the characteristic γ L value for each observed c1 value assuming that R C ≈ R LC . Then taking into account these γ L and R C values and assuming the RRL regime E BLC reads (see also Kalapotharakos et al. 2017).
where E acc is the accelerating electric field component. We also note that in Eq. (12), the radius of curvature has been considered to be equal to R LC 9 .
5. DISCUSSION AND CONCLUSIONS The discovery of the FP relation and the underlying simple theoretical description of it (K19) set strong diagnostics for the modeling of the magnetosphere structure and pulsar high-energy emission. Thus, a major 9 It is noted that Lyutikov et al. (2012)  goal of our study was the exploration of the robustness of the FP relation using a considerably larger data-set that has become available through the recent publication of the 4FGL-DR3 catalog. The motivation behind this exploration was the development of advanced kinetic particle-in-cell models that will be compared to the data evaluating with confidence their ability to reproduce the observed trends (e.g., the FP relation). We have already developed these kind of models and the results of this study will be presented in a forthcoming paper. The FP relation incorporates four parameters, i.e., E, B , cut , L γ . The former two are easily and more or less unambiguously derived by the the observed P,Ṗ values even though their derivation requires some additional assumptions regarding the moment of inertia value(s) and the spin-down power law (e.g., vacuum, force-free).
An accurate derivation of L γ is more problematic not only because of the related flux and distance uncertainties 10 , but also because it depends on the beaming factor, which measures the uniformity of emission on the sky (2PC). The beaming factor varies considerably be- Figure 14. The FP does not require CR emission in the RRL regime. However, assuming CR emission in the RRL regime, the accelerating electric field component, EBLC, in BLC units, is revealed. Here, we plot EBLC as a function oḟ E for all the YPs (black points) and MPs (gray points). The derivation of EBLC assumes that RC ≈ RLC and r = 10 6 cm (see Eq. 12).
tween different models while it may also vary within the same model depending, in principle, on the observer angle. However, the consideration of a different average beaming-factor value would collectively affect all the L γ values leaving unaffected the corresponding scaling, i.e., the exponents in the FP relation. Any L γ ambiguity that is the result of the variation of beaming factors around the mean value contributes just to the spread around the FP.
On the other hand, the derivation of the characteristic cut value is the most intriguing one for a variety of reasons. The exponent of cut in the FP relation is larger than all the other exponents, which implies that the luminosity function is more sensitive to cut . Moreover, the observations indicate that the spectral characteristic cut values have a rather limited range that barely exceeds one order of magnitude. Finally, the observed phase-averaged spectral shapes and the function models that seem to fit them provide a variety of characteristic energies whose exact physical meaning and their possible relation to the synchro-curvature emission, near the curvature radiation regime, is not always clear. The phase-averaged spectra consist of collections of emission from different magnetosphere regions reflecting, in general, different physical properties (e.g., accelerating electric field components, particle energies, radii of curvature of particle trajectories). The phase-resolved spectra corresponding to different phases may also trace different magnetosphere regions although, in this case, the variety of the physical properties should be significantly reduced compared to that in the phase-averaged spectra.
The most relevant cut values that should enter the FP relation are the ones corresponding to the peaks of the γ-ray light curves, which trace the γ-ray emis-sion corresponding to the core of the dissipative regions, i.e., the equatorial current sheet. However, despite the existence of tens of bright pulsars, only a handful of them have published phase-resolved spectra, which makes the phase-averaged spectra currently the only available source of analysis.
Single e ± radiating in single region of space produce spectra with exponential cutoffs, while the phaseaveraged spectra are better fitted by sub-exponential cutoffs. The contribution of different magnetosphere regions is able to explain the significant deviation from the pure exponential cutoff. Thus, we decided to use a toy-model approach to explore the phase-averaged spectra behavior assuming that these are simply superpositions of pure exponential cutoff spectra (i.e., the phaseresolved ones). Our considerations for constructing synthetic spectra were based on the properties of the few available phase-resolved spectra. However, these considerations taken from only a few cases were quite broad and therefore, they were probably not able to correctly describe the spectral shape statistics of the entire Fermi population. Nonetheless, our goal was not to reproduce the spectral statistics of the observed objects, but to explore how well various energy parameters probe the highest characteristic cut energy value between the two peaks, i.e., the X value that should be considered in the FP relation.
Thus, we examined different spectrum function models, focusing on the determination of the capacity of the various spectral photon energy parameters for probing X . Using the synthetic phase-averaged spectra, we found that indeed these have sub-exponential cutoffs.
However, we showed that the cutoff energy corresponding to the pure exponential cutoff function model is superior in probing X . Thus, our analysis implies that even though the sub-exponential model function may better describe the totality of the spectra, the cutoff energy that is more relevant to the cut in the FP relation is better captured by the one corresponding to the pure exponential model function, i.e., c1 .
We note that the shape of most of the synthetic spectra showed a sub-exponential cutoff corresponding to a b exponent (see Eq. 3) that was higher than those indicated by the most bright pulsars even though the fitting of several synthetic spectra provided quite low b values. The existence of the relatively high b values implies that either the spectral rules (e.g., relative powers and characteristic cutoff energies of the main emission components) adopted in our study do not always describe realistic spectra or the individual spectra that contribute to the synthetic phase-averaged ones are not purely exponential. As discussed above, the phase-resolved spectra may also probe different magnetosphere regions that contribute to the same phase (Bai & Spitkovsky 2010;Contopoulos & Kalapotharakos 2010;Kalapotharakos et al. 2014). Nonetheless, the variety of the physical properties (e.g., γ L , accelerating field components, R C ) of the regions contributing only to one phase is expected to be smaller than those in the phase-averaged spectra. This practically means that any deviation of the phase-resolved spectra from the pure exponential cutoff should be smaller (i.e., b closer to 1) than what is observed in the phase-averaged spectra. A contribution of phase-resolved spectra with smaller (than 1) b values would result in phase-averaged spectra with even smaller b values. Finally, there is a possibility that the b value at least for some pulsars is affected by the existence of an additional very high energy emission component (inverse Compton, IC; see Harding et al. 2021), which may form a softer transition between the two spectrum components (i.e., CR and IC).
In any case, a detailed description of the variety of the phase-resolved spectra will remain unclear until we get access to a larger amount of high-quality data. Moreover, any, so-far, available phase-resolved spectra have adopted the pure exponential model function, i.e., b = 1 (see DeCesar 2013). Hopefully, this work will further motivate observers to perform a detailed phase-resolved spectroscopy incorporating, if possible, more objects.
Taking into account our conclusions regarding the interpretation of c1 , we used the results of the spectral analysis in 4FGL-DR3 to reconstruct the phaseaveraged spectra for 190 pulsars, which allowed the calculation of the corresponding c1 values. Taking also into account the 4FGL-DR3 γ-ray fluxes as well as the distances and the spin-down powers from the ATNF catalog, we updated the observed FP relation by increasing more than two fold the number of Fermi pulsars compared to those we used in our original study K19. The extended sample showed results in agreement with those corresponding to the sample of 2PC indicating that the high-energy part of the Fermi spectra is consistent with CR.
Our study provides a preamble for the anticipated 3PC focusing mainly on the exploitation and further interpretation of the 4FGL-DR3 data that allow the derivation of meaningful parameter values that enter the FP relation. Nonetheless, the 3PC is expected to provide advanced data treatments and possibly improved fitting function models expanding also considerably the number of objects compiled in this study.
Finally, we revisited the FP theoretical background by relaxing the requirements that lead to the derivation of the theoretical FP relation. More specifically, we showed that the RRL regime assumption is redundant for the FP relation; however, it sets additional constraints that allow the determination of the underlying accelerating electric field component.

ACKNOWLEDGMENTS
We acknowledge helpful discussions with Isabelle Grenier and David Smith. We would like to thank the LAT group for providing valuable comments and suggestions, which helped us to improve the quality of the paper. We would also like to thank the International Space Science Institute (ISSI) for providing financial support for the organization of the meeting of the ISSI Team that was led by I. Contopoulos and D. Kazanas. C.K. and Z.W. are supported by the Fermi Guest Investigator program and the NASA Theory Program. D. K. is supported by the NASA Data Analysis Program. The material is based upon work supported by NASA under award number 80GSFC21M0002. Table 3. A set of parameters for 190 4FGL-DR3 young and millisecond pulsars.  Table 3 continued  Table 3 continued  Table 3 continued  Table 3 continued