The NANOGrav 15 yr Data Set: Search for Transverse Polarization Modes in the Gravitational-wave Background

Recently we found compelling evidence for a gravitational-wave background with Hellings and Downs (HD) correlations in our 15 yr data set. These correlations describe gravitational waves as predicted by general relativity, which has two transverse polarization modes. However, more general metric theories of gravity can have additional polarization modes, which produce different interpulsar correlations. In this work, we search the NANOGrav 15 yr data set for evidence of a gravitational-wave background with quadrupolar HD and scalar-transverse (ST) correlations. We find that HD correlations are the best fit to the data and no significant evidence in favor of ST correlations. While Bayes factors show strong evidence for a correlated signal, the data does not strongly prefer either correlation signature, with Bayes factors ∼2 when comparing HD to ST correlations, and ∼1 for HD plus ST correlations to HD correlations alone. However, when modeled alongside HD correlations, the amplitude and spectral index posteriors for ST correlations are uninformative, with the HD process accounting for the vast majority of the total signal. Using the optimal statistic, a frequentist technique that focuses on the pulsar-pair cross-correlations, we find median signal-to-noise ratios of 5.0 for HD and 4.6 for ST correlations when fit for separately, and median signal-to-noise ratios of 3.5 for HD and 3.0 for ST correlations when fit for simultaneously. While the signal-to-noise ratios for each of the correlations are comparable, the estimated amplitude and spectral index for HD are a significantly better fit to the total signal, in agreement with our Bayesian analysis.


Introduction
Einstein's theory of general relativity (GR) predicts the existence of gravitational waves (GWs) with two transverse polarization modes that propagate at the speed of light (Eardley et al. 1973a(Eardley et al. , 1973b)).Observations by the LIGO-Virgo-Kagra collaboration have shown that GR describes gravitational radiation from massive freely accelerating objects in the Universe (Abbott et al. 2016(Abbott et al. , 2023)).Although these observations have shown that GR describes observational data (Abbott et al. 2019a(Abbott et al. , 2019b(Abbott et al. , 2020(Abbott et al. , 2021a(Abbott et al. , 2021b;;The LIGO Scientific Collaboration et al. 2021), pulsar timing array (PTA) experiments offer a unique opportunity to probe other possible metric theories of gravity external to Einstein's GR (Lee et al. 2008;Yunes & Siemens 2013).
Modified theories of gravity are often introduced to resolve some of the current challenges facing fundamental physics, such as the nature of dark matter, and dark energy, and in attempts to reconcile quantum mechanics and gravity (see, e.g., Berti et al. 2015 and references therein).In metric theories of gravity, there can be up to six possible GW polarization modes (Eardley et al. 1973a(Eardley et al. , 1973b;;Will 1993).PTA searches for non-Einsteinian polarization modes may provide evidence for modified gravity theories by uncovering the different correlation patterns associated with such modes (Chamberlin & Siemens 2012;Yunes & Siemens 2013;Gair et al. 2015;Cornish et al. 2018;Bernardo & Ng 2023a, 2023b, 2023c;Afzal et al. 2023).
Millisecond pulsars (MSPs) emit radio beams from their magnetic poles and are extremely stable rotators.They appear to us as point sources of periodic radio bursts that arrive on 64 NASA Hubble Fellowship: Einstein Postdoctoral Fellow. 65NANOGrav Physics Frontiers Center Postdoctoral Fellow. 66Deceased. 67NSF Astronomy and Astrophysics Postdoctoral Fellow.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Earth with a consistency that rivals that of atomic clocks (Matsakis et al. 1997;Hobbs et al. 2012Hobbs et al. , 2019)).Pulsar timing experiments exploit the regularity of MSPs to search for lowfrequency (∼1-100 nHz) GWs by measuring deviations from the expected arrival time of radio pulses (Sazhin 1978;Detweiler 1979).Moreover, an array of these MSPs allows us to search for correlations between deviations of times of arrivals (TOAs) of pulses from MSP pairs (Hellings & Downs 1983;Foster & Backer 1990).
The North American Nanohertz Observatory for Gravitational Waves (NANOGrav), the European Pulsar Timing Array (EPTA), the Chinese Pulsar Timing Array, and the Parkes Pulsar Timing Array (PPTA) are the PTAs that possess the most sensitive data sets capable of measuring nHz GWs.NANOGrav, the EPTA, and the PPTA have seen strong evidence for a common red noise process (Arzoumanian et al. 2020;Chalumeau et al. 2021;Goncharov et al. 2021).Most recently in Agazie et al. (2023a; hereafter referred to as NG15), NANOGrav has found compelling evidence for quadrupolar correlations (Hellings & Downs 1983), while the EPTA and PPTA have seen these correlations at varying levels of significance (Antoniadis et al. 2023;Reardon et al. 2023).In this Letter, we complement our work in NG15 by searching for evidence for scalar-transverse (ST) correlations from the non-Einsteinian breathing polarization mode of gravity.Previous works have both set constraints and upper limits (Wu et al. 2022;Bernardo & Ng 2023d) as well as shown preference for ST- (Arzoumanian et al. 2021;Chen et al. 2021Chen et al. , 2022) ) and GW-like monopolar correlations (Arzoumanian et al. 2021).However, Arzoumanian et al. (2021) found these correlations were not significant as they were not robust to the solar system ephemeris and were associated with pulsar J0030 + 0451.In Section 2, we review the theoretical background required to identify and search for a general transverse polarization mode of gravity using PTAs.In Section 3, we then describe the analyses performed, both using Bayesian and frequentist approaches.Lastly, in Section 4, we present the evidence for/against the existence of ST correlations.

