Robust 1-norm Periodograms for Analysis of Noisy Non-Gaussian Time Series with Irregular Cadences: Application to VLBI Astrometry of Quasars (2024)

The Lomb–Scargle periodogram calculation is a powerful technique designed to reveal and characterize the periodic components in observational data sequences, which finds a wide scope of applications. For a review of its properties and underlying assumptions from the user's perspective, see VanderPlas (2018). The need for this technique arises from the character of astronomical data (observational measurements), which are practically never evenly sampled in time. This makes the standard Fourier power spectrum analysis inapplicable for astronomical time series. Detection of orbiting exoplanets from precision radial velocities of host stars is one of the well-known use cases for the Least-Squares (LS) periodogram method (Hara & Ford 2023). The periodic component of the measured radial velocity sequence is caused by the reflex orbital motion of the host star orbiting the system's barycentre. The period of the main sinusoidal mode in the computed periodogram in this case estimates the orbital period of the planet, which often cannot be directly observed.

The periodogram method finds a somewhat less known application in precision astrometry of celestial bodies' positions. Binary stars with unresolved or dim companions have periodic signals, which are the harmonics of the orbital frequency, in either of the sky coordinates referenced to a fixed celestial frame. Given a significantly long and precise cadence of position measurements covering at least one orbital period, the more general approach is to directly fit a set of Kepler elements of the emerging explicitly nonlinear 2D model, which proves a daunting and ambiguous task in the presence of even a small admixture of statistical outliers (Goldin & Makarov 2006). A robust and reliable periodogram decomposition is a welcome alternative when a large amount of observational data has to be processed with a low output of true positives. The need for a resilient periodogram algorithm, which can produce meaningful results outside of the normal distribution of data points, also emerges in the interpretation of high quality photometric time series. Magnetically active stars, for example, often manifest complex structures of signals in their light curves with periodic modulation mixed with stochastic, unpredictable bursts of radiation (Makarov & Goldin 2017).

Our main goal for this study is to develop and test a modification of the classical 2-norm periodogram algorithm (also known as the Lomb–Scargle periodogram) based on the principles of robust statistical estimation. This algorithm is intended to be used for processing of a massive data base that includes single-epoch astrometric measurements of thousands of radio-emitting quasars with the geodetic Very Long Baseline Interferometry (VLBI) world-wide facility. The system of accurate positions of these sources constitutes the fundamental International Celestial Reference Frame (ICRF3, Charlot et al. 2020), which underpins all other derivative celestial and geodetic reference frames. The astrometric stability of the most frequently observed quasars is of crucial importance for the overall accuracy and stability of ICRF3. We therefore develop a method to determine if some of the ICRF3 sources manifest periodic signals in their celestial positions, which could emerge from dual orbiting black holes in their centers, as well as a number of other effects in the extended structures and jets (Makarov et al. 2012).

The need for robust statistical estimation techniques generally arises in astronomical data processing when the available data are ridden with a large fraction of outliers outside of the commonly assumed Gaussian distribution of errors, representing a heavy-tailed sample distribution. Examples of critically important applications can be found in the mutual orientation alignment of different celestial reference frames (Malkin 2021; Frouard 2023; Lambert & Malkin 2023), where common object show a high rate of position offsets with extremely low formal probabilities.

In the most general setup of the problem, our task is to mathematically analyze a given time series (observations) d(t_i ), where periodic sinusoidal signals may be hidden. The data is discretized on a sequence of specific times of measurement t_i , i = 1, 2,..., N. When the measurements are taken on a regular, equally spaced grid with time step Δt, the problem is solved by the direct Fourier transform and subsequent computation of the Fourier power spectrum. The spectrum is quantified on a grid of angular frequencies f_k = 2 π/(kΔt), k = 2, 3,..., N, where the highest non-degenerate frequency f₂ = π/Δt is the angular Nyquist frequency. This set of frequencies is complete, because all other signals within the functional space spanned by the Fourier basis functions are not independent. In other words, the fitting function

