A nontrivial footprint of standard cosmology in the future observations of low-frequency gravitational waves

TL;DR

This paper explores how the Hubble constant influences gravitational wave propagation in the universe, proposing a new method to test cosmological models using pulsar timing arrays and local gravitational wave observations.

Contribution

It extends previous models by incorporating the Hubble constant into GW propagation analysis within the LCDM framework, providing a novel observational test.

Findings

Timing residuals show a distinguishable peak depending on the local Hubble constant.

Maximum timing residual correlates with the source's redshift.

Framework can improve gravitational wave measurements with PTAs.

Abstract

Recent research show that the cosmological components of the Universe should influence on the propagation of Gravitational Waves (GWs) and even it has been proposed a new way to measure the cosmological constant using Pulsar Timing Arrays (PTAs). However, these results have considered very particular cases (e.g. a de Sitter Universe or a mixing with non-relativistic matter). In this work we propose an extension of these results, using the Hubble constant as the natural parameter that includes all the cosmological information and studying its effect on the propagation of GWs. Using linearized gravity we considered a mixture of perfect fluids permeating the spacetime and studied the propagation of GWs within the context of the LCDM model. We found from numerical simulations that the timing residual of local pulsars should present a distinguishable peak depending on the local value of the Hubble constant. As a consequence, when assuming the standard LCDM model, our result predicts that the region of maximum timing residual is determined by the redshift of the source. This framework represents a new test for the standard cosmological model, and it can be used to facilitate the measurements of gravitational wave by ongoing PTAs projects.

Figures

Tables

Table 1. Table 1: List of values considered for the parameters in the numerical integration of the timing residual τ GW subscript 𝜏 GW \tau_{\textrm{GW}} in ( 26 ), according to current accuracy of PTAs.

Parameter	SI value
$Z$	$3\times {10}^{24}$ m	$\sim 100$ Mpc
$T$	$Z/c={10}^{16}$s	$\sim 300$ Myr
$L$	${10}^{19}$m	$\sim 1000$ ly
$\Omega$	${10}^{-8}\text{rad}/s$
$\epsilon$	$1.2\times {10}^9$ m

Table 2. Table 2: List of pulsars averaged for an hypothetical source at angular separation α 𝛼 \alpha . It is shown the data given in Manchester et al. ( 2005 ) , where ϕ italic-ϕ \phi is the galactic latitude –transformed to β i subscript 𝛽 𝑖 \beta_{i} – and L i subscript 𝐿 𝑖 L_{i} the distance between Earth and pulsar. We can note that this set simplify the computation of σ 𝜎 \sigma because all pulsars are very close to each other.

Pulsar Name	$\varphi$	$L_i$
J0024-7204E	$-{44.89}^{\circ }$	4.69 kpc
J0024-7204D	$-{44.88}^{\circ }$	4.69 kpc
J0024-7204M	$-{44.89}^{\circ }$	4.69 kpc
J0024-7204G	$-{44.89}^{\circ }$	4.69 kpc
J0024-7204I	$-{44.88}^{\circ }$	4.69 kpc

Equations

$$G_{\mu \nu }+\Lambda g_{\mu \nu }=\kappa T_{\mu \nu },$$

$$g_{\mu \nu }={\eta }_{\mu \nu }+h_{\mu \nu },\text{absolutevalue}{h}_{\mu \nu }\ll 1.$$

$$\square {\bar{h}}_{\mu \nu }=-2\Lambda {\eta }_{\mu \nu }-2\kappa T_{\mu \nu },$$

$${\bar{h}}_{\mu \nu }\equiv h_{\mu \nu }-\frac{1}{2}{\eta }_{\mu \nu }h,$$

$${\partial }_{\beta }{\bar{h}}^{\beta \alpha }=0.$$

$$\square h_{\mu \nu }^{(GW)}=0,\square h_{\mu \nu }^{(\Lambda )}=-2\Lambda {\eta }_{\mu \nu },h_{\mu \nu }^{(\text{bg})}=-2\kappa T_{\mu \nu }.$$

$$h_{\mu \nu }^{\text{(GW)}}=\frac{A_{\mu \nu }}{r}\sin [k_{\mu }{x}^{\mu }]+\text{order}{H}_0^2,$$

$$\text{differential}{s}^2=-\text{differential}{T}^2+a^2(T)\left(\frac{\text{differential}{R}^2}{1-k{R}^2}+R^2\text{differential}{\Omega }^2\right),$$

