: Continuous waves, LIGO, Hardware injections, Sensitivity analysis, Signal calibration, Template Generation Method, F-Statistic, Bayesian

1 Introduction↩︎

Since the first detection of gravitational waves from merging black holes in September 2015 [1] by the LIGO detector, there have been 90 published detections of compact star mergers by the LIGO-Virgo-KAGRA (LVK) Collaboration from the first three observing runs (O1–O3) [2]. In addition, there were nearly 300 public alerts [3] of LVK merger detections in low latency by real-time analyses during the recently completed O4 run, which began May 24, 2023. Despite this success in detection of transient signals, no persistent gravitational-wave signals, continuous-wave (CW) or stochastic, have been detected.

CW gravitational signals are expected to be emitted by rapidly spinning neutron stars with a small non-axisymmetry (mass or mass-current quadrupole asymmetry), but with signal amplitudes at the Earth several orders of magnitude weaker than the amplitudes of transient signals detected to date. Detecting such weak signals requires integration over long time intervals, with the best possible sensitivity accessible via fully coherent integrations over full observing spans [4].

Reliably extracting such signals requires an accurate understanding of detector response over the observing time, including absolute phase response over as many as 100 billion signal cycles. To help validate that understanding, the LIGO Scientific Collaboration uses a “hardware injection” system [5] that allows direct simulation of gravitational-wave signals via modulation of the positions of the interferometer mirrors¹.

This direct modulation is carried out via a “photon calibrator” consisting of an auxiliary laser directed at one of the interferometer mirrors at a known incident angle and known intensity, which is modulated quasi-sinusoidally [6]. The actuation response function is relatively stable (approximately that of a free mass for frequencies well above the mirror’s pendulum frequency) once the absolute scale factor and phase dependence are determined through meticulous pre-run calibration. This stability makes the actuation nearly independent of the bulk of rest of the LIGO calibration process, which takes into account frequency-dependent and time-dependent effects of the detector response function.

Hence the monitoring of CW injections, as measured in the gravitational-wave data channel constructed via that detector response, provides a key validation of the calibration. In addition, because additional cleaning of the data channel for both transient instrumental glitches and instrumental spectral lines is often carried out to improve sensitivity to astrophysical signals, the hardware injection recovery can verify the cleaning has done no harm, and where expected, has indeed improved sensitivity.

A secondary benefit of the hardware injections (but not one that by itself justifies the injection campaign) is allowing validation of search programs, one that is independent of the varying methods applied by different CW search programs for efficient large-scale “software injections” used in tuning algorithmic parameter choices and in determining signal sensitivities [7], [8].

This paper describes three different monitoring programs used both for “real time” monitoring (cadence: once per day) to catch potential unintended detector response changes and for longer-term analysis and validation of data cleaning techniques. The first method uses exact templates for the intended injections, allowing only amplitude and phase offsets to vary. The second method uses the venerable \(\mathcal{F}\)-statisticdetection statistic, which serves as the foundation of a multitude of CW searches for different potential sources [4] and allows a frequentist estimate of amplitude, phase offset and two angles describing the orientation of the simulated star. The third method applies a Bayesian statistic used in searches for CW radiation from known pulsars [9]. Collectively, these programs have helped validate both the nominal calibration of the gravitational strain channel and cleaning techniques used to remove deterministic instrumental contaminations during the O4 observing run,

1.1 Continuous Gravitational Waves↩︎

Continuous gravitational waves represent a distinct class of gravitational wave signals [4]. This radiation, expected from rapidly rotating galactic neutron stars, arises from slight stellar asymmetry in mass quadrupole moment or mass current quadrupole moment [10]–[19]. Those asymmetries lead to much weaker signal amplitudes, however, than what has been observed from compact binary mergers, despite the much larger distances of the merger sources. Long integrations (months to years) are expected to be necessary to detect CW sources, with the best sensitivity requiring phase agreement between the actual signal and the model used to extract the signal.

For a mass quadrupole moment asymmetry, the amplitude \(h_0\) of the continuous wave generated by a spinning neutron star is related to its equatorial ellipticity \(\epsilon\) and its distance \(d\) from Earth. The equation for the amplitude can be derived to be: [4]:

\[h_0 = \frac{16\pi^2 G \epsilon I_{zz} f^2}{c^4 d},\]

where \(G\) is the gravitational constant, \(c\) is the speed of light, \(I_{zz}\) is the moment of inertia of the neutron star about its rotational axis \(z\), and \(f\) is the rotational frequency of the neutron star.

The equatorial ellipticity \(\epsilon\), is defined as:

\[\epsilon = \frac{|I_{xx} - I_{yy}|}{I_{zz}}\]

where \(I_{xx}\) are \(I_{yy}\) are the other two moments of inertia along the principal axes of the neutron star.

Based on the absence of detections from previous searches and consistent with astrophysical knowledge, we expect plausible long-lived continuous gravitational waves to have a strain less than \(10^{-24}\) and to originate from galactic sources, at distances of O(kpc)[20], [21]. It is only the long-lived nature of these waves that offers hope for their detection as we can integrate data over long durations, to produce a statistically significant Signal-to-Noise Ratio (SNR). For a coherent integration and stationary detector noise, the SNR is expected to increase proportionally to the square root of the observation time as:

\[\text{SNR} \quad \propto \quad {h_0 \sqrt\frac{T_{\rm obs}}{S(f)}}\]

