The NANOGrav 12.5 yr Data Set: Search for Gravitational Wave Memory

Any system that radiates gravitational waves (GWs) will also cause a permanent change in the spacetime metric. This effect, which was first derived in Zel'dovich & Polnarev (1974), is called gravitational wave memory. Later on, it was discovered that the gravitational waves emitted by a system are themselves a source of memory (Christodoulou 1991; Wiseman & Will 1991; Blanchet & Damour 1992; Thorne 1992). This effect is known as nonlinear memory, because it arises from nonlinearities in the Einstein field equations. Much work has been done to estimate the size of the effects of nonlinear GW memory, and it has been shown that there is a reasonable chance that the GW memory effect is significant enough to be observed using modern GW detectors (Wiseman & Will 1991; Arun et al. 2004; Favata 2009a, 2009b, 2010).

One such GW detector is a pulsar timing array (PTA). A PTA is a collection of millisecond pulsars (MSPs) that have extremely stable rotational periods (Lorimer 2008). Because of their stability, it is expected that, by carefully observing the times of arrival (TOAs) of radio pulses from these MSPs, it is possible to observe timing residuals induced by the passage of GWs (Sazhin 1978; Detweiler 1979; Hellings & Downs 1983; Agazie et al. 2023). The combination of multiple MSPs into a PTA also offers boosted sensitivity when trying to detect GW signals that produce predictable correlations among multiple pulsars (Foster & Backer 1990; Lommen 2015). Currently, there are several PTA collaborations in operation, including the North American Nanohertz Observatory for Gravitational Waves (NANOGrav; McLaughlin 2013), the European Pulsar Timing Array (EPTA; Desvignes et al. 2016), the Parkes Pulsar Timing Array (PPTA; Manchester et al. 2013), and the Indian Pulsar Timing Array (InPTA; Paul et al. 2019). Together, these collaborations form the International Pulsar Timing Array (IPTA; Verbiest et al. 2016). In addition, the Chinese Pulsar Timing Array (CPTA; Xu et al. 2023) and the MeerTime Pulsar Timing Array (MPTA; Miles et al. 2022b) have recently released their first analyses.

In the absence of any exotic physics, PTAs are expected to first detect a gravitational wave background originating from an ensemble of supermassive black hole binary (SMBHB) systems, followed by continuous waves from particularly bright SMBHBs (Rosado et al. 2015). Much work has already been done to characterize the GW background and place limits on continuous GWs (e.g., Arzoumanian et al. 2016, 2018a, 2020; Chen et al. 2021; Antoniadis et al. 2022; Agazie et al. 2023; Arzoumanian et al. 2023; Falxa et al. 2023). During the final inspiral and merger of a SMBHB system, the SMBHB strongly emits GWs that are outside the frequency band detectable by PTAs. However, the accumulated memory from these GWs may be significant enough to be detected by PTAs (Seto 2009; Pshirkov et al. 2010; van Haasteren & Levin 2010; Cordes & Jenet 2012; Madison et al. 2014). In addition to these signals, it is also possible to detect or constrain more exotic GW sources using PTA data. For example, cosmic strings are expected to emit strong bursts of GWs (Damour & Vilenkin 2000; Siemens et al. 2006). Yonemaru et al. (2020) have placed limits on GW bursts from cosmic strings based on the second PPTA data release.

Several studies have already been done using PTA data to constrain GW memory. NANOGrav has published constraints on GW memory using their 5 yr and 11 yr data sets (Arzoumanian et al. 2015; Aggarwal et al. 2020, hereafter NG5mem and NG11mem, respectively). The PPTA has also published constraints in Wang et al. (2015). Madison et al. (2016) have published the results of a search for GW memory from five galaxy clusters using PPTA data.

Additionally, several studies have considered GW memory using ground- and space-based detectors like LIGO–Virgo–KAGRA and LISA. Lasky et al. (2016) suggest a method to detect accumulated memory of many individual mergers, each of which is too weak to see by itself. Hübner et al. (2021) find no evidence of memory using data from 50 detections of GWs made during the third observation run of LIGO and Virgo. Boersma et al. (2020) forecasted that GW memory could be detected with total S/N = 3 after approximately five years of aLIGO operation at design sensitivity. Favata (2009a) estimated that memory from SMBHB mergers may be detectable (S/N ∼ 5) out to z ≲ 2 with LISA, and Islo et al. (2019) estimated that LISA may see between 1 and 10 memory events in its lifetime.

