Compatibility of the Observed Rotation Parameters of Radio Pulsars on Long Time Scales ()

Arkady Efimovich Avramenko^{}, Boris Yakovlevich Losovsky^{}, Vyacheslav Dmitrievich Pugachev^{}, Tatyana Viktorovna Shabanova^{}

Pushchino Radio Astronomy Observatory, P. N. Lebedev Physical Institute, Russian Academy of Sciences, Pushchino, Russia.

**DOI: **10.4236/ijaa.2018.81003
PDF HTML XML
787
Downloads
1,346
Views
Citations

Pushchino Radio Astronomy Observatory, P. N. Lebedev Physical Institute, Russian Academy of Sciences, Pushchino, Russia.

We present the results of an analysis of the rotation regularities of pulsars from the current observations and retrospective data on the timing of pulsars at Pushchino Observatory, obtained since 1978 to 2017. Parametric approximation of numerical series of Times of Arrival (ToAs) the convergence polynomial series defined the combinations of numerical values, unique for each pulsar, identical in any reference system for coinciding epochs, irrespective of duration of observations that corresponds to the coherence of periodic radiation of pulsars. The braking index at all observed pulsars is n =−(0.9 ± 0.2) that corresponds monotony of slowing down of rotation (>0). The polynomial power series of the intervals calculated in observed parameters of rotation represent an analytical time scale on which variations of observed intervals, uncorrelated with rotation, are counted, as at pulsars B1919+21 and B0809+74, or the unpredictable variations of intervals defining the movement of a pulsar in the supernova remnants, as at PSR B0531+21 or at PSR B1822-09 which spontaneous movement also doesn’t exclude its belonging to the supernova remnants.

Share and Cite:

Avramenko, A. , Losovsky, B. , Pugachev, V. and Shabanova, T. (2018) Compatibility of the Observed Rotation Parameters of Radio Pulsars on Long Time Scales. *International Journal of Astronomy and Astrophysics*, **8**, 24-45. doi: 10.4236/ijaa.2018.81003.

1. Introduction

The most significant results in studying the physical properties of pulsars are connected with periodic regularity of a radio emission of the rotating magnetized neutron star. The energy losses accompanying star radiation lead to the gradual slowing down of its rotation, so the numerical value of the period of rotation specified in the catalog is fixed together with the relating observation epoch.

The rotation period, being observed parameter, represents the measured value, which is determined directly by the moments of arrival of impulses of radiation on the radio telescope. Already from the beginning of the 70th years, near only several years after opening and fast addition of the list of new pulsars, the numerical values of the period and its derivative were accurately measured for all the pulsars included in the first fundamental catalog of pulsars [1] , which have not lost meaning so far.

Many of results are obtained by the pulsar timing based on the precise measuring the topocentric Times-of-Arrival (TоAs) of the pulsar signals at the radio telescopes, denominated in the UTC time scale [2] . The process of pulsar timing is fundamentally dependent on an accurate description of everything that affects the ToAs of the pulsed radiation at the telescope. In addition to a local time standard and the Solar-System ephemerides pulsar timing requires knowledge of the pulsar’s pulse period and spin-down, its astrometric position and proper motion, its distance, parallax (where detectable) dispersion measure in the interstellar medium and (on the binary system) any orbital parameters. All of these parameters are included in a timing model, which can be used to predict the phase of the pulsar’s periodic signal at any point in time. The differences between the measured ToAs and those predicted by the timing are the “timing residuals”, which are the unmodelled difference between the observations and the theory [3] .

The continuing research of pulsars, especially in the process of accumulation long-term, within several decades, archives of timing, showed that braking of pulsars is not limited only to the first derivative. An analysis of timing noise has shown that these deviations are due to a random, continuous wandering of the pulse phase that produces long-term cubic polynomial components in the timing residuals, which are correlated with the period derivative and are quantified by a parameter based on the non-zero second derivative. The timing residuals cannot be explained by the models that are based on noise-like processes because such models are not consistent with observations over long time scales.

Crucial is an exact definition of the second derivative which contribution usually doesn’t exceed several microseconds within 2 - 3 years of observations that is several orders lower not than a threshold of detection against the background of unmodelled timing noise. As Shabanova notes in the generalizing paper [4] : “For ordinary pulsars, the values of P2dot due to slowdown are very small in comparison to the measurement uncertainties. In actuality, the long-term cubic polynomial components in the timing residuals of most ordinary pulsars lead to the large observed P2dot values and accordingly to anomalously high values of the braking indices (Cordes & Helfand 1980; Gullahorn & Rankin 1982; Cordes & Downs 1985)”.

Abnormally high, not giving in to a physical explanation of residuals and second derivative variations are also discussed in the extensive publication of Hobbs, Lyne and Kramer [5] : “For a large number of pulsars, the timing residuals deviate by more than one pulse period…The observed ν2dot values for the majority of pulsars are not caused by magnetic dipole radiation or by any other systematic loss of rotational energy, but are dominated by the amount of timing noise present in the residuals and the data span”.

Therefore, the accepted timing models quite often lead to considerable distortions in estimates of physical properties of the pulsars, which are strongly differing for the same objects on different sources. Unpredictable, abnormally high variations of residuals are interpreted as sudden failures of the period of rotation (glitches). In this consideration, especially significant is a reliable extraction of physical characteristics of rotation of pulsars, including the second derivative, on directly observed data of, in view of basic restrictions of reliability of statistical estimates on residuals.

Our approach, in general, is to find the analytical relation of the pulsar time intervals and the physical parameters so that the numerical values of these parameters should be determined and best matched with measured values of the observed intervals. The comprehensive fitting should be enough only those parameters whose exact values can be derived directly from observation.

2. Analytical Model of Pulsar Timing

The main point of the model is to convert the numerical series of the observed ToAs, measured with respect to the some initial event of pulsar radiation, into the analytical form of intervals, calculated from the observed rotation parameters that most closely correspond to the numerical series of the ToAs [6] .

Pulsar time (PT) intervals of any pulsed periodic radiation event N express by the function of its rotation parameters ${P}_{0},\stackrel{\dot{}}{P},\stackrel{\xa8}{P}$ in the time domain:

$PT\left({P}_{0},\stackrel{\dot{}}{P},\stackrel{\xa8}{P}\right)={P}_{0}N+\frac{1}{2}{P}_{0}\stackrel{\dot{}}{P}{N}^{2}+\frac{1}{6}\left({P}_{0}^{2}\stackrel{\xa8}{P}-2{P}_{0}{\stackrel{\dot{}}{P}}^{2}\right){N}^{3},N=1,2,3,\cdots $ (1)

