# Testing the polarization of gravitational wave background with LISA-TianQin
network

Yu Hu  Pan-Pan Wang  Yu-Jie Tan  Cheng-Gang Shao 11footnotetext: Corresponding
author.

###### Abstract

While general relativity predicts only two tensor modes for gravitational wave
polarization, general metric theories of gravity allows up to four additional
modes, including two vector and two scalar modes. Observing the polarization
modes of gravitational waves could provide a direct test of the modified
gravity. The stochastic gravitational wave background (SGWB), which may be
detected by space-based laser-interferometric detectors at design sensitivity,
will provide an opportunity to directly measure alternative polarization. In
this paper, we investigate the performance of the LISA-TianQin network for
detecting alternative polarizations of stochastic backgrounds, and propose a
method to separate different polarization modes. First, we generalize the
small antenna approximation to compute the overlap reduction functions for
SGWB with arbitrary polarization, which is suitable for any time-delay
interferometry combination. Then we analyze the detection capability of LISA-
TianQin for SGWB with different polarizations. Based on the LISA-TianQin
orbital characteristics, we propose a method to distinguish different
polarization modes from their mixed data. Compared with ground-based
detectors, the LISA-TianQin network is more capable of resolving polarizations
of SGWB. In particular, the LISA-TianQin network has the potential to resolve
two scalar modes that ground-based detectors cannot.

##  1 Introduction

Since the first gravitational wave (GW) signal GW150914 was detected by
Advanced LIGO in 2015 [1], Advanced LIGO and later Advanced Virgo have
detected more and more GWs sourced from compact binary coalescences[2, 3, 4,
5, 6]. At the same time, this also implies that there should be many weak GWs
that cannot be recognized by detectors. The combined weak signal from the
population of compact binary constitutes the stochastic gravitational wave
background (SGWB). In addition to the astrophysical sources, there are many
ways to generate SGWB in the early Universe, such as cosmological phase
transitions [7, 8], inflation [9, 10, 11, 12], cosmic strings [13, 14], etc.
The detection of the GWB can provide us with information on the astronomical
distribution [15, 16] and cosmology [17, 18, 19, 20, 21], and also provide an
opportunity to test the theory of gravity [22, 23].

Pulsar Timing Arrays (PTA) has accumulated more than a decade of data, and is
expected to detect SGWB in the nHz band in the near future [24, 25, 26, 27].
In the Hz hand, the ground-based laser interferometers give an upper limits of
fractional energy density of SGWB $\Omega_{\text{GW}}\leq 5.8\times 10^{-9}$
[28]. The future space-based interferometers such as LISA [29], TianQin [30],
Taiji [31] and DECIGO [32], will be able to detect SGWB with high sensitivity
in the mHz band [33, 34, 35]. Once we have detected the SGWB, we next need to
know what information can be obtained from it. An important role of SGWB is to
test the theory of gravity. General relativity (GR) predicts that there only
exists two tensor polarizations for gravitational wave, the plus and cross
modes. However, there are additional four polarizations allowed by the generic
metric theories of gravity, including two vector modes and two scalar
modes[36, 37, 38, 39, 40, 41, 42, 43]. The observation of vector or scalar
modes would cast doubt on general relatively, and the absence could also be
used to constrain modified gravity [44]. If parity is considered, there is
circular polarization in the SGWB, describing the asymmetry of the two tensor
polarizations [45, 46, 47, 48].

Research on the detection of polarization modes of the SGWB has become a hot
topic, as more and more modified theories of gravity are proposed. For PTA,
the detectability of non-GR polarizations of the SGWB was first investigated
in [49]. After that, more detailed extension can be found in [50, 51]. No
evidence of non-GR polarizations of SGWB has been found in more than a decade
of PTA data. The constrains of the amplitude or energy density for non-GR
modes in SGWB at frequency of 1/yr can be found in [52, 53, 54]. On the other
hand, the ground-based gravitational-wave detection network is gradually
forming, as more and more detectors join in. It has the potential to detect
the SGWB produced by compact binary mergers. Data from Advanced LIGO’s and
Advanced Virgo’s three observing run can constrain the fractional energy
density of different polarization modes at Hz band [55, 56, 57]. However, the
ability of the ground-based laser interferometer network to study the
polarizations of SGWB is limited. The current sensitivity of ground-based
detectors is limited enough to detect the SGWB. What more, the two kind of
scalar modes are completely different, one is the transverse and the other is
longitudinal. Unfortunately, the responses of the ground-based laser
interferometer to scalar-breathing and scalar-longitudinal modes are
completely degenerate, which means that the two modes cannot be distinguished,
no matter how sensitive the detectors are [58]. However, the two scalar
polarizations are of a completely different nature, and they cannot be
converted into one another via a rotational coordinate transformation along
the direction of propagation (whereas two tensor modes or two vector modes
can). The physical images behind the appearance of scalar breathing and
longitudinal modes are distinct. For example, in the scalar-tensor theory of
gravity, there exists only one extra breathing mode if the scalar field is
massless, or a mixture of breathing mode and longitudinal mode if the scalar
field is massive [59]. It is therefore necessary to distinguish between the
two scalar modes in the detector.

For the future space-based gravitational-wave detectors, the situation is
different. The degeneracy of response between the breathing and longitudinal
modes is broken at relatively high frequencies [60, 61], which implies that
the space-based detectors have potential to distinguish breathing and
longitudinal modes of gravitational waves. Besides, space-based detectors can
detect signals from a much larger range of frequencies, compared to ground-
based detectors. The abundance of sources in the sensitive band and the
sufficiently high sensitivity of the detectors ensure that the polarization
analysis can be performed. What more, the relative orientation of two space-
based detectors may change as they are in different positions in their
respective orbits. For example, for the LISA-TianQin network, the normal
vector of the constellation plane of LISA is vary with time and that of
TianQin points to a fixed direction [30]. This means that the overlap
reduction function (ORF), the transfer function between the spectrum of SGWB
and the power spectrum of cross correlation signal, will vary accordingly. In
this way, it is easier to distinguish between different polarization modes,
which can be inferred from the viewpoint that the detectors can be regarded as
different detectors at different positions.

For the space-based detectors, the laser phase noise is usually orders of
magnitude higher than other noises due to the mismatch of arm lengths, and
also much larger than the gravitational wave signal. Fortunately, the time
delay interferometry (TDI) can be used to suppress the laser phase noise [62,
63, 64, 65]. It is worth mentioning that three noise quadrature channels can
be constructed, called A, E and T. The T channel is relatively insensitive to
the gravitational-wave compared with other channels. So the readout of T
channel can be used to model the noise in order to subtract the noise from the
A and E channels to get the SGWB. This is called the null channel method [66,
67, 68]. However, the persuasiveness of the null channel is limit, it is
required that the cross-correlation signals are detected in the data of two
detectors to claim detection of SGWB. The cross-correlation analysis based on
TDI combination is a standard method for space-based detectors to detect the
SGWB. However, there is a lack of study on the performance of the space-based
detectors for SGWB with different polarizations, as well as the treatment of
the time varying effect. Among the promising space-based detectors, the
difference in the orbital configuration of LISA and TianQin is most
significant, so the time-varying effect of LISA-TianQin network is the most
pronounced. In this paper, we investigate the detection capability of the
LISA-TianQin network for SGWB with different polarization, and propose a
method inspired by the time-varying effect to distinguish different
polarization modes from their mixed data.

The outline of the paper is as follows. In section 2, we review the SGWB in
general metric theories of gravity and introduce the LISA-TianQin network. In
section 3, we review the correlation analysis for detecting SGWB and calculate
the overlap reduction function for LISA-TianQin network. In section 4, we
study the sensitivity and detectability for the SGWB of alternative
polarizations. Then, in section 5, we consider the method to separation the
polarizations for LISA-TianQin network. Finally, a discussion is presented in
in section 6.

##  2 SGWB statistic and LISA-TianQin network

###  2.1 Stochastic background of non-GR polarizations

The metric perturbations corresponding to SGWB can be expressed as a
superposition of plane waves of different frequencies from different
directions [68]:

$h_{ab}(t,\vec{x})=\int_{-\infty}^{\infty}df\int
d^{2}\Omega_{\hat{n}}h_{ab}(f,\hat{n})e^{i2\pi f(t+\hat{n}\cdot\vec{x}/c)}.$
(2.1)

The fourier coefficients $h_{ab}(f,\hat{n})$ are random variables, whose
statistic is significant. In generic metric theory, the coefficients can be
expanded in terms of the six spin-2 polarization tensors:

$h_{ab}(f,\hat{n})=\sum_{A}h_{A}(f,\hat{n})e^{A}_{ab}(\hat{n}),$ (2.2)

where $A=\\{+,\times,X,Y,B,L\\}$ represents different polarization modes,
where ${+,\times}$ represent tensor modes predicted by general relativity, and
${X,Y}$ and ${B,L}$ represent vector and scalar modes allowed by the general
metric theory of gravity. Explicitly, the six spin-2 polarization tensors are

$\displaystyle e^{+}_{ab}(\hat{n})=$
$\displaystyle\hat{\theta}_{a}\hat{\theta}_{b}-\hat{\phi}_{a}\hat{\phi}_{b},\quad
e^{\times}_{ab}(\hat{n})=\hat{\theta}_{a}\hat{\phi}_{b}+\hat{\phi}_{a}\hat{\theta}_{b},$
(2.3) $\displaystyle e^{X}_{ab}(\hat{n})=$
$\displaystyle\hat{\theta}_{a}\hat{n}_{b}+\hat{n}_{a}\hat{\theta}_{b},\quad
e^{Y}_{ab}(\hat{n})=\hat{\phi}_{a}\hat{n}_{b}+\hat{n}_{a}\hat{\phi}_{b},$
$\displaystyle e^{B}_{ab}(\hat{n})=$
$\displaystyle\hat{\theta}_{a}\hat{\theta}_{b}+\hat{\phi}_{a}\hat{\phi}_{b},\quad
e^{L}_{ab}(\hat{n})=\sqrt{2}\hat{n}_{a}\hat{n}_{b},$

and $\hat{\theta}$, $\hat{\phi}$ are the standard angular unit vectors tangent
to the sphere:

$\displaystyle\hat{n}$
$\displaystyle=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta),$ (2.4)
$\displaystyle\hat{\theta}$
$\displaystyle=(\cos\theta\cos\phi,\cos\theta\sin\phi,-\sin\theta),$
$\displaystyle\hat{\phi}$ $\displaystyle=(-\sin\phi,\cos\phi,0).$

The statistical properties of the SGWB are described by the probability
distribution of the metric perturbations. In this work, we assume that the
SGWB is Gaussian, stationary and isotropic. And without loss of generality we
can assume that the background has zero mean $\left\langle
h_{A}(f,\hat{n})\right\rangle=0$. So all the information is encoded in the
quadratic expectation:

$\left\langle
h_{A}(f,\hat{n})h_{A^{\prime}}^{*}\left(f^{\prime},\hat{n}^{\prime}\right)\right\rangle=\frac{1}{8\pi}S^{A}_{h}(f)\delta\left(f-f^{\prime}\right)\delta_{AA^{\prime}}\delta^{2}\left(\hat{n},\hat{n}^{\prime}\right),$
(2.5)

where $S^{A}_{h}(f)$ can be regarded as the component corresponding to the $A$
polarization of a one-sided gravitational-wave strain power spectral density
function. We further assume that both the tensor and vector modes are
unpolarized, which implies that

$\displaystyle S^{+}_{h}$ $\displaystyle=S^{\times}_{h}=S^{T}_{h}/2,$ (2.6)
$\displaystyle S^{X}_{h}$ $\displaystyle=S^{Y}_{h}=S^{V}_{h}/2.$

However, the two scalar modes should be considered as two independent
polarization modes, since one is the longitudinal and the other is transverse.
It is worth noting that some models predict the presence of parity violating
SGWB sources [69, 70], which can be recognized with the interferometers [71,
72, 73].

The function $S^{A}_{h}(f)$, which characterizes the spectral shape of the
SGWB within each polarization sector, can be detected directly without
assuming a model. However, the amplitude of SGWB for each polarization is
characterized by the fractional energy density [74],

$\Omega^{A}_{\mbox{gw}}(f)=\frac{1}{\rho_{c}}\frac{d\rho^{A}_{\mbox{gw}}}{d\ln
f},$ (2.7)

defined as the energy density per logarithmic frequency bin, normalized by the
critical energy density of the closed Universe $\rho_{c}\equiv
3c^{2}H_{0}^{2}/8\pi G$. Here $G$ is the gravitational constant, and
$H_{0}=67.4\text{km s}^{-1}\text{Mpc}^{-1}$ is the Hubble constant [75]. In
general relatively, the relation between $S^{A}_{h}(f)$ and
$\Omega^{A}_{\mbox{gw}}(f)$ is [74]

$\Omega^{A}_{\mbox{gw}}(f)=\frac{2\pi^{2}}{3H_{0}^{2}}f^{3}S^{A}_{h}(f).$
(2.8)

In alternative theories of gravity, eq. (2.8) may not hold unless the stress-
energy of gravitational waves also obeys Isaacson’s formula [76]:

$\rho_{\mbox{gw}}=\frac{c^{2}}{32\pi
G}\left\langle\dot{h}_{ab}(t,\vec{x})\dot{h}^{ab}(t,\vec{x})\right\rangle.$
(2.9)

In this case, $\Omega^{A}_{\mbox{gw}}(f)$ can be understood as a function of
the observable $S^{A}_{h}(f)$ rather than the fractional energy density.

Many theoretical models of SGWB predict that the shape of
$\Omega^{A}_{\mbox{gw}}(f)$ can be modeled as power laws [68], such that

$\Omega^{A}_{\mbox{gw}}(f)=\Omega^{\alpha_{A}}_{0}\left(\frac{f}{f_{0}}\right)^{\alpha_{A}}.$
(2.10)

Here $\Omega^{\alpha_{A}}_{0}$ is the amplitude of polarization $A$ at a
reference frequency $f_{0}$ and $\alpha_{A}$ is the corresponding spectral
index. For instance, the tensor polarization background from compact binary
coalescences is modeled by power law with index $\alpha_{T}=2/3$ [56] and for
the inflationary cosmic background is $\alpha_{T}=0$ [77].

###  2.2 LISA-TianQin network

The LISA and TianQin are proposed space GW missions targeting to detect the GW
in the frequency band 0.1 mHz – 1 Hz. The difference is that LISA is
heliocentric orbit and TianQin is geocentric. In addition, the relative angle
between their detecter plane will vary with time. TianQin has a geocentric
orbit and consists of three satellite formations in a nearly equilateral
triangle. Accurate to the first order, the coordinates of the three satellites
of TianQin are [78]

