EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 707 (March 2026) A&A, 707 (2026) L8 Full HTML
Open Access
Issue
A&A
Volume 707, March 2026
Article Number L8
Number of page(s) 4
Section Letters to the Editor
DOI https://doi.org/10.1051/0004-6361/202658930
Published online 06 March 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1. Introduction

The data on the pulsar-timing monitoring obtained by the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) (Agazie et al. 2023a, 2025) and other collaborations opened a new important window to trace the variety of ongoing processes in the Universe. The released pulsar timing array (PTA) data are linked to the stochastic gravitational-wave background (GWB) in the nanohertz band produced by binary supermassive black holes (Sesana et al. 2008; Agarwal et al. 2026; Sesana & Figueroa 2025), and they are also used to constrain several models such as those pertaining to the nature of dark matter and the evolution of early cosmology (e.g., Tiruvaskar & Gordon 2026; Ben-Dayan et al. 2026 and references therein).

The key issue in the pulsar timing is the analysis of the stochastic properties of the time of arrival (TOA) data of residuals for the studied 68 pulsars, for example the differences between the observed data and modeled parameters (Lam et al. 2025). Therefore, the study of tiny properties of the stochastic background by different methods and hence, of the residuals’ association to white noise, can be crucial for both revealing the details from the pulsar astrometric parameters of pulsars up those of the supermassive black hole binaries (e.g. of their orbital eccentricities) or of the contribution of other effects (see Fiscella et al. 2025 and references therein). Moreover, white noise can contain contributions from instrumental noise, calibration of the radio telescopes, in addition to interference from terrestrial sources. Thus, while the comprehensive study of white noise is crucial for revealing the origin of the nanoherz gravitational wave background, the use of different descriptors can be of particular importance for the separation of contributions for the abovementioned effects in white noise among others.