The right part of the Equation (1) represents a power series expansion Maclaurin of intervals PT, counted from the initial emitted pulsar event. The structure of intervals PT contains linear, square and cubic components, which are defined by fixed rotation period P_{0} on the initial epoch and its derivatives
$\stackrel{\dot{}}{P},\stackrel{\xa8}{P}$ .

The power series on the right side (1)is limited to only 3 parametric components,sufficient for its convergence.Indeed,for a typical pulsar(
${P}_{0}\approx 1s,\stackrel{\dot{}}{P}\approx {10}^{-15}\text{-}{10}^{-16}s\cdot {s}^{-1},\stackrel{\xa8}{P}\approx {10}^{-28}\text{-}{10}^{-29}{s}^{-1}$)within the two-year span of observations linear component reaches 10^{7} - 10^{8} s, square component is about 10 s, cubic component is not more than a few microseconds. Further continuation of the series shows that estimates the contribution of the third derivative is less than 1ns, and therefore cannot be detected on noise background, so higher derivatives of order 2 in a power series (1) cannot be considered at all.

In the relation (1) the intervals РТ are calculated from the parameters of rotation, which are known merely suspected a priori. However, by linking these parameters with the observed intervals PTi, we obtain an equation whose solutions are the observable parameters
${P}_{0}^{\ast},\stackrel{\dot{}}{P},\stackrel{\xa8}{P}$ , corresponding to the condition of convergence of the series (1). The equation of the observed intervals PT_{i} in accordance with the relation (1) is:

$P{T}_{i}=\left(1+{\alpha}_{i}\right){\left({P}_{0}^{\ast}N+\frac{1}{2}{P}_{0}^{\ast}\stackrel{\dot{}}{P}{N}^{2}+\frac{1}{6}\left({P}_{0}^{\ast 2}\stackrel{\xa8}{P}-2{P}_{0}^{\ast}{\stackrel{\dot{}}{P}}^{2}\right){N}^{3}\right)}_{i}$ (2)

Here are: PT_{i} is the numerical values of the observed intervals obtained from the planetary ephemerides;
${P}_{0}^{\ast},\stackrel{\dot{}}{P},\stackrel{\xa8}{P}$ are the observed pulsar rotation parameters obtained by solving Equation (2);
${\alpha}_{i}$ is the divergence measure of series (2) of the PT_{i} approximated by the rotation parameters of pulsar.

By parametric approximation of the intervals, PT_{i} (2) counted from the initial observed pulsar event, the fixed rotation period and its derivatives on the initial epoch in accordance with the power series expansion Maclaurin (1), are defined. The combination of parameters
${P}_{0}^{\ast},\stackrel{\dot{}}{P},\stackrel{\xa8}{P}$ is unique for each pulsar, wherein convergence of polynomial series (1) is reached. The 3-component series in the second brackets of the Equation (2) is an analytical form
$P{T}_{calc}$ of the observed intervals
$P{T}_{obs}$ , which corresponds to numerical value
$P{\u0422}_{i}/\left(1+{\alpha}_{i}\right)$ defining by the parametric approximation of PT_{i} by the calculated intervals with computing precision (inaccuracy less than 1 ns):

$P{T}_{obs}=P{T}_{i}$ ; (3)

$P{T}_{calc}={\left({P}_{0}^{\ast}N+\frac{1}{2}{P}_{0}^{\ast}\stackrel{\dot{}}{P}{N}^{2}+\frac{1}{6}\left({P}_{0}^{\ast 2}\stackrel{\xa8}{P}-2{P}_{0}^{\ast}{\stackrel{\dot{}}{P}}^{2}\right){N}^{3}\right)}_{i}$ (4)

In a generalized view of a ratio (3) and (4) for all events of radiation of a pulsar observed within any interval are expressed as follows:

$P{T}_{obs}=\left(1+{\alpha}_{PT}\right)P{T}_{calc}$ , (5)

The numerical value ${\alpha}_{PT}$ , as well as rotation parameters, is also the solution of the Equation (2). It characterizes unmodelled timing noise in the sequence of observed intervals $P{T}_{obs}$ , it isn’t correlated with an analytical pulsar time scale $P{T}_{calc}$ which is defined by the observed parameters of rotation only and, therefore, doesn’t depend on that unmodelled timing noise.

It has been shown that the Equation (2) is form-invariant under coordinate transformations, in which the numerical values of the observed rotation period are coinciding in the barycentric and topocentric coordinate systems at the same epochs of the local time. Left part of the Equation (2) consists the observed topocentric $T{T}_{obs}$ or barycentric $T{B}_{obs}$ intervals. The right part contains the intervals $T{T}_{calc}$ or $T{B}_{calc}$ , which are calculated by the observed rotation parameters obtained by approximation of $T{T}_{obs}$ or $T{B}_{obs}$ , respectively.

So, the analytical model of timing (2) determines the unique consistent values ${P}_{0}^{\ast},\stackrel{\dot{}}{P},\stackrel{\xa8}{P}$ corresponding to the coherence of periodic radiation of pulsar, on the observed ToAs only. Random variations of ToAs due to imprecise knowledge of physical variables and unmodelled noise of observations have no significant impact on the modelled intervals of coherent periodic radiation of pulsar, which are identified with a machine accuracy and sub-nanosecond resolution by the consistent rotation parameters of the pulsar.

As a result of dividing the observed intervals, the random variations are not contained detectable rotation components. They can be analyzed additionally to identify with affecting factors, for example, the distortion of an atomic time scale of the radio telescope, for possible physical interpretation. It is assumed also, that the pulsar rotation parameters are independent of any physical factors, which are involved in a parametrical fitting.

3. Consistency of the Observed Rotation Parameters of Pulsars

It is evident, that the value of the period ${P}_{0}^{\ast}$ depends on any choice of the epoch of the initial event, taking into account the nonzero derivatives $\stackrel{\dot{}}{P}$ , $\stackrel{\xa8}{P}$ . The corresponding settings of the rotation parameters also satisfy the convergence of the series expansion (2) for any extension in the vicinity the variable specified by $t={P}_{0}^{\ast}N$ . Included in the Equation (2) the numerical value ${P}_{0}^{\ast}$ represents the rotation period at the initial epoch of observations, and any current epoch t takes a value consistent with the derivatives $\stackrel{\dot{}}{P}$ , $\stackrel{\xa8}{P}$ :

