Frequency Evolution Behavior of Pulse Profile of PSR B1737+13 with the Inverse Compton Scattering Model

The radio radiation mechanism is one of the open questions in pulsar physics. Multiband observations are very important for constraining the pulsar radiation mechanism. In this paper, we investigate the pulse profiles of PSR B1737+13 and its evolution with the frequency. The integrated pulse profiles are obtained from the European Pulsar Network and the Australia Telescope National Facility data, together with recent observations from the largest dish Five-hundred-meter Aperture Spherical radio Telescope. The radiation components are separated with the squared hyperbolic secant functions, and the radiation altitudes of each radiation component at different frequencies are calculated. It is found that the radio radiation at different frequencies comes from different altitudes. The frequency evolutions of separations for the inner and outer cone components are studied. It is found that the separations of the inner and outer cone components have opposite frequency dependence. We simulate the RFM of PSR B1737+13 with the inverse Compton scattering (ICS) model and find that the RFM can be naturally described by the ICS model. Through the simulation, the radio radiation region of PSR B1737+13 is determined, and the result shows that the radio radiation of this pulsar may be generated in the annular gap region.


Introduction
Since pulsars were discovered more than 50 yr ago (Hewish et al. 1968), observations and theoretical studies of pulsars have achieved great progress. To date, thousands of pulsars have been discovered (Manchester et al. 2005;Han et al. 2021;Pan et al. 2021). Many theories have been proposed to understand their radiation mechanism (e.g., Sturrock 1971;Ruderman & Sutherland 1975;Qiao & Lin 1998;Gil & Sendyk 2000). Goldreich & Julian (1969) first proposed that a neutron star is surrounded a charge-separated magnetosphere. Ruderman & Sutherland (1975) suggested that there is a core gap (CG) accelerating region (the region surrounded by the critical magnetic field lines) above the polar cap of the pulsars, and high-energy particles accelerated in this region will flow outward along the lines of magnetic to generate observable radio radiation. Based on this model, some works have been made to study the pulsar radiation physics and great success has been achieved (e.g., Michel 1982;Gil et al. 2003). Qiao et al. (2004) proposed the annular region (AG), a region above the pole cap surrounded by the critical magnetic field line and the last open magnetic field line that can also accelerate high-energy particles (Du et al. 2011(Du et al. , 2012Qiao et al. 2013;Shang et al. 2020). They proposed an inverse Compton scattering (ICS) model and suggested that the radio emission of the pulsars is the result of the secondary particles scattering off "low-frequency" waves (Zhang & Qiao 1996;Qiao & Lin 1998;Qiao et al. 2001). They further studied the pulsar radiation mechanism emission and found that the ICS process may be the radiation mechanism of pulsars (Qiao 1988). The ICS model has been widely used to explain various observational phenomena, for example, the polarization properties, frequency evolution behavior, and mode changing of pulsars (Zhang et al. 1997b;Qiao et al. 1999;Xu et al. 2000;Qiao et al. 2001;Zhang et al. 2007;Shang et al. 2017Shang et al. , 2021. Despite the success, so far, there are still gaps between observations and theories, and a unified theory of the radiation mechanism of pulsars is still an unsolved mystery. Multiband observations are very important for understanding the emission physics of pulsars, especially the frequencydependent behaviors of pulse profiles and their evolution, and are very useful for testing and constraining pulsar radiation models. The pulse profile variations with respect to observing frequencies can provide very valuable information, such as the region, beam shape, and spectrum behavior of pulsar emission. Based on multiband observations, some studies on the radiation region, beam shape, spectrum, and width of pulsars have been conducted to clarify the mechanism and radiation geometry of pulsar radio radiation (e.g., Rankin 1983Rankin , 1983Lyne & Manchester 1988;Qiao et al. 2004;Pilia et al. 2016;Jankowski et al. 2018;Zhao et al. 2019;Shang et al. 2021). For the beam shape of pulsars, some models have been proposed to explain the various evolutionary behaviors of the pulse profiles, such as the "core-cone" beam, the "core-double cone" beam, the fan beam, the microbeam (Backer 1976;Rankin 1983;Lyne & Manchester 1988;Wang et al. 2014;Dyks et al. 2016), and so on. These models have their own advantages in explaining some observed characteristics; however, consensus on the radiation mechanism of the pulsars has not been fully achieved yet. The ICS model considers that the radiation beam may have two or three radiation components (Qiao & Lin 1998;Qiao et al. 2001), and it has been widely and successfully used to study the profile evolution with the frequency for various types of pulsars (Zhang et al. 2007;Shang et al. 2017Shang et al. , 2021; these types have mainly triple, core-single (Rankin 1983), and conal triple profiles (Rankin 1993b). The study of multiple-profile pulsars with the ICS mode is lacking and is of great importance to the model itself and the physics of pulsars.
In recent years, the observational data have been increasing; in particular, large numbers of high-frequency and lowfrequency data have been accumulated, which brings opportunities and challenges to the ICS model. B1737+13 is a multiple-profile pulsar (Rankin 1990), which has been observed by the Five-hundred-meter Aperture Spherical radio Telescope (FAST) very recently. In this paper, in the framework of the ICS model, we will study the evolution behavior of the pulse profile of PSR B1737+13 with the frequency with the aim to test the ICS model and to constrain the radiation physics of pulsars. This paper is structured as follows. In Section 2, the information of the collected data and the reduction methods are given. In Section 3, we separate the components of the pulse profile and give the frequency dependence of the different radiation components. In Section 4, the ICS model simulations of evolution behaviors for this pulsar are presented. The radiation region is calculated and given in this section. The discussions and conclusions are given in Section 5.

Data Information and Reduction Methods
In this work, the 19-beam receiver of FAST was used to observe PSR B1737+13 on 2018 August 30 (MJD 58360). The data were sampled at 49.152 s intervals. The observed frequency range was from 1000 and 1500 MHz with 4096 channels of 0.122 MHz. The data were recorded in the search mode PSRFITS file (Hotan et al. 2004). To analyze the data set, each observation was folded to a phase period by using the DSPSR (Hotan et al. 2004;van Straten & Bailes 2011) and the timing ephemeris provided by the pulsar catalog (PSRCAT; Manchester et al. 2005). In the observations, the signal from the target will be affected by radio-frequency interference (RFI). We first automatically flagged the RFI and removed 5% of the band edges using the PSRCHIVE software package (Hotan et al. 2004;van Straten et al. 2012). All RFI were removed manually by using the interactive software command pazi provided by PSRCHIVE. To study the evolution of the pulse with the frequency, each observation with a bandwidth of 400 MHz was divided into four subbands; the bandwidth of each subband was 100 MHz. The center frequencies of these subbands were 1100, 1200, 1300, and 1400 MHz, respectively. Finally, we obtained the pulse profile of each subband data according to the standard extraction process of the pulse profile of the pulsar. In addition to the observations of FAST, we also collected other multiband pulse profiles from the European Pulsar Network (EPN) 7 and the Australia Telescope National Facility (ATNF). 8 In this paper, the pulse profiles at center frequencies of 102.75, 143, 168, 408, 610, 910, 1369, 1420, 1642, and 4850 MHz (Seiradakis et al. 1995von Hoensbroech & Xilouris 1997;Gould & Lyne 1998;Kassim & Lazio 1999;Pilia et al. 2016;Johnston & Kerr 2018) were obtained from the EPN, while the pulse profiles at a frequency of 3100 MHz were obtained by reprocessing the archive data of the ATNF.

Component Separation and Frequency Evolution
Behavior of the Pulse Profile The beam shape and the acceleration regions are very important for constraining the pulsar emission model. Based on the "core-cone" model (Backer 1976) and "core-double cone" model (Rankin 1983(Rankin , 1993b(Rankin , 1993a, the emission beam contains two or three emission components, i.e., a center cone and one or two nested cones. The observed pulse profiles are the section of the emission beam sweeps across the line of sight. The core component and the two cone components of the emission beam correspond to the central peak component and the wing components of the pulse profiles, respectively. Usually, each component is assumed to have a Gaussian shape (Kramer 1994;Wu et al. 1998;Zhang et al. 2007). Lu et al. (2016) found that the square hyperbolic secant function can better fit the leading wings of the pulse profile compared with the Gaussian function. Following Lu et al. (2016), we assume that each component has a square hyperbolic secant shape. Therefore, the fitting function of the integrated profile can be written as follows: where I i , p i , and w i are the intensity, position of peak, and FWHM of the ith peak component, respectively.
As shown in Figure 1, when the frequency increases, the pulse profile of B1737+13 evolves from a triple profile to a profile with no obvious number of peaks. Compared to pulsars with more stable shapes and peak numbers, such as B0525+21, the pulse profile of B1737+13 is more complicated. According to the morphological classification of the pulse profile proposed by Rankin (1983), it is a classic multiple-profile pulsar (Rankin 1990(Rankin , 1993bOlszanski et al. 2019). Many previous studies have shown that this pulsar has five components (Hankins & Rickett 1986;Force & Rankin 2010;Skrzypczak et al. 2018). In this work, we use five square hyperbolic secant functions to separate the pulse profiles at each frequency, and the fitting results are shown in Figure 2. The fitting parameters for each frequency are listed in Table 1. The first and fifth peaks correspond to the outer cone radiation component, the second and fourth peaks correspond to the inner cone radiation component, and the middle peak corresponds to the core radiation component. Using the parameters in Table 1, we calculate the separation Δf of the inner cone and outer cone radiation components at each frequency and fit the frequency evolution of Δf with the power-law function: where A is a coefficient, ν is the separation index of the components, and ω is the observing frequency.
As shown in Figure 3, the separations of these two cone components have an opposite frequency dependence. The component separation index of the outer cone is ν = −0.13, and the component separation index of the inner cone is ν = 0.15. This phenomenon, where the separation of the inner cone component increases with increasing frequency, deviates from the radius-frequency mapping (RFM). Although the RFM model can explain the increase in the frequency of the component separation, some recent studies have found that the frequency dependence of the pulsar width does not always conform to this relationship (Chen & Wang 2014;Pilia et al. 2016;Zhao et al. 2019;Xu et al. 2021). This presents a challenge to this model.

Profile Evolution, Radiation Altitude, and Radiation
Region of PSR B1737+13

The Inverse Compton Scattering Model
The ICS model was first proposed by Qiao (1988). The basic picture of this model is as follows. It is assumed that the periodic breakdown of the inner gap near the polar cap region will generate low-frequency electromagnetic waves (Ruderman & Sutherland 1975). These low-frequency photons propagate outward to a finite height of interest and are scattered via inverse Compton scattering by the high-energy secondary particles moving along the magnetic field lines. The upward-scattered radio photons will allow us to observe radio emission. For this process, the frequency w¢ of the radio emission can be determined by the following equation: where γ is the Lorentz factor, ω 0 is the initial frequency of the photons generated by gap sparking, assumed to be 10 6 Hz (Qiao & Lin 1998;Qiao et al. 2001;Lu et al. 2016), and h = » u 1 c . θ i is the angle between the direction of the photons moving from the sparking point and the direction of motion of a particle.
where f is the azimuth angle of the magnetic field line, f c is the azimuth angle of the spark point, θ c is the polar cap angle of the spark point (Qiao & Lin 1998), θ is the polar angle of the emission point, R is the radius of pulsar, and r is the radiation altitude. The radiation altitude in the dipole field can be written as: where R e is the maximum radius of this magnetic field line and R e = λR lc (λ 1). R lc = Pc/2π is the radius of the light cylinder, and P is the pulse period of pulsar. The polar angle θ and beam angle θ μ are related by where α is the magnetic inclination angle, and β is the impact angle. Δf is the separation of the emission components, which can be obtained by observation. α and β are provided by Lyne & Manchester (1988) and Rankin (1993b). The results of  Lyne & Manchester (1988) show that α and β are 40°.8 and 2°. 5, respectively, but Rankin (1993b) provides values of α and β of 41°. and 1°. 9, respectively. In the calculations below, we use the results of the latest year, that is, the inclination angle and the impact angle are α = 41°and β = 1°.9. Of course, different α and β will only change the value of the simulation parameter, so this will not affect our results (Zhang et al. 2007).

Energy Loss of High-energy Particles
Since the secondary particles usually slow down due to various energy-loss mechanisms (Zhang et al. 1997a), the Lorentz factor of high-energy particles will gradually decrease when the highenergy particles flow along the magnetic field lines. Qiao & Lin (1998) pointed out that the energy loss of high-energy particles can be assumed to be: where γ 0 is the initial Lorentz factor, and k is a factor related to the energy loss of high-energy particles. Some studies assume that the Lorentz factor of high-energy particles varies is an exponentially manner (Qiao et al. 2001;Zhang et al. 2007;Shang et al. 2017Shang et al. , 2021: , 9 0 ( ) but the above equation ignores an energy term (Wang et al. 2013;Lu et al. 2016). Therefore, the change in the Lorentz factor of high-energy particles with the radiation altitude should be written as: where γ c is the corresponding value of the final secondary particles before and after the pair cascade. Based on Equations (5)-(7), we calculate the radiation height at each frequency of the outer cone radiation component and combine Equation (3) to obtain the value of the Lorentz factor for each radiation altitude. Here, we assume that the radiation comes from a field line of 0.98 times the radius of the polar cap (the reason for using this value is explained in Section 4.3). Based on the radiation altitude and Lorentz factor of the outer cone component at different frequencies, we compare Equations (8)-(10) to determine which function is more suitable for describing the change in the energy of the high-energy particles. In order to facilitate the fitting, we simplify Equations (8)-(10) into the following form: The fitting results are shown in Figure 4, and the fitting parameters are listed in Table 2. Among them, Equation (11) indicates that the energy of high-energy particles varies in a linear form, shown as a straight line in Figure 4, which obviously does not coincide with the data points. So we only give the fitting results of the black circles. As shown in Figure 4, Equation (13) (which is similar to Equation (10)) can better describe the decreasing trend of the Lorentz factor of high-energy particles. In the following simulations and calculations, we assume that the Lorentz factor of the high-energy particle varies with the radiation height in the form of Equation (10). Interestingly, if Equation (10) is the energy-loss function of high-energy particles, then the Lorentz factor will have a lower limit (∼210; see Table 2). This implies that the Lorentz factor of high-energy particles will not decrease after a certain altitude. This may be caused by the acceleration of high-energy particles in the flow process (Qiao et al. 2002). On the one hand, the energy of high-energy particles is continuously reduced due to the ICS process, and on the other hand, it is increased for some reason. When they are equal, the energy of high-energy particles tends to be constant. This requires further research.

Simulation of Frequency Behaviors of Pulse Profiles in ICS Model
The ICS process means that the low-frequency photons are scattered via inverse Compton scattering by the high-energy secondary particles moving along the magnetic field lines, thereby radiating the electromagnetic waves we observe. In the ICS model, for a certain open magnetic field line, different altitudes will have different polar angles and beam angles. Through Equation (3), we can obtain the dependence of the radiation altitude and frequency, for which the Lorentz factor of highenergy particles decays in the form of Equation (10). Similarly, we can obtain a relationship between the frequency and beam angle. Since the range of the azimuthal angle and the spark point are from 0 to 2π, we can determine the range of the beam angle and emission altitude at a certain frequency. Figure 5 shows our simulation results for B1737+13; panels (a) and (b) represent the evolution of the beam angle and radiation altitude with the frequency, respectively. The area enclosed by the red line and the blue line represents the range of the beam angle and radiation height at each frequency. The simulated curve in Figure 5 shows a drop-rise-drop pattern, corresponding to the core, inner cone, and outer cone components, respectively. It also reflects the frequency dependence of different components. The data points in the figure are obtained from Equations (5)-(7) with the fitting parameters in Table 1. Through the ICS simulation, the values of the initial Lorenz factor, the final secondary particle, and the energy-loss factors of particles are γ 0 = 1.9 × 10 4 , γ c = 210, and k = 0.195, respectively. The feet of the field lines and the locations of the sparking points above the polar cap are located at 0.98 R p and 0.93 R p , respectively, where R p = R 1.5 Ω 0.5 c −0.5 is the polar cap radius defined from the last open field line (Ruderman & Sutherland 1975). Ω is the rotation angular velocity.
Our results show that the ICS model can well describe the RFM of PSR B1737+13, as shown in Figure 5. The inner cone components gradually separate from the increase in the frequency and cannot be well explained by the traditional RFM model (Mitra & Rankin 2002;Shang et al. 2017;Zhao et al. 2019). Through calculations, we find that the inner cone radiation component mainly comes from an altitude of tens of kilometers. What is interesting is that its high-frequency radiation has a lower radiation height than low-frequency radiation. The evolution of the separation of the radiation component of the outer cone with the frequency is opposite to that of the inner cone. The outer cone radiation altitude is about several hundred kilometers, and high-frequency radiation is generated at low altitudes (∼120 km), while low-frequency radiation is generated at high altitudes (∼250 km). Compared with the outer cone component, the inner cone component has a lower radiation altitude, which is consistent with the results of some previous studies (Gupta & Gangadhara 2003;Srostlik & Rankin 2005;Zhang et al. 2007;Shang et al. 2017). Rankin (1990) pointed that the core radiation is thought to originate very close to the pulsar surface. As shown in Figure 5, our results confirm this view and indicate that the radiation altitude of the pulsar core is ∼20 km.

Radiation Region
In the dipole magnetic field model, the shape of the annular gap is a function of the inclination angle α (Qiao et al. 2004;Wang et al. 2006), and both AG and CG may obtain a potential drop to accelerate particles to produce observable radio radiation (Qiao et al. 2004;Lu et al. 2016;Shang et al. 2017Shang et al. , 2021. With the simulation results, we try to determine the location of the spark above the pole cap and the magnetic field lines where the radiation comes from. A section of the polar cap of this pulsar is shown in Figure 6. The black solid and blue dashed lines are the last open field and the critical field line, respectively. The section between the blue dashed and black solid lines is the AG region. The region surrounded by the blue dashed line is the CG region. The green solid line denotes the feet of the field line used for calculation in this simulation. The red solid line is the section of the spark point above the polar cap, and the red star dots A and ¢ A are the radiation points used in this simulation. The ICS simulation shows that the observed radio radiation of PSR B1737+13 may be generated in the AG region.

Conclusions and Discussions
We investigate the RFM of PSR B1737+13 by using the observations of FAST and the data archive obtained from the EPN and ATNF databases. Using the squared hyperbolic secant functions, we separate the radiation components of this pulsar and calculate the radiation altitude of each radiation component at different frequencies. Then, the frequency evolution of the   separation of the inner and outer cone components are studied. Finally, the RFM is simulated with the ICS model. Based on the simulated parameters, the radio radiation region is determined. We measure the separations between the inner cone components and the separations between the outer components, respectively. It is found that the separation evolutions with the frequency for the inner and outer cones are different. We fit the separation evolutions of the inner and outer cones with a powerlaw functions. It is found that the separation index of the inner and outer cones are −0.13 and 0.15, respectively. Mitra & Rankin (2002) presented a fitting result of the separation evolution for 10 pulsars with a function introduced by Thorsett (1991) and found that the constant separation appears to reflect the inner cone radiation. Some studies have found that the separation of the inner cone components does not approach zero (Zhang et al. 2007;Shang et al. 2017). Qiao & Lin (1998) suggested that this frequency behavior is the result of different frequency-mappings in the inner and outer cones, which need to be verified by a larger sample and a wider frequency of high-quality observations.
The Lorentz factor of high-energy particles flowing along the field lines will decrease due to various mechanisms. We have compared three different forms of energy-loss functions. The results showed Equation (10) can better describe the process of the Lorentz factor of high-energy particles decreasing with the increase in the radiation altitude. It is worth noting that the Lorentz factor of the high-energy particles of this pulsar decreases rapidly with the increase in altitude and gradually becomes flat after exceeding a certain altitude. However, it will not drop to zero, which is consistent with the conclusions of Lu et al. (2016). This implies that the Lorentz factor of highenergy particles will not decrease at higher altitudes. More observations of low frequency and high signal-to-noise ratio are needed to verify this result.
We have used the ICS model to simulate the evolution behavior of the beam angle and radiation altitude of the two cone components of PSR B1737+13 with the frequency (see Figure 5). The ICS model can naturally give the variation of the component separation and radiation altitude of the inner and outer cone components of the pulsar with the frequency. It can be clearly seen from Figure 5 that the radiation of the same frequency can come from three different altitudes. The core cone component is thought to originate very close to the pulsar surface (Rankin 1990). For the outer cone component, its highfrequency radiation comes from a low altitude, and its lowfrequency radiation comes from a high altitude. However, the inner cone radiation is just the opposite. Its high-frequency radiation originates from a higher altitude than that of the lowfrequency radiation. At the same time, the outer cone component has a higher radiation altitude than the inner cone component (Gupta & Gangadhara 2003;Srostlik & Rankin 2005). According to the inclination angle α (Rankin 1993b) and the ICS simulation results of PSR B1737+13, we determined the pulse radiation region, and the results show that the radiation of this pulsar can originate from the AG region.
Through ICS process simulation, the initial Lorentz factor γ 0 is found to be 1.9 × 10 4 . For PSR B1737+13, the value of the Lorentz factor is larger than that (γ ∼ 800) of secondary particles and lower than the energy of the primary particles γ ∼ 10 6 given by Ruderman & Sutherland (1975). Since not all particle pairs are produced at the bottom of the gap above the polar cap, not all primary particles can obtain such a large Lorentz factor (Zhang et al. 2007). The Lorentz factor does not exceed this value, which means that the value used in this article is reasonable.