Multifractal Analysis of Pulsar Timing Residuals: Assessment of Gravitational Wave Detection

TL;DR

This paper develops a multifractal analysis pipeline, including a new irregular sampling algorithm, to detect gravitational wave signatures in pulsar timing residuals, revealing multifractal behavior and setting bounds on gravitational wave background amplitude.

Contribution

The paper introduces a novel multifractal analysis pipeline with an irregular sampling algorithm and a new cross-correlation function to analyze pulsar timing residuals for gravitational wave detection.

Findings

All pulsar data are nonstationary with universal features.

Observed data exhibit multifractal behavior due to data correlations.

Upper bound on GWB amplitude is $oxed{< 2.0 imes 10^{-15}}$.

Abstract

We introduce a pipeline including multifractal detrended cross-correlation analysis (MF-DXA) modified by either singular value decomposition or the adaptive method to examine the statistical properties of the pulsar timing residual ($PTR$) induced by a gravitational wave (GW) signal. We propose a new algorithm, the so-called irregular-MF-DXA, to deal with irregular data sampling. Inspired by the quadrupolar nature of the spatial cross-correlation function of a gravitational wave background, a new cross-correlation function, $\bar{\sigma}_{\times}$, derived from irregular-MF-DXA will be introduced. We show that, this measure reveals the quadrupolar signature in the $PTRs$ induced by stochastic GWB. We propose four strategies based on the $y$-intercept of fluctuation functions, the generalized Hurst exponent, and the width of the singularity spectrum to determine the dimensionless amplitude and power-law exponent of the characteristic strain spectrum as $\mathcal{H}_c(f)\sim\mathcal{A}_{yr}(f/f_{yr})^{\zeta}$ for stochastic GWB. Using the value of Hurst exponent, one can clarify the type of GWs. We apply our pipeline to explore 20 millisecond pulsars observed by Parkes Pulsar Timing Array. The computed scaling exponents confirm that all data are classified into a nonstationary class implying the universality feature. The value of the Hurst exponent is in the range $H\in [0.56,0.87]$. The $q$-dependency of the generalized Hurst exponent demonstrates that the observed $PTRs$ have multifractal behavior, and the source of this multifractality is mainly attributed to the correlation of data which is another universality of the observed datasets. Multifractal analysis of available $PTRs$ datasets reveals an upper bound on the dimensionless amplitude of the GWB, $\mathcal{A}_{yr}< 2.0\times 10^{-15}$.

Figures

Tables

Table 1. Table 1 Hurst exponent, H 𝐻 H , width of singularity spectrum, Δ ​ α Δ 𝛼 \Delta\alpha , scaling exponent of temporal autocorrelation, γ 𝛾 \gamma , rms, total time span (TTS) of post-fit timing residuals, and the upper limit on dimensionless amplitude of GWB of 20 MSPs observed in PPTA project. The error-bar corresponds to 1 ​ σ 1 𝜎 1\sigma confidence interval.

PSR Number	PSR Name	$H$	$\Delta \alpha$	$\gamma$	rms $\left(\mu s\right)$	TTS (yr)	${𝒜}_{yr}^{up}\left(95\%\right)$
1	J0437-4715	0.78 $\pm$ 0.03	0.89 $\pm$ 0.06	-1.56 $\pm$ 0.06	0.08	4.76	$5.0\times {10}^{-15}$
2	J0613-0200	0.68 $\pm$ 0.06	1.22 $\pm$ 0.04	-1.37 $\pm$ 0.11	1.07	5.99	$7.0\times {10}^{-15}$
3	J0711-6830	0.56 $\pm$ 0.10	1.40 $\pm$ 0.08	-1.13 $\pm$ 0.19	0.89	5.99	$6.0\times {10}^{-15}$
4	J1022+1001	0.65 $\pm$ 0.06	1.04 $\pm$ 0.04	-1.30 $\pm$ 0.13	1.72	5.88	$8.5\times {10}^{-15}$
5	J1024-0719	0.87 $\pm$ 0.03	1.60 $\pm$ 0.03	-1.74 $\pm$ 0.07	1.13	5.99	-
6	J1045-4509	0.84 $\pm$ 0.02	1.29 $\pm$ 0.04	-1.68 $\pm$ 0.05	2.77	5.94	-
7	J1600-3053	0.75 $\pm$ 0.05	1.34 $\pm$ 0.04	-1.50 $\pm$ 0.09	0.68	5.93	-
8	J1603-7202	0.68 $\pm$ 0.04	1.29 $\pm$ 0.05	-1.37 $\pm$ 0.07	2.14	5.99	$2.5\times {10}^{-15}$
9	J1643-1224	0.83 $\pm$ 0.04	0.89 $\pm$ 0.02	-1.66 $\pm$ 0.08	1.64	5.87	-
10	J1713+0747	0.74 $\pm$ 0.04	1.20 $\pm$ 0.05	-1.48 $\pm$ 0.09	0.31	5.71	$2.0\times {10}^{-15}$
11	J1730-2304	0.60 $\pm$ 0.11	1.79 $\pm$ 0.04	-1.21 $\pm$ 0.23	1.47	5.93	-
12	J1732-5049	0.81 $\pm$ 0.03	1.56 $\pm$ 0.03	-1.62 $\pm$ 0.07	2.22	5.08	$2.0\times {10}^{-15}$
13	J1744-1134	0.85 $\pm$ 0.04	1.52 $\pm$ 0.03	-1.70 $\pm$ 0.09	0.32	5.87	-
14	J1824-2452A	0.70 $\pm$ 0.03	1.26 $\pm$ 0.05	-1.40 $\pm$ 0.07	2.44	5.75	$10.0\times {10}^{-15}$
15	J1857+0943	0.71 $\pm$ 0.05	1.45 $\pm$ 0.02	-1.42 $\pm$ 0.10	0.84	5.93	-
16	J1909-3744	0.76 $\pm$ 0.06	1.32 $\pm$ 0.06	-1.52 $\pm$ 0.11	0.13	5.75	$6.0\times {10}^{-15}$
17	J1939+2134	0.80 $\pm$ 0.02	1.25 $\pm$ 0.02	-1.61 $\pm$ 0.04	0.68	5.88	-
18	J2124-3358	0.65 $\pm$ 0.07	1.23 $\pm$ 0.04	-1.30 $\pm$ 0.13	1.90	5.99	$6.0\times {10}^{-15}$
19	J2129-5721	0.66 $\pm$ 0.07	1.54 $\pm$ 0.04	-1.32 $\pm$ 0.13	0.80	5.86	$7.0\times {10}^{-15}$
20	J2145-0750	0.69 $\pm$ 0.06	1.29 $\pm$ 0.05	-1.38 $\pm$ 0.11	0.78	5.99	-

Equations

$$PT{R}_a(i),PT{R}_b(i),i=1,...,N$$

$$X_{\diamond }(j)$$

$${\mathcal{E}}_{\times }(s,\nu )=\frac{1}{s}\sum _{i=1}^s$$

$$\widetilde{X_{\diamond }(j)}=\frac{1}{s}\sum _{k=-s_1}^{s_2}{X}_{\diamond }(j-k)$$

$${\epsilon }_{X_{\diamond }}(i)=X_{\diamond }(i)-\widetilde{X_{\diamond }(i)}$$

$${\mathcal{E}}_{\times }(s,\nu )=\frac{1}{s}\sum _{i=1}^s{\epsilon }_{X_a}(i+(\nu -1)s)\times {\epsilon }_{X_b}(i+(\nu -1)s)$$

$${\mathcal{F}}_{\times }(q,s)={\left(\frac{1}{2N_s}\sum _{\nu =1}^{2N_s}{\left|{\mathcal{E}}_{\times }(s,\nu )\right|}^{q/2}\right)}^{1/q}$$

$${\mathcal{F}}_{\times }(0,s)=\exp \left(\frac{1}{4N_s}\sum _{\nu =1}^{2N_s}\ln |{\mathcal{E}}_{\times }(s,\nu )|\right)$$

$${\mathcal{F}}_{\times }(q,s)\sim s^{h_{\times }(q)}$$

$${\mathcal{F}}_q(s)={\mathcal{G}}_{h(q)}{s}^{h(q)}$$

$$\mathcal{G}$$

$$\xi (q)=qh(q)-1$$

$$D(q)=\frac{\xi (q)}{q-1}=\frac{qh(q)-1}{q-1}$$

$$f(\alpha )=\alpha q-\xi (q)$$

$${\alpha }_{\min }=\underset{q\to +\infty }{\lim }\frac{\partial \xi (q)}{\partial q},{\alpha }_{\max }=\underset{q\to -\infty }{\lim }\frac{\partial \xi (q)}{\partial q}$$

$${\sigma }_{\times }({\Theta }_{ab})\equiv {\sum }_s\left(\frac{{\sum }_{\nu =1}^{2N_s}{\mathcal{E}}_{\times }(s,\nu )}{\sqrt{\left[{\sum }_{\nu =1}^{2N_s}{\mathcal{E}}_a(s,\nu )\right]\left[{\sum }_{\nu =1}^{2N_s}{\mathcal{E}}_b(s,\nu )\right]}}\right)$$

$${\bar{\sigma }}_{\times }(\Theta )=\frac{1}{4\pi }\int d\Omega {\sigma }_{\times }({\Theta }_{ab})$$

$${\mathcal{C}}_{\times }({\Theta }_{ab})$$

$${\overline{\mathcal{C}}}_{\times }(\Theta )$$

$$\overline{\Gamma }(\Theta )=\frac{3}{2}\psi \ln (\psi )-\frac{\psi }{4}+\frac{1}{2}$$

$$PT{R}_{reg}(t)=\int d{t}^{\prime }PT{R}_{irre}(t^{\prime })\mathcal{W}(t-t^{\prime })$$

$${\mathcal{E}}^2(s,\nu )=\frac{1}{s_{\nu }^{\prime }(s)}\sum _{i=1}^{s_{\nu }^{\prime }(s)}{\left[X(i+(\nu -1)s^{\prime })-{\tilde{X}}_{\nu }(i)\right]}^2$$

$${\mathcal{F}}_q(s)={\left(\frac{{\sum }_{\nu =1}^{2N_s}\frac{{\left[{\mathcal{E}}^2(s,\nu )\right]}^{q/2}}{{\sigma }_{\mathcal{E}}^2(s,\nu ,q)}}{{\sum }_{\nu =1}^{2N_s}\frac{1}{{\sigma }_{\mathcal{E}}^2(s,\nu ,q)}}\right)}^{1/q}$$

