An Improved Wiener Filtration Method for Constructing the Ensemble Pulsar Timescale

An improved Wiener filtration (WF) method for constructing the ensemble pulsar timescale (EPT) is proposed to solve the existing problems in the WF algorithm. The improvements consist of three parts: (i) adjusting the cross-power spectral density (PSD) and PSD estimated by the weighted average (WA) algorithm by the concept of the energy upper limit; (ii) combining the cross-PSD and WA PSD to form a different signal modulus for the filtration of residuals of different noise levels; and (iii) setting a weight for each residual by a more general algorithm based on the concept of effective power. We use this improved WF method to construct the EPT by both simulated data and observational data from the second data release of the International Pulsar Timing Array. The results from the simulated data indicate this improved WF can successfully suppress the noise level and reform the common signal. For observational data, this method also successfully detects the fluctuations of International Atomic Time (TAI) and local atomic time TA(NTSC). A stability analysis shows that an EPT will have the potential to reveal the instability of TAI on a scale longer than 7 yr. This improved WF method can also be used to detect and monitor the noise in pulsar timing residuals inversely, by using multiple atomic time standards, which can in turn improve the EPT.


Introduction
Millisecond pulsars (MSPs), with brilliant spinning stability, are the most stable natural frequency sources that provide a way to investigate the performance (especially in the long term) of timescales based on quantum theory in the solar system.In history, many authors (Petit & Tavella 1996;Rodin 2006Rodin , 2008;;Hobbs et al. 2012Hobbs et al. , 2020;;Rodin & Fedorova 2018;Píriz 2019) have proposed various methods to construct the ensemble pulsar timescale (EPT).An effective EPT can keep a lower noise level compared with the pulsar timescale (PT) established by a signal MSP and hence it behaves more stable.In addition to time metrology, an EPT is also used as a probe in astronomy and astrophysics, such as correcting solar system ephemeris, detecting nanohertz gravitational waves, and testing theories of gravity.
Today, timekeeping plays a significant role in scientific explorations and engineering applications.In this work, the atomic clock is at the center.The time laboratory of each country independently maintains its own atomic time TA(k), where k is the indicator of the laboratory.By collecting time data from different laboratories around the world, the Bureau International des Poids et Mesures (BIPM) calculates International Atomic Time (TAI), a paper timescale released monthly that conforms to the second definition.Apart from TAI, BIPM annually recalculates a terrestrial time (TT) named TT(BIPMXX)-where "XX" is the last two digits of the year this TT belongs to-which makes small corrections of TAI and presently represents the most precise atomic timescale realized on the Earth (Petit 2014).The common problem of atomic timescales is that their stabilities will degrade with device aging and error accumulation.An EPT, established by MSPs, provides one possible way of detecting these variations.
The key point of EPT construction is to extract the common signal among multiple pulsar timing residuals.In 1996, Petit & Tavella (1996) proposed the weighted average (WA) algorithm.Rodin (2008) pointed out the validity of WA will degrade when low-frequency noise exists in timing residuals.Apart from WA, cross-correlation is another commonly used method to search for a common signal.Based on this, Rodin (2006Rodin ( , 2008) ) and Rodin & Fedorova (2018) proposed a kind of brief algorithm to construct an EPT by using Wiener filtration (WF) theory.The algorithm is performed in the time domain in Rodin (2006Rodin ( , 2008) ) and in the frequency domain in Rodin & Fedorova (2018).However, practical calculations indicate that when some MSPs demonstrate significant noises, an EPT constructed by WF will fluctuate higher than the reference timescale.To suppress the effect of white noise, Yang et al. (2022) improved the WF algorithm in the time domain by using a power-law power spectral density (PSD) model and successfully reduced the fluctuations of the EPT.After carefully examining the problems of the WF algorithm in the frequency domain, we find that there exist some problems in extracting common signals.This paper is going to analyze these problems and propose significant improvements.The contributions of this paper include: (i) Indicating the problems of cross-PSD estimation in the classical WF algorithm and proposing the concept of an energy "upper limit," which corrects the amplitude of the signal recovered from pulsar timing residuals.This paper also make corrections to use the real part of the cross-PSD, which is the unbiased estimator of the common-signal PSD, to do filtration.
(ii) Combining cross-correlation and WA to search for the true PSD of common signals.This paper also corrects the Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
filtration method for residuals of different noise levels.As residuals of high-level noise generally cannot provide a true PSD estimation of the clock signal and may destruct the filtration of residuals of low-level noise, we need to treat the two conditions differently.
(iii) Proposing a new way for setting weight, based on the concept of effective power.
Using simulated data and observational data from the second data release of the International Pulsar Timing Array (IPTA), we demonstrate that this improved WF method successfully suppresses the noise of the EPT.The structure of this paper is as follows.Section 2 reviews the basis of pulsar timing, WA, and the WF algorithm.Section 3 analyzes the WF algorithm and proposes improvements to the filtration algorithm by using simulated data.In Section 4, we establish an EPT by using the observational data with the help of an autoregression (AR) model and conduct analysis of the results.Finally, Section 5 draws the conclusion.

