Constraining Screened Modified Gravity with Spaceborne Gravitational-wave Detectors

General relativity (GR) is always considered as the most successful theory of gravity. However, various difficulties of this theory are also well known. For instance, in the theoretical side, GR has the singularity and quantization problems (DeWitt 1967; Kiefer 2007). In the experimental side, all the observations in cosmological scale indicate the existence of so-called dark matter and dark energy, which might hint at the invalidity of GR at this scale (Sahni 2004; Cline 2013). For these reasons, since GR was proposed by Einstein in 1915, a large number of experimental tests have been performed on various scales, from submillimeter-scale tests in the laboratory to tests at solar system and cosmological scales (Adelberger 2001; Hoyle et al. 2001; Jain & Khoury 2010; Will 2014; Burrage et al. 2015; Koyama 2016; Bertoldi et al. 2019; Sabulsky et al. 2019). Unfortunately, most of these efforts have focused on the gravitational effects in weak fields. Since the observable gravitational-wave (GW) signals can only be generated in strong gravitational fields and are nearly freely propagating in spacetime once generated (Maggiore 2008), there is an excellent opportunity to experimentally test the theory of gravity in the strong-field regime (Abbott et al. 2016a, 2019a, 2019b). Recently, with the discovery of compact binary coalescence GW signals by the aLIGO and aVirgo collaborations (Abbott et al. 2016b, 2016c, 2016d, 2017a, 2017b, 2017c, 2017d, 2019c), testing GR in the strong gravitational fields becomes one of the key issues in GW astronomy (Kostelecký & Mewes 2016; Miller & Yunes 2019; Sathyaprakash et al. 2019).

The testing of GR by GW observations entails comparing the predictions of GW signals in GR and those in the alternative theories and constraining their differences by observations. Therefore, the choice of typical alternative gravitational theory and the calculation of GW waveforms in the theory have crucial roles (Yunes & Siemens 2013; Berti et al. 2018). A natural alternative to GR is the scalar–tensor theory, which invokes a conformal coupling between matter and an underlying scalar field (see, for instance, Brans–Dicke gravity; Brans & Dicke 1961), besides the standard spacetime metric tensor. The coupling between scalar field and matter leads to the scalar force (fifth force), and the tight experimental constraints (Adelberger et al. 2009; Williams et al. 2012) require that the fifth force must be screened in high-density environments. In our previous works (Zhang et al. 2016, 2017, 2019a, 2019b, 2019c; Liu et al. 2018b), we studied the general scalar–tensor gravity with screening mechanisms, which can suppress the fifth force in dense regions and allow theories to evade the solar system and laboratory tests, in a unified theoretical framework called screened modified gravity (SMG). In this framework, the chameleon, symmetron, dilaton, and f(R) models in the literature are the specific cases of this theory. We have calculated the parameterized post-Newtonian (PPN) parameters (Zhang et al. 2016), the post-Keplerian (PK) parameters (Zhang et al. 2019b), the effective cosmological constant (Zhang et al. 2016), the effective gravitational constant (Zhang et al. 2016), and the change in the orbital period of the binary system caused by the gravitational radiation (Zhang et al. 2017). Based on these results, we have derived the constraints on the model parameters by considering the observations in solar system, cosmological scale, binary pulsar, and lunar laser ranging measurements (Zhang et al. 2019a, 2019b, 2017, 2016). In addition, in Liu et al. (2018b), we calculated in detail the GW waveforms produced by the compact binary coalescences during the inspiraling stage, and we derived the deviations from that in GR, which are partly quantified by the parameterized post-Einsteinian (PPE) parameters. Utilizing these results, we also obtained the potential constraints on the theory by the future ground-based Einstein telescope.

In addition to the ground-based GW detectors, space-borne detectors are also proposed. In the near future, the missions Laser Interferometer Space Antenna (LISA), Taiji, and TianQin will be launched around the 2030s (Danzmann et al. 2016; Luo et al. 2016; Hu & Wu 2017). Due to the large arm lengths of these detectors, the sensitive frequency ranges become (10⁻⁴, 10⁰) Hz, lower than those of ground-based detectors, and the extreme-mass-ratio inspirals (EMRIs) are the important GW sources (Amaro-Seoane et al. 2007; Babak et al. 2017). The event rate of the EMRIs is difficult to estimate because it depends on factors that are poorly constrained by observation. According to the estimations of Babak et al. (2017), at least a few EMRIs per year can be detected by LISA, irrespective of the astrophysical model. For the most optimistic astrophysical assumptions, this number can reach a few thousand per year. An EMRI normally consists of a stellar compact object, such as a white dwarf (WD), neutron star (NS), or stellar-mass black hole (BH), and a massive BH, which is an excellent source for the test of gravity (Barack & Cutler 2007; Gair et al. 2013). In previous works (Scharre & Will 2002; Will & Yunes 2004), the authors have investigated the constraints on Brans–Dicke gravity, massive gravity, and so on, assuming the GW signals of BH–NS binaries observed by the LISA mission. Similarly, in this article, we will study the constraints on SMG theory by the GW signals produced by the BH–WD and BH–NS binaries. In our discussion, we will consider both the LISA and TianQin missions. Taiji is similar to LISA (Wu et al. 2019), so we suspect the potential constraint from Taiji is also similar to that from LISA. In the calculation, we consider three different cases for the detection. In case one, we constrain the SMG by detecting the extra GW modes. In case two, we constrain the theory by Fisher information matrix analysis, but consider only the restricted GW waveforms, and in case three, we do the same analysis but including the higher order amplitude corrections in the templates. In comparison with the results in these cases, we investigate the contributions of extra polarization modes and the higher order amplitude corrections in the model constraints.

This paper is organized as follows. In Section 2, a brief introduction to screened modified gravity is presented, and the Fourier transform of the GW waveforms in SMG is rewritten for convenient reference. Two aspects of detector information that are relevant to our analysis, noise curves and antenna pattern functions, are introduced in Section 3. In Section 4, the method employed in this work and the process used to get the constraints are shown in detail. The results are presented and discussed in Section 5, where we compare the constraints given by the forecasts of future space-borne detectors with the constraints obtained by the current experiments in the three specific SMG models. The full waveform of 2.5 PN in amplitude and 3.5 PN in phase with the corrections concerned with the SMG and the process to derive the antenna pattern functions are given in Appendixes A and B.

Throughout this paper we adopt the units where $c={\hslash }=1$. The reduced Planck mass is ${M}_{\mathrm{Pl}}=\sqrt{1/(8\pi G)}$, where G denotes the Newtonian gravitational constant. Because in this article we consider only the GW sources in the very low redshift range, the redshifts are not explicitly expressed in the formulae of this paper. The distance parameter denotes the luminosity distance, and the chirp mass and total mass in this paper denote the directly measured values in the detector frames.

The GW waveforms of binaries with circular orbits in general SMG have been calculated in previous work (Zhang et al. 2017; Liu et al. 2018b). In this section, we will present a brief introduction to screened modified gravity and rewrite the formulae of waveforms in the SMG. The action of a general scalar–tensor theory in the Einstein frame is given by

$$\begin{eqnarray}\begin{array}{rcl}S & = & \int {d}^{4}x\sqrt{-g}\left[\displaystyle \frac{{M}_{\mathrm{Pl}}^{2}}{2}R-\displaystyle \frac{1}{2}{\left({\rm{\nabla }}\phi \right)}^{2}-V(\phi )\right]\\ & & +\,{S}_{m}\left[{A}^{2}(\phi ){g}_{\mu \nu },{\psi }_{m}\right],\end{array}\end{eqnarray} \tag{ 1 }$$

where g_μν is the matrix in the Einstein frame, g is the determinant of the matrix, R is the Ricci scalar, ϕ is the scalar field, and ψ_m is the matter field. Here, V(ϕ) is a bare potential that characterizes the scalar self-interaction, and A(ϕ) denotes a conformal coupling function representing the interaction between scalar field and matter field. In scalar–tensor theory, the scalar field can affect the effective mass of a compact object. As suggested by Eardley (1975), the matter action takes the form

$$\begin{eqnarray}&&{S}_{m}=-\sum _{a}\int {m}_{a}(\phi )d{\tau }_{a},\end{eqnarray} \tag{ 2 }$$

where the constant inertial masses of the compact objects are substituted by a function of the scalar field. The field equations can be obtained by the variation of the action with respect to g_μν and ϕ,

$$\begin{eqnarray}&&{G}_{\mu \nu }=8\pi G({T}_{\mu \nu }+{T}_{\mu \nu }^{\phi }),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }{{\rm{\nabla }}}^{\mu }\phi =\displaystyle \frac{\partial }{\partial \phi }(V(\phi )-T),\end{eqnarray} \tag{ 4 }$$

where T_μν and ${T}_{\mu \nu }^{\phi }$ are the energy–momentum tensor of the matter field and scalar field, respectively, and T is the trace of T_μν. The behavior of the scalar field is controlled by both V(ϕ) and T, by which we define the effective potential

$$\begin{eqnarray}&&{V}_{\mathrm{eff}}=V(\phi )-T.\end{eqnarray} \tag{ 5 }$$

As shown in Zhang et al. (2016), in the SMG, the effective potential can be rewritten as

$$\begin{eqnarray}&&{V}_{\mathrm{eff}}=V(\phi )+\rho A(\phi ),\end{eqnarray} \tag{ 6 }$$

where ρ is the conserved energy density in the Einstein frame. In the wave zone, the metric and the scalar field can be expanded around the flat background ${\eta }_{\mu \nu }$ and the scalar background (the vacuum expectation value (VEV) of the scalar field) ϕ_VEV:

$$\begin{eqnarray}&&{g}_{\mu \nu }={\eta }_{\mu \nu }+{h}_{\mu \nu },\qquad \phi ={\phi }_{\mathrm{VEV}}+\delta \phi .\end{eqnarray} \tag{ 7 }$$

The bare potential V(ϕ) and the coupling function A(ϕ) can be expanded as

$$\begin{eqnarray}\begin{array}{rcl}V(\phi ) & = & {V}_{\mathrm{VEV}}+{V}_{1}\delta \phi +{V}_{2}\delta {\phi }^{2}+{V}_{3}\delta {\phi }^{3}+{ \mathcal O }(\delta {\phi }^{4}),\\ A(\phi ) & = & {A}_{\mathrm{VEV}}+{A}_{1}\delta \phi +{A}_{2}\delta {\phi }^{2}+{A}_{3}\delta {\phi }^{3}+{ \mathcal O }(\delta {\phi }^{4}).\end{array}\end{eqnarray} \tag{ 8 }$$

The effective mass of the scalar field is given by

$$\begin{eqnarray}&&{\left.{m}_{s}^{2}\equiv \displaystyle \frac{{d}^{2}{V}_{\mathrm{eff}}}{d{\phi }^{2}}\right|}_{{\phi }_{\mathrm{VEV}}}=2\left({V}_{2}+{\rho }_{b}{A}_{2}\right),\end{eqnarray} \tag{ 9 }$$

where ρ_b is the background matter density. We can find that the effective mass of the scalar field depends on the ambient matter density. In conditions such as the solar system, the matter density is high and the mass of the scalar field is large, so the range of the force corresponding to the scalar field is too short to have detectable effects. In this way, the effects of the scalar field can be screened in the high-density environment and evade the tight constraints presented by solar system experiments. However, in the large scale, the matter density is low, and the scalar field can have significant effects to accelerate the expansion of the universe. The mass of the compact object m_a(ϕ) can also be expanded as

$$\begin{eqnarray}&&{m}_{a}(\phi )={m}_{a}\left[1+{s}_{a}\left(\displaystyle \frac{\delta \phi }{{\phi }_{\mathrm{VEV}}}\right)+{ \mathcal O }{\left(\displaystyle \frac{\delta \phi }{{\phi }_{\mathrm{VEV}}}\right)}^{2}\right],\end{eqnarray} \tag{ 10 }$$

where ${m}_{a}={m}_{a}({\phi }_{\mathrm{VEV}})$. The sensitivity of the ath object s_a is defined as

$$\begin{eqnarray}&&{\left.{s}_{a}\equiv \displaystyle \frac{\partial (\mathrm{ln}{m}_{a})}{\partial (\mathrm{ln}\phi )}\right|}_{{\phi }_{\mathrm{VEV}}}.\end{eqnarray} \tag{ 11 }$$

In most cases, the deviations from GR are quantified by the sensitivity (Zhang et al. 2016, 2017, 2019a, 2019b; Liu et al. 2018b). In the SMG, it is proportional to the screened parameter

$$\begin{eqnarray}&&{s}_{a}=\displaystyle \frac{{\phi }_{\mathrm{VEV}}}{2{M}_{\mathrm{Pl}}}{\epsilon }_{a},\end{eqnarray} \tag{ 12 }$$

and the screened parameter of a uniform-density object is given by

$$\begin{eqnarray}&&{\epsilon }_{a}=\displaystyle \frac{{\phi }_{\mathrm{VEV}}-{\phi }_{a}}{{M}_{\mathrm{Pl}}{{\rm{\Phi }}}_{a}},\end{eqnarray} \tag{ 13 }$$

where ${{\rm{\Phi }}}_{a}={{Gm}}_{a}/{R}_{a}$ is the surface gravitational potential, and ϕ_a is the position of the minimum of the effective potential inside the object. In this paper, we will focus on the screened parameter and forecast how tight constraints can be placed on it by the future space-borne GW detectors.

In the wave zone, the linear field equations are given by

$$\begin{eqnarray}&&\square {\bar{h}}_{\mu \nu }=-16\pi G{\tau }_{\mu \nu },\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\left(\square -{m}_{s}^{2}\right)\delta \phi =-16\pi {GS},\end{eqnarray} \tag{ 15 }$$

where ${\bar{h}}_{\mu \nu }={h}_{\mu \nu }-\tfrac{1}{2}{\eta }_{\mu \nu }{h}_{\lambda }^{\lambda }$, τ_μν is the total energy–momentum tensor, and S is the source term of the scalar field. The solutions of these equations in the wave zone can be obtained using the Green's function method and are expressed in terms of the mass multipole moments and the scalar multipole moments. Based on the solutions, the GW waveforms in the SMG were calculated in previous work (Liu et al. 2018b).