$$\mathbf{\Gamma}\equiv\left(\begin{array}[]{cccc}{\it PTR}_{1}&{\it PTR}_{1+\tau}&...&{\it PTR}_{1+N-(d-1)\tau-1}\\ \vdots&\vdots&\vdots&\vdots\\ {\it PTR}_{i}&{\it PTR}_{i+\tau}&...&{\it PTR}_{i+N-(d-1)\tau-1}\\ \vdots&\vdots&\vdots&\vdots\\ {\it PTR}_{d}&{\it PTR}_{d+\tau}&...&{\it PTR}_{d+N-(d-1)\tau-1}\\ \end{array}\right)$$

$$\mathbf{\Gamma }=\mathbf{US}{\mathbf{V}}^{†}$$

$${\widetilde{PTR}}_{i+j-1}={\tilde{\Gamma }}_{ij}.$$

$${\mathcal{Y}}_{\nu }^{\mathrm{overlap}}(j)=\left(1-\frac{j-1}{n}\right){\mathcal{Y}}_{\nu }(j+n)+\frac{j-1}{n}{\mathcal{Y}}_{\nu +1}(j)$$

$$\mathcal{P}({\rm Y}|\mathcal{D})=\frac{\mathcal{L}(\mathcal{D}|{\rm Y})\mathcal{P}({\rm Y})}{\int \mathcal{L}(\mathcal{D}|{\rm Y})\mathcal{P}({\rm Y})d{\rm Y}}$$

$${\chi }^2({\rm Y})\equiv {\Delta }^{†}.C^{-1}.\Delta$$

Full text

1. Introduction

Pulsar timing has received extensive attention for astrophysical interests due to possessing a stable rotational mechanism (Verbiest et al., 2015; Lorimer & Kramer, 2005). The pulsar timing residual ($PTR$) which is an important observable, is defined by the difference between the measured time of arrival (TOA) and those anticipated by a timing model (Verbiest et al., 2015). The observed $PTR$ is a precise indicator to elucidate some interesting physical properties of pulsars and other cosmological and astrophysical foreground processes (Manchester et al., 2013). The influence of unknown physical phenomena on the variation of the pulsar’s spin and spin-down and the presence of foreground effects impose randomness on the $PTR$. Therefore, the $PTR$ is categorized in a $(1+1)$-dimensional stochastic process (where only one of the parameters is independent, while the other parameter is represented as a function of the independent parameter). Therefore, the stochastic nature of the data requires implying robust methods to extract reliable information from the $PTR$.

Millisecond pulsars (MSPs) were first suggested as detectors of gravitational waves (GWs) by Sazhin (1978) and Detweiler (1979) because of the high stability and predictability of their rotational behavior (see also (Hobbs et al., 2009)). Indeed, GWs can be produced by different mechanisms ranging from the early epoch to the present era. Continuous wave sources (Peters, 1964), burst sources (Thorne & Braginskii, 1976; Damour & Vilenkin, 2001; Kocsis et al., 2006) and stochastic backgrounds (Enoki, 2007; Maggiore, 2000) are the most well-known classes among the GW sources. As an illustration, we refer to relic GWs, including GWs by cosmic strings and primordial perturbations (Maggiore, 2000; Hobbs, 2011; Damour & Vilenkin, 2005; Pshirkov & Tuntsov, 2010). The GWs are also produced in the formation of supermassive black holes and binary black hole mergers (Rajagopal, M., & Romani, 1995; Taylor & Gair, 2012). Recently, the GWs of black hole mergers have been detected by LIGO instruments, which can be an evidence of dark matter in the early universe or can correspond to the binary black hole of stellar origin (Mandic et al., 2016). Other classes of GWs include the continuous, inspiral, burst, and stochastic types of GWs (Meadors, 2014; Coyne & Owen, 2016). Many approaches have been proposed and utilized during past decades to detect mentioned types of GWs (Coyne & Owen, 2016; Pai et al., 2001; Zhu et al., 2014; Jenet et al., 2005). The very low amplitude of GWs, different sources and mechanisms for GW production on one hand, and the extended range of frequency on the other hand lead to introducing various indirect approaches such as predictions of energy loss due to GW emission (Taylor, 1994) and direct approaches such as detecting the effect of GWs on pulsar timing residuals (Jenet et al., 2005). The two main methods for detection of GWs are known as interferometers (such as LISA and LIGO) and pulsar timing arrays (Roebber & Holder, 2017).

For the frequency interval $\nu\in[10^{-8},10^{-6}]$, there are several pulsar timing array projects that observe the imprint of GWs using pulsar timing detectors (Jenet et al., 2006). In the context of pulsar timing array approach, some famous projects have been proposed, namely the Parkes Pulsar Timing Array (PPTA) (Manchester et al., 2013; Manchester, 2008; Hobbs, 2013), the European Pulsar Timing Array (EPTA) (Janssen et al., 2008; Kramer & Champion, 2013), the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) (Demorest & Lommen, 2009; McLaughlin, 2013) and International Pulsar Timing Array (IPTA) (Verbiest et al., 2015). The Square Kilometer Array (SKA) (Cordes et al., 2004; Lazio, 2013) radio telescope would further improve the sensitivity of pulsar timing measurements to detect GWs. For a recent and more complete discussions on various experiments and methods to detect GWs, see Zhu et al. (2015, 2016); Lommen (2015); Ellis (2014); Manchester (2010); George & Huerta (2017) .

A pulsar timing array utilizing the Parkes radio telescope in Australia is an experiment to detect GWs by observing 20 bright MSPs (Manchester et al., 2013). The PPTA observations must be continued more than 5 yr in order to achieve a precision of 100 ns. Since pulsar timing residuals are induced by GWs (Jenet et al., 2005), therefore some authors used the statistical correlation of pulsars timing residuals to evaluate their capability of detecting GWs (Hellings & Downs, 1983; Jenet et al., 2005).

In order to elaborate the benchmark of different types of GWs based on $PTR$ datasets, we should consider four aspects. First of all, various sources of GWs encourage us to find deep insight regarding the properties of GWs emitted by different sources. Second, we should examine the performance of different statistical methods and their sensitivities. Third, we should take into account the quadrupolar signature of GWs. Finally, the upper limit on the amplitude of the gravitational wave background (GWB) should be computed. The main method to detect a stochastic GW employing pulsar timing arrays is to search for a correlation between $PTR$s and compare it with quadrupole spatial cross-correlation calculated by Hellings and Downs (Hellings & Downs, 1983; Jenet et al., 2005; Van Haasteren et al., 2009, 2011). If such a signature is not detected, one can set an upper bound on the amplitude of GWs using frequentist or Bayesian approaches. The upper limits provided by Van Haasteren et al. (2011) and Shannon et al. (2015) for a stochastic background produced by supermassive black hole binaries is $\mathcal{A}_{yr}\leq 6\times 10^{-15}$ and $\mathcal{A}_{yr}\leq 1\times 10^{-15}$, respectively, where the latter is the lowest claimed upper limit so far. In addition, Demorest et al. (2012) have used the high-precision pulsar timing data recorded as part of the NANOGrav project and finally provided an upper limit on the power-spectrum amplitude of the nHz-frequency stochastic supermassive black hole GWB.

Pulsar datasets are manipulated by trends and noises. Statistical models for noises, trends, and signals play crucial roles in any parametric GW detection approaches. Subsequently, it is necessary to implement robust and novel methods for removing destructive effects from desired parts of signals.

Our work in this paper has the following advantages and novelties.

i) Inspired by the properties of a self-similar process characterized by a scaling exponent called ”Hurst exponent” (Hurst, 1951; Tessier et al., 1996; Pandey et al., 1998; Pelletier & Turcotte, 1997), for the first time, we have used Multifractal Detrended Fluctuation Analysis (MF-DFA) (Peng et al., 1995; Kantelhardt et al., 2002), Multi-Fractal Detrended Moving Average Analysis (MF-DMA) (Alessio et al., 2002; Carbone et al.,, 2004; Arianos et al.,, 2007; Gu & Zhou, 2010; Shao et al., 2015) and Multifractal Detrended Cross-correlation Analysis (MF-DXA) (Podobnik & Stanley, 2008; Zhou, 2008) methods to analyze the observed (including 20 MSPs inferred from (Manchester et al., 2013)) and simulated pulsar timing residuals induced by GW signals (simulated by the TEMPO2 software package Hobbs et al. (2006)). We will evaluate the multiscaling behavior of the underlying data from statistical point of view.

ii) We modify MF-DFA, MF-DMA and MF-DXA by additional denoising algorithms, namely Adaptive Detrending (AD) (Hu et al., 2009) or Singular Value Decomposition (SVD) (Golub & Van Loan, 1996; Nagarajan & Kavasseri, 2005c, a) methods to exclude or at least to reduce the contribution of unknown trends and noises as much as possible. These methods are used as precomplementary denoising procedures.

iii) The standard version of multifractal analysis is a reliable algorithm when the input is a regular sampling series. Observed $PTR$s are unevenly sampled data sets; we therefore modify parts of the MF-DXA algorithm and call it the irregular MF-DXA method. In addition, noise modeling can be revealed by multifractal analysis.

iv) We check the multifractal nature of pulsar timing residuals. We also determine the sources of multifractality based on our statistical approaches.

v) According to the quadrupolar signature on the spatial cross-correlation function of PTRs, the detectability of the stochastic GWB is evaluated according to the MF-DXA of PTRs. The cross-correlation exponent will be determined. We also give a new spatial cross-correlation function for pulsar timing residues.

vi) We introduce some criteria not only for discrimination of the stochastic GWB footprint and single sources GWs on pulsar timing residuals but also for determining the dimensionless amplitude of GWB. Then, an upper bound on the amplitude of stochastic GWs will be computed. Finally, this view provides a new insight to use pulsar timing residuals for further astrophysical studies.