is the exact and unique representation of any sequence d_i with a zero mean. This is no longer true if the cadence {t_i } is irregular. The Fourier harmonics are not orthogonal if sampled on an irregular cadence. In principle, this difficulty can be bypassed by constructing an ad hoc orthogonal basis from the Fourier harmonics using, for example, the Gram-Schmidt process, but the practical value of such representation is dubious, because the emerging fitting functions do not find a simple interpretation. However, we can disregard the issue of nonorthogonality and seek a solution to Equation (1) for a chosen f_k . Under this generalization, the classical periodogram analysis is equivalent to the LS fitting of model (1) for a grid of trial periods p_k = 2 π/f_k (Scargle 1982). The emerging LS problems, for a specific p_k , can be written as

where the design matrix A has two columns with calculated sequences of and values, and the number of its rows is equal to the number of data points N. The right-hand part vector d is the vector of centralized observations, and the vector of unknowns x comprises the two coefficients c_k and s_k . Standard LS algorithms can quickly solve this system to obtain the solution vector

Any standard LS algorithms can be used to solve this system. This has to be done for each trial period p_k ≡ 2π/f_k separately, including the setup of the design matrix A . In the example used in our paper, the number of trial periods N_p = 1000, but in other applications, it can be up to O(10⁵). This is not a problem for modern computers and LS algorithms, but in the 1970-ies, when the periodogram method started to attract astronomers' attention, the speed of computation was a crucial consideration. This motivated Lomb (1976) to propose a modification where the design matrix A is orthogonalized by introducing a phase shift τ_k in each of the fitting functions in (1), to the effect that the normal matrix A ^T A becomes diagonal. This modification has little practical advantage now but it brings in additional restrictions precluding necessary extensions of the model, as we will now discuss. Lomb's modification is therefore not recommended. A more accommodating and rigorous way of orthogonalization, if such a technical action is deemed desirable, was proposed by Ferraz-Mello (1981). Scargle (1982) also mentions that the spectral power of a pure noise data has a more predictable statistical distribution, but this argument is only valid for Gaussian noise with equal variances, which is never the case in practice.

Astronomical data series d often include non-periodic components. Using the archetypical example of detection of exoplanet signals in precision radial velocity measurements of host stars, the expected additional terms include a constant offset (the systemic radial velocity of the exoplanet system) and possibly a linear trend from perspective acceleration or a distant binary companion. It is not recommended to estimate these terms separately and prior to the periodogram solution and subtract them from the original data, which, unfortunately, is often done in practice. The reason why this pre-processing leads to an error in the periodogram is that these terms are not orthogonal to the fitted sine functions on a non-uniform cadence of data points. Unlike the regular Fourier transform, the trial frequencies are not integer multiples (harmonics) of the time interval. If a sinusoidal signal is present in the data, its estimated amplitude or power will be affected by the biased estimate of the constant term. The only correct and consistent way of dealing with additional terms is to include them in the fitting model (Ferraz-Mello 1981; Cumming et al. 1999). For example, the fitting model suitable for exoplanet detection can be

for each trial period p_k . Note that a separate constant term x₀ and a linear slope x₁ are obtained for each trial period, and they are not equal for different trial periods. The variation of these terms with the trial period reflects the error introduced into the periodogram by subtracting the common terms a priori.

The periodogram estimate can be the amplitude of the fitted sinusoid

or the power

As discussed above, we are using the amplitude periodogram in this paper, which has a more intuitive interpretation as the amplitude of the periodic signal in the same units as the measurements. The linear condition equations still take the form (2), but the design matrix A now has four columns, and vector x includes four unknowns x₀, x₁, x₂, and x₃. The solution vector is obtained from the LS solution, Equation (3). If a significant signal a(p_k ) is detected, the corresponding values x₀(p_k ), x₁(p_k ) provide the best estimates of the constant term and the linear trend.

What is the confidence level of a detected signal in the LS periodogram? This is an estimate of crucial importance, because the probability of the null hypothesis (that the detected feature is just a random fluke), also known as the false alarm probability (FAP), determines if we can believe the result. Traditionally, a high formal confidence is desired in astronomical applications such as detection of exoplanet signals, the recommended value being 0.997 (the 3σ level for a normal distribution). A robust method of estimating the FAP is the bootstrap simulation, which is also extendable to non-Gaussian distributions of measurement error. This is the method of "last resort" when the signal-to-noise ratio (S/N) of the detected signal leaves room for a catastrophic false positive. It is computationally expensive, however, and requires a sufficiently large number of data points. Monte Carlo methods, which are also computationally expensive, can be efficient when the statistical distribution of the observational noise is known. A random number generator is used to construct a sequence of synthetic measurements on the given sequence of times t_i , then a periodogram solution is obtained for each realization of noise, and the signal amplitude is computed. Repeating this process multiple times (O(10³) is usually required for an accurate estimation) allows us to estimate the cumulative distribution function (CDF) of the posterior distribution of the periodogram amplitude at any trial period, and hence, the p-value of the null hypothesis (or FAP).