$\vec{r}^{TQ}_{n}=\vec{r}^{TQ}_{0}+\vec{R}^{TQ}_{n},$ (2.11)

where $\vec{r}^{TQ}_{0}=(x^{TQ}_{0},y^{TQ}_{0},z^{TQ}_{0})$ is the geocentric
coordinate,

$\displaystyle x^{TQ}_{0}$
$\displaystyle=R\cos\alpha_{TQ}+\frac{1}{2}e^{TQ}R(\cos 2\alpha_{TQ}-3),$
(2.12) $\displaystyle y^{TQ}_{0}$
$\displaystyle=R\sin\alpha_{TQ}+\frac{1}{2}e^{TQ}R\sin 2\alpha_{TQ},$
$\displaystyle z^{TQ}_{0}$ $\displaystyle=0,$

and $\vec{R}^{TQ}_{n}=(X^{TQ}_{n},Y^{TQ}_{n},Z^{TQ}_{n})$ are the coordinates
of the satellites in the geocentric coordinate system,

$\displaystyle X^{TQ}_{n}$
$\displaystyle=R_{1}\left(\cos\phi_{s}\sin\theta_{s}\sin\alpha_{n}+\cos\alpha_{n}\sin\phi_{s}\right),$
(2.13) $\displaystyle Y^{TQ}_{n}$
$\displaystyle=R_{1}\left(\sin\phi_{s}\sin\theta_{s}\sin\alpha_{n}-\cos\alpha_{n}\cos\phi_{s}\right),$
$\displaystyle Z^{TQ}_{n}$ $\displaystyle=-R_{1}\sin\alpha_{n}\cos\alpha_{n}.$

Here $R=1\text{AU}$, $e^{TQ}=0.0167$, $\alpha_{TQ}=2\pi f_{m}t-\alpha_{0}$,
$f_{m}=1/\text{yr}$, $\alpha_{0}=102.9^{\circ}$, $\alpha_{n}=2\pi
f_{sc}t+\kappa_{n}$, $\kappa_{n}=2/3(n-1)\pi$, $R_{1}=1\times 10^{8}m$,
$\theta_{s}=-4.7^{\circ}$, $\phi_{s}=120.5^{\circ}$,
$f_{sc}=1/(3.64\text{days})$ and $n=1,2,3$ respectively denotes one of the
three satellites. The detector plane orientation is fixed as
$(\cos\theta_{s}\cos\phi_{s},\cos\theta_{s}\sin\phi_{s},\sin\theta_{s})$ and
the arm length is $L^{TQ}=\sqrt{3}\times 10^{8}\text{m}$. The displacement
measurement noise $\sqrt{S^{TQ}_{x}}=1\times
10^{-12}\text{m}/\text{Hz}^{-1/2}$ and the residual acceleration noise
$\sqrt{S^{TQ}_{a}}=1\times 10^{-15}\text{m s}^{-2}/\text{Hz}^{-1/2}$. LISA has
a heliocentric orbit at $20^{\circ}$ behind the Earth. The satellite formation
consists of three satellites to form an approximate equilateral triangle [79],
and the coordinates of the three LISA satellites are
$r^{LS}_{n}=(x^{LS}_{n},y^{LS}_{n},z^{LS}_{n})$, namely,

$\displaystyle x^{LS}_{n}$
$\displaystyle=R\cos\alpha_{LS}+\frac{1}{2}e^{LS}R\left(\cos\left(2\alpha_{LS}-\kappa_{n}\right)-3\cos\kappa_{n}\right),$
(2.14) $\displaystyle y^{LS}_{n}$
$\displaystyle=R\sin\alpha_{LS}+\frac{1}{2}e^{LS}R\left(\sin\left(2\alpha_{LS}-\kappa_{n}\right)-3\sin\kappa_{n}\right),$
$\displaystyle z^{LS}_{n}$
$\displaystyle=-\sqrt{3}e^{LS}R\cos(\alpha_{LS}-\kappa_{n}),$

where $\alpha_{LS}=\alpha_{TQ}-20^{\circ}$ and $e^{LS}=0.0048$. The detector
plane is inclined to the orbit plane by $60^{\circ}$ and the arm length is
$L^{LS}=2.5\times 10^{9}\text{m}$. The displacement measurement noise
$\sqrt{S^{LS}_{x}}=1.5\times 10^{-11}\text{m}/\text{Hz}^{-1/2}$ and the
residual acceleration noise $\sqrt{S^{LS}_{a}}=3\times 10^{-15}\text{m
s}^{-2}/\text{Hz}^{-1/2}$.

###  2.3 Noise and response for time delay interferometry

Figure 1: The configuration of the space-based laser interferometry consists
of laser sources and links.

There are six laser links between the three satellites. Laser noise can be
effectively reduced by constructing time delay interferometry combinations.
Any TDI combination can be expressed in terms of a polynomial of the delay
operator acting on the six received signals,

$s_{TDI}=\sum_{i}P_{i}s_{i}.$ (2.15)

Here $i=1,2,3,1^{\prime},2^{\prime},3^{\prime}$ represents a link
respectively: $2\rightarrow 1$, $3\rightarrow 2$, $1\rightarrow 3$,
$3\rightarrow 1$, $1\rightarrow 2$, $2\rightarrow 3$. The detector
configuration is illustrated in figure 1. The time delay operator is defined
as $D_{i}s_{j}(t)=s_{j}(t-L_{i}/c)$, and converted to the frequency domain to
$\tilde{D}_{i}=e^{-i2\pi fL_{i}/c}$. For example, the coefficients of the
first-generation TDI Michelson combination $X$ are given by

$\displaystyle P_{1}$
$\displaystyle=D_{2^{\prime}2}-1,P_{2}=0,P_{3}=D_{2^{\prime}}-D_{33^{\prime}2^{\prime}}$
(2.16) $\displaystyle P_{1^{\prime}}$
$\displaystyle=1-D_{33^{\prime}},P_{2^{\prime}}=D_{2^{\prime}23}-D_{3},P_{3^{\prime}}=0,$

which represents a laser interferometry link:

$[1\rightarrow 2\rightarrow 1\rightarrow 3\rightarrow 1]-[1\rightarrow
3\rightarrow 1\rightarrow 2\rightarrow 1].$

The equivalent expressions for Michelson channel $Y$ and $Z$ can be obtained
by permuting the label $\\{1,2,3\\}$. Three noise-orthogonal channels can be
constrained by the linear combinations of $X,Y,Z$,

$A=\frac{Z-X}{\sqrt{2}},E=\frac{X-2Y+Z}{\sqrt{6}},T=\frac{X+Y+Z}{\sqrt{3}}.$
(2.17)

After removing the phase noise, there are residual acceleration noise and
displacement measurement noise in the TDI combination. The power spectral
densities (PSDs) of the remaining noises is

$\displaystyle P_{n,TDI}=$
$\displaystyle\frac{1}{L^{2}}\left[C_{1}S_{x}\left(1+(2\text{mHz}/f)^{4}\right)\right.$
(2.18)
$\displaystyle\left.+(2C_{1}+C_{2}\cos\beta)\frac{S_{a}\left(1+(0.4\text{mHz}/f)^{2}\right)\left(1+(f/8\text{mHz})^{2}\right)}{(2\pi
f)^{4}}\right],$

where $\beta=2\pi fL/c$ and the coefficients

$\displaystyle C_{1}$ $\displaystyle=\sum_{i}|\tilde{P}_{i}|^{2},$ (2.19)
$\displaystyle C_{2}$
$\displaystyle=\Re(\tilde{P}_{1}\tilde{P}^{*}_{2^{\prime}}+\tilde{P}_{2}\tilde{P}^{*}_{3^{\prime}}+\tilde{P}_{3}\tilde{P}^{*}_{1^{\prime}}).$

For example, the PSD of $X$ channel is

$\displaystyle P_{n,X}=$
$\displaystyle\frac{16\sin^{2}\beta}{L^{2}}\left[S_{x}\left(1+(2\text{mHz}/f)^{4}\right)\right.$
(2.20)
$\displaystyle\left.+\frac{S_{a}\left(1+(0.4\text{mHz}/f)^{2}\right)\left(1+(f/8\text{mHz})^{2}\right)}{(2\pi
f)^{4}}(3+\cos 2\beta)\right].$

The PSDs of X and A channel for LISA and TianQin are shown in figure 2.

Figure 2: The PSD $P_{n}(f)$ for different channels. The solid and dashed
lines denote the X and A channels respectively. And we adopt the red/blue to
label the LISA/TianQin.

Since the gravitational waves are weak, it is accurate enough to calculate the
detecter response to the linear order. In frequency domain, the gravitational
waves signal can be expressed as [68]

$\tilde{h}(f)=\int
d^{2}\Omega_{\hat{n}}\sum_{A}F^{A}(f,\hat{n})h_{A}(f,\hat{n}),$ (2.21)

the response function $F^{A}(f,\hat{n})=R^{ab}(f,\hat{n})e^{A}_{ab}(\hat{n})$.
The response for any TDI channel (2.15) is

$R_{TDI}^{ab}=\sum_{i}P_{i}R^{ab}\left(f,\hat{n},\hat{u}_{i},\vec{r}_{i}\right),$
(2.22)

where $R^{ab}\left(f,\hat{n},\hat{u}_{i},\vec{r}_{i}\right)$ is the impulse
response of a single arm, $\hat{u}_{i}$ is the direction unit vector of the
arm and $\vec{r}_{i}$ is the midpoint of arm. Here the choose of $\vec{r}_{i}$
is different from the literature for the convenience of calculation, such that
the impulse response become [80]

$R^{ab}(f,\hat{n},\hat{u},\vec{r})=\frac{1}{2}u^{a}u^{b}\mathcal{T}_{\hat{u}}(f,\hat{n}\cdot\hat{u})e^{i2\pi
f\hat{n}\cdot\vec{r}/c},$ (2.23)

where

$\displaystyle\mathcal{T}(f,\hat{n}\cdot\hat{u})$
$\displaystyle\equiv\frac{c}{i2\pi
fL}\frac{1}{1+\hat{n}\cdot\hat{u}}\left[e^{\frac{i\pi
fL}{c}(\hat{n}\cdot\hat{u})}-e^{-\frac{i\pi
fL}{c}(2+\hat{n}\cdot\hat{u})}\right]$ (2.24) $\displaystyle=e^{-\frac{i\pi
fL}{c}}\operatorname{sinc}\left(\frac{\pi
fL}{c}[1+\hat{n}\cdot\hat{u}]\right)$

is the transfer function.

##  3 Correlation analysis

###  3.1 Cross-correlation signal

Usually, the SGWB is very weak, which is masked by the noise of the detector,
and the characteristics are close to the noise. So it is difficult to
distinguish the noise and the GWB signal in a single detector. The correlation
analysis is a powerful method to detect SGWB [74]. We review this method and
apply it to the LISA-TianQin network in this section.

We start with the output signals of the two detectors,

$\displaystyle s_{I}(t)=h_{I}(t)+n_{I}(t),$ (3.1) $\displaystyle
s_{J}(t)=h_{J}(t)+n_{J}(t),$

where $I,J$ respect TianQin and LISA in this paper. And the correlation signal
is

$S_{IJ}=\int_{-T/2}^{T/2}dts_{I}(t)s_{J}(t),$ (3.2)

or in frequency domain

$S_{IJ}=\int_{-\infty}^{\infty}df\int_{-\infty}^{\infty}df^{\prime}\delta_{T}\left(f-f^{\prime}\right)\tilde{s}_{I}(f)\tilde{s}_{J}^{*}\left(f^{\prime}\right),$
(3.3)

where $\delta_{T}(f)=\int_{-T/2}^{T/2}dte^{-i2\pi ft}=\frac{\sin(\pi fT)}{\pi
f}$. Assume that the two detector noises are uncorrelated, the mean of the
correlation signal is

$\mu=\langle
S_{IJ}\rangle=\int_{-\infty}^{\infty}df\int_{-\infty}^{\infty}df^{\prime}\delta_{T}\left(f-f^{\prime}\right)\left\langle\tilde{h}_{I}(f)\tilde{h}_{J}^{*}\left(f^{\prime}\right)\right\rangle.$
(3.4)

Combining eq. (2.5) and eq. (2.21),

$\left\langle\tilde{h}_{I}(f)\tilde{h}_{J}^{*}\left(f^{\prime}\right)\right\rangle=\frac{1}{2}\delta\left(f-f^{\prime}\right)\sum_{A=\\{T,V,B,L\\}}\Gamma_{IJ}^{A}(f)S_{h}^{A}(f),$
(3.5)

where the overlap reduction functions are

$\displaystyle\Gamma^{T}_{IJ}(f)=$ $\displaystyle\frac{1}{8\pi}\int
d^{2}\Omega_{\hat{n}}\sum_{A={+,\times}}F^{A}_{I}(f,\hat{n})F^{A*}_{J}(f,\hat{n}),$
(3.6) $\displaystyle\Gamma^{V}_{IJ}(f)=$ $\displaystyle\frac{1}{8\pi}\int
d^{2}\Omega_{\hat{n}}\sum_{A={X,Y}}F^{A}_{I}(f,\hat{n})F^{A*}_{J}(f,\hat{n}),$
$\displaystyle\Gamma^{B}_{IJ}(f)=$ $\displaystyle\frac{1}{4\pi}\int
d^{2}\Omega_{\hat{n}}F^{B}_{I}(f,\hat{n})F^{B*}_{J}(f,\hat{n}),$
$\displaystyle\Gamma^{L}_{IJ}(f)=$ $\displaystyle\frac{1}{4\pi}\int
d^{2}\Omega_{\hat{n}}F^{L}_{I}(f,\hat{n})F^{L*}_{J}(f,\hat{n}).$

We will calculate the ORFs in next subsection. So the mean of signal is

$\mu=\frac{T}{2}\int_{-\infty}^{\infty}df\sum_{A}\Gamma_{IJ}^{A}(f)S_{h}^{A}(f).$
(3.7)

Assume that the signal is much smaller than the noise, such that

$\left\langle\tilde{s}_{i}^{*}(f)s_{i}^{*}\left(f^{\prime}\right)\right\rangle\approx\left\langle\tilde{n}_{i}^{*}(f)n_{i}^{*}\left(f^{\prime}\right)\right\rangle=\frac{1}{2}\delta\left(f-f^{\prime}\right)P_{i}(f).$
(3.8)

The variance is

$\sigma^{2}=\left\langle S_{IJ}^{2}\right\rangle-\langle
S_{IJ}\rangle^{2}\approx\frac{T}{4}\int_{-\infty}^{\infty}dfP_{I}(f)P_{J}(f).$
(3.9)

So the signal-to-noise ratio (SNR) is