$\u0420\left(t\right)={P}_{0}^{\ast}+{\stackrel{\dot{}}{\u0420}}_{0}\cdot t+\frac{1}{2}\stackrel{\xa8}{P}\cdot {t}^{2},\text{}t={P}_{0}^{\ast}N,\text{}1<N<\infty $ (6)

Relation (6) corresponds to a property of the coherence of pulsar periodic radiation and hence indicates monotonous slowing down of the pulsar rotation ( $\stackrel{\dot{}}{P}$ > 0, $\stackrel{\xa8}{P}$ > 0), which is an observational sign of convergence of the polynomial power series (1). Monotonic slowing down of the pulsar caused the negative values of derivatives $\stackrel{\dot{}}{\upsilon}$ < 0, $\stackrel{\xa8}{\upsilon}$ < 0 in the frequency domain, respectively.

The approximation of the observed intervals (2) also detects the random variations of the observed emission period ${\u0420}^{\ast}\left(t\right)$ . The magnitude of variation is expressed as follows:

$\Delta \u0420\left(t\right)=\alpha \left(t\right)\left({P}_{0}^{\ast}+{\stackrel{\dot{}}{\u0420}}_{0}\cdot t+\frac{1}{2}\stackrel{\xa8}{P}\cdot {t}^{2}\right)$ (7)

By integration of function in an interval {tj, ti}, (here i > j), the deviation of observed intervals is calculated:

$\Delta P{T}_{j,i}=\frac{1}{{P}_{0}^{\ast}}{\displaystyle \underset{{t}_{j}}{\overset{{t}_{i}}{\int}}\Delta P\left(t\right)\text{d}t}$ (8)

Replacing continuous function with its discrete approximation on the Equation (2), and assuming that change of the observed period in an interval {tj, ti} happens evenly, we obtain deviations of observed intervals in the form of the final differences between any two observed events of radiation of a pulsar within this interval:

$\Delta P{T}_{j,i}=\frac{\left({\alpha}_{i}+{\alpha}_{j}\right)}{2}\left(P{T}_{i}-P{T}_{j}\right)$ . (9)

The main characteristic of the observed slowing down of the pulsar rotation is the braking index n describing dependence of angular speed $\Omega =2\pi \upsilon =2\pi /P$ of rotation of a neutron star on time. The braking index is the only direct test connected with a slowing down mechanism of radio pulsars now. On the braking index it is possible to compare pulsars, the period and which derivatives differ on several orders, but being observed parameters, they are available to comparison according to catalogs and other sources. They it is equivalent are expressed in a frequency and time domain:

$n=\frac{\Omega \stackrel{\xa8}{\Omega}}{{\stackrel{\dot{}}{\Omega}}^{2}}=\frac{\upsilon \stackrel{\xa8}{\upsilon}}{{\stackrel{\dot{}}{\upsilon}}^{2}}$ , or $n=2-\frac{P\stackrel{\xa8}{P}}{{\stackrel{\dot{}}{P}}^{2}}$ (10)

Considering that the numerical values of the period P and its derivative $\stackrel{\dot{}}{P}$ were correctly measured for all the pulsars contained in the fundamental catalogs, there are followed the evident discrepancy for a pulsar slowing down evaluation. Instead of inspecting n = 3 for the magnetic dipole braking model accepted by the authors, there were obtained the measured braking indices for the non-recycled pulsars range from −287,986 to +36,246 with a mean of −1713 and a median of 22. If to restrict the sample of the 30-years data spanning, there were obtained braking indices ranging from −1701 to +36,246 with a mean of +3750 and median of +29. Out of a sample of 366 pulsars, 193 (53 per cent) have a positive $\stackrel{\xa8}{\upsilon}$ value and the remaining 173 (47 per cent) have negative $\stackrel{\xa8}{\upsilon}$ [5] .

Revealed abnormal dispersion of a braking index of pulsars in different sources demonstrate first of all discrepancies exactly of the second derivative challenging additional studying and specification exactly about it.

Parametrical approximation of observed intervals was applied to timing of pulsars В0809+74, В1919+21, В0834+06, J1509+5531, B2217+47, В0329+54, В0531+21 in 2007-2015 on the BSA radio telescope, which operated at frequencies close to 111.3 MHz. A 64-channel radiometer with a channel bandwidth of 20 kHz was used for the pulsar observations. The data were sampled at intervals of 2.56 or 1.28 ms. The BSA radio telescope, made up of a linearly polarized transit antenna with a beam size of about (3.5/cosδ) min, provided the duration of the observing session ranging from 3 to 11 minutes at different pulsar declination δ. During this time, individual pulses were summed synchronously with a predicted topocentric pulsar period to form the mean pulse profile in each 20 kHz channel. After dispersion removal, all the channel profiles were summed to produce a mean pulse profile for the given observing session. The topocentric arrival times of the pulses for each observing session were calculated by cross-correlating the mean pulse profile with a standard low-noise template.

By approximation of observed intervals according to the Equation (2) were obtained the consistency parameters of rotation $P,\stackrel{\dot{}}{P},\stackrel{\xa8}{P}$ which are given in Table 1 together with corresponding braking indexes n [7] .

As a component of the second derivative
$\stackrel{\xa8}{P}$ in the PT intervals, it does not exceed a few microseconds, even at the pulsars with relatively large values
$\stackrel{\xa8}{P}$ = 10^{−28} - 10^{−29} s^{−1}, it is detected not earlier than after two years of observations. The

Table 1. Rotation parameters and braking indexes of the observed pulsars.

pulsar B1919+21 (
$\stackrel{\xa8}{P}$ near 10^{−30} s^{−1}) requires at least 7 - 8 years of observations to detect
$\stackrel{\xa8}{P}$ . Even after 9 years of observations Pulsar B0809+74 within 2007-2009 and 2012-2017 it has not been found the meaningful sign of the second derivative, i.e.
$\stackrel{\xa8}{P}$ < 10^{−30} s^{-1}, which does not exceed a specified detection threshold.

Calculated from the observed parameters $P,\stackrel{\dot{}}{P},\stackrel{\xa8}{P}$ the value of braking index is a numerical invariant n = -(0.9 ± 0.2), corresponding to coherent radiation at monotonous slowing down of these pulsars.

4. Observations and Timing Analysis