In this paper, we present our analysis of the NANOGrav 12.5 yr data set (Alam et al. 2020) for GW memory. We find that there is no significant evidence for GW memory in the data set. The model including a memory signal and a common spatially uncorrelated red noise (CURN) process is only very marginally favored, with a Bayes factor of 2.8, when compared to a model including only a CURN process. The posteriors from a full PTA analysis show that there is a very weak hint for GW memory at three different epochs: MJDs 54000, 55400, and 57300. However, a more detailed analysis shows that these three features are spurious. Each event is only supported by one or two pulsars, and one even lies inside a data gap in which there are no TOAs. When using more flexible red-noise models for three of the pulsars in the data set, the Bayes factor drops to BF = 1.3.

Thus, finding no GW memory events, we present upper limits constraining the amplitudes of any GW memory as functions of trial burst epoch and sky location. In addition, we use the constraints as a function of burst epoch to set constraints on rates of all astrophysical events which produce GW memory.

In Section 2, we will describe the NANOGrav 12.5 yr data set. Then, in Section 3, we will discuss how a GW memory wave front affects TOAs from a PTA. Next, in Section 4, we will summarize the mathematical techniques and software used in this search. Finally, we will discuss the results in Section 5.

In this paper, we analyze the NANOGrav 12.5 yr narrowband data set (Alam et al. 2020). We will briefly summarize some key points about this data set, but more details may be found in Alam et al. (2020) (hereafter NG12).

This data set contains TOAs from observations of 47 pulsars made between 2004 July and 2017 June. However, in this analysis, we used only the 45 pulsars with at least three years of data. These observations were performed using the Arecibo Observatory (AO) and Green Bank Telescope (GBT). All pulsars in the decl. range 0° < δ < + 39° were observed at Arecibo, with the remaining pulsars observed at GBT. In addition, PSRs B1937+21 and J1713+0747, which lie within the aforementioned decl. range, were also observed at GBT. Based on work by Burt et al. (2011) and Christy et al. (2014), six pulsars were also observed in a high-cadence program: PSRs J0030+0451, J1640+2224, J1713+0747, J1909−3744, J2043+1711, and J2317+1439. These six pulsars were observed weekly starting from 2013 at GBT and 2015 at AO. The remaining pulsars were observed monthly.

Each pulsar was, where possible, observed using two different receivers at different frequency ranges, to help understand and model out interstellar medium (ISM) and dispersion measure (DM) effects. At Arecibo, observations were performed with the 1.4 GHz receiver, and one of either the 430 MHz or 2.1 GHz receivers, depending on the noise characteristics of the observed pulsar. At the GBT, monthly observations used the 1.4 GHz and 820 MHz receivers. Weekly observations, however, only used the 1.4 GHz receiver. The observations were initially recorded using the ASP/GASP backends at AO and GBT, respectively (Demorest 2007). Later, between 2010 and 2012, these backends were replaced by the wideband backends PUPPI/GUPPI at AO and GBT, respectively (DuPlain et al. 2008; Ford et al. 2010).

Each pulsar's timing model was fitted using TEMPO ⁶⁴ (Nice et al. 2015) and checked for consistency using TEMPO2 ⁶⁵ (Hobbs et al. 2006) and PINT ⁶⁶ (Luo et al. 2019).

In this section, we will discuss the effects of GW memory on TOA residuals in pulsar timing data and summarize all of the components of the signal model and noise model used in this search.

Qualitatively, a memory wave front passing over a single pulsar will cause the observed rotational frequency of the pulsar to suddenly increase or decrease by a constant amount. A memory wave front passing over the Earth will cause the observed rotational frequencies of all pulsars in the PTA to suddenly increase or decrease by a constant amount. In either case, this sudden change to the observed rotational frequency will introduce timing residuals because of the difference between the pulsar's timing model fitted rotational frequency and the observed rotational frequency. Because the difference between the expected frequency and observed frequency is a constant, this will cause residuals to accumulate linearly over the course of the observation. In the case of a memory wave front passing over just one pulsar, we will only see residuals in the TOAs of that pulsar. If a wave front passes over the Earth, we will see the residuals begin accumulating in the TOAs of every single observed pulsar.

For a GW memory wave front propagating in the direction $\hat{{\boldsymbol{k}}}$ with polarization angle ψ passing over the line of sight to a pulsar at sky position p , the residuals induced may be calculated (Estabrook & Wahlquist 1975; Hellings & Downs 1983):

$$\begin{eqnarray}&&\delta {t}_{\mathrm{mem}}(t)=B(\hat{{\boldsymbol{k}}},\hat{{\boldsymbol{p}}},\psi ){h}_{\mathrm{mem}}(t).\end{eqnarray} \tag{ 1 }$$