As shown in Liu et al. (2018b), in addition to + and × polarization modes in GR, the massive scalar field induces two polarizations: breathing polarization h_b and longitudinal polarization h_l. The response of an interferometric detector is given by

$$\begin{eqnarray}&&h(t)={F}_{+}{h}_{+}+{F}_{\times }{h}_{\times }+{F}_{b}{h}_{b}+{F}_{l}{h}_{l}\end{eqnarray} \tag{ 16 }$$

where ${F}_{+,\times ,b,l}$ denotes the antenna pattern functions depending on the direction of GW sources (θ, φ), detector configuration, polarization angle ψ, and the frequency of GWs for space-borne detectors, and ${h}_{+,\times ,b,l}$ denotes gravitational waveforms for the plus, cross, breathing, and longitudinal polarization modes, respectively. Besides the same parameters in waveforms of GR (which are total mass m = m₁+m₂, symmetric mass ratio η = m₁m₂/(m₁+m₂)², chirp mass ${M}_{c}={\eta }^{3/5}m$, distance D, inclination angle ι between the line of sight and the binary orbital angular momentum, the time of coalescence t_c, and the orbital phase of coalescence Ψ_c), there are five extra parameters peculiar to SMG. They are the effective mass of the scalar field m_s, the expansion coefficients of the coupling function A₀ and A₁, and the screened parameters of binary ₁, ₂ (see Zhang et al. 2017; Liu et al. 2018b for more details). The Fourier transform can be obtained by using the stationary phase approximation. The constraint $| {A}_{0}-1|$ is less than 10⁻¹⁰ according to the solar system experiments (Zhang et al. 2016). Therefore, similar to Liu et al. (2018b), we can safely adopt A₀ = 1 in our calculation. As shown in Equation (74) of that paper (Liu et al. 2018b), the difference between the parameters in the Einstein frame and those in the Jordan frame is only a factor A₀. The parameters of the waveform in the Einstein frame and in the Jordan frame are the same when we adopt A₀ = 1. In addition, because the Compton wavelength ${m}_{s}^{-1}$ is roughly cosmological scale (${m}_{s}^{-1}\sim 1\,\mathrm{Mpc}$), as in Zhang et al. (2017), we set m_s = 0 in the waveforms, which makes the "l" polarization vanish (${\tilde{h}}_{l}(f)=0$). The results can be rewritten as follows. Harmonic number one is given by

$$\begin{eqnarray}\begin{array}{rcl}{F}_{b}{\tilde{h}}_{b}^{\left(1\right)}(f) & = & {\left(\displaystyle \frac{5}{48}\right)}^{\tfrac{1}{2}}{\pi }^{-\tfrac{1}{2}}\displaystyle \frac{{\left({{GM}}_{c}\right)}^{5/6}}{D}{\left(2f\right)}^{-\tfrac{7}{6}}\\ & & \times \,\left[-\displaystyle \frac{5}{384}E{\epsilon }_{d}^{2}{\left(2\pi {fGm}\right)}^{-1}+E{\left(2\pi {fGm}\right)}^{-\tfrac{1}{3}}\right]\\ & & \times \,\exp \left\{i\left[2\pi {{ft}}_{c}-\displaystyle \frac{\pi }{4}+{\rm{\Psi }}(f)\right]\right\},\end{array}\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{F}_{l}{\tilde{h}}_{l}^{\left(1\right)}(f)=0,\end{eqnarray} \tag{ 18 }$$

where

$$\begin{eqnarray}&&E=-{F}_{b}{A}_{1}{M}_{\mathrm{Pl}}{\epsilon }_{d}\sin \iota {\left(1+\displaystyle \frac{1}{2}{\epsilon }_{1}{\epsilon }_{2}\right)}^{1/3},\end{eqnarray} \tag{ 19 }$$

with _d = ₁ − ₂. Here, Ψ(f) takes the form

$$\begin{eqnarray}&&{\rm{\Psi }}(f)=-{{\rm{\Psi }}}_{c}+\displaystyle \frac{3}{256{\left(2\pi {{fGM}}_{c}\right)}^{5/3}}\sum _{i=-2}^{7}{{\rm{\Psi }}}_{i}{\left(2\pi {fGm}\right)}^{i/3},\end{eqnarray} \tag{ 20 }$$

where Ψ_c is the orbital phase of coalescence, and the coefficients Ψ_i are presented in Appendix A. Here, Ψ_i(i ≥ 0) is the coefficient in the 3.5 PN phase function of the Fourier domain waveform, and the coefficient Ψ₋₂ is concerned with the correction of dipole radiation. Harmonic number two is given by

$$\begin{eqnarray}\begin{array}{rcl}{F}_{b}{\tilde{h}}_{b}^{\left(2\right)}(f) & = & {\left(\displaystyle \frac{5}{96}\right)}^{\tfrac{1}{2}}{\pi }^{-\tfrac{2}{3}}\displaystyle \frac{{\left({{GM}}_{c}\right)}^{\tfrac{5}{6}}}{D}{f}^{-\tfrac{7}{6}}T\left[{F}_{b}{S}_{-1}{\left(\pi {fGm}\right)}^{-\tfrac{2}{3}}+{F}_{b}\right]\\ & & \times \,\exp \left\{i[2\pi {{ft}}_{c}-\displaystyle \frac{\pi }{4}+2{\rm{\Psi }}(f/2)]\right\},\end{array}\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{F}_{l}{\tilde{h}}_{l}^{\left(2\right)}(f)=0,\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&\begin{array}{l}{F}_{+}{\tilde{h}}_{+}^{\left(2\right)}(f)+{F}_{\times }{\tilde{h}}_{\times }^{\left(2\right)}(f)={\left(\displaystyle \frac{5}{96}\right)}^{\tfrac{1}{2}}{\pi }^{-\tfrac{2}{3}}\displaystyle \frac{{\left({{GM}}_{c}\right)}^{\tfrac{5}{6}}}{D}{f}^{-\tfrac{7}{6}}\\ \qquad \qquad \times \,\left[Q+{{QS}}_{-1}{\left(\pi {fGm}\right)}^{-\tfrac{2}{3}}\right]{e}^{-i{\varphi }_{\left(\mathrm{2,0}\right)}}{P}_{\left(\mathrm{2,0}\right)}\\ \qquad \qquad \times \,\exp \left\{i\left[2\pi {{ft}}_{c}-\displaystyle \frac{\pi }{4}+2{\rm{\Psi }}(f/2)\right]\right\},\end{array}\end{eqnarray} \tag{ 23 }$$

where

$$\begin{eqnarray}&&Q={\left(1+\displaystyle \frac{1}{2}{\epsilon }_{1}{\epsilon }_{2}\right)}^{2/3},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{S}_{-1}=-\displaystyle \frac{5}{384}{\epsilon }_{d}^{2},\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&T=-{A}_{1}{M}_{\mathrm{Pl}}{\rm{\Gamma }}{\left(1+\displaystyle \frac{1}{2}{\epsilon }_{1}{\epsilon }_{2}\right)}^{2/3}{\sin }^{2}\iota ,\end{eqnarray} \tag{ 26 }$$

with Γ = (₁m₂+₂m₁)/m. And in the expression of ${\tilde{h}}_{+,\times }^{\left(2\right)}$, we have adopted conventions similar to Van Den Broeck & Sengupta (2006), where ${e}^{-i{\varphi }_{\left(\mathrm{2,0}\right)}}{P}_{\left(\mathrm{2,0}\right)}=-(1+{\cos }^{2}\iota ){F}_{+}\,-i(2\cos \iota ){F}_{\times }$.

We also would like to investigate whether the constraints can be improved if the higher order amplitude corrections of the PN gravitational waveform are taken into consideration. In the stage of adiabatic inspiraling, the analytic waveforms can be obtained by using PN approximation where the waveforms can be expanded in terms of the orbital velocity. Thanks to great efforts over the past few decades, the PN waveforms have been calculated to very high orders and are sufficiently precise to extract small signals buried in the large noise by matched filtering in GW experiments. More details can be found in the review article by Blanchet (2014). Since the matched filtering method used in the GW detections is more sensitive to the phase of the templates than the amplitude, the restricted waveform is the most commonly used waveform model, in which only the dominant harmonic is taken into account, other than the leading order all amplitude corrections are discarded, but all the available orders of phases are included. However, some works (Van Den Broeck & Sengupta 2006; Arun et al. 2007; Trias & Sintes 2008a, 2008b) have shown that there can be considerable consequences if higher order amplitude corrections are included in the templates. Here we consider the full PN waveform in which amplitude terms are included up to 2.5 PN order and phase terms are included up to 3.5 PN order. The full waveforms are shown in Appendix A, where we adopt conventions similar to Van Den Broeck & Sengupta (2006).

In this work, we consider two proposed space-borne GW detectors, LISA and TianQin, to forecast the constraints on SMG. LISA is a mission led by the European Space Agency that can detect GWs in the milli-Hertz (0.1 mHz–1 Hz) range (Danzmann et al. 1996, 2016). LISA consists of three identical spacecraft that maintain an equilateral triangular configuration in an Earth-trailing heliocentric orbit between 50 and 65 million km from Earth. The distance between the two spacecraft is 2.5 million km according to the new LISA design (Danzmann et al. 2016). The line connecting the Sun and the center of mass of the detector keeps a 60° angle with respect to the plane of the constellation. In addition to the revolution around the Sun, the detector rotates clockwise (viewed from the Sun) around its center of mass with a period of one year. Pictures of this orbit configuration can be found in Figure 4.8 of Danzmann et al. (1996) or in Figure 4 of Danzmann et al. (2016). TianQin has a similar equilateral triangular configuration and is sensitive to the same frequency range. Different from LISA, TianQin is in a geocentric orbit with a period of 3.65 days. The distance between each pair of spacecraft is about 1.7 × 10⁵ km. The normal vector of the detector plane is fixed and points toward the reference source J0806.3+1527, which is a candidate ultracompact WD binary in the Galaxy (longitude = 1205, latitude = −47 in the ecliptic coordinate system) and is a strong periodic GW source in the milli-Hertz range. Illustrations of TianQin's configuration and orbit can be found in Figure 1 of Luo et al. (2016) or in Figures 1 and A1 of Hu et al. (2018). Two aspects of the detectors are relevant to our analysis, the noise spectrum and the antenna pattern functions, which will be introduced respectively in the following subsections.

Zoom In Zoom Out Reset image size

3.1. Antenna Beam-pattern Functions

In the proper detector frame, the response of a ground-based laser interferometer to GWs can be calculated using the equation of geodesic deviation, which is (Maggiore 2008; Poisson & Will 2014)

$$\begin{eqnarray}&&{\ddot{\xi }}^{i}=\displaystyle \frac{1}{2}{\ddot{h}}_{{ij}}{\xi }^{j}.\end{eqnarray} \tag{ 27 }$$

However, it is not completely correct to straightforwardly extend the same treatment to the situation of space-borne GW detectors. In the derivation of the equation of the geodesic deviation (27), the approximation that the distance between two test masses is much smaller than the typical scale over which the gravitational field changes significantly has be adopted. This means that Equation (27) can be used to derive the response function when the detector arm length is shorter than the reduced wavelength of the GWs (see Section 1.3, 9.1 of Maggiore 2008, and Cornish & Rubbo 2003; Hu et al. 2018 for more details). The corresponding frequency, which is called the transfer frequency, is given by

$$\begin{eqnarray}&&{f}_{* }=\displaystyle \frac{c}{2\pi L},\end{eqnarray} \tag{ 28 }$$

where L is the arm length of the detector. This condition is satisfied for the ground-based GW detectors because the sensitive bands of ground-based detectors are far below their transfer frequency, while it is not always satisfied for space-borne detectors. For instance, this critical frequency is 0.019 Hz for LISA, and it is 0.28 Hz for TianQin, which are similar to the sensitive frequency bands of the detectors. The transfer frequencies of LISA and TianQin are illustrated by vertical dotted blue and orange lines, respectively, in Figure 1. For science objectives such as supermassive black hole binaries, Equation (27) can be safely used (Klein et al. 2016; Feng et al. 2019). But in the case of BH–NS binaries considered in this work, the rough estimation of last stable frequency is much higher than the upper limit of the detector's sensitive band. It is not proper to use the response function derived by extending the approach that is used for ground-based interferometers.

The propagation of GWs in the time during which the photons travel from laser source to photodetector is neglected if one uses Equation (27) to derive the response function. When frequencies are higher than the transfer frequency, there may be GWs of a few wavelengths passing through the path of photons during the time between emission and reception of the photons, which makes the effect of GWs cancel out itself and deteriorates the response of detectors to the GWs. In order to get the exact response of detectors, the integration along the null geodesic of photons between two test masses should be calculated. The response functions of LISA-like detectors for two GR polarizations are given by Cornish & Larson (2001), Cornish & Rubbo (2003), and Rubbo et al. (2004), and the same process can be extended to other polarizations (Liang et al. 2019).

We present the general form of antenna pattern functions here, while the explicit expression and the details of the process can be found in Appendix B. The response of LISA or TianQin to GWs has been shown in Equation (16), where the antenna beam-pattern functions ${F}_{+,\times ,b,l}$ are given by

$$\begin{eqnarray}&&{F}_{A}=\displaystyle \frac{1}{2}{\epsilon }_{{ij}}^{A}\left[{\hat{l}}_{1}^{i}{\hat{l}}_{1}^{j}\,{ \mathcal T }(f,{\hat{{\boldsymbol{l}}}}_{1}\cdot \hat{{\boldsymbol{\Omega }}})-{\hat{l}}_{2}^{i}{\hat{l}}_{2}^{j}\,{ \mathcal T }(f,{\hat{{\boldsymbol{l}}}}_{2}\cdot \hat{{\boldsymbol{\Omega }}})\right],\end{eqnarray} \tag{ 29 }$$