Now, in view of the results of the pulsar timing expressed as analytical power series, we will expand a scope of this model, having extended it also to other observation data and objects for comparison of results and their physical interpretations. For this purpose, we will address the list from 27 pulsars, observed from 1978 to 2012 years on the BSA FIAN radio telescope by an identical technique in the same conditions [4] . All pulsars from Table 1, except PSR B0531+21 contain in the list 27 PSR of Shabanova as well. In turn, all pulsars from Shabanova’s list contain in the extensive list from 366 pulsars observed over the past 36 years using the 76-m Lovell radio telescope at the Jodrell Bank observatory, which are investigated in details regarding the existence of irregularities of rotation in a well-known work of Hobbs et al. [5] . Observations on the BSA and Lovell radio telescopes were carried out almost in parallel at the same years, authors were rather well mutually informed and actively discussed results, paying special attention to the available divergences, analyzed their possible reasons.

Timing parameters of 27 pulsars were determined by the initial topocentric ToAs performed on regular observations on the BSA radio telescope in the automatic mode. Timing parameters were determined using the program TEMPO 2 and the JPL DE200 ephemeris. The topocentric arrival times collected at the PRAO and the geocentric arrival times obtained from the archival JPL timing data were all corrected to the barycenter of the solar system at infinite frequency. The position and the proper motion required for this correction were taken from [5] and were held constant during the fitting procedure. The positions of most of the pulsars from our sample do not require correction since their residuals (defined as the observed minus predicted arrival times) obtained after the fitting procedure did not show a sinusoid with a period of 1 yr within the precision of the measurements. Position corrections were needed for pulsars B1508+55, B2020+28, and B2224+65, which are fast moving pulsars, but a timing solution did not provide improved positions for these pulsars because of the large timing noise. As the timing behavior of these pulsars was analyzed over a long interval, the presence of a weak annual sine wave at the end of this interval does not affect the large-scale structure of the timing residuals.

Further are given the results of timing of 4 pulsars, two of which (B1919+21 and B0809+74) contain both in Table 1, and in the list of 27 pulsars of Shabanova, a pulsar B0531+21 is only in Table 1, pulsar B1822-09 is given in the list of Shabanova only. Timing of all listed pulsars is performed on analytical model and observed rotation parameters of pulsars corresponding to it.

4.1. Pulsars B1919+21 and В0809+74

At these pulsars, apparently from Table 1, there are the smallest values of the second derivative, at PSR B0809+74 it doesn’t reach a detecting threshold even for several years of observations. Therefore, it is possible to hope that retrospective archives of observations of Shabanova lasting several decades will allow to extract the reliable value of the second derivative. On the example of PSR B1919+21 we will compare results of detecting of second derivatives (P2dot near 10 - 30 s^{−1}) on previous data sets of long-term observations of Shabanova and later observations on the BSA FIAN radio telescope.

1) B1919+21.

The results of our timing processing of PSR B1919+21 are presented in Figures 1-4. We get topocentric and by means of TIMAPR [2] barycentric ToAs from the observations of Shabanova during 1994-2012 that partially coincides in time with 8 years’ observations 2008-2009 and 2012-2017 by Losovsky of that pulsar in Table 1. We set a task to compare values of the second derivative coordinated with other parameters of rotation on both of these selections.

On analytical model we looked for equivalent decisions for the second derivative both in topocentric and in the barycentric coordinates. As duration of both selections is sufficient for sure detection of a cubic component, we can directly compare the results of detecting of the second derivative and test its consistency with other parameters of rotation. When detecting the second derivative we operated with absolute values of ToAs both in the topocentric system (TT), and in the barycentric system (TB). On graphics of Figure 1 their differences are shown.

The second derivative was determined according to ratios (3), (4) and (5) direct detecting by observed intervals of PT_{obs}. Considering that parameters P

Figure 1. Difference of topocentric (TT) and barycentric (TB) Times of Arrival [s] during 1994-2012 observations of PSR B1919+21.

and Pdot are known, the component of the second derivative has been just subtracted from PT_{obs} of a component of this interval calculated in the known parameters P and Pdot:

$P{T}_{obs}\left(P2dot\right)=P{T}_{obs}-P{T}_{calc}\left(P,Pdot\right)$ . (11)

Detecting of the second derivative was made equal in both systems, observed ToAs were counted is independent of the general initial impulse defining an initial epoch of observations of all subsequent impulses of TTobs and TBobs in each reference system.

The result of detecting of the P2dot component in topocentric and barycentric systems is given in Figure 2, and the result of cubic approximation is shown on Figure 3. On Figure 4 not modelled noise variations of observed intervals which aren’t correlated with rotation of a pulsar are shown. Their values (RMS) within almost 20-year interval doesn’t exceed 0.1 mks, and, therefore, signal to noise ratio (SNR) of the detected cubic P2dot component on Figure 2 appears not less than 250, it more than is enough for sure detecting of the second derivative confirmed by the coincidence of observed components of the second derivative and their cubic approximation in both reference frames.

The numerical value of the second derivative establishes approximation is 3.99573∙10^{−30} s^{−1}. It coincides with the P2dot value obtained in later observations in Table 1 and is coordinated with the first derivative and the period defined for any epoch, according to a ratio (6) with the same braking index n = -0.94. It is necessary to denote to Figure 3 on distinction of values of the period, which correspond to different initial epoch of MJD: in the topocentric reference frame P_{0}(TT) = 1.337302232070039186 s (MJD 49616.69854), in the barycentric reference frame P_{0}(TB) = 1.337302232070342499 s (MJD 49616.70114). The specified distinctions of the period are defined by the difference of an initial epoch of observations (TT-TB)_{0} = -224.4 s on Figure 1 due to ToAs difference in topocentric and barycentric reference frames which, according to the expression (6), define the corresponding values of the period for an initial epoch of observations in each reference frame.

It is easy to notice that the value of the period is defined much more precisely, than it is specified in the catalog. And the longer duration of the observations, the higher the accuracy corresponding to an approximation of observed intervals. With a duration of observations in several decades the inaccuracy of approximation

(a) (b)

Figure 2. Direct detection of second derivative P2dot in the observed intervals TTobs (a) and TBobs (b) [mks]. The cubic deviation of observed intervals reaches about 25 mks and surely exceeds a threshold of detection of the second derivative.

(a) (b)