The projection factor $B(\hat{{\boldsymbol{k}}},\hat{{\boldsymbol{p}}},\psi )$ accounts for the fact that the GW memory has a quadrupolar antenna pattern. The effect of the GW memory on the TOAs from a pulsar depends on where that pulsar lies inside the antenna pattern. For a pulsar at p and a wave front propagating in the direction $\hat{{\boldsymbol{k}}}$ separated by an angle α, we can write the projection factor as

$$\begin{eqnarray}&&B(\hat{{\boldsymbol{k}}},\hat{{\boldsymbol{p}}},\psi )=\displaystyle \frac{1}{2}\cos (2{\psi }_{\hat{{\boldsymbol{p}}}})(1-\cos \alpha ),\end{eqnarray} \tag{ 2 }$$

where α is the angle between $\hat{{\boldsymbol{p}}}$ and $\hat{{\boldsymbol{k}}}$, and ${\psi }_{\hat{{\boldsymbol{p}}}}$ is the angle between the principal polarization vector (defined by ψ) and the pulsar line of sight projected onto a plane perpendicular to $\hat{{\boldsymbol{k}}}$.

The second factor in Equation (1) carries the strength and time dependence of the burst. For a wave front with a characteristic strain of h₀, we may write the time dependence as

$$\begin{eqnarray}&&{h}_{\mathrm{mem}}(t)={h}_{0}[(t-{t}_{0}){\rm{\Theta }}(t-{t}_{0})-(t-{t}_{i}){\rm{\Theta }}(t-{t}_{i})],\end{eqnarray} \tag{ 3 }$$

where t₀ is the time at which the memory wave front passes over the Earth, ${t}_{i}={t}_{0}+(| {{\boldsymbol{p}}}_{i}| /c)\hspace{1mm}[1+\cos ({\theta }_{i})]$ is the retarded time at which the same wave front passed over the pulsar, and Θ is the Heaviside function. The left- and right-hand terms in Equation (3) are called the "Earth term" and "pulsar term," respectively. In reality, because the distance to each pulsar in our PTA is on the order of thousands of light-years, and the observation baselines of ongoing PTA experiments is tens of years, we expect only one nonzero term in Equation (3).

The characteristic strain of GW memory h₀ depends on the amount of energy radiated as GWs, ΔE_rad, and the orientation and distance of the source relative to the observer. For a binary merger

$$\begin{eqnarray}&&{h}_{0}=\displaystyle \frac{{\rm{\Delta }}{E}_{\mathrm{rad}}}{24r}{\sin }^{2}\iota \left(17+{\cos }^{2}\iota \right),\end{eqnarray} \tag{ 4 }$$

where r is the comoving distance to the source, ι is the orbital inclination angle of the binary, and ΔE_rad is a function of the individual masses and spins. As our signal model includes only the memory portion, we are agnostic to the particulars of the signal's origin and parameterize our model with h₀ directly.

Because of the significance of the detection of a CURN process in NG12gwb, we also include a CURN with a fixed spectral index as part of our model. The CURN is modeled as a power law, with a power spectrum characterized by two hyperparameters (A, γ) (Phinney 2001):

$$\begin{eqnarray}&&P(f)={A}^{2}{\left(\displaystyle \frac{f}{{\mathrm{yr}}^{-1}}\right)}^{-\gamma },\end{eqnarray} \tag{ 5 }$$

where f is the frequency of the spectral component, A is the characteristic amplitude of the red-noise process at the reference frequency yr⁻¹, and γ is the spectral index of the process. A stochastic gravitational wave background generated by an ensemble of SMBHB is expected to have a spectral index of γ_SMBHB = 13/3 ≈ 4.33. However, the maximum a posteriori value found for the spectral index of the CURN in NG12gwb was approximately γ_MAP = 5.5. In this paper, we present two sets of results using both of these fixed spectral indices for the CURN.

In addition to these signals, we include Gaussian white noise and Gaussian red noise on a per-pulsar basis. The Gaussian red noise accounts for long-timescale changes in the pulsar's rotational frequency. Some examples of processes that can cause these changes include spin noise (Shannon & Cordes 2010; Lam et al. 2016), stochastic variations in dispersion measure (Demorest et al. 2013; Keith et al. 2013; Jones et al. 2017), and mode changing (Lyne et al. 2010; Miles et al. 2022a).

The Gaussian white noise is parameterized by three parameters known as EQUAD, EFAC, and ECORR. EQUAD and EFAC modify the measured TOA uncertainty: EQUAD adds additional white noise in quadrature, and EFAC multiplies the total TOA uncertainty after including EQUAD. ECORR describes white noise that is correlated between TOAs gathered in the same observation epoch but uncorrelated between different observations. This term nominally accounts for pulse jitter noise. For this analysis, the white noise parameters are fixed to their median values as determined by single pulsar noise analyses, for the sake of computational efficiency.