Here we describe a computationally efficient and direct method of confidence estimation for LS periodograms in the extended form (Equation (4)). Noting that the standard periodogram solution obtained from the LS adjustment per Equation (3) is already based on the assumption that the measurement error is normally distributed, a direct computation of the periodogram CDF can be performed. If the covariance matrix C _d ≡ E[ d · d ^T ] is known or assumed, the corresponding covariance of the solution vector is

where . It is often assumed that the measurements d are statistically independent, in which case the matrix C _d is diagonal. If the errors also have the same variance σ², this equation further simplifies to C _x = σ² C .

The solution covariance C _x is a 4 × 4 symmetric matrix, which is computed for each trial period p_k . We are mostly interested in the a(p_k ) statistics per Equation (5). The two involved statistics x₂ and x₃, in accordance with the assumed normal distribution for each data point, are binormal variates, whose covariance matrix C _a is the corresponding 2 × 2 block of C _x . The estimated vector can be standardized to obtain

so that is a binormal uncorrelated variate of unit variance. This is equivalent to determining the error ellipse for binormal variates. The components of can be interpreted as the upper and lower S/N of the given periodogram result. Consequently, the quadratic form

is a χ²-distributed variate with 2 degrees of freedom. The corresponding confidence of rejecting the null hypothesis can be computed from the CDF of the distribution χ²[2] for each periodogram value. For a graphical representation of perodogram results, it is convenient to compare the confidence levels to specific points, which correspond to the ±1σ, ±2σ, and ±3σ intervals of the normal distribution, which have the cumulative probabilities of 0.683, 0.955, and 0.997, respectively. The corresponding levels of ψ (computed as CDF) are 2.296, 6.180, and 11.829. Periodogram amplitudes with ψ-values above 11.829 can then be regarded as highly confident positive detections at a confidence level above 0.997.

In the exoplanet detection literature, an alternative method of FAP-estimation is often used, developed by Baluev (2008). It is also based on the assumption that the signal contains only a finite set of model (base) functions, and the random component of the data vector is pure Gaussian noise with the known standard deviations σ_i . The statistical significance of a single periodogram value can then be naturally estimated from the properly normalized difference of the reduced χ² statistic of residuals with and without the corresponding harmonic terms (e.g., Equation (6) in Cumming et al. 1999), which follows the F-distribution (or beta-distribution if two or more specific periodogram frequencies are considered). However, while we have in practice a large number of periodogram value realizations, only the maximum value and the corresponding trial period are of interest. Even in the absence of detectable signal, given a large number of trials, it is probable that the highest significance value exceeds the threshold confidence level. This probability can be estimated within the extreme value statistic of an F-distributed hom*oscedastic random process assuming that the periodogram values are independent. This method is vulnerable to aliasing, which is caused by the limited spectral window of the given data series. The adjacent periodogram values are not independent, and periodogram features become increasingly wider "window functions" toward the longest trial periods. Non-uniform cadences with long gaps can also generate aliasing, spectral leakage, and spurious periodogram peaks. The extreme-value distribution fitting method (Süveges 2014) is more general for non-Gaussian processes, but it still refers to the null hypothesis of uncorrelated white noise in the data, which is inaccurate for the specific applications in this paper, or the spectroscopic detection of exoplanets (Makarov et al. 2009).

Astronomical time series often have unequal formal errors of individual data points. The formal error represents the expected standard deviation of the measurement, which can vary in a wide range because of observational conditions, instrument setup, etc. The LS solution in Equation (3), on the other hand, is unweighted, because it does not involve the estimated formal errors. The standard way of dealing with processing data of non-uniform precision is to use weighted LS fitting. It can be applied to LS periodogram analysis too in the framework of weighted periodogram solution. The basic equation replacing (2) becomes