The rest of this paper is organized as follows. In section 2, we will explain MF-DFA, MF-DMA, and MF-DXA, dealing with irregular sampled data, AD, and SVD in detail. A new measure for the spatial cross-correlation of $PTR$ is presented in this section. Noise and trend modeling and posterior analysis to obtain scaling exponents are also discussed in section 2. Section 3 is devoted to the theoretical notions of the GWB and data description for observed as well as synthetic datasets. We will implement the multifractal methods on simulated timing residuals series in section 4. We will also study a new spatial cross-correlation function derived by the MF-DXA method in the search of the footprint of the stochastic GWB in the sensitive range of pulsar timing residual series. Four strategies to reveal the imprint of GWs on the residual time series in a noiseless observation will also be explained in section 4. We will examine the multifractality of the observed pulsar timing residuals in section 5. In that section, we will also look for an upper bound on the amplitude of stochastic GWs using some observed $PTR$s. Section 6 is devoted to summary and conclusion.

2. Methodology: Multifractal analysis

Nonstationary sources such as trends and artificial noises usually influence the observed time series. To infer reliable results, these spurious effects should be well characterized and distinguished from the intrinsic fluctuations. Concerning trends, Wu et al. (2007) stated that, in principle, there is no universal definition for trends, and any proper algorithm for evoking trends from underlying series should remove the contribution of trends without destroying fluctuations. One of the well-studied methods for this purpose is MF-DFA (Peng et al., 1995; Kantelhardt et al., 2002), used in various areas, such as economical time series (Mantegna & Stanley, 2000; Liu & Stanley, 1999; Vandewalle et al., 1999; Ivanov et al., 2004; Ferreira et al., 2017), river flow (Movahed & Hermanis, 2008; Hajian & Movahed, 2010), sunspot fluctuations (Movahed et al., 2006; Hu et al., 2009), cosmic microwave background radiations (Movahed et al., 2011), music (Jafariet al., 2007; Jennings et al., 2004), plasma fluctuations (Kimiagar et al., 2009), identification of a defective single layer in two-Dimensional material (Shidpour and Movahed, 2018), traffic jamming (Xiao-Yan et al., 2007), image processing, medical measurements (Soares et al., 2009, 2010), and astronomy (Zunino et al., 2014). Cross-correlation has also been introduced and applied in some cases (Podobnik & Stanley, 2008; Podobnik et al., 2009, 2011; Qian et al., 2015; Zebende, 2011; Zebende et al., 2013; Kristoufek, 2015). The MF-DXA examining higher-order detrended covariance was introduced by Zhou (2008). Although the approaches in multifractal detrended analysis, such as the MF-DFA and MF-DXA methods, diminish polynomial trends, previous researche demonstrated that sinusoidal and power-law trends are not completely removed (Hu et al., 2001; Chen et al., 2002). Mentioned trends make some crossovers in fluctuation functions (Hu et al., 2001; Chen et al., 2002; Kantelhardt et al., 2001; Nagarajan & Kavasseri, 2005a, b, c). Several robust methods have been proposed to eliminate cross-overs produced by sinusoidal and power-law trends: Fourier Detrended Fluctuations Analysis (F-DFA) (Chianca et al., 2005; Nagarajan & Kavasseri, 2005b), Singular Value Decomposition (SVD) (Golub & Van Loan, 1996; Nagarajan & Kavasseri, 2005c, a), Adaptive Detrending method (AD) (Hu et al., 2009), and Empirical Mode Decomposition (EMD)(Huang et al., 1998). In this work, we implement the AD and SVD methods to reduce the contribution of noise and magnify the effect of GWs in our results for further cleaning preprocessors.

2.1. Multifractal-based analysis

Finding scaling exponents in the context of auto-correlation and cross-correlation analysis has many inaccuracies due to nonstationarity, noises, and undesired trends. To resolve the mentioned difficulties, a well-known method based on decomposing the original signal into its positive and negative fluctuation components has been proposed by Jun et al. (2006). Motivated by such a decomposition method, Podobnik et al. introduced the cross-correlation between two non-stationary fluctuations by means of the DFA method (Podobnik & Stanley, 2008). A modification of detrended cross-correlation analysis (DCCA) is known as MF-DXA was invented by Zhou (2008). The pipeline of MF-DXA is considered as follows (Podobnik & Stanley, 2008; Zhou, 2008)111If both signals are identical, we have the MF-DFA/MF-DMA method..

$(1)$: We consider two typical $PTR$ series named by $PTR_{a}$ and $PTR_{b}$, located at $\hat{n}_{a}$ and $\hat{n}_{b}$ with respect to the line of sight, respectively, as the input data sets to study their mutual multifractal property:

[TABLE]

The pulsar timing observations are almost unevenly sampled datasets. We need equidistant sampling series. A trivial but not essentially optimum way is to interpolate between two successive data. Different methods to reconstruct regular series will be explained in subsection 2.2. Therefore, here we assume that the input data are regular and ready for further tasks. Moreover, the observed data have variable error bars, and, to take into account heteroskedasticity, we use error-propagation formalism in all statistical analysis, such as averaging, fitting, and computing fluctuation functions throughout this paper.

$(2)$: To magnify the hidden self-similarity property, we make profile series according to:

[TABLE]

Here the subscript $\diamond$ can be replaced by ”$a$” or ”$b$”.

$(3-a)$: The above profile series must be divided into $N_{s}=\mathrm{int}(N/s)$ nonoverlapping segments of length $s$. The range of nonoverlapping window values is $N_{s}\in[N_{s}^{\rm min},N_{s}^{\rm max}]$. To take into account the remaining unused part of the data from the opposite end of the data, the enumeration must to be repeated from the mentioned part. In this case, we will have $2N_{s}$ segments. In the framework of the MF-DCCA method, we should compute the following fluctuation function in each segment as follows:

[TABLE]

for segments $\nu=1,...,N_{s}$. For the opposite end, we have:

[TABLE]

where $\nu=N_{s}+1,\cdot\cdot\cdot,2N_{s}$ and $\tilde{X}_{\diamond}^{(\nu)}(i)$ is a weighted fitting polynomial function in the $\nu$th segment with an arbitrary order describing the local trend for data with variable error bars. Usually a linear function for modeling local trends is considered (Bunde et al., 2000). The MF-DCCA$m$ denotes that the order of the polynomial function used in the MF-DCCA is ”$m$”. Throughout this paper, we take $m=1$ unless stated otherwise. To reduce the statistical uncertainties in the computed fluctuation functions, we set $s>m+2$ (Kantelhardt et al., 2002). On the other hand, this method becomes unreliable for very large window sizes, i.e. $s>\frac{N}{4}$. There is a discontinuity for fitting a polynomial at the boundary of each partition in the MF-DCCA method; to resolve this discrepancy, MF-DMA has been introduced (Alessio et al., 2002; Carbone et al.,, 2004; Arianos et al.,, 2007; Gu & Zhou, 2010; Shao et al., 2015). Accordingly, instead of doing item $(3-a)$, we carry out the following procedure:

$(3-b)$: For each moving window with size $s$, we calculate the moving average function:

[TABLE]

where $s_{1}=\lfloor(s-1)\theta\rfloor$ and $s_{2}=\lceil(s-1)(1-\theta)\rceil$. The symbol $\lfloor a\rfloor$ represents the largest integer value not greater than $a$ and $\lceil a\rceil$ is devoted to the smallest integer value not smaller than $a$. In the above equation, $\theta$ plays a crucial role. The $\theta=0$ refers to the backward moving average, while $\theta=1$ is the so-called forward moving average; finally $\theta=0.5$ is related to the centered moving average (Xu et al., 2005; Gu & Zhou, 2010). Therefore, detrended data are constructed by subtracting the calculated moving average function from the cumulative series, $X_{\diamond}$ as:

[TABLE]

where $s-s_{1}\leq i\leq N-s_{1}$. Now $\varepsilon_{X_{\diamond}}(i)$ values are divided into $N_{s}={\rm int}[N/s]$ nonoverlapping windows with the same size of $s$ and we calculate the fluctuation function:

[TABLE]

$(4)$: Using Eqs. (3) and (2.1) for the MF-DCCA (MF-DFA) based method (Peng et al., 1992, 1994; Buldyrev et al., 1995; Kantelhardt et al., 2002; Shao et al., 2012) and Eq. (7) for the MF-DMA algorithm, the corresponding $q$th-order fluctuation function can be computed by:

[TABLE]

For $q=0$, we have:

[TABLE]

$(5)$: The scaling behavior of the fluctuation function according to:

[TABLE]

gives the cross-correlation exponent $h_{\times}(q)$. The $q$-parameter enables us to quantify the contribution of different values of fluctuation functions in Eqs. (8) and (9). The small fluctuations play a major role in summation for $q<1$, while large fluctuations become dominant for $q\geq 1$. We emphasize that for heteroskedastic data, the summation in Eqs. (8) and (9) should incorporate variable errorbars, and weighted fitting polynomials must be considered. It turns out that for $a=b$, the usual generalized Hurst exponent, $h(q)$, is retrieved. In this case we have:

[TABLE]

for $q=2$, the $\mathcal{G}$ is

[TABLE]

and $\sigma^{2}=\langle PTR^{2}\rangle$ for zero mean data. Any $q$-dependency of $h(q)$, confirms that the underlying data set is a multifractal process. For the class of the nonstationary series (corresponding to a fractional Brownian motion; fBm) the exponent derived by using MF-DFA is $h(q=2)>1$. Therefore, in this case, the Hurst exponent is given by $H=h(q=2)-1$. In the stationary case, $h(q=2)<1$ (corresponding to a fractional Gaussian noise; fGn) and $H=h(q=2)$. For completely stationary random data, $H=0.5$, while for a persistent data set, $0.5<H<1.0$. For an anticorrelated data set, $H<0.5$ (Taqqu et al., 1995; Peng et al., 1994; Ossadnik et al., 1994). When the Hurst exponent is determined, the scaling exponents of autocorrelation for an fGn process read as $\mathcal{C}(\tau)=\langle x(t)x(t+\tau)\rangle\sim\tau^{-\gamma}$ for $\tau\gg 0$ with $\gamma=2-2H$, while for a fBm signal, we have $\mathcal{C}(t_{i},t_{j})=\langle x(t_{i})x(t_{j})\rangle\sim t_{i}^{-\gamma}+t_{j}^{-\gamma}-|t_{i}-t_{j}|^{-\gamma}$ for $|t_{i}-t_{j}|\gg 0$ with $\gamma=-2H$. The associated power spectrum is $S(f)\sim f^{-\beta}$ with $\beta=2H-1$ and $\beta=2H+1$ for the fGn and fBm processes, respectively. The relation between the generalized Hurst exponent and the scaling exponent of the partition function known as the multifractal scaling exponent based on the standard multifractal formalism becomes (Kantelhardt et al., 2002):

