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

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].

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}$.

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.

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.

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)

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.