$\rho=\frac{\mu}{\sigma}=\sqrt{T}\frac{\int_{-\infty}^{\infty}df\sum_{A}\Gamma_{IJ}^{A}(f)S_{h}^{A}(f)}{\sqrt{\int_{-\infty}^{\infty}dfP_{I}(f)P_{J}(f)}}.$
(3.10)

A prerequisite for the cross-correlation analysis is that the ORFs remain
constant for the duration of observation. LISA-TianQin does not meet this
condition, so some improvements need to be made. The total observation time
$T$ is divided into $N$ segments, with each segment $\Delta T$. $\Delta T$
need to be greater than the light travel time between the two detectors and
less than the time scale over which the ORF will vary. Therefore, ORFs in each
segment can be treated as constants, labeled as $\Gamma_{IJ,k}^{A}(f)$, and
the cross-correlation analysis could be performance. The SNR for each segment
is

$\rho_{k}=\sqrt{\Delta
T}\frac{\int_{-\infty}^{\infty}df\sum_{A}\Gamma_{IJ,k}^{A}(f)S_{h}^{A}(f)}{\sqrt{\int_{-\infty}^{\infty}dfP_{I}(f)P_{J}(f)}},$
(3.11)

and the total SNR for observation time $T$ reads

$\rho_{total}=\sqrt{\sum_{k=1}^{N}\rho_{k}^{2}}.$ (3.12)

###  3.2 Overlap reduction function

The overlap reduction function can be interpretation as the response of
correlation of two detectors to the isotropic SGWB [68]. The overlap reduction
function mainly depends on three factors: detector similarity, separation and
orientation relative to one another. In the past literature, one usually
considers two identical detectors placed in difference location. In fact, for
detectors like LISA-TianQin, the different arm lengths result in slightly
different frequency bands for their respective responses, which is an
important factor leading to the reduction of ORF. On the other hand, changes
in orbits cause ORFs to change over time. In addition, the small-antenna limit
that applies to ground-based detectors is no longer always applicable to
space-based detectors in the detection band. For example, the characteristic
frequencies $f=c/(2L)$ for LISA and TianQin are $0.06\text{Hz}$ and
$0.86\text{Hz}$ respectively. the small-antenna limit required that the
frequency much smaller than the characteristic frequencies. However, the most
sensitive frequency band of LISA-TianQin is $10^{-3}$Hz – $10^{-1}$ Hz, so the
small antenna limit is not always satisfied. Based on the above
considerations, the ORF of LISA-TianQin network deserves a careful discussion.

For any TDI channel, the ORF for $A$ polarization is

$\displaystyle\Gamma_{TDI}^{A}(f)$ $\displaystyle=\frac{1}{8\pi}\int
d^{2}\Omega_{\hat{n}}\sum_{A}F_{I}^{A}(f,\hat{n})F_{J}^{A^{*}}(f,\hat{n})$
(3.13) $\displaystyle=\frac{1}{8\pi}\sum_{i,j}P_{I,i}P_{J,j}^{*}\int
d^{2}\Omega_{\hat{n}}\sum_{A}R_{I}^{ab}\left(f,n,\hat{u}_{I,i},\vec{r}_{I,i}\right)R_{J}^{cd^{*}}\left(f,n,\hat{u}_{J,j},\vec{r}_{J,j}\right)e_{ab}^{A}e_{cd}^{A}$
$\displaystyle=\frac{e^{\frac{i}{2}\left(\beta_{J}-\beta_{I}\right)}}{32\pi}\sum_{i,j}P_{I,i}P_{J,j}^{*}\hat{u}_{I,i}^{a}\hat{u}_{I,i}^{b}\hat{u}_{J,j}^{c}\hat{u}_{J,j}^{d}\Gamma^{A}_{abcd}\left(\alpha_{ij},\beta_{I},\beta_{J},\hat{u}_{I,i},\hat{u}_{J,j},s_{ij}\right),$

where $\beta_{I}=2\pi fL_{I}/c$ , $\beta_{J}=2\pi fL_{J}/c$ and

$\displaystyle\Gamma^{A}_{abcd}$
$\displaystyle\left(\alpha_{ij},\beta_{I},\beta_{J},\hat{u}_{I,i},\hat{u}_{J,j},\hat{s}_{ij}\right)=\int
d^{2}\Omega_{\hat{n}}e^{-i\alpha_{ij}\hat{n}\cdot\hat{s}_{ij}}$ (3.14)
$\displaystyle\times\sum_{A}\operatorname{sinc}\left(\frac{\beta_{I}}{2}\left[1+\hat{n}\cdot\hat{u}_{I,i}\right]\right)\operatorname{sinc}\left(\frac{\beta_{J}}{2}\left[1+\hat{n}\cdot\hat{u}_{J,j}\right]\right)e_{ab}^{A}e_{cd}^{A}.$

Here $\alpha_{ij}=2\pi fs_{ij}/c$,
$s_{ij}\equiv\left|\Delta\vec{x}_{ij}\right|=\left|\vec{r}_{J,j}^{\prime}-\vec{r}_{I,i}\right|$
and $\hat{s}_{ij}=\Delta\vec{x}_{ij}/s_{ij}$. To keep the definition
consistent, the sum means
$e_{ab}^{+}e_{cd}^{+}+e_{ab}^{\times}e_{cd}^{\times}$ for $A=T$,
$e_{ab}^{X}e_{cd}^{X}+e_{ab}^{Y}e_{cd}^{Y}$ for $A=V$, $2e_{ab}^{B}e_{cd}^{B}$
for $A=B$ and $2e_{ab}^{L}e_{cd}^{L}$ for $A=L$. In this way, the ORF of any
TDI channel can be disassembled and calculated between two separate arms. In
general, the integral in eq. (3.14) can not be calculated analytically. In the
short antenna limit $\beta,\beta^{\prime}\ll 1$, it can be calculated
analytically [81, 74]. The result is

$\displaystyle\Gamma_{abcd}^{A(0)}(\alpha,\hat{s})$
$\displaystyle=A^{A(0)}(\alpha)\delta_{ab}\delta_{cd}+B^{A(0)}(\alpha)\left(\delta_{ac}\delta_{bd}+\delta_{bc}\delta_{ad}\right)$
(3.15)
$\displaystyle+C^{A(0)}(\alpha)\left(\delta_{ab}s_{c}s_{d}+\delta_{cd}s_{a}s_{b}\right)$
$\displaystyle+D^{A(0)}(\alpha)\left(\delta_{ac}s_{b}s_{d}+\delta_{ad}s_{b}s_{c}\right.$
$\displaystyle\left.+\delta_{bc}s_{a}s_{d}+\delta_{bd}s_{a}s_{c}\right)+E^{A(0)}(\alpha)s_{a}s_{b}s_{c}s_{d},$

where the coefficients is

$X^{A(0)}=(M^{(0)})^{-1}Y^{A(0)}.$ (3.16)

Here,

$X^{A(0)}=\left[\begin{array}[]{c}A^{A(0)}\\\ B^{A(0)}\\\ C^{A(0)}\\\
D^{A(0)}\\\
E^{A(0)}\end{array}\right],M^{(0)}=\left[\begin{array}[]{rrrrr}9&6&6&4&1\\\
6&24&4&16&2\\\ 6&4&8&8&2\\\ 4&16&8&24&4\\\ 1&2&2&4&1\\\ \end{array}\right],$
(3.17)

and

$Y^{T(0)}=32\pi\left[\begin{array}[]{c}0\\\ j_{0}(\alpha)\\\ 0\\\
2j_{1}(\alpha)/\alpha\\\ j_{2}(\alpha)/\alpha^{2}\end{array}\right],$ (3.18)
$Y^{V(0)}=32\pi\left[\begin{array}[]{c}0\\\ j_{0}(\alpha)\\\ 0\\\
j_{0}(\alpha)-j_{1}(\alpha)/\alpha\\\
j_{1}(\alpha)/\alpha-4j_{2}(\alpha)/\alpha^{2}\end{array}\right],$ (3.19)
$Y^{B(0)}=32\pi\left[\begin{array}[]{c}j_{0}(\alpha)\\\ j_{0}(\alpha)\\\
2j_{1}(\alpha)/\alpha\\\ 2j_{1}(\alpha)/\alpha\\\
2j_{2}(\alpha)/\alpha^{2}\end{array}\right],$ (3.20)
$Y^{L(0)}=16\pi\left[\begin{array}[]{c}j_{0}(\alpha)\\\ 2j_{0}(\alpha)\\\
2j_{1}(\alpha)/\alpha-2j_{0}(\alpha)\\\
4j_{1}(\alpha)/\alpha-4j_{0}(\alpha)\\\
\left(8/\alpha^{2}-1\right)j_{2}(\alpha)-j_{1}(\alpha)/\alpha\end{array}\right].$
(3.21)

So the ORF of any TDI channel in the small antenna limit is

$\displaystyle\Gamma_{TDI}^{A0}(f)$
$\displaystyle=\frac{e^{\frac{i}{2}\left(\beta_{J}-\beta_{I}\right)}}{32\pi}\sum_{i,j}P_{I,i}P_{J,j}\left\\{A^{A(0)}(\alpha_{ij})\right.$
(3.22)
$\displaystyle+2B^{A(0)}(\alpha_{ij})\left(1+\hat{u}_{I,i}\cdot\hat{u}_{J,j}\right)$
$\displaystyle+C^{A(0)}(\alpha_{ij})\left(\left(\hat{u}_{I,i}\cdot\hat{s}_{ij}\right)^{2}+\left(\hat{u}_{J,j}\cdot\hat{s}_{ij}\right)^{2}\right)$
$\displaystyle+4D^{A(0)}(\alpha_{ij})\left(\hat{u}_{I,i}\cdot\hat{u}_{J,j}\right)\left(\hat{u}_{I,i}\cdot\hat{s}_{ij}\right)\left(\hat{u}_{J,j}\cdot\hat{s}_{ij}\right)$
$\displaystyle\left.+E^{A(0)}(\alpha_{ij})\left(\hat{u}_{I,i}\cdot\hat{s}_{ij}\right)^{2}\left(\hat{u}_{J,j}\cdot\hat{s}_{ij}\right)^{2}\right\\}.$

Notice that there is an extra phase factor
$e^{\frac{i}{2}\left(\beta_{J}-\beta_{I}\right)}$, due to the difference in
the arm lengths of LISA and TianQin. For LISA-TianQin network, the phase
factor $e^{i2\pi f/65\text{mHz}}$ can not ne ignored in the detection band.

Since the small antenna limit may not work well, we extend it slightly. We
expand the eq. (3.14) as a Taylor series of the frequency $f$. To the zeroth
order term is the expression for the small antenna approximation above.
Expanding to the next order, we obtain that

$\displaystyle\Gamma_{abcd}^{A2}$
$\displaystyle\left(\alpha_{ij},\beta_{I},\beta_{J},\hat{u}_{I,i},\hat{u}_{J,j},\hat{s}_{ij}\right)=\Gamma_{abcd}^{A(0)}\left(\alpha_{ij},\hat{s}_{ij}\right)$
(3.23)
$\displaystyle+\Gamma_{abcd}^{A(2)}\left(\alpha_{ij},\beta_{I},\hat{u}_{I,i},\hat{s}_{ij}\right)+\Gamma_{abcd}^{A(2)}\left(\alpha_{ij},\beta_{J},\hat{u}_{J,j},\hat{s}_{ij}\right),$

where

$\displaystyle\Gamma_{dbcd}^{A(2)}(\alpha,\beta,\hat{u},\hat{s})$
$\displaystyle=-\frac{1}{6}\int
d^{2}\Omega_{\hat{n}}\left(\frac{\beta}{2}[1+\hat{n}\cdot\hat{u}]\right)^{2}$
(3.24)
$\displaystyle\times\sum_{A}e_{ab}^{A}e_{cd}^{A}e^{-i\alpha\hat{n}\cdot\hat{s}}.$

Following the same method, we construct it with $\delta_{ab}$, $s_{a}$ and
$u_{a}$ under the premise of ensuring its symmetry. Then solve the linear
equations for the coefficients, we can get the ORF of the second order. The
details of calculation are provided in appendix A.

Figure 3: The tensor ORF of the $X$ channel for LISA-TianQin network. The real
and imaginary parts of ORF are the left and right respectively, and the time
is $t=0$.  Figure 4: The different polarization ORF of the $A$ channel for
LISA-TianQin network. The ORF curve is plotted for that $t=0$. And the real
and imaginary parts are represented by different types of curves. Figure 5:
The comparison of the A channel ORF of two scalar modes for $t=0$. The left-
hand side is the real part, and the right-hand side is the imaginary part.
Figure 6: The ORF at different times. The left and right sides from top to
bottom are the real and imaginary parts of the ORF at $t=0$, $t=1\text{day}$
and $t=\text{yr}/2$ respectively.

In order to quantify the accuracy of the expand ORF, we compared the zero-
order (small antenna approximation) and second-order with numerical
integration, as shown in figure 3, taking the $X$ channel as an example . For
frequencies below the characteristic frequency $f<c/(2L^{LS})=0.06\text{Hz}$,
the second-order expression is more accurate than the zero-order expression,
which is in good agreement with the numerical integration. If the phase factor
$e^{\frac{i}{2}\left(\beta_{J}-\beta_{I}\right)}$ is ignored, the zero-order
accuracy will be even worse. For frequencies above the characteristic
frequency $f>c/(2L^{LS})=0.06\text{Hz}$, the accuracy of second-order is worse
compared with zero-order, since the error of second-order grows faster than
zero-order’s. To improve accuracy, we choose to splice the two at the
characteristic frequency,

$\Gamma^{A}_{TDI}(f)=\begin{cases}\Gamma_{TDI}^{A0}(f)&f\geq c/(2L^{LS})\\\
\Gamma_{TDI}^{A2}(f)&f<c/(2L^{LS})\end{cases}.$ (3.25)

They are smooth at the connection point, both equal to 0. This concatenated
expression has high enough precision for SGWB data analysis. The detector is
insensitive to frequency bands above the characteristic frequency due to loud
noise and low response for SGWB. Also, even if the numerical integration
method is used, the calculation for frequencies above the characteristic
frequency is very slow and the result is inaccurate, due to the rapid
oscillation caused by the high frequency.

