Dissecting the stochastic gravitational wave background with astrometry

Astrometry, the precise measurement of star motions, offers an alternative avenue to investigate low-frequency gravitational waves through the spatial deflection of photons, complementing pulsar timing arrays reliant on timing residuals. Upcoming data from Gaia, Theia, and Roman can not only cross-check pulsar timing array findings but also explore the uncharted frequency range bridging pulsar timing arrays and LISA. We present an analytical framework to evaluate the feasibility of detecting a gravitational wave background, considering measurement noise and the intrinsic variability of the stochastic background. Furthermore, we highlight astrometry's crucial role in uncovering key properties of the gravitational wave background, such as spectral index and chirality, employing information-matrix analysis. Finally, we simulate the emergence of quadrupolar correlations, commonly referred to as the generalized Hellings-Downs curves.


Introduction
Recent breakthroughs in pulsar timing arrays (PTAs) have heralded a transformative era for the observation of gravitational waves (GWs) at nano-Hertz (nHz) frequencies.These advancements leverage the precise timing information obtained from pulsars, enabling a galactic-scale gravitational wave detector.This innovative approach first led to the detection of a common-spectrum process by the Nanohertz Observatory for Gravitational Waves (NANOGrav) collaboration [1], a finding subsequently corroborated by the Parkes PTA, European PTA, and International PTA [2][3][4].More recently, multiple PTA collaborations have found evidence for the quadrupolar angular correlation of these signals on the sky [5][6][7][8].These collective findings consistently support the presence of a quadrupolar correlation function known as the Hellings-Downs curve [9], a pivotal characteristic of gravitational-wave induced timing residuals.Furthermore, the amplitude and spectral index of the inferred stochastic gravitational wave background (SGWB) align broadly with predictions derived for a cosmological population of supermassive black hole binaries (SMBHBs) as the gravitational wave sources [10,11].Nonetheless, while not definitively ruled out, certain cosmological sources may still provide potential explanations for the observed SGWB.
On a parallel front, it is essential to note that GWs not only perturb photon travel times but also deflect their paths.This phenomenon serves as the basis for astrometric detection of GWs, which relies on measuring correlated wobbling movements of stars on the sky [12][13][14][15][16]. Astrometry represents another promising avenue for the detection and characterization of gravitational waves, providing a complementary perspective to the PTA approach [17][18][19][20][21][22][23][24][25][26][27][28][29][30].While astrometry currently has less constraining power than PTAs [23,27,29], complete datasets from space-borne astrometry missions like Gaia [31], including full time series, have the potential to provide sensitivity comparable to PTAs.The next-generation astrometric observations have the potential to exceed the sensitivity of PTAs [26].
Compared to PTAs, astrometry boasts several noteworthy complementary advantages.The distinct response functions of PTA and astrometry make them sensitive to incoming GWs from different directions as pulsars or stars are more concentrated toward the direction of the Galactic Center [17].Additionally, the strain sensitivity of astrometry remains nearly constant across the frequency spectrum [17], unlike the linear decrease in PTA sensitivity toward higher frequencies.This opens a new window for GWs, exploring uncharted frequency ranges lying between PTAs and the Laser Interferometer Space Antenna (LISA) band [32].For instance, the Nancy Grace Roman Space Telescope (Roman), with its significantly higher observing cadence, can indeed extend the frequency range to above 10 −4 Hz [28,[33][34][35].Moreover, the presence of parity-odd correlations among astrometric observables or PTAastrometry cross-correlations can reveal the existence of a chiral component of SGWB [21,30,36].Hence, it is imperative to explore in detail the prospects of astrometry.
This study is dedicated to forecasting the potential of astrometry in discovering SG-WBs and characterizing their properties with future data releases, in conjunction with PTAs or without.We establish a framework for predicting the feasibility of astrometric SGWB detection and the resolution of key SGWB parameters.These parameters, including the normalization of characteristic strain, spectral index, and chirality, are crucial for comprehending the distribution of SMBHBs, their potential interaction with the environment, and the presence of any sub-leading cosmological sources.The vector nature of astrometry observables introduces various options for cross-correlations [19,21], in addition to the redshift-only correlation of PTA.Each of these correlations possesses unique quadrupolar correlation functions [19,21] and corresponding variances due to the stochastic nature of GWs.Identifying them will not only provide a cross-check of the PTA Hellings-Downs curve at nHz but may also uncover a GW signal at higher frequencies.The paper's structure is organized as follows: In Sec. 2, we review the basics of both PTA and astrometry for the detection of the SGWB.Moving to Sec. 3, we analytically calculate the sensitivities for various cross-correlation choices associated with PTA and astrometry and evaluate the resolution of key SGWB properties.Section 4 delves into the simulation of spatial correlations, specifically exploring the generalized Hellings-Downs curve, and offers a theoretical insight into intrinsic variance of SGWBs.Finally, in Sec. 5, we draw our conclusions and discuss our findings.

Responses and Angular Correlations for PTA and Astrometry
GWs induce perturbations in the paths of photons along geodesics.These perturbations give rise to two distinct categories of observable phenomena: shifts in the temporal arrival of photons and the proper motion of their sources across the celestial sphere [16].PTAs, functioning as an exceptional network of cosmic clocks, possess the remarkable capability to precisely measure the arrival times of radio pulses from distant pulsars.On the other hand, astrometry is dedicated to the precise determination of the positions and motion of stars across the celestial sphere.The shifts depend on both the metric perturbations at the observation point (Earth term) and on the emission sources (pulsar or star terms).However, the latter can usually be disregarded when the distance between the two points significantly exceeds the wavelength of the GWs [19], or treated as noise when correlations among different baselines result in only the Earth term being coherently summed up.For the purposes of this study, we will focus solely on the Earth term.
The received GW strain at a specific location can be represented as a sum of frequency modes: (2.1) In the above equation, f , Ω and P label the frequency, incoming direction, and polarization mode of the GW, respectively.In this work, we will only consider polarization modes within the framework of Einstein's gravity.The strain amplitude in the frequency domain is given by h P (f, Ω), and ϵ P ij ( Ω) is the polarization basis tensor satisfying ϵ P ij ( Ω)ϵ ij P ′ ( Ω) = 2δ P P ′ .The time-domain strain imposes a real condition, necessitating that h P (f, Ω) = h P (−f, Ω) * and ϵ P ij ( Ω) * = ϵ P ji ( Ω) for linear polarization basis P = +/×.The frequency-domain signals from PTAs and astrometry can be universally expressed as [36]: Here, we introduce X ≡ {δz, δx} to encompass both the photon redshift δz from PTAs and the proper motion on the celestial sphere δx from astrometry.The subscript a designates the a-th pulsar/star, with the line-of-sight direction denoted as na .The redshift and astrometric response functions are elucidated in Ref. [16] as follows: Here, we utilize a Cartesian coordinate system where the components are labeled by i, j and l.Notably, the timing residual signal represents the time integral of the redshift δz, introducing an additional factor of 1/(2πf ).For a Gaussian, stationary GW background, the two-point correlation function of h P is as follows: where P P P ′ (f, Ω) represents the power spectrum of correlation between the P mode and the P ′ mode.When the SGWB exhibits isotropy, the power spectrum matrix can be parameterized as [37]: for P/P ′ ∈ {+, ×}, where the real quantities I(f ) and V (f ) represent the total intensity and the circular polarization, respectively.The isotropic SGWB, as defined in Eqs.(2.4) and (2.5), results in correlations between two received signals from a pair of pulsars, stars, or pulsar-star pairs, as follows: (2.6)For PTA with redshift correlations, Eq. (2.6) leads to the well-known Hellings-Downs curve [9]: where the function Γ z (θ ab ) depends solely on θ ab ≡ na • nb due to the rotational invariance of an isotropic SGWB.
Correlations involving astrometric motions of a pair of stars or a star and a pulsar can be categorized into directions that are parallel (ê a || and êb || ) and perpendicular (ê ⊥ ) to their great arc [19], defined as: (2.8) Using these definitions, the correlations in Eq. (2.6) can be simplified to the following expressions [19,21,30,36]: where the dimensionless correlation functions satisfy (2.10) These functions are commonly known as generalized Hellings-Downs curves.

Gravitational Wave Signal in Spherical Harmonic Space
An alternative representation of SGWB signals in PTAs and astrometric observation is achieved through the use of spherical harmonic space, as demonstrated in prior works such as Refs.[21,28,[38][39][40].The key advantage of this representation lies in its ability to diagonalize both the signals and the SGWB-induced variances, simplifying the definition of estimators.In this formalism, both the GW-induced redshift and angular deflection can be expressed as discrete summations over the harmonic basis: where Y ℓm represents spherical harmonic functions, and Y E ℓm and Y B ℓm are the E-and Bcomponents of vector spherical harmonic functions, respectively.Correspondingly, z ℓm (f ), E ℓm (f ), and B ℓm (f ) are their respective expansion coefficients.Employing the orthogonality of the spherical harmonic basis, we can determine these components as follows: (2.12) It is important to note that in realistic observations, a finite number of pulsars and nonuniform distributions of pulsars/stars should be taken into account when reconstructing these components.These factors can lead to a mixture of modes and consequently introduce noise [41].However, recent simulations using mock data have shown that the influence of this mixture is negligible [40].
To simplify notation, we introduce X ℓm ≡ {z ℓm , E ℓm , B ℓm } to represent all sphericalharmonic components.With this representation, we can construct the rotationally invariant power spectra: These power spectra satisfy Here, T represents the total observation time, which is used to account for the factor δ(0) → T in the discrete frequency domain.
A convenient way to calculate X ℓm is by using the total angular momentum (TAM) decomposition of GW strain [42], which includes both transverse polarization modes, denoted as α ∈ {T E, T B}.Consequently the GW-induced (h) spherical harmonic components can be expressed as [21]: for ℓ ≥ 2. Here, F X,α ℓ are projection factors derived for various combinations of X and α [21].The respective strain amplitudes, denoted as h α ℓm (f ), exhibit correlations that can be parameterized as: This assumes an isotropic SGWB with vanishing linear polarization components.The power spectra matrix P αβ is structured as [37,43]: for α/β ∈ {T E, T B}.
The ensemble average of C XX ′ ℓ (f ) comprises an independent sum of GW-induced signals (h) and measurement noise (n): The signal part is directly derived from Eq. (2.14) and (2.15): More explicitly, each component of Eq. (2.18) is directly related to either the total intensity I(f ): or the circular polarization V (f ) [25,30,36]: (2.20) Here, we define On the contrary, the Gaussian noise inherent in each measurement results in an ℓindependent noise spectrum [39].In the case of a uniform distribution of N X pulsars or patches (N E = N B ≡ N δx ) on the celestial sphere, the noise spectra in harmonic space can be expressed as: Here, S (n) X (f ) represents the noise spectra from each pulsar/patch.For astrometry, these spectra satisfy S (n) δx (f ) is the measurement noise of the average proper motion of the patch.The rotational invariant noise spectra directly follow from Eq. (2.22): (2.23)

Dissecting the Stochastic Gravitational Wave Background
In this section, we assess the sensitivity of both PTA and astrometry utilizing rotational invariant power spectra C XX ′ ℓ .As highlighted earlier, these observables derive considerable advantage from their diagonal nature in harmonic space.Our exploration commences with the formulation of estimators built upon C XX ′ ℓ , followed by a thorough analysis of their properties.Subsequently, we employ the information matrix [44][45][46] to gauge the sensitivity towards total intensity, spectral index, and chirality.

Rotational Invariant Estimators
In the frequency domain, discrete frequencies f k span from 1/T to 1/(2∆t), where ∆t is the cadence of observation.We use the integer k to label measurements in a certain frequency bin such that The total intensity I k can be estimated from XX ′ = zz/EE/BB/zE, while the circular polarization V k arises from the parity-odd observables zB/EB.Each estimator is directly constructed from Eq. (2. 19) and (2.20): whose ensemble averages are either I k or V k .Given that each power spectrum in harmonic space is independent, the signal-to-noise ratio (SNR) of these estimators in a given frequency bin receives contributions from all achievable ℓ-modes: Here, ℓ max ∼ N X /2 is the highest observable ℓ for the constructed power spectra.The denominator corresponds to the variance of the estimator in Eq. (3.1), derived using Isserlis' theorem [47] for Gaussian fields X ℓm .
To gain a quantitative understanding of Eq. (3.2), we depict the SNR 2 k as a function of 1 for the I k estimators.These estimators include PTA-only (zz) in blue, astrometry-only (EE or BB) in orange, and PTA-astrometry cross-correlation (zE) in green.Two scenarios are presented with ℓ max = 10 (solid lines) and 50 (dashed lines).The PTA-only case (zz) aligns with results previously obtained using cross-correlations in separation-angle space [48].To better understand the behaviour of the results plotted in Fig. 1, we consider the plot categorized into three regions: The measurement noise is subdominant, and the data variance is predominantly determined by the stochastic variance of GW.The saturated value of SNR 2 k is ℓmax ℓ=2 (2ℓ+1) = (ℓ 2 max + 2ℓ max − 3)/2, which is universal for all cases.
• Intermediate signal: Between the weak and strong signal limits, the contribution of the lowest few ℓ-modes of (C XX ′ (h)ℓ,k ) 2 surpasses that of the measurement noise.As C XX ′ (h)ℓ,k decreases with higher X,k /N X ) for I k estimators, encompassing PTA-only correlation zz (blue), astrometry-only correlation EE or BB (orange), and PTAastrometry cross-correlation zE (green).Two options for the highest ℓ-mode constructed, ℓ max = 10 (solid lines) and 50 (dashed lines), are considered.In the case of zE correlation, we assume ℓ, the transition commences at lower values of ℓ ≥ 2 until the ℓ max -mode becomes dominant.In this region, we introduce ℓ eff ≤ ℓ max to label the highest ℓ-mode whose power spectra dominates over the measurement noise.Consequently, the SNR 2 k becomes The two cases with different ℓ max exhibit nearly identical behavior in the weak signal regions, start to deviate in the middle of the intermediate signal regions, and eventually reach saturation at their respective maximum SNR values.
For estimators of circular polarization V k , the variance incorporates C XX (h)ℓ,k , which is proportional to I k .Consequently, their maximum SNR is typically suppressed by the ratio of circular polarization V k /I k .

Sensitivity and Parameter Estimation
In addition to exploring the total intensity and circular polarization, many aspects of the SGWB remain to be unraveled.In this section, we establish an analytical framework to dissect the SGWB.This includes sensitivity estimation and the assessment of various pertinent parameters, all based on the information matrix [44].
We introduce a vector X ℓm,k defined as follows: encompassing both redshift and astrometric observables.Assuming Gaussianity and stationarity for both the SGWB and the measurement noise, each X ℓm,k follows an independent complex multivariate normal distribution: Here, Σ ℓ,k is the 3 × 3 covariance matrix: where all the components are defined across Eqs.(2.19), (2.20) and (2.23).
Given the statistical independence of X ℓm,k for different m, ℓ, and k, the joint probability can be expressed as the product of individual probabilities.Consequently, the ln-likelihood function is given by: where O represents the model parameters, and Σ ℓ,k (O) is the covariance matrix constructed from the corresponding model parameters.Importantly, Σ ℓ,k (O) adheres to the condition Σ ℓ,k (O truth ) = Σ ℓ,k , as defined in Eq. (3.5), where O truth represents the true parameters.The ability to estimate the parameters is gauged by the information matrix [44]: Here, O i represents the i-th parameter, and ⟨. . .⟩ denotes the ensemble average over X ℓm,k .The inverse of the information matrix I provides the uncertainties for parameter estimation: In cases where the only model parameter of interest is the GW spectrum intensity I k , the SNR of this amplitude is given by: It is crucial to note that the information-matrix analysis is accurate primarily when the SNR is robust [45,46,49].Therefore, we consistently choose SNR ≥ 1 as the threshold to ensure the validity of the subsequent discussion.

Power-law Stochastic Gravitational Wave Background
In this study, we examine power-law SGWB models and investigate the possible presence of a circularly polarized component.The power-law model is a prevalent approach to characterize the power spectrum of SGWB.Many SGWB models arising from astrophysical or primordial origins can be effectively approximated by power laws in the nHz frequencies.The intensity spectrum of this model is expressed as follows: where f ref represents the reference frequency, commonly chosen as f ref = 1/yr in PTA observations, and α denotes the spectrum index.The dimensionless characteristic strain can then be defined as The spectral index of the SGWB can deviate from a constant value.At the low-frequency end of the nHz spectrum, potential interactions with the environment [52][53][54] or orbital eccentricities of SMBHBs [55] lead to a turning of the spectrum slope [10].The high-frequency turning, occurring around ∼ 10 −6 Hz for SMBHBs with masses of ∼ 10 9 M ⊙ , happens as SMBHBs approach the merger phase [56].However, in the frequency range we consider, the low-frequency deviation is expected to be small, while the high-frequency contribution is sub-leading in sensitivity.Thus, we focus on a constant spectral index in the following.
Another feature of the SGWB is chirality, parameterized by macroscopic circular polarization.Cosmological models, such as those described in Refs.[57][58][59][60][61], can generate chirality through parity-violating interactions.A finite sum of nearby SMBHBs can also produce a random fraction of chirality [62][63][64].As C zz (h)ℓ,k in Eq. (2. 19) is dependent solely on I k and lacks any V k dependence, measuring the isotropic circular polarization map using PTA-only observations is not possible.It is noteworthy that PTA can still probe anisotropic circular polarization [65][66][67][68], which is beyond the scope of this study.
Fiducial SGWB Model from PTA Observation Recent observations by NANOGrav [5], PPTA [7], EPTA [6], and CPTA [8] suggest that the SGWB signal is consistent with that produced by SMBHBs.Assuming that the nearly-circular orbit evolution is predominantly driven by GW emission, the SGWB spectrum can be well-described by a power-law model with α = −2/3, despite a potential deviation at the low-frequency end due to environmental effects or eccentric orbits [10,69,70].With α = −2/3 fixed, the strain found by NANOGrav (NG) is given by where A NG ≃ 2.4 × 10 −15 [5].Utilizing Eq. (3.11), the power spectrum of our fiducial SGWB model has a reference intensity In the subsequent sections, the fiducial model will be employed to compare sensitivities across different observational channels.Any deviation from α = −2/3 can be interpreted as environmental effects or primordial origins.

Power-law Model Parameters
The power-law model involves two parameters, denoted as O = {log 10 A, α}.The derivatives with respect to these parameters are obtained through the chain rule: where we define the dimensionless factor ∆ f ≡ 1/(T f ref ).The information matrix is expressed as 10) ln(k∆ f ) ln (10) ln where σ k ≡ σ(I k ).Note that the matrix inside the summation may appear singular, but the overall matrix after summation is not.The uncertainties of {log 10 A, α} are given by respectively.An alternative parameterization using O = {log 10 A, γ}, where α ≡ (3 − γ)/2, results in a straightforward rescaling of σ(γ) = 2σ(α).The assessability of both A and α is significantly dependent on the total SNR.However, the weighted summation of individual frequency bins adds complexity to the equations.In the following, we will directly calculate σ −2 (α) and σ −2 (log 10 A) and examine their behavior in different observation channels.