where the weight matrix . The covariance of the right-hand part is now the identity. The formal covariance of the periodogram coefficients of interest transforms from Equation (7) into

The subsequent analysis of periodogram uncertainties is the same as described in Section 4.

We performed limited experiments using the formal weights on the example described in Section 6 to estimate the impact of this additional modification. We found rather limited changes in the computed periodogram amplitudes with the unweighted and weighted LS options. The most prominent features indicating possible signals have approximately the same shape and location. The greatest difference is found in the estimation of the 1σ and 3σ confidence intervals. The weighted covariance of the periodogram coefficients C _a is generally much smaller for the weighted solution than for the unweighted solution. This is caused by a large spread of individual formal errors, and the fact that the weighted LS solution is optimal. If all the formal errors are equal, the periodogram covariances and the derived amplitudes become equivalent in the two solutions. Thus, the weighted covariance C _a is the global minimum of all possible unweighted counterparts. The lower covariances result in narrower confidence intervals, and the net result is that most of the periodogram solution becomes a highly confident positive detection. This result is completely misleading for the given example, because, as we will see in the next Section, the formal errors of the data points have little bearing on the actual dispersion and statistical distribution of the data.

Let us summarize the implicit assumptions involved in the LS periodogram method.

• 1.

  The data vector is a composition of random uncorrelated noise and a single monochromatic sinusoidal signal, whose amplitude and period are to be determined.

• 2.

  The measurement noise is Gaussian.

• 3.

  The data sequence is centralized, i.e., has a zero mean—unless the extended version of the method is employed.

An example when the second assumption is violated can be found in (Makarov et al. 2010). The astrometric position (photocenter) of the Sun as measured by a distant observer is subject to stochastic variations caused by the presence of sunspot groups and bright plage areas on the rotating surface. Each photometric feature generates a time-variable shift of the unresolved disk on the timescale of days, which is not periodic because of the phase scrambling. The composition of such stochastic signals is an unpredictable "jitter." The measured shifts from the mean photocenter show an utterly non-Gaussian distribution because the intrinsic distributions of the sunspot sizes, lifetimes, and positions within the disk are not normal. In this case, the nominal LS periodogram, as well as the traditional FAP estimation, are likely to produce misleading and inaccurate results.

Figure 1 shows the observed time series used in this paper to illustrate the application of the proposed 1-norm periodogram analysis. It shows the high-accuracy astrometric data collected by geodetic VLBI for the ICRF3 source IERS B0642+449 for nearly 40 yr of continuous observations. Each data point corresponds to a one-day "session" with multiple delay measurements of this source together with a number of other ICRF3 sources. The observational data are represented as coordinate offsets (left panel) and y = δ_obs − δ_mean (right panel) in mas, where {α_mean, δ_mean} are the weighted mean coordinates for this source in the equatorial coordinate system. The formal errors for each observation are shown as error bars. This enigmatic high-redshift (z = 3.41) gamma-ray blazar is obviously one of the astrometrically unstable ICRF3 sources with shifting position mostly in the R.A. component. The origin of the position variations is outside of the topic of this paper, but we briefly note the study by Xu et al. (2016), who detected a dual structure of IERS B0642+449 from a closed delay analysis of a high-intensity geodetic VLBI session. The detected separation of the dual components is approximately 0.46 mas and the position angle is 2622. One of the interesting applications of sub-mas astrometry with VLBI is the possibility of detection of orbiting dual AGNs. With an estimated scale of 7.6 pc mas⁻¹ at this redshift, a binary black hole with a total mass of 10¹⁰ M_☉ and a period of 40 yr may have an angular separation of about 10 μas, which may be within the reach with this type of data. Central engine binarity may be one of the explanations for the observed quasi-periodic modulation of some gamma-ray blazars' light curves (Ackermann et al. 2015).