In figure 4, we show the ORF of $A$ channel for different polarizations. We
choose not to normalize the ORF, because normalization is thankless for the
TDI combination of the LISA-TianQin network. Unlike ground-based detectors,
TDI combinations are insensitive to low frequencies because of the virtual
equal arm interferometric measurements. Actually, the ORF should be
proportional to $\beta_{I}\beta_{J}$. On the other hand, it is not trivial to
define them co-located for any TDI channel, since LISA and TianQin have
different arm lengths. Furthermore, the ORF of different polarizations also
share some common characteristics. Their zero point distributions are similar:
$f\approx nc/(2|\Delta\vec{x}|)$, $f=nc/(2L_{I})$ and $f=nc/(2L_{J})$. Due to
the different arm lengths of the two detectors, TianQin is still sensitive in
the frequency band $f>c/(2L^{LS})\approx 60\text{mHz}$, so the ORF does not
decay rapidly to 0 beyond $f=c/(2L^{LS})\approx 60\text{mHz}$. The peak in the
frequency band of 60mHz to 860mHz is mainly contributed by the response of
Tianqin, although it has virtually no contribution to the SNR due to the
intense noise. As the frequency continues to increase, the ORF will continue
to decay, which is determined by the properties of the spherical Bessel
function.

It is well acknowledged that the angular antenna response of the Michelson
interferometer-type detector to breathing and longitudinal modes is degenerate
at low frequency approximations, namely,
$F^{L}(f,\hat{n})=\sqrt{2}F^{B}(f,\hat{n})$. This results in the degeneracy of
the ORFs for the two scalar polarizations in the low frequency range, with a
quantitative relationship of $\Gamma^{L}=2\Gamma^{B}$. However, as the
frequency increases, the degeneracy between them is removed as shown in figure
5. This implies that it is possible to resolve two scalar modes through the
LISA-TianQin network. The key point is that the degeneracy is broken, not how
different the ORFs of the two scalar modes are, as distinguishing between the
two modes requires different correlation channels to provide more information.
As the quantitative relationship of $\Gamma^{L}=2\Gamma^{B}$ always remains in
the low frequency range, it does not provide any information to distinguish
between the two scalar modes, regardless of the number of correlated channels
added. In the high frequency range, the relationship of
$\Gamma^{L}=2\Gamma^{B}$ is no longer maintained, which means that the ratio
of $\Gamma^{L}$ and $\Gamma^{B}$ for different correlation channels may vary.
This provides the possibility of distinguishing between breathing and
longitudinal modes.

What more, we plot the ORF for different polarizations at different times in
figure 6. The time-varying property makes data analysis more difficult, but
also increases the chance of resolving different polarizations. Correlations
at different locations can be used as multi-correlation channels to
distinguish between the four polarizations. Based on this idea, an algorithm
can be developed to achieve the separation of different polarizations. And the
accuracy of a polarization identification depends on how different the rate of
change in time of its ORF is from other modes’. Details are provided in
section 5.

##  4 Sensitivity for the background of alternative polarizations

Figure 7: The SNR of power-law SGWB with different amplitude $\Omega^{A}$ and
index $\alpha^{A}$ for different polarizations. The reference frequency is
choose to be $f_{0}=1\text{mHz}$. The A channel is chosen for calculation here
and the rest of the article.  Figure 8: The normalized SNR as a function of
the cutoff frequency.

###  4.1 Optimal filter

In addition to the naive method of directly correlating the output of two
detecters in section 3, the idea of matched filter is often used to improve
the SNR. Matched filtering starts by multiplying the correlation signal (3.3)
by a filter function,

$S_{m}=\int^{\infty}_{-\infty}df\int^{\infty}_{-\infty}df^{\prime}\delta_{T}(f-f^{\prime})\tilde{s}_{I}(f)\tilde{s}^{*}_{J}(f^{\prime})Q(f).$
(4.1)

The SNR is [74]

$\rho_{m}^{2}=T\left(\frac{3H_{0}^{2}}{2\pi^{2}}\right)^{2}\frac{(Q,\frac{\sum_{A}\Gamma_{IJ}^{A}(f)\Omega^{A}_{h}(f)}{f^{3}P_{I}(f)P_{J}(f)})^{2}}{(Q,Q)},$
(4.2)

where the inner product of $A$ and $B$ is defined as
$(A,B)\equiv\int_{-\infty}^{\infty}dfA^{*}(f)B(f)P_{I}(f)P_{J}(f)$. And the
optimal filter is

$Q(f)\propto\frac{\sum_{A}\Gamma_{IJ}^{A}(f)\Omega^{A}_{h}(f)}{f^{3}P_{I}(f)P_{J}(f)},$
(4.3)

when mean that the model and the actual energy-density spectrum match. The
resulting optimal SNR is given by

$\rho_{o}\approx\sqrt{T}\frac{3H_{0}^{2}}{2\pi^{2}}\left[\int^{\infty}_{-\infty}df\frac{\left|\sum_{A}\Gamma_{IJ}^{A}(f)\Omega^{A}_{h}(f)\right|^{2}}{f^{6}P_{I}(f)P_{J}(f)}\right]^{1/2}.$
(4.4)

The premise of the above expression is that the noise amplitude is much larger
than the signal. In general, it will be modified to [82]

$\rho_{o}=\sqrt{T}\frac{3H_{0}^{2}}{2\pi^{2}}\left[\int^{\infty}_{-\infty}df\frac{\left|\sum_{A}\Gamma_{IJ}^{A}(f)\Omega^{A}_{h}(f)\right|^{2}}{f^{6}M(f)}\right]^{1/2},$
(4.5)

where

$\displaystyle M(f)$
$\displaystyle=P_{I}(f)P_{J}(f)+S_{h}(f)\left(P_{I}(f)\Gamma_{JJ}(f)+P_{J}(f)\Gamma_{II}(f)\right)$
(4.6)
$\displaystyle+S^{2}_{h}(f)\left(\Gamma^{2}_{IJ}(f)+\Gamma_{II}(f)\Gamma_{JJ}(f)\right).$

Once again, considering the time-varying effect of LISA-TianQin, the total SNR
after segmented processing is

$\rho_{o,total}=\sqrt{T}\frac{3H_{0}^{2}}{2\pi^{2}}\left[\int^{\infty}_{-\infty}df\frac{\left|\sum_{A}\tilde{\Gamma}_{IJ}^{A}(f)\Omega^{A}_{h}(f)\right|^{2}}{f^{6}M(f)}\right]^{1/2},$
(4.7)

where the average ORFs is defined as

$\tilde{\Gamma}_{IJ}^{A}(f)=\frac{1}{N}\sqrt{\sum_{k=1}^{N}|\Gamma_{IJ,k}^{A}(f)|^{2}}.$
(4.8)

Assuming that the spectrum of SGWB is power-law, the optimal SNR calculated by
eq. (4.4) can be used to evaluate the detection capability of LISA-TianQin
network for SGWB. We show the SNR of SGWB with different polarizations in
figure 7, where the observation time is chosen as $T=1\text{yr}$ and the TDI
channel is selected as channel A here and the rest of the article. The result
shows that the LISA-TianQin network is more sensitive to tensor and vector
polarizations than the two scalar polarizations. For example, for the same
spectral index $\alpha^{A}=2/3$, in order to achieve a SNR of 10, the
amplitudes at reference frequency for different polarizations are required to
be: $\Omega_{0}^{T}=8.59\times 10^{-12}$, $\Omega_{0}^{V}=7.58\times
10^{-12}$, $\Omega_{0}^{B}=1.71\times 10^{-11}$ and $\Omega_{0}^{L}=8.61\times
10^{-12}$ . In addition, in order to find the most sensitive frequency band of
LISA-TianQin, we define a normalized SNR as a function of the cutoff
frequency,

$\hat{\rho}_{o}(f_{\text{cut}})=\frac{1}{\rho_{o}}\sqrt{T}\frac{3H_{0}^{2}}{2\pi^{2}}\left[2\int^{f_{\text{cut}}}_{0}df\frac{\left|\sum_{A}\Gamma^{A}(f)\Omega^{A}_{h}(f)\right|^{2}}{f^{6}P_{I}(f)P_{J}(f)}\right]^{1/2}.$
(4.9)

We show the function $\hat{\rho}_{o}(f_{\text{cut}})$ in figure 8 where the
power-law SGWB with spectrum index $\alpha_{A}=2/3$ are assumed. The results
show that the frequency band that contributes the most to the SNR of
background from compact binary coalescences is $1\text{mHz}-10\text{mHz}$. The
contribution of high-frequency bands to the SNR is negligible, with each
polarization SNR of more than 0.06Hz accounting for $8.5\times 10^{-8}$,
$6.5\times 10^{-8}$, $7.7\times 10^{-8}$ and $1.1\times 10^{-7}$,
respectively. This reflects the fact that strong noise prevents high-frequency
bands from being accurately detected at present. The error of the ORF in high
frequency region does not matter. Of course, as detector technology advances,
noise levels decrease, the development of high-precision ORFs in high
frequency region is necessary to detect potential high-frequency sources. Also
note that the curves for the scalar-breathing and the scalar-longitudinal mode
are the same, since their responses are the same at low frequency
approximations. If we increase the spectrum index, which means that the high
frequency part of the SGWB is louder, the curves for the scalar-breathing and
the scalar-longitudinal mode will be differentiated.

###  4.2 Parameter estimation accuracy

Figure 9: The uncertainties of the amplitude $\Omega^{A}$ and index
$\alpha^{A}$ for different polarizations. The uncertainties of the index
$\alpha^{A}$ are calculated with suitable amplitude parameters such that the
SNR is 10.

Under the premise that the model is accurate, the optimal filter can not only
improve the signal-to-noise ratio, but also provide an estimate of the model
parameters. When the SNR is high enough, the Fisher information matrix (FIM)
can be used to evaluate parameter estimation accuracy [83]. The FIM is defined
by

$F_{ab}=T\int^{\infty}_{-\infty}df\frac{(\sum_{A}\Gamma^{A}(f)\frac{\partial
S^{A}_{h}(f)}{\partial\theta_{a}})(\sum_{A}\Gamma^{A*}(f)\frac{\partial
S^{A}_{h}(f)}{\partial\theta_{b}})}{P_{I}(f)P_{J}(f)},$ (4.10)

and the estimate error of a parameter $\Delta\theta_{a}$ is estimated from the
inverse of FIM

$\Delta\theta_{a}=\sqrt{F_{aa}^{-1}}.$ (4.11)

There are a total of eight parameters we need to estimate for a power law
model that includes all polarizations. And the parameter estimation errors are
given by

$\displaystyle\Delta\Omega^{A}_{0}$
$\displaystyle=\frac{1}{\sqrt{T}}\frac{2\pi^{2}}{3H_{0}^{2}}\left[\int^{\infty}_{-\infty}df\frac{|\Gamma^{A}(f)|^{2}(f/f_{0})^{2\alpha_{A}}}{f^{6}P_{I}(f)P_{J}(f)}\right]^{-1/2},$
(4.12) $\displaystyle\Delta\alpha_{A}$
$\displaystyle=\frac{1}{\sqrt{T}}\frac{2\pi^{2}}{3H_{0}^{2}}\frac{1}{\Omega^{A}_{0}}\left[2\int^{\infty}_{-\infty}df\frac{|\Gamma^{A}(f)|^{2}(f/f_{0})^{2\alpha_{A}}\ln(f/f_{0})}{f^{6}P_{I}(f)P_{J}(f)}\right]^{-1/2}.$

The uncertainties of the amplitude $\Omega^{A}$ and index $\alpha^{A}$ for
different polarizations are shown in figure 9. In general, the accuracies of
the amplitude parameters for tensor and vector are higher than the two scalar
polarizations. This means they are easier to detect. The uncertainties of the
index are plotted with suitable amplitude parameters such that the SNR is 10.
In such a case, the spectral parameter estimation accuracy of the two scalar
modes is higher. Another point worth noting is that the uncertainty curves of
the scalar-breathing mode and the scalar-longitudinal mode almost overlap
because their ORFs are degenerate at low frequencies. And when the exponent
increases, which means that the high frequency signal is stronger, the curves
are not overlapping, because the degeneracy breaks at high frequency.

##  5 Polarization separation

If the alternative polarization really exists, different polarization modes
will be mixed in the cross-correlation signal. We need suitable methods to
distinguish difference modes in the cross-correlation signal. One possible
approach is Bayesian model selection proposed in [23], which is adopted to
detect alternative polarization background in the cross-correlation data of
ground-based detecter [55, 56, 57]. This approach is efficient and easy to
implement. However, it may be prone to bias if the model does not fit well
with the true background.

If one want to get away from being bound by model assumptions, there is a
intuitive method to separate polarizations using multiple detecters pairs
[58]. Different detector pairs have different ORFs, which breaks the
degeneracy between different polarizations and provides enough degrees of
freedom to separate different polarizations. Notably, this approach can be
implemented independently for each frequency bin. The accuracy of the results
for each frequency bin depends on the corresponding signal and noise strength.
However, this method requires a higher signal-to-noise than the previous
method. Another disadvantage is that multiple detector pairs naturally mean
that more than two detectors are required.

Since the LISA-TianQin network varies with orbits, their ORFs for difference
polarizations varies accordingly as see in section 3. So their correlation
signals at different positions respond differently to different polarizations,
which means that it can be equivalently regarded as different detector pairs
at different times. There is an opportunity to extract different polarization
patterns from the data from the two detectors, using a similar approach to the
multi-detector based approach [58].

The key point is to break the coupling of different polarizations in the data
through multiple correlation channels. If data for $n$ detectors is available,
a total of $N=n(n-1)/2$ correlation channels can be constructed as

$Z_{ij}(f)\propto\tilde{s}_{i}(f)\tilde{s}^{*}_{j}(f),$ (5.1)

where $1\leq i<j\leq n$. These statistics satisfy

$\left\langle\bm{Z(f)}\right\rangle=\bm{M}\begin{bmatrix}\Omega^{T}(f)\\\
\Omega^{V}(f)\\\ \Omega^{B}(f)\\\ \Omega^{L}(f)\end{bmatrix},$ (5.2)

where $\bm{Z(f)}=[Z_{12}(f),Z_{13}(f),...,Z_{(n-1)n}(f)]^{T}$ and

$\bm{M}=\begin{bmatrix}\Gamma^{T}_{12}&\Gamma^{V}_{12}&\Gamma^{B}_{12}&\Gamma^{L}_{12}\\\
\Gamma^{T}_{13}&\Gamma^{V}_{13}&\Gamma^{B}_{13}&\Gamma^{L}_{13}\\\
...&...&...&...\\\
\Gamma^{T}_{(n-1)n}&\Gamma^{V}_{(n-1)n}&\Gamma^{B}_{(n-1)n}&\Gamma^{L}_{(n-1)n}\\\
\end{bmatrix}$ (5.3)

is the detecter correlation matrix. If the rank of the matrix is greater than
or equal to 4, $\text{rank}(\bm{M})\geq 4$, the polarization components can be
separated by inverting eq. (5.2), namely,

$\begin{bmatrix}\hat{\Omega}^{T}(f)\\\ \hat{\Omega}^{V}(f)\\\
\hat{\Omega}^{B}(f)\\\
\hat{\Omega}^{L}(f)\end{bmatrix}=\bm{M(f)}^{-1}\left\langle\bm{Z(f)}\right\rangle.$
(5.4)