Pulsar Timing Arrays
In PTA-only observations, the ln-likelihood defined in Eq. (3.6) exclusively involves the redshift z ℓm,k , obtained by marginalizing over all E ℓ,m and B ℓ,m : As discussed previously, this likelihood is not sensitive to V k .According to Eq. (3.8), the uncertainty of each {I k } is given by (3.17) The total SNR follows Eq. (3.9) which aligns with the estimator using C zz ℓ,k in Eq. (3.2).Given various ξ z k value ranges, the SNR can be categorized into weak, intermediate, and strong signal regions, as discussed in Sec.3.1.

PTA Sensitivity to the Power-law SGWB
We proceed to estimate the sensitivity to the power-law SGWB in realistic PTA observations.The Gaussian noise for timing residual consists of a white noise component and a red noise one, relatively well-fitted by the powerlaw spectrum [71].By transitioning from timing residual (TR) to redshift, the noise spectrum for a single pulsar can be expressed as Here, S (r) z and S (w) z are the corresponding noise components, σ TR is the timing residual uncertainty of each measurement for a pulsar, ∆t is the cadence of the observation.
We consider recent NANOGrav and future SKA observations.The red noise for NANOGrav is fit to be S (r) z = 1.3 × 10 −28 yr and γ r = 0 [5,71], while for SKA, we assume there is only white noise.The benchmark parameters for the two observations are taken from [71] and [72,73], respectively: Here, T obs is the total observation time.Comparing our noise model for NANOGrav with the fiducial SGWB model defined in Eq. (3.10), we find that the SGWB dominates the noise in the lowest 5 bins, consistent with NANOGrav's result [5].With N z = 50 pulsars for NANOGrav, this leads to ℓ max = 5, while for SKA, ℓ max = 10.
In the left part of Fig. 2, we depict the SNR distribution of the power-law SGWB as a function of log 10 A and α for both PTA observations.The fiducial SGWB model yields an SNR ≃ 4.0 for NANOGrav (NG) and SNR ≃ 34.1 for SKA.Notably, we observe that above SNR ≈ 80, the NANOGrav sensitivity reaches a saturation phase, where a stronger A no longer improves the SNR due to the intrinsic variation in SGWB, consistent with the result in Ref. [48].On the other hand, for SKA with a larger value of ℓ max and more sensitive frequency bins, the threshold for saturation is significantly higher.
We also compute the estimation uncertainty of the power-law model parameters, log 10 A and α, using Eq.(3.8), as illustrated in Fig. 3.For the fiducial SGWB model, the NANOGrav case yields σ(log 10 A) ≃ 0.25 and σ(α) ≃ 0.24, consistent with their data analysis [5].In comparison, the SKA demonstrates superior resolution over NANOGrav by factors of 21 and 14 for the two parameters, respectively.Akin to the total SNR observations, the uncertainty ceases to improve beyond the saturation phase for NANOGrav.