where A denotes different polarizations (A = +, × , b, l), ${\epsilon }_{{ij}}^{A}$ denotes polarization tensors, and ${\hat{{\boldsymbol{l}}}}_{1}$ and ${\hat{{\boldsymbol{l}}}}_{2}$ denote the unit vectors of two arms. Comparing the antenna beam-pattern functions of ground-based detectors, we find the differences are quantified by the transfer functions ${ \mathcal T }(f,{\hat{{\boldsymbol{l}}}}_{1}\cdot \hat{{\boldsymbol{\Omega }}})$ and ${ \mathcal T }(f,{\hat{{\boldsymbol{l}}}}_{2}\cdot \hat{{\boldsymbol{\Omega }}})$, which are given by

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal T }(f,{\hat{{\boldsymbol{l}}}}_{1}\cdot \hat{{\boldsymbol{\Omega }}}) & = & \displaystyle \frac{1}{2}\left\{\mathrm{sinc}\left[\displaystyle \frac{f}{2{f}_{* }}(1-{\hat{{\boldsymbol{l}}}}_{1}\cdot \hat{{\boldsymbol{\Omega }}})\right]\right.\\ & & \times \,\exp \left[-i\displaystyle \frac{f}{2{f}_{* }}(3+{\hat{{\boldsymbol{l}}}}_{1}\cdot \hat{{\boldsymbol{\Omega }}})\right]\\ & & +\,\mathrm{sinc}\left[\displaystyle \frac{f}{2{f}_{* }}(1+{\hat{{\boldsymbol{l}}}}_{1}\cdot \hat{{\boldsymbol{\Omega }}})\right]\\ & & \times \,\left.\exp \left[-i\displaystyle \frac{f}{2{f}_{* }}(1+{\hat{{\boldsymbol{l}}}}_{1}\cdot \hat{{\boldsymbol{\Omega }}})\right]\right\},\end{array}\end{eqnarray} \tag{ 30 }$$

where $\mathrm{sinc}(x)=\tfrac{\sin x}{x}$, and $\hat{{\boldsymbol{\Omega }}}$ denotes the unit vector of the GW propagation direction. In the low-frequency limit f ≪ f_*, the transfer functions ${ \mathcal T }(f,{\hat{{\boldsymbol{l}}}}_{1}\cdot \hat{{\boldsymbol{\Omega }}})$ and ${ \mathcal T }(f,{\hat{{\boldsymbol{l}}}}_{2}\cdot \hat{{\boldsymbol{\Omega }}})$ approach unity, which returns to the case of ground-based detectors.

For the equilateral triangular configuration, there are two independent output signals. The second output signal, following the previous work (Cutler 1998), is equivalent to the response of a two-arm detector rotated by π/4 with respect to the first one, in the assumptions that the noise is Gaussian, stationary, and totally symmetric. As shown in Appendix B, (θ', φ', ψ') are employed to denote the GW source direction and the polarization angle in the detector coordinate system. In terms of (θ', φ', ψ'), the two output signals can be expressed as

$$\begin{eqnarray}\begin{array}{rcl}{h}^{I}(t) & = & \displaystyle \sum _{A}{F}_{A}(\theta ^{\prime} ,\varphi ^{\prime} ,\psi ^{\prime} ){h}_{A},\\ {h}^{{II}}(t) & = & \displaystyle \sum _{A}{F}_{A}(\theta ^{\prime} ,\varphi ^{\prime} -\pi /4,\psi ^{\prime} ){h}_{A},\end{array}\end{eqnarray} \tag{ 31 }$$

where A = (+, × , b, l).

3.2. Noise Spectra of GW Detectors

The noise of a GW detector can be characterized by the one-side noise power spectral density (PSD) S_n(h). We employ the noise curve of LISA from Belgacem et al. (2019):

$$\begin{eqnarray}&&{S}_{n}(f)=\displaystyle \frac{4{S}_{\mathrm{acc}}(f)+{S}_{\mathrm{other}}}{{L}_{\mathrm{LISA}}^{2}}\left[1+{\left(\displaystyle \frac{f}{1.29{f}_{* }}\right)}^{2}\right]+{S}_{\mathrm{conf}}(f).\end{eqnarray} \tag{ 32 }$$

Here, f_* = 0.019 Hz is the transfer frequency of LISA, and L_LISA is the arm length, which is 2.5 million km according to the new LISA design (Danzmann et al. 2016). The acceleration noise is given by

$$\begin{eqnarray}\begin{array}{rcl}{S}_{\mathrm{acc}}(f) & = & \displaystyle \frac{9\times {10}^{-30}\,{{\rm{m}}}^{2}\,{\mathrm{Hz}}^{3}}{{\left(2\pi f\right)}^{4}}\left[1+{\left(\displaystyle \frac{6\times {10}^{-4}\mathrm{Hz}}{f}\right)}^{2}\right.\\ & & \left.\times \,\left(1+{\left(\displaystyle \frac{2.22\times {10}^{-5}\mathrm{Hz}}{f}\right)}^{8}\right)\right].\end{array}\end{eqnarray} \tag{ 33 }$$

The other noise is given by

$$\begin{eqnarray}&&{S}_{\mathrm{other}}=8.899\times {10}^{-23}\,{{\rm{m}}}^{2}\,{\mathrm{Hz}}^{-1}.\end{eqnarray} \tag{ 34 }$$

In addition to the noise from instruments, the numerous compact WD binaries in the Galaxy can emit GWs of a few milli-Hertz and produce confusion noise. The confusion noise from unresolved binaries is approximated by

$$\begin{eqnarray}&&{S}_{\mathrm{conf}}(f)=\displaystyle \frac{A}{2}{e}^{-{s}_{1}{f}^{\alpha }}{f}^{-7/3}\left\{1-\tanh \left[{s}_{2}(f-\kappa )\right]\right\},\end{eqnarray} \tag{ 35 }$$

with A = (3/20)3.2665 × 10⁻⁴⁴ Hz^4/3, s₁ = 3014.3 Hz^−α, α = 1.183, s₂ = 2957.7 Hz⁻¹, and κ = 2.0928 × 10⁻³ Hz.

For TianQin, we employ the noise curve provided by Luo et al. (2016), Hu et al. (2018), and Feng et al. (2019),

$$\begin{eqnarray}\begin{array}{rcl}{S}_{n}(f) & = & \left[\displaystyle \frac{{S}_{x}}{{L}_{\mathrm{TianQin}}^{2}}+\displaystyle \frac{4{S}_{a}}{{\left(2\pi f\right)}^{4}{L}_{\mathrm{TianQin}}^{2}}\left(1+\displaystyle \frac{{10}^{-4}\,\mathrm{Hz}}{f}\right)\right]\\ & & \times \,\left[1+{\left(\displaystyle \frac{f}{1.29{f}_{* }}\right)}^{2}\right],\end{array}\end{eqnarray} \tag{ 36 }$$

where f_* = 0.28 Hz is the transfer frequency of TianQin, L_TianQin = 1.73 × 10⁸ m is the arm length, and S_x = 10⁻²⁴ m² Hz⁻¹ and S_a = 10⁻³⁰ m² s⁻⁴ Hz⁻¹ are the position noise and acceleration noise, respectively. The sensitivity curves $\sqrt{{S}_{n}(f)}$ of LISA and TianQin are presented in Figure 1.

4.1. Fisher Information Matrix

The Fisher matrix approach is widely used to estimate the precision of future experiments. Compared with techniques like Monte Carlo analysis, the Fisher matrix is a simpler way to efficiently estimate errors of parameters in GW detection with sufficient accuracy in the high signal-to-noise ratio (S/N) cases (Finn 1992; Finn & Chernoff 1993; Cutler & Flanagan 1994). The elements of a Fisher matrix are given by

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{{ij}}=\left\langle \displaystyle \frac{\partial \tilde{h}(f)}{\partial {p}_{i}},\displaystyle \frac{\partial \tilde{h}(f)}{\partial {p}_{i}}\right\rangle ,\end{eqnarray} \tag{ 37 }$$

where $\tilde{h}(f)$ is the Fourier transform of the output h(t) of the detectors, and p_i are the parameters to be estimated. The angle brackets denote the detector-dependent inner product,

$$\begin{eqnarray}&&\langle \tilde{a}(f),\tilde{b}(f)\rangle =4{\int }_{{f}_{1}}^{{f}_{2}}\displaystyle \frac{\tilde{a}(f){\tilde{b}}^{* }(f)+{\tilde{a}}^{* }(f)\tilde{b}(f)}{2}\displaystyle \frac{{df}}{{S}_{n}(f)},\end{eqnarray} \tag{ 38 }$$

where S_n(f) is the PSD of the detector. The upper limit of integral interval f₂ is determined by min(${f}_{\mathrm{LSO}}$, f_up), where f_up is the upper limit of the detector's sensitive band (1 Hz), and f_LSO is the last stable orbital frequency of the binary, which will be discussed in the Section (4.2). The lower limit of integral interval f₁ is given by max(f_low, f_obs). Here, f_low is the lower limit of the detector's sensitive band (10⁻⁴ Hz). The f_obs corresponds to the orbital frequency at T_obs earlier from the time corresponding to f₂. Approximately, T_obs can be regarded as the designed mission duration of the detector. For LISA, T_obs = 4 yr (Danzmann et al. 2016), and for TianQin, T_obs = 5 yr (Luo et al. 2016). Using the formula of orbital decay to leading order (Equation (43) with _d = 0), the relation between orbital frequency at the beginning and end of any time interval can be given by

$$\begin{eqnarray}&&{f}_{\mathrm{obs}}={f}_{2}{\left[1+\displaystyle \frac{256}{5}{\left({{GM}}_{c}\right)}^{\tfrac{5}{3}}{T}_{\mathrm{obs}}{\left(2\pi {f}_{2}\right)}^{\tfrac{8}{3}}\right]}^{-\tfrac{3}{8}}.\end{eqnarray} \tag{ 39 }$$

The Fisher matrix for the combination of the two independent output signals is given by

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{{ij}}={{\rm{\Lambda }}}_{{ij}}^{I}+{{\rm{\Lambda }}}_{{ij}}^{{II}},\end{eqnarray} \tag{ 40 }$$

where ${{\rm{\Lambda }}}_{{ij}}^{I}=\langle \tfrac{\partial {\tilde{h}}^{I}}{\partial {p}_{i}},\tfrac{\partial {\tilde{h}}^{I}}{\partial {p}_{i}}\rangle$, ${{\rm{\Lambda }}}_{{ij}}^{{II}}=\langle \tfrac{\partial {\tilde{h}}^{{II}}}{\partial {p}_{i}},\tfrac{\partial {\tilde{h}}^{{II}}}{\partial {p}_{i}}\rangle$. Using the definition of inner product, the combined S/N of two independent signals is

$$\begin{eqnarray}&&{\rho }^{2}={\left({\rho }^{I}\right)}^{2}+{\left({\rho }^{{II}}\right)}^{2},\end{eqnarray} \tag{ 41 }$$

with ${({\rho }^{I})}^{2}=\langle {\tilde{h}}^{I},{\tilde{h}}^{I}\rangle$ and ${({\rho }^{{II}})}^{2}=\langle {\tilde{h}}^{{II}},{\tilde{h}}^{{II}}\rangle$. The covariance matrix Σ can be derived by taking the inverse of Fisher matrix Λ, that is, Σ = Λ⁻¹, and the estimation of the rms error of a parameter p_i is given by ${\rm{\Delta }}{p}_{i}=\sqrt{{{\rm{\Sigma }}}_{{ii}}}$. The correlation coefficients between parameters p_i and p_j are given by

$$\begin{eqnarray}&&{c}_{{ij}}=\displaystyle \frac{{{\rm{\Sigma }}}_{{ij}}}{{\left({{\rm{\Sigma }}}_{{ii}}{{\rm{\Sigma }}}_{{jj}}\right)}^{1/2}}.\end{eqnarray} \tag{ 42 }$$

4.2. Constraining Screened Modified Gravity

As shown in Zhang et al. (2017) and Liu et al. (2018b), the orbital decay of a compact binary, due to gravitational radiation, is given by

$$\begin{eqnarray}&&\dot{\omega }(t)=\displaystyle \frac{96}{5}{\left({{GM}}_{c}\right)}^{\tfrac{5}{3}}{\omega }^{\tfrac{11}{3}}\left[1+\displaystyle \frac{5}{192}{\left({Gm}\omega \right)}^{-\tfrac{2}{3}}{\epsilon }_{d}^{2}\right],\end{eqnarray} \tag{ 43 }$$

where ω is the orbital angular frequency and _d = ₁ − ₂ is the difference between the screened parameters of two objects. We find that the asymmetric binary systems can induce more phase corrections, which induces the foremost difference between SMG and GR. Therefore, in this work, we consider the asymmetric systems, BH–NS binaries and BH–WD binaries, as the GW sources.

Because the stationary phase approximation, which is used to get the Fourier transform of a detector's response, is maintained only in the stage where the change in orbital frequency is negligible in the period of a single circle, and the PN waveforms are not accurate enough in the late stage, a specific frequency must be chosen where the waveforms are truncated. As rough estimations, we employ the Roche radius of a rigid spherical body as the last stable distance between the two objects (Scharre & Will 2002), which is given by

$$\begin{eqnarray}&&d={R}_{\mathrm{comp}}{\left(2\displaystyle \frac{{M}_{\mathrm{BH}}}{{M}_{\mathrm{comp}}}\right)}^{\tfrac{1}{3}}.\end{eqnarray} \tag{ 44 }$$

Here, ${R}_{\mathrm{comp}}$ denotes the radius of the companion, and M_BH and M_comp denote the mass of the black hole and the companion, respectively. The corresponding orbital frequency is given by

$$\begin{eqnarray}&&f=\displaystyle \frac{1}{2\pi }\sqrt{\displaystyle \frac{G({M}_{\mathrm{BH}}+{M}_{\mathrm{comp}})}{{d}^{3}}}.\end{eqnarray} \tag{ 45 }$$