Basic Knowledge
In this section, we review the basic knowledge of pulsar timing and the EPT construction method-WA and the WF algorithm, which are the cornerstones of EPT construction in this paper.

Basis of Pulsar Timing
Before taking the rotation of an MSP as frequency source, a procedure named pulsar timing needs to be run to correct time delays caused by various physical phenomena.Briefly, time corrections can be described as (Hobbs et al. 2006)

VP
where t obs and t refer to the time of arrival (TOA) at the observation site and the emission time of the pulses, respectively.In the right-hand side of Equation (1), each term corresponds to the clock correction, atmospheric propagation delay, solar system Einstein delay, solar system Roemer delay, solar system Shapiro delay, dispersion delay, excess vacuum propagation delay, and time delay caused by binary motion, respectively.The rotation of an MSP is generally expressed by a part Taylor series as where f is the rotation phase and f 0 is the phase at a certain epoch, and ν and  n are the rotation frequency and its derivative at the epoch.For MSPs, higher-order derivatives of ν are treated as noises, so are not included in Equation (2).The difference between the observed rotation phase f i , obtained by Equations (1) and (2), and the theoretical one I i , which should be an integer nearest to f i , forms a residual: Then the timing parameters can be updated by minimizing the weighted timing residuals.Lorimer & Kramer (2005) and Edwards et al. (2006) describe pulsar timing in more detail.As the computation concerned is complicated, the Australia Telescope National Facility developed a tool, tempo2,5 which is quite popular in pulsar timing analysis.In this paper, we use tempo2 to process pulsar timing data.

WA Algorithm for Constructing EPT
WA was proposed by Petit & Tavella (1996).The basic principle is very simple.Consider the pulsar timing residual r is clock offset PT-AT, where AT denotes atomic time, which is the reference timescale (also called the reference clock) of pulsar timing.Then the EPT can be described as where μ k is the weight for the pulsar k.In this paper, we use n to denote the number of pulsars in EPT construction and μ to represent the weight in EPT construction.Generally, the weight in WA is set by the rms (or weighted rms) σ rms : where σ rms,i is the rms (or weighted rms) of the residuals for the ith pulsar.
The principle of WA is under the assumption that if the noises induced by the PTs of different pulsars are independent, then they can be canceled out by each other by addition, and the common signal induced by AT will become clear in EPT-AT.

WF Algorithm for Constructing EPT
WF is a kind of optimal filtration dealing with wide stationary stochastic processes (SPs).In the construction of an EPT, all the data are post-processed, so noncausal filtration is involved.The frequency response of a noncausal WF is (Ding et al. 2007) Starting from here, we use x and s to denote the observational data and expected signal, respectively (in certain cases, x and s refer to the timing residuals and reference clock signal, respectively), ε to denote the noise in timing residuals, S to denote the PSD, and  and -  1 to denote the Fourier transform and its inverse transform.e i ω or e i2 π f is a special notation used to mark the PSD symbol in equations, where e, i, ω, and f refer to the Napierian base, imaginary unit, angular frequency, and frequency, respectively.ω and f satisfy ω = 2πf.In the frequency domain, WF can be described as We can understand Equation (7) as this: consider S xs and S xx are calculated by a periodogram, so the multiplication of S xs (e i ω ) and { }  x produce the multiplication of S xx (e i ω ) and { }  s , then s is recovered by the inverse Fourier transform.
In EPT construction, the expected signal is the reference clock, which is unknown and is estimated by the cross-correlation method.Rodin & Fedorova (2018) proposed the filtration method as where x a and x b refer to two different timing residuals, ¯( ) is the arithmetic average of all the cross-PSDs, and s k is the expected signal waveform recovered from the residuals x k .Then the EPT is formed by the WA of each signal s k (Rodin & Fedorova 2018): where μ k is the normalized weight of s k .Equation (8) can also be understood intuitively: consider S is evaluated by a periodogram and the divisor (data amount) has been canceled in the frequency response, then is a unit circle in the complex plane, so only provides phase information in the inverse Fourier transform, while ¯( ) is the estimation of the signal modulus.
We will offer a brief review of the cross-correlation method used in WF.
In this paper, we use N to donate the length of the residuals' time series and * to donate complex conjugation.the noises ε a and ε b are produced by different physical mechanisms and therefore should be independent, while s is generated by the reference timescale and therefore should be the same for all MSPs.If we take E to represent the mathematical expectation in this paper, and assume E(ε) = 0, we have e e e e ee just represents the PSD of the expected signal.