where \(h_0\) is the strain amplitude, \(T_{\rm obs}\) is the total observation time, and \(S(f)\) represents the power spectral noise density [22], [23] at the frequency \(f\) of the signal. As energy and angular momentum are lost by the star to gravitational waves, the spin frequency and the CW signal frequency (double the spin frequency) should decrease slowly, leading to a small negative first derivative \(df/dt\). In the LIGO injections described in the next section, only the first derivative is simulated, for simplicity, and signal amplitude is held constant within a given observing run, although a true CW source would have higher-order frequency derivatives that become important over long integration times and would have slight amplitude reduction as the spin frequency declines [4].

The eventual detection of continuous gravitational waves from a galactic neutron star should yield insight into the neutron star equation of state, and in the event of successful identification of electromagnetic radiation from the same star, should provide a splendid multimessenger source, one that might well be studied for years and decades to come, to understand both the dynamics of these exotic objects and to probe fundamental physics in a regime inaccessible to laboratory experiments [24].

1.2 LIGO and its Hardware Injection Monitoring: an overview↩︎

Hardware injections [5] make up one element of the LIGO detector’s [25]–[29] calibration and validation framework [30], [31], confirming our ability to identify and reconstruct gravitational wave signals. These injections use physical modulations of the detector’s test mass mirrors to simulate the effects of a gravitational wave, introducing controlled signals into the interferometric data stream. The objectives of hardware injections are to verify the fidelity of signal reconstruction pipelines, ensure the accuracy of amplitude and phase calibrations, and evaluate the robustness of data analysis methods under varying signal conditions [5].

Continuous wave injections pose a technical challenge because of their extended duration [4] and the need to maintain absolute phase coherence. To address this challenge, the LIGO injection system uses the same barycentering software [32] used in published CW searches and in software injections, along with a data buffering mechanism [5] to ensure smooth, continuous insertion of signal waveforms without unintended discontinuities in injected waveforms. Because the CW injections are extremely narrowband, it has proven convenient to use a frequency-domain technique in which a Fourier transform of injected strain is computed for short intervals of data with a frequency-domain inverse actuation response function applied to the transform before inverting to create a time-domain waveform for injection into the photon calibrator system. The inverse actuation function includes a correction for time delays in the injection system itself. Although the software infrastructure exists for direct time-domain injection of CW signals [5], it was found in the O1 observing run that such injection of CW signals leads to high-frequency noise injection.

Daily reconstructions of hardware injections, based on the previous day’s data and based on the cumulative data, are carried out to ensure there have been no significant unintended changes in detector response. This check is most important each week, following regular Tuesday maintenance periods. Graphs are created and inspected to track trends in amplitude recovery, phase stability, and the expected growth in SNR over time. Additionally, cumulative analyses over the observation period to date allow tracking of detector stability, aiding in identifying drifts, systematic biases or spectral line contaminations that can affect sensitivity [33].

In CW searches, the hardware injections provide a set of sources for validating analysis programs. Because most CW search programs are designed to analyze data stretches of months to years in duration using narrowband spectral bands at a time, it is typically not feasible to carry out software injection campaigns by superposing time-domain strain waveforms on the data and then computing Fourier transforms. Instead, more computationally efficient techniques, tailored to each search program, are used in which Fourier coefficients or derived quantities are superposed. The associated variation in the technical implementations of injections used for tuning and validating different pipelines leaves open the possibility of inconsistency across different programs. The existence of a set of hardware injections already present in the data means that performance can be evaluated and compared without ambiguity across different programs. Validation of results with hardware injections is a routine requirement in LIGO-Virgo-KAGRA CW analyses.

As noted above, however, despite this utility, one would not implement hardware injections for this purpose alone. One could, in principle, simply create a data set copy with software injections to allow unambiguous cross-pipeline comparisons. The primary benefit and justification for hardware injections is to validate assumptions about detector response, including the effects of calibration and data cleaning, especially the self-gating of large transients [34], [35], to provide a true end-to-end validation of the search.

2 LIGO CW Hardware Injections↩︎

The CW injections are meant to simulate neutron stars in both isolated and binary systems [32]. These signals span a range of frequencies and amplitudes, ensuring coverage of the detector’s operational bandwidth. Key parameters such as nominal source frequency (at a particular reference time), strain amplitude, phase, and polarization are defined, enabling comparisons between the intended injected signals and the corresponding reconstructed values [5]. Below we describe the details of the signal model and the parameter choices made for the 18 different stars simulated during the O4 observing run.

2.1 Signal Model for Continuous Wave Detection↩︎

Carrying out CW searches in LIGO data with optimum sensitivity requires correcting for signal frequency and amplitude modulations. The frequency (or phase) modulations result from the detector’s motion with respect to an approximately inertial frame of reference (taken to be the solar system barycenter). For a star in a binary system, there are also frequency modulations from the orbital motion of the source. Amplitude modulations arise from the changing orientations of the LIGO interferometers with respect to the line of sight to the star. The relative frequency modulations due to the Earth’s daily rotation are O(10\(^{-6}\)), and those due to the Earth annual orbit around the solar system barycenter (SSB) are O(10\(^{-4}\)). The amplitude modulations are O(1) and have periodicities of 12 hours and 24 hours (sidereal) [23]. Here we provide more detail on some of these modulations in the form of the signal model. The barycentering code used to account for these modulations is based on the same model used in the pulsar astronomy community [36] and includes General Relativistic effects.

We treat a continuous gravitational wave emitter as an isolated, spinning, rigid, triaxial ellipsoid. The strain waveform perceived by an interferometer can be modeled as:

\[h(t) = F_+ (t; \psi) h_+ \cos \Phi(t) + F_\times (t; \psi) h_\times \sin \Phi(t),\]

where \(h_+\) and \(h_\times\) are the amplitudes of two orthogonal, linearly polarized signals defined with respect to the spin axis of the star:

\[\label{eqn:polarizations} h_+ = h_0 \frac{1 + \cos^2 \iota}{2}; \qquad {\rm and} \qquad h_\times = h_0 \cos \iota,\tag{1}\]

where \(h_0\) is the amplitude of the gravitational wave, \(\psi\) denotes the polarization angle, \(\iota\) is the inclination angle of the neutron star’s spin axis relative to the line of sight, and \(\Phi(t)\) describes the time-dependent phase of the signal [23]. The detector response functions, \(F_+ (t; \psi)\) and \(F_\times (t; \psi)\), capture the directional sensitivity of the interferometer to the plus and cross polarizations of the gravitational wave [23].

The phase evolution, \(\Phi(t)\), plays a pivotal role in accurate signal modeling. It can be expressed as:

\[\label{eq:phase95evolution} \Phi(t) = \Phi_0 + 2\pi \sum_{m=0}^N \frac{f^{(m)}(\tau_{\text{ref}}) \, \tau(t)^{m+1}}{(m+1)!}\tag{2}\]

where \(\tau(t)\) is the SSB arrival time of the signal, which includes relativistic effects such as Shapiro delay and Einstein delay. The term \(\tau_{\text{ref}}\) represents a reference time, at which the frequency \(f\) and its higher derivatives \(f^{(m)}\) are defined. The precision of this phase model is crucial for retaining coherence over long observation periods.

Modulations arise on daily and annual timescales due to the Earth’s rotation and orbital motion around the Sun. These effects can be described as:

\[\tau(t) = t + \frac{\mathbf{r}(t) \cdot \hat{\mathbf{n}}}{c} - \tau_{\text{ref}}\]

where \(\mathbf{r}(t)\) is the position vector of the detector relative to the SSB, \(\hat{\mathbf{n}}\) is the unit vector pointing from the SSB toward the source, and \(c\) is the speed of light.

2.2 Parameters of O4 Injections↩︎

A campaign of 18 CW hardware injections (indexed from 0 to 17) has been carried out during the O4 run, comprising 16 isolated and 2 binary systems. These injections have the same source parameters as those used in previous observing runs [21], except for amplitudes which have been reduced in keeping with frequency-dependent reductions in detector noise levels. Six of the injections spread across the band are given a high enough amplitude to be detected in coherent integration of a single day of data. These six loud injections are used for verifying on a daily basis that no unintended large changes have occurred in the calibration, such as might occur from the addition of an inadvertent time delay or even sign flip in the digitally controlled electronics.

The sources selected for the CW O4 hardware injections span a frequency range of 12 to 2991 Hz. The O4 campaign includes 6 loud simulations (injections 1, 4, 5, 6, 10 and 14), i.e., having relatively high SNRs for daily monitoring and 12 weaker injections for longer-term monitoring at other nominal frequencies. Table ¿tbl:tab:parameters1? lists the parameters of all 18 injections. Additional orbital parameters for the two binary simulations (injections 16 and 17) are listed in Table ¿tbl:tab:parameters2?.

3 Methodology↩︎

The methods used to monitor hardware injection analysis are presented in the following subsections: a template generation method, a \(\mathcal{F}\)-statisticapproach, and a Bayesian inference. Each method has specific strengths, to be discussed below.

3.1 Template Generation Method↩︎

The template generation method of injection recovery exploits the expected phase and amplitude modulations of the signal and a matched filter approach to extract estimates of true injection amplitude and phase offset. Although all three monitoring methods presented in this article include some degree of templating, since all assume the correct sky location and at least approximately correct frequency evolution, the template generation method does so to the greatest degree.

The observation span is split into segments of duration \(T_{\mathrm{SFT}}\) = 1800 s for which discrete Fourier transforms are computed, known as “short Fourier transforms” or SFTs. The SFTs for the data for each detector, used for spectral monitoring and astrophysical figures of merit, are already available at the start of each UTC day’s injection monitoring run. That run begins several hours after the corresponding day’s end, to allow for data transfer delays from the observatories to the Caltech computing cluster, where monitoring is run centrally. Parallel sets of narrowband SFTs for each day and for each injection are then computed for a pure signal template identical to what is supposed to have been injected in the data, using LIGO Analysis Library software [32], including a python-code wrapper [37]. The band of each SFT for a given injection is centered on the injected signal’s nominal detector-frame frequency at the start of that SFT’s time interval.

In the following, denote the \(k\)th data set SFT as \(d^k\) and the complex Fourier coefficient contained in its \(i\)th bin as \(d^k[i]\), with the corresponding signal template quantities \(s_{\rm T}^k\) and \(s_{\rm T}^k[i]\). The (real) powers in these bins are simply the absolute squares of these values: \(D^k[i]= |d^k[i]|^2\) and \(S_{\rm T}^k[i]= |s_{\rm T}^k[i]|^2\).

The analysis for each SFT \(k\) is performed in two steps. First, a discrete frequency band (selection window) is determined that captures at least 99% of the total signal template power. Second, an inner product of the data and the complex conjugate of the signal template bins within the selection window is computed, from which estimates of signal amplitude and phase offset can be computed, both for SFT \(k\) alone and for the cumulative set of SFTs from 1 to \(k\) in the observing run. Uncertainties and SNRs are also computed based on empirical estimates taken from neighboring frequency windows of the same width.