$${\rho }_i=\{\begin{array}{ll}\frac{4}{3({\omega }_i+1{)}^2\kappa T^2}& \text{if}{\omega }_i\ne -1\\ \Lambda /\kappa & \text{if}{\omega }_i=-1\end{array},$$

$$\text{differential}{s}^2=-\left(1-\frac{\Lambda }{3}{r}^2\right)\text{differential}{t}^2+\frac{\text{differential}{r}^2}{1-\frac{\Lambda }{3}{r}^2}+r^2\text{differential}{\Omega }^2.$$

$$r(T,R)$$

$$t(T,R)$$

$$r(T,R)$$

$$t(T,R)$$

$$\text{differential}{s}^2=-\frac{\text{differential}{t}^2}{\left(1-\frac{\kappa {\rho }_i{r}^2}{3}\right){\left(1+\frac{\kappa {\rho }_i{r}^2(3{\omega }_i+1)}{6}\right)}^{\frac{1-3{\omega }_i}{1+3{\omega }_i}}}+\frac{\text{differential}{r}^2}{1-\frac{\kappa {\rho }_i{r}^2}{3}}+r^2\text{differential}{\Omega }^2,$$

$$t$$

$$r$$

$$t$$

$$r$$

$$a(T)=1+H_0\Delta T+\text{order}{H}_0^2,$$

$$H_0=\sqrt{\frac{\kappa {\rho }_{\text{eff}}(T_0)}{3}}=\sqrt{\frac{\Lambda }{3}+\frac{\kappa {\rho }_{d0}}{3}+\frac{\kappa {\rho }_{r0}}{3}}.$$

$$t$$

$$r$$

$${h^{\prime }}_{\mu \nu }^{\text{(GW)}}=\frac{\left(1+R{H}_0\right)}{R}{A}_{\mu \nu }^{\prime }\sin \left[-w_{\text{eff}}T+k_{\text{eff}}R\right]+\text{order}{H}_0^2,$$

$$w_{\text{eff}}\equiv \Omega \left(1-R{H}_0\right),k_{\text{eff}}\equiv \Omega \left(1-\frac{R}{2}{H}_0\right).$$

$$\varphi (T)={\varphi }_0\left[T-\frac{L}{c}-{\tau }_0(T)-{\tau }_{\text{GW}}(T)\right]$$

$${\tau }_{\text{GW}}(T)=-\frac{1}{2}{\hat{n}}^i{\hat{n}}^j{H}_{ij}(T),$$

$$H_{ij}(T)=\frac{L}{c}{\int }_{-1}^0{h^{\prime }}_{\mu \nu }^{\text{(GW)}}\left(T+\frac{L}{c}x,\text{absolutevalue}\vec{R}(x)\right)\text{differential}x.$$

$${\tau }_{\text{GW}}^{\Lambda \text{CDM}}=-\frac{L\epsilon }{2c}({\sin }^2\alpha {\cos }^2\beta +$$

$$\times {\int }_{-1}^0\frac{1+H_0\left[T_e+\frac{xL}{c}\right]}{Z+xL\cos \alpha }[\cos \Theta +\sin \Theta ]\text{differential}x,$$

$$\Theta =\Omega \left\{\left(1-\frac{Z+xL\cos \alpha }{c}{H}_0\right)\left(T+\frac{xL}{c}\right)-\left(1-\frac{Z+xL\cos \alpha }{2c}{H}_0\right)\left(\frac{Z+xL\cos \alpha }{c}\right)\right\}.$$

$$\sigma =\sqrt{\frac{1}{N_p{N}_t}\sum _{i,j=1}^{N_p,N_t}{\left(\frac{{\tau }_{\text{GW}}}{{\sigma }_t}\right)}^2},$$

$$\sigma (\alpha )=\sqrt{\frac{1}{5\cdot 101}\sum _{i=1}^5\sum _{j=1}^{101}{\left(\frac{{\tau }_{\text{GW}}({\beta }_i)}{{\sigma }_t}\right)}^2}$$

Full text

I Introduction

In the last two decades a lot of astrophysical evidence has been found suggesting that our Universe is expanding at an accelerated rate (Perlmutter et al., 1997; Riess et al., 1998; Planck Collaboration et al., 2018; Riess et al., 2018; Hotokezaka et al., 2019; Wong et al., 2019). This observation is the cornerstone of modern cosmology and represents an ideal setting in which the large scale aspects of gravitation can be tested. A century after the formulation of General Relativity (Einstein, 1916), it remains as the most successful description of gravity, and after many attempts, the Lambda–Cold Dark Matter ($\Lambda$CDM) model has become the standard theoretical framework in which we can study the cosmological phenomena, describing our flat, isotropic, homogeneous and expanding Universe.

Within the $\Lambda$CDM model, the main quantity, among others, that describes the rate of expansion is the Hubble parameter and its value at the present day, the Hubble constant, is denoted by $H_{0}$. The value of this constant represents the current rate of expansion and contains information of the composition of the Universe. Recently there has been great controversy about its actual value, in particular, from the tension in the data obtained from the early (Planck Collaboration et al., 2018) and the late Universe (Riess et al., 2018; Wong et al., 2019). This fact has called the local behavior of $H_{0}$ into question, and therewith many potential explanations to the phenomenon have been raised: From possible new physics to as-yet unrecognized uncertainties from the observations (Odderskov et al., 2014; Ko and Tang, 2016; Freedman, 2017; Bringmann et al., 2018; Camarena and Marra, 2018; Mörtsell and Dhawan, 2018; Di Valentino et al., 2018; Feeney et al., 2019). Nevertheless, no consensus has been reached so far and an intense debate is still going on, waiting for improved empirical data or more robust and complete theories of gravitation.

On the other hand, the measurement of Gravitational Waves (GWs) carried out by the ground-based detector LIGO (Abbott et al., 2016), is considered as one of the last experimental verifications of general relativity. Actually, GWs have been used as a new alternative for the measurement of cosmological parameters: A gravitational–wave standard siren was used to measure the Hubble constant independently (LIGO Collaboration et al., 2017). This kind of observations show us that the efforts involved in the accurate measurement of gravitational waves could be useful in the analysis of cosmological parameters within the next years. Nevertheless, ground-based detectors are not the unique way to detect GWs. Another type of ongoing gravitational-wave experiment is the Pulsar Timing Array (PTA), which uses the residual time of the arriving electromagnetic emissions from different millisecond-pulsars located in the Milky Way and their correlations to determine the presence of gravitational radiation (Barke et al., 2015; Hobbs and Dai, 2017; Cordes et al., 2019; Burke-Spolaor et al., 2019). Many projects have been developed, as NANOGrav (McLaughlin, 2013), Parkes (Hobbs, 2013) or the European PTA (Ferdman et al., 2010), members of the International Pulsar Timing Array collaboration (Hobbs et al., 2010) (IPTA), which together, seek to identify and measure low-frequency (i.e. $\sim 10^{-9}$ to $10^{-8}$ Hz) gravitational waves coming from astrophysical sources as an isotropic and stochastic background (Lentati et al., 2015; Verbiest et al., 2016; Arzoumanian et al., 2018) or Continuous Gravitational Waves (CGWs) from individual Supermassive Binary Black Holes (SMBBHs) (Babak et al., 2016; Mingarelli et al., 2017; Kelley et al., 2019), in which PTA experiments can detect timing differences of $\sim 100$ ns and the dimensionless strain amplitude of CGWs is expected to be $h\sim 10^{-15}$ for sources at redshift $z\lesssim 0.5$. The second data release of the IPTA collaboration was published recently Perera et al. (2019), but the results are not conclusive yet

Considering the previous facts, the work of Bernabeu et al. (2011) shows the effects produced by a non–zero cosmological constant $\Lambda$ on the propagation of gravitational waves, which includes corrections of order $\sim\sqrt{\Lambda}$ to the phase and the amplitude. Later, it was found Espriu and Puigdomènech (2013); Espriu (2014) that certain modifications to the frequency (i.e. the usual redshift) and a non–trivial correction to the wave number can be found. Furthermore, in those papers was proposed that the magnitude of residual time in a PTA experiment could change due to the action of the cosmological constant. A wide explanation and a review of this phenomenon can be found in the work of Alfaro et al. (2019), where the action of non-relativistic matter was included into the phenomenon, showing that Dark Matter increases the effect of $\Lambda$ on the propagation of gravitational waves. In that sense, we expect that each cosmological component of the Universe affects the propagation of GWs similar to the case of $\Lambda$, reason of why we will follow the idea of using the Hubble constant as the main control parameter.

The manuscript is structured as follows: In Section II we develop the theoretical formalism of our framework, using the linearization of General Relativity and the $\Lambda$CDM model to show how the Hubble constant affects the propagation of Gravitational Waves. The entire Section III is devoted to analyze the possibility of using PTA experiments to observe the studied phenomenon: In subsection III.1 we shortly explain the working of PTA and in subsection III.2 we show how our framework can be tested using the previous setup. Furthermore, in subsection IV.3 we present a reasonable approximation that give us a relationship between the location of the maximum timing residual and the redshift of the source. Finally, in Section V we present our conclusions and an outlook for future research.

II Gravitational Waves in an expanding Universe

II.1 Linearized Theory of General Relativity and Standard Cosmology

The classical theory of General Relativity, finally consolidated by Einstein (1916), predicts that matter and the curvature of spacetime are related through the Field Equations,

[TABLE]

where $\kappa=8\pi G/c^{4}$, $G_{\mu\nu}$ are the components of the Einstein tensor, $T_{\mu\nu}$ are the components of the stress-energy tensor and $\Lambda$ is the cosmological constant, which plays an important role within the context of cosmology since it can be interpreted as the energy related to the vacuum, called commonly Dark Energy. The prediction of the existence of gravitational radiation (Einstein, 1916), came from the linearization of (1), which is nothing else than doing perturbation theory around a flat spacetime,

[TABLE]

The computation of the linearized Field Equations can be found in any text of General Relativity (Cheng, 2010), giving the following expression,

[TABLE]

where $\bar{h}_{\mu\nu}$ is the trace-reversed perturbation defined by

[TABLE]

which satisfies the Lorenz gauge condition

[TABLE]

Note that there is still gauge freedom, but we can completely fix it using the Transverse–Traceless gauge, in which $\bar{h}_{\mu\nu}=h_{\mu\nu}$. Therefore, the metric perturbation can be decomposed (Alfaro et al., 2019) as $h_{\mu\nu}=h_{\mu\nu}^{(\textrm{GW})}+h_{\mu\nu}^{(\Lambda)}+h_{\mu\nu}^{(\textrm{bg})}$, which satisfy the following field equations,

[TABLE]

The linearization is well justified because the perturbations from the Minkowski spacetime are very small when we are far away from the source, which will be our case of study. As explained by Bernabeu et al. (2011) and Espriu (2014), by using a time–independent –and at order $\sqrt{\Lambda}$– coordinate transformation, it can be shown that the solution of $h_{\mu\nu}^{(\Lambda)}$ corresponds to the linearization of a static spherically symmetric metric (i.e. the Schwarzschild-de Sitter metric) in the Lorenz gauge. In this work the relevant corrections will be of order $H_{0}$ (i.e. $\sim\sqrt{\Lambda}$), and therefore, the analysis of GWs and the decomposition (6) will be simpler in the Lorenz gauge. As the contribution of different terms turns to be additive, the discussion seems straightforward.

If the source is located at a large (but non cosmological) distance, the perturbation will be described approximately by a combination of harmonic functions (Espriu, 2014; Alfaro et al., 2019)

[TABLE]

where $A_{\mu\nu}$ are the components of the polarization tensor, $k^{\mu}$ is the 4–wavevector and $x^{\mu}$ is the 4–position. We will denote $\Omega$ as the angular frequency of a monochromatic gravitational wave. It is important to note that the coordinates $\{t,r\}$ corresponds to the spherically symmetric context explained before and they will represent the coordinates which origin is placed at the –usually spherical shaped– source of GWs.

On the other hand, it is commonly accepted that we live in a mostly isotropic and homogeneous Universe. As a first approximation, we can use the FLRW metric to describe the geometry of spacetime (Cervantes-Cota and Smoot, 2011), which is of the form

[TABLE]

where $a(T)$ is the adimensional scale factor and the spacetime is described by comoving coordinates $\{T,R\}$. In this work we will use $k=0$, which represents a globally flat geometry of the Universe, as it is currently observed (Planck Collaboration et al., 2018). If we consider different perfect fluids as material components of the Universe, and label them with the subscript $i$, having each of them an energy density $\rho_{i}$, isotropic pressure $p_{i}$ and a equation of state $p_{i}=\omega_{i}\rho_{i}$, the Friedmann equations will give us the following expression for the energy density of the i–th fluid,

[TABLE]

where $T$ is the time coordinate in the FLRW metric. An extensive treatment of this computation can be found in the Appendix A.

II.2 Coordinate transformation and an application to the $\Lambda$CDM model

Let us consider the following situation: Two remote galaxies are merging so their central super-massive black holes are orbiting around a common center of mass and slowly approaching to each other. This is basically a Keplerian problem with –approximate– spherical symmetry. For this reason, near the source of gravitational waves, the set of coordinates $\{t,r,\theta,\phi\}$ with spherical symmetry is very useful in order to describe spacetime and their perturbations. A very detailed discussion about these considerations are given in Espriu (2014), where also it is established that GWs in their simplest form and expressed in these source-centered coordinates will have the form of the equation (7).

However, these coordinates are not useful in cosmology because the cosmological measurements are described in comoving coordinates, i.e. $\{T,R,\theta,\phi\}$. Thus, the main objective of this research is to find the coordinate transformation between $\{t,r\}$ and $\{T,R\}$ and thus, taking advantage of the principle of covariance and relativity, we will able to describe the propagation of gravitational waves as seen by a cosmological observer. For a compendium with all the results found before the writing of this paper, including the coordinate transformations of a de Sitter Universe, besides other implications, see Alfaro et al. (2019).

The easiest example that we can give to illustrate the situation is showing the de Sitter case, where only the action of the cosmological constant is taken into account. We note that in a vacuum background, an approximately spherical source of GWs would produce a Schwarzschild metric. Thus, if we take $\Lambda\neq 0$, then when we are far from the source (i.e. neglecting the mass term at the cosmological horizon, where $\Lambda r^{3}\gg 6M$), the geometry of spacetime will be described approximately by a de Sitter (dS) metric,

[TABLE]

On the other hand, the FLRW metric –expressed in comoving coordinates– will be given by (8) and the scale factor is of the form $a(T)=a_{0}\exp(\sqrt{\Lambda/3}\Delta T)$ (see appendix A). The main idea is to express the coordinates of the SdS metric in terms of the comoving coordinates of the FLRW metric, as they are two equivalent representations of the same spacetime. It was found that the coordinate transformation is given by,

[TABLE]

If we expand them at order $\sqrt{\Lambda}$, we get,

[TABLE]

This result was obtained previously Bernabeu et al. (2011); Espriu and Puigdomènech (2013), and it was the starting point of this line of work. From these transformations it was found how the cosmological constant affects on the propagation of GWs and it was explained how this effect could be measured using PTAs Espriu (2014). These expressions show us how comoving coordinates are related to the coordinates in the dS metric, and replacing the coordinates in (7) would show how GWs are seen by a cosmological observer. That is the main idea and it is what we are going to exploit next.

In order to develop a more general discussion of the phenomenon, we will first consider a Universe filled by a single fluid with an arbitrary equation of state, i.e. $p_{i}=\omega_{i}\rho_{i}$. The methodology to be used is basically build a diagonal, spherically symmetric and asymptotically flat metric (that we will denote by SS$\omega_{i}$), described in the coordinates $\{t,r\}$, that recovers the corresponding FLRW metric in comoving coordinates, then find the coordinate transformation between both frames and, finally, replace the coordinates in (7), showing how it affects on the propagation of gravitational radiation.

In the appendix B is available the full derivation of the exact expression of the the SS$\omega_{i}$ metric, which is of the form

[TABLE]

and the coordinate transformations between $\{t,r\}$ and $\{T,R\}$ in terms of $\rho_{i}$ and $\rho_{0}=\rho_{i}(T_{0})$, which are

[TABLE]

where $A$ and $n$ are constants (see appendix B). The expansion of (16a) and (16b) in terms of the energy density of the fluid at the present day, i.e. $\rho_{0}$, becomes

[TABLE]

These results agree with the previous works Espriu and Puigdomènech (2013); Espriu (2014); Alfaro et al. (2019), and they show that, regardless the equation of state, the first term in the expansion is always at order $\sqrt{\rho_{0}}$. Using this fact, we can expand the First Friedmann equation at first order in $H_{0}$,

[TABLE]

where $H_{0}$ is the Hubble constant, given by

[TABLE]

As the $r$ coordinate has to transform as $r\to a(T)R$ to preserve spherical symmetry, a comparison between (17a), (17b) and (18) shows that the cosmological components are added inside the square root, as was discussed by Alfaro et al. (2019) in the case of non-relativistic matter. Thus, in order to obtain the correct limits for the previous models, the most general linearized coordinate transformations must be of the form,

[TABLE]

Note that these are linear in $H_{0}$, so only small effects will be considered. By replacing (20a) and (20b) into (7), we obtain an expression in terms of the comoving coordinates,

[TABLE]

where $A_{\mu\nu}^{\prime}$ are the transformed components of the polarization tensor, and the effective angular frequency and wave number are given by

[TABLE]

From the last two expressions we can infer how the Hubble constant affects the propagation of GWs when a cosmological observer (e.g. laboratories in the surface of the Earth or local celestial bodies as pulsars) is measuring them, using comoving coordinates. These results show that the previous findings Bernabeu et al. (2011); Espriu and Puigdomènech (2013); Espriu (2014); Alfaro et al. (2019) were merely approximated and, therefore, incomplete. At the same time, we discard the idea of the possibility of measuring the cosmological constant $\Lambda$ separately from the other components of the Universe, since all of them are coupled within $H_{0}$. It is important to note that the expression for the effective frequency $w_{\textrm{eff}}$ in the equation (22) reproduces the usual cosmological redshift expected from the expansion of the Universe. However, the effect on the wave number $k_{\textrm{eff}}$ cannot be derived from other simpler considerations, e.g. time dilation or redshift, and represents an additional feature of this framework. On the other hand, as it is discussed by Alfaro et al. (2019), even when the phase velocity of the GW is not exactly $1$ (in natural units), if it is computed with respect to the ruler distance traveled, it can be shown that its value is exactly equal to $1$.

III Pulsar Timing Arrays

III.1 Timing residual and the working of PTAs

The results obtained in the last section, e.g. equations (21) and (22), show that the Hubble constant should influence on the propagation of Gravitational Waves. Now we will set an experimental framework in which this effect can be eventually measured. For this, we will use the light coming from a local pulsar and the shift in the time of arrival of the electromagnetic (EM) pulse due to the pass of GWs. In the following picture we show the simplest configuration, which will guide our discussion.

From figure 1, we note that the Earth and the pulsar are gravitationally bounded to the Milky Way, so they do not feel the Universe’s expansion. However, the source of GWs and the system Earth-pulsar are not bounded, so they do feel accelerated expansion and, therefore, the discussion and results of the previous chapter apply. The pulsar emits light with a particular EM field. Denoting the time-dependent phase of this field at the pulsar as $\phi_{0}$, then the phase of the EM pulse measured from Earth can be expressed as

[TABLE]

where $c$ the speed of light, $\tau_{0}(T)$ is the timing correction associated to the motion of the Earth respect to the Solar system and $\tau_{\textrm{GW}}(T)$ is the timing correction due to the action of GWs passing through the system. The correction due to the action of Gravitational Waves is given (Finn, 2009; Deng and Finn, 2011; Espriu, 2014; Alfaro et al., 2019) by the following expression

[TABLE]

where $\hat{n}=(-\sin\alpha\cos\beta,-\sin\alpha\cos\beta,\cos\alpha)$ is a unit vector pointing from Earth to the pulsar and $H_{ij}$ is the integral of the metric perturbation along the null geodesic in the path pulsar–Earth, which could be parameterized by $\vec{R}(x)=\vec{P}+L(1+x)\hat{n}$ with $x\in[-1,0]$. Using this path, $H_{ij}(T)$ takes the following form

[TABLE]

III.2 Including the $\Lambda$CDM model

In our framework (see figure 1), the source of GWs is far away from Earth, though the pulsar is at a local distance. Therefore, we can consider the reasonable approximation $L/Z\ll 1$ and using it we can show that $R(x)\approx Z+xL\cos\alpha$. Now we compute the equation (24) using (25). However, in the TT-Lorenz gauge, for a GW propagating through the Z axis, the only non-zero values of the components of $A^{\prime}_{\mu\nu}$ are in the X,Y components (Bernabeu et al., 2011). We can also additionally assume, in order to simplify the computation of the integral, that $\absolutevalue{A^{\prime}_{\mu\nu}}\equiv\varepsilon\;\forall\mu,\nu$. Thus, the full timing residual in the arrival time of the pulsar due to the pass of GWs in the $\Lambda$CDM model becomes,

[TABLE]

where

[TABLE]

IV Results and discussion

IV.1 Simulation of the timing residual of an individual pulsar

From a geometric argument, we can always fix the pulsar, the Earth and the source of gravitational waves in the same plane, thus we can set $\beta=0$, and the geometrical parameters involved in $\tau_{\textrm{GW}}$ are the angle $\alpha$ between the Pulsar and the GW-Source, the distance Earth–Pulsar, $L$, and the distance Earth–GW-source, $Z$. In order to perform a numerical analysis we can choose some reasonable values of the parameters that appear in the equation (26) and fix them to visualize the behavior of the timing residual $\tau_{\textrm{GW}}$. Thus, the setup described in the figure 1 can be approximately modeled with the values that appear in the Table 1.

For the source of GWs, we choose a typical distance $Z$ where supermassive black holes are present and although is large, it is not a cosmological distance near to the Big Bang. On the other hand, the distance between Earth and the pulsar is within the margin of a local galactic scale. It can be seen, from table 1, that $L\ll Z$, as required from the previous considerations. The angular frequency is of the expected order for future PTA projects and the same argument is used to fix $\varepsilon$, due to that it satisfies $\absolutevalue{h}\sim\varepsilon/R\sim 10^{-15}$, where $\absolutevalue{h}$ and $\Omega$ are within the expected accuracy of PTA projects, e.g. the EPTA (Babak et al., 2016) or the NANOGrav collaboration (Kelley et al., 2019). Employing these parameters, the numerical integration of $\tau_{\textrm{GW}}$ gives the results that are shown in figure 2a. As we can see, the value of the timing residual can be positive or negative. Since the meaningful physical magnitude is the amount of time, rather than the direction of the shift, we can also plot the absolute value of $\tau_{\textrm{GW}}$, but now changing the value of $H_{0}$ within a region of parameters, obtaining the plot in figure 2b.

IV.2 Statistical significance and the timing residual analysis

The most important feature of these figures is the presence of a considerable peak in the value of the timing residual for a certain angle $\alpha$. Moreover, we note that this peak changes its angular position with the value of the Hubble constant. This is our first clue of the existence of a distinguishable signal coming from the cosmological effects on the propagation of GWs.

In order to make an analysis of the possible signal shown in the previous figures, we will use some pulsars from the ATNF catalog (Manchester et al., 2005) and follow the ad hoc analysis of Espriu (2014). As we know, pulsars are stable clocks whose periods are known with great accuracy. Assuming a modest precision of $\sigma_{t}=9.6\times 10^{-7}s\approx 10^{-6}s$ which is obtained by averaging the precision achieved of best pulsars in the IPTA collaboration, we can define a statistical significance of the timing residual, of the form

[TABLE]

where index $i$ running from 1 to $N_{p}$ (number of pulsars averaged) and $j$ running from 1 to $N_{t}$ (number of observations). Assuming we perform measurements every 11 days through 3 years, then $N_{t}=101$. The pulsars belong to the considered cluster are shown in table 2.

We will keep $\alpha$ as a free parameter and suppose that an hypothetical GW source is located at $\alpha$ radians between Earth and pulsars. Thus, the statistical significance is given by

[TABLE]

The result of the simulation can be observed in the Figure 3a, showing the characteristic peak as we expected. However, we can develop a more realistic simulation. In figure 3a, only a cluster of 5 pulsars were considered and all of them were averaged at the same angle $\alpha$. However, one can expect that all the pulsars are located at different angles (in galactic coordinates) and basically being randomly located. Therefore, we have considered 11 randomly distributed groups of 5 pulsars each (see appendix D), two test clusters of pulsars with a suitable location (65 pulsars in total) and a source of gravitational waves located at galactic coordinates $\theta_{S}=20^{\circ}$ and $\phi_{S}=15^{\circ}$.

Then, we averaged them using the statistical significance given by

[TABLE]

and plot it as a function of the average angle of the group, namely $\bar{\sigma}_{k}=\sum_{i=1}^{5_{k}}\alpha_{i}/5$. From this simulation we obtained figure 3b, where we can note how the randomly distributed pulsars mostly do not show any signal, except for those that are located very close to the maximal angle, namely $\alpha_{m}$, where the value of the timing residual is maximized. For a pulsar located within the vicinity of the angle $\alpha_{m}$, the effect of the Hubble constant on the propagation of GWs and the capability of measure them using PTAs is greatly increased.

This fact indicates that only the pulsars placed near the angle $\alpha_{m}$ with respect to the source of GWs will show the characteristic peak in the timing residual with great statistical significance, which implies a major obstacle when trying to observe this effect. Nevertheless, as more pulsars are observed and studied, it is more likely to measure the existence of this peak, which could represent a challenge for the future astrophysical research of PTA experiments.

First, in order to understand the role of the angle $\alpha_{m}$, in which the maximum of $\tau_{\textrm{GW}}$ is reached, we analyze the dependency on the original angular frequency of the incoming gravitational wave, namely $\Omega$, e.g. see eq. (27). From figure 4a we note that for the region $10^{-6}\;\textrm{rad/s}<\Omega<10^{2}\;\textrm{rad/s}$, the value of $\absolutevalue{\tau_{\textrm{GW}}}$ is practically zero. However, in the region $10^{-8}\;\textrm{rad/s}<\Omega<10^{-6}\;\textrm{rad/s}$ the value starts to rise. This is the main reason why the other type of detectors as LISA or LIGO are useless in this context: Only PTA works in the proper range of the frequency spectrum (McLaughlin, 2013; Barke et al., 2015; Verbiest et al., 2016).

In figure 4b we note a lack of angular dependence on the maximum values of $\tau_{\textrm{GW}}$, and moreover, the same behavior is observed for the distance $L$, in figure 5a.

However, in figure 5b it is clear that for different values of the distance Earth–GW-Source, the angle $\alpha_{m}$ changes dramatically and therefore it depends on the distance $Z$ in an explicit but unknown way. For the Hubble constant, figure 6a shows a similar situation. In fact, it represents how the dependency on $H_{0}$ is actually quite similar to the case of $Z$ in figure 5b. From this behavior we can speculate about a possible relationship between these parameters. Zooming up, we obtain the figure 6b, where we can see that the maximum value of the timing residual $\tau_{\textrm{GW}}$ (i.e. the white spots) has a slightly oscillatory structure around a characteristic maximal angle $\alpha_{m}$. It can be noticed that the full white strip has a very small slope, showing a slow variation of the angular position, in agreement with the figure 2b.

We can summarize the previous analysis as follows: The numerical integration of the equation (26), with the parameters given by the Table 1, gives us the density plots shown above. From figure 4a we can establish the crucial role of PTAs in the eventual measurement of the effects caused by the cosmological components of the Universe on the propagation of gravitational waves: The kind of effect that Dark Matter, Dark Energy and others cosmological fluids induce to the propagation of GWs is sensitive to the frequency spectrum of current Pulsar Timing Arrays experiments.

Furthermore, the figures 5b and 6a support the hypothesis of the existence of an implicit relationship between the $H_{0}$, $Z$, and $\alpha_{m}$. From the figures 2a and 2b we note that it could be difficult to find the best value of $H_{0}$ that fit the data unequivocally, due to the complicated dependencies of the parameters within the equation (26) and because of the background (i.e. the null test for $H_{0}$, the flat Minkowski spacetime) has a not negligible amplitude for a large part of the angular values of $\alpha$. However, this situation does not happen within the vicinity of the characteristic peak of $\tau_{\textrm{GW}}$ when $H_{0}\neq 0$, which occurs at a certain angle $\alpha_{m}$, where the difference with the background is considerable. In that sense, if this tiny effect is ever observed, it is very likely that must be the case of a pulsar located in the vicinity of the angle $\alpha_{m}$, in such a way that the value of its timing residual is maximized, and that the signal is strong enough to rule out the experimental noise. In the section IV.4 the relationship between $H_{0}$, $\alpha_{m}$ and $Z$ will be clarified, and we could state the importance of it.

IV.3 A relationship between PTA observables and the Hubble constant

Using the stationary phase approximation for $\tau_{\textrm{GW}}$ and considering reasonable asymptotic expansions, we obtain the following expression,

[TABLE]

A derivation of this equation can be found in appendix C. The expression in (31) provides a precise relationship between $H_{0}$ and the observables $\alpha_{m}$ and $Z$. Under ideal assumptions, this formula can be used to estimate the local value of the Hubble constant knowing the two main observables: $Z$ and $\alpha_{m}$. In figure 7 we show the behavior of the approximation formula with respect to the numerical analysis, included some values of current observations.

However, there are some experimental obstacles when we want to measure $H_{0}$ using this framework, in particular because we need to determine the value of $Z$ by independent astrophysical methods. Nevertheless, a PTA experiment could determine, considering the experimental uncertainties, the value of the timing residual $\absolutevalue{\tau_{\textrm{GW}}}$, the value of the frequency of the source $\omega_{\textrm{eff}}$ and the direction of the incoming Gravitational Wave (and therefore, the value of $\alpha$). In the ongoing PTA experiments, several pulsars are observed at once. Each of them will be located at different angular positions. If a GW is passing through the set of pulsars, then each of them will have, according to our model, a different associated $\absolutevalue{\tau_{\textrm{GW}}}$. Evidently, it is more likely that the pulsars have a relative angular position with the source in which the $\absolutevalue{\tau_{\textrm{GW}}}$ is not maximum (i.e. $\alpha\neq\alpha_{m}$). Clearly, these cases are very difficult to measure experimentally, because the value of $\absolutevalue{\tau_{\textrm{GW}}}$ for these cases is indistinguishable from the background. However, there is a tiny possibility that one or some of these pulsars are located just at an angular position where the value of $\absolutevalue{\tau_{\textrm{GW}}}$ is clearly distinguishable with respect to the background (i.e. $\alpha\approx\alpha_{m}$). If that were the case, then the value of $\absolutevalue{\tau_{\textrm{GW}}}$ should be large enough to be measured by the expected accuracy for PTA experiments.

IV.4 A test of the $\Lambda$CDM model using PTAs

With the spirit of testing the standard cosmological model, we know that for a small redshift, i.e. $z\ll 1$, it holds that $z\approx(Z/c)H_{0}$ Ryden (2003). Thus, when we compare this known formula with (31), we obtain an expression that could be interpreted as the redshift of the source of GWs in terms of the maximal angle $\alpha_{m}$,

[TABLE]

This relationship between the redshift of the source and $\alpha_{m}$ had not been identified in any of the previous works of Espriu and Puigdomènech (2013); Espriu (2014) and Alfaro et al. (2019), hence it represents a new and interesting finding in this line. This result implies that the $\Lambda$CDM model predicts large timing residuals for pulsars located approximately at the same direction of the (local) sources of gravitational waves. This effect could be useful for PTAs research and the measurement of Continuous Gravitational Waves through this method, a task that is expected to be accomplished in the next years (Babak et al., 2016; Mingarelli et al., 2017; Kelley et al., 2019): The angular separation between the sources of GWs and the monitored pulsars is strongly constrained in order to measure distinguishable differences in the timing residual. Therefore, if a PTA experiment seeks to measure a strong signal of gravitational radiation coming from single sources (e.g. mergers of supermassive black holes), the monitored pulsars must be located near to the angle $\alpha_{m}$ with respect to the sources. If that is not the case, the measurement of a powerful signal is very unlikely. In that sense, this prediction of the standard $\Lambda$CDM model represents an alternative way to improve the detection of GWs using PTA experiments by taking advantage of this local cosmological effect on the propagation of gravitational waves within our expanding Universe.

V Conclusions

In this manuscript we have extended the line of research developed in Bernabeu et al. (2011); Espriu and Puigdomènech (2013); Espriu (2014); Alfaro et al. (2019) on the action of cosmological parameters in the propagation of gravitational waves. This effect, an additional and different contribution besides the redshift in the frequency, is completely caused by the coordinate transformation between a frame centered onto the source of GWs and the comoving frame of the FLRW metric. The main result is that harmonic gravitational waves become anharmonic when a cosmological observer measures them, i.e. the Earth or PTAs. In particular, we have generalized the previous attempts to a case where all the cosmological contributions (e.g. Dark Energy, non-relativistic matter and radiation) are naturally included by the Hubble constant, in contrast to the previously developed idea of analyzing each of them separately, which seems to be wrong at first order in the perturbations.

We have shown that all the coordinate transformations involved are linear in terms of $H_{0}$, which turned out to be the parameter that governs the behavior of local cosmological fluctuations in the propagation of gravitational waves. From the expression (31) an unexplored possibility of study the local behavior of $H_{0}$ opens, and besides, from (32) we shown a direct relationship between the redshift of the source (for $z\ll 1$) and $\alpha_{m}$, in which the standard $\Lambda$CDM model predicts the observation of a strong signal in the time of arrival of electromagnetic emission of a monitored pulsar, when it is located at an angle close to $\alpha_{m}$ with respect to the source of continuous gravitational waves.

This local cosmological action appears to be very sensitive to PTA observations (basically due to the frequency spectrum of gravitational waves that PTAs could see). Although the process of measurement could be quite tough due to the experimental difficulties in the measurement of the astrophysical parameters involved (e.g. distances, redshift and timing residuals), some very privileged millisecond pulsars can be analyzed with an incredible accuracy (Manchester, 2017) and they could be useful in order to test the $\Lambda$CDM model within the near future. For this reason we believe it is worthwhile to perform a more detailed study of this phenomenon, for example, considering different sources or not–monochromatic GWs, in order to have a more realistic astrophysical analysis.

Acknowledgments

The authors thank D. Espriu and L. Gabbanelli for many interesting conversations. M. Gamonal and J. Alfaro are partially supported by Fondecyt 1150390 and CONICYT-PIA-ACT14177.

Appendix A On the equations of standard Cosmology

The components of the stress-energy tensor that appears in (3) can be found by considering a perfect fluid (i.e. a fluid that does not have viscosity and does not conduce heat), with energy density $\rho$ and isotropic pressure $p$, filling the whole Universe. The components of $T^{\mu\nu}$ take the following form,

[TABLE]

where $U^{\mu}$ are the components of the 4-velocity of the fluid. When this expression is inserted into the Einstein’s Field equations we obtain the two Friedmann equations: The first equation obtained from the 00 component of (1) and the second from the combination between the trace of the field equations and the first Friedmann equation, giving,

[TABLE]

where $\rho_{i}$ and $p_{i}$ are the energy density and the isotropic pressure of the i-th fluid respectively and $H(T)$ is the Hubble parameter, which value at the present day, i.e. $H(T_{0})\equiv H_{0}$, is known as the Hubble constant $H_{0}$. In order to obtain the time evolution of the scale factor an equation of state must be provided. In Alfaro et al. (2019) was used $p=0$, which corresponds to the equation of state of non-relativistic dust. However, in this work we will use $p_{i}=\omega_{i}\rho_{i}$, with $\omega_{i}$ constant, in order to develop a general discussion of the phenomenon. Using the Friedmann equations we can find that

[TABLE]

where $\rho_{0}=\rho(T_{0})$ is the current energy density of the i-th fluid and $a_{0}=a(T_{0})$ is the current scale factor (usually taken as 1), and both are integration constants. Replacing this expression into (34) provides a solution of the scale factor in terms of the comoving time and the equation of state,

[TABLE]

In the $\omega_{i}\neq-1$ case, when we combine (36) with (37) we can obtain the general form of the energy density of the i–th fluid,

[TABLE]

Currently, the standard cosmological model is the $\Lambda$CDM: Includes a positive cosmological constant $\Lambda$ (which represents the so-called Dark Energy) and Cold Dark Matter (which is the union of baryonic matter and non-relativistic dark matter). In the $\Lambda$CDM model, when a global flat geometry is considered, we can use (36) to write an effective energy density in terms of the scale factor and the currently evaluated energy densities,

[TABLE]

where $\rho_{\Lambda}=\Lambda/\kappa$, $\rho_{d0}$ is the current density of non-relativistic matter (i.e. Cold Dark Matter and baryonic matter, $\omega_{d}=0$) and $\rho_{r0}$ is the current radiation density ($\omega_{r}=1/3$). These expressions will be used in order to construct a spherically symmetric metric which reproduces the corresponding geometry of a FLRW metric for a perfect fluid with an arbitrary $\omega_{i}$.

Appendix B On the derivation of the SS$\omega_{\scriptsize{i}}$ metric

As we have to impose a spherically symmetric geometry we will have the transformation $r^{2}\,\differential\Omega^{2}\to a(T)^{2}R^{2}\,\differential\Omega^{2}$. Using the second rank tensor property of the metric tensor when we perform coordinate transformations,

[TABLE]

and the requirement that the new metric must be diagonal, we obtain the relation

[TABLE]

By computing the partial derivatives we obtain the expressions

[TABLE]

and from (41) we find that

[TABLE]

Thus, from the last equation, $\partialderivative{T}{r}$ becomes

[TABLE]

and using (40) we can obtain the components of the metric,

[TABLE]

From (36) we can write the SS$\omega_{i}$ metric as

[TABLE]

but, using (36) and (45), we get

[TABLE]

If we properly redefine $\tilde{\rho}_{i}\equiv\kappa\rho_{i}$, the last expression becomes

[TABLE]

but it can be noticed from (50) that we can form the expression

[TABLE]

where $c$ and $n$ are unknown constants that we suppose exist. Unfolding the last expression and using the linear independence of $r$, we obtain that the constants are

[TABLE]

Therefore, we can integrate (51) and write

[TABLE]

where $F(t)$ is a function of $t$. By a dimensional analysis, we note that in natural units $[\tilde{\rho}_{i}]=L^{-2}$ and therefore $[F(t)]=L^{2n}$. As there is no other parameter involved apart from $t$, and also as $[t]=L$ in natural units, then we set $F(t)=At^{2n}$, with $A$ as a dimensionless arbitrary constant. For any fluid we can expect that at later stage it will be diluted homogeneously, which implies that for $t\to\infty$ the metric (48) is almost flat. Then,

[TABLE]

On the other hand, (53) can be written as

[TABLE]

but when we take the derivative with respect to $t$ and solving for $\partial_{t}\rho_{i}$, we obtain

[TABLE]

and if we square, divide by $3\kappa\rho_{i}^{3}$ and replace the previous results, we can found the following equality

[TABLE]

Computing the limit $\rho_{i}\to 0$ as the fluid dilutes at distant times, we can set $A$,

[TABLE]

and using that the metric is asymptotically flat, which implies that the previous limit is equal to one, we get the value of $A$,

[TABLE]

Finally, with the constant $A$ known, we can provide an exact expression for the the SS$\omega_{i}$ metric, which becomes

[TABLE]

and from (55) we can express the coordinate transformation between the SS$\omega_{i}$ and the FLRW frames in terms of $\rho_{i}$ y $\rho_{0}$,

[TABLE]

Appendix C On the accuracy in the approximation of $H_{0}$

In order to simplify the computation, we can omit the geometrical prefactor that appears in (26), because it is common to every observation and is $H_{0}$–independent. Therefore, we define a reduced timing residual,

[TABLE]

with $\absolutevalue{R_{1}}\leq\frac{LH_{0}}{cZ}\sim 10^{-31}$ s. Then, we take the reduced timing residual from (63) and note that $R_{1}$ is given by

[TABLE]

Thus we can bound the value of $R_{1}$ by

[TABLE]

Then we can reasonable neglect $R_{1}$ in the equation (63). Now we can express $\tau_{\textrm{GW}}^{\textrm{red}}$ in terms of the imaginary part of the complex exponential and write, since $\Theta(x,\alpha)$ is quadratic in $x$:

[TABLE]

where $B(\alpha)$ is defined as $B(\alpha)\equiv\int_{-1}^{0}\differential x\;e^{i\lambda(x-x^{\ast})^{2}}$, $x^{*}$ satisfies $\evaluated{\partialderivative*{\Theta(x,\alpha)}{x}}_{x=x^{*}}=0$, thus

[TABLE]

and $\lambda$ is given by

[TABLE]

The integral $B(\alpha)$ can be written in terms of the error function, giving

[TABLE]

where $u^{*}\equiv\sqrt{\lambda}x^{*}$. Using the asymptotic expansion of the error functions for $u^{*}\gg 1$ (see Abramowitz and Stegun, 1972), we can write

[TABLE]

Inserting the last expression into (66), $\tau_{\textrm{GW}}^{\textrm{red}}$ becomes

[TABLE]

where $C=H_{0}Z^{2}\Omega/2c^{2}$. From this expression we can see that the maximum of $\tau_{\textrm{GW}}^{\textrm{red}}$ clearly happens for $u^{*}\to 0$. This condition implies, from (68), that the angle corresponding to the maximum absolute value of $\tau_{\textrm{GW}}$ satisfies $x^{*}=0$, or, rearranging the terms, the approximation formula (31). In order to justify the validity of the asymptotic expansion, we can explore around $u^{*}=0$, finding that for a variation in the angle $\Delta\alpha$, then $u^{*}\sim i\sqrt{\frac{Z\Omega}{c}}\Delta\alpha\sim 10^{4}\Delta\alpha$. Thus, the expansion is well defined for $\Delta\alpha\gg 10^{-4}$.

Appendix D Table of pulsars of the ATNF catalog