Analysis of and Improvements to WF Method
This section conducts analyses of and discussions about WF for the EPT, based on simulated data, and contains four parts.The first part introduces the properties of the simulated data in use and their generating principles.The second part shows the results of WF and analyzes the reasons leading to them.In the third part, the improvements to WF are then proposed.The last part exhibits the effects of this improved WF method, by comparing with the results of classical WF.

Generating Simulated Data
The reason for using simulated data for the analysis is that we can generate different equispaced timing residuals with uncorrelated noise and a precisely equal common signal.The common signal added to the simulated data is TT(TAI)-with the quadratic trend removed-in the Modified Julian Date (MJD) interval 50209-55579, and the result can be compared to that obtained by observational data later on.To make the simulated data satisfy the condition for comparison, the following principles are followed.
(i) To obtain equispaced residuals that are strictly aligned to the point of time of the clock offset TT(BIPM15)-TT(TAI), in the step of generating the simulated data, the observational cadence is set to be about 5 days per TOA.Then we use piecewise linear interpolation to get the residuals at each time point.
(ii) Apart from white noise, we also consider adding different red noise to the simulated data.White noise and red noise are added to the residuals by the fake plugin and simRedNoise plugin of tempo2.The properties of the red noise are characterized by the PSD where A p , γ, and f c , respectively, refer to the amplitude, index, and corner frequency.
(iii) We make six sets of simulated data and denote the Xth one to be SDX.The noise parameters of the simulated data are listed in Table 1.SD1 is purely white, whose rms is 100 ns.SD2-SD4 have the same strength white noise as SD1, but consist of different types of red noise.SD2 consists of a weak red noise, whose PSD parameters come from the noise model of J0437-4715, in the first data release of the IPTA Version B data set.SD3 and SD4 both have a stronger red noise, with the PSD parameters coming from the noise models of J1744-1134 and J1713 + 0747, respectively, but the noise amplitudes are both 100 times stronger than the original values.SD5 is purely white, with 1 μs rms.SD6 has the same strength white noise as SD5, but consists of a weak red noise of the same PSD as SD2.
(iv) In generating the simulated data purely with noise, we set the reference clock to be TT(BIPM15).Then the pre-fit timing residuals are TT(BIPM15)-PT.The quadratic trends in each set of residuals are fitted and removed to keep each SD satisfying wide stationary SPs, and a common signal (clock signal) is added to the timing residuals by the operation TT TAI TT BIPM15 PT PT TT TAI .14 Figure 1 depicts the six sets of timing residuals obtained by Equation (14).

Problems and Discussion of WF Method
Figure 2 depicts EPT-TT(TAI) established by two different data sets: (i) only SD1 and SD2; and (ii) all the SDs.In constructing the EPT, weight is set to be inversely proportional to the rms square of each SD.The error bars in Figure 2 are calculated using the method described in Yang et al. (2022), and the same goes for other places.From the figure, we can conclude that the inclusion of SDs with high-level noise totally deactivates the effectiveness of the EPT.We explain this phenomenon in two aspects: (i) By comparing the PSD of the clock offset TT(BIPM15)-TT(TAI) and ¯( ) estimated by the SD, we find the latter exhibits very different energy distributions when all the SDs are used.In this condition, ¯( ) is much higher in the lowfrequency interval that covers the energy of TT(BIPM15)-TT (TAI) and also has a large energy peak in a higher-frequency interval.This difference can be understood as follows.That the cross-PSD of different residuals equals the PSD of the clock signal is theoretically right in the sense of mathematical expectation.In practical calculations, Equation (10) is the estimation of the cross-PSD by time average under the assumption of ergodicity.For an SP with finite length, Equation (10) may not give a reliable estimation, but rather provide an energy upper limit of the cross-PSD by the following truth (Feng & Bo 2012): e S e S e . 1 5 Hence, with the inclusion of residuals with strong noise, ¯( ) x x i a b will give a much higher energy estimation of the clock signal.
(ii) The phase of signal s k is the phase of . It is also disturbed by noise.Consider that pulsar timing residual x is the addition of signal s and noise ε, so we have Hence, when the energy of ε is more weak than that of s, the phase of { }  x can reflect the phase of { }  s without too much distortion; while the energy of ε is much stronger, { }  x will not reflect the phase of { }  s at all.This effect is shown in Figure 3. Therefore, when using pulsar timing residuals having very strong noise to do the filtration, even when we have the PSD of s exactly, the waveform of s still cannot be recovered from the residuals.