To determine the selection window for each SFT, the following algorithm is used. First, the template bin \(i_0\) with maximum expected signal power is found, and the frequency window of minimum / maximum indices \(i_{\rm min}\) / \(i_{\rm max}\) are initialized to \(i_{\rm min}=i_{\rm max}=i_0\). In an iterative process, \(i_{\rm min}\) and \(i_{\rm max}\) are adjusted until the fraction \(R_k\) of total window power is at least 99%of the total signal power in the SFT, that is, until

\[R_k= \frac{\sum_{i=i_{\rm min}}^{i_{\rm max}} S_{\rm T}^k[i]}{\sum_{{\rm all}\, i} S_{\rm T}^k[i]} \ge 0.99.\]

At each iteration \(i_{\rm min}\) is decremented by one or \(i_{\rm max}\) incremented by 1, depending on which bin yields a greater increase in signal template power. For later convenience, we also define a Fourier transform conversion factor \(\Lambda_k\):

\[\Lambda_k = \frac{\sum_{i=i_{\rm min}}^{i_{\rm max}} S_{\rm T}^k[i]}{h_{\rm Inj}^2}.\]

where \(h_{\rm Inj}\) is the template signal amplitude used in the time series from which \(s_{\rm T}^k[i]\) is computed. Once the selection window is defined for SFT \(k\), a (complex) product sum \(P[k]\) is defined from the inner product of the data with the template in that window:

\[P[k]= \sum_{i=i_{\rm min}}^{i_{\rm max}} d^k[i](s_{\rm T}^k[i])^*,\]

where, ideally, \(\Re\{P[k]\}\) contains the injected signal and has an expectation value proportional to \(h_0\), and the expectation value of \(\Im\{P[k]\}\) is zero. If, however, the signal injection has an incorrect phase offset (from imperfect actuation or from an imperfect calibration), then the expectation value of \(\Im\{P[k]\}\) can be non-zero.

To normalize these quantities to yield direct injection amplitude, we apply the scale factor \(\Lambda_k\) and correct for the template amplitude contribution to the product sum:

\[h_{\cos}[k]= \frac{\Re\{P[k]\}}{\Lambda_k h_{\rm Inj}};\]

\[h_{\sin}[k]= \frac{\Im\{P[k]\}}{\Lambda_k h_{\rm Inj}}.\]

to define”cosine” and “sine” projections of the data onto the template.

For each SFT \(k\), an estimate of the injection amplitude \(h_0(k)\) and phase mismatch \(\Delta\Phi_{\rm est}[k]\) is found from

\[h_{\rm est}[k]= \sqrt{(h_{\cos}[k])^2 + (h_{\sin}[k])^2};\]

\[\Delta\Phi_{\rm est}[k]= \arctan(h_{\sin}[k]/h_{\cos}[k])\]

where, ideally, \(\Delta\Phi_{\rm est}[k]\) has an expectation value of zero.

As an empirical estimate of the noise contamination, the same computation is carried out for additional, adjacent control windows of the same width as the selection window, with different bands of data multiplied against the same signal template. Specifically for a control index offset \(m\), the (complex) control product sum is \[P_{\rm ctrl}[m,k]= \sum_{i=i_{\rm min}}^{i_{\rm max}} d[n_{\rm win}m+i] (s_{\rm T}^k[i])^*,\]

where \(n_{\rm win}= i_{\rm max}-i_{\rm min}+1\), that is, the width in bins of the selection window.

We determine a sample of 10 control values using \(m = -5\ldots5\), excluding \(m=0\), and from them compute the variances of the cosine and sine components \(V_{\cos}[k]\equiv{\rm Var}\{\Re\{P_{\rm ctrl}[m,k]\}\}\) and \(V_{\sin}[k]\equiv{\rm Var}\{\Im\{P_{\rm ctrl}[m,k]\}\}\) separately. The variances are then used to compute weights:

\[w_{\cos}[k]= \frac{1}{V_{\cos}[k]}; \qquad w_{\sin}[k]= \frac{1}{V_{\sin}[k]}.\]

These weights are used, in turn, to compute cumulative best estimates of cosine and sine projections from all SFTs up to and including SFT \(k_{\rm max}\):

\[\hat{h}_{\cos}= \frac{\sum_{k=1}^{k_{\rm max}} w_{\cos}[k]h_{\cos}[k]}{\sum_{k=1}^{k_{\rm max}} w_{\cos}[k]};\]

\[\hat{h}_{\sin}= \frac{\sum_{k=1}^{k_{\rm max}} w_{\sin}[k]h_{\sin}[k]}{\sum_{k=1}^{k_{\rm max}} w_{\sin}[k]}.\]

which have estimated uncertainties:

\[\sigma^{\hat{h}}_{\cos}= \sqrt{\frac{1}{\sum_{k=1}^{k_{\rm max}} w_{\cos}[k]}};\]

\[\sigma^{\hat{h}}_{\sin}= \sqrt{\frac{1}{\sum_{k=1}^{k_{\rm max}} w_{\sin}[k]}};\]

As for individual SFTs, one can compute cumulative injection amplitude and phase mismatch from these quantities:

\[h_0 = \sqrt{(\hat{h}_{\cos})^2 + (\hat{h}_{\sin})^2};\]