The techniques used in this search are documented in Sun et al. (2023). As such, in this section, we will give only a brief overview of the techniques. The residuals in a single pulsar's TOAs may be written as the sum of multiple stochastic and deterministic processes:

$$\begin{eqnarray}&&{\boldsymbol{\delta }}{\boldsymbol{t}}={\boldsymbol{\delta }}{{\boldsymbol{t}}}_{\mathrm{mem}}+M{\boldsymbol{\epsilon }}+F{\boldsymbol{a}}+{F}_{\mathrm{gw}}{{\boldsymbol{a}}}_{\mathrm{gw}}+{\boldsymbol{n}}.\end{eqnarray} \tag{ 6 }$$

Above, δ t are the residual time series for the pulsar. The term δ t _mem represents the residuals induced by GW memory; M is the design matrix accounting for small errors in the linearized pulsar timing model ; F is the design matrix for a pulsar-intrinsic Gaussian red-noise process modeled as a Fourier series with coefficients a ; similarly, F_gw and a _gw are the design matrix and Fourier coefficients for the CURN; finally, n represents the uncertainties in the TOAs from Gaussian white noise.

Given estimations of the timing model parameters, GW memory signal, and Gaussian process parameters, we can construct residuals r :

$$\begin{eqnarray}&&{\boldsymbol{r}}={\boldsymbol{\delta }}{\boldsymbol{t}}-{\boldsymbol{\delta }}{{\boldsymbol{t}}}_{\mathrm{mem}}-M{\boldsymbol{\epsilon }}-F{\boldsymbol{a}}-{F}_{\mathrm{gw}}{{\boldsymbol{a}}}_{\mathrm{gw}}.\end{eqnarray} \tag{ 7 }$$

Because the residuals r are expected to arise only from Gaussian white noise (having subtracted out all other effects), we can compute the likelihood of any set of model parameters as

$$\begin{eqnarray}\begin{array}{rcl}L({\boldsymbol{\epsilon }},{\boldsymbol{a}},{{\boldsymbol{a}}}_{\mathrm{gw}},{\boldsymbol{\delta }}{{\boldsymbol{t}}}_{\mathrm{mem}}) & = & p({\boldsymbol{\delta }}{\boldsymbol{t}}| {\boldsymbol{\epsilon }},{\boldsymbol{a}},{{\boldsymbol{a}}}_{\mathrm{gw}},{h}_{\mathrm{mem}},{t}_{0},\hat{{\boldsymbol{k}}},\hat{{\boldsymbol{p}}},\psi )\\ & = & \displaystyle \frac{\exp \left(-\tfrac{1}{2}{{\boldsymbol{r}}}^{T}{N}^{-1}{\boldsymbol{r}}\right)}{\sqrt{2\pi \det N}}.\end{array}\end{eqnarray} \tag{ 8 }$$

It is also possible to analytically marginalize this likelihood over the timing model and red-noise parameters. For a full description, see Lentati et al. (2013), van Haasteren & Vallisneri (2014), and van Haasteren & Vallisneri (2015). The final marginalized likelihood is

$$\begin{eqnarray}&&p({\boldsymbol{\delta }}{\boldsymbol{t}}| {h}_{\mathrm{mem}},{t}_{0},\hat{{\boldsymbol{k}}},\hat{{\boldsymbol{p}}},\psi )=\displaystyle \frac{\exp \left(-\tfrac{1}{2}{{\boldsymbol{q}}}^{T}{C}^{-1}{\boldsymbol{q}}\right)}{\sqrt{2\pi \det C}},\end{eqnarray} \tag{ 9 }$$

where we have the definitions:

$$\begin{eqnarray}&&{\boldsymbol{q}}=\delta {\boldsymbol{t}}-\delta {{\boldsymbol{t}}}_{\mathrm{mem}},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&C=N+{\mathrm{TDT}}^{T},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}T=\left[\begin{array}{cc}M & F\end{array}\right],\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}D=\left[\begin{array}{cc}\infty & 0\\ 0 & \phi \end{array}\right],\end{eqnarray} \tag{ 13 }$$

where ∞ above is a diagonal matrix of infinities (which we can understand as unconstrained priors on timing model parameters), and ϕ is a covariance matrix for the individual red noise and CURN Fourier coefficients. Because we are using a CURN, these ϕ matrices are also diagonal, and they simply contain the red-noise power at each frequency bin given by Equation (5).