In this way, the spectrum of each polarization component can be reconstructed
from the detector data by eq. (5.4). However, the degeneracy of ORFs of two
scalar polarizations for ground-based causes the rank of the detector
correlation matrix to be at most 3. So the SGWB detected by the ground-based
detector network can only be separated into tensor, vector and scalar
components.

Two space-based detectors may achieve the capabilities of the above-mentioned
ground-based detector network to distinguish polarizations in SGWB. It is only
necessary to replace the different correlation channels with the correlation
of the two space-based detectors at different positions. Below we use another
description of maximum likelihood estimation mentioned in [68] to obtain the
same result. First, the data is divided into $N$ segments (indexed by $i$),
and the duration of each segment of data $\Delta T$ satisfies the condition
that it is greater than the light travel time between the two detectors and
less than the time scale of the ORF changes. A cross-correlation statistic can
be constructed as

$\hat{C}_{i}(f)=\frac{2}{\Delta
T}\frac{2\pi^{2}}{3H_{0}^{2}}f^{3}\tilde{s}_{1,i}(f)\tilde{s}^{*}_{2,i}(f),$
(5.5)

normalized such that the statistic’s mean is

$\left\langle\hat{C}_{i}(f)\right\rangle=\sum_{A}\Gamma_{i}^{A}(f)\Omega^{A}(f),$
(5.6)

where $\Gamma_{i}^{A}(f)$ is the ORF for i-th segment data. The statistic
$\bm{\hat{C}(f)}=[\hat{C}_{1}(f),\hat{C}_{2}(f),...,\hat{C}_{N}(f)]^{T}$ acts
as the role of different correlation channels $\bm{Z(f)}$. And the variance is

$\sigma_{i}^{2}(f)=\frac{1}{2\Delta T\Delta
f}\left(\frac{2\pi^{2}}{3H_{0}^{2}}\right)^{2}f^{6}P_{1,i}(f)P_{2,i}(f)$ (5.7)

where $\Delta f$ is the frequency bin width, $P_{1,i}(f)$ and $P_{2,i}(f)$ are
the noise power spectral density of detectors.

The likelihood function for $\bm{\hat{C}(f)}$ is

$\mathcal{L}[\bm{\hat{C}(f)}|\mathcal{A}]\propto\text{exp}\left[-\sum_{i}^{N}\sum_{f}\frac{|\hat{C}_{i}(f)-\sum_{A}\Gamma_{i}^{A}(f)\Omega^{A}(f)|^{2}}{2\sigma_{i}^{2}(f)}\right],$
(5.8)

where
$\mathcal{A}=[\Omega^{T}(f),\Omega^{V}(f),\Omega^{B}(f),\Omega^{L}(f)]^{T}$
represent the fractional energy density of different polarizations.
Equivalently, we can write it as

$\mathcal{L}[\bm{\hat{C}(f)}|\mathcal{A}]\propto\exp\left[-\frac{1}{2}(\bm{\hat{C}}-\bm{M}\mathcal{A})^{\dagger}\mathcal{N}^{-1}(\bm{\hat{C}}-\bm{M}\mathcal{A})\right],$
(5.9)

where the detector correlation matrix becomes

$\bm{M}=\begin{bmatrix}\Gamma^{T}_{1}&\Gamma^{V}_{1}&\Gamma^{B}_{1}&\Gamma^{L}_{1}\\\
\Gamma^{T}_{2}&\Gamma^{V}_{2}&\Gamma^{B}_{2}&\Gamma^{L}_{2}\\\
...&...&...&...\\\
\Gamma^{T}_{N}&\Gamma^{V}_{N}&\Gamma^{B}_{N}&\Gamma^{L}_{N}\\\ \end{bmatrix}.$
(5.10)

And the noise covariance matrix is
$\mathcal{N}_{ij}=\delta_{ij}\sigma_{i}^{2}(f)$.

Figure 10: The efficient overlap functions for different polarizations.
Figure 11: Equivalent sensitivity curves of different polarization components.
Figure 12: The reconstructed spectrum for each polarizations by the equivalent
multi-detector-based method. The data contains tensor and breathing
polarisations: $\Omega^{T}=8.59\times 10^{-12}(f/1\text{mHz})^{2/3}$ and
$\Omega^{B}=1.71\times 10^{-11}(f/1\text{mHz})^{2/3}$. The red dashed line
represents the actual injected spectrum, and the light blue solid line
represents the reconstructed spectrum.

The maximum-likelihood estimator for the fractional energy density of
different polarizations are [68]

$\hat{\mathcal{A}}=F^{-1}X,$ (5.11)

where $F=M^{\dagger}\mathcal{N}^{-1}M$, $X=M^{\dagger}\mathcal{N}^{-1}\hat{C}$
are the Fisher matrices and ’dirty’ map for this analysis. The maximum-
likelihood estimator $\hat{\mathcal{A}}$ is equivalent to eq. (5.4), which is
easier to calculate in practice and avoids finding the inverse of a huge
matrix. After the data is available, the reconstructed spectrum of each
polarization component can be obtained by calculating eq. (5.11) for each
frequency. The error will be determined by the detector noise at that
frequency and the degeneracy of ORFs of different polarizations. Specifically,
the marginalized uncertainties of the maximum-likelihood estimates are given
by the inverse of the Fisher matrices,

$\sigma^{2}_{\Omega^{A}}=\left\langle\mathcal{A}\mathcal{A}^{\dagger}\right\rangle_{AA}=(F^{-1})_{AA}.$
(5.12)

Without loss of general, we assume that the noise power spectrum for different
time segments are equal. The Fisher matrix is

$F_{AB}=\sum_{i}\Gamma_{i}^{A}\Gamma_{i}^{B*}/\sigma_{i}^{2}.$ (5.13)

And the efficient overlap function can be defined by

$\Gamma_{\text{eff}}^{A}(f)\equiv\sigma^{-1}_{\Omega^{A}}\sigma_{i}=\frac{\sigma_{i}}{\sqrt{(F^{-1})_{AA}}},$
(5.14)

which can reflect the intensity of the response to different polarization
components. The efficient overlap functions for LISA-TianQin network are shown
in figure 10. The effective overlap functions of the two scalar components are
significantly lower than that of tensors and vectors, due to the similar ORFs
of scalar-breathing and scalar-longitudinal modes. In order to quantify the
sensitivity of this method to different polarization components, detector
noise should also be included. This is the role of uncertainties
$\sigma_{\Omega^{A}}$, which are shown in figure 11. For a certain frequency
$f$, the intensity of $A$ component must meet
$\Omega^{A}(f)>\sigma_{\Omega^{A}}(f)$ to be recognized.

To visually measure the method’s ability to distinguish between different
polarizations, an overall SNR can be defined as

$\displaystyle\text{SNR}_{\text{eff}}^{A}$
$\displaystyle=\left[\int_{0}^{\infty}df\frac{|\Gamma^{A}_{\text{eff}}(f)\Omega_{h}^{A}(f)|^{2}}{\sigma_{i}^{2}}\right]^{1/2}$
(5.15) $\displaystyle=\sqrt{\Delta
T}\frac{3H_{0}^{2}}{2\pi^{2}}\left[2\int^{\infty}_{0}df\frac{\left|\Gamma_{\text{eff}}^{A}(f)\Omega^{A}_{h}(f)\right|^{2}}{f^{6}P_{I}(f)P_{J}(f)}\right]^{1/2}.$

By convention, a SNR of 10 is considered to be a signal extract from the data.
For an astronomical background spectrum with index $\alpha=2/3$, the
amplitudes at reference frequency, $f_{0}=1\text{mHz}$, for different
polarizations are required to be: $\Omega_{0}^{T}=1.27\times 10^{-11}$,
$\Omega_{0}^{V}=9.55\times 10^{-12}$, $\Omega_{0}^{B}=1.43\times 10^{-10}$ and
$\Omega_{0}^{L}=7.23\times 10^{-11}$ . Comparing with the results of optimal
SNR, we find that the ability of the method to resolve vector polarization is
the best, followed by tensor. The ability to distinguish two scalars is the
worst, and the SNR is equivalent to 12% of the optimal SNR. For comparison,
for the flat spectral index $\alpha=0$, the amplitudes for different
polarizations need to exceed: $\Omega_{0}^{T}=2.90\times 10^{-11}$,
$\Omega_{0}^{V}=2.59\times 10^{-11}$, $\Omega_{0}^{B}=3.61\times 10^{-10}$ and
$\Omega_{0}^{L}=1.84\times 10^{-10}$ .

To verify the ability of the above method in resolving polarization modes, we
apply it to the simulated data of different situations. Notably, we simulated
a mixture of tensor and breathing modes to examine whether it can distinguish
between breathing and longitudinal modes. We simulate a year of data and
divided it into 3600 segments, and the noise is generated according to eq.
(5.7). The spectrum injected into each piece of data are
$\Omega^{T}=8.59\times 10^{-12}(f/1\text{mHz})^{2/3}$ and
$\Omega^{B}=1.71\times 10^{-11}(f/1\text{mHz})^{2/3}$. The parameters are
choose such that individual tensor and scalar-breathing components both
correspond to an optimal SNR of 10, mixed together with an optimal SNR of 15.
The spectrum for each polarization reconstructed with eq. (5.11) is shown in
figure 12. The spectral accuracy of the reconstructed spectrum in the
frequency domain $10^{-3}\text{Hz}-10^{-2}\text{Hz}$ is relatively high. It
can be concluded that in the most sensitive frequency band, the different
modes can indeed be distinguished. Even two scalar modes are capable of
distinguishing, but the accuracy is worse than tensor or vector modes.

Figure 13: The reconstructed spectrum for each polarizations by the equivalent
multi-detector-based method in the case of all four modes. The spectra of each
polarization are: $\Omega^{T}=1.42\times 10^{-11}(f/1\text{mHz})^{2/3}$,
$\Omega^{V}=1.62\times 10^{-11}$, $\Omega^{B}=1.71\times
10^{-11}(f/1\text{mHz})^{2/3}$ and $\Omega^{L}=8.61\times
10^{-12}(f/1\text{mHz})^{2/3}$.  Figure 14: The comparison of the
reconstructed spectrum with observation time $T=1\text{years}$ and
$T=5\text{years}$. The injected spectra are the same as in figure 13.

As discussed above, the separate breathing mode of the two scalar modes will
not be mistaken for longitudinal mode. Furthermore, some alternative theories
of gravity would predict more complex mixed polarizations. On the other hand,
the spectral index for different polarizations may be different. Therefore, we
simulate a more complex case with all four polarizations at the same time, in
which the spectral index of the vector mode is set to distinguish it from the
other three. The special spectrum injected into the data are:
$\Omega^{T}=1.42\times 10^{-11}(f/1\text{mHz})^{2/3}$, $\Omega^{V}=1.62\times
10^{-11}$, $\Omega^{B}=1.71\times 10^{-11}(f/1\text{mHz})^{2/3}$ and
$\Omega^{L}=8.61\times 10^{-12}(f/1\text{mHz})^{2/3}$. The indices for tensor,
breathing and longitudinal modes are set to $2/3$, while the spectra of vector
mode is set to be flat, namely, $\alpha_{V}=0$. As shown in figure 13, the
spectra of each polarization are correctly reconstructed. The noise of the
breathing and longitudinal modes is significantly greater than that of the
other two modes, but each can be identified independently. The reconstructed
polarization spectrum fall exactly near the injected true value, and the
difference between the polarizations is only the strength of noise.

figure 14 shows the result of increasing the observation time for the same
signal, where the parameters of the signal are exactly the same as figure 13,
and the observation time is increased to 5 years. The SNR increases with
observation time, specifically $\text{SNR}\propto\sqrt{T}$. Equivalently, five
times the observation time would reduce the uncertainty by a factor of 2.2.
This behavior is clearly presented in figure 14. The noise of the
reconstruction spectrum observed over five years is significantly reduced for
the same signal.

The simulation results verify the ability of the equivalent multi-detector-
based method to separate different polarizations. The reconstructed spectrum
for each polarization obtain by this method falls near the true value, and the
uncertainty depends on the noise of the detector and the correlation between
one component and other polarization component in the correlation matrix in
eq. (5.10). It may be misunderstood that this correlation between different
polarized components refers to their ORF similarity. In fact, it does not
matter how many times the ORF for one mode is that for the other, what matters
is how different the values of each of them change at different positions in
orbit. For example, the degeneracy of breathing and longitudinal modes in the
low frequency band results in the ORF for longitudinal mode being twice that
of the breathing mode. But the key point is that this double relationship is
always maintained in the low frequency band, regardless of how the detectors
move. For high frequency bands, the degenerate relation is broken, and the
ORFs of breathing and longitudinal modes have the opportunity to vary with
time differently. The time-varying effect has been encapsulated in the
equivalent ORFs $\Gamma_{\text{eff}}^{A}(f)$ constructed by the inverse of
correlation matrix, which could quantifies the ability to separate
polarization.

##  6 Discussion

In this paper, we studied the detectability of alternative polarizations of
SGWB with LISA-TianQin network. The different orbits of LISA and TianQin make
them very different from ground-based detectors. The relative orientation of
their orbital planes changes with time, which means that the ORF will vary
accordingly. In other words, the SGWB signal in their cross-correlation signal
varies over time. This will pose certain challenges for data analysis. On the
other hand, LISA-TianQin has advantages for resolving polarization modes in
SGWB in principle.

TDI technology is applied to space gravitational wave detectors to suppress
laser phase noise. Based on the small antenna approximation, we obtain the ORF
of the second-order expansion of the frequency for any TDI channel. This
method can be applied to any laser interferometric detector, even with
different configurations and different arm lengths. For the LISA-TianQin
network, the accuracy can be effectively improved in its most sensitive
frequency band. Then we study the detectability of alternative polarizations
in SGWB with LISA-TianQin network. We calculated the signal-to-noise ratio and
parameter estimation accuracy of power-law SGWB with different polarization
modes.

Once the SGWB is detected, it is necessary to distinguish different
polarization modes from it. An equivalent multi-detector-based approach can be
applied on the LISA-TianQin network, thanks to its special orbital variation.
The ORFs of the LISA-TianQin network are different at different times, so they
can be equivalently regarded as different detectors. In fact, the cross-
correlated signals at different times form a system of linear equations for
different polarization modes. In turn, the system of equations can be solved
to obtain different polarization modes, including the influence of noise of
course. Although this method does not rely on model assumptions, it requires a
high signal-to-noise ratio. Its resolution at a certain frequency bar is
affected by the degeneracy of the ORFs for different polarizations in addition
to noise at that frequency bar. what’s more, it has more advantages in
distinguishing scalar-breathing mode and scalar-longitudinal mode.