\[\Delta\Phi= \arctan(\hat{h}_{\sin}/\hat{h}_{\cos}).\]

which have estimated uncertainties:

\[\sigma_{h_0} = \sqrt{\frac{(\hat{h}_{\cos})^2(\sigma^{\hat{h}}_{\cos})^2 + (\hat{h}_{\sin})^2(\sigma^{\hat{h}}_{\sin})^2}{(\hat{h}_{\cos})^2 + (\hat{h}_{\sin})^2}};\]

\[\sigma_{\Delta\Phi} = \frac{\sqrt{(\hat{h}_{\cos})^2(\sigma^{\hat{h}}_{\cos})^2 + (\hat{h}_{\sin})^2(\sigma^{\hat{h}}_{\sin})^2}}{(\hat{h}_{\cos})^2 + (\hat{h}_{\sin})^2};\]

The SNR is computed from:

\[{\rm SNR} = \frac{h_0}{\sigma_{h_0}}.\]

Example graphs of recovered injection properties vs.time are presented in section 4.1, along with summary graphs for all injections.

3.2 \(\mathcal{F}\)-statistic Method↩︎

The \(\mathcal{F}\)-statistic-based monitoring exploits well established LIGO Analysis Library (lalsuite) [32] software to compute both daily and cumulative \(\mathcal{F}\)-statistic [23] values for each day of the run for each interferometer. In addition to computing an \(\mathcal{F}\)-statisticvalue (denoted by \(2\!\mathcal{F}\)) for each injection, interferometer and time span, the program (lalpulsar_ComputeFstatistic_v2 [38]) computes best estimates and uncertainties for the injection amplitude \(h_0\), cosine of the inclination angle \(\iota\), phase offset \(\phi_0\) and polarization angle \(\psi\).

In brief, the \(\mathcal{F}\)-statisticis constructed from a combination of functions computed numerically from the SFTs of the observation span. These functions take into account the antenna pattern responses of the detector at the midpoints of the SFTs, and their combination is maximized with respect to \(h_0\), \(\cos(\iota)\), \(\phi_0\) and \(\psi\). In the \(2\!\mathcal{F}\)implementation, the source spindown is fixed at the intended injected value, but the frequency is allowed to vary over a narrow band around the intended value. The original \(\mathcal{F}\)-statisticdefinition may be found in [23], and the algorithm implementation in the lalsuite package may be found in [38].

Analytically maximizing \(2\!\mathcal{F}\)with respect to the four intrinsic source parameters allows solving for the best estimates for those parameters, with the Fisher matrix used in the maximization process implying estimated uncertainties. These quantitities are used in daily produced graphs.

Specifically, for each injection monitored and for each interferometer, graphs are produced for 1) the \(2\!\mathcal{F}\)for that day’s data; 2) the cumulative \(2\!\mathcal{F}\)from the start of the run through that day; and 3-6) cumulative graphs of the four estimated signal parameters vs.time. Examples are shown in Section 4.2.

3.3 Bayesian Method↩︎

The Bayesian method [9], [39] starts by greatly compressing the calibrated gravitational-wave strain data and then uses this compressed data, and the assumption of stationary noise with varying timescales, to estimate posterior probabilities on the unknown signal parameters: \(h_0\), \(\iota\), \(\psi\), and \(\Phi_0\). The method is able to coherently analyze any length of observed data provided the phase evolution can be accurately tracked over the observation period, which is the case for many known pulsars and certainly the case for the hardware injections. For the compression, for each hardware injection the phase evolution from Eqn. 2 is calculated over the observation period. The strain data is multiplied by \(e^{-i\Phi(t)}\) (with the value of \(\Phi_0\) not included), i.e., heterodyned, which produces a complex time series in which the signal power is equally split between a component now shifted to \(\sim 0\) Hz (D.C.; with small modulations left in from the diurnal antenna patterns), and one shifted to twice the signal frequency. The noise at the signal frequency is also equally split in the same way, meaning the signal-to-noise ratio for each component is the same as the original signal. This complex time series is then heavily low-pass filtered using a 9\(^{\rm th}\) order Butterworth IIR filter with a knee-frequency of 0.25 Hz, and down-sampled (via averaging) from a sample rate of 16 384 Hz to 1/60 Hz. The resulting heterodyned complex time series contains approximately white Gaussian noise over a bandwidth of \(\pm 1/120\) Hz about the signal, which now has the form

