Rotations of the polarization of a gravitational wave propagating in Universe

In this paper, we study the polarization of a gravitational wave (GW) emitted by an astrophysical source at a cosmic distance propagating through the Friedmann-Lema\^itre-Robertson-Walk universe. By considering the null geodesic deviations, we first provide a definition of the polarization of the GW in terms of the Weyl scalars with respect to a parallelly-transported frame along the null geodesics, and then show explicitly that, due to different effects of the expansion of the universe on the two polarization modes, the so-called"+"and"$\times$"modes, the polarization angle of the GW changes generically, when it is propagating through the curved background. By direct computations of the polarization angle, we show that different epochs, radiation-, matter- and $\Lambda$-dominated, have different effects on the polarization. In particular, for a GW emitted by a binary system, we find explicitly the relation between the change of the polarization angle $|\Delta \varphi|$ and the redshift $z_s$ of the source in different epochs. In the $\Lambda$CDM model, we find that the order of $|\Delta \varphi|{\eta_0 F}$ is typically $O(10^{-3})$ to $O(10^3)$, depending on the values of $z_s$, where $\eta_0$ is the (comoving) time of the current universe, and $F\equiv\Big(\frac{5}{256}\frac{1}{\tau_{obs}}\Big)^{3/8}\left(G_NM_c\right)^{-5/8}$ with $\tau_{obs}$ and $M_c$ being, respectively, the time to coalescence in the observer's frame and the chirp mass of the binary system.