[TABLE]

For a monofractal data set, $\xi(q)$ is a linear function (Kantelhardt et al., 2002). The generalized multifractal dimension is also given by:

[TABLE]

where $D(q=0)=D_{f}$ is the fractal dimension of the time series and $D(q=1)$ is related to the so-called entropy of the underlying system (Halsey et al., 1986). A more complete quantitative measure of multifractality is the singularity spectrum and indicates how the box probability of standard multifractal formalism behaves at small scales. It is defined by the Legendre transformation of $\xi(q)$ as (Feder, 1989):

[TABLE]

and the H${\rm\ddot{o}}$lder exponent is $\alpha\equiv d\xi(q)/dq$. In the case of multifractality, a spectrum of the H${\rm\ddot{o}}$lder exponent is obtained instead of a single exponent. The domain of the H${\rm\ddot{o}}$lder spectrum, $\alpha\in[\alpha_{\rm min},\alpha_{\rm max}]$, becomes (Muzy et al., 1994; Arneodo et al., 1995):

[TABLE]

Subsequently, the width $\Delta\alpha\equiv\alpha_{\rm max}-\alpha_{\rm min}$ is a reliable measure for quantifying the multifractal nature of the underlying data. The higher value of $\Delta\alpha$ is associated with the higher multifractal nature reflecting the complexity of the signal. As other complexity measures, one can point to the $q$-order Lyapunov exponent (Eckmann & Procaccia, 1986), and the Lempel-Ziv complexity (Lempel & Ziv, 1976). Inspired by the common cross-correlation definition, relying on Eq. (7), we define the new cross-correlation function (Zebende, 2011; Zebende et al., 2013):

[TABLE]

here $\Theta_{ab}=\arccos|\hat{n}_{a}.\hat{n}_{b}|$. Averaging on all available pairs separated by $\Theta$ leads to:

[TABLE]

The $\bar{\sigma}_{\times}$ introduced by Eq. (18) based on fluctuation functions computed in the context of detrended cross-correlation contains the quadrupolar signature if $PTR$s are modified by the GWB signal. Therefore, this is a new criterion that enables us to assess the footprint of GWs more precisely.

Now we turn to the spatial cross-correlation function for $PTR$s taking into account stationarity as:

[TABLE]

In the presence of an isotropic GWB, by averaging the cross-correlation on all available pairs separated by $\Theta$ leads to:

[TABLE]

The $\overline{\Gamma}(\Theta)$ is given by the Hellings and Downs equation (Hellings & Downs, 1983; Jenet et al., 2005):

[TABLE]

where $\psi\equiv{[1-\cos(\Theta)]}/{2}$. We should notice that the Hellings and Downs curve is only a function of the angular separation between pulsar pairs separated by $\Theta$, and it is independent of the frequency (Romano & Cornish, 2017).

The new cross-correlation coefficient defined by Eq. (18) is related to the traditional cross-correlation $\mathcal{C}_{\times}$ in a complex way, and the relation is not analytically tractable without any approximation, and we will evaluate it numerically in the next section. However, according to Eq. (7), the mapping between $\mathcal{C}_{\times}$ and $\sigma_{\times}$ does not change the sign of $\sigma_{\times}$. Thus, the quadrupolar signature of the Hellings and Downs function is preserved. It is worth mentioning that, besides probable GW signal superimposed in the $PTR$s, the following fluctuations can be existed in the recorded data: the correlated red (fractal) noise; clock errors, which are the same in all pulsars (i.e., monopolar); and ephemeris errors (which are dipolar). There are no known noise sources other than GWs that are quadrupolar (Tiburzi et al., 2015).

Applying MF-DXA on $PTR$s determines the value of the temporal scaling exponent, $h_{\times}$. We expect to find constant $h_{\times}(q)$ with respect to different separation angles ($\Theta$) for an isotropic GWB, while for the other local source of GWs, the $h_{\times}(q)$ depends on $\Theta_{ab}$ in an arbitrary manner.

2.2. Dealing with irregularly sampled data

The pulsar timing observations are unevenly sampled; i.e. they are not a set of equidistant sampling values, and the underlying series is nonuniform, requiring some sort of interpolation technique. The Lomb-Scargle periodogram proposed a least-squares pipeline to resolve this problem (Lomb, 1976; Scargle, 1982). Radon transformations have also been used for irregular sampling analysis (Ronen et al., 1991; Duijndam & Schonewille, 1999; Duijndam et al., 1999); see also (Gulati & Ferguson, 2009) and references therein. Extrapolation of irregularly recorded data onto a regular grid was introduced by Ferguson (2006). For constructing Fourier expansion, nonuniform discrete Fourier transform was introduced by Gulati & Ferguson (2009); Anholm et al. (2009). A trivial but not necessarily optimum method with less computational burden is to interpolate between two successive data points in recorded series. A more robust method is to apply kernel functions on the irregular data, as (see also (Monaghan & Lattanzio, 1985)):

[TABLE]

where $PTR_{reg}$ and $PTR_{irre}$ are regular and nonuniform sampled data, respectively. Here $\mathcal{W}$ is a normalized window function. A typical functional form for this window function can be Gaussian. In general, the choice of the window function, $\mathcal{W}$, depends on the smoothness, accuracy requirements, and computation efficiency (Monaghan & Lattanzio, 1985).

Here we propose a new approach to find robust scaling properties for irregular sampled data. If there is no a priori information for the smoothing procedure, we suggest applying a gaussian kernel to the data followed by a linear interpolation to regularize datasets. Subsequently, we can construct the profile using such regular data (Eq. (2)). To reduce the contribution of artificial data points produced in this interpolation, we introduce tge irregular MF-DXA method. In this new algorithm, we modify the fluctuation function procedure given by Eqs. (3) and (2.1) for identical $PTR$s as:

[TABLE]

In the above equation, only the data points recorded during observation in each segment with size $s$ will be considered for further computations. Therefore, the number of data in the $\nu$th window with size $s$ is represented by $s_{\nu}^{\prime}(s)$ which in general is not equal to $s$. Now Eq. (8) becomes a weighted average:

[TABLE]

where $\sigma_{\mathcal{E}}^{2}(s,\nu,q)$ is the variance of $\left[{\mathcal{E}}^{2}(s,\nu)\right]^{q/2}$. We similarly replace the averaging procedure in any relevant parts with the weighted averaging.

Recently, Ma et al. (2010) showed that the global scaling exponents of long-correlated signals remain unchanged for up to 90% of data loss, while for anticorrelated series, even less than 10% of data loss creates a significant modification in the original scaling exponents. This research shows that one can compute the scaling exponents for long-range correlated irregularly sampled data points if one regularizes the data set through linear interpolation and then applies DFA. But for an anticorrelated signal, the DFA method does not lead to reasonable results. Our new proposal demonstrates that for synthetic series with known Hurst exponents, our modification leads to more reliable estimations for scaling exponents, not only for correlated series but also for anticorrelated datasets. Our simulations show that the $PTR$ can be considered as long-range correlated fluctuation. Therefore, our results are almost are not affected by the type of regularization.

2.3. SVD

It is important to find trends and noise sectors in data analysis, especially in the astronomical data. When we use MF-DFA, MF-DMA, and MF-DXA, an essential demand corresponding to presenting a scaling behavior must be satisfied, as represented by Eqs. (10) and (11). In some cases, there exist one or more crossovers corresponding to different correlation behaviors of the pattern in various scales (Hu et al., 2001; Chen et al., 2002; Kantelhardt et al., 2001; Nagarajan & Kavasseri, 2005a, b, c). The MF-DFA and MF-DXA methods cannot remove the effect of all undesired parts of the underlying signal; therefore, we implement complementary tasks to properly recover the scaling behavior of fluctuation functions properly and to obtain the reliable scaling exponents. There are some preprocessing methods for denoising in the literature; for instance, The EMD method (Huang et al., 1998), the Fourier-detrended (Fourier-based filtering) method (Chianca et al., 2005; Nagarajan & Kavasseri, 2005b), the SVD method (Golub & Van Loan, 1996; Nagarajan & Kavasseri, 2005c, a) and the AD algorithm (Hu et al., 2009). In this paper, we utilize the SVD method and AD algorithm. The main part of the SVD method can be described in the following steps (Nagarajan & Kavasseri, 2005a, c; Hajian & Movahed, 2010):

(I): Construct a matrix whose elements are $PTR$s in the following order:

[TABLE]

where $d$ is the embedding dimension, $\tau$ is the time delay, and $1\leq i\leq d$. Considering a time series of size $N$, the maximum value of the embedding dimension $d$ is equal to $d\leq N-(d-1)\tau+1$ (Nagarajan & Kavasseri, 2005c, b; Shang et al., 2009).

(II): Decompose the matrix $\mathbf{\Gamma}$ to left (${\mathbf{U}}_{d\times d}$) and right (${\mathbf{V}}_{(N-(d-1)\tau)\times(N-(d-1)\tau)}$) orthogonal matrices:

[TABLE]

where ${\mathbf{S}}_{d\times(N-(d-1)\tau)}$ is a diagonal matrix and its elements are the desired singular values. If we are interested in examining the fluctuations with high frequency, we should remove dominant wavelengths. In this case, for removing trends containing $p$ dominant wavelengths, we set $2p+1$ largest eigenvalues of matrix $\mathbf{S}$ to zero; therefore, long periods or short frequencies are eliminated. In other words, $p$ dominant eigenvalues and associated eigenvectors correspond to long wavelength (short-frequency part) subspace, while $d-p$ eigenvalues and the corresponding eigen-decomposed vectors represent short-wavelength (high-frequency part) subspace.

In this paper, we look for the footprint of GWs superimposed on the $PTR$s signals. As shown in Fig. 1, the GW part behaves as a dominant trend in $PTR$s; consequently, we essentially need to do denoising using the SVD method to magnify the contribution of superimposed GWs. To this end, we should remove small eigenvalues corresponding to a low-pass filter. In this paper, we eliminate the high-frequency part of the signal by keeping the $2p+1$ largest eigenvalues of the matrix $\mathbf{{S}}$.

Finally, the new eigenvalues matrix, $\mathbf{\tilde{S}}$, is determined. According to the filtered matrix, $\mathbf{\tilde{\Gamma}}={\mathbf{U\tilde{S}V^{\dagger}}}$, the cleaned time series is constructed by:

[TABLE]

Here $1\leq i\leq d$ and $1\leq j\leq N-(d-1)\tau$. Now the cleaned $\widetilde{{\it PTR}}$ datasets will be used as input for the MF-DFA or MF-DXA discussed in previous subsections.

2.4. AD algorithm

Another robust algorithm to examine trends is the AD method introduced by Hu et al. (2009). The implementation of the AD algorithm is a complementary method for determining local and global trends. Therefore, after applying the AD method on observed pulsar timing series, the corresponding dominant trend output data will be used as an input for the MF-DFA or MF-DXA methods. The AD method includes the following steps (Hu et al., 2009). A discrete series, ${\it PTR}(i)$ with $i=1,\cdots,N$ is partitioned with overlapping windows of length $2n+1$ and, accordingly, each neighboring segment has $n+1$ overlapping points. An arbitrary polynomial ${\cal Y}$ is constructed in each window of length $2n+1$. In order to have the continuous trend function avoid a typical sharp jump in it, the following weighted function for the overlapping part of the $\nu$th segment is considered (Hu et al., 2009):

[TABLE]

where $j=1,2,\cdots,n+1$. The two free parameters, namely $n$ and the order of the fitting polynomial, should be determined properly (Hu et al., 2009).

The size of each segment was calculated by $2n+1=2\times{\rm int}\left[(N-1)/(w_{\rm adaptive}+1)\right]+1$. It turns out that by increasing the value of $w_{\rm adaptive}$ and the order of the fitting polynomial, the fluctuations disappear, and, consequently, the fluctuations are suppressed. For the nonoverlapping segments, the AD data are given by ${\it PTR}(i)-{\cal Y}_{\nu}(i)$, while for the overlap part it is ${\it PTR}(i)-{\cal Y}_{\nu}^{\rm overlap}(i)$. Since the GW, as the dominant part of the signal, is our desired part of the signal, we instead use $\widetilde{{\it PTR}}(i)={\cal Y}_{\nu}(i)$, while for the overlap part, we consider $\widetilde{{\it PTR}}(i)={\cal Y}_{\nu}^{\rm overlap}(i)$. Now $\widetilde{{\it PTR}}(i)$ is used for further analysis in MF-DFA or MF-DXA.

2.5. Trend and noise modeling

In real observational data to carry out parametric detection, reliable statistical models of the noise and signal should be well established. A proposal for noise modeling is based on the denoiseing procedure carried out by the SVD or AD algorithms. Previously, we were interested in removing the contribution of undesired noise modulated on real data. Now we concentrate on the $\widetilde{{\it PTR}}$ given by Eq. (27) in the context of SVD analysis as a model of trends and $PTR-\widetilde{{\it PTR}}$ for noise. Also, if we use the AD approach, the global variation part of the signal corresponds to both $\mathcal{Y}$ and $\mathcal{Y}^{\rm overlapp}$ (Eq. (2.4)). For the noise part, we should consider ${\it PTR}(i)-{\cal Y}_{\nu}(i)$, while for the overlap part, it is ${\it PTR}(i)-{\cal Y}_{\nu}^{\rm overlap}(i)$. Therefore, SVD or AD, as well as the internal part of the MF-DFA and MF-DXA algorithms, are able to give a robust model for trends and noise. Also, extracting intrinsic functions based on EMD can be a good proposal for this purpose (Huang et al., 1998).

2.6. Posterior Analysis

In this paper, we turn to Bayesian statistics (Colistete et al., 2004) to compute the reliable value of the generalized Hurst exponent (Eqs. (10) and (11)). Let $\{{\mathcal{D}}\}:\{{\mathcal{F}}_{q}(s)\}$ and $\{\Upsilon\}:\{h(q)\}$ represent the measurements and model parameters, respectively. The posterior function is defined by:

[TABLE]

where $\mathcal{L}$ is the likelihood and $\mathcal{P}(\Upsilon)$ is the prior probability function including all information concerning model parameters. Here we adopt the top-hat function for ${\mathcal{P}}(h(q))$ in the interval $h(q)\in[0,4]$. According to the central limit theorem, the functional form of likelihood becomes multivariate Gaussian, i.e. ${\mathcal{L}}({\mathcal{D}}|\Upsilon)\sim\exp(-\chi^{2}/2)$. The $\chi^{2}$ for determining the best-fit value for the scaling exponent coordinated by multifractal formalism reads as:

[TABLE]

where $\Delta\equiv[{\mathcal{F}}_{q}^{{\rm obs.}}-{\mathcal{F}}_{q}^{{\rm the.}}]$ and ${C}$ is the covariance matrix. The ${\mathcal{F}}_{q}^{{\rm obs.}}(s)$ and ${\mathcal{F}}_{q}^{{\rm the.}}(s;h(q))$ are fluctuation functions computed directly from the data and determined by Eqs. (10) or (11), respectively. In the case of the diagonal covariance matrix, the $\chi^{2}$ becomes:

[TABLE]

Here $\sigma_{{\rm obs.}}(s)=\left\langle\left[\delta\mathcal{F}_{q}^{{\rm obs.}}(s)\right]^{2}\right\rangle$, which is related to the diagonal elements of ${C}$ and can be computed using a standard statistical error propagator from primary uncertainties on $PTR$ datasets (Eq. (1) to Eqs. (8) and (9)). The $1\sigma$ error bar of $h(q)$ is determined by:

[TABLE]

Subsequently, we report the best value of the scaling exponent at a $1\sigma$ confidence interval as $h(q)_{-\sigma^{-}_{h(q)}}^{+\sigma^{+}_{h(q)}}$.

3. Data Description

In this section, we will describe theoretical models for GW signals. The observational datasets, synthetic series for pure timing residuals, and GWs, in order to examine the multiscaling behavior of $PTR$s as an indicator of GWs, will be described in this section.

3.1. Theoretical notions of the GWB on PTRs

The potential sources of GWs could be massive accelerated objects (Rajagopal, M., & Romani, 1995; Taylor & Gair, 2012), burst sources (Thorne & Braginskii, 1976; Damour & Vilenkin, 2001) or stochastic background sources (Maggiore, 2000; Hobbs, 2011; Damour & Vilenkin, 2005; Pshirkov & Tuntsov, 2010; Hobbs, 2011; Hobbs et al., 2009). Isotropic stochastic GWB produced by coalescing supermassive binary black hole mergers is the strongest potentially detectable signal of GWs (Hobbs et al., 2009). Therefore, we use the GWB model to produce synthetic data. The characteristic strain spectrum, $\mathcal{H}_{c}(f)$, for a stochastic GWB can be described by the power-law relation (Hobbs, 2011):

[TABLE]

where $f$ is the frequency of GWs, $f_{1yr}\equiv\frac{1}{1yr}$; $\mathcal{A}_{yr}$ is the dimensionless amplitude of the GWB; and $\zeta$ is a scaling exponent and for almost all expected GWs is $\zeta<0$. The corresponding $\zeta$ exponent takes the following values for different mechanisms: $\zeta=-\frac{2}{3}$ for coalescing black hole binaries, $\zeta=-1$ for cosmic strings, and $\zeta=-\frac{7}{6}$ for primordial GWs from the Big Bang (Hobbs, 2011). We should mention that the power-law relation obtained in Equation (33) is not unique and there is another framework represented by (Yardley, 2004; Sesana et al., 2008). The dimensionless amplitude of GWs has been predicted by most authors in the range of $\mathcal{A}_{yr}\in[10^{-15},10^{-14}]$; however, according to Refs. (Yardley, 2004; Sesana et al., 2008) the expected range of $\mathcal{A}_{yr}$ for a stochastic GWB is $\mathcal{A}_{yr}\in[10^{-16},3\times 10^{-15}]$.

3.2. Synthetic Data Sets for GWB

To simulate synthetic series, we use the TEMPO2 software package that carries out the fitting procedure of TOA (Hobbs et al., 2006). This package is used to simulate pure timing residuals (Hobbs et al., 2009). To simulate the GWB, the ”GWbkgrd” plug-in of TEMPO2 will be used (Taylor & Gair, 2013). In the absence of GW signal, we have pure pulsar timing residual represented by ${\it PTR}_{\rm pure}$, while signal induced by GWB is indicated by ${\it PTR}(t)$.

In order to test the effect of GWs on the $PTR$s, we simulate $100$ timing residuals with 1076 data points that are separated by 13 days with an rms of 100 ns. Then we add the effect of GWB on the simulated pure $PTR$ using different seeds for a given $\mathcal{A}_{yr}$. The chosen accuracy for simulation has been used in other work as a level at which a GWB might be detected (Jenet et al., 2005); however, it should be noted that only two of the PPTA pulsars (J0437-4715 and J1909-3744) have rms noise of this order (Table 1).

The GWB introduces two terms for each polarization, one set of which is referred to as the Earth terms. These Earth terms are correlated. However, the other set, referred to as the pulsar terms, has equal amplitude but a long and unknown time delay, so these terms are effectively uncorrelated noise with the same red spectrum as the Earth terms. Our simulations include both the Earth and the pulsar terms. We simulate 20 pure PTRs for pulsars separated in the sky according to the ephemeris of 20 MSPs observed in the PPTA project (Table 1). An isotropic GWB induces a particular spatial cross-correlation in $PTR$s leading to a quadrupolar signature (Hellings and Downs curve) (Hellings & Downs, 1983; Jenet et al., 2005). Subsequently, to examine the GWB, we will examine the cross-correlation property of the simulated data.

The upper panel of Fig. 1 indicates a typical pure timing residual simulated by TEMPO2 with zero mean uncorrelated series. We also depict the superposition of pure timing residuals with the GW model introduced in (Hobbs et al., 2009), in the middle panel of Fig. 1.

3.3. Observed Data

We use the timing residual data of 20 MSPs observed by the PPTA project at three bandwidths, namely $10$, $20$, and $50$cm, by implementing the Parkes 64 m radio telescope (PTA) (Manchester et al., 2013). The PTA telescope is located in Australia at an altitude of -33*∘* and can observe all of the inner Galaxy. Due to the higher stability of the short-period MSPs, the observed pulsars have short periods and are selected from bright ones. Also, these MSPs have narrow pulse widths in order to reduce uncertainties in the corresponding TOA. Finally, isolated wide-binary MSPs have been selected to avoid the effects of the companion star.