\[\begin{align} h'(t) = \frac{h_0}{2}\left(\frac{1}{2}F_{+}(t, \psi)(1+\cos{}^2\iota)e^{i\Phi_{0}} -iF_{\times}(t, \psi)\cos{\iota}e^{i\Phi_{0}} \right). \end{align}\]

where \(\Phi_{0}\) is defined in Eqn. 2 . To infer the signal parameters for each injection, their joint posterior probability distribution is calculated via Bayes’ rule

\[p(h_0, \cos{\iota}, \Phi_0, \psi|\mathbf{D}) \propto p(\mathbf{D}| h_0, \cos{\iota}, \Phi_0, \psi) p(h_0, \cos{\iota}, \Phi_0, \psi),\]

where \(\mathbf{D}\) is the heterodyned time series from a given detector and for a particular injection, \(p(\mathbf{D}| h_0, \cos{\iota}, \Phi_0, \psi)\) is the likelihood function of the data, and \(p(h_0, \cos{\iota}, \Phi_0, \psi)\) is the joint prior probability for the parameters. In this analysis, the priors are treated as a constant, i.e., uniform distributions over the allowed parameter ranges. The likelihood function assumes Gaussian noise with an unknown variance over a series of stationary chunks. The length of each stationary chunk is calculated via change-point analysis using a method based on BayesianBlocks [39], [40]. The unknown variance in each chunk can be marginalized over, to produce a Student’s t-likelihood [9]. The full likelihood is then calculated as the product of the likelihoods for each chunk.

In practice, the posterior probability distribution is sampled from using the Nested Sampling method [39], [41], [42], giving a series of samples \(\{h_0, \cos{\iota}, \psi, \Phi_0\}\) that are distributed according to the posterior and from which credible intervals from the marginal distributions for each parameter can be calculated.

4 Results↩︎

Here we present example and summary results from the three complementary monitoring programs. For simplicity, most results are shown for strain data not subjected to self-gating of loud transients, but similar and consistent results are observed for corresponding self-gated data (see below).

4.1 Template Generation Method↩︎

4.1.1 Differential Analysis of Weak and Strong Injections↩︎

The analysis focuses on recovering and comparing key parameters, including amplitude, phase, and frequency evolution, for both weak and strong injections. This approach involves assessing the consistency of recovered strain amplitudes across signals of different strengths, evaluating the precision of phase recovery to identify biases between weak and strong signals, and examining variations in the SNR as a function of signal strength to gain insights into the detector’s performance and reliability [33] .

Diagnostic plots shown below are used to visualize trends and discrepancies in these metrics. For weak signals, the proximity to the noise floor may result in greater variability in recovered parameters, whereas strong signals exhibit higher recovery accuracy. These visualizations help identify performance boundaries and calibration inconsistencies.

Sample results from the cumulative SNR monitoring of Injection 6 (\(\sim\)​144 Hz) and Injection 8 (\(\sim\)​189 Hz) injections in both H1 and L1 are shown in Figure 1. In each panel the deep pink curve is the expected \(\sqrt{T}\) growth of the injected SNR for the intended injection amplitude, the dark blue curve is the daily recovered SNR from our narrow-band search, and the gray bars show the \(\pm1\sigma\) uncertainty across the frequency band.

One sees in H1 for Injection 6 (upper left panel) excellent agreement between expected and actual rise in SNR, but L1 (lower left) there is some discrepancy by the end of the O4a period. In these cumulative graphs, adjacent SNR values and their error bars shown are highly correlated since nearly the same data is used to evaluate them. For the weaker Injection 8 (upper and lower right), both detectors exhibit larger early fluctuations as the estimator initially latches onto noise excursions, but by the end of O4a, there is good agreement between expectation and signal recovery.

Figure 2 presents the SFT‐weighted estimates of the gravitational‐wave amplitude \(h_0\) for Injection 6 (\(\sim\)​144 Hz) and Injection 8 (\(\sim\)​189 Hz) in the H1 (top row) and L1 (bottom row) detectors. In each plot the deep pink curve marks the true injected \(h_0\), the dark blue trace the recovered weighted average, and the gray error bars the \(\pm1\sigma\) spread across individual SFT estimates. For the strong Injection 6, both H1 and L1 show an initial noise‐driven overshoot in early June, but by mid‐June the recovered \(h_0\) converges to within a few percent of the intended injected value ( in both H1 and L1), and the uncertainties shrink as more SFTs accumulate. A small discrepancy is seen in the reconstructed L1 amplitude for Injection 6, but consistent with a statistical fluctuation from the collective mean reconstructed / intended amplitude (see Fig. 6). By contrast, the weaker Injection 8 exhibits large early scatter—error bars exceeding 30% of the injected amplitude — and a pronounced bias as the maximization briefly latches onto statistical fluctuations. Over the next three months, however, both detectors’ weighted averages drift steadily toward the true \(h_0\) and their \(1\sigma\) bands contract, achieving agreement within the error envelope by September 2023 in H1 and by October in L1. The fluctuations early in the run for weak injections merely

Figure 3 illustrates the cumulative phase offset reconstruction for Injection 6 (\(\sim\)​144 Hz) and Injection 8 (\(\sim\)​189 Hz) in the H1 (top row) and L1 (bottom row) detectors. In each column the upper panel shows the SFT‐weighted average of the cosine component of \(h_0\), the middle panel the sine component, and the lower panel the recovered signal phase, all compared against the injected references (deep pink curves) with their \(\pm1\sigma\) spreads (gray bars). For the strong Injection 6, both H1 and L1 exhibit an early noise‐driven overshoot in the cosine term that settles to the true \(h_0\) within a few weeks, while the sine component decays to near zero and the recovered phase (dark blue) rapidly locks onto the injected phase with shrinking uncertainties.

In contrast, the weaker Injection 8 shows larger initial scatter in both cosine and sine components – particularly in L1 – accompanied by broader phase error bars; over the course of O4a the weighted averages converge toward the intended injected values, and the recovered phase offset agrees with the intended value.

Figure 4 presents the (magnified) evolution of the recovered signal relative phase offset (magenta trace with gray \(\pm\)​1-\(\sigma\) error bars) compared to the intended value (zero, by definition – red curve) for Injection 6 (\(\sim\)​144 Hz) and Injection 8 (\(\sim\)​189 Hz) in H1 (top row) and L1 (bottom row). In H1 for Injection 6 (upper left), the relative phase offset estimate exhibits a large initial scatter—spanning \(\pm\)​150\(^\circ\)—as the maximization latches onto noise fluctuations, but by mid-June it rapidly converges to within a few degrees of the injected phase and remains stable thereafter. L1 for Injection 6 (lower left) shows a similar pattern with a larger early bias and broad uncertainties, converging toward the intended value by late June.

For the weaker Injection 8 in H1 (upper right), the recovered relative phase offset initially undershoots (to –100°) then overshoots (+20\(^\circ\)) before drifting toward zero over the next three months; in L1 (lower right), the phase begins with extreme fluctuations (\(\pm\)​200\(^\circ\)) and a pronounced negative bias, but gradually aligns to the injected phase by July and holds within \(\pm\)​5\(^\circ\) for the remainder of O4a.

4.1.2 Integrated Summary of Detector Response to Injections↩︎

The daily summary graphs provide a more global snapshot of H1 and L1 calibration fidelity over the full search band for CW searches, giving insight into potential broadband discrepancies not easily seen for an individual injection.

Figure 5 compares the injected and recovered gravitational-wave amplitudes \(h_0\) for all 18 hardware injections, plotted as a function of signal frequency for the H1 (left) and L1 (right) detectors. In each panel the deep pink “×” marks the true injected \(h_0\), the dark blue square the SFT-weighted average recovered by our narrow-band search, and the adjacent label the pulsar index. The vertical axis spans nearly three decades of amplitude—from \(\sim10^{-26}\) to \(10^{-23}\), reflecting the range of both intended SNRs and frequency-dependent detector noise levels. Very good agreement between intended and recovered amplitudes is visible in these graphs, but a more discriminating global comparison comes from measuring the fractional differences, namely the differences between recovered and intended amplitudes divided by the intended values.

These relative amplitudes are shown in Figure 6, with H1 values and uncertainties in the left panel and L1 values in the right panel. For H1, one sees excellent agreement with the ideal value of zero with a weighted mean of -0.0083 \(\pm\) 0.0025. For L1, one sees a slight discrepancy of -0.0112 \(\pm\) 0.0023, with a corresponding bias toward negative values visible in the graph. This discrepancy is consistent with the estimated uncertainty for the actuation functions determined shortly before the start of O4a running.

Similarly, combined phase offset reconstruction graphs provide a validation of the understanding of timing delays, both in the detector response and in the data stream. Figure 7 shows the reconstructed phases (radians) for all injections in H1 (left panel) and L1 (right panel). One sees relative phase values near the ideal of zero. If there were an uncorrected fixed time delay (or advance) in the detector response, one would expect a phase difference that grows linearly with frequency. To check for this potential discrepancy, we convert reconstructed relative phase into an equivalent time offset, as shown (\(\mu\)s) in Figure 8. One sees that the weighted mean H1 time delay is about 3 \(\mu\)s, and that for L1 it is about 6 \(\mu\)s. Both values lie within the actuation function timing uncertainty estimated prior to the start of O4a. For reference, the clock cycle for the controls and data acquisition system is about 61 \(\mu\)s.

LIGO data collected during both the O3 and O4 oberving runs have been characterized by infrequent but extremely loud transient instrumental glitches of uncertain origin, glitches loud enough in their strain power to raise the noise floor below \(\sim\)​500 Hz appreciably even for spectra meausered over time scales of 30 minutes or more. Identification and self-gating of these transients has been carried out for the O3 [34] and O4 [35] data. The CW hardware injection provide a handy validation of this gating in establishing that SNR for gated data improves relative to that of ungated data at low frequencies and that no statistically significant bias in reconstructed parameters is created by the gating. Figure 9 shows, for example, reconstructed \(h_0\) values for H1 and L1 for all 18 injections. As expected, central values change little, but error bars are somewhat reduced for gated data.

Collectively, these amplitude and phase offset reconstructions validate the fidelity of our understanding of the detector response to continuous gravitational waves over as many as tens of billions of signal cycles in the O4a run period.

4.2 F-statistic Method↩︎

Sample results from the \(\mathcal{F}\)-statisticmonitoring of injections are presented here. Figure 10 shows graphs of cumulative H1 and L1 \(2\!\mathcal{F}\)values for Injection 6 (\(\sim\)​144 Hz), one of the loud signals, through the O4a period, along with estimated \(h_0\) values and uncertainties. One sees a monotonically increasing cumulative \(2\!\mathcal{F}\)value, as expected. For stationary detection noise and stable data collecting livetime, one would expect an approximately linear rise in \(2\!\mathcal{F}\)with time. Downward deviations in slope can arise from poor livetime and/or degraded noise during a given day or days. The estimated \(h_0\) is biased upward early in the run because the maximization over a narrow band can initially sieze upon an upward statistical fluctuation. Ideally, as the run progresses, the program converges upon the correct frequency bin, as the estimated amplitude drifts toward the intended value.

Figure 10 also shows the cumulative estimated \(\cos(\iota)\), \(\phi_0\) and \(\psi\) for the same Injection 6. Again, one sees deviations from the intended values that decrease with time, as expected. For this nearly linearly polarized signal (\(|\cos(\iota)| \ll 1\)), the trends behave as one would naively expect. For nearly circularly polarized signal (\(|\cos(\iota)|\approx 1\)), on the other hand, there is a strong anti-correlation between reconstructed \(h_0\) and \(|\cos(\iota)|\) values, and near degeneracy in the definitions of the phase offset \(\phi_0\) and the polarization angle \(\psi\). Figure 11 shows corresponding results for an example (Injection 0) with \(|\cos(\iota)|\) = 0.795. One sees early in the run a strong anti-correlation between deviations of the reconstructed \(h_0\) and \(\cos(\iota)\) values from their injected values. Similarly, one sees a strong anti-correlation between deviations of the reconstructed \(\phi_0\) and \(\psi\) values from their injected values.

A rough measure of the detected strength of a signal can be obtained via an effective strain defined by the quadrature sums of the magnitudes of the plus- and cross-polarized signals defined in Eqn. 1 : \(h_{\rm eff} \equiv \sqrt{h_+^2+h_\times^2}\). From Eqn. 1 , one sees that \(h_{\rm eff}\) increases monotonically with \(h_0|\cos(\iota)|\), hence the anti-correlation seen in Fig. 11.

Figure 12 shows the \(2\!\mathcal{F}\)values recovered for all of the isolated-star injections for both H1 and L1 by the end of O4a. The values obtained for the two detectors differ because of different noise levels and integrated livetimes. For most injections, the recovered L1 values are higher. The loudest six injections used for daily monitoring of unintended calibration changes are apparent, especially for L1. Figure 12 also shows the \(h_0\) values recovered for the same injections in both detectors, along with the intended values. One sees good agreement for most injections (albeit with large uncertainties for Injection 13 12.34 Hz which is in a band marked by highly non-Gaussian artifacts. The region below 20 Hz has not traditionally been used for CW searches.

Notable discrepanices, however, are seen for Injections 2, 7 and 10. These discrepancies in reconstructed \(h_0\) are reflected in Fig. 13 which shows corresponding reconstructed \(\cos(\iota)\) values. Combining \(h_0\) and \(\cos(\iota)\) into effective strain \(h_{\rm eff}\) yields much better agreement, as expected, which is shown in Fig. 13, where excellent consistency is seen for every injection. Nonetheless, the disentanglement of the four source parameters reconstructed in single-detector \(\mathcal{F}\)-statisticanalysis highlights the value of the templated monitoring method for clean validation of amplitude and phase response.

4.3 Bayesian Method↩︎

The hardware injections for both H1 and L1 were recovered using data between 2023-05-24 15:00:00 and 2024-01-16 16:00:00. For each injection and detector, the data was heterodyned as described in Section 3.3 and the posterior probability distributions for \(\{h_0, \cos{\iota}, \psi, \Phi_0\}\) were estimated. Of particular interest are the recovered amplitude and initial phase. In Figure 14, the recovered posteriors on \(h_0\) and \(\cos{\iota}\) are recast as posteriors on \(h_+\) and \(h_{\times}\), as in Eqn. 1 , and the distributions are shown as violin plots of the fractional offset from the expected values and the true absolute values. The narrower the distribution, the more precisely the amplitude could be recovered. For the fractional offset, the closer to zero means that the recovered injection amplitude more closely matched the expected value. We see no specific trend in amplitude offsets as a function of frequency and the posteriors are broadly consistent with the expected amplitudes.

In Figure 15, we similarly see violin plots of posterior distributions representing the difference between the recovered \(\Phi_0\) distribution and the expected injection value. There is no trend in the offset as a function of injection signal frequency and the posteriors are broadly consistent with the injected value. The broadest posteriors generally reflect the high degree of correlation between the \(\Phi_0\) and \(\psi\) which occurs as the signal becomes closer to circularly polarized, i.e., as \(h_{\times}\) comes to dominate the amplitude.

5 Conclusion↩︎

In summary, this paper presents three techniques used in the real-time monitoring and analysis of hardware injections in LIGO data, highlighting their role in validation of searches for astrophysical sources of continuous gravitational waves.

The first, templated technique applies all of the expected frequency and amplitude modulations, in order to achieve the best precision on absolute amplitude and phase offset. The second technique, based on the long-established \(\mathcal{F}\)-statistic, allows for unknown stellar orientation angles and hence unknown amplitude modulations. The third, Bayesian inference technique also allows for unknown orientations and presents posterior distribution functions. The complementary approaches and the mostly independent code bases give further confidence in the validation.

As we move closer to the detector sensitivity necessary for the first-ever detection of an astrophysical CW source, it is critical to verify our search and signal reconstruction techniques, including understanding detector response well enough to phase-track a source over billions of signal cycles.

6 Acknowledgments↩︎

We gratefully acknowledge useful discussions and long collaboration with current and former colleagues in the LIGO-Virgo-KAGRA continuous waves working group and the LIGO calibration working group. We thank, in particular, David Barker, Joseph Betzweiser, Evan Goetz, Michael Thomas and Keith Thorne for ensuring reliable, nearly 24/7 LIGO hardware injections during the O4 run. We also thank Evan Goetz for helpful suggestions and comments on this manuscript. This work was supported in part by National Science Foundation Awards PHY-2110181 and PHY-2408883.

This material is based upon work supported by NSF’s LIGO Laboratory which is a major facility fully funded by the National Science Foundation. LIGO was constructed and is operated by the California Institute of Technology and Massachusetts Institute of Technology with funding from the U.S. National Science Foundation under grant PHY- 0757058. The authors also gratefully acknowledge the support of the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO. The authors are grateful for computational resources provided by the LIGO Laboratory and supported by National Science Foundation Grants PHY-0757058 and PHY-0823459.

References↩︎

---

1. The Virgo and KAGRA Collaborations rely on unmodulated actuations for similar validation↩︎