Improvements to WF Method
According to the analysis above, it is clear that the inclusion of strong noise timing residuals will greatly affect the EPT.Rodin (2006) developed the concept of inverse WF, which uses multiple atomic clocks to filter out fluctuations in pulsar timing residuals, as a very good noise cancellation method.In this paper, we do not dabble in the noise cancellation problem, but propose a method to reduce the noise level.We have known  that different MSPs play different roles in EPT construction.For clarity, rewrite Equation (8) as where { } s modulus k is the estimation of the signal modulus for the residuals x k .The relation between { } s modulus k and the PSD is connected by Equation (10).In the following, we propose the improvements in three aspects.

Decision of Signal Energy
The result obtained by using six SDs in Figure 2 implies the estimation of signal energy is incorrect.Based on the discussions in the previous subsection and Equation (15), it is reasonable to put forward an energy-estimating principle: the modulus of each cross-PSD should not exceed the square root of multiplication of the auto-PSD of the two sets of timing residuals having the lowest noise levels.Thus, under the condition there are two sets of timing residuals that can provide the correct estimate of the magnitude of the signal PSD, we can use this estimation (which we call the upper limit) to adjust each cross-PSD, especially those obtained by timing residuals having very strong noise.
To speak in detail, the first step is to calculate the auto-PSD of each set of residuals and get the square root of the product of each PSD pair.Then we can make a rough check by observing the energy level.When there are several possible upper limits, and we cannot confirm which to choose by direct observations, a further comparison of the main energy area should be made.As we know, the energy of the clock signal concentrates in a low-frequency interval, hence for residuals with relatively weak noise, the form of the upper limit should be like Figure 4 depicts.The comparison of the main energy region is quantified by two types of inflection points: ω l , ω h and ω lh , ω hh .ω l , ω h is defined as exceeding the frequency point in the low-frequency (l) or high-frequency (h) direction, so the PSD can be regarded as spectral white noise.ω lh , ω hh is the low-and high-halfheight-frequency point of the main peak of the PSD in the main energy region.For the form in Figure 4(a), the main energy region can be defined as [0, ω hh ] or [0, ω h ]; while for the form in Figure 4(b), [ω lh , ω hh ] or [ω l , ω h ] can be used.Precisely speaking, how to define the main energy region depends on the extent to which the cross-PSD is constrained.In practice, it is not difficult to try and compare.

PSD Estimation
Here, PSDs estimated by cross-correlation and WA are both used.
(i) PSD obtained by cross-correlation.
In general, a cross-PSD is a complex function of ω.The exact cross-PSD estimated by periodogram is Rubiola & Vernotte (2010) pointed out the unbiased estimator for the PSD of the common signal is where  refers to the real part of a complex number.Hence, in our analysis, we use Equation (19) but not Equation (10) to estimate the cross-PSD.
(ii) PSD obtained by WA.
Apart from cross-correlation, WA is another way to extract the common signal, hence, it can be used to make an estimation of the signal PSD.When the common signal is rather large compared to the noise of several sets of pulsar timing residuals, WA can provide a rather good EPT.However, as the noise in the timing residuals increases, this PSD obtained by WA will be distorted and enlarged greatly.In this case, this PSD should be adjusted with other cross-PSDs together.
(iii) Construct the PSD for WF.
For each PSD obtained in (i) and (ii), it should be adjusted generally.The method to confine the PSD in the main energy region is dividing it by a factor.However, the upper limit also has error.Practically, if the PSD obtained in (ii) keeps a very similar level, or has the same order of magnitude as the upper limit, it can also be reserved with no adjustment.After this, each PSD in (i) and (ii) will have the correct level under the assumption.Then the PSD used for filtration is formed by  where S WA (e i ω ) is obtained by WA and S s,i (e i ω ) is obtained by cross-correlation, and m is the total number of PSDs used to form S fil (e i ω ).In theory, the values of modulus{s k } for different residuals to establish a best EPT should be different.But practically, it is more feasible to divide the PSD for filtration into two categories: one to conduct filtration for the residuals forming the upper limit and another to conduct filtration for other residuals.Residuals forming the upper limit have lower noise levels and can reserve the common signal more factually.
The cross-PSD obtained by these residuals will be more like the true PSD of the common signal.For these residuals, S fil (e i ω ) is obtained by setting S s,i (e i ω ) to take the values obtained by these residuals.For other residuals, S s,i (e i ω ) takes all values of the cross-PSD.Apart from the above principle, there are two other aspects we need to take care of.The first one is if there exist some negative values in S fil (e i ω ), these values should be set to zero.The other is if some values of S fil (e i ω ) are pretty large outside the main energy region, these values can also be further suppressed by the upper limit to acquire a smoother curve of the EPT.