The PTR series for these MSPs as observed datasets are publicly available222 https://datanet.csiro.au/dap/. We have used the TEMPO2 software to extract post-fitted PTRs from timing model data presented by Manchester et al. (2013).

The spectralModel333http://www.atnf.csiro.au/research/pulsar/tempo2 plug-in is utilized for temporal smoothing and making an equally spaced grid of observed data (Coles et al., 2011). Then, we applied our analysis on post-fitted data.

The names of 20 MSPs with the corresponding rms and total time span are reported in Table 1. It is worth noting that several phenomena, such as atmospheric delays, vacuum retardation due to observatory motion, Einstein delay, and Shapiro delay, can affect the TOA (Edwards et al., 2006) and they should be dismissed to have a post-fitted timing residual that is called ${\it PTR}$. The lower panel of Fig. 1 illustrates a typical post-fit pulsar timing residual of PSR J0437-4715 observed by the PPTA project (Manchester et al., 2013). The fitting procedure has been done with the TEMPO2 software.

4. Multifractal analysis of synthetic PTR series

In this section we will evaluate the multifractal nature of synthetic datasets. The capability of our analysis as a detector of gravitational waves and a pipeline for determining the type of GWB will be explained in this section.

4.1. Multifractal nature of synthetic Data

At first, we examine the multifractal nature of synthetic ${\it PTR}_{\rm pure}$ and its superposition with simulated GWB. Since, in simulation, our data are regular, therefore we apply common assessment algorithms. Fig. 2 illustrates the fluctuation functions versus scale computed by DMA for $PTR_{\rm pure}$. The results derived by the DFA method are in agreement with the DMA algorithm. The average value of the Hurst exponent for all simulated pure $PTR$s is $\langle H\rangle=0.51\pm 0.02$ at a $1\sigma$ level of confidence, confirming that $PTR_{\rm pure}$ is an uncorrelated data set (Hobbs, 2011). Now we superimpose the synthetic $PTR_{\rm pure}(t)$ with simulated GWB with a given set of free parameters.

We apply DFA and DMA on simulated $PTR(t)$ for various GWB amplitudes. Fig. 3 illustrates $\mathcal{F}_{2}(s)$ as a function of $s$ for the simulated series. These results confirm that there is at least one crossover in fluctuation function versus $s$. We should eliminate the crossover in fluctuation function to determine the generalized Hurst exponent. To this end, we apply either SVD or AD to the datasets, and the clean series are used for further analysis by either the DFA or DMA methods. For SVD, we consider $p=1$ and $d=40$; therefore, the three largest eigenvalues are set to zero, and the new eigenvalues matrix ($\mathbf{\tilde{S}}$), filtered matrix ($\mathbf{\tilde{\Gamma}}$), and cleaned data ($\widetilde{{\it PTR}}$) are constructed.

Fig. 3 indicates $\mathcal{F}_{2}(s)$ computed by the DFA and DMA algorithms after applying either the SVD or AD method. Gu & Zhou (2010) demonstrated that DMA with $\theta=0$ (backward) has the best performance; therefore, we use the backward DMA method throughout this paper. We deduce that applying an SVD preprocess can efficiently remove the crossover, and we are able to assign a scaling exponent for fluctuation function versus $s$. The situation for AD preprocessing is somehow different, but it is consistent with the SVD results. The generalized Hurst exponent and $\xi$ versus $q$ for three types of $PTR$s superimposed by different values of GWB amplitudes are depicted in Fig. 4. The upper panels of Fig. 4 illustrate the $h(q)$ and $\xi(q)$ for synthetic $PTR$s affected by GWB with different amplitudes with the same $\zeta$. As we expect, the value of $h(q=2)$ that is related to $\zeta$ for all samples is almost same. The lower panel shows $h(q)$ and $\xi(q)$ for simulated $PTR$s with different $\zeta$.

4.2. Irregular MF-DXA of simulated $PTR$s

The quadrupolar signature of the spatial cross-correlation function of $PTR$s is considered as a particular measure for detecting the imprint of the GWB (Taylor et al., 2017). Previously, the Hellings and Downs curve has been examined for detection of stochastic GWB (Hellings & Downs, 1983; Jenet et al., 2005; Taylor et al., 2017).

Implementation of Irregular MF-DXA on $PTR$s provides a reliable cross-correlation exponent and coefficient in the presence of nuisance trends and noises. Irregular MF-DXA is indeed a crucial part of our pipeline for searching the significance of GWB. Here, due to the regularity of the simulated data, we consider the usual MF-DXA. To show the validity of this idea, we simulate 20 pure $PTR$s for pulsars separated in the sky according to the ephemeris of 20 MSPs observed in the PPTA project given in Table 1. Then, we add the effect of GWB to each pure $PTR$.

In Fig. 5, we show $\bar{\sigma}_{\times}(\Theta)$ for simulated $PTR$s. Here we have simulated 50 realizations for 20 pulsars. The points plotted in Fig. 5 are the average of these 50 realizations. As indicated in this figure, when synthetic $PTR$s are affected by GWB with $\mathcal{A}_{yr}=50\times 10^{-17}$ and $\zeta=-2/3$, we can recognize a quadrupolar feature in $\bar{\sigma}_{\times}(\Theta)$ which is a benchmark for existing GWB. This behavior is similar to the Hellings and Downs curve indicated by the dashed line in Fig. 5. One of the advantages of this new measure is that, when undesired parts exist in the observed series, we are able to infer the contribution of the GWB signal robustly. Eqs. (17) and (18) also confirm that $\bar{\sigma}_{\times}$ is almost insensitive to the value of $\mathcal{A}_{yr}$. To make a more conservative pipeline for assessing the GWB signal, it is necessary to compute the cross-correlation coefficient, $\bar{\sigma}_{\times}(\Theta)$, in addition to the usual spatial cross-correlation function known as the Hellings and Downs curve. After obtaining the feature, we carry out the rest part of the MF-DMA analysis to determine the type and amplitude of GWB signal.

4.3. Strategies for Searching GWs

According to the results presented in the previous sections, the randomness of pure $PTR$s exhibits that deviations from uncorrelated behavior can be considered as additional features presented in the recorded data. Unfortunately, the observed PTRs may include intrinsic fractal noise, interstellar plasma, uncertainties in the Earth’s motion, master clocks, and receiver signals. It has been demonstrated that the noise from some of these sources is wavelength dependent and has spatial correlation, either monopole or dipole in nature. Subsequently, relying on multifractal analysis modified by preprocessing algorithms such as the AD or SVD methods of individual PTRs probably gives rise spurious results in the framework of GW searching. To get rid of the effect of undesired components, we rely on the quadrupole structure of the GWB and carry out the irregular MF-DXA approach.

Therefore, we begin with Irregular-MF-DXA on all available $PTR$s distributed over all directions and then compute $\bar{\sigma}_{\times}(\Theta)$ as a function of separation angle, $\Theta$. The existence of a feature similar to Fig. 5 in observed $PTR$s would imply detection of a GWB. Note that Fig. 5 is the average of 50 realizations. One observation with these parameters would have error bars almost 7 times larger, so the GWB would be detected but the significance would be much less. Thereafter, we will turn to the multifractal behavior of the $PTR$ series to determine the type and amplitude of the GWB. In order to determine the type of stochastic GWB with a strain spectrum modeled by Eq. (33), after preprocessing to remove noise and foreground, we apply multifractal methods to compute a reliable Hurst exponent. This exponent is related to the power-spectrum exponent. Finally, the best-fit value of $\zeta$ and its associated error bar are determined (Hobbs et al., 2009). However, there are many complications in the real data sets, making the inference procedure less straightforward to assess GWs. We therefore introduce four criteria as follows:

$I)$ According to Eqs. (11) and (12), the intercept of fluctuation function for PTRs contains the intensity of superimposed GWs. Therefore, after recognizing a quadrupolar signature in analyzing pairs of PTRs, the following quantity is able to indicate the intensity of GWB: $\Delta h_{1}(\mathcal{A}_{yr},\zeta)\equiv\sum_{q=q_{min}}^{q=q_{max}}|\mathcal{G}_{h(q)}(\mathcal{A}_{yr},\zeta)-\mathcal{G}_{h(q)}(\mathcal{A}_{yr}=0)|$. In practice, we find a robust mathematical relation between $\Delta h_{1}(\mathcal{A}_{yr},\zeta)$ and $\mathcal{A}_{yr}$ for any given $\zeta$ (or, equivalently, $H$) and rms of white noise, as follows. We do many simulations for a given value of $\zeta$ with different $\mathcal{A}_{yr}$ values. Then, we apply either SVD or AD to make clean data. The clean data are used for further analysis. According to our simulation for $\zeta=-2/3$ and rms=100 ns, the mathematical relation between $\mathcal{A}_{yr}$ and $\Delta h_{1}$ in the range of $\mathcal{A}_{yr}\in[10^{-17},10^{-15}]$ reads as:

[TABLE]

where $a=(-1.15\pm 0.40)\times 10^{12}$, $b=(2.84\pm 0.54)\times 10^{7}$ and $c=-74.45\pm 16.88$. This fitting function is not unique, and here we select one with a high goodness of fit before going further. Also, for any other rms dictated by experiment, the above analysis should be repeated again to find the corresponding fitting function.

$II)$ For pure PTRs, we found that the Hurst exponent is almost $0.5$, while there will be deviations in the generalized Hurst exponent for $PTR$ signals affected by GWs (Eq. (33)) for a given amplitude $\mathcal{A}_{yr}$, and $\zeta$. Therefore, another powerful measure to quantify the intensity of the GWB would be $\Delta h_{2}(\mathcal{A}_{yr},\zeta)\equiv\sum_{q=q_{min}}^{q_{max}}|h(q;\mathcal{A}_{yr},\zeta)-h_{\rm shuf}(q;\mathcal{A}_{yr},\zeta)|$. Where $h_{\rm shuf}(q;\mathcal{A}_{yr},\zeta)$ is for completely randomized $PTR$ and ”shuf” refers to shuffled. In practice, we find a robust mathematical relation between $\Delta h_{2}(\mathcal{A}_{yr},\zeta)$ and $\mathcal{A}_{yr}$ for any given $\zeta$ (or, equivalently, $H$) and rms of white noise. The corresponding shuffled series are produced using original series. Now by calculating the generalized Hurst exponent for original and shuffled data, one can compute $\Delta h_{2}$. We find that the following function is a good fit to our simulations for $\mathcal{A}_{yr}$ in the range of $\mathcal{A}_{yr}\in[10^{-17},10^{-15}]$ versus $\Delta h_{2}$ for $\zeta=-2/3$ and rms=100 ns:

[TABLE]

where $a=0.19\pm 0.06$, $b=-1.57\pm 0.92$, and $c=7.40\pm 3.30$. This fitting function is not unique, and here we select a high goodness of fit. Before going further, it is worth noting that the whitened noise generation is serious in many simulations. An optimal algorithm to evaluate noise quality in many simulations, especially in data generation by the TEMPO2 software, can be carried out by the shuffling procedure explained here. Subsequently, our proposal in this regard can be straightforwardly implemented as a new plug-in.

$III)$ Since GWs may induce non-Gaussianity in $PTR$, it is interesting to take into account $\Delta h_{3}(\mathcal{A}_{yr},\zeta)\equiv\sum_{q=q_{min}}^{q_{max}}|h(q;\mathcal{A}_{yr},\zeta)-h_{\rm sur}(q,\mathcal{A}_{yr},\zeta)|$. In the mentioned criterion, $h_{\rm sur}(q;\mathcal{A}_{yr},\zeta)$ is the generalized Hurst exponents computed for Gaussian datasets with the same correlation function as the original series. Here ”$\rm sur$” represents surrogated data or phase-randomized surrogated series, including the multiplication of Fourier-transform data by a random phase with a uniform distribution function (Prichard & Theiler, 1994). We simulated the $PTR$ accompanying the GWB with different amplitudes, and the following fitting function is determined for $\mathcal{A}_{yr}$ in the same range as above versus $\Delta h_{3}$ for $\zeta=-2/3$ and rms=100 ns:

[TABLE]

where $a=68.03\pm 11.73$ and $b=-321.50\pm 65.10$.

$IV)$ The width of the singularity spectrum, which quantifies the nature of multifractality, is another benchmark for determining the amplitude of GWs superimposed on the $PTR$s. This measure is defined by $\Delta h_{4}(\mathcal{A}_{yr},\zeta)\equiv|\Delta\alpha(\mathcal{A}_{yr},\zeta)-\Delta\alpha(\mathcal{A}_{yr}=0)|$. According to our simulations, we find:

[TABLE]

for $\zeta=-2/3$ and rms=100 ns in the range of $\mathcal{A}_{yr}\in[10^{-17},10^{-15}]$. Here $a=106.30\pm 7.80$, $b=1.62\pm 0.42$, and $c=1.52\pm 9.74$.

Let us summarize our strategy based on the above criteria for searching GWs in observation. As explained in section 2, in the case of the proper value of signal-to-noise (S/N) for each observed PTR, we remove all known contributions from foreground contamination. Therefore, we make regular series according to methods explained in subsection 2.3. Now we are ready to apply either AD or SVD method to extract the dominant part of the signal (the trend part) from the noise. Then, we apply the MF-DXA method to compute $\bar{\sigma}_{\times}$, and we compute the spatial cross-correlation to identify the probable quadrupolar signature. In the case of finding the mentioned signature, we go through the detection of GWs. Otherwise, we can only carry out the upper-limit approach. We also apply irregular MF-DXA on the proper part of the series for all available pairs of observed $PTR$s to examine the temporal part of the cross-correlation function and deduce the temporal scaling exponent. In the case of the homogeneous and isotropic source of the GWB, $h_{\times}$ is independent from the angular separation of $PTR$s, while for anisotropic or different single sources of GWs, the scaling exponent of the temporal part of the cross-correlation gets various values for different pairs. Utilizing either irregular MF-DFA or irregular MF-DMA on cleaned data leads to computing $h(q)$. The best-fit value of $\zeta$ is then determined by using the power-spectrum exponent. Following the benchmarks, we compute $\Delta h_{1}$, $\Delta h_{2}$, $\Delta h_{3}$, and $\Delta h_{4}$ for the observed $PTR$s. The GWB amplitude can be conservatively read from the corresponding plots, as indicated in Fig. 6 or stated by Eqs (34), (35), (36) and (37). It is worth noting that the functional form of $\Delta h$ should be determined for each value of $\zeta$ and given rms of white noise associated with observed data. Finally, we are able to compute the upper limit on $\mathcal{A}_{yr}$ using posterior analysis (see section 5). Fig. 7 is a schematic representation of the pipeline.

Here we emphasize some important considerations for dealing with observed $PTR$s. First of all, we define a relative difference between the scaling exponent computed for the observed $PTR$s and that computed for the $PTR$s without GWB to reduce the contribution of noise and trends. Finally, in our approach, the level of noise is almost no longer serious when we focus on the scaling exponent.

5. Implementation of multifractal methods on observed PTR data

Here we use the MF-DFA and MF-DMA methods modified by either AD or SVD detrending procedures to examine the multifractal and complexity behavior of observed pulsar timing residuals.

5.1. Implementation on Observed Data

As discussed in subsection 2.3, observed $PTR$s datasets are in the form of irregularly sampled series, and here we use the spectralModel plug-in for the temporal smoothing algorithm to construct equidistant regular series for further analysis (Coles et al., 2011). The size of the current observed data is not large enough to use the irregular version of MF-DFA and MF-DMA introduced by Eqs. (23) and (24).

Fig. 8 illustrates the MF-DMA results for various observed PTRs. These results confirm that there is a crossover in fluctuation functions versus $s$, corresponding to $s_{\times}\sim 70$ days. For the scaling exponent for $s<s_{\times}$, we have $h(q=2)\in[1.03,1.82]$, demonstrating that datasets have a nonstationary nature, while for $s>s_{\times}$, we find $h(q=2)\in[0.07,1.55]$.

In order to get rid of these crossovers and have a scaling behavior in fluctuation functions, we apply either AD or SVD separately on modified observed datasets. Then, the cleaned data will be used as input for the MF-DFA and MF-DMA algorithms. Fig. 9 illustrates a typical observed $PTR$ (red line) and the trend (black line) determined by AD (upper panel) and SVD (lower panel). The corresponding residual between the observed data and trend is indicated in the bottom of this figure. Fig. 10 represents the fluctuation functions computed for a typical observed $PTR$ by DFA and DMA applied on cleaned data provided by AD and SVD separately. The slope of the fluctuation functions for $q=2$ in reliable scales is $h(q=2)\in[1.56,1.87]$, demonstrating that all underlying series are categorized in the nonstationary class. The corresponding Hurst exponent, $H=h(q=2)-1$, belongs to $H\in[0.56,0.87]$. The value of the Hurst exponents for all observed PTRs at the $68\%$ level of confidence is depicted in Fig. 11. This result confirms that the dominant part of observed PTRs belongs to the long-range correlated signal. The lower panel of Fig. 11 shows the $q$-dependency of the generalized Hurst exponent after applying SVD on observed data and determined by MF-DMA. The results for MF-DFA are consistent with those determined by MF-DMA. Since $h$ depends on $q$, we conclude that all observed PTRs are multifractal. Singularity spectra of some observed PTRs are plotted in the upper panel of Fig. 12. The strength of the multifractal nature of PTRs is determined by the width of the singularity spectrum, $\Delta\alpha=\alpha_{\rm max}-\alpha_{\rm min}$. This value for observed data is reported in Table 1 and is also shown in the lower panel of Fig. 12. The range of the mentioned singularity spectra is $\Delta\alpha\in[0.89,1.79]$. Other relevant exponents are reported in Table 1.

An interesting question is, what are the sources of multifractality of observed PTRs? As explained in more detail by Kantelhardt et al. (2002), in principle, different correlation functions at small and large fluctuations can be considered as a source of multifractality. In addition, heavy-tailed probability distribution contributes to the multifractal behavior. In order to distinguish the two mentioned types of multifractality, we follow the method introduced in (Kantelhardt et al., 2002). By shuffling the series, the scaling behavior of the ratio of fluctuation functions, $\mathcal{F}_{q}(s)/\mathcal{F}_{q}^{\rm shuf}(s)$, is represented as:

[TABLE]

where $h_{\rm{shuf}}(q)$ is the generalized Hurst exponent for shuffled data. The case of $h_{\rm{cor}}(q)\equiv h(q)-h_{\rm{shuf}}(q)=0$ refers to multifractality sourced by the distribution function. In this case, we can compute $h_{\rm PDF}(q)\equiv h(q)-h_{\rm sur}(q)$. If both $h_{\rm{cor}}(q)$ and $h_{\rm{PDF}}(q)$ depend on $q$, both sources are playing roles in the multifractality of the data. In our samples, all PTRs have $h_{\rm{shuf}}(q)=0.50$ at a $1\sigma$ confidence interval, confirming that the correlation in datasets is almost the main source of multifractality. This property is a universal feature of all observed PTRs investigated in this paper.

The multifractality responsible for observed PTRs can also be examined by our method. To this end, we have used different models for the noise component according to the SimRedNoise plug-in of TEMPO2 and applied the MF-DMA method on those series. The upper panel of Fig. 13 indicates that the width of the singularity spectrum computed by the MF-DMA method is almost independent of the amplitude of the red-noise model. The lower panel illustrates the dependency of $\Delta\alpha$ on the exponent of the red-noise power spectrum considered as the $P_{red}(f)=A_{red}(1+f^{2}/f_{c}^{2})^{-Q/2}$ model, where $A_{red}$, $f_{c}$, and $Q$ are the amplitude of the power spectrum, corner frequency, and power-law index, respectively (Hobbs et al., 2006; Archibald et al., 2016). In this equation, $Q=0$ corresponds to white noise, and $Q=2,4$, and $6$ are related to phase noise, frequency noise and spin-down noise. Subsequently, we can deduce that the red noise can be responsible for multifractality of observed PTRs as well as GWs. Therefore, a part of our reported multifractality is related to red noise.

In Fig. 14, we indicate $\bar{\sigma}_{\times}$ as function of $\Theta$ for 20 MSPs observed in the PPTA project (listed in Table 1). We have not obtained an obvious quadrupolar signature for the mentioned observed series due to the high value of rms, short length in the size of the data, unresolved foreground contamination, and systematic noise. In the next subsection, we will go through the finding upper limit for the amplitude of the probable GWB superimposed in observed $PTR$s.