Figure 3. Cubic parametric approximation of second derivative P2dot in topocentric (a) and barycentric (b) reference frames. Cubic component of P2dot, in accordance with Equation (4), is calculated with an achievable computer accuracy (uncertainty less than 1 ns). Second derivative P2dot = 3.99573∙10 s^{−1} - 30 s^{−1}. It is agreed with the rotation period P and derivative P2dot according to the catalog [1] , at braking index n = -0.94.

(a) (b)

Figure 4. Random variations of the observed intervals of PSR B1919+21 in topocentric (a) and barycentric (b) reference frames. The variations are uncorrelated with periodic rotation of pulsar and don’t exceed 0.1 mks within 1994-2012 observations.

reaches a limit of discrete resolution of the approximating calculator that in this case corresponds to picosecond permission of observed intervals, thus that not correlated noise of timing don’t exert impact on the result of the approximation.

Thus, by means of analytical models of timing (2) the compatibility of observed intervals of radiation in a topocentric and a barycentric reference frame is reached. Observed intervals, due to the coherence of radiation of a pulsar, are expressed through the rotation parameters identical for coinciding epoch and coordinated within any duration of the observations. Therefore, the periods and its derivatives received as a result of the approximation of observed intervals stay mutually coordinated within any duration, irrespective of the choice of an initial epoch of observations. It is only necessary to provide the duration of observations sufficient for sure detecting of the second derivative taking into account the required accuracy of the parametric approximation of observed intervals.

2) В0809+74

As it was shown upper, even after 9 years of observations of pulsar B0809+74 on the BSA FIAN radio telescope it has not been found the meaningful sign of the second derivative, i.e.
$\stackrel{\xa8}{P}$ < 10 - 30 s^{-1}, which does not exceed a specified detection threshold, as noted in Table 1. The absence of a significant second derivative at relatively short intervals of 2 - 3 years simplified the processing of the timing data when only the first derivative could be limited. However, the increasing demands on the accuracy of timing, raises the question of the exact definition of the second derivative, especially sharply, taking into account the anomalous variations of its existing estimates. Therefore, the unique 34-year observations pulsar B0809+74, performed by Shabanova et al. [4] , are of particular interest as a chance to overcome the limiting detection threshold and to obtain the actual measured value of the cubic component of the second derivative by using long-term archival data.

The method of processing observational data is completely analogous to that used for the pulsar B1919+21. In Figures 5-8 shows the same structure of graphs, calculations and estimates as in Figure 1 for the pulsar B1919+21. In the graph on Figure 5 shows the difference between the initial topocentric and barycentric ToAs within 34 years of observation. We note that two breaks in observations of 1 - 2 or more years in the initial stage had no effect on the timing results, since, firstly, all the intervals of the observed events are counted from the general beginning and, secondly, the intervals are approximated observable

Figure 5. Difference of topocentric (TT) and barycentric (TB) ToAs [s] during 1978-2012 observations of PSR B0809+74.

(a) (b)

Figure 6. Direct detection of P2dot in the observed intervals TTobs (left) and TBobs (right) [mks]. The cubic deviation of observed intervals only after 34 yrs of observations overcomes a threshold of detection of the second derivative, which is about 2, 5 mks.

(a) (b)

Figure 7. Cubic parametric approximation of second derivative P2dot in the topocentric (left) and barycentric (right) reference frames. P2dot = 5.89∙10 - 32 s^{−1}. The braking index n = -0.71, taking into account the agreed values of the rotation period and derivative Pdot in the Cataloque [1] .

(a) (b)

Figure 8. Random variations of the observed intervals of PSR B0809+74 in the topocentric (left) and barycentric (right) reference frames. The random variations of the observed intervals are within the range 0.2 - 0.3 mks and do not have a significant effect on the accuracy of the measurement of the second derivative. (a) RMS = 0.275 mks; (b) RMS = 0.159 mks.

parameters of rotation, agreed upon at any epoch, no matter what events were chosen for observations.

Compared to the PSR P1919+21, the timing of the PSR B0809+74 reveals the following features:

• Higher level of unmodelled timing noise and about an order of magnitude smaller the value of the detected component of the second derivative, relatively low value of S/N ≈ 10;

• The more accurate installation of the initial MJD epoch of observations in topocentric and barycentric reference frames as a condition for the equivalence of the cubic components of P2dot and the consistency of the initial period Po at the 34-year interval is defined;

• The value of P2dot is approximately 2 orders of magnitude smaller than that of the B1919+21 pulsar and 3 - 4 orders of magnitude lower than that of the pulsars B0834+06, J1509+5531, B2217+47, B0329+54 from the Table 1;

These features of the timing of the B0809+74 pulsar, which are directly related to its physical properties, are manifested in the conditions of data processing in the accepted parametric time processing model.

It should be specially noted that the braking index of the PSR B0809+74 is n = -0.71 and, unlike the other pulsars in Table 1, is close to the lower limit of the range n = -(0.9 ± 0.2) established by observation. So far, according to observations of only one pulsar, it is unclear whether this quantity is random or it reflects the general pattern characteristic of pulsars with a small deceleration, that is, pulsars with relatively small derivatives similar to PSR B0809+74. In general, such pulsars with a period of about a second, in which
$\stackrel{\dot{}}{P}$ £ 10 - 16 c∙с^{−1}, are not known so much. So, in the ATNF catalog, containing more than 2600 pulsars, there are only about two dozen of them. Several such pulsars are given in Table 2.

From the observational data on the timing of these pulsars, a direct measurement of the cubic component of the observed intervals can be used to measure the exact value of the second derivative to refine the lower limit of the braking index. It should only take into account that for a reliable detection of P2dot observations will be required for 40 years or more.

4.2. On the Consistency of the Rotation Parameters of Millisecond Pulsars

For millisecond pulsar, whose periodic radiation is also coherent, the derivative of the period
$\stackrel{\dot{}}{P}$ usually does not exceed 10^{−18} - 10^{−21} s∙s^{−1}. The expected magnitude of the second derivative
$\stackrel{\xa8}{P}$ , by the consistency of the rotation parameters according to the criterion of monotony slowing down (6), has to be within 10^{−35} - 10^{−39} s^{−1}. This is a lot of order below the detection threshold, and
$\stackrel{\xa8}{P}$ can be evaluated theoretically only, by analogy with calculation of the braking index of second pulsars. Therefore, proceeding from the numerical invariance of the braking index with an average value of the braking index n_{avrg} = -0.94, the numerical values of
$\stackrel{\xa8}{P}$ as their upper limit for millisecond pulsars have been calculated in the known rotation parameters
$P,\stackrel{\dot{}}{P}$ .