Setting Weights for Residuals
Weight reflects the significance of the residuals in EPT construction.The commonly used methods to set weight include: (i) set μ k to be equal with each other; (ii) set μ k inversely proportional to the rms square of the residuals; and (iii) set μ k to be inversely proportional to the stability value of the residuals.However, we must be careful that the EPT should be sensitive to the changing of the reference clock.Method (i) will make the weight of the residuals of strong noise too large.Methods (ii) and (iii) are both valid for the case in which the signal is much weaker than the noise in the residuals.
Take note of the fact that residuals of different noise strengths have different sensitivities to a certain signal, so the following truths are the basic principles in setting weight: (i) when one set of residuals is pure and noiseless, the weight of this set residuals should be 1 and we do not need to consider other residuals any more; (ii) when the noise strength of one set of residuals is approaching infinity, the weight should be approaching 0; and (iii) when the signal strength is becoming more and more strong and approaching infinity, the unnormalized weight should be approaching 1. Denote P s and e P k as the effective power of the signal and noise for the kth set of residuals x k , respectively, and we can express the weight as In practice, we do not know the exact values of P s and e P k .So we first make a guess of the signal amplitude, which to a great extent relies on a priori knowledge.When the signal strength can be treated as constant from the start to the end, we can use a simple signal, such as a sinusoidal signal, to estimate P s .Here, P s is calculated by where ŝ is the assumed signal.To calculate e P k , we need to subtract the signal from the residuals.Hence, e P k can be expressed by Here, + or − represents the direction of the signal, and x k,i is the ith value for the residuals x k .To avoid underestimating the weight, we can choose the sign that makes e P k smaller.The last column of Table 1 lists the weights for the simulated data.The signal added to each residual is a whole-period sine wave of 200 ns amplitude.

Effects of Improved WF Method
Figure 5 depicts the result of the improved WF.We can see the inclusion of SD3-SD6 not only does not result in distortion, but also lowers the stochastic noise of the EPT.In EPT construction, the upper limit is decided by SD1 and SD2 and the main energy region is defined as [ω lh , ω hh ].The energy of the PSD outside the main energy region is also suppressed by the upper limit.The rms of the EPT residuals ([EPT-TT(TAI)]-[TT(BIPM15)-TT(TAI)]) are: (i) the classical WF method using SD1 and SD2-75.8ns; and (ii) the improved WF method-55.7ns.

Improved WF Method Applied to Observational Data
In this section, we apply the improved WF method to observational data from the second data release of IPTA.We first introduce the data pre-processing method and stress the difference between observational data and simulated data.Then we use a new method to estimate the upper limit for some data sets.Finally, we apply this improved WF method to observational data to construct the EPT and conduct some discussions.