For the NS, the typical values of mass and radius are M_NS = 1.4 M_⊙ and R_NS = 10 km. For the WD, we employ the parameters of the target in PSR J1738+0333, which are given by M_WD = 0.181 M_⊙ and R_NS = 0.037R_⊙ (Antoniadis et al. 2012). For the massive BH in the binary systems, we consider two cases in this paper, M_BH = 1000 M_⊙ and M_BH = 10,000 M_⊙. The last stable orbital frequencies of BH–WD binaries are approximately 0.0042 Hz for both binaries with different BH masses, which are illustrated by the vertical dashed green line in Figure 1. Since the f_LSO of BH–NS binaries exceeds the upper limit of the detectors' sensitive band, it is not shown in Figure 1. Similar to previous work (Will & Yunes 2004), the locations of the GW sources are set to the Virgo cluster, where the distance D = 16.5 Mpc and the celestial position in the ecliptic coordinate system is longitude 18104 and latitude 1433.

Because the screened parameter of BH is zero (Liu et al. 2018b), the terms related to ₁₂ vanish, and $| {\epsilon }_{d}|$ becomes the screened parameter of the companion _NS or _WD because the GW waveforms contain the term ${\epsilon }_{\mathrm{NS},\mathrm{WD}}^{2}$ rather than _NS,WD. For this reason, similar to Liu et al. (2018b) and Zhang et al. (2019b), we constrain the parameter ${\epsilon }_{\mathrm{NS},\ \mathrm{WD}}^{2}$ instead of _NS,WD in our analysis. In addition, because A₁ and _d always appear together in the GW waveforms, to evade a singular matrix in the practical computations of the Fisher matrix, we consider A_1d as a combination and constrain it, instead of the parameter A₁.

In the Fisher matrix analysis, there are 11 free parameters in the full response functions, which are

$$\begin{eqnarray}&&\left(\iota ,\mathrm{ln}{M}_{c},\mathrm{ln}\eta ,\mathrm{ln}r,{t}_{c},{{\rm{\Psi }}}_{c},\theta ,\varphi ,\psi ,{A}_{1}{M}_{\mathrm{Pl}}{\epsilon }_{d},{\epsilon }_{d}^{2}\right).\end{eqnarray} \tag{ 46 }$$

As mentioned above, the rms errors of the parameters can be estimated by the Fisher matrix method for future GW experiments. The GW waveforms return to those of GR when A₁ = 0 and _NS,WD = 0. So, in this article, we set A₁ = _NS,WD = 0 in the fiducial waveforms, and their rms errors can be considered as the upper limits of A₁ and _NS,WD by the potential observations. The values of the other parameters in fiducial waveforms are set as $\left({t}_{c}=0,{{\rm{\Psi }}}_{c}=0,\psi =0\right)$. The fiducial waveforms are presented in Figure 2, where we have set the inclination angle ι = 45° and the mass of BH is chosen to be 1000M_⊙.

Zoom In Zoom Out Reset image size

In our analysis, we consider three different cases to investigate the capabilities of LISA and TianQin:

• 1.  
  In the first case, we assume the GW detectors will constrain SMG only by observing the extra polarization modes of GWs. As mentioned in the previous discussion, it includes only the breathing polarization mode h_b, because we have set the mass of the scalar field ${m}_{s}\to 0$ in our calculation. The expansion coefficient of the coupling function A₁ is set to a specific value. The f(R) gravity can be cast into the form of scalar–tensor theory, and the scalar degree of freedom can be suppressed in high-density regions by the chameleon mechanism. The coupling function in f(R) gravity is given by $A(\phi )\,=\tfrac{1}{\sqrt{f^{\prime} (R)}}=\exp \left(\tfrac{\xi \phi }{{M}_{\mathrm{Pl}}}\right)$ with $\xi =1/\sqrt{6}$ (Liu et al. 2018a). Here, we choose ${A}_{1}{M}_{\mathrm{Pl}}=1/\sqrt{6}$ as a characteristic value. Thus, the amplitude of h_b is quantified by the screened parameter _NS,WD, which will be constrained by the potential observation. We consider the values of _NS,WD that can make the S/N reach 10 as the constraints of screened parameters of NSs or WDs.
• 2.  
  In the second case, we constrain the SMG by analyzing the deviation of the GW waveform in SMG from that in GR. The Fisher matrix technique will be employed for the analysis. Because the GW detectors are sensitive to the GW phases, rather than the amplitudes, in this case only the phase correction induced by dipole radiation is taken into consideration. We consider a waveform model including the restricted waveforms where all phase corrections are included while all amplitude corrections except the leading order are discarded and the phase correction induced by dipole radiation. As shown in Liu et al. (2018b), the standard PPE framework can be applied to the waveforms in SMG when we consider the two tensor polarizations h₊ and h_×. The general form of the detector's response function is given by

  $$\begin{eqnarray}&&\tilde{h}(f)={\tilde{h}}_{\mathrm{GR}}(f)\left[1+\alpha {\left(\pi {M}_{c}f\right)}^{\tfrac{a}{3}}\right]{e}^{i\beta {\left(\pi {M}_{c}f\right)}^{\tfrac{b}{3}}},\end{eqnarray} \tag{ 47 }$$

  where $\alpha ,a,\beta ,b$ are the four PPE parameters and ${\tilde{h}}_{\mathrm{GR}}(f)$ denotes the response function in GR. In this case, we only consider the non-GR correction in phase. The PPE parameters are taken to be $\alpha =a=0,\beta \,=-\tfrac{5}{14336}{\epsilon }_{d}^{2}{\eta }^{2/5},b=-7$, and the restricted waveforms are employed in ${\tilde{h}}_{\mathrm{GR}}(f)$. Since the terms related to A₁ only present in the amplitude corrections are discarded in this case, there are 10 parameters remaining in the Fisher matrix.
• 3.  
  In the third case, the higher order amplitude corrections of PN waveforms are included, and all available correction terms of GW waveforms are taken into account. The full response functions are presented in Appendix A. In order to investigate the influence of higher order amplitude corrections on the constraints of SMG, we use the waveform model that includes 3.5 PN phase corrections, 2.5 PN amplitude corrections, and corrections concerned with the SMG both in amplitude and phase as the input signals of the Fisher matrix. The 11 parameters in Equation (46) all exist in the Fisher matrix, and we can obtain the constraints of both A₁M_{Pl d} and ${\epsilon }_{\mathrm{NS},\mathrm{WD}}$ from the Fisher matrix analysis.

5.1. Constraints on the Screened Parameters of Neutron Star and White Dwarf

5.2. Comparison with Other Observational Constraints

In this section, we would like to compare the above results, which are forecasts for the future space-borne GW detectors, with the constraints placed by the present experiments, including pulsar timing observations, lunar laser ranging (LLR), and the Cassini experiment. We find that the constraint given by GW observations is complementary to the constraint from the Cassini experiment, but weaker than those from LLR and binary pulsars. Due to the strong surface gravitational potentials of neutron stars, although the screened parameter of a neutron star can be constrained quite well, the constraint of the scalar background ϕ_VEV is worse than that given by LLR. As for pulsar timing experiments of binary pulsars, the constraints on ϕ_VEV are actually from the constraints on the orbital period decay caused by energy loss through gravitational radiation. This is similar to GW observations that constrain ϕ_VEV by using the GW waveform. As we can see from the waveform (Equations (71), (72), and (82)), the terms concerned with SMG include the terms relevant to (2πfGm)⁻¹, ${(2\pi {fGm})}^{-\tfrac{1}{3}}$, ${(\pi {fGm})}^{-\tfrac{2}{3}}$ in amplitude (Equations (71), (72)) and ${(2\pi {fGm}/n)}^{-\tfrac{2}{3}}$ in phase (Equation (82)), which means that the SMG effects are more obvious in the lower frequency range. The most sensitive frequency of space-borne GW detectors is about 10⁻² Hz, while the orbital period of binary pulsars is on the order of 0.1 day. It is reasonable that the constraints given by pulsar timing are better than the constraints given by GW observations. This result implies that, at least for the three SMG models considered in this work, the GW observations by space-borne detectors may be not a good tool for the task of testing SMG theories.

Pulsar binary systems provide very useful tools for testing gravity theories. The first indirect evidence of the existence of GWs was given by the measurement of binary pulsar orbital period decay (Taylor & Weisberg 1982). By monitoring the orbital period change, the deviation from GR can be constrained. During the Apollo program and the Lunokhod missions, laser reflectors were installed on the moon. The laser pulses emitted on Earth can be reflected by the reflectors. By measuring the round-trip time, the Earth–Moon distance can be measured with extreme accuracy. The constraints on the Nordtvedt parameter and time variation of the gravitational constant can be given by LLR experiments (Hofmann et al. 2010). In this paper, we adopt the constraints in our previous work (Zhang et al. 2019a), which gave the upper bound on the scalar background ϕ_VEV (the VEV of the scalar field in SMG) as

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\phi }_{\mathrm{VEV}}}{{M}_{\mathrm{Pl}}}\right)}_{\mathrm{pulsar}}\leqslant 4.4\times {10}^{-8},\end{eqnarray} \tag{ 48 }$$

by pulsar observations of PSRs J1738+0333 and J0348+0432 at 95.4% confidence level (CL), and the constraints by LLR at 95.4% CL of

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\phi }_{\mathrm{VEV}}}{{M}_{\mathrm{Pl}}}\right)}_{\mathrm{LLR}}\leqslant 7.8\times {10}^{-15}.\end{eqnarray} \tag{ 49 }$$

The Cassini satellite was in solar conjunction in 2002. The Shapiro time-delay measurements using the Cassini spacecraft yielded a very tight constraint on the PPN parameter γ (Bertotti et al. 2003):

$$\begin{eqnarray}&&| {\gamma }_{\mathrm{obs}}-1| \leqslant 2.3\times {10}^{-5}.\end{eqnarray} \tag{ 50 }$$

These constraints will be compared with the potential constraint from future GW observations in this subsection.

In SMG, the screened parameter of an NS or WD can be approximated by (Zhang et al. 2017)

$$\begin{eqnarray}&&{\epsilon }_{a}=\displaystyle \frac{{\phi }_{\mathrm{VEV}}}{{M}_{\mathrm{Pl}}{{\rm{\Phi }}}_{a}},\end{eqnarray} \tag{ 51 }$$

where a denotes NS or WD, ${{\rm{\Phi }}}_{a}={{Gm}}_{a}/{R}_{a}$ is the surface gravitational potential of the a object, and ϕ_VEV is the scalar background in SMG. The constraints on the screened parameter _a can be converted to the constraints on the scalar background ϕ_VEV, and vice versa. Here, we consider the best constraint on the screened parameter given by TianQin, which is _NS ≤ 3.3 × 10⁻⁵, and compare with other observational constraints on SMG. The corresponding constraint on ϕ_VEV given by TianQin is

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\phi }_{\mathrm{VEV}}}{{M}_{\mathrm{Pl}}}\right)}_{\mathrm{GW}(\mathrm{NS})}\leqslant 7.6\times {10}^{-6}.\end{eqnarray} \tag{ 52 }$$

Similarly, we also consider the best result of A₁,

$$\begin{eqnarray}&&{A}_{1}{M}_{\mathrm{Pl}}{\epsilon }_{\mathrm{NS}}\leqslant 2.6\times {10}^{-3},\end{eqnarray} \tag{ 53 }$$

as a typical value of constraint on A₁ to compare with other constraints. In the following, we will compare these constraints in three specific SMG models: chameleon, symmetron, and dilaton theories.

5.2.1. Chameleon

The chameleon model was proposed by Khoury & Weltman (2004a, 2004b), who introduced the screening mechanism by making the mass of the scalar field depend on the environment density. The original chameleon model has been ruled out by the combined constraints from the solar system and cosmology (Hees & Füzfa 2012; Zhang et al. 2016). The idea of the chameleon can be revived by introducing a potential and coupling function that has an exponential form (Brax et al. 2004)

$$\begin{eqnarray}&&V(\phi )={{\rm{\Lambda }}}^{4}\exp \left(\displaystyle \frac{{{\rm{\Lambda }}}^{4}\alpha }{{\phi }^{\alpha }}\right),\,\,\,A(\phi )=\exp \left(\displaystyle \frac{\beta \phi }{{M}_{\mathrm{Pl}}}\right).\end{eqnarray} \tag{ 54 }$$

Here, α, β are positive dimensionless constants, and Λ denotes the energy scale of the theory, which is required by the cosmological constraints to be close to the dark energy scale $2.24\times {10}^{-3}$ eV (Hamilton et al. 2015; Zhang et al. 2016). The scalar background ϕ_VEV in the chameleon model is given by (Zhang et al. 2016, 2017)

$$\begin{eqnarray}&&{\phi }_{\mathrm{VEV}}={\left[\displaystyle \frac{\alpha {M}_{\mathrm{Pl}}{{\rm{\Lambda }}}^{4+\alpha }}{\beta {\rho }_{b}}\right]}^{\tfrac{1}{1+\alpha }},\end{eqnarray} \tag{ 55 }$$

where ρ_b is the background matter density corresponding to the galactic matter density ρ_gal ≃ 10⁻⁴² GeV⁴ (Zhang et al. 2017). The PPN parameter γ in the chameleon model is given by (Zhang et al. 2016)

$$\begin{eqnarray}&&\gamma =1-\displaystyle \frac{2\beta {\phi }_{\mathrm{VEV}}}{{M}_{\mathrm{Pl}}{{\rm{\Phi }}}_{\mathrm{Sun}}},\end{eqnarray} \tag{ 56 }$$

where ${{\rm{\Phi }}}_{\mathrm{Sun}}$ denotes the surface gravitational potential of the Sun. The expansion coefficient of coupling function A₁ is given by Zhang et al. (2016), which is ${A}_{1}/{A}_{0}=\beta /{M}_{\mathrm{Pl}}$. Recall that what we actually get by the process discussed above is the constraints on A₁M_PlNS, which takes the form

$$\begin{eqnarray}&&{A}_{1}{M}_{\mathrm{Pl}}{\epsilon }_{\mathrm{NS}}=\beta \displaystyle \frac{{\phi }_{\mathrm{VEV}}}{{M}_{\mathrm{Pl}}{{\rm{\Phi }}}_{\mathrm{NS}}},\end{eqnarray} \tag{ 57 }$$