Our proposed method for resolving polarization relies on the change of
relative orientation of the two detectors over time, with LISA-TianQin being
the most pronounced due to differences in orbital configuration. However, our
analysis can be extended to include any space-based detector. We will compare
the resolution capabilities of different detector pairs and evaluate the
ability improvement of multi-detector network cooperation as future work.

##  Appendix A Second order overlap reduction function of any TDI channel

The ORF frequency second-order expansion of any TDI channel requires the
calculation of coefficients

$\Gamma_{abcd}^{A(2)}(\alpha,\beta,\hat{u},\hat{s})=-\frac{1}{6}\int
d^{2}\Omega_{\hat{n}}\left(\frac{\beta}{2}[1+\hat{n}\cdot\hat{u}]\right)^{2}\sum_{A}e_{ab}^{A}e_{cd}^{A}e^{-i\alpha\hat{n}\cdot\hat{s}}.$
(A.1)

According to the symmetry of the index, we construct it with $\delta_{ab}$,
$s_{a}$ and $u_{a}$:

$\displaystyle\Gamma_{abcd}^{A(2)}$
$\displaystyle(\alpha,\beta,\hat{u})=A^{A(2)}\delta_{ab}\delta_{cd}+B^{A(2)}\left(\delta_{ac}\delta_{bd}+\delta_{bc}\delta_{ad}\right)+C_{1}^{A(2)}\left(\delta_{ab}s_{c}s_{d}+\delta_{cd}s_{a}s_{b}\right)$
(A.2)
$\displaystyle+C_{2}^{A(2)}\left(\delta_{ab}u_{c}u_{d}+\delta_{cd}u_{a}u_{b}\right)+D_{1}^{A(2)}\left(\delta_{ac}s_{b}s_{d}+\delta_{ad}s_{b}s_{c}+\delta_{bc}s_{a}s_{d}+\delta_{bd}s_{a}s_{c}\right)$
$\displaystyle+D_{2}^{A(2)}\left(\delta_{ac}u_{b}u_{d}+\delta_{ad}u_{b}u_{c}+\delta_{bc}u_{a}u_{d}+\delta_{bd}u_{a}u_{c}\right)+E_{1}^{A(2)}s_{a}s_{b}s_{c}s_{d}+E_{2}^{A(2)}u_{a}u_{b}u_{c}u_{d}$
$\displaystyle+E_{3}^{A(2)}\left(u_{a}u_{b}s_{c}s_{d}+s_{a}s_{b}u_{c}u_{d}\right)+E_{4}^{A(2)}\left(u_{a}u_{c}s_{b}s_{d}+u_{a}u_{d}s_{b}s_{c}+u_{b}u_{c}s_{a}s_{d}+u_{b}u_{d}s_{a}s_{c}\right).$

Contracting the expression with $\delta^{ab}\delta^{cd}$,
$\left(\delta^{ac}\delta^{bd}+\delta^{bc}\delta^{ad}\right)$,$\dots$,
$\left(u^{a}u^{c}s^{b}s^{d}+u^{a}u^{d}s^{b}s^{c}+u^{b}u^{c}s^{a}s^{d}+u^{b}u^{d}s^{a}s^{c}\right)$,
the linear system of equations for the coefficients $A^{A(2)}$, $B^{A(2)}$,
$\dots$, $E_{4}^{A(2)}$ is obtained, which can be expressed as matrix
equations as

$M^{(2)}X^{A(2)}=Y^{A(2)}.$ (A.3)

Here

$M^{(2)}=\left[\begin{array}[]{cccccccccccc}9&6&6&6&4&4&1&1&2&4\eta^{2}\\\
6&24&4&4&16&16&2&2&4\eta^{2}&4\left(1+\eta^{2}\right)\\\
6&4&8&2+6\eta^{2}&8&8\eta^{2}&2&2\eta^{2}&2\left(1+\eta^{2}\right)&8\eta^{2}\\\
6&4&2+6\eta^{2}&8&8\eta^{2}&8&2\eta^{2}&2&2\left(1+\eta^{2}\right)&8\eta^{2}\\\
4&16&8&8\eta^{2}&24&4+20\eta^{2}&4&4\eta^{2}&8\eta^{2}&4+12\eta^{2}\\\
4&16&8\eta^{2}&8&4+20\eta^{2}&24&4\eta^{2}&4&8\eta^{2}&4+12\eta^{2}\\\
1&2&2&2\eta^{2}&4&4\eta^{2}&1&\eta^{4}&2\eta^{2}&4\eta^{2}\\\
1&2&2\eta^{2}&2&4\eta^{2}&4&\eta^{4}&1&2\eta^{2}&4\eta^{2}\\\
2&4\eta^{2}&2\left(1+\eta^{2}\right)&2\left(1+\eta^{2}\right)&8\eta^{2}&8\eta^{2}&2\eta^{2}&2\eta^{2}&2\left(1+\eta^{4}\right)&8\eta^{2}\\\
4\eta^{2}&4\left(1+\eta^{2}\right)&8\eta^{2}&8\eta^{2}&4+12\eta^{2}&4+12\eta^{2}&4\eta^{2}&4\eta^{2}&8\eta^{2}&4\left(1+\eta^{2}\right)^{2}\end{array}\right],$
$X^{A(2)}=\left[A^{A(2)},B^{A(2)},C_{1}^{A(2)},C_{2}^{A(2)},D_{1}^{A(2)},D_{2}^{A(2)},E_{1}^{A(2)},E_{2}^{A(2)},E_{3}^{A(2)},E_{4}^{A(2)}\right]^{-1},$

where $\eta=\hat{u}\cdot\hat{s}$. And the coefficients in $Y^{A(2)}$ can be
calculated by contracting eq. (A.1) with $\delta^{ab}\delta^{cd}$,
$\left(\delta^{ac}\delta^{bd}+\delta^{bc}\delta^{ad}\right)$,$\dots$ For
example,
$Y^{T(2)}(1)=-\frac{1}{6}\int
d^{2}\Omega_{\hat{n}}\left(\frac{\beta}{2}[1+\hat{n}\cdot\hat{u}]\right)^{2}\sum_{A=+,\times}e_{ab}^{A}e_{cd}^{A}\delta^{ab}\delta^{cd}e^{-i\alpha\hat{n}\cdot\hat{s}}=0$.
The solution of coefficient $X^{A(2)}$ is

$X^{A(2)}=(M^{(2)})^{-1}Y^{A(2)}.$ (A.4)

The coefficients $Y^{A(2)}$ for tensor, vector, scalar-breathing and scalar-
longitudinal modes are

$Y^{T(2)}=\frac{4\pi\beta^{2}}{3}\left[\begin{array}[]{c}0\\\
-j_{0}(\alpha)+(2i\eta-1/\alpha)j_{1}(\alpha)+\eta^{2}j_{2}(\alpha)\\\ 0\\\
0\\\
\left(\frac{3}{\alpha^{2}}+\frac{4i\eta}{\alpha}+1\right)j_{2}(\alpha)+\eta^{4}j_{4}(\alpha)+\left(-\frac{6\eta^{2}}{\alpha}+2i\eta\left(\eta^{2}-1\right)\right)j_{3}(\alpha)-\frac{3j_{1}(\alpha)}{\alpha}\\\
\left(\frac{2\left(\eta^{2}-2\right)}{5\alpha^{2}}+\frac{4i\eta}{\alpha}\right)j_{2}(\alpha)+\frac{2\eta^{2}j_{3}(\alpha)}{\alpha}-\frac{2j_{1}(\alpha)}{\alpha}\\\
\left\\{\begin{aligned}
&\left(-\frac{8}{\alpha^{3}}+\frac{2i\eta}{\alpha^{2}}-\frac{\eta^{2}}{\alpha}\right)j_{3}(\alpha)+\left(\frac{11}{\alpha^{2}}+\frac{2i\eta\left(\eta^{2}-1\right)}{\alpha}\right)j_{4}(\alpha)\\\
&+\left(\frac{-\eta^{4}+4\eta^{2}-5}{2\alpha}+\frac{1}{4}i\eta\left(\eta^{2}-1\right)^{2}\right)j_{5}(\alpha)+\frac{1}{8}\left(\eta^{2}+1\right)\left(\eta^{2}-1\right)^{2}j_{6}(\alpha)\end{aligned}\right\\}\\\
\frac{\left(\eta^{2}+1\right)j_{4}(\alpha)}{\alpha^{2}}+\left(\frac{2\left(\eta^{2}-5\right)}{\alpha^{3}}+\frac{2i\eta}{\alpha^{2}}\right)j_{3}(\alpha)\\\
\frac{\left(15\eta^{4}-6\eta^{2}-1\right)j_{4}(\alpha)}{2\alpha^{2}}+\left(\frac{8-24\eta^{2}}{\alpha^{3}}+\frac{2i\eta\left(5\eta^{2}-3\right)}{\alpha^{2}}\right)j_{3}(\alpha)\\\
\left(\frac{2\left(3\eta^{4}-6\eta^{2}+7\right)}{\alpha^{2}}+\frac{2i\eta\left(\eta^{2}-1\right)}{\alpha}\right)j_{4}(\alpha)+\left(\frac{-6\eta^{2}-26}{\alpha^{3}}+\frac{2i\left(3\eta^{3}+\eta\right)}{\alpha^{2}}\right)j_{3}(\alpha)+\frac{\left(\eta^{4}-1\right)j_{5}(\alpha)}{\alpha}\end{array}\right],$

$Y^{V(2)}=\frac{4\pi\beta^{2}}{3}\left[\begin{array}[]{c}0\\\
\eta^{2}j_{2}(\alpha)+\left(2i\eta-\frac{1}{\alpha}\right)j_{1}(\alpha)-j_{0}(\alpha)\\\
0\\\ 0\\\
\left(-\frac{14}{\alpha^{2}}+\frac{8i\eta}{\alpha}\right)j_{2}(\alpha)-\frac{1}{2}\left(\eta^{2}+1\right)^{2}j_{4}(\alpha)+\left(\frac{10\eta^{2}+6}{\alpha}-i\left(\eta^{3}+\eta\right)\right)j_{3}(\alpha)\\\
\left(\frac{-\eta^{2}-13}{\alpha^{2}}+\frac{8i\eta}{\alpha}\right)j_{2}(\alpha)+\left(\frac{6\eta^{2}+10}{\alpha}-2i\eta\right)j_{3}(\alpha)+\left(-\eta^{2}-1\right)j_{4}(\alpha)\\\
\left\\{\begin{aligned}
&\left(-\frac{10}{\alpha^{3}}+\frac{6i\eta}{\alpha^{2}}\right)j_{3}(\alpha)+\left(\frac{9\eta^{2}+7}{\alpha^{2}}+\frac{i\eta\left(\eta^{2}-3\right)}{\alpha}\right)j_{4}(\alpha)\\\
&+\left(\frac{15\eta^{4}-18\eta^{2}-1}{2\alpha}-i\eta^{3}\left(\eta^{2}-1\right)\right)j_{5}(\alpha)+\frac{1}{2}\left(\eta^{2}-\eta^{6}\right)j_{6}(\alpha)\end{aligned}\right\\}\\\
\left(\frac{-\eta^{2}-9}{\alpha^{3}}+\frac{6i\eta}{\alpha^{2}}\right)j_{3}(\alpha)+\left(\frac{4\left(\eta^{2}+3\right)}{\alpha^{2}}-\frac{2i\eta}{\alpha}\right)j_{4}(\alpha)+\frac{\left(-\eta^{2}-1\right)j_{5}(\alpha)}{\alpha}\\\
\left(\frac{10-30\eta^{2}}{\alpha^{3}}+\frac{4i\eta\left(4\eta^{2}-1\right)}{\alpha^{2}}\right)j_{3}(\alpha)+\left(\frac{15\eta^{2}\left(\eta^{2}+2\right)-13}{\alpha^{2}}-\frac{4i\eta\left(2\eta^{2}-1\right)}{\alpha}\right)j_{4}(\alpha)+\frac{\left(1-5\eta^{4}\right)j_{5}(\alpha)}{\alpha}\\\
\left\\{\begin{aligned}
&\left(\frac{-11\eta^{2}-29}{\alpha^{3}}+\frac{4i\eta\left(\eta^{2}+5\right)}{\alpha^{2}}\right)j_{3}(\alpha)+\left(\frac{-3\eta^{4}+4\eta^{2}-17}{2\alpha}-i\eta\left(\eta^{2}-1\right)\right)j_{5}(\alpha)\\\
&+\left(\frac{3\eta^{4}+35\eta^{2}+26}{\alpha^{2}}-\frac{i\eta\left(3\eta^{2}+5\right)}{\alpha}\right)j_{4}(\alpha)+\frac{1}{2}\left(1-\eta^{4}\right)j_{6}(\alpha)\end{aligned}\right\\}\end{array}\right],$

$Y^{B(2)}=\frac{4\pi\beta^{2}}{3}\left[\begin{array}[]{c}\left(\eta^{2}+1\right)j_{2}(\alpha)+\left(2i\eta-\frac{4}{\alpha}\right)j_{1}(\alpha)\\\
\left(\eta^{2}+1\right)j_{2}(\alpha)+\left(2i\eta-\frac{4}{\alpha}\right)j_{1}(\alpha)\\\
\left(-\frac{12}{\alpha^{2}}+\frac{4i\eta}{\alpha}\right)j_{2}(\alpha)+\left(\eta^{4}-1\right)j_{4}(\alpha)+\left(\frac{10-6\eta^{2}}{\alpha}+2i\eta\left(\eta^{2}-1\right)\right)j_{3}(\alpha)\\\
\frac{2\left(\eta^{2}+1\right)j_{3}(\alpha)}{\alpha}+\left(\frac{2\left(\eta^{2}-7\right)}{\alpha^{2}}+\frac{4i\eta}{\alpha}\right)j_{2}(\alpha)\\\
\left(-\frac{12}{\alpha^{2}}+\frac{4i\eta}{\alpha}\right)j_{2}(\alpha)+\left(\eta^{4}-1\right)j_{4}(\alpha)+\left(\frac{10-6\eta^{2}}{\alpha}+2i\eta\left(\eta^{2}-1\right)\right)j_{3}(\alpha)\\\
\frac{2\left(\eta^{2}+1\right)}{\alpha}j_{3}(\alpha)+\left(\frac{2\left(\eta^{2}-7\right)}{\alpha^{2}}+\frac{4i\eta}{\alpha}\right)j_{2}(\alpha)\\\
\left\\{\begin{aligned}
&\left(-\frac{16}{\alpha^{3}}+\frac{4i\eta}{\alpha^{2}}\right)j_{3}(\alpha)+\left(\frac{22-18\eta^{2}}{\alpha^{2}}+\frac{4i\eta\left(\eta^{2}-1\right)}{\alpha}\right)j_{4}(\alpha)\\\
&+\frac{1}{4}\left(\eta^{2}+1\right)\left(\eta^{2}-1\right)^{2}j_{6}(\alpha)+\left(\frac{-\eta^{4}+6\eta^{2}-5}{\alpha}+\frac{1}{2}i\eta\left(\eta^{2}-1\right)^{2}\right)j_{5}(\alpha)\end{aligned}\right\\}\\\
\frac{2\left(\eta^{2}+1\right)}{\alpha^{2}}j_{4}(\alpha)+\left(\frac{4\left(\eta^{2}-5\right)}{\alpha^{3}}+\frac{4i\eta}{\alpha^{2}}\right)j_{3}(\alpha)\\\
\left(\frac{2\left(3\eta^{4}-6\eta^{2}+7\right)}{\alpha^{2}}+\frac{2i\eta\left(\eta^{2}-1\right)}{\alpha}\right)j_{4}(\alpha)+\left(\frac{2i\left(3\eta^{3}+\eta\right)}{\alpha^{2}}-\frac{6\eta^{2}+26}{\alpha^{3}}\right)j_{3}(\alpha)+\frac{\left(\eta^{4}-1\right)}{\alpha}j_{5}(\alpha)\\\
\left(\frac{21\eta^{4}-18\eta^{2}+13}{\alpha^{2}}+\frac{2i\eta\left(\eta^{2}-1\right)}{\alpha}\right)j_{4}(\alpha)+\left(-\frac{2\left(27\eta^{2}+5\right)}{\alpha^{3}}+\frac{2i\eta\left(13\eta^{2}-5\right)}{\alpha^{2}}\right)j_{3}(\alpha)+\frac{\left(\eta^{4}-1\right)}{\alpha}j_{5}(\alpha)\end{array}\right]$