For this study we applied the Kolmogorov stochasticity parameter (KSP) to analyze the properties of the white noise (cf. Barenboim et al. 2025) of TOA associated with both the TOA data and the pulsar timing models. The Kolmogorov stochasticity parameter is a sensitive quantitative descriptor of the degree of randomness of sequences (Kolmogorov 1933; Arnold 2008a). The KSP has been applied to the analysis of various signals (e.g., Arnold 2008a; Atto et al. 2013; Hoffman et al. 2023, including human genomic sequences (Gurzadyan et al. 2015), as well as those of an astrophysical origin (Gurzadyan & Kocharyan 2008; Gurzadyan et al. 2009, 2010; Gurzadyan & Durret 2011; Galikyan et al. 2025, see also Gurzadyan 1999). For example, the KSP enabled us to trace the properties of the cold spot as of a non-Gaussian region in the cosmic microwave background sky (Planck data) and to conclude on its void origin in the large-scale matter distribution (Gurzadyan et al. 2014); the void nature of the cold spot was soon after confirmed by an galaxy survey study (Szapudi et al. 2015).

Thus, we used the KSP for the following: (a) to quantify how independent the TOA residuals are; (b) to reveal how homogeneous the randomness is regarding individual pulsars, i.e., if there are pulsars attributed to low and high randomness signals; (c) to test if nonstationary effects can be traced in the white noise; and (d) to determine how closely the timing models for individual pulsars match their observations.

2. Kolmogorov stochasticity parameter

First, let us define the Kolmogorov stochasticity parameter as in (Kolmogorov 1933; Arnold 2008a,b,c, 2009a,b). Let X be a real-valued random variable. Its cumulative distribution function (CDF) is

F ( x ) ≡ Pr ( X ≤ x ) , Mathematical equation: $$ \begin{aligned} F(x) \equiv \Pr (X \le x), \end{aligned} $$(1)

where Pr(⋅) denotes the probability under the assumed probabilistic model. Throughout this section, F refers to the theoretical (model/null) CDF.

Given an observed sample x1, ..., xn (not necessarily distinct), let x(1) ≤ x(2) ≤ ⋅ ⋅ ⋅ ≤ x(n) denote the order statistics. The corresponding empirical CDF (ECDF) is defined by

F n ( x ) ≡ 1 n ∑ i = 1 n 1 { x i ≤ x } , Mathematical equation: $$ \begin{aligned} F_{n}(x) \equiv \frac{1}{n}\sum _{i=1}^{n}\mathbf 1 \{x_{i}\le x\}, \end{aligned} $$(2)

which is a right-continuous step function,

F n ( x ) = { 0 , x < x ( 1 ) , k / n , x ( k ) ≤ x < x ( k + 1 ) , k = 1 , ⋯ , n − 1 , 1 , x ( n ) ≤ x . Mathematical equation: $$ \begin{aligned} F_{n}(x) = {\left\{ \begin{array}{ll} 0,&x < x_{(1)}, \\ k/n,&x_{(k)}\le x < x_{(k+1)}, \quad k=1,\dots ,n-1, \\ 1,&x_{(n)}\le x. \end{array}\right.} \end{aligned} $$(3)

Our basic goal is to quantify the discrepancy between the ECDF of the data and a reference (theoretical) CDF F. In applications, F represents a null/model hypothesis (e.g., a specified distribution, or a distribution implied by a simulated noise model). When F is fully specified and continuous, the resulting statistic is distribution-free under the null.

When the (one-sample) Kolmogorov distance is defined,

D n ≡ sup x ∈ R | F n ( x ) − F ( x ) | , Mathematical equation: $$ \begin{aligned} D_{n} \equiv \sup _{x\in \mathbb{R} } |F_{n}(x) - F(x)|, \end{aligned} $$(4)

and for its scaled form, the Kolmogorov stochasticity parameter,

λ n ≡ n D n = n sup x ∈ R | F n ( x ) − F ( x ) | . Mathematical equation: $$ \begin{aligned} \lambda _{n} \equiv \sqrt{n}D_{n} = \sqrt{n} \sup _{x\in \mathbb{R} } |F_{n}(x) - F(x)|. \end{aligned} $$(5)

Because the sample { X i } i = 1 n Mathematical equation: $ \{X_{i}\}_{i=1}^{n} $ is random, the ECDF Fn and thus λn are random variables. Under the null hypothesis H0 : X1, ..., Xn ∼ F (with F continuous), λn has a distribution function

Φ n ( λ ) ≡ Pr ( λ n ≤ λ ) , Mathematical equation: $$ \begin{aligned} \Phi _{n}(\lambda ) \equiv \Pr (\lambda _{n}\le \lambda ), \end{aligned} $$(6)

which converges to a universal limiting distribution Φ(⋅) independent of F:

Φ ( λ ) = ∑ k = − ∞ + ∞ ( − 1 ) k e − 2 k 2 λ 2 , Φ ( 0 ) = 0 . Mathematical equation: $$ \begin{aligned} \Phi (\lambda ) = \sum _{k=-\infty }^{+\infty }{(-1)^k e^{-2k^2\lambda ^2}},\\ \Phi (0) = 0. \end{aligned} $$(7)

This is the statement of Kolmogorov’s theorem (Kolmogorov 1933).

Note that, according to this theorem and the property of the limiting function Φ, the probable values of λ are within 0.3 and 2.2. In other words, at its lower values, Φ(λ) tends rapidly to 0 and to 1 at higher values.

An important feature of the KSP criterion is its applicability to relatively small sequences (orbits of dynamical systems) as is explicitly shown in Arnold (2008a), which compared two sequences of 15 two-digit numbers and concluded that the stochasticity probability of the one sequence is approximately 300 times higher than for other one. Then, for random sequences, the stochasticity parameter λn has a distribution close to the function Φ(λ), while the distribution is different from that function for nonrandom sequences.

The numerical experiments performed in Gurzadyan et al. (2011), using a single scaling of the ratio of stochastic to regular components of the generated sequences including a uniform distribution, revealed the efficiency of the KSP as an indicator of the degree of randomness of composite signals.

3. Pulsar timing

Pulsar timing is based on measuring TOAs of radio pulses with high precision and comparing them to the predictions of a timing model. The timing model encodes all deterministic contributions to the TOAs–such as the pulsar’s spin evolution, astrometric parameters, binary motion (when present), and propagation effects. Any mismatch between the measured TOAs and the timing model – i.e., the timing residuals – contains information about additional physical or instrumental processes not captured by the deterministic model.

These stochastic contributions are commonly separated into low-frequency, temporally correlated “red noise” (e.g., intrinsic spin noise or a stochastic gravitational-wave background), which are the type of signals that researchers usually look for. It is important to note that not all pulsars exhibit a red-noise component, and the contribution to the red noise can come from various effects (Fiscella et al. 2025), which can introduce characteristic structures in timing residuals. High-frequency, uncorrelated “white noise” mainly arises from radiometer noise and measurement uncertainties, representing random fluctuations in the time series. The white-noise component can be calculated as WhiteNoise = TOA − TimingModel − RedNoise (ifpresent).

For our analysis, we used PINT (Luo et al. 2019, 2021; Susobhanan et al. 2024), a high-precision open-source package for pulsar timing data analysis written in Python. It is developed in collaboration with NANOGrav and is largely used for analysis (Agazie et al. 2023a,b; Smith et al. 2023), among other things. It allows one to read TOA data and load timing models with their parameters, for example, provided by NANOGrav (Agazie et al. 2025). Moreover, in PINT, one can simulate the white and red noise components of pulsar signals separately and estimate the white-noise contribution. The simulated white noise and the observed white noise should have the same statistical properties, i.e., the signal that we refer to as white noise in the observed data is indeed a white-noise signal. To test this, we used the KSP as a statistical characteristic of a signal and compared the observed signal with the simulated one as described in Sect. 4.1.

4. Analysis

Consider a time series of length N that represents the residual component of the TOAs. To obtain the KSP characteristic histogram for the series, we calculated λ inside a sliding window of length 10−1 ⋅ N with a sliding step of 10−3 ⋅ N. Inside the window, we randomly selected 80% of the points and calculated λ using the generalized normal distribution as the theoretical distribution. Its probability density is given via Eq. (8). The final λ for that window is the mean of the 80%s which was repeated five times. Namely, the algorithm estimated the distribution of λ, and its running allowed the λ distribution’s stability to be estimated four more times, as 80% of the subsamplings each yielded different distributions. This enabled us to ensure the stability of the procedure, to define the errors bars, and to decrease the influence of outliers, if any. The specific choice of parameters – the 80% subsampling and the 5 × 5 repetitions – were made in order to adequately represent the distributions for comparison, while avoiding both over-averaging and undersampling:

f β ( x ) = 1 2 α Γ ( 1 / β ) e ( | x − μ | α ) β α = Γ ( 1 / β ) Γ ( 3 / β ) ; μ = 0 . Mathematical equation: $$ \begin{aligned} \begin{aligned}&f_\beta (x) = \frac{1}{2\alpha \Gamma (1/\beta )}e^{\left(\frac{|x-\mu |}{\alpha }\right)^\beta } \\&\alpha = \sqrt{\frac{\Gamma (1/\beta )}{\Gamma (3/\beta )}}; \quad \mu =0 \end{aligned} .\end{aligned} $$(8)

In Eq. (8) we have only a single free parameter β that describes the sharpness of the distribution. The scale parameter α was chosen such that the theoretical distribution had a unit variance. We considered two scenarios of data normalization and a choice of β. First, when β was a global parameter of a time series as described in Sect. 4.1, and, second, where β was local and had a time dependence (Sect. 4.2).

4.1. β = const

Consider a time series of residuals of a pulsar, { R i } i = 1 N Mathematical equation: $ \{R_i\}_{i=1}^{N} $, where N is the length of the series. We assumed that they are identically distributed. In this scenario, we first normalized our series by subtracting the mean of the series R ¯ Mathematical equation: $ \overline{R} $ and dividing by its standard deviation ΔR. Hence, we obtain normalized series r i = R i − R j ¯ Δ R Mathematical equation: $ r_i = \frac{R_i-\overline{R^j}}{\Delta R} $.

We performed this procedure on observed and simulated data using PINT (Luo et al. 2019, 2021; Susobhanan et al. 2024). Simulations used the same models, and the parameters within them were used to extract residuals from the TOA. We observed that some of the pulsars’ KSP distributions are not similar, as shown in Fig. 1, while for others they are, as shown in Fig. 2.

Thumbnail: Fig. 1. Refer to the following caption and surrounding text. Fig. 1.

KSP distributions for the simulated/observed data of pulsars’ (J1022+1001, J1125+7819, J1903+0327) white noise, when the histograms do differ. Orange and blue histograms correspond to simulations and observations, respectively.

Thumbnail: Fig. 2. Refer to the following caption and surrounding text. Fig. 2.

Same as in Figure 1, but for sample pulsars (J1745+1017, J1910+1256) when the KSP distributions for the simulated/observed data are similar.

For a quantitative picture, we created five epochs of simulations for each pulsar presented, and calculated the reduced χ2 of the KSP histograms between observed and simulated data. To calculate χ2, we linearly interpolated histogram values and their corresponding errors into the same bins. Table 1 shows the minimum, maximum, mean, median, and standard deviation of χ2 for each pulsar. Thus, the key point of our analysis is that the same numerical scheme was performed both for the data of observations and simulations.

Table 1.

Reduced-χ2 data for each pulsar for β = const.

4.2. β = β(t)

Now we assume that the properties of the distribution may change over time, hence β has a time dependence. We calculated β around each residual using maximum likelihood estimation. The algorithm for finding βi for a residual time series Ri is the following:

  1. We weighed each residual value by weight,

    w ik ∝ exp [ − 1 2 ( t k − t i h ) 2 ] , k = 1 ; N , ¯ Mathematical equation: $$ \begin{aligned} w_{ik} \propto \exp \left[-\frac{1}{2} \left(\frac{t_k - t_i}{h}\right)^2\right], \qquad k= \overline{1;N,} \end{aligned} $$

    where ti, tk are the time moments for Ri and Rk, respectively, and h is a hyperparameter that was chosen via cross-validation.

  2. We calculated the weighted mean and standard deviation and obtained normalized time series for each time moment ti: r ik = R k − R i ¯ Δ R i Mathematical equation: $ r_{ik} = \frac{R_k - \overline{R_i}}{\Delta R_i} $.

  3. We estimated βi as

    β i = β [ − ∑ k w ik log f β ( r ik ) ] . Mathematical equation: $$ \begin{aligned} \beta _i = _\beta \left[ -\sum _k w_{ik}\log f_\beta (r_{ik})\right]. \end{aligned} $$

The following procedure of calculating the KSP histogram is the same, with the exception that we normalized the distribution inside the sliding window similarly to Section 4.1, i.e., without the weights. Table 2 shows the χ2 values for this scenario.

Table 2.

Same as in Table 1, but for β = β(t).

Figures 3, 4, 5 illustrate examples of three differentiable cases: non-similar distributions became similar; and similar and non-similar KSP distributions remained as they were.

Thumbnail: Fig. 3. Refer to the following caption and surrounding text. Fig. 3.

Case of pulsar J1125+7819, when the KSP distributions became similar after introducing β(t).

Thumbnail: Fig. 4. Refer to the following caption and surrounding text. Fig. 4.

KSP distributions for pulsar J1910+1256. Introducing the time dependence in β does not significantly affect χ2.

Thumbnail: Fig. 5. Refer to the following caption and surrounding text. Fig. 5.

KSP distributions for the pulsar J1022+1001 after introducing the time dependence in β. Although ⟨χ2⟩ is smaller in the β(t) approach (Table 2) than in the β = const case (Table 1), they are of the same order.

5. Conclusions

We analyzed NANOGrav pulsar timing data using Kolmogorov’s method. Namely, we concentrated on the white-noise component of the observations, as it is the remainder after subtracting the known physical contributions from the TOAs. Comparing the observed white noise with simulated white noise and using the Kolmogorov stochasticity parameter of the signal as a metric enabled us to probe the possible physical components in the TOAs. These physical components may be deterministic, and thus they might be included in the timing models, or stochastic, as part of the red noise (e.g., the GWB). The analysis was performed blindly with the same algorithm applied to all pulsars considered.

Assuming a generalized normal distribution – with a single free β parameter describing the sharpness of the distribution – as the underlying distribution for the white noise, we find both cases, when the white noise stochasticity agrees with the observations and simulations and cases where it does not. The introduction of a time-dependent β(t) improved the similarity; however, instances remain in which the observational and simulated stochasticity still disagree. This can indicate the existence of nonstationary physical processes influencing the pulsar timing.

The KSP method can be used to test different models affecting pulsar timing, including for example modified theories of gravity or models for the dark sector of the universe (e.g., Gurzadyan et al. 2023, 2025, 2026). Namely, then the task will be to reproduce TOAs with the observed KSP properties, especially depending upon the availability of higher precision and more observational data.

Acknowledgments

We are thankful to the referee for many valuable comments. We acknowledge the use of The NANOGrav 15-Year Data Set; https://zenodo.org/records/16051178 and ASCL Code Record of Astrophysics Source Code Library ASCL.net; ascl:1902.007.

References

  1. Agarwal, N., Agazie, G., Anumarlapudi, A., et al. 2026, ApJ, 988, L11 [Google Scholar]
  2. Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023a, ApJ, 951, L8 [NASA ADS] [CrossRef] [Google Scholar]
  3. Agazie, G., Alam, M. F., Anumarlapudi, A., et al. 2023b, ApJ, 951, L9 [CrossRef] [Google Scholar]
  4. Agazie, G., Kaplan, D. L., Susobhanan, A., et al. 2025, arXiv e-prints [arXiv:2510.16668] [Google Scholar]
  5. Arnold, V. I. 2008a, Nonlinearity, 21, T109 [Google Scholar]
  6. Arnold, V. I. 2008b, Uspekhi Mat. Nauk, 63, 5 [CrossRef] [Google Scholar]
  7. Arnold, V. I. 2008c, IC/2008/001, ICTP, Trieste [Google Scholar]
  8. Arnold, V. I. 2009a, Trans. Moscow Math. Soc., 70, 31 [CrossRef] [Google Scholar]
  9. Arnold, V. I. 2009b, Funct. An. Other Math., 2, 139 [Google Scholar]
  10. Atto, A. M., Berthoumieu, Y., & Megret, R. 2013, Entropy, 15, 4782 [Google Scholar]
  11. Barenboim, G., Ireland, A., & Stebbins, A. 2025, arXiv e-prints [arXiv:2511.13866] [Google Scholar]
  12. Ben-Dayan, I., Kumar, U., & Verma, A. 2026, J. High Energy Astrophys., 50, 100510 [Google Scholar]
  13. Fiscella, S. V. S., Lam, M. T., Agazie, G., et al. 2025, arXiv e-prints [arXiv:2509.21203] [Google Scholar]
  14. Galikyan, N., Kocharyan, A. A., & Gurzadyan, V. G. 2025, A&A, 696, L21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Gurzadyan, V. G. 1999, EPL, 46, 114 [Google Scholar]
  16. Gurzadyan, V. G., & Durret, F. 2011, EPL, 95, 69001 [Google Scholar]
  17. Gurzadyan, V. G., & Kocharyan, A. A. 2008, A&A, 492, L33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Gurzadyan, V. G., Allahverdyan, A. E., Ghahramanyan, T., et al. 2009, A&A, 497, 343 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Gurzadyan, V. G., Kashin, A. L., Khachatryan, H. G., et al. 2010, EPL, 91, 19001 [Google Scholar]
  20. Gurzadyan, V. G., Ghahramanyan, T., & Sargsyan, S. 2011, EPL, 95, 19001 [Google Scholar]
  21. Gurzadyan, V. G., Kashin, A. L., Khachatryan, H., et al. 2014, A&A, 566, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Gurzadyan, V. G., Yan, H., Vlahovic, G., et al. 2015, Roy. Soc. Open Sci., 2, 150143 [Google Scholar]
  23. Gurzadyan, V. G., Fimin, N. N., & Chechetkin, V. M. 2023, A&A, 677, A161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Gurzadyan, V. G., Fimin, N. N., & Chechetkin, V. M. 2025, A&A, 694, A252 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Gurzadyan, V. G., Fimin, N. N., & Chechetkin, V. M. 2026, A&A, 705, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Hoffman, C., Cheng, J., Ji, D., & Dabaghian, Y. 2023, PNAS, 120, e2218245120 [Google Scholar]
  27. Kolmogorov, A. N. 1933, G. Ist. Ital. Attuari, 4, 83 [Google Scholar]
  28. Lam, M. T., Kaplan, D. L., Agazie, G., et al. 2025, arXiv e-prints [arXiv:2506.03597] [Google Scholar]
  29. Luo, J., Ransom, S., Demorest, P., et al. 2019, Astrophysics Source Code Library [record ascl:1902.007] [Google Scholar]
  30. Luo, J., Ransom, S., Demorest, P., et al. 2021, ApJ, 911, 45 [NASA ADS] [CrossRef] [Google Scholar]
  31. Sesana, A., & Figueroa, D. G. 2025, arXiv e-prints [arXiv:2512.18822] [Google Scholar]
  32. Sesana, A., Vecchio, A., & Colacino, C. N. 2008, MNRAS, 390, 192 [NASA ADS] [CrossRef] [Google Scholar]
  33. Smith, D. A., Abdollahi, S., Ajello, M., et al. 2023, ApJ, 958, 191 [NASA ADS] [CrossRef] [Google Scholar]
  34. Susobhanan, A., Kaplan, D. L., Archibald, A. M., et al. 2024, ApJ, 971, 150 [Google Scholar]
  35. Szapudi, I., Kovacs, A., Granet, B. R., et al. 2015, MNRAS, 450, 11 [Google Scholar]
  36. Tiruvaskar, S., & Gordon, C. 2026, Phys. Rev., 113, 043501 [Google Scholar]

All Tables

Table 1.

Reduced-χ2 data for each pulsar for β = const.

In the text
Table 2.

Same as in Table 1, but for β = β(t).

In the text

All Figures

Thumbnail: Fig. 1. Refer to the following caption and surrounding text. Fig. 1.

KSP distributions for the simulated/observed data of pulsars’ (J1022+1001, J1125+7819, J1903+0327) white noise, when the histograms do differ. Orange and blue histograms correspond to simulations and observations, respectively.
In the text
Thumbnail: Fig. 2. Refer to the following caption and surrounding text. Fig. 2.

Same as in Figure 1, but for sample pulsars (J1745+1017, J1910+1256) when the KSP distributions for the simulated/observed data are similar.
In the text
Thumbnail: Fig. 3. Refer to the following caption and surrounding text. Fig. 3.

Case of pulsar J1125+7819, when the KSP distributions became similar after introducing β(t).
In the text
Thumbnail: Fig. 4. Refer to the following caption and surrounding text. Fig. 4.

KSP distributions for pulsar J1910+1256. Introducing the time dependence in β does not significantly affect χ2.
In the text
Thumbnail: Fig. 5. Refer to the following caption and surrounding text. Fig. 5.

KSP distributions for the pulsar J1022+1001 after introducing the time dependence in β. Although ⟨χ2⟩ is smaller in the β(t) approach (Table 2) than in the β = const case (Table 1), they are of the same order.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.