where we have adopted ${A}_{0}=1$. Using the above formulae (Equations (55)–(57)), the constraints on the parameters of the chameleon model can be obtained by the constraints on ϕ_VEV, A₁, and γ (Equations (48)–(50), (52), (53)). We find that the expression of A₁ in Equation (57) is similar to that of the PPN parameter γ in Equation (56). The comparison can be glimpsed by comparing β ϕ_VEV/M_Pl, which is

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\beta {\phi }_{\mathrm{VEV}}}{{M}_{\mathrm{Pl}}}\right)}_{{A}_{1}}\leqslant 5.4\times {10}^{-4},\,\,{\left(\displaystyle \frac{\beta {\phi }_{\mathrm{VEV}}}{{M}_{\mathrm{Pl}}}\right)}_{\mathrm{Cassini}}\leqslant 2.4\times {10}^{-11}.\end{eqnarray} \tag{ 58 }$$

Since the constraint given by A₁ is much weaker than that derived from other observations, the allowed range will fill the full region of Figure 3. For this reason, the constraint corresponding to A₁ is not shown in Figure 3.

Zoom In Zoom Out Reset image size

The other four constraints are illustrated in Figure 3, where the dashed line denotes the forecast for the GW constraint, the solid lines denote the constraints of real experiments (pulsar and LLR), their allowed regions the areas to the right of the corresponding lines, and the region allowed by the Cassini experiment is illustrated by the yellow area. Although the GW observations can give the tight constraint on the screened parameter of the NS, the constraint on the scalar background ϕ_VEV cannot be improved simultaneously because the surface gravitational potential of NS is much larger than that of a WD or the solar system. We find that the most stringent bound on chameleon is still given by the combined constraint of LLR and Cassini (Zhang et al. 2019a), which gives α ≥ 0.35.

5.2.2. Symmetron

The symmetron models are characterized by a Mexican hat potential and a quadratic coupling function (Hinterbichler & Khoury 2010; Hinterbichler et al. 2011; Davis et al. 2012),

$$\begin{eqnarray}&&V(\phi )={V}_{0}-\displaystyle \frac{1}{2}{\mu }^{2}{\phi }^{2}+\displaystyle \frac{\lambda }{4}{\phi }^{4},\,\,\,A(\phi )=1+\displaystyle \frac{{\phi }^{2}}{2{M}^{2}},\end{eqnarray} \tag{ 59 }$$

where μ and M are mass scales, λ is a positive dimensionless coupling constant, and V₀ is the vacuum energy of the bare potential V(ϕ). In the symmetron model, the VEV ϕ_VEV is given by (Zhang et al. 2016, 2017)

$$\begin{eqnarray}&&{\phi }_{\mathrm{VEV}}=\displaystyle \frac{{m}_{s}}{\sqrt{2\lambda }},\end{eqnarray} \tag{ 60 }$$

which is proportional to the scalar mass. Similar to the chameleon model, the constraints on the scalar background ϕ_VEV can be interpreted as the constraints on the parameters m_s and λ of the symmetron model, which are shown in the upper panel of Figure 4. The m_s is the effective mass of the scalar field background. The scalar field background plays the role of dark energy, which should have effects on large scales to accelerate the expansion of the universe. So the ${m}_{s}^{-1}$ is considered as roughly cosmological scale (∼1 Mpc; Zhang et al. 2017).

Zoom In Zoom Out Reset image size

The PPN parameter γ in the symmetron model is (Zhang et al. 2017, 2016)

$$\begin{eqnarray}&&\gamma =1-2\displaystyle \frac{{\phi }_{\mathrm{VEV}}^{2}}{{M}^{2}{{\rm{\Phi }}}_{\mathrm{Sun}}}.\end{eqnarray} \tag{ 61 }$$

Thus, we can obtain the constraint on the scalar background ϕ_VEV with mass scale M from the Cassini experiment, which is presented in the bottom panel of Figure 4. The expansion coefficient of coupling function A₁ is given by (Zhang et al. 2016)

$$\begin{eqnarray}&&{A}_{1}=\displaystyle \frac{{\phi }_{\mathrm{VEV}}}{{M}^{2}}.\end{eqnarray} \tag{ 62 }$$

The term A₁M_{Pl NS}, which is treated as a parameter in the computations, takes the form

$$\begin{eqnarray}&&{A}_{1}{M}_{\mathrm{Pl}}{\epsilon }_{\mathrm{NS}}=\displaystyle \frac{{\phi }_{\mathrm{VEV}}^{2}}{{M}^{2}{{\rm{\Phi }}}_{\mathrm{NS}}}.\end{eqnarray} \tag{ 63 }$$

Thus, we can compare the constraints of A₁ with the constraints of the Cassini experiment by comparing ${\phi }_{\mathrm{VEV}}^{2}/{M}^{2}$ directly:

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\phi }_{\mathrm{VEV}}^{2}}{{M}^{2}}\right)}_{{A}_{1}}\leqslant 5.4\times {10}^{-4},\,\,\,{\left(\displaystyle \frac{{\phi }_{\mathrm{VEV}}^{2}}{{M}^{2}}\right)}_{{Cassini}}\leqslant 2.4\times {10}^{-11}.\end{eqnarray} \tag{ 64 }$$

Since the surface gravitational potential of NS is much larger than that of the Sun, the constraint given by A₁ is much weaker than that given by the Cassini experiment. The constraint given by A₁ is also not plotted in Figure 4.

The other four constraints are illustrated in Figure 4. The upper panel shows the upper bound on the mass of scalar field m_s (or the lower bound on ${m}_{s}^{-1}$) with the coupling constant λ. The dashed line denotes the forecast for GW constraints, and the solid lines denote the constraints of real experiments (pulsar and LLR). The most stringent constraint is given by LLR. If ${m}_{s}^{-1}\leqslant 1$ Mpc, the constraint on λ given by LLR is λ ≥ 10^−85.3. The bottom panel shows the constraints in parameter space (ϕ_VEV, M). The yellow area denotes the allowed region of constraint given by the Cassini experiment. The dashed and solid vertical lines represent the constraints on ϕ_VEV given by the forecast of GWs (Equation (52)) and the real experiments (Equations (49) and (48)), respectively. The corresponding allowed regions are the areas to the left of the lines. The most stringent constraint is still given by the combined constraints of LLR and Cassini.

5.2.3. Dilaton

In the dilaton model, the potential and coupling functions take the forms (Damour & Polyakov 1994a, 1994b; Brax et al. 2010)

$$\begin{eqnarray}&&V(\phi )={V}_{0}\exp (-\displaystyle \frac{\phi }{{M}_{\mathrm{Pl}}}),\,\,\,A(\phi )=1+\displaystyle \frac{{\left(\phi -{\phi }_{\star }\right)}^{2}}{2{M}^{2}},\end{eqnarray} \tag{ 65 }$$

where V₀ is a constant that has the dimension of energy density, M denotes the energy scale of the theory, and ϕ_⋆ represents the approximate value of the scalar field today. The scalar background ϕ_VEV and PPN parameter γ are given by (Zhang et al. 2016, 2017)

$$\begin{eqnarray}&&{\phi }_{\mathrm{VEV}}={\phi }_{\star }+\displaystyle \frac{{M}^{2}{\rho }_{{{\rm{\Lambda }}}_{0}}}{{M}_{\mathrm{Pl}}{\rho }_{b}},\,\,\,\gamma =1-2\displaystyle \frac{{\left({\phi }_{\mathrm{VEV}}-{\phi }_{\star }\right)}^{2}}{{M}^{2}{{\rm{\Phi }}}_{\mathrm{Sun}}}.\end{eqnarray} \tag{ 66 }$$

Here, ρ_b is the background matter density, which is the galactic matter density ρ_gal ≃ 10⁻⁴² GeV⁴ in our calculation, and ${\rho }_{{{\rm{\Lambda }}}_{0}}$ denotes the density of dark energy, which is ${\rho }_{{{\rm{\Lambda }}}_{0}}\simeq 2.51\,\times {10}^{-47}\,{\mathrm{GeV}}^{4}$ (Zhang et al. 2016). The screened parameter of a object takes the form ${\epsilon }_{a}=\tfrac{{\phi }_{\mathrm{VEV}}-{\phi }_{a}}{{M}_{\mathrm{Pl}}{{\rm{\Phi }}}_{a}}$, where ${\phi }_{a}={\phi }_{\star }+\tfrac{{M}^{2}{\rho }_{{{\rm{\Lambda }}}_{0}}}{{M}_{\mathrm{Pl}}{\rho }_{a}}$ is the minimum of the effective potential inside the object, and ρ_a is the matter density inside the object. Because the matter density in compact objects is much larger than that in the cosmological background, we can drop the term $\tfrac{{M}^{2}{\rho }_{{{\rm{\Lambda }}}_{0}}}{{M}_{\mathrm{Pl}}{\rho }_{a}}$ in the relation between the screened parameter _a and the parameter M of the dilaton model, which has the form

$$\begin{eqnarray}&&{\epsilon }_{a}{{\rm{\Phi }}}_{a}=\displaystyle \frac{{M}^{2}}{{M}_{\mathrm{Pl}}^{2}}\displaystyle \frac{{\rho }_{{{\rm{\Lambda }}}_{0}}}{{\rho }_{b}}.\end{eqnarray} \tag{ 67 }$$

Another parameter for which we can get the constraints from GW observation takes the form

$$\begin{eqnarray}&&{A}_{1}{M}_{\mathrm{Pl}}{\epsilon }_{\mathrm{NS}}={\left(\displaystyle \frac{M}{{M}_{\mathrm{Pl}}}\right)}^{2}{\left(\displaystyle \frac{{\rho }_{{{\rm{\Lambda }}}_{0}}}{{\rho }_{b}}\right)}^{2}\displaystyle \frac{1}{{{\rm{\Phi }}}_{\mathrm{NS}}}.\end{eqnarray} \tag{ 68 }$$

As in the previous models, the constraints on screened parameters _a, A₁, and PPN parameter γ can be switched to the constraints on the parameter M of the dilaton model. The results are shown in Table 5. As mentioned before, although the future observations of GWs from BH–NS binaries can constrain the screened parameter of NSs very stringently, the constraint on the mass scale M by GW observations is not stringent, due to the large gravitational surface potential of NS. The tightest constraint on the dilaton model is still given by LLR observation.

Table 5.  Constraints of the Dilaton Model Derived from Various Observations

	GW(NS)	GW(A1)	Pulsar	LLR	Cassini
M/MPl	≤0.55	≤920	≤0.042	≤1.8 × 10−5	≤0.20

Download table as:  ASCIITypeset image

Gravitational waves provide an excellent opportunity to test GR, which is always considered as the most successful theory of gravity, in the strong gravitational fields. In this issue, the calculation of GW waveforms in alternative gravitational theory is important. The SMG is one of the simplest extensions of GR in the scalar–tensor framework, which naturally explains the acceleration of cosmic expansion by introducing the scalar field. In addition, in this theory, the fifth force caused by the scalar field can be suppressed in the dense regions to satisfy various tests in the solar system and laboratories. For these reasons, the SMG theory and its specific models, including chameleon, symmetron, dilaton, and f(R), have been widely studied in the literature. Based on the GW waveforms produced by the coalescence of compact binaries in general SMG derived in Liu et al. (2018b), in this article we investigate the potential constraints on the general SMG theory from future GW observations. In our calculations, we focus on future space-borne missions, including LISA, TianQin, and Taiji, and assume EMRIs, including BH–NS and BH–WD in the Virgo cluster, as the GW targets. By comparing three different cases, we find that the extra polarization modes, the breathing mode and the longitude mode, have little contribution to the constraining of model parameters. The modifications of GW waveforms in the plus and cross modes, in particular the correction terms in the GW phases, dominate the constraint of SMG parameters. If a GW signal produced by the coalescence of a BH–NS system is detected by LISA, Taiji, or TianQin, the screened parameter _NS can be constrained at the level of $\lt { \mathcal O }({10}^{-5})$. On the other hand, limited by the durations and the sensitive frequency bands of the GW detectors, we find that the GW signals produced by the coalescence of BH–WD systems are difficult to detect by LISA, Taiji, or TianQin. For three specific SMG models (chameleon, symmetron, and dilaton), we compare this potential constraint with the other existing constraints derived by the Cassini experiment, LLR observations, and binary pulsars. We find that the constraint from GW observation is complementary to that from the Cassini experiment, but weaker than those from LLR observations and binary pulsars.

This work is supported by NSFC grants No. 11773028, No. 11633001, No. 11653002, No. 11421303, No. 11903030, No. 11903033, the Fundamental Research Funds for the Central Universities, the China Postdoctoral Science Foundation grant No. 2019M652193, and the Strategic Priority Research Program of the Chinese Academy of Sciences grant No. XDB23010200.

As shown in previous works (Van Den Broeck & Sengupta 2006; Arun et al. 2007; Trias & Sintes 2008a, 2008b), for space-borne detectors, considerable consequences can result if higher order amplitude corrections are taken into consideration. In order to investigate whether including higher PN order amplitude corrections can affect the constraints on SMG or not, we adopt a waveform that combines the SMG corrections which were derived in our previous work (Liu et al. 2018b) and the full GR waveform, which is up to 3.5 PN in the phase and 2.5 PN in the amplitude. In this appendix, we present explicitly this waveform in which the corrections caused by SMG and the Doppler modulation peculiar to space-borne detectors are included. We adopt conventions similar to Van Den Broeck & Sengupta (2006), where the full waveforms in GR have been presented.

In the PN approximation, the waveforms can be expressed as expansions in the typical internal speed of the source. The general forms of two GR polarizations are written as