and

$Y^{L(2)}=\frac{4\pi\beta^{2}}{3}\left[\begin{array}[]{c}\frac{1}{2}\left(\eta^{2}+1\right)j_{2}(\alpha)+\left(i\eta-\frac{2}{\alpha}\right)j_{1}(\alpha)\\\
\left(\eta^{2}+1\right)j_{2}(\alpha)+\left(2i\eta-\frac{4}{\alpha}\right)j_{1}(\alpha)\\\
\left(-\frac{8}{\alpha^{2}}+\frac{6i\eta}{\alpha}\right)j_{2}(\alpha)-\left(\eta^{2}+1\right)\eta^{2}j_{4}(\alpha)+\left(\frac{13\eta^{2}+1}{\alpha}-2i\eta^{3}\right)j_{3}(\alpha)\\\
\left(\frac{-2\eta^{2}-6}{\alpha^{2}}+\frac{6i\eta}{\alpha}\right)j_{2}(\alpha)+\left(\frac{5\eta^{2}+9}{\alpha}-2i\eta\right)j_{3}(\alpha)-\left(\eta^{2}+1\right)j_{4}(\alpha)\\\
\left(-\frac{16}{\alpha^{2}}+\frac{12i\eta}{\alpha}\right)j_{2}(\alpha)-2\left(\eta^{2}+1\right)\eta^{2}j_{4}(\alpha)+2\left(\frac{13\eta^{2}+1}{\alpha}-2i\eta^{3}\right)j_{3}(\alpha)\\\
\left(-\frac{4\left(\eta^{2}+3\right)}{\alpha^{2}}+\frac{12i\eta}{\alpha}\right)j_{2}(\alpha)+2\left(\frac{5\eta^{2}+9}{\alpha}-2i\eta\right)j_{3}(\alpha)-2\left(\eta^{2}+1\right)j_{4}(\alpha)\\\
\left\\{\begin{aligned}
&\left(-\frac{18}{\alpha^{3}}+\frac{15i\eta}{\alpha^{2}}\right)j_{3}(\alpha)+\left(\frac{3\left(33\eta^{2}+1\right)}{2\alpha^{2}}-\frac{10i\eta^{3}}{\alpha}\right)j_{4}(\alpha)\\\
&+\frac{1}{2}\left(\eta^{6}+\eta^{4}\right)j_{6}(\alpha)+\left(-\frac{\eta^{2}\left(13\eta^{2}+3\right)}{\alpha}+i\eta^{5}\right)j_{5}(\alpha)\end{aligned}\right\\}\\\
\left\\{\begin{aligned}
&\left(\frac{3\left(13\eta^{2}+21\right)}{2\alpha^{2}}-\frac{10i\eta}{\alpha}\right)j_{4}(\alpha)+\left(-\frac{6\left(\eta^{2}+2\right)}{\alpha^{3}}+\frac{30i\eta}{2\alpha^{2}}\right)j_{3}(\alpha)\\\
&+\left(-\frac{7\eta^{2}+9}{\alpha}+i\eta\right)j_{5}(\alpha)+\frac{1}{2}\left(\eta^{2}+1\right)j_{6}(\alpha)\end{aligned}\right\\}\\\
\left\\{\begin{aligned}
&2\left(-\frac{13\eta^{2}+5}{\alpha^{3}}+\frac{3i\eta\left(2\eta^{2}+3\right)}{\alpha^{2}}\right)j_{3}(\alpha)+\left(-\frac{9\eta^{4}+22\eta^{2}+1}{\alpha}+2i\eta^{3}\right)j_{5}(\alpha)\\\
&+\left(\frac{12\eta^{4}+77\eta^{2}+13}{\alpha^{2}}-\frac{2i\eta\left(7\eta^{2}+3\right)}{\alpha}\right)j_{4}(\alpha)+\left(\eta^{4}+\eta^{2}\right)j_{6}(\alpha)\end{aligned}\right\\}\\\
\left\\{\begin{aligned}
&\left(-\frac{4\left(13\eta^{2}+5\right)}{\alpha^{3}}+\frac{12i\eta\left(2\eta^{2}+3\right)}{\alpha^{2}}\right)j_{3}(\alpha)+2\left(-\frac{9\eta^{4}+22\eta^{2}+1}{\alpha}+2i\eta^{3}\right)j_{5}(\alpha)\\\
&+2\left(\frac{12\eta^{4}+77\eta^{2}+13}{\alpha^{2}}-\frac{2i\eta\left(7\eta^{2}+3\right)}{\alpha}\right)j_{4}(\alpha)+2\left(\eta^{4}+\eta^{2}\right)j_{6}(\alpha)\end{aligned}\right\\}\end{array}\right]$

respectively.

Combining eq. (3.13), eq. (A.1) and eq. (A.2), we get the ORF frequency
second-order expansion,

$\displaystyle\Gamma_{IJ}^{A2}($ $\displaystyle
f)=\frac{e^{\frac{i}{2}(\beta_{J}-\beta_{I})}}{32\pi}\sum_{i,j}P_{I,i}P_{J,j}\bigg{\\{}A^{A(0)}+A^{A(2)}_{I}+A^{A(2)}_{J}+2\left(B^{A(0)}+B^{A(2)}_{I}+B_{J}^{A(2)}\right)\left(\hat{u}_{I,i}\cdot\hat{u}_{J,j}\right)^{2}$
(A.5)
$\displaystyle+\left(C^{A(0)}+C_{I,1}^{A(2)}+C_{J,1}^{A(2)}\right)\left(\left(\hat{u}_{I,i}\cdot\hat{s}_{ij}\right)^{2}+\left(\hat{u}_{J,j}\cdot\hat{s}_{ij}\right)^{2}\right)+\left(C_{I,2}^{A(2)}+C_{J,2}^{A(2)}\right)\left(\left(\hat{u}_{I,i}\cdot\hat{u}_{J,j}\right)^{2}+1\right)$
$\displaystyle+4\left(D^{A(0)}+D_{I,1}^{A(2)}+D_{J,1}^{A(2)}\right)\left(\hat{u}_{I,i}\cdot\hat{u}_{J,j}\right)\left(\hat{u}_{I,i}\cdot\hat{s}_{ij}\right)\left(\hat{u}_{J,j}\cdot\hat{s}_{ij}\right)+4\left(D_{I,2}^{A(2)}+D_{J,2}^{A(2)}\right)\left(\hat{u}_{I,i}\cdot\hat{u}_{J,j}\right)^{2}$
$\displaystyle+\left(E^{A(0)}+E_{I,1}^{A(2)}+E_{J,1}^{A(2)}\right)\left(\hat{u}_{I,i}\cdot\hat{s}_{ij}\right)^{2}\left(\hat{u}_{J,j}\cdot\hat{s}_{ij}\right)^{2}+\left(E_{I,2}^{A(2)}+E_{J,2}^{A(2)}\right)\left(\hat{u}_{I,i}\cdot\hat{u}_{J,j}\right)^{2}$
$\displaystyle+E_{I,3}^{A(2)}\left(\left(\hat{u}^{\prime}_{j}\cdot\hat{s}_{ij}\right)^{2}+\left(\hat{u}_{I,i}\cdot\hat{s}_{ij}\right)^{2}\left(\hat{u}_{I,i}\cdot\hat{u}_{J,j}\right)^{2}\right)+E_{J,3}^{A(2)}\left(\left(\hat{u}_{I,i}\cdot\hat{s}_{ij}\right)^{2}+\left(\hat{u}_{J,j}\cdot\hat{s}_{ij}\right)^{2}\left(\hat{u}_{I,i}\cdot\hat{u}_{J,j}\right)^{2}\right)$
$\displaystyle+4\left(E_{I,4}^{A(2)}+E_{J,4}^{A(2)}\right)\left(\hat{u}_{I,i}\cdot\hat{u}_{J,j}\right)\left(\hat{u}_{I,i}\cdot\hat{s}_{ij}\right)\left(\hat{u}_{J,j}\cdot\hat{s}_{ij}\right)\bigg{\\}}.$

## Acknowledgments

This work is supported by the National Key R&D Program of China under Grant
No. 2022YFC2204602, the Natural Science Foundation of China (Grants No.
12247154, No. 11925503), the Postdoctoral Science Foundation of China (Grant
No. 2022M711259), Guangdong Major project of Basic and Applied Basic Research
(Grant No. 2019B030302001).