Both coordinate trajectories in Figure 1 appear to include long-term variations and, possibly, periodic components on the timescale of a few hundred days. Are they statistically significant? We begin with the standard LS periodogram analysis using the extended model Equation (4). The need to include the linear terms, in particular, comes from the possibility of a "secular" proper motion in the data, which is not part of the astrometric model used in the VLBI data reduction pipeline. We compute the periodogram fitting coefficients {x₀,..., x₃} for an exponential grid of 1000 trial periods, p_k = p₀ dex(q k), with the exponent step q = (3652.5 − p₀)/1000 in days. The longest trial period is then 10 yr, which is practically limited by the time span of the available data.

The results are shown in the upper row plots of Figure 2 for the two coordinate components. The periodogram amplitude estimates are connected with a black line to aid the eye. The significance of each periodogram point is also computed for this unweighted solution according to Section 4. We color-coded the significance by the normalized confidence level, so that estimates below the 1σ level are marked with blue dots, and estimates above the 3σ level are marked with red dots. A large number of values appear to be highly significant with periods across the entire range, including some short periods below 100 d, which obviously cannot be physical. This result, with a jungle of sharp peaks in the short-period domain and a few prominent features in the long-period domain, is typical of LS periodograms for "noisy" data. The inference is completely false, and we will now reveal why that happens.

The single-epoch positions measured with VLBI are two-dimensional, and each position determination {x, y} comes with a formal covariance G , which is a 2 by 2 matrix. It is convenient to consider the normalized and centralized single-epoch position offset

because it is a scalar variate, which is expected to follow a Rayleigh distribution with scale 1, reducing the dimensionality of statistical analysis to 1. The true coordinates are not known, but they can be separately estimated as the weighted mean position. Note that even the well-known formula for covariance G is based on the underlying assumption of normally distributed random errors. The variate D allows us to test this basic assumption. The histogram of D values computed for the example data set in Figure 1 is shown in Figure 3. For reference, the expected Rayleigh[1] distribution (normalized to the same area) is shown with the blue line. We can see that the actual distribution of astrometric offsets is very far from the expectation, and the difference cannot be fixed just by scaling the formal errors. Although the mode of the empirical distribution is approximately where it is expected to be (at 1), a long and powerful tail stretching far beyond the Rayleigh[1] curve indicates that nearly half of the available measurements have values associated with nil probabilities of occurrence within the assumed statistical model.

The heavy-tailed nature of the data distribution invalidates the LS periodogram method. The data points with large deviations from the mean should not be called outliers in this case, because they represent a large part, if not the majority, of the population. Simple fixes such as clipping the data outside the 3σ threshold are not justified. The numerous deviant data points corrupt any LS estimation, and generate bogus signals in this periodogram analysis. Methods of robust estimation are designed to handle heavy-tailed data in a more consistent way. In particular, the 1-norm estimation seeks to minimize the sum of absolute values of residuals rather than the sum of their squares:

This merit function diminishes the impact of large deviants and permits a meaningful solution for any intrinsically symmetric populations. It is robust with respect to the subsample containing high normalized offsets, because each individual data point has an effectively lower weight in the solution, irrespective of its value.

The same periodogram models (Equations (1) and (4)) can be used as in the classical LS method. The main differences in implementation are of the technical character. The main optimization problem is no longer linear, and it cannot be formalized as Equation (2). Consequently, there is no direct calculation of the associated covariance matrix of the periodogram coefficients. This can still be done numerically using Monte Carlo simulations. The solution itself is implemented with one of the existing global nonlinear optimization methods with vector-valued arguments, such as the Nelder-Mead (simplex downhill), differential evolution, or simulated annealing methods. ³ These methods are computationally much more expensive than the regular LS periodogram. However, we achieved a computing time of about 1 minute on a regular laptop for the given example with 1668 data points and 1000 trial periods for the two time series.

The resulting 1-norm periodograms for the given data sets are shown in Figure 2, lower row. They are expressed in the same values (amplitudes, per Equation (5)) and units as the LS periodograms in the upper row, so that they can be directly compared. Formal confidence levels cannot be directly computed for the 1-norm solutions, because the population distribution is non-Gaussian. We reproduce, however, the re-normalized confidence intervals 1σ and 3σ from the LS solution to emphasize the significance of the results.