$$\begin{eqnarray}\begin{array}{rcl}{h}_{+,\times }(t) & = & \displaystyle \frac{2m\eta }{D}x\left\{{H}_{+,\times }^{\left(0\right)}+{x}^{1/2}{H}_{+,\times }^{\left(1/2\right)}+{x}^{1}{H}_{+,\times }^{\left(1\right)}\right.\\ & & +\,\left.{x}^{3/2}{H}_{+,\times }^{\left(3/2\right)}+{x}^{2}{H}_{+,\times }^{\left(2\right)}+{x}^{5/2}{H}_{+,\times }^{\left(5/2\right)}\right\},\end{array}\end{eqnarray} \tag{ 69 }$$

where x is the expansion parameter, which is defined as x = [2πmF(t)]^2/3 with F(t) the orbital frequency. The expansion coefficients ${H}_{+,\times }^{(s)}$ consist of linear combinations of $\cos [n{\rm{\Psi }}(t)]$ and $\sin [n{\rm{\Psi }}(t)]$, where Ψ(t) is the orbital phase and the number of harmonics n = 7 for 2.5 PN order in amplitude. The explicit expressions of ${H}_{+,\times }^{(s)}$ can be found in Arun et al. (2004, 2005). As mentioned in Equation (16), the response function depends not only on the waveforms but also on the antenna pattern functions. The analytic expressions for the Fourier transform of a detector's response function can be obtained by using a stationary phase approximation. The expressions are the sum of seven harmonics:

$$\begin{eqnarray}&&\tilde{h}(f)=\sum _{k=1}^{7}{\tilde{h}}^{\left(k\right)}(f).\end{eqnarray} \tag{ 70 }$$

Taking into account the corrections caused by SMG, the explicit expressions of ${\tilde{h}}^{(k)}(f)$ are as follows:

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{h}}^{\left(1\right)}(f) & = & {\left(\displaystyle \frac{5}{48}\right)}^{\tfrac{1}{2}}{\pi }^{-\tfrac{2}{3}}\displaystyle \frac{{\left({{GM}}_{c}\right)}^{5/6}}{D}{\left(2f\right)}^{-\tfrac{7}{6}}\left\{-\displaystyle \frac{5}{384}{\pi }^{\tfrac{1}{6}}E{\epsilon }_{d}^{2}{\left(2\pi {fGm}\right)}^{-1}\right.\\ & & +\,{\pi }^{\tfrac{1}{6}}E{\left(2\pi {fGm}\right)}^{-\tfrac{1}{3}}\\ & & +\,{e}^{-i{\varphi }_{\left(\mathrm{1,1}/2\right)}}{P}_{\left(\mathrm{1,1}/2\right)}{\left(2\pi {fGm}\right)}^{\tfrac{1}{3}}\\ & & +\,\left[{e}^{-i{\varphi }_{\left(\mathrm{1,3}/2\right)}}{P}_{\left(\mathrm{1,3}/2\right)}+{e}^{-i{\varphi }_{\left(\mathrm{1,1}/2\right)}}{P}_{\left(\mathrm{1,1}/2\right)}{S}_{1}\right](2\pi {fGm})\\ & & +\,\left[{e}^{-i{\varphi }_{\left(\mathrm{1,2}\right)}}{P}_{\left(\mathrm{1,2}\right)}+{e}^{-i{\varphi }_{\left(\mathrm{1,1}/2\right)}}{P}_{\left(\mathrm{1,1}/2\right)}{S}_{3/2}\right]{\left(2\pi {fGm}\right)}^{\tfrac{4}{3}}\\ & & +\,\left[{e}^{-i{\varphi }_{\left(\mathrm{1,5}/2\right)}}{P}_{\left(\mathrm{1,5}/2\right)}+{e}^{-i{\varphi }_{\left(\mathrm{1,3}/2\right)}}{P}_{\left(\mathrm{1,3}/2\right)}{S}_{1}\right.\\ & & \left.\left.+{e}^{-i{\varphi }_{\left(\mathrm{1,1}/2\right)}}{P}_{\left(\mathrm{1,1}/2\right)}{S}_{2}\right]{\left(2\pi {fGm}\right)}^{\tfrac{5}{3}},\}\right\}\\ & & \times \,{\rm{\Theta }}({f}_{\mathrm{LSO}}-f)\exp \left\{i\left[2\pi {{ft}}_{c}-\displaystyle \frac{\pi }{4}+{\rm{\Psi }}(f)\right]\right\},\end{array}\end{eqnarray} \tag{ 71 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{h}}^{\left(2\right)}(f) & = & {2}^{-\tfrac{1}{2}}{\left(\displaystyle \frac{5}{48}\right)}^{\tfrac{1}{2}}{\pi }^{-\tfrac{2}{3}}\displaystyle \frac{{\left({{GM}}_{c}\right)}^{5/6}}{D}{\left(f\right)}^{-\tfrac{7}{6}}\\ & & \times \,\left\{\left[{{TF}}_{b}{S}_{-1}+{{QS}}_{-1}{e}^{-i{\varphi }_{\left(\mathrm{2,0}\right)}}{P}_{\left(\mathrm{2,0}\right)}\right]{\left(\pi {fGm}\right)}^{-\tfrac{2}{3}}\right.\\ & & +\,\left[{{TF}}_{b}+{{Qe}}^{-i{\varphi }_{\left(\mathrm{2,0}\right)}}{P}_{\left(\mathrm{2,0}\right)}\right]\\ & & +\,\left[{e}^{-i{\varphi }_{\left(\mathrm{2,1}\right)}}{P}_{\left(\mathrm{2,1}\right)}+{e}^{-i{\varphi }_{\left(\mathrm{2,0}\right)}}{P}_{\left(\mathrm{2,0}\right)}{S}_{1}\right]{\left(\pi {fGm}\right)}^{\tfrac{2}{3}}\\ & & +\,\left[{e}^{-i{\varphi }_{\left(\mathrm{2,3}/2\right)}}{P}_{\left(\mathrm{2,3}/2\right)}+{e}^{-i{\varphi }_{\left(\mathrm{2,0}\right)}}{P}_{\left(\mathrm{2,0}\right)}{S}_{3/2}\right](\pi {fGm})\\ & & +\,\left[{e}^{-i{\varphi }_{\left(\mathrm{2,2}\right)}}{P}_{\left(\mathrm{2,2}\right)}+{e}^{-i{\varphi }_{\left(\mathrm{2,1}\right)}}{P}_{\left(\mathrm{2,1}\right)}{S}_{1}+{e}^{-i{\varphi }_{\left(\mathrm{2,0}\right)}}\right.\\ & & \left.\times \,{P}_{\left(\mathrm{2,0}\right)}{S}_{2}\right]{\left(\pi {fGm}\right)}^{\tfrac{4}{3}}\\ & & +\,\left[{e}^{-i{\varphi }_{\left(\mathrm{2,5}/2\right)}}{P}_{\left(\mathrm{2,5}/2\right)}+{e}^{-i{\varphi }_{\left(\mathrm{2,3}/2\right)}}{P}_{\left(\mathrm{2,3}/2\right)}{S}_{1}+{e}^{-i{\varphi }_{\left(\mathrm{2,1}\right)}}\right.\\ & & \times \,\left.\left.{P}_{\left(\mathrm{2,1}\right)}{S}_{3/2}+{e}^{-i{\varphi }_{\left(\mathrm{2,0}\right)}}{P}_{\left(\mathrm{2,0}\right)}{S}_{5/2}\right]{\left(\pi {fGm}\right)}^{\tfrac{5}{3}}\right\}\\ & & \times \,{\rm{\Theta }}(2{f}_{\mathrm{LSO}}-f)\\ & & \times \,\exp \left\{i\left[2\pi {{ft}}_{c}-\displaystyle \frac{\pi }{4}+2{\rm{\Psi }}(f/2)\right]\right\},\end{array}\end{eqnarray} \tag{ 72 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{h}}^{\left(3\right)}(f) & = & {3}^{-\tfrac{1}{2}}{\left(\displaystyle \frac{5}{48}\right)}^{\tfrac{1}{2}}{\pi }^{-\tfrac{2}{3}}\displaystyle \frac{{\left({{GM}}_{c}\right)}^{5/6}}{D}{\left(2f/3\right)}^{-\tfrac{7}{6}}\\ & & \times \,\left\{{e}^{-i{\varphi }_{\left(\mathrm{3,1}/2\right)}}{P}_{\left(\mathrm{3,1}/2\right)}{\left(2\pi {fGm}/3\right)}^{\tfrac{1}{3}}\right.\\ & & +\,\left[{e}^{-i{\varphi }_{\left(\mathrm{3,3}/2\right)}}{P}_{\left(\mathrm{3,3}/2\right)}+{e}^{-i{\varphi }_{\left(\mathrm{3,1}/2\right)}}{P}_{\left(\mathrm{3,1}/2\right)}{S}_{1}\right]\left(2\pi {fGm}/3\right)\\ & & +\,\left[{e}^{-i{\varphi }_{\left(\mathrm{3,2}\right)}}{P}_{\left(\mathrm{3,2}\right)}+{e}^{-i{\varphi }_{\left(\mathrm{3,1}/2\right)}}{P}_{\left(\mathrm{3,1}/2\right)}{S}_{3/2}\right]{\left(2\pi {fGm}/3\right)}^{\tfrac{4}{3}}\\ & & +\,\left[{e}^{-i{\varphi }_{\left(\mathrm{3,3}/2\right)}}{P}_{\left(\mathrm{3,3}/2\right)}{S}_{1}\right.\\ & & \left.\left.+{e}^{-i{\varphi }_{\left(\mathrm{3,1}/2\right)}}{P}_{\left(\mathrm{3,1}/2\right)}{S}_{2}\right]{\left(2\pi {fGm}/3\right)}^{\tfrac{5}{3}}\right\}\\ & & \times \,{\rm{\Theta }}(3{f}_{\mathrm{LSO}}-f)\exp \\ & & \times \,\left\{i\left[2\pi {{ft}}_{c}-\displaystyle \frac{\pi }{4}+3{\rm{\Psi }}(f/3)\right]\right\},\end{array}\end{eqnarray} \tag{ 73 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{h}}^{\left(4\right)}(f) & = & {4}^{-\tfrac{1}{2}}{\left(\displaystyle \frac{5}{48}\right)}^{\tfrac{1}{2}}{\pi }^{-\tfrac{2}{3}}\displaystyle \frac{{\left({{GM}}_{c}\right)}^{5/6}}{D}{\left(f/2\right)}^{-\tfrac{7}{6}}\\ & & \times \,\left\{{e}^{-i{\varphi }_{\left(\mathrm{4,1}\right)}}{P}_{\left(\mathrm{4,1}\right)}{\left(\pi {fGm}/2\right)}^{\tfrac{2}{3}}\right.\\ & & +\,\left[{e}^{-i{\varphi }_{\left(\mathrm{4,2}\right)}}{P}_{\left(\mathrm{4,2}\right)}+{e}^{-i{\varphi }_{\left(\mathrm{4,1}\right)}}{P}_{\left(\mathrm{4,1}\right)}{S}_{1}\right]{\left(\pi {fGm}/2\right)}^{\tfrac{4}{3}}\\ & & +\,\left[{e}^{-i{\varphi }_{\left(\mathrm{4,5}/2\right)}}{P}_{\left(\mathrm{4,5}/2\right)}+{e}^{-i{\varphi }_{\left(\mathrm{4,1}\right)}}{P}_{\left(\mathrm{4,1}\right)}{S}_{3/2}\right]\\ & & \times \,\left.{\left(\pi {fGm}/2\right)}^{\tfrac{5}{3}}\right\}\\ & & \times \,{\rm{\Theta }}\left(4{f}_{\mathrm{LSO}}-f\right)\\ & & \times \,\exp \left\{i\left[2\pi {{ft}}_{c}-\displaystyle \frac{\pi }{4}+4{\rm{\Psi }}\left(f/4\right)\right]\right\},\end{array}\end{eqnarray} \tag{ 74 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{h}}^{\left(5\right)}(f) & = & {5}^{-\tfrac{1}{2}}{\left(\displaystyle \frac{5}{48}\right)}^{\tfrac{1}{2}}{\pi }^{-\tfrac{2}{3}}\displaystyle \frac{{\left({{GM}}_{c}\right)}^{5/6}}{D}{\left(2f/5\right)}^{-\tfrac{7}{6}}\\ & & \times \,\left\{{e}^{-i{\varphi }_{\left(\mathrm{5,3}/2\right)}}{P}_{\left(\mathrm{5,3}/2\right)}\left(2\pi {fGm}/5\right)\right.\\ & & +\,\left[{e}^{-i{\varphi }_{\left(\mathrm{5,5}/2\right)}}{P}_{\left(\mathrm{5,5}/2\right)}+{e}^{-i{\varphi }_{\left(\mathrm{5,3}/2\right)}}{P}_{\left(\mathrm{5,3}/2\right)}{S}_{1}\right]\\ & & \times \,\left.{\left(2\pi {fGm}/5\right)}^{\tfrac{5}{3}}\right\}\\ & & \times \,{\rm{\Theta }}\left(5{f}_{\mathrm{LSO}}-f\right)\\ & & \times \,\exp \left\{i\left[2\pi {{ft}}_{c}-\displaystyle \frac{\pi }{4}+5{\rm{\Psi }}\left(f/5\right)\right]\right\},\end{array}\end{eqnarray} \tag{ 75 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{h}}^{\left(6\right)}\left(f\right) & = & {6}^{-\tfrac{1}{2}}{\left(\displaystyle \frac{5}{48}\right)}^{\tfrac{1}{2}}{\pi }^{-\tfrac{2}{3}}\displaystyle \frac{{\left({{GM}}_{c}\right)}^{5/6}}{D}{\left(f/3\right)}^{-\tfrac{7}{6}}\\ & & \times \,\left\{{e}^{-i{\varphi }_{\left(\mathrm{6,2}\right)}}{P}_{\left(\mathrm{6,2}\right)}{\left(\pi {fGm}/3\right)}^{\tfrac{4}{3}}\right\}\\ & & \times \,{\rm{\Theta }}\left(6{f}_{\mathrm{LSO}}-f\right)\\ & & \times \,\exp \left\{i\left[2\pi {{ft}}_{c}-\displaystyle \frac{\pi }{4}+6{\rm{\Psi }}\left(f/6\right)\right]\right\},\end{array}\end{eqnarray} \tag{ 76 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{h}}^{\left(7\right)}\left(f\right) & = & {7}^{-\tfrac{1}{2}}{\left(\displaystyle \frac{5}{48}\right)}^{\tfrac{1}{2}}{\pi }^{-\tfrac{2}{3}}\displaystyle \frac{{\left({{GM}}_{c}\right)}^{5/6}}{D}{\left(2f/7\right)}^{-\tfrac{7}{6}}\\ & & \times \,\left\{{e}^{-i{\varphi }_{\left(\mathrm{7,5}/2\right)}}{P}_{\left(\mathrm{7,5}/2\right)}{\left(2\pi {fGm}/7\right)}^{\tfrac{5}{3}}\right\}\\ & & \times \,{\rm{\Theta }}\left(7{f}_{\mathrm{LSO}}-f\right)\\ & & \times \,\exp \left\{i\left[2\pi {{ft}}_{c}-\displaystyle \frac{\pi }{4}+7{\rm{\Psi }}\left(f/7\right)\right]\right\},\end{array}\end{eqnarray} \tag{ 77 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{S}_{1} & = & \displaystyle \frac{1}{2}\left(\displaystyle \frac{743}{336}+\displaystyle \frac{11}{4}\eta \right),\\ {S}_{3/2} & = & -2\pi ,\\ {S}_{2} & = & \displaystyle \frac{7266251}{8128512}+\displaystyle \frac{18913}{16128}\eta +\displaystyle \frac{1379}{1152}{\eta }^{2},\\ {S}_{5/2} & = & -\pi \displaystyle \frac{4757}{1344}-\displaystyle \frac{3}{16}\left(-63+44\pi \right)\eta ,\end{array}\end{eqnarray} \tag{ 78 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{e}^{-i{\varphi }_{\left(n,s\right)}}{P}_{\left(n,s\right)} & = & \left[{F}_{+}{C}_{+}^{\left(n,s\right)}+{F}_{\times }{C}_{\times }^{\left(n,s\right)}\right]\\ & & +\,i\left[{F}_{+}{D}_{+}^{\left(n,s\right)}+{F}_{\times }{D}_{\times }^{\left(n,s\right)}\right].\end{array}\end{eqnarray} \tag{ 79 }$$

Here, ${C}_{+,\times }^{\left(n,s\right)}$ and ${D}_{+,\times }^{\left(n,s\right)}$ denote the prefactors of $\cos (n{\rm{\Psi }})$ $\sin (n{\rm{\Psi }})$ in ${H}_{+,\times }^{(s)}$, respectively, which can be found in Arun et al. (2004, 2005), and we will not repeat them here. Note that f_LSO is the last orbital frequency where the waveforms are truncated. We employ the Roche radius of rigid spherical bodies in Equation (44) as a rough estimation of the last stable distance between two objects of a binary. The SMG corrections enter the waveform in harmonic one (71) and harmonic two (72). The terms relevant to (2πfGm)⁻¹ and ${(2\pi {fGm})}^{-\tfrac{1}{3}}$ in harmonic one (71) are directly added into the waveform. These two terms correspond to Equation (17). The term relevant to ${(\pi {fGm})}^{-\tfrac{2}{3}}$ in harmonic two (72), which corresponds to the terms having the same power of (πfGm) in Equations (21) and (23), is also directly added into the waveform. The above three terms are zero in the waveform of GR. The term relevant to (πfGm)⁰ in harmonic two (72) is the modified leading-order term, which can return to the case of GR when TF_b = 0 and Q = 1. The rest terms are all from higher order amplitude corrections in GR, which can also be found in Van Den Broeck & Sengupta (2006). Different from ground-based detectors, the antenna pattern functions of space-borne detectors depend on time. In the Fourier transform of a detector's response function obtained using stationary phase approximation, the time t in F₊, F_×, and F_b is replaced by function t(f), which is given by

$$\begin{eqnarray}&&t(f)={t}_{c}-\displaystyle \frac{5}{256{\left({{GM}}_{c}\right)}^{5/3}}{\left(2\pi f\right)}^{-8/3}\sum _{i=-2}^{7}{\tau }_{i}{\left(2\pi {fGm}\right)}^{i/3},\end{eqnarray} \tag{ 80 }$$

with the coefficients

$$\begin{eqnarray}\begin{array}{rcl}{\tau }_{-2} & = & -\displaystyle \frac{1}{48}{\epsilon }_{d}^{2},\\ {\tau }_{-1} & = & 0,\\ {\tau }_{0} & = & 1,\\ {\tau }_{1} & = & 0,\\ {\tau }_{2} & = & \displaystyle \frac{743}{252}+\displaystyle \frac{11}{3}\eta ,\\ {\tau }_{3} & = & -\displaystyle \frac{32}{5}\pi ,\\ {\tau }_{4} & = & \displaystyle \frac{3058673}{508032}+\displaystyle \frac{5429}{504}\eta +\displaystyle \frac{617}{72}{\eta }^{2},\\ {\tau }_{5} & = & -\left(\displaystyle \frac{7729}{252}-\displaystyle \frac{13}{3}\eta \right)\pi ,\\ {\tau }_{6} & = & -\displaystyle \frac{10052469856691}{23471078400}+\displaystyle \frac{128{\pi }^{2}}{3}\\ & & +\,\displaystyle \frac{6848\gamma }{105}+\left(\displaystyle \frac{3147553127}{3048192}-\displaystyle \frac{451{\pi }^{2}}{12}\right)\eta \\ & & -\,\displaystyle \frac{15211}{1728}{\eta }^{2}+\displaystyle \frac{25565}{1296}{\eta }^{3}+\displaystyle \frac{3424}{105}\mathrm{ln}\left[16{\left(2\pi {mf}\right)}^{2/3}\right],\\ {\tau }_{7} & = & \left(-\displaystyle \frac{15419335}{127008}-\displaystyle \frac{75703}{756}\eta +\displaystyle \frac{14809}{378}{\eta }^{2}\right)\pi .\end{array}\end{eqnarray} \tag{ 81 }$$

Here, γ = 0.5772 is the Euler–Mascheroni constant, τ₋₂ is concerned with the corrections induced by SMG, and τ_i(i ≥ 0) is concerned with the frequency evolution at 3.5 PN in phase (Buonanno et al. 2009). The phase Ψ(f) is given by

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Psi }}(f) & = & -{{\rm{\Psi }}}_{c}+\displaystyle \frac{3}{256{\left(2\pi {{fGM}}_{c}\right)}^{5/3}}\\ & & \times \,\sum _{i=-2}^{7}{{\rm{\Psi }}}_{i}{\left(2\pi {fGm}\right)}^{i/3}-{{\rm{\Psi }}}_{D}[t(f)],\end{array}\end{eqnarray} \tag{ 82 }$$