NANOGrav N o C o n s t r a i n t Gaia
N o C o n s t r a i n t

Astrometry
For astrometry, we have two independent measurements, X ℓm,k = {E ℓm,k , B ℓm,k }, resulting in a covariance matrix: Here, We insert the covariance matrix Eq. (3.21) into the likelihood function Eq. (3.6), labeling it as ln L (ast) .
We first calculate the SNR of the total intensity The total SNR is derived from Eq. (3.9): This expression is equivalent to the sum of estimators using C EE ℓ,k and C BB ℓ,k in Eq. (3.2).The SNR is very similar to that of PTA in Eq. (3.17), despite the additional factor B 2 ℓ ∝ 1/(ℓ(ℓ + 1)), which makes higher ℓ-modes more suppressed.Consequently, while astrometry can accommodate a significantly higher ℓ max due to the larger number of stars compared to pulsars, redistributing stars into an appropriate number of patches will not diminish sensitivity.
Astrometric Sensitivity to the Power-law SGWB Due to the large number of observed stars in astrometry, an efficient strategy, without sacrificing sensitivity, involves dividing the celestial sphere into various patches [17,18].The distribution of these patches can be realized through HEALPix [74].Each patch combines the spatial deflections of stars on the celestial sphere into an average δx.For N S stars evenly distributed in N δx patches, the two-sided noise power spectrum of each patch is given by Here, σ θ represents the uncertainty of each measurement for a star.The SNR contribution from a specific ℓ-mode is characterized by ξ δx k in Eq. (3.23), which is independent of the number of patches N δx .On the other hand, N δx determines ℓ max , influencing SNR in the strong signal region.Thus, a reasonable choice of N δx is necessary in astrometric observation.
Gaia [31] and upcoming missions such as Roman [28,[33][34][35] and Theia [75,76] play pivotal roles in astrometric observations.The Gaia mission, having measured the proper motion of over ∼ 10 9 stars and ∼ 10 6 quasi-stellar objects (QSOs) for more than 10 years, has provided a rich dataset.The QSO data from Gaia has been leveraged to constrain ultralow (≪ nHz) frequency gravitational waves [23,27,29].In this study, we consider the full dataset from Gaia upon its release, including its comprehensive measurements.The proposed mission Theia boasts significantly improved resolution but with a small field of view [75,76].We anticipate that the next-generation upgrade of Gaia, which we abbreviate as 'XG-Gaia', will achieve µas-resolution while maintaining its cadence and the number of observed stars [77].The benchmark parameters for the two astrometric missions are listed as follows [17,77]: Note that we adopt a conservative estimate for the number of stars compared to Ref. [26].
We evaluate the SNR for various measurements of the power-law SGWB model, assuming uniform distribution of stars across the celestial sphere, each possessing identical measurement properties.While a total of N S = 1.5 × 10 9 stars would allow for an ℓ max ∼ 10 4 , we opt for N δx = 40000 for Gaia and XG-Gaia with ℓ max = 200, as further increasing the number no longer significantly enhances the SNR within our parameter space of interest.
In the right panel of Fig. 2, we depict the SNR distribution of Gaia and XG-Gaia using Eq.(3.23) concerning power-law model parameters.Generally, PTA's sensitivity exhibits a more pronounced change in terms of α, attributed to its more tentative decrease in sensitivity at higher frequencies.The fiducial SGWB signal yields an SNR ≃ 1.5 for Gaia, compared to 4.0 in NANOGrav, indicating Gaia resides in the marginal region for cross-checking PTA discoveries.On the other hand, XG-Gaia achieves an SNR ≃ 43.9, higher than SKA, benefiting from contributions at higher ℓ and frequency modes.
In Fig. 3, we compare the resolution of power-law parameters in astrometry with those in PTA.The resolution approximately follows the distribution of 1/SNR.XG-Gaia, with resolutions below 10 −2 for both amplitude and spectral index, can provide an exceptionally precise dissection of the SGWB spectrum.This high precision is pivotal for gaining insights into the distribution and evolution of SMBHBs, including potential environmental effects.
The Roman telescope, distinguished by its much higher cadence, demonstrates sensitivity to the SGWB in the frequency band situated between PTA and LISA [28,35].The associated frequency range of Roman offers valuable insights into whether the SGWB displays a frequency turning point higher than the typical PTA band.This aspect is crucial for elucidating the mass distribution of SMBHBs and investigating potential cosmological components of the SGWB.Additionally, it allows for cross-correlation with other proposed measurements in the same frequency band, such as binary neutron star resonance [78,79].

Identifying Circular Polarization with Astrometry
Astrometry holds the potential to explore circular polarization through EB correlation [21,30,36].We estimate the resolution of the circular polarization fraction v k ≡ V k /I k , considering parameters on each frequency mode O k ≡ {I k , v k }.The corresponding information matrix, as per Eq.(3.7), is given by where Due to the complexity of the expression in Eq. (3.26), we consider a small portion of circular polarization in the limit |v k | ≪ 1.In this limit, the off-diagonal terms I vI k are subleading, implying that I k and v k are uncorrelated.We find the inverse of the uncertainties as (3.27) Here, the uncertainty for I k aligns with the expression in Eq. (3.22), as expected.
Considering a power-law model with α fixed at the fiducial value and assuming a constant circular polarization fraction v k across all frequencies (i.e., v k = v), the uncertainty of v converges to the total SNR: Note that this assumption applies to cosmological sources rather than a finite sum of nearby SMBHBs.
In the top panel of Fig. 4, we display the posterior distribution on the {log 10 A, v} plane for Gaia and XG-Gaia, assuming a truth value of v = 0.1.The degeneracy between the two parameters turns out to be negligible.We also present the marginalized distribution of both parameters, finding that the uncertainty for XG-Gaia, σ(v) ≃ 0.022, is very close to 1/SNR.Thus, XG-Gaia serves as a powerful test for circular polarization that can reach a fractional value as low as 2.2%.In contrast, the marginal detection of SGWB for Gaia (SNR ≃ 1.5) makes it challenging to resolve the circular polarization component.

Synergistic Analyses between PTAs and Astrometry
In this section, we explore the synergistic potential of joint observations involving both PTAs and astrometry, considering all possible correlations.Assuming an ideal scenario where all measurements share the same cadence and total observation time, we consider independent measurements X ℓm,k = {z ℓm,k , E ℓm,k , B ℓm,k } with a covariance matrix Σ ℓ,k defined in Eq. (3.5) are proportional to V k .Inserting Σ ℓ,k into the likelihood Eq. (3.6), the SNR in the limit where where δx,k /N δx ), and we omit the index ℓ and k on the right-hand side for simplicity.We keep ℓ (syn) max to be the one of PTA, as ℓ max for astrometry is typically much higher.
In the weak signal region, where Here, the β −2 , β −1 , and β 0 terms correspond to the estimation from PTA, PTA-astrometry correlation, and astrometry only, respectively.Thus, in addition to the sum of each type of observation, the synergistic observation has the additional contribution from the crosscorrelation among the two.
As astrometry has a higher ℓ max than PTAs, the remaining astrometric-only ℓ-modes can be included separately into the analysis: Here, we define ℓ (syn) max and ℓ (ast) max to distinguish the maximal ℓ from the synergistic analysis and that from the astrometry-only.