Table 3 summarizes the consistent rotation parameters of 7 millisecond pulsars. The parameters $P,\stackrel{\dot{}}{P}$ of the PSRs J1640+2224, J2145-0750, J0613-0200, J1713+0747, J1643-1224 were taken from the catalog ATNF [8] , these parameters of the PSRs B1937+21 and B1855+09 were published in the fundamental research Kaspi et al. [9] , special devoted to these pulsars, together with ToAs, which were kindly provided via the Internet by the authors.

The observations of Kaspi et al. 1994 were performed on observatory in Arecibo at frequencies of 1.4 GHz of PSR B1855+09 ((J1857+0943), and, respectively 1.4 and 2.4 GHz of PSR B1937+21 (J1939+2134). The method of precision timing the moments of observed impulses of PSR B1855+09 with an accuracy of 0.8 mks and PSR B1937+21 with an accuracy of 0.2 mks were measured. In Table 3 at these pulsars the format of an epoch differs in a large number of decimal signs of a fractional part whereas at all other pulsars of Table 3 the integer MJD only

Table 2. Second pulsars with small rotation deceleration.

Table 3. Consistency rotation parameters of the millisecond pulsars.

are specified. At PSR B1937+21, besides, the period format is expanded to 18 decimal signs. For these pulsars the second derivatives in frequency domain are specified:

for PSR B1937+21 F2dot = (13.2 ± 0.3) E-27 s^{−3} or P2dot= -(3.2 ± 0.7) E-32 s^{−1}

for PSR B1855+09 F2dot = - (1.0 ± 0.9) E-27 s^{−3} or P2dot= (2.87 ± 2.58) E-32 s^{−1}.

The calculation for a ratio (10) showed that the numerical values of the second derivative
$\stackrel{\xa8}{P}$ of the period of millisecond pulsars В1937+21, В1855+09, J1640+2224, J2145-0750, J0613-0200, J1713+0747, J1643-1224 are in the range 10^{−35} - 10^{−38} s^{−1}. It is at least 7 - 8 orders less, than at second pulsars, and is much lower than a detection threshold that excludes a possibility of measurement
$\stackrel{\xa8}{P}$ directly by the timing of millisecond pulsars.

In work of Hobbs et al. (2010) [5] from total number 366 PSR 25 millisecond pulsars are analyzed. They are selected in Table 4 together with their frequency parameters of rotation FO, Fdot, F2dot. Some of pulsars, namely J1939+2134, J1857+0943, J0613-0200 contain in Table 3 in which parameters of rotation it is equivalent in a time domain are expressed.

It is obvious that for millisecond pulsars at which the second derivative is a priori 7 - 8 orders less, than at second pulsars to find its numerical value in the observation way, even a parametric approximation of observed intervals, it is impossible in principle. The abnormal discrepancies in the definition of the second derivative are especially inherent in millisecond pulsars at which the relation of a cubic component in observed intervals to unmodelled variance of residuals objective is several orders less, than at ordinary pulsars.

Therefore the only opportunity to define the second derivative at millisecond pulsars is to establish it in the settlement way as the top limit coordinated with known for the period (frequency) of rotation and the first derivative on a numerical invariant of an braking index of n = -(0.9 ± 0.2) according to a ratio (10).

Table 4. Calculated F2dot (top limit) of millisecond pulsars in Table 1 (Hobbs at al., 2010).

Such computing procedure was performed in relation to all millisecond pulsars in Table 4. In the right column the table the calculated values of the top limit of F2dot received at the fixed value n_{avrg} = -0.94 are placed. Thereby the coordinated values of the second derivative in according to the sign, corresponding to monotonous slowing down, are obtained.

4.3. The Observed Spontaneous Movement of Pulsars in the Supernova Remnants

The specific features of pulsars associated with the supernova remnants, are larger magnetic fields, shorter spin periods, than the typical pulsars. These pulsars differ much bigger, than at usual second ones, in energy losses therefore the period derivative at them is more at least on 1 - 2 orders. Even before opening and the beginning of active observations of pulsars, F. Pacini [10] predicted that after explosion supernova has to be a very dense star kernel. Developing these representations, F. Pacini [11] noted that the rotation of neutron stars, which can be formed during the explosion of supernova, has to be connected with a surrounding gaseous cover.

V. Kaspi selected about 30 pulsars, which can be associated with the supernova remnants [12] . Important criterion of accessory are speeds of the movement of a pulsar, is emphasized that appropriate measurements of the movement are extremely important for confirmation of authenticity of many associations of pulsars listed in the table with supernova. Proper motion measurements for young pulsars are crucial for deciding whether many of the associations are genuine. The most striking conclusion following examination of the associations between neutron stars and supernova remnants is that the paradigm that every young neutron star is like the Crab pulsar requires reconsideration.

Meanwhile, it should be noted that the pulsar B0531+21, in particular, was not included in the number of pulsars 366 analyzed by Hobbs et al. [5] : “Although PSRs B0531+21,

3) PSR B0531+21

The observations of pulsar B0531+21 in the radio range 111 MHz detect an anomalous deviation the observation period P and intervals PT that do not fit into a convergent power series (2) of parameterized components observed PT intervals (Figures 9-12). However, the observed rotation parameters obtained from the Equation (2), matched with the criteria of coherence and correspond to the typical monotone spin down with braking index n = -(0.9 ± 0.2). Meanwhile, the relative magnitude of the observed deviations of the period P during the two-year observations reached the order of 10^{−9} (Figure 9(a)), which is 5 - 7 orders of magnitude greater than the other pulsars in Table 1. In addition, 21.02.2012 (MJD 55978), it has inverted the sign of the deviation of P, and further observation showed a sharp increase in its order of magnitude up to 10^{−6} [7] .

These spontaneous deviations of the observed emission period that is not associated with changes in the pulsar rotation parameters, can be interpreted as a Doppler shift due to the motion of the pulsar because of the effects of unbalanced forces in radial directions. The radial velocity is then determined by the deviation ${\alpha}_{i}$ of the period P obtained from Equation (2) and the light velocity c: ${V}_{i}={\alpha}_{i}\cdot c$ (Figure 9(a)). Figure 9(b) shows the deviation PT intervals and the corresponding displacement of the pulsar L along the line of sight.