where the coefficients Ψ_i are given by

$$\begin{eqnarray*}\begin{array}{rcl}{{\rm{\Psi }}}_{-2} & = & -\displaystyle \frac{5}{336}{\epsilon }_{d}^{2},\\ {{\rm{\Psi }}}_{-1} & = & 0,\\ {{\rm{\Psi }}}_{0} & = & 1,\\ {{\rm{\Psi }}}_{1} & = & 0,\\ {{\rm{\Psi }}}_{2} & = & \displaystyle \frac{20}{9}\left[\displaystyle \frac{743}{336}+\displaystyle \frac{11}{4}\eta \right],\\ {{\rm{\Psi }}}_{3} & = & -16\pi ,\\ {{\rm{\Psi }}}_{4} & = & 10\left[\displaystyle \frac{3058673}{1016064}+\displaystyle \frac{5429}{1008}\eta +\displaystyle \frac{617}{144}{\eta }^{2}\right],\\ {{\rm{\Psi }}}_{5} & = & \pi \left[\displaystyle \frac{38645}{756}+\displaystyle \frac{38645}{756}\mathrm{ln}\left(\displaystyle \frac{f}{{f}_{\mathrm{LSO}}}\right)\right.\\ & & -\,\left.\displaystyle \frac{65}{9}\eta \left(1+\mathrm{ln}\left(\displaystyle \frac{f}{{f}_{\mathrm{LSO}}}\right)\right)\right],\\ {{\rm{\Psi }}}_{6} & = & \left(\displaystyle \frac{11583231236531}{4694215680}-\displaystyle \frac{640{\pi }^{2}}{3}-\displaystyle \frac{6848\gamma }{21}\right)\\ & & +\,\eta \left(-\displaystyle \frac{15737765635}{3048192}+\displaystyle \frac{2255{\pi }^{2}}{12}\right)\\ & & +\,\displaystyle \frac{76055}{1728}{\eta }^{2}-\displaystyle \frac{127825}{1296}{\eta }^{3}-\displaystyle \frac{6848}{21}\mathrm{ln}\left[4{\left(2\pi {fGm}\right)}^{1/3}\right],\\ {{\rm{\Psi }}}_{7} & = & \pi \left(\displaystyle \frac{77096675}{254016}+\displaystyle \frac{378515}{1512}\eta -\displaystyle \frac{74045}{756}{\eta }^{2}\right),\end{array}\end{eqnarray*}$$

with γ = 0.5772 the Euler–Mascheroni constant. Here, Ψ_i(i ≥ 0) is the coefficient in the 3.5 PN phase function of the Fourier domain waveform, and the coefficient Ψ₋₂ is concerned with the correction of dipole radiation in SMG. Note that Ψ_D denotes the Doppler modulation, which is the difference between the phase of the wavefront at the detector and at the barycenter. The expression of Ψ_D is given by (Cutler 1998; Hu et al. 2018)

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{D}=2\pi {fR}\sin \theta \cos \left[\displaystyle \frac{2\pi t(f)}{T}+{b}_{0}-\varphi \right],\end{eqnarray} \tag{ 83 }$$

where θ and φ are the ecliptic colatitude and longitude of the GW source, $R=1\,\mathrm{au}$, T is one year, and b₀ is the ecliptic longitude of the detector at t = 0. Finally, the other parameters, E, Q, S₋₁, and T, are defined by Equations (19), (24), (25), and (26).

The response functions of LISA-like detectors can be found in Cornish & Larson (2001), Cornish & Rubbo (2003), Rubbo et al. (2004), and Liang et al. (2019). In this appendix, we will give the expressions of antenna pattern function in a specific coordinate system.

The general form of antenna pattern function for LISA-like detectors is given in Equations (29) and (30). In the low-frequency range, the transfer functions ${ \mathcal T }(f,{\hat{{\boldsymbol{l}}}}_{1}\cdot \hat{{\boldsymbol{\Omega }}})$ and ${ \mathcal T }(f,{\hat{{\boldsymbol{l}}}}_{2}\cdot \hat{{\boldsymbol{\Omega }}})$ approach unity, and the antenna pattern functions return to the cases that are similar to the ground-based detectors. Here we will write explicitly the polarization tensor ${\epsilon }_{{ij}}^{A}$ and the vectors ${\hat{{\boldsymbol{l}}}}_{1}$, ${\hat{{\boldsymbol{l}}}}_{2}$, $\hat{{\boldsymbol{\Omega }}}$ in a specific coordinate system. This process is similar to the case of ground-based detectors (Maggiore 2008; Poisson & Will 2014) or the case of space-borne detectors in the low-frequency approximation (Cutler 1998; Hu et al. 2018).

We choose the coordinate system tied to the detector, which is denoted by $\hat{{\boldsymbol{x}}}\hat{{\boldsymbol{y}}}\hat{{\boldsymbol{z}}}$. The interferometer arms are put in this coordinate system, as shown in Figure 5. The unit vectors of two arms can be expressed as

$$\begin{eqnarray}&&{\hat{{\boldsymbol{l}}}}_{1}=\left(\begin{array}{c}\cos \tfrac{\pi }{12}\\ \sin \tfrac{\pi }{12}\\ 0\end{array}\right),\,\,\,\,{\hat{{\boldsymbol{l}}}}_{1}=\left(\begin{array}{c}\cos \tfrac{5\pi }{12}\\ \sin \tfrac{5\pi }{12}\\ 0\end{array}\right)\end{eqnarray} \tag{ 84 }$$

in this coordinate system.

Zoom In Zoom Out Reset image size

In a general metric theory, there are up to six possible polarization modes. Besides the h₊ and h_× modes in GR, there are the purely transverse h_b mode, the purely longitudinal h_l mode, and two mixed modes h_x and h_y. In the coordinates $({\hat{{\boldsymbol{x}}}}_{1},{\hat{{\boldsymbol{x}}}}_{2},{\hat{{\boldsymbol{x}}}}_{3})$ where the GW travels along the ${\hat{{\boldsymbol{x}}}}_{3}$ direction, these polarizations can be expressed as

$$\begin{eqnarray}{h}_{{ij}}=\left(\begin{array}{ccc}{h}_{b}+{h}_{+} & {h}_{\times } & {h}_{x}\\ {h}_{\times } & {h}_{b}-{h}_{+} & {h}_{y}\\ {h}_{x} & {h}_{y} & {h}_{l}\end{array}\right).\end{eqnarray} \tag{ 85 }$$

In SMG, there are four modes: h₊, h_×, h_b, and h_l (Liu et al. 2018b). We can use polarization tensors ${\epsilon }_{{ij}}^{A}$ to expand the metric perturbation as

$$\begin{eqnarray}&&{h}_{{ij}}(t)=\sum _{A}{\epsilon }_{{ij}}^{A}{h}_{A}(t),\end{eqnarray} \tag{ 86 }$$

where $A=+,\times ,\ b,l$ labels the polarization modes. By using the unit vector $\hat{{\boldsymbol{\Omega }}}$ (which points in the propagation direction of the GW) and the unit vectors $\hat{{\boldsymbol{u}}}$ and $\hat{{\boldsymbol{v}}}$ (which are orthogonal to $\hat{{\boldsymbol{\Omega }}}$ and orthogonal to each other), the polarization tensors can be rewritten as