Data Pre-processing
Pulsar timing observations provide a quite uneven TOA distribution.To allow discrete Fourier analysis to be applied, it is necessary to transform uneven data so that it is equispaced.In this section, we use the average method with an average span of 30 days to achieve this task.Compared with interpolation, the average method can suppress some high- frequency noise.Though it disturbs the noise properties of the timing residuals, the spectrum of the clock signal will remain unchanged (see the Appendix) if the frequency corresponding to the average span is much higher than the bandwidth upper limit of the clock signal.Figure 6 depicts the frequency spectrum of the original clock offsets TT(BIPM15)-TT(TAI) and TT(BIPM15)-TA(NTSC) we use in this paper.The frequency corresponding to an average span of 30 days is about 6.08 yr −1 ( » log 6.08 0.78), which we can see is high enough to clock offset.The upper panel of Figure 7 depicts the timing residuals of J0437-4715 in data set 2 (DS2) before and after using the average method, while the lower panel shows the corresponding stabilities evaluated by σ z (τ) (briefly σ z ), a statistic described in detail by Matsakis et al. (1997), with τ being the measured interval length.In Figure 7, the last values of σ z curves are omitted for rather large estimation errors.It is clear the stabilities of the residuals basically remain unchanged over the long timescale.
To balance the data-span coverage and the amount of MSPs for the construction of the EPT, we select seven MSPs from the second data release of IPTA Version A. The seven MSPs have a common data span of 14.7 yr.The data input to WF are generated as follows: (i) Getting a precise timing model by fitting all the TOAs with a reference clock set to TT(BIPM15), which is the default reference timescale for the second data release of IPTA.
(ii) Setting the reference clock with which we are going to construct the EPT.If the reference clock is TT(BIPM15), this step is skipped.
(iii) Obtaining pre-fit timing residuals and selecting residuals in the longest common interval for all MSPs in use.
(iv) For the residuals selected above, using the average method over successive intervals of 30 days in length and placing each mean residual in the center of the corresponding interval.The successive intervals for each MSP are aligned in order.This step generates equispaced residuals that preserve the clock signal.
(v) If there are no residuals in one interval, piecewise linear interpolation is used to generate the mean residual from adjacent mean residuals points.
(vi) The quadratic trend is fitted and removed from each evenly spaced residual.This will make the residuals stationary, so WF can be used correctly.
The basic parameters of the observational data are listed in Table 2.

PSD Estimation by AR Model
Before searching for the upper limit of each DS, we need first to look into the problem of PSD estimation.As a PSD estimator, a periodogram has very large variance, which means the periodogram value may deviate from the true PSD greatly.If the upper limit is not precise in magnitude, the finally constructed EPT will not reveal the fluctuation of the reference clock correctly.Therefore, we should first confirm the validity of the PSD we are going to use.
The AR model is a widely used method for describing a class of SP.Compared with the moving average (MA) or autoregression moving average (ARMA) models, the solution of the AR parameters only involves solving linear equations.Apart from this fact, the MA and ARMA models can also be approximated by the AR model to sufficient order.Therefore, the AR model can provide a way of obtaining a PSD estimation.Generally, the expression of the AR model for the SP Y t is described as (He 2003) where A(z) is a complex polynomial, z is a complex number, { } e¢ t is a white-noise process WN(0, σ 2 ) with variance σ 2 , and a and p are the coefficient and order of the model.The PSD of Y t described by Equation ( 24) is Different from a periodogram, the PSD given by Equation (25) has a much higher resolution and a very smooth spectral curve.It is also independent of the PSD provided by the periodogram, so we can use it to help verify the upper limit.However, the AR model is not suitable for estimating the PSDs of signals with a strong single periodicity, such as a sinusoidal signal (Ding et al. 2007).In this case, the energy peak of the PSD would drift from the true value.In practice, we must watch carefully to avoid this case before deciding the upper limit.
The validity of the AR model can be verified by a whitenoise test (He 2003).From a χ 2 test, approximately follows a χ 2 distribution, where rj is the autocorrelation coefficient of the sample y 1 , y 2 ,...,y N1 and N1 is the number of sample points.y 1 , y 2 ,...,y N1 are produced by the model and observational data, so they should be independent of white noise.Therefore, from , there would be a large possibility that the model is incorrect.
By carefully examining the difference of the PSDs obtained by the two algorithms, we find for DS1, DS2, and DS3 the upper limit should be calculated by the AR model.Figure 8 depicts the χ 2 tests for DS1, DS2, and DS3.