Parameters of the movement―the speed, acceleration and movement of a pulsar, as well as rotation parameters, are defined by the Equation (2). Only,

(a) (b)

Figure 9. Observed movement of the pulsar В0531+21 during 01.2010-11.2011: $\Delta P/P$ ―deviation of the rotation period, V―velocity along the line of sight (a); $\Delta PT$ ―deviation of the observed intervals, L―observed motion of pulsar (b). Period P = 0.033403347409400 s (48743.0), Pdot = 4.2096∙10-32 [1] . Observed velocity smoothly increases up to value about 0.1 m/s, then it is stabilized near this level within the next year.

Figure 10. Observed velocity V of the PSR В0531+21 during 01.2010-09.2017. Observed changes of speed of the movement of a pulsar are unpredictable, they have sign-variable character and vary in range from -500 m/s to 3000 m/s and more.

Figure 11. Observed motion L of the PSR В0531+21. The pulsar moved in space more than 750 million kilometers within Crab Nebula during our observations from 01.2010 till 09.2017.

Figure 12. Observed acceleration ΔV/ΔPT of the PSR В0531+21. On the measured changes of velocity ΔV within intervals ΔPT between observations, the numerical values of acceleration of moving pulsar are calculated. They have spontaneous character and are in range approximately ± 3.5∙10^{−4} m/s^{−2}.

unlike other observed pulsars from Table 1, for PSR B0531+21 it is necessary to consider that the dimensionless parameter $\alpha $ of approximation of intervals (2) is associated with the movement of a pulsar here and is interpreted as the Doppler shift of the observed period of radiation owing to the movement of a pulsar under the influence of unbalanced forces in the radial directions. The relative shift of the period is expressed as follows:

${\alpha}_{i}=\Delta {P}_{i}/{P}_{0}^{\ast}$ .

The beam velocity of the movement of a pulsar is defined on an epoch of observation as follows:

${V}_{i}={\alpha}_{i}\cdot c$ , here сis light speed.

Displacement ${L}_{i}$ of a pulsar within some interval counted from an initial epoch in an assumption of constant speed is defined by the following expression:

${L}_{i}={V}_{i}\cdot P{T}_{i}={\alpha}_{i}\cdot c\cdot P{T}_{i}$ .

Displacement ${L}_{j,i}$ of a pulsar between randomly chosen observed events in an interval {РТj,РТi} is calculated on the ratio (9).

Over the span of nearly 6-year observations of pulsar B0531+21 within 2010-2015, its velocity changed in range, from 0.1 m/s or less up to 500 m/s or more; displacement of pulsar along the line of sight changed from a few hundred meters up to 100 million km and more. Calculated from deviation in velocity within the intervals of 1 - 7 days between observations, acceleration changed in range, from 10^{−8} m∙s^{−2} and less up to 10^{−4} m∙s^{−2} or more.

Note that the implied pulsar transverse velocity V_{t}, which is determined by the angular pulsar displacement from the SNR Centre, is specified by Kaspi [12] as a constant value 125 km/s. This is several orders of magnitude greater than the observed value, and does not consider it significant and unpredictable variations.

The observed deviation of the emission period in the pulsar-driven Crab Nebula shows that part of the rotational energy of PSR B0531+21 is converted into energy of the observed accelerated motion, without any violating the monotony its slowing down and coherence of the periodic radiation.

4) PSR B1822-09

This pulsar from a list of 27 PSR of Shabanova et al. (2013) [4] , drew the attention of a derivative Pdot = 5.2432E−14 s∙s^{−1} which 1 - 2 orders more than value, usual for second pulsars, and there is closer to Pdot = 4.096E−13 s∙s^{−1} of a pulsar B0531+21 which is located in the supernova remnant.

Shabanova, as she notes, was the first to detect slow glitches of the pulsar B1822-09. Slow glitches, according to Shabanova et al. [4] , represent a unique glitch phenomenon that results in regular oscillations in the rotation frequency. A full pattern of the timing and frequency residuals for the pulsar B1822-09 over the 26 year span between 1985 and 2012 show that the pulsar experienced three discrete glitches and suffered a sequence of five slow glitches that occurred over the 1995-2004 interval. Hobbs et al. [5] , concerning this subject, comment on it in connection with others, including own works: “Zou et al. (2004) reported a phenomenon known as ‘slow glitches’ for PSRs B1822-09 and J1825-0935. ‘Slow glitches’ are characterized by a permanent increase in frequency, but no significant change in the slow-down rate. More recently, Shabanova (2007) identiﬁed a further ﬁve slow glitches in the timing of PSR B1822-09. The observed timing residuals are not dissimilar to those described in our sample of pulsars in Section 3.2.3. We therefore suggest that slow glitches are not a unique phenomenon, but are caused by the same process as the timing noise”.

We carried out the processing of observation data of a pulsar of B1822-09 on the analytical model of time processing, similar to PSR B0531+21, having taken as initial the topocentric ToAs PSR B1822-09 of Shabanova within an interval of 1994-2005. Besides, to this series, we added data of later observations, performed on BSA FIAN in 2015-2016. Results of processing are given in the Figures 13-15.

On the coordinated values of parameters of rotation P_{0} = 0.6897223949 s (MJD 46800.0) and Pdot = 5.432E−14 s∙s^{−1} from the PSR Catalogue [1] , P2dot = 3.66E−27 s^{−1} that corresponds to an braking index n = - 0.94 and to monotonous slowing down of rotation were found and measured the spontaneous deviations of the observed emission period which are interpreted in this case by analogy with the PSR B0531+21 as a Doppler shift due to the motion of the pulsar. Observed spontaneous deviations of the period fall just on an interval 1994-2005

Figure 13. Observed velocity V of the PSR В1822-09. For processing data of observations of Shabanova during 1994-2005 together with later data of observations on the BSA radio telescope in 2015-2016, were taken. The intensive growth of velocity within the first two years was replaced by gradual subsequent reduction of growth up to almost constant value V near 350 m/s after 5 - 6 years.

Figure 14. Observed motion within 09.1994-02.2016. The pulsar moving after sharp acceleration with actually constant speed more than 15 years overcame distance about 250 million kilometers during this time.

(a) (b)

Figure 15. Observed acceleration of the pulsar. The measured numerical value of acceleration on the active state (a) reaches an order 10^{−4} m/s^{−2}, whereas in rather quiet state (b) the value of acceleration is two orders less.

to which the estimated glitch or noise of timing specified by the author’s of this pulsar belongs.