Parameter Estimation in Synergistic Analyses of PTA and Astrometry
We explore two pairs of joint observations: the ongoing NANOGrav + Gaia and the future SKA + XG-Gaia.A challenge in these synergistic analyses arises from differences in cadence and total observation time, as detailed in Eq. (3.20) and (3.25).In our approach, we adopt a conservative strategy by selecting the longer cadence, the shorter observation time, and the smaller value among N z and N δx from each pair.The corresponding benchmark parameters are as follows: NANOGrav + Gaia : The measurement uncertainties and the total number of observed stars remain consistent with Eq. (3.20) and (3.25).The total SNR and parameter resolution for the power-law SGWB model are comparable between the more sensitive of the pair, namely NANOGrav and XG-Gaia.Thus, no additional figures similar to Fig. 2 and 3 are presented for further comparison.
The measurement of circular polarization in PTA-astrometry cross-observations can be realized through the zB correlation, in addition to the EB correlation in the astrometry-only correlations.Similar to the process outlined in Sec.3.2.3,we first calculate the information matrix I ij k and derive the corresponding uncertainty of v. Due to the complexity of the expressions, we present the result for the case when |v| ≪ 1: This expression is the linear sum of contributions from both the zB correlation: and the EB correlation in Eq. (3.28).
In the bottom panel of Fig. 4, we present the posterior distribution for {log 10 A, v} for the two joint observations, along with the two astrometry-only observations.The SKA + XG-Gaia pair exhibits slightly better resolution σ(v) ≃ 0.022 compared to the XG-Gaia-only observation, attributed to XG-Gaia's overall better sensitivity.On the other hand, the NANOGrav + Gaia pair can achieve a much superior resolution with σ(v) ≃ 0.39 compared to the Gaia-only observation, benefiting from NANOGrav's higher sensitivity.Thus, we conclude that to effectively constrain circular polarization using current-generation observations, a cross-correlation between PTA and Gaia is necessary.