EPT Construction and Discussion of Results
Similar to the case of simulated data, we first calculate the PSD of each DS and find the upper limit.The results show that all the upper limits are provided by MSP pair J0437-4715 and J1713 + 0747.Then the upper limits are used for PSD adjustment.Figure 9 depicts the adjusted PSDs and the upper limit for DS2.For each DS, the main energy region is defined as [ω lh , ω hh ], which is varying with the DS.After calculating S(e i ω ) in Equation (20), the energy outside the main energy region is suppressed by the upper limit.
To determine the weights of different MSPs, we add sine waves with amplitudes of 50-400 ns to DS1-DS5.The weights determined by Equation (21) are depicted in Figure 10.It is clear that as the signals increase, the differences of the weight values are reducing.However, for signals under 400 ns, the EPT is mainly influenced by J0437-4715 and J1713 + 0747.Here, we choose the weights by an intuitive and reasonable assumption: the signal amplitudes of TT(BIPM15) and TT(TAI) are 50 ns and 200 ns, and TT(TAI1) and TA(NTSC) are 400 ns.The weights for DS1 to DS6 are listed in Table 3.  4): the properties of the released data.Columns (5)-( 8): the weighted rms (wrms) of the residuals after the six pre-processing steps.(i) In the first line, the fluctuation levels for DS2 and DS3 are basically consistent with those of clock signals TT(TAI) and TT(TAI1), which means the pulsar pair determining the upper limit provides basically correct energy estimations.The variation trends of the EPT for DS2 and DS3 are consistent with TT(TAI) and TT(TAI1).It is difficult to say if there exist any systematical errors in TT(BIPM15), because the general noise level of the EPT for DS1 is larger than 50 ns.
(ii) In the second line, the EPT residuals of DS2 are more in accordance with EPT-TT(BIPM15).For DS3, as the clock signal doubles, the EPT residuals differ from EPT-TT(BIPM15) more significantly.The third line is another reflection of the changes of the EPT.We can see the EPT varies linearly at many points.
(iii) There exist some discrepancies between the EPT and the reference clock for DS2 and DS3 in the first line, and also for the EPT variations in the third line.The significant discrepancies of the EPT in the first line mainly concentrate on the beginning of the DS and the MJD interval 51526-52276.
As J0437-4715 and J1713 + 0747 play the key role, with a sum weight over 0.8 in EPT construction, we will conduct an inspection of the residuals to see what leads to these discrepancies.Figure 12 depicts the equispaced timing residuals (PT-TT(TAI)) of the two MSPs (with the quadratic trend fitted and removed) and the corresponding signals recovered from them by WF.At the beginning of the timing residuals, we can see clear bulge structures in both timing residuals.This kind of correlated noise is also seen in MJD interval 51526-52276.Therefore, in crosscorrelation and WA, this structure can hardly be canceled out.Now let us consider that the noise of the timing residuals of different pulsars is not uncorrelated.Denote the residuals x i and x j as in which ε is the common noise of the residuals x i and x j ; ε, ε i , and ε j are independent from each other; and ε i and ε j are zero mean.Then there exists ( )   would be significantly nonzero in most cases.Hence, the change of correlation from Equations (29) to (30) is no longer linear.This can explain the nonlinear variations of the EPT when the reference clock is changed.If the noise in the timing residuals is pure and low, as in the simulated data, a result like Figure 5 will be obtained.
The procedure proposed above can also be used to construct the EPT referred to other timescales, such as TA(NTSC).TA(NTSC) is the local atomic timescale maintained by NTSC. Figure 13 depicts the EPT constructed by DS4 and DS5.From the figure, we can see the EPT successfully indicates the variation trend of TA(NTSC), though the EPT residuals in the second line and EPT variations in the third line show some discrepancies.For DS4, we cannot give a confirmation of the existence of systematic error in TT(BIPM15).The fluctuations of EPT-TT(BIPM15) are generally larger than 50 ns.
Figure 14 depicts the stability of the EPT obtained above and of the corresponding atomic timescales.On a short scale, the stability of the EPT is corrupted by high-frequency noise, so it is not as good as TAI and TA(NTSC).The stability advantage of the EPT starts to appear after 4-5 yr.The instability of TAI can be slightly reflected on a scale of 7 yr, but the error bars of σ z are too large to make a precise assertion.In the right panel, it is hardly possible to get valid stability information of the EPT in EPT-TA(NTSC) on a scale of 6-7 yr.We must point out there may be some inaccuracies in the stability estimation by σ z here.If the variation rate of EPT-TA(NTSC) is nearly quadratic  in both the former and latter half-intervals, the last value of the σ z curve for EPT-TA(NTSC) will hardly reveal the EPT performance.
The most significant application of pulsar timing in astrophysics is searching for isotropic and stochastic gravitational-wave backgrounds (GWBs).The characteristic strain spectrum can be expressed as a power-law form where A is the dimensionless amplitude of the GWB, A log 14.7 1.9 0.9 at the 95% confidence level) given by the Chinese Pulsar Timing Array (CPTA; Xu et al. 2023) and -+ -