Theoretical Background
In this section, we will first review the basics of GW polarization modes in Section 2.1, and proceed to outline the theoretical considerations needed to predict the signature of such modes in a PTA GWB signal in Section 2.2.Finally, in Section 2.3, we will explicitly describe the model for a general transverse GWB signal, which we will later search for in Section 3 using the NANOGrav 15 yr data set.

Generalized Polarization Modes in Metric Theories of Gravity
In metric theories of gravity, there can be between two and six independent polarization modes for GWs (Eardley et al. 1973a).These modes are the "electric" components of the Riemann tensor R 0i0j , where i and j are the spatial components.These components were originally found by Newman & Penrose (1962) making use of tetrad and spinor calculus.
For the purposes of this work, we assume a coordinate system such that a null-plane GW travels along the +z-axis at the speed of light (c), where the components of the Riemann tensor only depend on the retarded time u = t − z/c.The assumptions lead to the following coefficients, which depend on combinations of the independent electric components of the Riemann tensor: We may relate these to a matrix of the GW polarization modes by Here, A + and A × represent the two tensor modes of GWs, the only two allowed by GR.The shear modes are given by A V1 and A V2 , while the scalar breathing and scalar longitudinal modes are A B and A L , respectively.Searching for the coefficients in Equation (2) allows for a theory-independent way to perform a test of gravity, without the need to be concerned with the specifics of any metric theory of gravity.We will utilize this technique in searches for a GWB using PTAs.

Pulsar Timing and Isotropic Gravitational Wave Background
A GW propagating through the Earth-pulsar line of sight will induce a change in the expected time of arrival for the pulsar's radio pulse.These perturbations were first calculated in the late 1970s (Sazhin 1978;Detweiler 1979) and have since been used to predict the GWB signature.For pulsar timing, the measured variation in the pulse TOAs can be used to calculate GW-induced residuals, R a GW , of pulsar a following the relation which is the quantity measured directly in PTAs, where z a is the GW-induced redshift.For a detailed explanation of Equation (4), refer to Arzoumanian et al. (2021) and Chamberlin & Siemens (2012).
The fractional energy density of the background is given to be where ρ GW is the energy of the gravitational wave, f is the frequency, and ρ c is the critical density necessary for a closed Universe.For the purposes of this analysis, we assume the GWB is produced by a large number of independent, weak, unresolvable sources isotropically distributed throughout the sky.Hence, the correlation between the strain functions is written as where S h ( f ) is the one-sided power spectral density of the GWB; related to The spectral characteristics of the GWB are often described via the characteristic strain, This quantity is useful, as it includes the effects of the number of cycles during the GW source in-spiral throughout the frequency band f as discussed in Taylor (2021).While several models exist for describing the nature of h c ( f ) (NG15), in this work, we will restrict ourselves to that of a power-law model for each polarization g such that where A g is a dimensionless amplitude, f yr is the reference frequency, and α g is the spectral index.Using Equations ( 6), (8), and (9), we can find the cross-correlation estimator between pulsars a and b where f L and f H are lower and upper frequencies and is the overlap reduction function (ORF), which is related to the spatial geometry of the two pulsars in relation to the Earth and we have introduced P g defined as In the above, to align with the more widely used terminology for spectral index, we have made the reparameterization γ g = 3 − 2α g .