## References

  * [1] LIGO Scientific Collaboration and Virgo Collaboration collaboration, _Observation of Gravitational Waves from a Binary Black Hole Merger_ , _Phys. Rev. Lett._ 116 (2016) 061102.
  * [2] LIGO Scientific Collaboration and Virgo Collaboration collaboration, _GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral_ , _Phys. Rev. Lett._ 119 (2017) 161101.
  * [3] R. Abbott, T. Abbott, S. Abraham, F. Acernese, K. Ackley, A. Adams et al., _Observation of gravitational waves from two neutron star–black hole coalescences_ , _Astrophys. J. Lett._ 915 (2021) L5.
  * [4] LIGO Scientific Collaboration and Virgo Collaboration collaboration, _Binary Black Hole Mergers in the First Advanced LIGO Observing Run_ , _Phys. Rev. X_ 6 (2016) 041015.
  * [5] LIGO Scientific Collaboration and Virgo Collaboration collaboration, _GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs_ , _Phys. Rev. X_ 9 (2019) 031040.
  * [6] LIGO Scientific Collaboration and Virgo Collaboration collaboration, _GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo during the First Half of the Third Observing Run_ , _Phys. Rev. X_ 11 (2021) 021053 [arXiv:2010.14527].
  * [7] C. Hogan, _Gravitational radiation from cosmological phase transitions_ , _Mon. Not. R. Astron. Soc._ 218 (1986) 629.
  * [8] C. Caprini, R. Durrer and X. Siemens, _Detection of gravitational waves from the QCD phase transition with pulsar timing arrays_ , _Phys. Rev. D_ 82 (2010) 063511 [arXiv:1007.1218].
  * [9] V. Mukhanov, H. Feldman and R. Brandenberger, _Theory of cosmological perturbations_ , _Phys. Rep._ 215 (1992) 203.
  * [10] M.S. Turner, _Detectability of inflation-produced gravitational waves_ , _Phys. Rev. D_ 55 (1997) R435 [astro-ph/9607066].
  * [11] S. Vagnozzi, _Implications of the NANOGrav results for inflation_ , _Mon. Not. R. Astron. Soc._ 502 (2021) L11 [arXiv:2009.13432].
  * [12] M. Benetti, L.L. Graef and S. Vagnozzi, _Primordial gravitational waves from NANOGrav: A broken power-law approach_ , _Phys. Rev. D_ 105 (2022) 043520 [arXiv:2111.04758].
  * [13] T. Damour and A. Vilenkin, _Gravitational radiation from cosmic (super)strings: Bursts, stochastic background, and observational windows_ , _Phys. Rev. D_ 71 (2005) 063510 [hep-th/0410222].
  * [14] X. Siemens, V. Mandic and J. Creighton, _Gravitational-Wave Stochastic Background from Cosmic Strings_ , _Phys. Rev. Lett._ 98 (2007) 111101 [astro-ph/0610920].
  * [15] N. Mazumder, S. Mitra and S. Dhurandhar, _Astrophysical motivation for directed searches for a stochastic gravitational wave background_ , _Phys. Rev. D_ 89 (2014) 084076 [arXiv:1401.5898].
  * [16] T. Callister, L. Sammut, S. Qiu, I. Mandel and E. Thrane, _Limits of Astrophysics with Gravitational-Wave Backgrounds_ , _Phys. Rev. X_ 6 (2016) 031018.
  * [17] C. Ungarelli and A. Vecchio, _High energy physics and the very early universe with LISA_ , _Phys. Rev. D_ 63 (2001) 064030 [gr-qc/0003021].
  * [18] C. Caprini, M. Hindmarsh, S. Huber, T. Konstandin, J. Kozaczuk, G. Nardini et al., _Science with the space-based interferometer eLISA. II: Gravitational waves from cosmological phase transitions_ , _JCAP_ 2016 (2016) 001 [arXiv:1512.06239].
  * [19] N. Bartolo, C. Caprini, V. Domcke, D.G. Figueroa, J. Garcia-Bellido, M.C. Guzzetti et al., _Science with the space-based interferometer LISA. IV: probing inflation with gravitational waves_ , _JCAP_ 2016 (2016) 026 [arXiv:1610.06481].
  * [20] C. Caprini, M. Chala, G.C. Dorsch, M. Hindmarsh, S.J. Huber, T. Konstandin et al., _Detecting gravitational waves from cosmological phase transitions with LISA: an update_ , _JCAP_ 2020 (2020) 024 [arXiv:1910.13125].
  * [21] P. Auclair, J.J. Blanco-Pillado, D.G. Figueroa, A.C. Jenkins, M. Lewicki, M. Sakellariadou et al., _Probing the gravitational wave background from cosmic strings with LISA_ , _JCAP_ 2020 (2020) 034 [arXiv:1909.00819].
  * [22] A. Maselli, S. Marassi, V. Ferrari, K. Kokkotas and R. Schneider, _Constraining Modified Theories of Gravity with Gravitational-Wave Stochastic Backgrounds_ , _Phys. Rev. Lett._ 117 (2016) 091102 [arXiv:1606.04996].
  * [23] T. Callister, A.S. Biscoveanu, N. Christensen, M. Isi, A. Matas, O. Minazzoli et al., _Polarization-Based Tests of Gravity with the Stochastic Gravitational-Wave Background_ , _Phys. Rev. X_ 7 (2017) 041058.
  * [24] Z. Arzoumanian, P.T. Baker, H. Blumer, B. Bécsy, A. Brazier, P.R. Brook et al., _The NANOGrav 12.5 yr data set: search for an isotropic stochastic gravitational-wave background_ , _Astrophys. J. Lett._ 905 (2020) L34.
  * [25] S. Chen, R. Caballero, Y. Guo, A. Chalumeau, K. Liu, G. Shaifullah et al., _Common-red-signal analysis with 24-yr high-precision timing of the European Pulsar Timing Array: Inferences in the stochastic gravitational-wave background search_ , _Mon. Not. R. Astron. Soc._ 508 (2021) 4970 [arXiv:2110.13184].
  * [26] B. Goncharov, R. Shannon, D. Reardon, G. Hobbs, A. Zic, M. Bailes et al., _On the evidence for a common-spectrum process in the search for the nanohertz gravitational-wave background with the Parkes Pulsar Timing Array_ , _Astrophys. J. Lett._ 917 (2021) L19 [arXiv:2107.12112].
  * [27] J. Antoniadis, Z. Arzoumanian, S. Babak, M. Bailes, A. Bak Nielsen, P. Baker et al., _The International Pulsar Timing Array second data release: Search for an isotropic gravitational wave background_ , _Mon. Not. R. Astron. Soc._ 510 (2022) 4873 [arXiv:2201.03980].
  * [28] LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration collaboration, _Upper limits on the isotropic gravitational-wave background from Advanced LIGO and Advanced Virgo’s third observing run_ , _Phys. Rev. D_ 104 (2021) 022004 [arXiv:2101.12130].
  * [29] P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender et al., _Laser interferometer space antenna_ , arXiv:1702.00786.
  * [30] J. Luo, L.-S. Chen, H.-Z. Duan, Y.-G. Gong, S. Hu, J. Ji et al., _TianQin: a space-borne gravitational wave detector_ , _Class. Quant. Grav._ 33 (2016) 035010 [arXiv:1512.02076].
  * [31] W.-R. Hu and Y.-L. Wu, _The Taiji Program in Space for gravitational wave physics and the nature of gravity_ , _Natl. Sci. Rev._ 4 (2017) 685.
  * [32] S. Kawamura, M. Ando, N. Seto, S. Sato, T. Nakamura, K. Tsubono et al., _The Japanese space gravitational wave antenna: DECIGO_ , _Class. Quant. Grav._ 28 (2011) 094011.
  * [33] Z.-C. Liang, Y.-M. Hu, Y. Jiang, J. Cheng, J.-d. Zhang and J. Mei, _Science with the TianQin Observatory: Preliminary results on stochastic gravitational-wave background_ , _Phys. Rev. D_ 105 (2022) 022001 [arXiv:2107.08643].
  * [34] G. Wang and W.-B. Han, _Alternative LISA-TAIJI networks: Detectability of the isotropic stochastic gravitational wave background_ , _Phys. Rev. D_ 104 (2021) 104015 [arXiv:2108.11151].
  * [35] N. Seto, _Gravitational wave background search by correlating multiple triangular detectors in the mHz band_ , _Phys. Rev. D_ 102 (2020) 123547 [arXiv:2010.06877].
  * [36] C. Brans and R.H. Dicke, _Mach’s Principle and a Relativistic Theory of Gravitation_ , _Phys. Rev._ 124 (1961) 925.
  * [37] D.M. Eardley, D.L. Lee and A.P. Lightman, _Gravitational-Wave Observations as a Tool for Testing Relativistic Gravity_ , _Phys. Rev. D_ 8 (1973) 3308.
  * [38] T. Jacobson and D. Mattingly, _Einstein-aether waves_ , _Phys. Rev. D_ 70 (2004) 024003 [gr-qc/0402005].
  * [39] E. Sagi, _Propagation of gravitational waves in the generalized tensor-vector-scalar theory_ , _Phys. Rev. D_ 81 (2010) 064031 [arXiv:1001.1555].
  * [40] D. Liang, Y. Gong, S. Hou and Y. Liu, _Polarizations of gravitational waves in $f(R)$ gravity_, _Phys. Rev. D_ 95 (2017) 104034 [arXiv:1701.05998].
  * [41] F. Moretti, F. Bombacigno and G. Montani, _Gauge invariant formulation of metric $f(R)$ gravity for gravitational waves_, _Phys. Rev. D_ 100 (2019) 084014 [arXiv:1906.01899].
  * [42] F. Moretti, F. Bombacigno and G. Montani, _Gravitational Landau damping for massive scalar modes_ , _Eur. Phys. J. C_ 80 (2020) 1 [arXiv:2005.08010].
  * [43] F. Moretti, F. Bombacigno and G. Montani, _The Role of Longitudinal Polarizations in Horndeski and Macroscopic Gravity: Introducing Gravitational Plasmas_ , _Universe_ 7 (2021)  [arXiv:2111.11342].
  * [44] C.M. Will, _The confrontation between general relativity and experiment_ , _Living Rev. Rel._ 17 (2014) 1 [arXiv:1403.7377].
  * [45] N. Seto and A. Taruya, _Measuring a Parity-Violation Signature in the Early Universe via Ground-Based Laser Interferometers_ , _Phys. Rev. Lett._ 99 (2007) 121101 [arXiv:0707.0535].
  * [46] N. Seto and A. Taruya, _Polarization analysis of gravitational-wave backgrounds from the correlation signals of ground-based interferometers: Measuring a circular-polarization mode_ , _Phys. Rev. D_ 77 (2008) 103001 [arXiv:0801.4185].
  * [47] R. Kato and J. Soda, _Probing circular polarization in stochastic gravitational wave background with pulsar timing arrays_ , _Phys. Rev. D_ 93 (2016) 062003 [arXiv:1512.09139].
  * [48] Y.-K. Chu, G.-C. Liu and K.-W. Ng, _Observation of a polarized stochastic gravitational-wave background in pulsar-timing-array experiments_ , _Phys. Rev. D_ 104 (2021) 124018 [arXiv:2107.00536].
  * [49] K. Lee, F.A. Jenet and R.H. Price, _Pulsar timing as a probe of non-Einsteinian polarizations of gravitational waves_ , _Astrophys. J._ 685 (2008) 1304.
  * [50] S.J. Chamberlin and X. Siemens, _Stochastic backgrounds in alternative theories of gravity: Overlap reduction functions for pulsar timing arrays_ , _Phys. Rev. D_ 85 (2012) 082001 [arXiv:1111.5661].
  * [51] J.R. Gair, J.D. Romano and S.R. Taylor, _Mapping gravitational-wave backgrounds of arbitrary polarisation using pulsar timing arrays_ , _Phys. Rev. D_ 92 (2015) 102003.
  * [52] N.J. Cornish, L. O’Beirne, S.R. Taylor and N. Yunes, _Constraining Alternative Theories of Gravity Using Pulsar Timing Arrays_ , _Phys. Rev. Lett._ 120 (2018) 181101 [arXiv:1712.07132].
  * [53] Z. Arzoumanian, P.T. Baker, H. Blumer, B. Bécsy, A. Brazier, P.R. Brook et al., _The NANOGrav 12.5-year Data Set: Search for Non-Einsteinian Polarization Modes in the Gravitational-wave Background_ , _Astrophys. J. Lett._ 923 (2021) L22 [arXiv:2109.14706].
  * [54] Y.-M. Wu, Z.-C. Chen and Q.-G. Huang, _Constraining the Polarization of Gravitational Waves with the Parkes Pulsar Timing Array Second Data Release_ , _Astrophys. J._ 925 (2022) 37 [arXiv:2108.10518].
  * [55] LIGO Scientific Collaboration and Virgo Collaboration collaboration, _Search for Tensor, Vector, and Scalar Polarizations in the Stochastic Gravitational-Wave Background_ , _Phys. Rev. Lett._ 120 (2018) 201102 [arXiv:1802.10194].
  * [56] LIGO Scientific and Virgo Collaboration collaboration, _Search for the isotropic stochastic background using data from Advanced LIGO’s second observing run_ , _Phys. Rev. D_ 100 (2019) 061101 [arXiv:1903.02886].
  * [57] LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration collaboration, _Upper limits on the isotropic gravitational-wave background from Advanced LIGO and Advanced Virgo’s third observing run_ , _Phys. Rev. D_ 104 (2021) 022004 [arXiv:2101.12130].
  * [58] A. Nishizawa, A. Taruya, K. Hayama, S. Kawamura and M.-a. Sakagami, _Probing nontensorial polarizations of stochastic gravitational-wave backgrounds with ground-based laser interferometers_ , _Phys. Rev. D_ 79 (2009) 082002 [arXiv:0903.0528].
  * [59] D. Liang, Y. Gong, S. Hou and Y. Liu, _Polarizations of gravitational waves in $f(R)$ gravity_, _Phys. Rev. D_ 95 (2017) 104034.
  * [60] D. Liang, Y. Gong, A.J. Weinstein, C. Zhang and C. Zhang, _Frequency response of space-based interferometric gravitational-wave detectors_ , _Phys. Rev. D_ 99 (2019) 104027.
  * [61] C. Liu, W.-H. Ruan and Z.-K. Guo, _Constraining gravitational-wave polarizations with Taiji_ , _Phys. Rev. D_ 102 (2020) 124050.
  * [62] M. Tinto, J.W. Armstrong and F.B. Estabrook, _Discriminating a gravitational wave background from instrumental noise in the LISA detector_ , _Phys. Rev. D_ 63 (2000) 021101.
  * [63] M. Tinto, F.B. Estabrook and J.W. Armstrong, _Time-delay interferometry for LISA_ , _Phys. Rev. D_ 65 (2002) 082003.
  * [64] T.A. Prince, M. Tinto, S.L. Larson and J.W. Armstrong, _LISA optimal sensitivity_ , _Phys. Rev. D_ 66 (2002) 122002 [gr-qc/0209039].
  * [65] M. Tinto and S.V. Dhurandhar, _Time-delay interferometry_ , _Living Rev. Rel._ 24 (2020) 1.
  * [66] C.J. Hogan and P.L. Bender, _Estimating stochastic gravitational wave backgrounds with the Sagnac calibration_ , _Phys. Rev. D_ 64 (2001) 062002 [astro-ph/0104266].
  * [67] M.R. Adams and N.J. Cornish, _Discriminating between a stochastic gravitational wave background and instrument noise_ , _Phys. Rev. D_ 82 (2010) 022002 [arXiv:1002.1291].
  * [68] J.D. Romano, N. Cornish et al., _Detection methods for stochastic gravitational-wave backgrounds: a unified treatment_ , _Living Rev. Rel._ 20 (2017) 1 [arXiv:1608.06889].
  * [69] S.H.S. Alexander, M.E. Peskin and M.M. Sheikh-Jabbari, _Leptogenesis from Gravity Waves in Models of Inflation_ , _Phys. Rev. Lett._ 96 (2006) 081301 [hep-th/0403069].
  * [70] L. Sorbo, _Parity violation in the Cosmic Microwave Background from a pseudoscalar inflaton_ , _JCAP_ 2011 (2011) 003 [arXiv:1101.1525].
  * [71] N. Seto, _Prospects for Direct Detection of the Circular Polarization of the Gravitational-Wave Background_ , _Phys. Rev. Lett._ 97 (2006) 151101 [astro-ph/0609504].
  * [72] V. Domcke, J. García-Bellido, M. Peloso, M. Pieroni, A. Ricciardone, L. Sorbo et al., _Measuring the net circular polarization of the stochastic gravitational wave background with interferometers_ , _JCAP_ 2020 (2020) 028 [arXiv:1910.08052].
  * [73] G. Orlando, M. Pieroni and A. Ricciardone, _Measuring parity violation in the Stochastic Gravitational Wave Background with the LISA-Taiji network_ , _JCAP_ 2021 (2021) 069 [arXiv:2011.07059].
  * [74] B. Allen and J.D. Romano, _Detecting a stochastic background of gravitational radiation: Signal processing strategies and sensitivities_ , _Phys. Rev. D_ 59 (1999) 102001 [gr-qc/9710117].
  * [75] Planck Collaboration, Aghanim, N., Akrami, Y., Ashdown, M., Aumont, J., Baccigalupi, C. et al., _Planck 2018 results - VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [arXiv:1807.06209].
  * [76] R.A. Isaacson, _Gravitational Radiation in the Limit of High Frequency. II. Nonlinear Terms and the Effective Stress Tensor_ , _Phys. Rev._ 166 (1968) 1272.
  * [77] N.J. Cornish and S.L. Larson, _Space missions to detect the cosmic gravitational-wave background_ , _Class. Quant. Grav._ 18 (2001) 3473 [gr-qc/0103075].
  * [78] X.-C. Hu, X.-H. Li, Y. Wang, W.-F. Feng, M.-Y. Zhou, Y.-M. Hu et al., _Fundamentals of the orbit and response for TianQin_ , _Class. Quant. Grav._ 35 (2018) 095008 [arXiv:1803.03368].
  * [79] S.V. Dhurandhar, K.R. Nayak, S. Koshti and J.-Y. Vinet, _Fundamentals of the LISA stable flight formation_ , _Class. Quant. Grav._ 22 (2005) 481 [gr-qc/0410093].
  * [80] F.B. Estabrook and H.D. Wahlquist, _Response of Doppler spacecraft tracking to gravitational radiation_ , _Gen. Rel. Grav._ 6 (1975) 439.
  * [81] E.E. Flanagan, _Sensitivity of the Laser Interferometer Gravitational Wave Observatory to a stochastic background, and its dependence on the detector orientations_ , _Phys. Rev. D_ 48 (1993) 2389.
  * [82] N.J. Cornish, _Detecting a stochastic gravitational wave background with the Laser Interferometer Space Antenna_ , _Phys. Rev. D_ 65 (2001) 022004 [gr-qc/0106058].
  * [83] C. Cutler and M. Vallisneri, _LISA detections of massive black hole inspirals: Parameter extraction errors due to inaccurate template waveforms_ , _Phys. Rev. D_ 76 (2007) 104018 [arXiv:0707.2982].