Following the procedure of Sun et al. (2023), which is based on methods in Section 6 of NG5mem, we then compute the pulsar-term likelihoods (marginalized over intrinsic pulsar red noise and fixed spectral index CURN) on a grid of trial parameters {h_i , t_B } where h_i is the post-projection strain of the memory signal in the ith pulsar, and t_B is the burst epoch. The post-projection strain is given by the product of the projection factor (Equation (2)) and the intrinsic strain of the memory signal h₀. We can only use the post-projection strain because the residuals of one pulsar cannot break the degeneracy between the location of the signal's origin and the intrinsic strain. Additionally, Sun et al. (2023) do not include a CURN, but we choose to include this additional noise process because of the results of NG12gwb, in which it was shown that there is significant evidence for a CURN in the NANOGrav 12.5 yr data set (NC12). These pulsar-term likelihood tables may be used to set upper limits on pulsar-term GW memory.

Then, we can combine the pulsar-term likelihoods to compute Earth-term likelihoods by making use of the factorizability of the signal model:

$$\begin{eqnarray}\begin{array}{rcl}p({\boldsymbol{\delta }}{\boldsymbol{t}}| \hat{{\boldsymbol{k}}},\psi ,{t}_{B},{h}_{0}) & = & \displaystyle \prod _{i=1}^{{N}_{\mathrm{psr}}}{p}_{i}({\boldsymbol{\delta }}{\boldsymbol{t}}| \hat{{\boldsymbol{k}}},\psi ,{t}_{B},{h}_{0})\\ & = & \displaystyle \prod _{i=1}^{{N}_{\mathrm{psr}}}{p}_{i}({\boldsymbol{\delta }}{\boldsymbol{t}}| {h}_{i},{t}_{B}),\end{array}\end{eqnarray} \tag{ 14 }$$

where above we have implicitly used Equation (2) to combine the burst parameters $\hat{{\boldsymbol{k}}}$, ψ, and h₀ into the post-projection, pulsar-term GW memory strain h_i . In this way, it becomes very computationally inexpensive to compute the red-noise-marginalized Earth-term likelihoods of any GW memory events on a full grid of trial parameters $\{{\mathrm{log}}_{10}{h}_{0},{t}_{B},{\theta }_{B},{\phi }_{B},{\psi }_{B}\}$, where h₀ is the intrinsic strain of the memory event, t_B is the event epoch, (θ_B , ϕ_B ) are the polar and azimuthal angles of the sky location of the event source, and ψ_B is the polarization angle of the memory wave front. Once we have the likelihoods on a grid of trial parameters, we can simply numerically marginalize over any of the trial parameters to obtain marginalized likelihoods or posterior probability distributions.

These signal models and this likelihood calculation are implemented in enterprise ⁶⁷ (Ellis et al. 2020) and enterprise_extensions ⁶⁸ (Taylor et al. 2021)

5.1. Earth-term Memory Search

We began by performing an Earth-term Bayesian search for GW memory using MCMC sampling. We compared two models: (1) a noise-only model and (2) a noise and GW memory model. The noise-only model included intrinsic pulsar red noise, white noise, and a common red-noise process. The signal model included the same noise processes with an additional GW memory signal. The two models were simultaneously sampled using the product-space sampling method (Carlin & Chib 1995; Godsill 2001), allowing us to determine the posterior probability for the memory signal and compute the Bayes factor for the signal model compared to noise only. The resulting Bayes factor of 2.8, shows the GW memory model is marginally favored over noise only. However, this Bayes factor is too small to be considered a detection. The posterior probability distributions for the memory signal and global spatially uncorrelated red-noise process are shown in Figure 1. Based on the posterior probability of the burst epoch, we can identify three "hot spots" near MJDs 54000, 55400, and 57300.

Zoom In Zoom Out Reset image size

The features near MJD 54000 and 55400 were both present in the analysis of NG11mem, where they were the most significant GW memory false-alarm events in NG11 and NG9, respectively. The feature near MJD 54000 lies near the start of our observations and at a time where there were large data gaps for several pulsars. At early times in our data set, there were fewer pulsars being observed, and the observations were less regular. The sparsity of data makes it harder to constrain any signal in these times. Events that occur early in the data are also more degenerate with the quadratic pulsar timing model fit to the pulsar rotational frequency and frequency derivative. This means that the signal model can be consistent with a high-amplitude memory event that is effectively removed by the marginalization of the timing model.

Using a dropout analysis, we identified three pulsars in particular that supported each of the three aforementioned features: PSRs J0030+0451, J1744−1134, and J2043+1711. We performed another Bayesian search using a free spectral noise model, which treats the power in each frequency bin in the power spectral density as an independent parameter, rather than requiring a power-law red-noise power spectral density for these three pulsars. We found that using a more flexible noise model in these pulsars completely removes the features at MJD 55400 and MJD 57300. Because each of these features has no support from any other pulsars, we conclude that these events are related to noise in individual pulsars and not actual GW memory events. The analysis using the free spectral noise model in three pulsars results in a Bayes factor for GW memory of 1.3. In general, more complex noise models like those used in J. Simon et al. (2023, in preparation) should help prevent noise features from contaminating searches for GW memory in the future.