Restriction to Transverse Modes
As discussed in Section 2.1, there exist between two and six possible independent polarization modes for a GW in metric theories of gravity.Calculating the effects of longitudinal modes requires additional steps and assessments, such as accurate knowledge of distances to the pulsars, handling the frequency dependence of the ORF, as well as having a significant number of pulsars at small-angular separations to capture the unique ORF signature of such polarization modes (Arzoumanian et al. 2021).Thus, we will restrict ourselves to the three transverse modes, A + A × , and A B , for the purposes of this Letter.
Given only the transverse tensor and scalar mode, we may generalize Equation ( The ORFs for the tensor-transverse (TT) and ST modes have been calculated previously in Chamberlin & Siemens (2012) and Gair et al. (2015), with ξ ab being the angular distance on the sky between pulsars a and b and A plot of the transverse ORFs as a function of angular separation is shown in Figure 1.Where ab TT G is represented by the more widely known Hellings and Downs (HD) curve (Hellings & Downs 1983), and henceforth the TT mode will be represented by HD.With this structure in hand, we may proceed to the analysis of the NG15yr data and the investigation of the general transverse modes.

Searches for a General Transverse GWB in the
NANOGrav 15 yr Data Set In this section, we complement our previous work in NG15 by analyzing the NANOGrav 15 yr data set for statistical significance of HD plus ST correlations.We first describe the pulsar noise modeling in Section 3.1, and then we present the results of the Bayesian and frequentist analyses in Sections 3.2 and 3.3, respectively.We will take an agnostic approach to the mixing between the HD and the ST polarization modes of gravity by allowing each mode to possess its own independent power spectral density as suggested by Equation (13).

Noise Modeling Details
Through individual pulsar analyses, we obtain posteriors for both the red and white noises intrinsic to each pulsar.Red noise has more power at lower frequencies.We model the intrinsic pulsar red noise as a power law with variable characteristic amplitude and spectral index following Equation (12).We model the power spectra using frequency bins from 1/T obs to 30/T obs to cover a frequency range in which pulsar noise transitions from red-noise-dominated to white-noise-dominated and for T obs being the longest observational baseline among the considered pulsars in the data set.The white noise is described by three parameters: a linear scaling of TOA uncertainties, noise added in quadrature to the TOA uncertainties, and noise common to a given epoch at all frequency subbands.These parameters are called EFAC, EQUAD, and ECORR, respectively, and are set to their fixed values in NG15.For detailed explanations of these parameters, refer to Agazie et al. (2023b).
In addition to pulsar-intrinsic noise, we also include a common red noise, which is a red noise process that is shared among all pulsars.The GWB is expected to appear as a common red noise process that is spatially correlated across pulsars.The common red noise process is also modeled as a power law but uses three different overlap-reduction functions corresponding to three different kinds of common red noise processes: CURN, HD, and ST.The common red noise is modeled using frequency bins from 1/T obs to 14/T obs following NG15.The CURN treats the common red noise process as spatially uncorrelated (i.e., ab ).The HD and ST models include additional spatial correlations, which are shown in Figure 1 and Equations (14a) and (14b).
The upper and lower bounds of the model parameters we use are shown below.Note that the subscript "int" refers to the intrinsic red noise processes, while the subscript g refers to the common red noise process (i.e., CURN, HD, or ST).

Bayesian Analyses
Our Bayesian analyses follow Section 3.1 as well as NG15.In short, in terms of a likelihood function, all the various noise modelings follow from (Johnson et al. 2023) and we then use the Woodbury matrix identity to invert this covariance matrix.We find In the above, F is a matrix with alternating columns of sine and cosine components representing a discrete Fourier transform of the red noise processes, D is covariance matrix for the white noise parameters, and f is the covariance matrix of the red noise components.We use Bayesian analyses to compare several models of interest via Bayes factor estimation (Figure 2), and to obtain posterior distributions for A log g 10 and γ g for HD and ST signals (Figure 3).While we include an ST-only model in our analysis, it is important to note that all metric theories of gravity must include the Einsteinian polarization modes of gravity.In Figure 2, we observe that correlated Bayesian models are preferred over the uncorrelated model.The most favored model is a GWB with HD correlations with a Bayes factor of 200.When ST is modeled alongside HD, Bayes factors are uninformative given they are on the order of unity when compared to each correlation alone.
We can use the transitive nature of Bayes factors as a consistency check of our results.For instance, going around the bottom half of Figure 2 we can take the Bayes factor of ST/ CURN and multiply it by the Bayes factor of HD/ST to obtain the Bayes factor of HD/CURN.This results in 90 × 2.2 = 198, which is consistent with the Bayes factor for HD/CURN of ∼200 we obtained by directly comparing those two models.
We note that in Figure 3 when fitting for one correlation signature both HD correlations and ST correlations are able to explain the total signal.This agrees with the large Bayes factors favoring these models over CURN.However, the recovered power spectral estimates for ST are poor when modeled alongside HD.To check the consistency of the power spectral estimates we see that A log 14.17  .We see that values for CURN and HD are more consistent with each other.While the ST spectral estimates do overlap with the median of the CURN spectral estimates, we observe the 68% credible region for γ ST and A log 10 ST expand over about 43% and 31% of the prior region, respectively.Therefore, the addition of the ST correlation yields no additional information and we see that the HD signal in this model explains most of the total signal.

Optimal Statistic Analyses
The optimal statistic (Anholm et al. 2009;Chamberlin et al. 2015) allows for a robust and computationally inexpensive analysis of the correlation content of a PTA data set.The amplitude and the uncertainty of the pairwise cross correlations are estimated by maximizing the ratio of the likelihood of the fiducial GWB over the noise-only model.The fiducial model contains a GWB signal along with intrinsic red and white noise components while the noise-only model includes the intrinsic noises and a common uncorrelated red noise process.We have employed the noise-marginalized version of the optimal statistic technique in which 10 4 random draws from the posteriors of all of the model parameters of a CURN model are used to estimate the required power spectra.Additionally, since our goal is to search for a general transverse GWB signal in which two nonorthogonal types of correlations might simultaneously exist in the data set, a chi-squared statistic fitting for both HD and ST correlations is used to find the optimal estimators of the signal-to-noise ratio (S/N) and the amplitude of the correlated signal (A g ˆ) for each polarization mode.Note that A g ˆdiffers from A g as the former is an optimal estimator of the latter.See Vigeland et al. (2018) and Sardesai & Vigeland (2023) for more details.
The estimated amplitudes for single-component (SCOS) and multicomponent (MCOS) noise-marginalized optimal statistics are shown in Figure 4. We see that for SCOS the amplitude reconstruction of the HD correlations is in excellent agreement with the CURN amplitude posterior, whereas the estimated amplitude for ST correlations is only consistent with CURN.For MCOS, the correlations are fit for simultaneously so the total power is split and the fit correlations are shifted toward smaller amplitude values and less consistent with CURN.
However, the HD correlations explain most of the total CURN signal when both correlations are present.
The S/Ns for SCOS and MCOS are shown in Figure 5.We note that the median S/Ns for HD correlations, 5.0 and 3.5 for SCOS and MCOS, respectively, are larger than the median S/ Ns for ST correlations, 4.6 and 3.0 for SCOS and MCOS, respectively.The difference in median S/N values is not significant as the medians lie within the interquartile ranges of each other.While the S/N values are similar for HD and ST correlations, as noted before, the consistency of the estimated HD amplitude with CURN suggests that HD, not ST, correlations make up most of the common red noise process.

Discussion
NANOGrav's 15 yr data set shows compelling evidence for quadrupolar HD interpulsar correlations.In this work, we explored the possibility of deviations from the HD curve caused by the presence of an additional ST mode.
Our Bayesian analyses show the Bayes factor for HD over ST is ∼2, and the Bayes factor for a model with both correlations compared to a model with just HD is ∼1.These results are largely consistent with a similar study by Chen et al. (2023), in which they searched NANOGrav's 15 yr data set for nontensorial GWBs on a similar timescale to the work presented here.Taking the spectral parameter recovery into account, as in Figure 3, we found each correlation, when fit for individually, is in agreement with CURN.We also found more informative A log g 10 and γ g recovery for HD than ST, and HD parameters show better agreement with CURN spectral parameters when correlations are included together.The analyses in this Letter, as well as those in Bernardo & Ng (2023c) and Chen et al. (2023), do not rule out the possibility of ST correlations in our data.However, our analysis also shows no statistical need for an additional stochastic process with ST correlations.This is also the case for our frequentist analyses.When fitting the interpulsar correlation data for a single correlation signature, we find that HD correlations completely account for the total signal due to the amplitudes consistency with CURN, but ST correlations are only somewhat consistent.When we fit for both correlations simultaneously, we still see that HD correlations are able to explain most of the total signal.For the Einsteinian polarization modes with HD correlations are present in all metric theories of gravity.Thus, even though we do not find a convincing Bayes factor favoring only HD correlations over only ST correlations and they have similar S/ N values in our frequentist analyses, there is no metric theory that predicts only ST GWs.Metric theories do allow ST correlations to be present alongside HD correlations, but we did not find strong evidence in favor of HD plus ST correlations over only HD correlations.In addition, we no longer report higher S/Ns and Bayes factors for ST correlations as we did in Arzoumanian et al. (2021).We have seen a larger increase in favor of HD correlations than ST correlations in both S/Ns and Bayes factors.These changes are consistent with simulations in Arzoumanian et al. (2021) and, with no evidence indicating otherwise, we expect this trend to continue with additional data.
We also performed dropout analysis tests, similar to what was done in Arzoumanian et al. (2021), to determine if particular pulsars play a role in the observed ST significance.We found that J0030 + 0451 and J0613−0200 are responsible for a majority of the ST significance.When we remove these two pulsars from the analysis, we find that the Bayes factor for HD/CURN increases to ∼600, while the Bayes factor for ST/ CURN is reduced to ∼30.We suspect improved noise modeling (as used in Falxa et al. 2023 andAgazie et al. 2023b) on these and other pulsars will shed some light on this, and we leave this for future work.
Other recent work (Allen & Romano 2023; Allen 2023) has shown that the HD correlation signature has a cosmic variance, and accounting for the cosmic variance improves the ability to see correlations for the GWB in noisy data (Bernardo & Ng 2023e).While searches have included cosmic variance (Bernardo & Ng 2022, 2023b, 2023e, 2023f), this concept is not addressed within this manuscript.It introduces another aspect to consider when searching for the GWB as we will also need to account for the variance of the HD curve.Impacts of the cosmic variance will be considered in future work.
Future prospects for performing tests of gravity using PTA data are compelling.Large observational baselines as well as the addition of more MSPs to the observing array will enable more robust and sensitive searches for additional GW polarization modes.In addition to PTAs, future experiments in other regions of the GW frequency spectrum may also provide insight into alternative polarizations.Callister et al. (2017) and Wang & Han (2021) have shown the capabilities to make detections of alternative GW polarizations in their respective parameter spaces.
In this work, we reported on one test of gravity in which we searched for evidence for an ST polarization mode.While we did not find substantial evidence for or against this mode, the situation may change in the future due to the nature of the PTA data sets.It is also important to note that a number of the observed pulsars are dominated by white and intrinsic red noise processes, which could be suppressing a GW-sourced signal.For the case of a GWB, as we obtain more data on these pulsars, we will be able to provide more definitive answers about the possibility of the existence or absence of additional polarization modes of gravity.

Figure 1 .
Figure 1.A plot of transverse ORFs as a function of angular separation.The blue curve describes the Hellings and Downs curve, which is produced by the TT polarization mode, while the orange curve describes the shape of the correlations induced by the ST polarization mode of gravity.
Agazie et al. (2023b) for a more detailed explanation of the noise modeling adapted for the analyses of the NANOGrav 15 yr data set.

Figure 2 .
Figure 2. Bayes factors for various model comparisons between ST, HD, and CURN.Overall, the HD model is preferred over CURN and ST.Modeling ST alongside HD gives about equal odds over HD and ST only.All model comparisons are agnostic with respect to the spectral index of each model.See Section 3.1 for more details.The uncertainties are estimated using bootstrapping and Markov model techniques of Heck et al. (2019).An ST-only model is made for comparisons even though a metric theory must have HD present.All Bayes factors are presented as the model at the end of an arrow over the starting model of an arrow.For example, for the arrow pointing from CURN to ST, the values are Bayes Factors for ST/CURN.

Figure 3 .
Figure 3. (Left) Bayesian probability posterior distributions of A log g 10 and γ g from an HD correlated model (blue) and an ST correlated model (orange) showing the 1σ/2σ/3σ credible regions.Plotted for comparison are CURN (gray) posterior distributions for A log g 10 and γ g parameters.Each correlated signal is able to explain the total signal.(Right) Bayesian probability posterior distributions of A log g 10 and γ g for HD (blue) and ST (orange) from an HD+ST correlated model showing 1σ/2σ/ 3σ credible regions.Plotted for comparison are CURN (gray) posterior distributions for A log g

Figure 4 .
Figure 4. Distributions of the recovered amplitudes from single-component noise-marginalized optimal statistic (SCOS; left) and multicomponent noise-marginalized optimal statistic (MCOS; right) for HD (blue) and ST (orange) correlations.Additionally, CURN (gray) is plotted for comparison to determine the consistency with the common red noise process.

Figure 5 .
Figure 5. Distributions of the S/N for HD (blue) and ST (orange) correlations from the multicomponent-(MCOS) and single-component (SCOS) noisemarginalized optimal statistic techniques.The dashed lines correspond to the median values and the first and third quartiles.The median S/N values for HD are greater than ST but are within their interquartile range.
11) (O'Beirne et al. 2019; Arzoumanian et al. 2021) as It is worth pointing out that the effect of dipole radiation of binary sources in non-GR metric theories of gravity (O'Beirne et al. 2019) is accounted for by treating the spectral index γ g as a free parameter in our statistical models.