10
0.6 0.7 15 (at the 90% confidence level) given by the North American Nanohertz Observatory for Gravitational Waves (Agazie et al. 2023).In both panels of Figure 14, we show the theoretical curve of σ z induced by the GWB with h c ( f ) = 1.995 × 10 −15 at f = 1 yr −1 .Moreover, the upper bounds of h c ( f ) at f = 1 yr −1 with a 68% confidence level are 3.966 × 10 −15 and 6.547 × 10 −15 , respectively, given by σ z (τ = T span /2) of EPT-TT(BIPM15) for different data spans.It is clear that, due to this method, longer pulsar observational data would give more rigorous limits on GWBs.It is interesting to note that the value of σ z induced by the GWB estimated from the CPTA data is well below the upper bounds given by our results and is in accordance with our results at the 68% confidence level.It is worth pointing out that σ z cannot provide a precise detection for GWBs.A more realistic meaning is that σ z could be a tool to estimate a robust upper limit of the intensity of GWBs when the pulsar timing residuals are dominated by GWBs.To make the detection at higher confidence, one should find the quadripolar signals based on the spatial correlations between different pulsar pairs and fit the "Hellings and Downs (HD) curve."Tiburzi et al. (2016) pointed out that though the correlation of the common clock offset is different from the HD curve, this kind of noise with some others, such as the dipolar signal caused by the error of ephemeris and the fluctuations induced by solar wind, would also significantly affect GWB searching and lead to false detections.The improved WF is a good method to suppress the monopolar signals caused by clock errors, so it should potentially promote the research in GWB detection by PTAs.

Conclusion
In order to construct the EPT under the inclusion of timing residuals with high-level noise, we propose an improved method based on the WF algorithm in the frequency domain.This improved method can effectively reduce the energy and phase distortion in the EPT caused by high-level noise, so as to recover the common reference clock signal to the maximum extent.
The improved WF method successfully detects a common signal in simulated data under the setting of noise conditions.Based on the seven MSPs in the second data release of IPTA Version A, the EPT also indicates a global variation trend with TAI and TA(NTSC).For observational data, some obvious discrepancies exist in the EPT.This is due to the existence of correlated noise and the rather higher noise level in the timing residuals.It can be resolved by noise detection and elimination first.
From a probabilistic point of view, as the number of pulsars involved in EPT construction increases and the noise of the different residuals is truly uncorrelated, the effect of the noise from different MSPs will gradually be canceled out.But if most residuals keep relatively very strong noise, the weights and the effectiveness of these residuals would be very small.Thus, it is more effective to use random noise sources to detect and eliminate the existing noise in pulsar timing residuals so as to improve the current EPT level.This is a direction in EPT research.
If we denote the frequency spectrum of x(t) and x[kT] as X(ω) and X 1 (ω), respectively, and simplify å =-¥ +¥ k as ∑ k , we can get the expected signal PSD, because in the different timing residuals x a , x b , indicates that the cross-correlation function of x a , x b is equal to the autocorrelation function of s.

Figure 1 .
Figure1.Simulated data for testing the WF method.For SD3 and SD4, the error bars are too small and therefore fall within the residual marks.

Figure 3 .
Figure 3.The distortion of the signal phase under the interference of noise of different strengths in the complex plane.

Figure 4 .
Figure 4. Different PSDs type for determining the main energy region.

Figure 10 .
Figure 10.The changes of weight with respect to the amplitude of the sinusoidal wave.The residual used to obtain the weight is TT(BIPM15)-PT (MJD: 50206-55576).

Figure 11 .
Figure 11.EPT referred TT(BIPM15), TT(TAI) and artificial timescale TT(TAI1).Each column corresponding to a timescale.The first line shows the EPT and clock offset.The second line shows the EPT residuals.The third line is the EPT variations compared to TT(BIPM15)-reference clock.

Figure 13 .
Figure 13.EPT constructed by DS4 and DS5.The first line shows the EPT and clock offset.The second line shows the EPT residuals.The third line shows the EPT variations compared to the TT(BIPM15) reference clock.

Table 1
Noise Parameters of Simulated Data The last column lists the weight for each SDX in EPT construction.

Table 2
Data Used from the Second Data Release of IPTA Version A Figure9.Adjusted cross-PSDs, the PSD obtained by WA, and the upper limit for DS2.The upper limit is drawn by the black thick line.The cross-PSD determined by PSR J0437-4715 and J1713 + 0747 after adjustment is drawn by the red thick line.The PSD obtained by WA after adjustment is drawn by the blue thick line.The thin lines are 20 other cross-PSDs after adjustments.The main energy region is [0, 0.67](yr −1 ).

Table 3
Weights of Each MSP in EPT Construction (Verbiest et al. 2009)index of the GWB.Different origins generate GWBs with different values of α.As an example, we show the σ z (τ) induced by the GWB generated by supermassive black hole mergers in Figure14.The σ z (τ) caused by this type of GWB in the log-log plot has a slope of 2/3(Verbiest et al. 2009).The latest estimations of h c ( f ) with α = − 2/3 at f = 1 yr −1 are -