5.2. Pulsar-term Upper Limits

Because we make no detection, we report upper limits on GW memory strain amplitude.

Figure 2 shows the pulsar-term upper limits on GW memory using both γ_MAP and γ_SMBHB. Because the pulsar-term upper limits are computed one pulsar at a time, we lose all information relating to the sky location of the signal. This amplitude upper limit is a constraint on the product of $B(\hat{{\boldsymbol{k}}},\hat{{\boldsymbol{p}}},\psi )$ and h₀, because these two terms are fully degenerate in the pulsar-term search. In other words, it is impossible to differentiate between a weak memory event or one that originated in the sky such that the antenna pattern is weak at the pulsar's sky location. We see that the choice of spectral index does not affect the pulsar-term upper limits very much in most cases. Some pulsars (e.g., PSRs B1937+21, J0613−0200, J0645+5158, and J1713+0747) show small but significant differences.

Zoom In Zoom Out Reset image size

5.3. Earth-term Upper Limits

Figures 3 and 4 show the upper limits on GW memory strain amplitude in the NANOGrav 12.5 yr data set as a function of burst epoch and sky location, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To compute the upper limits as a function of burst epoch (Figure 4), we must compute amplitude posteriors that have uniform priors over the sky and polarization. Thus, we started by splitting up the sky into 48 HEALPix ⁶⁹ (Gorski et al. 2005) sky pixels using healpy ⁷⁰ (nside=2) (Zonca et al. 2019). Then, for each sky pixel, we computed likelihood tables for global GW memory events using Equation (14) and the pulsar-term likelihood tables. Finally, for each trial burst epoch, we took an equal number of samples from the amplitude posteriors from each source-orientation bin at that trial epoch. We then concatenated the samples taken from each of these amplitude posteriors together to form a sky-averaged strain amplitude posterior.

We must sample each source-orientation bin independently to construct our sky-averaged posteriors, because of the nature of the memory signal. Our PTA does not have uniformly distributed pulsars, and as such, there are parts of the sky in which we have little to no sensitivity. If we are not careful about sampling, and instead search over the entire sky simultaneously, our amplitude would be dominated by samples taken from source orientations to which our PTA is completely insensitive. Furthermore, because there is much more prior volume at high amplitudes, these samples would all heavily bias our amplitude posteriors toward higher amplitudes, which our PTA has no way of ruling out. This sampling scheme, in which we concatenate samples from different source-orientation bins, guarantees that our posterior is marginalized uniformly over the prior (Malmquist 1922).

Figure 3 shows the upper limits on memory strain as a function of sky location and fixed common red-noise spectral index. To obtain these upper limits, we first computed Earth-term lookup tables for 768 HEALPix sky pixels (nside=8) marginalized over the polarization of the memory wave front (see the upper limits as a function of burst epoch). For these Earth-term lookup tables, we limit the prior on the burst epoch to the last three years because some pulsars do not have more than three years of data. We can see from this comparison that the upper limits differ slightly depending on the choice of fixed spectral index for the common red noise. While we do not believe these differences are significant, there is a difference pattern that is very similar to the antenna response of a GW memory event around PSR J0613−0200. After repeating the same analysis, but omitting this pulsar, the nearby differences largely disappear. This suggests that this pulsar contains a noise feature that is difficult to model accurately using only a red-noise power law and white noise. When mismodeled, the excess noise is conflated with a GW memory signal, thus causing the upper limit differences in Figure 3.

Figure 4 shows the upper limits on GW memory in the NANOGrav 12.5 yr data set plotted as solid curves. We also show the results of the NANOGrav 11 yr search for GW memory plotted as a dashed green curve. The upper limits computed from the NANOGrav 12.5 yr data set do not improve significantly upon those computed from NG11mem in the overlapping epochs. However, the increased volume of timing data and number of pulsars do clearly result in improvements on the upper limits of approximately half an order of magnitude when compared to early upper limits. Continued observation and growth of this PTA will cause the upper limits in the future to be even lower and thus give much more stringent limits on GW memory.