Further observations of a pulsar of B1822-09 in 2015-2016 showed that the beam speed of a pulsar remained without significant changes and after a 10-year break in observations. Respectively acceleration became several orders less in comparison with an interval of 1994-2000 (Figure 15(b)), that is after about 6 years of the expressed dynamics at the beginning of observations there came a rather quiet pause during which almost unvarying motion of a pulsar observed in the beam direction continues.

5. Summary

The analytical timing model, based on the parametric approximation of the numerical series ToAs by the convergent power series (2), is an effective and reliable tool for direct extraction and testing of the consistent physical parameters of the pulsar rotation. For each pulsar, there is a single combination of numerical values $P,\stackrel{\dot{}}{P},\stackrel{\xa8}{P}$ ( $\upsilon ,\stackrel{\dot{}}{\upsilon},\stackrel{\xa8}{\upsilon}$ ) that are the same in any arbitrarily chosen observing reference system within an unlimited duration of observations, both current and retrospective.

The consistency of the rotation parameters does not depend on the choice of the initial epoch; in any case, the initial period value corresponding to the chosen sample of the initial epoch, taking into account the derivatives of the period (6), is determined by approximating the intervals PTobs.

A prior knowledge of the parameters of pulsar rotation from previous observations, confirmed by subsequent observations, indicates the coherence of periodic pulsar radiation and the necessary condition for the convergence of the observed PTobs series.

The main problems in approximating the observed intervals were related to the exact determination of the second derivative, since the cubic component of the polynomial power series of intervals PTcalc, due to a very small, only a few microseconds or less, reaches the threshold of reliable detection only after several tens of years of observations.

A direct consequence of the consistency of parameters ${P}_{0}^{*}$ , $\stackrel{\dot{}}{P}$ , $\stackrel{\xa8}{P}$ of coherent radiation is the monotonicity of the deceleration of pulsar rotation, expressed by the coincidence of the sign of the derivatives-positive in the time domain ( $\stackrel{\dot{}}{P},\stackrel{\xa8}{P}$ > 0) and negative in the frequency domain ( $\stackrel{\dot{}}{\upsilon},\stackrel{\xa8}{\upsilon}$ < 0). Accordingly, the braking index (10) for the observed pulsars is close to the negative unit: n = -(0.9 ± 0.2) and varies in a rather narrow range of values, while the rotation parameters they differ by several orders of magnitude. In essence, this is a numerical invariant for all pulsars, and propagates to those about them, in which the second derivative does not reach the detection threshold in observations and is calculated from the condition of consistency of the three rotation parameters in the time or frequency domain. The analytical model (2) excludes the anomalous distortions of the second derivative, inherent in its estimates from the barycentric residual deviations,

The initial numerical series ToAs in the polynomial power series of intervals PTcalc, transformed with machine accuracy, is an analytical time scale from which the random, unmodified variations of the observed intervals PTobs are measured, as in the pulsars B1919+21 and B0809+74, as well as the pronounced spontaneous variations of the intervals defining the movement of PSR B0531+21 pulsar in the supernova remnant and PSR B1822-09, the spontaneous movement of which also suggests that it belongs to the supernova remnant.

All observation data of 27 pulsars kindly allowed T.V. Shabanova for the staff of Pushchino Observatory became available for analysis in scientific research after considerable work on their systematization, which is executed by her colleague V.D. Pugachev.

Acknowledgements

Authors are thankful to the technical and scientific staff of the Pushchino Radio Astronomy Observatory for the active collaboration and comprehensive supporting in caring out of the regular long-term observations of the pulsars.

Conflicts of Interest

The authors declare no conflicts of interest.

[1] |
Taylor, J.H., Manchester, R.N. and Lyne, A.G. (1993) Catalog of 558 Pulsars. Astrophysical Journal Supplement Series, 88, 529-568. https://doi.org/10.1086/191832 |

[2] | Doroshenko, O.V. and Kopeikin, S.M. (1990) High Precision Pulse Timing for Single Pulsars. The Astronomical Journal, 67, 986-992. |

[3] | Lorimer, D.R. and Kramer, M. (2005) Handbook of Pulsar Astronomy. Cambridge University Press, Cambridge. |

[4] |
Shabanova, T.V. Pugachev, V.D. and Lapaev, K.A. (2013) Timing Observations of 27 Pulsars at the Pushchino Observatory from 1978 to 2012. The Astrophysical Journal, 775, 1-24. https://doi.org/10.1088/0004-637X/775/1/2 |

[5] |
Hobbs, G., Lyne, A.G. and Kramer, M. (2010) An Analysis of the Timing Irregularities for 366 Pulsars. Monthly Notices of the Royal Astronomical Society, 402, 1027-1048. https://doi.org/10.1111/j.1365-2966.2009.15938.x |

[6] | Avramenko, A.E. (2015) Parametric Invariance of the Relativistic Coordinate Pulsar Time Scales. Proceedings Journées Systèmes de Référence Spatio-temporels, St. Petersburg, 81-82. |

[7] | Avramenko, A.E. and Losovsky, B.Ya. (2016) A Continuous Observing Rotational Stability of Pulsars. Proceedings of Russian Astrometric Conference, Pulkovo, 11-16. |

[8] | Manchester, R.N., Hobbs, G.B., Teoh, A. and Hobbs, M. (2005) The ATNF Pulsar Catalogue. Astronomical Journal, 129, 1993-2006. |

[9] |
Kaspi, V.M., Taylor, J.H. and Ryba, M.F. (1994) High-Precision Timing of Millisecond Pulsars. III. Long-term Monitoring of PSRs B1885+09 and B1937+21. The Astrophysical Journal, 428, 713-728. https://doi.org/10.1086/174280 |

[10] |
Pacini, F. (1967) Energy Emission from a Neutron Star. Nature, 216, 567-568. https://doi.org/10.1038/216567a0 |

[11] |
Pacini, F. (1968) Rotating Neutron Stars, Pulsars and Supernova Remnants. Nature, 219, 145-146. https://doi.org/10.1038/219145a0 |

[12] |
Kaspi, V.M. (1998) Neutron Star/Supernova Remnant Associations. http//arxiv.org/pdf/astro-ph/9803026.pdf |

Journals Menu

Contact us

customer@scirp.org | |

+86 18163351462(WhatsApp) | |

1655362766 | |

Paper Publishing WeChat |

Copyright © 2022 by authors and Scientific Research Publishing Inc.

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.