Quite clearly, the robust 1-norm periodograms paint a different picture about the temporal variations of the given data. The amplitude values dropped by half or more, and most of the estimated values are now below the 1σ-interval. The largest reduction is seen in the high-frequency domain. Given the nature of the object under investigation, the low-frequency features are of special interest. We find the main features at different locations than with the LS method, and for the R.A. component, they clearly dominate the spectral power distribution. Intriguingly, there is a compact location around 1730 days with periodogram amplitudes above 3σ, which was completely insignificant in the LS solution. To test if this point is indeed associated with a high level of confidence, a much more extensive bootstrapping or Monte Carlo simulations are required.

We performed a non-parametric bootstrap by producing 100 data samples from the original time series by randomly permuting its elements but keeping the fixed cadence of epochs, thus inheriting the same distribution of the uncorrelated noise component as the original data. Then, for each data sample, we computed the LS and 1-norm periodograms. The bootstrap distribution appears reasonably symmetric. For each period, the N% confidence interval (0 ≤ N ≤ 100) is given by the non-parametric percentile bootstrap interval contained between the (N/2) th and (100 − N/2) th percentiles. For verification purposes, the periodogram analysis and the bootstrapping estimation was independently implemented by two authors using different computer languages, on the same data set. Figure 4 displays the obtained results where the grey lines represent the difference between the upper and lower bounds of the confidence intervals, respectively for N = 68.3% (dotted line), N = 95.5% (dashed line), and N = 99.7% (solid line). In the case of the 2-norm implementation, the bootstrap provides a confidence interval consistent with results shown in Figure 2. For the 1-norm spectrum, the confidence interval is lower than for the 2-norm, reflecting the robustness of the 1-norm estimation, which is less sensitive to the data points in the extended tail of the distribution. The amplitudes reach highest significant values, with a significance level well higher than 99.7%, for the right ascension suggesting that the bump observed at a period of ∼1700 days is not due to chance and the secondary peak at ∼2800 days should also be considered for further investigations. We note, however, that the bootstrap-estimated confidence only refers to the random uncorrelated noise in the data (of arbitrary PDF). The data may include a time-correlated component of physical or instrumental origin. Broad periodogram features seen in Figure 4 may require additional analysis using, e.g., structure functions of first or second order, which also incorporate periodic components with time-variable phase (Rutman 1978; Simonetti et al. 1985; Rutman & Walls 1991).

Diurnal geodetic VLBI sessions have been regularly scheduled over nearly 40 yr, using networks of stations separated by baselines of hundreds to thousands of kilometers long. Within each daily session, a number of widely separated radio sources are observed multiple times over the course of 24 hours. The resulting data are processed in a few data analysis centers, including the U.S. Naval Observatory. In this paper, we use a global solution for two-dimensional coordinates of epoch calculated at USNO (2022a) in the standard S/X band setup. This data product includes astrometric time series from more than 6000 diurnal sessions. The total number of sources is 5153 in this data set, but here we only consider 259 of them with more than 200 single-epoch measurements.

The 1-norm periodogram computation was uniformly applied to each of the frequently observed ICRF3 sources, separately for R.A. and decl. offsets from the weighted mean positions. These mean positions are specifically computed for the given data set and the solution version, so they may slightly differ from the published ICRF3 mean positions. The purpose of this numerical experiment was to identify sources with possible sinusoidal variations in the observed positions on the sky. The results are presented in a compact form in Table 1 for 49 quasars where the formal significance criterion ψ > 11.829, from Equation (9), is triggered in either of the coordinates for at least one trial period. The table provides IERS names of the sources (which should be prepended with letter B to match Simbad identification), our computed mean R.A. and decl. coordinates in degrees, the total number of diurnal sessions, and the significant trial periods.

The main result of this computation is that there seems to be no isolated single-frequency sinusoidal signals in the observational data similar to those that are found for astrometric binaries. Instead of well defined peaks in the periodograms, we find "packages" of trial periods with elevated amplitudes and significance levels. The emerging picture is more consistent with an ensemble of small vaguely periodic modes of non-commensurate frequencies. One-fifth of the detections have a rising periodogram amplitude toward the upper limit of this analysis (10 yr). The significant periods are mostly longer than 2 yr, although much shorter trial periods have been tested. Only 19% of the sample show periodic variations above the 3σ level. The typical peak amplitudes are in the range 100–200 μas.