Emergence of Generalized Hellings-Downs Correlation Patterns
An essential feature for identifying the quadrupolar nature of SGWB lies in the spatial correlation pattern, which can be revealed through either the coefficients of C XX ′ (h)ℓ as a function of ℓ, as discussed in Sec.2.2, or in the separation angle space, such as the Hellings-Downs curve recently explored by various PTA collaborations [5][6][7][8].While we extensively discussed observables in spherical harmonic space in Sec. 3, this section focuses on correlation functions in the separation angle space.
Instead of using the information matrix, we generate random realizations of both SGWBinduced redshift/deflections and measurement noises.These realizations ultimately lead to predictions of spatial correlation patterns with uncertainties in each separation angle bin.

Realization of Correlations in PTA and Astrometry
In this section, we elaborate on the detailed methodology employed for simulating SGWB-induced signals in both PTA and astrometric observations in configuration space, presenting various illustrative examples of results.
To commence, we partition the celestial sphere into N evenly distributed patches using HEALPix [74], denoting their central locations as {n a }.In the frequency domain, the signals are complex variables.Each patch is attributed a stochastic complex dimension-3 vector denoted as (δz a , δx a ), where δx a represents complex dimension-2 vectors on planes perpendicular to each na .The generation of δz a and δx a follows a probability distribution given by Here, C is the 3N × 3N complex covariance matrix, representing the linear sum of the SGWB-induced correlation C (h) and the noise C (n) .The SGWB covariance matrix, following Ref.[19], is defined as: where I and V are the true parameters of total intensity and circular polarization of the SGWB at a given frequency (with the frequency label k omitted for simplicity).The matrices C δz , C δzδx , and C δx are of dimensions N × N , N × 2N , and 2N × 2N , respectively.Their definitions mirror those of Eq. (2.7) and (2.9) but lack the δ-functions: Here, êa || , êb || , and ê⊥ are defined in Eq. (2.8), and the Γ functions are defined in Eqs.(2.7) and (2.10).The symmetry (C ab δx ) † = C ba δx ensures that C δx is Hermitian, and so is the covariance matrix C (h) .
The noise matrix is purely diagonal due to spatially uncorrelated measurement noise: where S (n) z and S (n) δx are defined in Eq. (3.19) and (3.24), respectively.In practice, to generate data from a complex normal distribution, we decompose both the signals and the covariance matrix into real and imaginary parts: Here, ℜ[C] contains signals proportional to I and measurement noise C (n) , while ℑ[C] only contains the term proportional to V .In Fig. 5, we present three cases of realizations for both redshift δz and angular deflection δx with V /I = −1, 0, and 1, respectively.The background circle colors represent the real part of δz, ranging from red (δz > 0) to blue (δz < 0).The real and imaginary components of the angular deflections are depicted with black and white arrows, respectively, with lengths proportional to their magnitudes.In these examples, noise is assumed to vanish.Consequently, in the two maximally polarized cases, the real and imaginary parts of δx are always perpendicular, with the relative phases having different signs for the two cases.On the other hand, for V = 0, they exhibit random behavior without any correlations.
Another consistency check, aiming to connect with the previous discussion in the spherical harmonic space in Sec. 3, involves the reconstruction of the spherical harmonic observable from {z a } and {δx a } using Eqs.(2.12) and (2.13), in the absence of measurement noise.Examples of the corresponding C XX ′ ℓ are depicted in Fig. 6.The violins at each ℓ-mode showcase the statistical distribution of the estimated C XX ′ ℓ values for all ℓ ≤ 6, based on a total of 10 4 realizations.The modes for ℓ = 0 and 1 turn out to be negligible, as expected.We normalize the remaining modes by their expected average values, involving C XX ′ (h)ℓ defined in Eqs.(2.19) and (2.20).Consequently, the red dots representing the average values consistently hover around 1. The variance for I estimators, including zz, EE, BB, and zE correlations, shown in orange, all exhibit uncertainties of approximately 1.This is again consistent with the theoretical predictions in the denominator of Eq. (3.2).For V estimators involving EB and zB correlations, their variance, shown in blue, is larger than 1, attributed to our assumption of V /I = 0.3, and the fact that the variance for V includes a contribution from I.
Figure 6: Reconstructions of rotationally invariant power spectra in the spherical harmonic space for ℓ ranging from 2 to 6.Each reconstruction is derived from a random realization of a map of 108 patches of {δz a } and {δx a } using Eqs.(2.12) and (2.13), in the absence of measurement noise.Total intensity I estimators are depicted in orange, while circular polarization V estimators are shown in blue, assuming V /I = 0.3.Each spectrum is normalized by its theoretical average value.Violins represent the variance of reconstructions from 10 4 realizations, and the red dots denote the average values of realizations.Both are consistent with the variance calculation in the denominator of Eq. (3.2).