Figure 5 shows the upper limits on the rate of SMBHB mergers that produce GW memory computed using the results shown in Figure 4. We do this by counting the number of epochs that have lower strain upper limits than a given fixed strain. From this, we can then constrain the rate of events that have strains at or below this fixed strain. In addition, this figure shows the predicted rate in Islo et al. (2019). From the right-hand-side plot, we can see that our rate upper limits do not improve much when compared to the NANOGrav 11 yr results. We also include the sky-marginalized pulsar-term upper limits. Notably, in this analysis, the Earth-term rate upper limits are more constraining than the combined pulsar-term upper limits. This indicates that the PTA contains enough pulsars that the sensitivity at low strain amplitudes is no longer dominated by a few pulsars. In addition, the NANOGrav 12.5 yr pulsar-term rate upper limits are worse than the 11 yr rate upper limits. This is due to the additional red-noise model used in the analysis of the 12.5 yr data set.

Zoom In Zoom Out Reset image size

In this paper, we have shown that there is no significant detection of GW memory in the NANOGrav 12.5 yr data set. We have therefore set upper limits on the strain amplitude of any GW memory events in the NANOGrav 12.5 yr data set in the presence of the CURN detected in NG12gwb. The addition of a CURN to the noise model does not significantly affect the upper limits, but it does have some covariance with the GW memory signal. We also see from Figure 3 that PSR J0613−0200 gives significantly different strain upper limits for sources in its vicinity, depending on the choice of spectral index for a CURN process. Furthermore, these differences have a quadrupolar shape, similar to the antenna response of a GW memory signal. This indicates the presence of some excess low-frequency noise in this pulsar that can be conflated with GW memory. From Figures 4 and 5, we see that the additional data in the 12.5 yr data set continue to increase our sensitivity to GW memory, especially later in the data set. This in turn will allow us to continue placing more stringent limits on the rates of memory-producing events.

It is important to remember that the predicted event rate upper limits shown in Figure 5 are only for SMBHB mergers. While the prospects for detecting GW memory from SMBHB mergers are low, there are many exotic sources that may be expected to emit GWs and produce GW memory (Cutler et al. 2014). Furthermore, it has been recently discovered that there is some evidence of a GW background in the NANOGrav 15 yr data set (Agazie et al. 2023). The spectrum detected in Agazie et al. (2023) may be used to update estimations of SMBHB systems and recompute merger event rates. This work is critical for continued development of future searches for GW memory. Additionally, pulsar glitches, which are instantaneous changes in the rotational frequencies, produce a signal almost identical to that of a pulsar-term GW burst with memory (Cordes & Jenet 2012). Pulsar glitches have been observed in two millisecond pulsars, PSRs B1821−24 and J0613−0200, thus far (Cognard & Backer 2004; McKee et al. 2016). Of these two, J0613−0200 is in NANOGrav's timing data set. However, this glitch occurred before NANOGrav began timing the pulsar, and it should therefore have no effect on this pulsar's timing model or residuals. The pulsar-term upper limits presented in this analysis may be used to set upper limits on glitches in every other pulsar as well. In general, this analysis may be used to cross-validate any detection of any loud GW-producing event.

Finally, the search for GW memory can reveal interesting noise features of a PTA's constituent pulsars. For example, the analysis presented in NG11mem shows that PSRs J1909−3744 and J0030+0451 had some excess, unmodeled noise. This analysis shows that there is some excess noise in PSRs J1744−1134 and J2043+1711 that conspires to give support for a memory event at MJD 57300. This makes them good candidates for any future studies of noise characteristics, like those presented in Lam et al. (2016) and Hazboun et al. (2020). In addition, Figure 3 suggests that PSR J0613−0200 may have a noise transient that is highly covariant with red noise. Previous work has shown that scattering variations may result in excess correlated noise in pulsar timing data sets (Keith et al. 2013; Chalumeau et al. 2021; Goncharov et al. 2021). In particular, Main et al. (2020, 2023) have shown that data from PSR J0613−0200 show significant evidence of scattering variations. These scattering variations may be the source of the differences in GW memory upper limits in the vicinity of this pulsar when using different CURN spectral indices. As pulsar timing baselines become longer and PTA sensitivity to red-noise increases, it will be critically important to explore how strong red noise and components of each pulsar's timing models affect detection prospects of GW memory.

Overall, the search for GW memory remains a critical part of the GW analysis pipeline because of its use in cross-validation of any potential detections of loud GWs, its ability to reveal unmodeled noise, and the possibility that a GW memory event may reveal exotic GW sources. Continued methods development, as applied to both GW memory and intrinsic pulsar noise, will allow us to perform more robust searches for these sources using future data sets.