We note that the robust periodograms were computed separately for the R.A. and decl. coordinates of epoch positions. In the case of Keplerian motion in a binary system, the detectable signal may be present in both coordinates with the same principal period (and its harmonics) but with different phase and amplitude. The probability density of the angle i between the line of sight and the vector of orbital angular momentum is proportional to . Therefore, nearly face-on projected orbits are less likely than nearly edge-on orbits. For marginally detectable trajectories, the detactable signal is mostly present in one dimension, which is uniformly distributed with respect to the local north direction. The largest extent of the projected orbit can be aligned with one of the coordinate axes with the same probability as a tilt of 45° or 135°. In the latter case, the detectable signal is split between the coordinates, and it should be harder to find it with confidence. Possible ways to deal with this problem include rotating the R.A.-decl. measurements on a grid of position angles to find a preferred direction maximizing the signal amplitude from a 1D periodogram. Technically, if a significant single-period signal is detected in both R.A. and decl. coordinates, a 2D version of robust periodogram can be implemented. The available coordinate measurements are combined in a single LS adjustment, but the number of unknown terms per trial period increases to a minimum of 8 because of the unknown phase. This may erode the confidence level of the signal, if the projected orbit is strongly elongated due to the geometric orientation or large eccentricity. Quasars IERS B1451−375, 2234+282, and 2318+049 are attractive targets for further investigation, because they show coherent periodicities in both coordinates from our results in Table 1.

We have shown in this paper that the classical LS periodogram method is firmly based on strong and restrictive assumptions about the distribution of post-fit residuals (which is assumed to be Normal) and the character of physical signals in the data. It provides an optimal, unbiased, and unique solution for periodogram power or amplitude only under these conditions. Whenever the sample distribution shows significant departures from the Gaussian PDF, or more complex signals are present that are not captured in the model, the LS method becomes corrupted and can produce absolutely misleading results.

We have considered a specific observational data set for a moderately variable ICRF3 source collected over >30 yr by the global geodetic VLBI system. The distribution of astrometric positional offsets with respect to the mean position on the sky is explicitly non-Gaussian when scaled with the given 2D formal covariances or in absolute values. The normalized offsets are well represented by a log-normal distribution with a tighter mode and a heavy tail extending to high values. Nearly half of the measurements are way outside of the expected distribution. As a result, the traditional 2-norm (LS) periodogram produces a complex structure with multiple features that are formally above the 3σ confidence interval across the spectrum of trial periods. This result is completely bogus. The robust 1-norm periodogram method, when applied to the same data, produces amplitudes that are smaller by half or more. Ranked by the same previously estimated single-point confidence, the 1-norm values are all insignificant except for a single point in the R.A. component with a period of 1730 days and amplitude 72 μas, which appears to be above 3σ. Is this periodic signal real? The best way to find out is to continue taking high-precision measurements of this source with VLBI for a few yr. A stable sinusoidal signal, which could be produced by an orbiting binary black hole, for example, would emerge more strongly on the longer timescale. Alternatively, physical models could be tested, where transient periodic signals wax and wane in segments of the data due to phase scrambling. This new method provides the opportunity to more reliably and extensively search for periodic signals in non-Gaussian time series at the margin of available accuracy.

The 1-norm periodogram computation was performed for 259 ICRF3 sources with more than 200 diurnal sessions collected over nearly 40 yr. These measurements are characterized by heavy-tailed sample distributions of residuals. We identified 49 objects (19%), which have at least one statistically significant periodogram value in either coordinate component. Short periods are never found, indicating a possible physical mechanism of these signals in the transient structure of the radio-emitting sources. The signals are not consistent with clean sinusoidal variation at a specific frequency, which would emerge for an orbiting binary black hole. Rather, the pattern is that of "vague periodicity" represented by packages of sine waves with a distribution of frequencies. A possible physical model is a source that moves in loops on the sky returning to the vicinity of the initial position after some characteristic time, which may also vary with time. The estimated amplitude of these vaguely periodic excursions is 70 μas and higher. About one-fifth of the detected signals are truncated by the upper boundary of our periodograms (10 yr). Further investigation of these astrometric wobbles and continuous daily measurements will refine the models and allow us to understand the nature of the phenomenon.

Facility: VLBI - .

Software: Wolfram Mathematica.