Generalized Hellings-Downs Correlations
In this section, we employ the developed simulation to project how spatial correlations in the configuration space manifest for various PTAs and astrometric observations, accounting for a more realistic scenario where both measurement noise and intrinsic variance of SGWB are present.
We consider distinct cases of observations: NANOGrav and SKA as PTAs, Gaia and its next-generation upgrade as astrometry missions, and their cross-correlations.Each observation channel differs in terms of the number of pulsars or patches, and noise levels, with corresponding benchmark parameters listed in Eqs.(3.20), (3.25), and (3.32).For the simulations, we use the fiducial SGWB model defined in Eq. (3.12), assuming V /I = 0 throughout the analysis.
We conduct a total of 100 simulations, each comprising realizations for all frequency bins, resulting in distributions of {δz a , δx a } maps for each frequency.For every pair of patches, characterized by their separation angle θ, we categorize them into 11 evenly distributed bins spanning from 0 to 180 represents the projection of δx a along êa || or ê⊥ .This computation yields a sky-averaged two-point function for a single realization.From the 100 simulations, we determine the standard deviation at the θ i -th separation angle and frequency mode k as σ θ i ,k , incorporating contributions from both the SGWB and measurement noises.
In the left part of Fig. 7, violins are employed to illustrate the distribution from 100 simulations of sky-averaged two-point functions for various observation channels, presenting only the first frequency bin of the observations.The red solid lines represent the theoretical prediction for the average correlations, as defined in Eqs.(2.7) and (2.10).In the right part, we conduct an average across all frequency bins for each separation angle, incorporating a weight factor 1/σ 2 θ i ,k .The frequency-averaged variance σ 2 θ i , defined as 1/σ 2 θ i ≡ Σ k 1/σ 2 θ i ,k , becomes narrower, and the average values, shown in blue, closely align with the theoretical predictions.
Figure 8 illustrates predictions for spatial correlations across various observation channels.Each gray line represents one realization obtained at a specific frequency bin.The opacity of each line corresponds to the weight factor of the frequency, defined as 1/σ 2 k ≡ Σ i 1/σ 2 θ i ,k .The black dashed lines represent the averages of these gray lines, accounting for the weight factor 1/σ 2 θ i ,k across all frequencies.For current observations involving NANOGrav or Gaia, the gray lines fluctuate around the generalized Hellings-Downs curves in red, while the dashed lines exhibit some deviation from it.In contrast, next-generation observations like SKA or XG-Gaia have averaged values that align well with the red.The gray lines are predominantly localized in specific regions, a result of the stochastic nature of the SGWB rather than measurement noise, owing to the high sensitivity of these measurements.In the next subsection, we will delve into a detailed discussion of the origin of uncertainty regions, commonly referred to as cosmic variances.