$$\begin{eqnarray}\begin{array}{rcl}{\epsilon }_{{ij}}^{+} & = & {\hat{u}}_{i}{\hat{u}}_{j}-{\hat{v}}_{i}{\hat{v}}_{j},\,\,\,\,{\epsilon }_{{ij}}^{\times }={\hat{u}}_{i}{\hat{v}}_{j}+{\hat{v}}_{i}{\hat{u}}_{j},\\ {\epsilon }_{{ij}}^{b} & = & {\hat{u}}_{i}{\hat{u}}_{j}+{\hat{v}}_{i}{\hat{v}}_{j},\,\,\,\,{\epsilon }_{{ij}}^{l}={\hat{{\rm{\Omega }}}}_{i}{\hat{{\rm{\Omega }}}}_{j}.\end{array}\end{eqnarray} \tag{ 87 }$$

Thus, F_A can be derived straightforwardly once the unit vectors $\hat{{\boldsymbol{\Omega }}}$, $\hat{{\boldsymbol{u}}}$, $\hat{{\boldsymbol{v}}}$, and ${\hat{{\boldsymbol{l}}}}_{1}$, ${\hat{{\boldsymbol{l}}}}_{2}$ are expressed in the same coordinate system. We employ (θ', φ') to represent the direction of the GW source in the detector coordinate system, where φ' is the azimuth angle and θ' is the altitude angle (we use the definition of θ', which is the angle between the direction of the GW source and the direction of $\hat{{\boldsymbol{z}}}$ in this work). Here, ψ' denotes the polarization angle in the detector coordinate system. Therefore, the vectors $\hat{{\boldsymbol{\Omega }}}$, $\hat{{\boldsymbol{u}}}$, $\hat{{\boldsymbol{v}}}$ in the detector coordinate system can be given by

$$\begin{eqnarray}\begin{array}{rcl}\hat{{\boldsymbol{u}}} & = & \left(\begin{array}{c}\cos \theta ^{\prime} \cos \varphi ^{\prime} \cos \psi ^{\prime} -\sin \varphi ^{\prime} \sin \psi ^{\prime} \\ \cos \theta ^{\prime} \sin \varphi ^{\prime} \cos \psi ^{\prime} +\cos \varphi ^{\prime} \sin \psi ^{\prime} \\ -\sin \theta ^{\prime} \cos \psi ^{\prime} \end{array}\right),\\ \hat{{\boldsymbol{v}}} & = & \left(\begin{array}{c}\cos \theta ^{\prime} \cos \varphi ^{\prime} \sin \psi ^{\prime} +\sin \varphi ^{\prime} \cos \psi ^{\prime} \\ \cos \theta ^{\prime} \sin \varphi ^{\prime} \sin \psi ^{\prime} -\cos \varphi ^{\prime} \cos \psi ^{\prime} \\ -\sin \theta ^{\prime} \sin \psi ^{\prime} \end{array}\right),\\ \hat{{\boldsymbol{\Omega }}} & = & \left(\begin{array}{c}-\sin \theta ^{\prime} \cos \varphi ^{\prime} \\ -\sin \theta ^{\prime} \sin \varphi ^{\prime} \\ -\cos \theta ^{\prime} \end{array}\right).\end{array}\end{eqnarray} \tag{ 88 }$$

Different from the ground-based detectors, the timescale of a GW signal detected by space-borne detectors is comparable to the timescale of the detector motion. The motion of space-borne detectors cannot be neglected, so (θ', φ', ψ') are considered to be time-dependent. We need to find the relationships between the direction of the GW source, as well as the polarization angle in the detector coordinate system, and those in the heliocentric coordinate system, which are denoted by (θ, φ, ψ). The relationships depend on the motion of the detectors. We will discuss LISA first.

We employ $\hat{{\boldsymbol{i}}}\hat{{\boldsymbol{j}}}\hat{{\boldsymbol{k}}}$ to denote the heliocentric coordinates tied to the ecliptic and (θ, φ, ψ) to represent the direction of the GW source and the polarization angle in the heliocentric coordinate system, respectively. Recalling the orbital configuration of LISA, which was introduced in Section 3, the unit vectors $\hat{{\boldsymbol{x}}}$, $\hat{{\boldsymbol{y}}}$, and $\hat{{\boldsymbol{z}}}$ of the detector coordinate system can be written in terms of the heliocentric coordinate system as

$$\begin{eqnarray}\begin{array}{rcl}\hat{{\boldsymbol{x}}} & = & \left[\displaystyle \frac{1}{2}\cos a(t)\cos b(t)+\sin a(t)\sin b(t)\right]\hat{{\boldsymbol{i}}}\\ & & +\,\left[\displaystyle \frac{1}{2}\cos a(t)\sin b(t)-\sin a(t)\cos b(t)\right]\\ & & \times \,\hat{{\boldsymbol{j}}}+\left[\displaystyle \frac{\sqrt{3}}{2}\cos a(t)\right]\hat{{\boldsymbol{k}}}\\ \hat{{\boldsymbol{y}}} & = & \left[\displaystyle \frac{1}{2}\sin a(t)\cos b(t)-\cos a(t)\sin b(t)\right]\hat{{\boldsymbol{i}}}\\ & & +\,\left[\displaystyle \frac{1}{2}\sin a(t)\sin b(t)+\cos a(t)\cos b(t)\right]\\ & & \times \,\hat{{\boldsymbol{j}}}+\left[\displaystyle \frac{\sqrt{3}}{2}\sin a(t)\right]\hat{{\boldsymbol{k}}}\\ \hat{{\boldsymbol{z}}} & = & \left[-\displaystyle \frac{\sqrt{3}}{2}\cos b(t)\right]\hat{{\boldsymbol{i}}}+\left[-\displaystyle \frac{\sqrt{3}}{2}\sin b(t)\right]\,\hat{{\boldsymbol{j}}}+\left(\displaystyle \frac{1}{2}\right)\hat{{\boldsymbol{k}}},\end{array}\end{eqnarray} \tag{ 89 }$$

where $a(t)={a}_{0}+\tfrac{2\pi t}{{T}_{\mathrm{LISA}}}$ is the phase of rotation around the detector's center, and $b(t)={b}_{0}+\tfrac{2\pi t}{{T}_{\mathrm{LISA}}}$ is the phase of revolution around the Sun. For the motion of LISA, the periods of rotation around the detector's center and the revolution around the Sun are both one year. The initial phases a₀ and b₀ are constant. We can take a₀ = 0 and b₀ = 0 without loss of generality. The direction of the GW source $\hat{{\boldsymbol{r}}}$ can be given in $\hat{{\boldsymbol{i}}}\hat{{\boldsymbol{j}}}\hat{{\boldsymbol{k}}}$ coordinates as

$$\begin{eqnarray}&&\hat{{\boldsymbol{r}}}=(\sin \theta \cos \varphi )\hat{{\boldsymbol{i}}}+(\sin \theta \sin \varphi )\hat{{\boldsymbol{j}}}+(\cos \theta )\hat{{\boldsymbol{k}}}.\end{eqnarray} \tag{ 90 }$$

Using the geometry relationships of those vectors and angles, θ' and φ' can be obtained:

$$\begin{eqnarray}&&\cos \theta ^{\prime} =\hat{{\boldsymbol{r}}}\cdot \hat{{\boldsymbol{z}}},\,\,\,\tan \phi ^{\prime} =\displaystyle \frac{\hat{{\boldsymbol{r}}}\cdot \hat{{\boldsymbol{y}}}}{\hat{{\boldsymbol{r}}}\cdot \hat{{\boldsymbol{x}}}}.\end{eqnarray} \tag{ 91 }$$

As for polarization angle ψ', we follow the definition in Cutler (1998). An ellipse can be obtained by projecting the binary's circular orbit on the plane of the sky (i.e., the plane orthogonal to the GW propagation direction). The major axis of this ellipse is defined as the vector $\hat{{\boldsymbol{u}}}$ mentioned above. The polarization angle is defined as the angle between the vector $\hat{{\boldsymbol{u}}}$ and the vector pointing in the direction of increasing θ'. A useful figure can be referred to in the literature (Poisson & Will 2014, Figure 11.5). According to this definition, the polarization is given by

$$\begin{eqnarray}&&\tan \psi ^{\prime} =\displaystyle \frac{\left[\hat{{\boldsymbol{L}}}-(\hat{{\boldsymbol{L}}}\cdot \hat{{\boldsymbol{r}}})\hat{{\boldsymbol{r}}}\right]\cdot \hat{{\boldsymbol{z}}}}{(\hat{{\boldsymbol{r}}}\times \hat{{\boldsymbol{L}}})\cdot \hat{{\boldsymbol{z}}}},\end{eqnarray} \tag{ 92 }$$

where $\hat{{\boldsymbol{L}}}$ denotes the unit vector parallel to the orbital angular momentum vector of the binary. The vector $\hat{{\boldsymbol{L}}}$ in the detector coordinates $\hat{{\boldsymbol{x}}}\hat{{\boldsymbol{y}}}\hat{{\boldsymbol{z}}}$ is time-dependent, so we prefer to express $\hat{{\boldsymbol{L}}}$ in the heliocentric coordinates $\hat{{\boldsymbol{i}}}\hat{{\boldsymbol{j}}}\hat{{\boldsymbol{k}}}$. The vector $\hat{{\boldsymbol{L}}}$ in the coordinates $\hat{{\boldsymbol{i}}}\hat{{\boldsymbol{j}}}\hat{{\boldsymbol{k}}}$ can be given by

$$\begin{eqnarray}\begin{array}{rcl}\hat{{\boldsymbol{L}}} & = & \left[\cos \iota \sin \theta \cos \varphi \right.\\ & & +\,\left.\sin \iota \left(\cos \psi \sin \varphi +\cos \theta \cos \varphi \sin \psi \right)\right]\hat{{\boldsymbol{i}}}\\ & & +\,\left[-\sin \iota \cos \psi \cos \varphi \right.\\ & & +\,\left.\sin \varphi \left(\cos \iota \sin \theta +\cos \theta \sin \iota \sin \psi \right)\right]\hat{{\boldsymbol{j}}}\\ & & +\,\left[\cos \iota \cos \theta -\sin \theta \sin \iota \sin \psi \right]\hat{{\boldsymbol{k}}}.\end{array}\end{eqnarray} \tag{ 93 }$$

In the above equation, inclination angle ι is the angle between $\hat{{\boldsymbol{L}}}$ and $\hat{{\boldsymbol{r}}}$, and the polarization angle ψ in heliocentric coordinates $\hat{{\boldsymbol{i}}}\hat{{\boldsymbol{j}}}\hat{{\boldsymbol{k}}}$ has the definition similar to ψ' in the coordinates $\hat{{\boldsymbol{x}}}\hat{{\boldsymbol{y}}}\hat{{\boldsymbol{z}}}$, which is the angle between the major axis of the projection ellipse and the vector pointing to the direction of increasing θ.

The parameters (θ', φ', ψ') in the detector coordinates can be eventually expressed in terms of the parameters (θ, φ, ψ, ι) in the heliocentric coordinates, which are considered to be time-independent, and the variables a(t), b(t) describing the motion of detectors in the heliocentric coordinates, which have simple relationships with time. Substituting Equation (88) into Equation (87), we can obtain the polarization tensors in the detector coordinates. The transfer functions can also be obtained by putting the expressions of ${\hat{{\boldsymbol{l}}}}_{1}$, ${\hat{{\boldsymbol{l}}}}_{2}$, $\hat{{\boldsymbol{\Omega }}}$ (Equations (88), (84)) into Equation (30). Substituting these results and Equation (84) into Equation (29), the antenna pattern functions can finally be assembled. The final results are straightforward but cumbersome. To avoid redundancy, the results are not presented here.

A similar calculation is also applicable for TianQin. Compared with LISA, TianQin will run in a geocentric orbit, and the orientation of TianQin is fixed to the reference source J0806.3+1527 instead of varying with time. A brief introduction to TianQin's orbit was given in Section 3. The base vectors of TianQin's detector coordinates ${(\hat{{\boldsymbol{x}}},\hat{{\boldsymbol{y}}},\hat{{\boldsymbol{z}}})}_{\mathrm{TianQin}}$ can be given by

$$\begin{eqnarray}\begin{array}{rcl}\hat{{\boldsymbol{x}}} & = & \left[\cos a(t)\cos {\theta }_{0}\cos {\varphi }_{0}-\sin a(t)\sin {\varphi }_{0}\right]\hat{{\boldsymbol{i}}}\\ & & +\,\left[\cos a(t)\cos {\theta }_{0}\sin {\varphi }_{0}+\sin a(t)\cos {\varphi }_{0}\right]\hat{{\boldsymbol{j}}}\\ & & +\,\left[-\cos a(t)\sin {\theta }_{0}\right]\hat{{\boldsymbol{k}}},\\ \hat{{\boldsymbol{y}}} & = & \left[-\sin a(t)\cos {\theta }_{0}\cos {\varphi }_{0}-\cos a(t)\sin {\varphi }_{0}\right]\hat{{\boldsymbol{i}}}\\ & & +\,\left[-\sin a(t)\cos {\theta }_{0}\sin {\varphi }_{0}+\cos a(t)\cos {\varphi }_{0}\right]\hat{{\boldsymbol{j}}}\\ & & +\,\left[\sin a(t)\sin {\theta }_{0}\right]\hat{{\boldsymbol{k}}},\\ \hat{{\boldsymbol{z}}} & = & \left(\sin {\theta }_{0}\cos {\varphi }_{0}\right)\hat{{\boldsymbol{i}}}+\left(\sin {\theta }_{0}\sin {\varphi }_{0}\right)\hat{{\boldsymbol{j}}}+\left(\cos {\theta }_{0}\right)\hat{{\boldsymbol{k}}},\end{array}\end{eqnarray} \tag{ 94 }$$

in the heliocentric coordinates, where (θ₀, φ₀) denote the direction of the reference source, and $a(t)={a}_{0}+\tfrac{2\pi t}{{T}_{\mathrm{TianQin}}}$ represents the rotation phase of the detector. Here, T_TianQin is the period of TianQin's rotation, which is about 3.65 days. Following the same process as discussed above, the antenna pattern functions of TianQin can be derived.