I. INTRODUCTION
With the detection of gravitational waves (GWs) by the LIGO collaboration [1,2], a new era, the gravitational wave astronomy, began, after exactly 100 years since Einstein first predicted the existence of GWs [3] by using his brand new theory of general relativity, established only one year earlier. The masses of the two binary black holes (BBHs) in the event GW151226 were, respectively, 14M and 7M [2], well in the range of stellar-mass black holes known so far [4]. On the other hand, in the event GW150914 [1], the masses of the two black holes were, respectively, 36M and 29M , which are much more massive than the known ones [4]. This has already stimulated lots of interest and various scenarios have been proposed, including their formation from very massive stars (50−100M ) [5], PopIII stars [6], and primordial black holes [7](see also [8]). Despite of the difference among the masses of BBHs in these two events, the distances of them to Earth are almost the same, about 400Mpc, which corresponds to a redshift of z 0.09 [1,2]. After that, LIGO/Virgo scientific collaborations detected dozen more GWs [9][10][11][12][13][14][15][16], including possibly the coalescence of neutron-star (NS)/BH binary, although some details of these detections have not been released yet [17]. Over the next few years, the advanced LIGO and Virgo will continuously increase their sensitivities, and once at their designed goals, together with other detectors of the second generation, such as KAGRA [18], they could be able to detect heavy BBHs up to redshifts of unity. The third generation detectors, both ground and space based, such as the Einstein Telescope [19], Cosmic Explorer [20], LISA [21], Taiji [22], Tianqin [23], and DECIGO [24], will substantially increase their sensitivities and allow us to detect GWs from sources up to redshifts of z 20, well within the epoch of reionization [25].
Therefore, the studies of such GWs and their sources, among other things, will open a new window to explore the early universe.
Such studies are normally divided into two different zones, the generation and propagation zones [26]. The propagation zone is further divided into three different phases, the inspiral, plunge/merger and ringdown phases. In the inspiral phase, the post-Newtonian (PN) approximations are usually used, while in merger phase, the gravitational field is very strong and the field equations (of Einstein as well as of other theories), become highly nonlinear, and heavy numerical calculations are often inevitably used. In the ringdown phase, quasi-normal modes (QNM) analysis is found sufficient.
When GWs are far from their sources (in the far propagation zone), they become very weak, and can be considered as linear perturbations. Previous studies were mainly focused on perturbations in the Minkowski background [27], but recently studies of perturbations in the Friedmann-Lemaître-Robertson-Walk (FLRW) background have started to attract attention [28][29][30][31][32] (See also [33]).
In particular, for a GW generated by an astrophysical source with a distribution of matter in a finite region, such as an inspiraling compact binary system, it can be considered as perturbations on the FLRW backgroundĝ µν ≡ a 2 η µν 1 , g µν = a 2 (η) (η µν + h µν ) , (1.1) where η µν = diag. (1, −1, −1, −1), η denotes the conformal time of the FLRW universe, and | h µν | 1. In this paper, we shall use the notations and conventions of [34], and in particular, we shall set the speed of light to unity, c = 1. In the transverse-traceless (TT) gauge, with δ ij D ij = 0 = ∂ i D ij , the linear perturbations D ij 's equations are given by [29], where ≡ĝ µν∇ µ∇ν ,∇ µ denotes the covariant derivative with respect toĝ µν , and Π ij denotes the part of matter perturbations δT i j , subjected to the TT constraints, δ ij Π ij = 0 = ∂ i Π ij . When far from the source, | x| D, where D denotes the spatial extent of the source, and | x| ≡ x 2 + y 2 + z 2 is the comoving distance between the center of the source and the observer.
In the FLRW background (1.1), when GWs propagate through this curved space-time, an additional term appeared due to the scattering off the curved backgorund, the "backscattering waves" or "tails" 2 D ij consists of two parts: one travels along the light-cone, denoted by D (γ) ij , which will be referred to as the propagation (or direct) part, and the other travels inside the light-cone, denoted by D (tail) ij , the so-called tails as shown in Fig. 1. They are given, respectively, by, a(η) , (1.4) where v ≡ 1 √ 2 (η − | x|), and where R 3 is the region of source, and Θ(x) denotes the Heaviside function, satisfies the equation, with where a (η) ≡ d 2 a(η)/dη 2 , etc. When a(η) is a constant, for which the background becomes the Minkowski spacetime, the last term of Eq.(1.8) vanishes, so Eq.(1.7) reduces to (∂ 2 η − ∂ 2 z )J(η, η ;σ) = 0, which represents waves propagating along the light-cone η = ±z, that is, the tail part reduces to the direct part, and the corresponding GWs all move exactly along the light-cones.
In this paper, we shall study the polarizations of GWs produced by an astrophysical source, which propagate through our universe over a cosmic distance. As a first step, we shall assume that our universe is homogeneous and isotropic. The efforts of the inhomogeneity will be considered somewhere else [36].
The rest of this paper is organized as follows. In Section II, assuming that the GWs are far from their sources, we show how to define the polarization angle of the GW propagating through the flat FLRW background at a cosmic distance. In Section III, different epochs, in which the universe is dominated by different components of matter field, are considered. In particular, we calculate explicitly the changes of the polarization angle in each of these epochs. Our main conclusions are summarized in section IV.
There are also four appendices, A, B, C and D, in which some mathematical calculations involved in this paper are given in detail. In particular, in Appendix C we consider the timelike geodesics deviations, from which we define the polarization of a GW propagating through the flat FLRW universe, and show that it is the same as that obtained by considering null geodesic deviations studied in Section II.