Cosmic Variances
As depicted in Fig. 8, even in the limit of SNR ≫ 1 for high-sensitivity observations that encompass a large number of patches, the sky-averaged spatial correlations still exhibit an uncertainty envelope.The discussion of these inherent uncertainties for PTAs, referred to as cosmic variance, has been the subject of recent studies in Refs.[80][81][82][83][84][85][86].In this section, we present a theoretical derivation of cosmic variances for general PTA and astrometric observations.
Cosmic variance manifests itself in patch pairs separated by a specific angular separation on the celestial sphere, denoted as the variance intrinsic to the sky-averaged spatial correlations at a given separation angle θ.To calculate cosmic variance, we begin by introducing the sky-averaged two-point functions [82,87]: where X can represent δz, δx || , or δx ⊥ .This expression can be further simplified in the spherical harmonic space, as seen in the case of PTAs [21, 38-41, 82, 88]: Here, P ℓ are the Legendre polynomials.
A parallel set of steps can be applied to astrometric observations, taking into account their expansion in terms of vector spherical harmonics as follows [21]: are functions of associated Legendre polynomials P m ℓ [21].The expression for δx ⊥ δx ⊥ * S θ only differs from Eq. (4.8) by switching G (ℓ1) ↔ G (ℓ2) .The expression for the PTA-astrometry cross-correlation can be derived similarly: The cosmic variance (CV) at a separation angle precisely corresponds to the variance of the sky-averaged two-point correlations defined in Eq. (4.6): which arises from Isserlis' theorem [47].
A similar result can be obtained for the astrometry-only correlations.We begin by calculating the first term in Eq. (4.11), focusing on the parallel directions: where we omit the θ-dependence in G (ℓ1) (θ) and G (ℓ2) (θ) for simplicity.The first two terms within the [• • • ] can be computed directly from Eq. (4.13), while the last two terms represent where we utilized the relationship C EE (h)ℓ = C BB (h)ℓ .The expressions for correlations involving δx ⊥ are derived analogously, with the result given by the above equation and G (ℓ1) ↔ G (ℓ2) swapped.
Note that in Figs.7 and 8, the astrometric observations involve a linear combination of parallel and perpendicular correlations, specifically (δx

Conclusion
In this study, we utilize a joint likelihood that incorporates both astrometric and PTA observations to predict the detection of SGWB and the resolution of SGWB parameters.Our analysis takes advantage of the diagonal structure inherent in both measurement noise and intrinsic SGWB variance in the spherical harmonic space.This results in an analytical framework that facilitates the comparison of various PTA and astrometric observations, providing intuitive insights.Astrometry showcases its advantages in the harmonic space: the abundance of stars allows for a high value of l max , yielding more independent estimators and, consequently, higher SNR in the strong signal regions.This SNR, in turn, translates into the resolution of SGWB parameters, enabling the precise dissection of SGWB properties.
Our findings yield several predictions regarding astrometric observations.Firstly, the upcoming full data release of Gaia is expected to marginally detect the SGWB, which has already been observed by current PTA observations.Consequently, Gaia can serve as a valuable cross-check for PTA results, providing an independent demonstration of spatial correlation distinct from PTA's Hellings-Downs curve.Gaia's individual observations can also contribute to checking the chirality nature of SGWB.However, a significantly improved resolution of chirality can be achieved through cross-correlations between current PTAs and Gaia.The next-generation upgrade of Gaia is poised to deliver the best-ever sensitivity to the SGWB, particularly for its spectrum and chirality.Precise measurements of these quantities are crucial for understanding the evolution of SMBHBs, including eccentricity distribution and environmental effects that may influence the low-frequency end of the spectrum, as well as any potential cosmological signatures.The exploration of the high-frequency turning point in the SGWB spectrum remains an intriguing area, where high-cadence observations by Roman can provide valuable insights.
Looking forward, the prospects for the field appear promising.The anisotropy in the SGWB, a facet beyond the scope of our current analysis, represents another crucial aspect that astrometric observations could potentially illuminate.Our assumption of uniformly distributed stars on the celestial sphere, with uniform noise levels and cadences, can be refined by considering the realistic star distribution found in Gaia datasets.Incorporating this distribution along with the response functions of astrometry and PTA further strengthens the complementarity between these two types of observations, as they exhibit sensitivity to different incoming directions of the GWs.Such angular-dependent sensitivity will play an essential role in resolving individual SMBHBs, likely marking the next significant milestone in gravitational astronomy.

10 AFigure 2 :Figure 3 :
Figure 2:The SNR distribution as a function of the power-law SGWB model parameters, the reference strain log 10 A and spectral index α as defined in Eq.(3.11), for the four considered observations.The left two correspond to PTAs, while the right two represent astrometric observations.The white regions indicate SNR < 1, where the information matrix is not applicable.Yellow stars mark the fiducial parameters from Eq. (3.12), with SNR values of 4.0, 1.5, 34.1, and 43.9 for NANOGrav, Gaia, SKA, and XG-Gaia, respectively.The current exclusion region, defined as SNR > 10 for NANOGrav, is highlighted by yellow dashed lines.

Figure 4 :
Figure4: Posteriors of the reference strain amplitude log 10 A and the fraction of circular polarization v ≡ V k /I k , for astrometric observations (left) and PTA-astrometry synergistic analyses (right).The red stars correspond to the true parameters, including the fiducial value for log 10 A and the assumed v = 0.1.The marginalized distribution of each parameter is presented next to the posterior distribution.The dark and light blue regions represent the 1σ and 2σ regions, respectively.The 1σ uncertainties for v are σ(v) = 0.57, 0.39, 0.022, and 0.022 for Gaia, NANOGrav + Gaia, XG-Gaia, and SKA + XG-Gaia, respectively.

Figure 5 :
Figure 5:Examples of realizations on the celestial sphere for V /I = −1 (top), 0 (middle), and 1 (bottom), featuring a total of 432 patches.Each patch contains information on the real part of the redshift δz represented on a circle, spanning from red (> 0) to blue (< 0).The real (black arrow) and imaginary (white arrow) components, with lengths proportional to their magnitudes, are also displayed.The measurement noise is not included in these examples.In cases with V /I = ±1, the real and imaginary arrows are always perpendicular, while in the unpolarized case, they are uncorrelated.

Figure 7 :
Figure 7: Distribution of sky-averaged two-point functions from 100 simulations presented in gray violins for various PTA and astrometric observation channels, displayed for the first frequency bin (left) and weight-averaged across all frequency bins (right).The red lines illustrate the generalized Hellings-Downs curves defined in Eqs.(2.7) and (2.10).The variance from the simulations encompasses both SGWB variance and measurement noises.The central values of the violins, shown in blue, align with the red lines.

Figure 8 :
Figure 8: Sky-averaged two-point functions from a single simulation for various PTA and astrometric observation channels.Each gray line represents a realization in a specific frequency bin, with opacity inversely proportional to the variance associated with that frequency bin.The black dashed lines depict the weighted averages of all frequency bins, aligning with the generalized Hellings-Downs curves (red) defined in Eqs.(2.7) and (2.10) for observations involving SKA or XG-Gaia.

Finally, turning our
attention to the PTA-astrometry cross-correlation, its CV can be obtained through analogous procedures, yielding: CV(ℜ δz δx || * ) θ = ℓ that the CVs presented in Eqs.(4.12), (4.17), and (4.18) share a common characteristic: each ℓ-mode is the square of the coefficients in their corresponding average values of Eqs.(4.7), (4.8), and (4.10), divided by (2ℓ + 1).The numerical values of these CVs align with the variance envelope observed in the right part of Fig.8.
• .Within each separation angle bin, we compute the averages of the products δz a δz * b , (δx [82]Here, ⟨• • • ⟩ denotes the ensemble average over the SGWB, as the definition of CV does not include measurement noises.The calculation of CV can be simplified using correlations in the spherical harmonic space[82].For the PTA-only observation, the CV becomes:CV(δz δz * ) θ =