5.2. Upper bound on GWB amplitude

Multifractal assessment of individual $PTR$s series is not adequate to make a decision on the significance of the stochastic GWB. Therefore, inspired by the unique signature of the GWB, i.e. the quadrupolar feature induced on the spatial correlation function of $PTR$s fluctuations, we apply multifractal cross-correlation analysis. This is a generalized function including spatial-temporal cross-correlation function and has some novelties compared to the standard spatial cross-correlation analysis. Our algorithm is a proper method for denoising and detrending.

The irregular MF-DXA applied to the observed irregular $PTR$s did not yield reliable results for detecting GWB due to the limited size and low S/N of the data. Different criteria introduced in this paper will enable us to detect the footprint of possible GWs with a future generation of surveys with high-S/N observations. Now we turn to assigning an upper bound on probable GWB amplitude.

Previous studies have mainly considered a model for the power spectrum of the PTR signal modulated by GWB, including the amplitude and scaling exponent of GWB. According to priors associated with the model parameters, the Bayesian method has been adopted (Lentati et al. (2013); Shannon et al. (2015) and references therein). In our approach, we proceed with our strategies for searching the GWB (subsection 4.3). The posterior probability function, ${\mathcal{P}}_{\diamond}(\mathcal{A}_{yr}|{\mathcal{D}})$, reads as:

[TABLE]

Here symbol ”$\diamond$” corresponds to one of four measures proposed for determining the amplitude of the stochastic GWB, and $\delta_{D}$ is the Dirac delta function. The $\Phi_{\mathcal{D}}(\Delta h_{\diamond})$ represents the functional form presented in Fig. 6. The integral form of Eq. (39) is given by:

[TABLE]

in which $|\mathcal{J}|$ is the Jacobian computed for $\Delta h^{\prime}_{\diamond}=\Phi_{\mathcal{D}}^{-1}(\mathcal{A}_{yr})$. Finally, the upper bound on $\mathcal{A}_{yr}^{up-\diamond}$ can be determined by:

[TABLE]

where $C.L.^{\diamond}$ and $\mathcal{A}_{yr}^{up-\diamond}$ are the confidence interval and upper limit associated with one of our strategies, respectively. According to the posterior function defined by Eq. (29), considering $\{\mathcal{D}\}=\{\Delta h_{\diamond}^{PTR}\}$ for a given observed pulsar called by $PTR$ and $\{\Upsilon\}=\mathcal{A}_{yr}$, we compute:

[TABLE]

where $\Delta_{PTR}\equiv\left[\Delta h^{PTR}-\langle\Delta h(\mathcal{A}_{yr})\rangle\right]$ and ${\mathcal{C}}_{\mathcal{A}_{yr}}$ is the $4\times 4$ covariance matrix of the four statistical features defined by $\Delta h_{1}$, $\Delta h_{2}$, $\Delta h_{3}$ and $\Delta h_{4}$ (see Eqs. (34), (35), (36) and (37)). The $\langle\Delta h({\mathcal{A}}_{yr})\rangle$ is the average of $\Delta h$ over 1000 synthetic datasets for a given $\mathcal{A}_{yr}$, where $\mathcal{A}_{yr}\in[10^{-16},10^{-14}]$ and with a step size of $5\times 10^{-16}$. According to the likelihood function, ${\mathcal{L}}({\Delta h^{PTR}}|\mathcal{A}_{yr})\sim\exp(-\chi^{2}(\mathcal{A}_{yr})/2)$, the $95\%$ upper bound on $\mathcal{A}_{yr}^{up}$ using the observed $PTR$s is defined by:

[TABLE]

We report the computed upper bound for some observed pulsar timing residuals at a $95\%$ confidence level in Table 1. One may note that the upper bound on $\mathcal{A}_{yr}$ has not been reported for some observed $PTR$s. This is because, for such cases, the upper value is not in the range of $\mathcal{A}_{yr}\in[10^{-16},10^{-14}]$ considered in this research. Our results are consistent with other reports (Shannon et al., 2015).

6. Summary and conclusion

The $PTR$ is a good indicator to examine relevant physical phenomena from the interior of pulsars, as well as cosmological events. In spite of high stability in some types of pulsars, $PTR$s are classified as stochastic processes due to superimposed unknown trends and noises. The GWs produced by either primordial or late events affect the $PTR$s. Therefore, quantifying the fluctuations of $PTR$s can be a proper measure for GW detection.

In this paper, for the first time, we utilized a multifractal approach in order to examine the statistical properties of synthetic and observed $PTR$s affected by trends and noises. In the presence of trends and unknown noises, only robust methods are able to recover the correct multifractal nature of underlying series. In this research, we used MF-DFA, MF-DMA, and MF-DXA modified by the preprocessors, so-called AD or SVD algorithms. The pulsar timing observations are unevenly sampled datasets. To mitigate this property, we modified some internal parts of the multifractal analysis and proposed the irregular MF-DXA method and examined its accuracy. Our results demonstrated that computed scaling exponents for anticorrelated and long-range-correlated irregular signals are consistent with the expectations.

We used synthetic $PTR$s simulated by the TEMPO2 pulsar timing package. A template proposed by Hobbs et al. (2009) was used to take into account the contribution of GWs. We simulated 1000 synthetic $PTR$s, and the MF-DFA, MF-DMA, and MF-DXA methods were implemented on the simulated series. Our results demonstrated that the ensemble average of the Hurst exponent of the simulated data is $\langle H\rangle=0.51\pm 0.02$, confirming that the pure $PTR$s belong to monofractal uncorrelated stationary processes. There is no crossover in fluctuation functions versus scale determined by MF-DFA and MF-DMA (Fig. 2). Adding mock GWB signal on pure $PTR$s leads to crossovers in the log-log plot of $\mathcal{F}_{2}$ as a function of $s$ as indicated in Fig. 3. To examine the scaling behavior of $PTR$s induced by GWs, we carried out either the SVD or AD method on the data. We found that SVD can remove the crossover on fluctuation function for any $q$. The time scale for crossover depends on the intensity of the GW signal. In the presence of GWs, $PTR$s belong to a multifractal process due to the $q$-dependency of the generalized Hurst exponent, $h(q)$, (Fig. 4). Therefore, we were able to classify the mentioned data in the universal class of the multifractal process. The value of multifractality increased by increasing the intensity of GWs.

Various components of a recorded PTR may behave as a scaling fluctuation. Therefore, applying a multifractal algorithm on individual PTRs may give spurious results in exploring GWs. We relied on quadrupolar structure associated with the impact of GWB on the spatial cross-correlation of $PTR$s. We carried out cross-correlation analysis by the irregular MF-DXA introduced in this paper on all available $PTR$s distributed in all directions. To this end, we defined a new cross-correlation function (Eq. (18)) and accordingly, we computed the ensemble average of $\langle\sigma_{\times}(\Theta_{ab})\rangle_{\rm pair}$ for all synthetic $PTR$s as a function of separation angle, $\Theta$. We obtained an analogous behavior as a quadrupolar signature in $\bar{\sigma}_{\times}$. According to a model for GWB, obviously, the temporal part must be independent from the separation angle of the $PTR$ pairs affected by isotropic GWB, while the amplitude of cross-correlation defined by the DXA method illustrates the Hellings and Downs curve (Fig. 5) similar to the usual spacial crosscorrelation.

We proposed four criteria to quantify the footprint of GWs on pulsar timing residuals. Comparing the $y$-intercept of fluctuation functions with the one computed for pure $PTR$s is our first measure. The second measure is devoted to the generalized Hurst exponent with the one computed for pure $PTR$s. Comparison between $h(q)$ and the generalized Hurst exponent computed for the Gaussian signal is the third criterion. The fourth criterion corresponds to the width of the singularity spectrum.

The strategy for GWB detection in observations is as follows. After removing foreground and systematic noise by applying either SVD or AD on datasets, cleaned data that are associated with the dominant part of the signal (the trend part) will be used as input for irregular MF-DXA. Having observed relevant features for GWB on PTRs, irregular MF-DFA or irregular MF-DMA methods are applied exclusively. The type of superimposed GWs can be recognized by determining the Hurst exponent. Finally, the dimensionless amplitude of expected GWB ($\mathcal{A}_{yr}$) can be determined by inserting relevant quantities extracted by our four measures given by Eqs. (34), (35), (36) and (37) for a given $\zeta$ and rms of white noise determined in observations.

There is a crossover in the log-log plot of fluctuation function versus window length of observed $PTR$s. For $s<s_{\times}$ and $s>s_{\times}$, the exponents $h(q=2)$ are $h(2)\in[1.03,1.82]$ and $h(2)\in[0.07,1.55]$, respectively. After applying SVD, the corresponding Hurst exponent is $H\in[0.56,0.87]$.

The $q$-dependency of $h(q)$ confirmed that all observed MSPs behave as multifractal fields. The relevant exponents for observed MSPs have been reported in Table 1. The source of multifractality is mainly the correlation in small and large scales and is a universal property of all observed pulsars examined in this paper. The contribution of red-noise model indicated the extra multifractality on observed MSPs. Consequently, the degree of multifractality reported for PPTA data sets is the upper value, and a part of this value is associated with the noise model.

To infer the statistical significance of the GWB impact on the $PTR$s, we computed $\bar{\sigma}_{\times}(\Theta)$ for 20 MSPs observed in the PPTA project. Due to a high value of rms and a short length in the size of recorded data, we have not found a quadrupolar signature. Thereafter, we computed the upper bound for PSRs reported in Table 1.

Final remarks are as follows. The observed $PTR$s are affected by noises classifying into intrinsic and extrinsic categories Hobbs et al. (2006); Caballero et al. (2016). Reliable statistical models for noise and signal were introduced. The shuffling procedure and its evaluation by multifractal detrended analysis can also be implemented in TEMPO2 and other subroutines for simulation of $PTR$s. It could be interesting to simulate various kinds of GWs and to consider timing noise. Evaluation of different noise models and sensitivity to frequency is beyond the scope of this paper and will be considered elsewhere.

The authors thank M. Farhang for her useful comments on the manuscript. Also, the authors appreciate R. Monadi for his useful discussion. We also appreciate W. Coles for his comments on noise models in pulsar timing residual data sets. Thanks to the anonymous referee for the very extremely useful comments and for helping us to improve this paper extensively. SMSM is grateful to the HECAP section of ICTP, where some parts of this work have been finalized.