II. POLARIZATIONS OF GRAVITATIONAL WAVES
In the rest of this paper, we shall assume that the GWs are far from their sources, so they are well-described by Eqs.(1.4)-(1.8). In addition, we shall choose the coordinates so that the z-axis is passing through the observer and the center of the source. Then, the transverse-traceless conditions lead to, where σ = +, ×, u ≡ 1 √ 2 (η + | x|), and R ≡ a(η)| x| is the physical distance from the center of the source to the observer, n is a constant, and In the radiation-dominated epoch we have T σ (u, v) = 0, while in the de Sitter spacetime, or the matter-dominated epoch, we have T σ (u, v) = 0, but with n = 0, 3/2, respectively [29].
In writing the above expressions, we expanded the direct part D ij (t, x) only to the first-order of 1/| x|. P σ (v) is related to the integration of the direct part of the source via the relation, Then, the analysis of propagation and polarizations of GWs in the Universe can follow the one given in [37][38][39][40][41][42] 3 , which is crucially based on the following decompositions of the Riemann tensor, where R µν and R denote, respectively, the Ricci tensor and scalar, and C µνλβ is the Weyl tensor. The latter is considered as representing pure gravitational fields, while R µν is connected to the energymomentum tensor T µν through the Einstein field equations R µν = 8πG N (T µν − g µν T /2), which is usually considered as the spacetime curvature produced by matter.
In terms of C µνλβ and R µν , the Bianchi identities take the form, which represents the interaction between matter and pure gravitational fields (represented by the is directly related to the energy-momentum tensor T µν through the Einstein field equations, with To study the phenomena, it is found very convenient to use the Newman-Penrose formula [43], a tetrad formula but with a particular choice of the tetrad, where an overline/bar denotes the complex conjugate, and It can be shown that, to the leading order of , each of l µ and n µ defines a null geodesic congruence, where ∇ ν denotes the covariant derivative with respect to g µν , and A ,η ≡ ∂A/∂η, etc. Projecting the Weyl and Ricci tensors onto the above null tetrad, we find that the ten independent components of the Weyl tensor are given by the five Weyl (complex) scalars, Ψ A , while the ten independent components of the Ricci tenor are given by the seven Ricci scalar, Φ AB , given explicitly in Appendix A for the metric Eq.(2.5). In particular, we have where The importance of the above decompositions is that each component of these scalars has its own physical interpretation. For example, Ψ 0 and Ψ 4 represent the gravitational waves propagating along the null geodesics, defined by l µ and n µ , respectively, while Ψ 1 and Ψ 3 represent the longitudinal components, and Ψ 2 the Coulomb component [44].

A. Polarizations of GWs
To study the polarization of GWs propagating in the FRLW background, let us consider null geodesic deviations. The main motivation of studying null geodesic deviations, instead of those timelike which are related to observations directly, is that GWs are traveling along their null geodesics before catching by the detectors. As shown in [37,38], the considerations of null geodesic congruences will lead to the same definition for the polarization of a given GW. However, in [37,38] only plane gravitational wave (Petrov Type N)spacetimes were considered, and in what follows we will show that this is also true for the spacetimes considered here, by considering null geodesic deviations.
In fact, in Appendix C we show that the definition of the polarization of a GW should be valid in any given spacetime, independent of its symmetry and the nature of the geodesic deviations, null or timelike.
To our purpose, in the following let us consider the null geodesic congruences formed by l µ . We construct the following four unity vectors, Let η µ be the geodesic deviation between two neighbor geodesics and η µ l µ = 0. Then, after some tedious but straightforward calculations, we find that the null geodesic deviation is given by, where e µν + ≡ e µ 2 e ν 2 − e µ 3 e ν 3 , e µν × ≡ e µ 2 e ν 3 + e µ 3 e ν 2 , e µν l,2 ≡ e µ 2 l ν + e ν 2 l µ , e µν l,3 ≡ e µ 3 l ν + e ν 3 l µ , e µν o ≡ e µ 2 e ν 2 + e µ 3 e ν 3 . (2.19) The first term of Eq.(2.18) appears due to the background that is curved, in which h σ is the function of u and v. If the background is flat, then h σ is a function of u or v only, and, as a result, this term will vanish. For matter that satisfies the energy conditions [45], we always have Φ 00 > 0. Thus, the second term always makes the geodesic congruence uniformly contracting in the plane orthogonal to the propagation of the GWs, spanned by e µ 2 and e µ 3 . In contrast to matter, the GW has different effects. In particular, the Re (Ψ 0 ) part stretches the e 2 -direction and meantime squeezes the e 3 -direction or vice versa, depending on the sign of Re (Ψ 0 ).
The Im (Ψ 0 ) part does the same, but now along the axes e 2 and e 3 , which are obtained by a 45 0 rotation in the (e 2 , e 3 )-plane. On the other hand, making the following rotation, Eq.(2.18) becomes, where e µν + = cos 2ϕ 0ê µν In [37], the angle ϕ 0 defined in Eq.(2.21) was referred to as the polarization angle of the GW with respect to the (e 2 , e 3 )-frame. Note that the definition of the polarization angle is gauge invariant as shown explicitly in Appendix A. In addition, to the first-order of h σ , this frame is of parallel transport along the null geodesics defined by l µ (See Appendix B for details), Similarly, if we consider the null geodesic congruence defined by n µ , we will find [37], with ζ µ n µ = 0. If we make rotation of the kind given by Eq.(2.20) but now with an angle ϕ 4 , we shall obtain an expression similar to Eq.(2.21), and in particular we have Thus, ϕ 4 defines the polarization angle of the Ψ 4 GW, moving in the opposite direction of Ψ 0 . It can be shown that the (e 2 , e 3 )-frame is also of parallel transport along the null geodesics defined by n µ , so a such defined polarization angle has an invariant interpretation. The definition of the rotation angle in the FLRW background is consistent with the one in the Minkowski background [37], and it is independent of the nature of the geodesic deviations, null or timelike, as shown in Appendix C.
Note that when h × = 0, we have P × (v) = 0, and that is, the polarization angles of Ψ 0 and Ψ 4 become zero, and all along the e 2 -direction.
In addition, in the Minkowski background, we have H = 0, H = 0 and a = 1, Then from Eqs.(2.13) and (2.14), we find that Hence the polarization angle is given by Thus, along the path of the GW wave, we have that is, the polarization of the wave does not change along its path in the Minkowski background, as it is usually expected.

B. Rotation of Polarization Angles in Curved Universe
However, the polarization angles ϕ 0 and ϕ 4 will change with time along their wave paths once the background is curved, as can be seen from Eqs.(2.21) and (2.27). To study the rotations in more details, let us first introduce the "scale-invariant" quantities via the relations [38][39][40], where Ψ 0 is given by Eq.(A.1), and as can be seen from Eqs.(A.2) in Appendix A. From these expressions we find that which implies that the "×" modes cannot be created from the "+" modes, or vice versa, where a, b = u, v and f, g are certain functions whose explicit forms are not important. To be more specific, let us assume that h + = 0, h × = 0 at a given moment, say, t = t 0 . Then, since the function g does not depend on h + , it vanishes at t = t 0 , too. As a result, the imaginary part of Ψ (0) 0 shall remain constant along the wave path, as one can see from Eq.(2.37). On the other hand, if the function g depends on h + , in general it will not vanish at t 0 even when h × = 0, as now h + = 0. So, the imaginary part of Ψ (0) 0 now can no longer remain constant along the wave path, and the polarization angle defined above will change, too. As shown explicitly in Appendix B, this is possible only when the second-order effects are taken into account, and to the linear order of such effects do not exist.
However, due to the presence of the matter field of the FLRW background, represented by Φ 02 , both of the "×" and "+" modes can get amplified or diluted, depending on their amplitudes and signs, so in principle the effects to the "×" modes are different from those to the "+" modes. As a result, the polarization angle ϕ 0 defined above will change along the wave path. The above analysis was for the GW represented by A similar analysis can be carried out for the Ψ 4 wave, too, and the same conclusion will be obtained, that is, the "×" modes cannot be created from the "+" modes, or vice versa, but the polarization angle of the Ψ 4 wave can still get changed, due to the different effects of the matter presented in our universe to the two different kinds of polarization modes.
Therefore, to the linear order of , there is no transfer between the two different kinds of modes.
However, due to the fact that the effects of the curved background on each of the two independent modes depend on their waveforms, the polarization angle can be still changed along the wave path.
These can be seen more clearly from the following studies of the polarization angles in various epochs of the universe.
Before proceeding to the next section, it should be emphasized that the conclusion that there is no transfer between the two different modes holds only in the linear level (of ), as shown explicitly in [37,38] that the nonlinear interaction of the GW wave with the background will lead to the phenomena of the Faraday rotation. holds for any GW sources. However, in order to have more Intuitive understanding, let us consider a binary system as a concrete example.

A. Radiation-Dominated Epoch
The scale factor in general is given by a power law a = (η/η 0 ) p in terms of the conformal time.
In particular, in the radiation-dominated epoch we have p = 1. Then, the tail vanishes T = 0 [29], and from Eqs.(2.13) and (2.14) together with z x, y, we find that Substituting the above equations into Eq.(2.21) and then performing the series expansion in terms of η, we find that In the second step of Eq.(3.3), we have substituted the result of P σ for the binary system, which can be found explicitly in Appendix D. And it is obvious that g(v) only depends on GWs frequency Suppose that the radiation-dominated epoch is through the whole cosmic history, so that η 0 is the age of the universe, then the change of the polarization angle from a moment η e to the Earth along the wave path (v = Constant) can be obtained from Eq.(3.2), which is given by 4 , where v e is the corresponding retarded time, which is a constant along the GW path, and a(η 0 ) = 1, a(η e ) = 1/(1 + z s ). Then, we find that Thus, in the source frame, we have where τ s denotes the time to the coalescence measured by the source's clock, and f (s) gw (τ s ) is the frequency in the source frame. However, it is more convenient to express the polarization angle in the observer's frame, which is given by where f (obs) gw (τ obs ) and τ obs are the frequency and time to the coalescence in the observer's frame. In Fig. 2, we plot |∆ϕ 0 | as a function of z s . we set 3600 ≤ z s ≤ 4000.

B. Matter-Dominated Epoch
In the matter-dominated epoch, we have a = (η/η 0 ) 2 and, (3.10) Then, from Eq.(2.2) we find that (v). Hence, substituting n = 3 2 and in the limit x, y << z we obtain It is obvious that C(v) is highly suppressed by P σ (v)/R while taking η 0 as the age of the universe, where a(η 0 ) = 1, ∆t peakwidth dηa(η), and η * is the peak time of the source's strength. Then, we find Assume that the matter-dominated epoch is through the whole cosmic history, then the accumulation of the polarization angle from the source to the earth is given by, (3.15) Thus, the rotation angle in observer's frame is (3.16) We can also plot |∆ϕ 0 | as a function of z s , which is given in Fig. 3.

C. de Sitter Background
In this case, the tail D (tail) ij (η, x) takes the form [29], (3.17) In the present case, we have (v) is related to the integration of the source, as that given by Eq. (3.17). With n = 0 and a(η) = (η 0 /η), from Eqs.(2.13) and (2.14), we find that Similar to the matter-dominated case, now D σ (v) is bounded by where ∆t peakwidth dηa(η), and η * is the peak time of the source's strength. Then, we find that, where η 0 = − 1 H . Substituting Eq.(3.19) into Eq.(2.21) and then carrying out the series expansion, we finally find Similarly, supposing that the de sitter epoch is through the whole cosmic history, we have (3.23) Thus, the polarization angle in observer's frame can be expressed as, To see this effect more explicitly, we plot |∆ϕ 0 | as a function of z s in Fig. 4.

D. More realistic considerations
In the last subsections, we have assumed that there is only one epoch throughout the whole cosmic history. The real universe is obviously not the case. In this subsection, we will consider a more realistic case by combining all these epochs together within the framework of the ΛCDM model.
We consider a binary system emitted GW signals at some time η e . Then, we should integrate the polarization rotations from η e to η Λ , the receiving time of detectors. The whole picture is sketched in Fig. 5. Obviously, we have (3.25) The fucntion g(v e ) in observer's frame is given by where z s is the redshift of the source at η e , and z r 3600 at η r [46], z m 0.4 at η m [46]. Ω Λ , Ω M and Ω R are the density parameters for dark energy, matter and radiation, respectively. H 0 is the Hubble constant.
Depending on the values of z s , the total rotation angle of the whole propagation process can be divided into the following three cases.
• z s > 3600, that is, η e < η r . The accumulated polarization angle from η e to η Λ is where |∆ϕ 0 | r is the polarization angle from η e to η r , |∆ϕ 0 | m is the polarization angle from η r to η m , and |∆ϕ 0 | Λ is the polarization angle from η m to η Λ .
At the end, we use the Plank 2015 data [47]: Ω Λ = 0.6935, Ω m = 0.3065, Ω R = 0. Thus, the accumulated polarization angle in the whole propagation process are shown in Fig. 6 and Fig. 7. To show the rotation angle more intuitively, let us consider the following two examples: • The binary system at redshift z s = 60. The range of frequency is from 10 Hz to 1000 Hz.
We plot ∆ϕ 0 as a function of frequency f (obs) gw in Fig. 8 where the source redshift is z s = 60.
• The binary system at redshift z s = 10. The range of frequency is from 10 −4 Hz to 10 −2 Hz.
We plot ∆ϕ 0 as a function of frequency f (obs) gw in Fig. 8 where the source redshift is z s = 10. typical sources of LIGO/Virgo, we find the typical value for |∆ϕ| is 10 −21 as shown in Fig. (8), which is very small and is far from being detected by the current detectors. However, for lower frequency, the observation effect increases linearly. Hence, it is possible for future GW detectors sensitive in extremely low frequency.
It should be noted that there are two timescales: One is related to the propagation of GWs, the other refers to the period of detections of GWs. The former can be comparable with the age of the universe for astrophysical GW sources with high redshift, while the latter is of the order of seconds for the current ground-based detectors. Therefore, the changes of polarizations of GWs cannot be directly related to the current detections of GWs. However, the effects of the background on the propagations and polarizations of such GWs cannot be ignored as shown in section III, especially for sources with low frequency and high redshift.
It would be very interesting to find the observational signature of the rotations of the polarization.
However, since the changes are due to the curvature of the background, they become significant only over cosmic scales. Therefore, to detect such effects, combinations of cosmic distant objects with the ground-and/or space-based detectors might be needed. Interestingly, analysis about the amplitude and phase of GWs traveling a long cosmological distance has been discussed a lot in [53].
In addition, in this paper we assumed that the GWs are propagating on the homogeneous and isotropic background. This implies that the wavelengths of such GWs are much shorter than the scale, over which the changes of the universe including its perturbations are negligible. This automatically ignores the gravitational lensing of such GWs by massive objects, such as galaxies and supermassive black holes. It would be very interesting and important to study such effects, specially their observational evidences.

Acknowledgements
We would like very much to express our gratitude to Prof. Shinji Mukohyama for his long time collaboration on the subjects, and valuable comments and suggestions, which lead us to sharp the content considerably, and make the current version more readable. We are very grateful to the anonymous referee for her/his valuable comments and suggestions. This work was partially In this appendix, we present some detailed computations of the Weyl and Ricci scalars in terms of the null tetrad. In particular, we find, When | z | x and | z | y, we find, Then from (2.2) we obtain Inserting them into Eqs.(A.1) and (A.2), we get Eqs.(2.13), (2.14) and When it comes to gauge invariance, there are two kinds of gauge transformations: transformations of the first kind as defined in [48].
In summary, Ψ 0 and Ψ 4 are invariant under any kind of the gauge transformations.

Appendix B: Parallel-transported Polarization Angle
It was shown in [37] that the above definition of the polarization angle is in general not paralleltransported along the wave paths. This makes it impossible to compare the polarization of a gravitational wave at two different points along its propagation path. One way to overcome this difficulty is to find a parallel-transported basis carried by the wave. Once we have this basis, polarization angle defined based on this basis should be parallel-transported. As suggested in [37], a paralleltransported basis can be found by rotating the coordinate with a proper angle in the (e µ 2 , e µ 3 ) plane. Similar considerations can be generalized to the current case. In particular, for GWs along the null geodesics defined by n µ , we let ϕ (0) 0 be the angle of rotating, and λ µ (2) and λ µ (3) be the new basis. Then, we find Then, the angle determines the polarization direction of the Ψ 0 wave relative to the (λ µ (2) , λ µ (3) ) basis. Similarly, for the Ψ 4 wave, the related angle is given by In this appendix, let us consider timelike geodesic deviations, as they are directly related to observations. To our purpose, in the following we consider timelike geodesic congruences formed by t µ , from Eqs.(2.16) we find Therefore, choosing A and B so that we find t µ ;ν t ν = 0, that is, t µ defines a timelike and affinely parametrized geodesic congruence. Recall that 2AB = a −2 . Thus, the general solution of the above equation is where A 0 is an integration constant. Without loss of the generality, we can always choose A 0 = 0, so that In what follows, we shall adopt the "+" sign.
On the other hand, the four spatial unity vectors we constructed in section II, form an orthogonal base. Thus, the component Ψ 0 (Ψ 4 ) represents the GW moving along the positive (negative) s µdirection, as shown in Fig. 1.
Let η µ be the geodesic deviation between two neighbor geodesics and η µ t µ = 0. Then, the timelike geodesic deviation is given by, where e µν ×2s ≡ e µ 2 s ν + e ν 2 s µ , e µν ×3s ≡ e µ 3 s ν + e ν 3 s µ , (C.6) The right hand side of Eq.(C.5) includes six terms, each of them represents a sort of polarization mode, and has the following physical interpretations [38,44]. The first terms is transverse, laying on the (e 2 , e 3 )-plane only, and makes the test particles contract or expand uniformly, depending on sign of the coefficient. It is remarkable to note that the Weyl scalars have no contributions to this term.
The second term has the effect of distorting a sphere of particles about the observer into an ellipsoid, which has s µ as principal axis, while in the perpendicular (e 2 , e 3 )-plane it is uniformly contracting or expanding, depending the signs of the coefficient, Ψ 2 +Ψ 2 + 2 (Φ Λ − Φ 11 ). This is a behavior quite similar to particles falling in towards a central attracting body with the inverse square law, so it is often to call this term as the Coulomb part of the field, and the strength of the Coulomb field of the pure gravitational field is given by the real part of Ψ 2 , while the strength of the matter field is represented by (Φ Λ − Φ 11 ). The third term represents a force that acts only on the (e µ 2 , s µ )-plane. we find that Therefore, if it is initially a circle in the (e µ 2 , s µ )-plane, this term will make the circle into an ellipse with its major axis along the 45 0 with respect to the s µ -direction. The strength of this force depends on the real part of the combination, ( In contrast to the third term, the fourth term represents a force that acts only on the (e µ 3 , s µ )plane, and the strength of this term depends on the imaginary part of the combination, (Ψ 1 − Ψ 3 ) − (Φ 01 − Φ 12 ). Note the difference between this term and the last one in the pure gravitational part, represented by Ψ 1 and Ψ 3 , respectively. In the third term, it is the real part of (Ψ 1 + Ψ 3 ), while in the fourth term, it is the imaginary part of (Ψ 1 − Ψ 3 ). This is different from the matter field, which is always proportional to the combination of (Φ 01 − Φ 12 ) in both terms.
The fifth term is also purely transverse and deforms a circle laying in the (e 2 , e 3 )-plane into ellipse with its major axis along the e 2 -direction. To see the effects of the last term, we can make a similar rotation as that given by Eq.(C.7) but now in the (e 2 , e 3 )-plane, so we find that e µν × = −e µν + .
(C.9) Therefore, the last term is also purely transverse and deforms a circle laying in the (e 2 , e 3 )-plane into ellipse with its major axis along a line that has a 45 0 angle with respect to the e 2 -direction.
To consider the polarization of the right-moving (along the positive direction of s µ ) GW represented by Ψ 0 , let us first make the rotation as Eq.(2.20), then we find that the Ψ 0 parts in the right-hand side of Eq.(C.5) become, Re (Ψ 0 ) e µν + + Im (Ψ 0 ) e µν × = Ψ 0Ψ0 1/2êµν That is, the Ψ 0 GW makes a circle in (e 2 , e 3 )-plane into ellipse with its major axis along theê 2direction. Such a defined angle is referred as to the polarization angle of the Ψ 0 GW [37,38].
Note that, to the first-order of h σ , the two spatial unit vectors e 2 and e 3 are of parallel transport along the time geodesics defined by t µ , so that the angle ϕ 0 defined above has an absolute meaning.
It must be noted that the matter component Φ 02 has similar effects on the circle laying in the (e 2 , e 3 )-plane, as we can see from the last two terms in Eq.(C.5). Therefore, this terms also makes the circle into an ellipse, with the same effects precisely as the Ψ 0 GW component, and in principle the detectors, such as LIGO and Virgo, cannot distinguish these effects, and will measure them together. However, from Eqs.(2.13)-(2.15) and (A.5)-(A.6), we find that, relative to Ψ 0 , the matter component Φ 02 is either suppressed by a factor H or (x 2 + y 2 )/z 2 . So, comparing with Ψ 0 , the effects of Φ 02 are negligible, and can be safely ignored for the current and forthcoming detectors, including the third generation of detectors.
With the above in mind, for the Ψ 4 wave, if we make a similar rotation as that given by Eq.(C. 10) but now with the angle given by, tan 2ϕ 4 ≡ + Im (Ψ 4 ) Re (Ψ 4 ) , (C.14) we find that [37,38], Re (Ψ 4 ) e µν + − Im (Ψ 4 ) e µν × = Ψ 4Ψ4 1/2êµν + , (C. 15) and such defined angle ϕ 4 is referred to as the polarization angle of the Ψ 4 wave. Again, in the current case it has an absolute meaning, as the frame defined by e 2 and e 3 are of parallel transport along the time geodesics.
In summary, the definition of the polarization of a GW is valid in any given spacetime, independent of its symmetry and the nature of the geodesic deviations, null or timelike.
Appendix D: The expression of P σ To calculate P σ , we first note that Π ij = Λ ij,kl δT kl , (D.1) where Λ ij,kl ( n) = P ik P jl − 1 2 P ij P kl and P ij ( n) = δ ij − n i n j . Then we get where S ij = R 3 d 3 x a 3 (v)δT ij (v, x ). We also define It is not difficult to prove that the following relation holds where a dot denotes the derivative with respect to (w.r.t.) the cosmic time t, with dt ≡ a(η)dη.
Following [51], with n = (sin θ sin φ, sin θ cos φ, cos θ), (D. 6) we find Particularly, let us consider a binary system with masses m 1 and m 2 , and also assume that the motion is circular. For simplification, we also assume that the orbit lies in the (x, y) plane, such that the relative coordinates x i 0 are x 0 = R cos(ω s t + π 2 ), (D.9) y 0 = R sin(ω s t + π 2 ), (D.10) where R, ω s are the orbital radius and frequency, respectively. Then, we find that M ij is given by M ij = µx i 0 x j 0 with µ = m 1 m 2 /(m 1 + m 2 ). By using Kepler's law, we finally get that where t ret , f gw and M c are the retarded time, frequency of GWs and the chirp mass respectively, ι is the angle between the normal to the orbit and the line-of-light which equals to θ. On a quasi-circular orbit with R = R(t), ω s = ω s (t), we define [51,52] Φ(t) ≡ 2 of the source frame and (x,ȳ,z) is the coordinate system of observer's frame. 5 ηret ≡ η − | x| = √ 2v is the retarded time, ηc is the value of η at the coalescence while vc is equal to the corresponding retarded time, andτ denotes the time interval to the coalescence.