The NANOGrav collaboration receives support from National Science Foundation (NSF) Physics Frontiers Center award numbers 1430284 and 2020265, the Gordon and Betty Moore Foundation, NSF AccelNet award number 2114721, an NSERC Discovery Grant, and CIFAR. The Arecibo Observatory is a facility of the NSF operated under cooperative agreement (AST-1744119) by the University of Central Florida (UCF) in alliance with Universidad Ana G. Méndez (UAGM) and Yang Enterprises (YEI), Inc. The Green Bank Observatory is a facility of the NSF operated under cooperative agreement by Associated Universities, Inc. The National Radio Astronomy Observatory is a facility of the NSF operated under cooperative agreement by Associated Universities, Inc. L.B. acknowledges support from the National Science Foundation under award AST-1909933 and from the Research Corporation for Science Advancement under Cottrell Scholar Award No. 27553. P.R.B. is supported by the Science and Technology Facilities Council, grant No. ST/W000946/1. S.B. gratefully acknowledges the support of a Sloan Fellowship, and the support of NSF under award #1815664. The work of R.B., R.C., D.D., N.La., X.S., J.P.S., and J.T. is partly supported by the George and Hannah Bolinger Memorial Fund in the College of Science at Oregon State University. M.C. and S.R.T. acknowledge support from NSF AST-2007993. M.C. and N.S.P. were supported by the Vanderbilt Initiative in Data Intensive Astrophysics (VIDA) Fellowship. Support for this work was provided by the NSF through the Grote Reber Fellowship Program administered by Associated Universities, Inc./National Radio Astronomy Observatory. Support for H.T.C. is provided by NASA through the NASA Hubble Fellowship Program grant #HST-HF2-51453.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. M.E.D. acknowledges support from the Naval Research Laboratory by NASA under contract S-15633Y. T.D. and M.T.L. are supported by an NSF Astronomy and Astrophysics Grant (AAG) award number 2009468. E.C.F. is supported by NASA under award number 80GSFC21M0002. G.E.F., S.C.S., and S.J.V. are supported by NSF award PHY-2011772. The Flatiron Institute is supported by the Simons Foundation. A.D.J. and M.V. acknowledge support from the Caltech and Jet Propulsion Laboratory President's and Director's Research and Development Fund. A.D.J. acknowledges support from the Sloan Foundation. N.La. acknowledges the support from Larry W. Martin and Joyce B. O'Neill Endowed Fellowship in the College of Science at Oregon State University. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). D.R.L. and M.A.M. are supported by NSF #1458952. M.A.M. is supported by NSF #2009425. C.M.F.M. was supported in part by the National Science Foundation under Grants No. NSF PHY-1748958 and AST-2106552. A.Mi. is supported by the Deutsche Forschungsgemeinschaft under Germany's Excellence Strategy—EXC 2121 Quantum Universe—390833306. The Dunlap Institute is funded by an endowment established by the David Dunlap family and the University of Toronto. K.D.O. was supported in part by NSF grant No. 2207267. T.T.P. acknowledges support from the Extragalactic Astrophysics Research Group at Eötvös Loránd University, funded by the Eötvös Loránd Research Network (ELKH), which was used during the development of this research. S.M.R. and I.H.S. are CIFAR Fellows. Portions of this work performed at NRL were supported by ONR 6.1 basic research funding. J.D.R. also acknowledges support from startup funds from Texas Tech University. J.S. is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2202388, and acknowledges previous support by the NSF under award 1847938. Pulsar research at U.B.C. is supported by an NSERC Discovery Grant and by CIFAR. S.R.T. acknowledges support from an NSF CAREER award #2146016. C.U. acknowledges support from BGU (Kreitman fellowship), and the Council for Higher Education and Israel Academy of Sciences and Humanities (Excellence fellowship). C.A.W. acknowledges support from CIERA, the Adler Planetarium, and the Brinson Foundation through a CIERA-Adler postdoctoral fellowship. O.Y. is supported by the National Science Foundation Graduate Research Fellowship under grant No. DGE-2139292.

Author contributions. An alphabetical-order author list was used for this paper in recognition of the fact that a large, decade-timescale project such as NANOGrav is necessarily the result of the work of many people. All authors contributed to the activities of the NANOGrav collaboration leading to the work presented here, and reviewed the manuscript, text, and figures prior to the paper's submission. Additional specific contributions to this paper are as follows. Z.A., H.B., P.R.B., H.T.C., M.E.D., P.B.D., T.D., J.A.E., R.D.F., E.C.F., E.F., N.G.-D., P.A.G., D.C.G., M.L.J., M.T.L., D.R.L., R.S.L., J.L., M.A.M., C.N., D.J.N., T.T.P., N.S.P., S.M.R., K.S., I.H.S., R.S., J.K.S., and S.J.V. developed the NANOGrav 12.5 yr data set by performing observations, computing TOAs, checking data, and creating and refining timing models. P.T.B., A.D.J., D.R.M., X.S., and J.P.S. performed the various analyses. J.P.S. and P.T.B. coordinated paper writing.