Full analytic expression of overlap reduction function for gravitational wave background with pulsar timing arrays

The metric perturbations corresponding to GWB can be expressed as

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

where $e^{A}_{ab}(\hat{n})$ is the spin-2 polarization tensors. $A=\{+,\times,X,Y,B,L\}$ represent different polarization mode, where ${+,\times}$ represent tensor mode predicted by general relativity, ${X,Y}$ and ${B,L}$ represent vector mode and scalar allowed by the general metric theory of gravity. Explicitly,

	$\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)
	$\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},$

where

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

For the isotropic, unpolarized and stationary GWB, the quadratic expected values satisfy that

$$\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),$$

(4)

where $S^{A}_{h}(f)$ can be regarded as the component corresponding to the $A$ polarization of one-sided gravitational-wave strain power spectral density function. For tensor mode and vector mode, it can be considered that $S^{+}_{h}=S^{\times}_{h}=S^{T}_{h}/2$, $S^{X}_{h}=S^{Y}_{h}=S^{V}_{h}/2$. However, the two scalar modes should be considered as two independent polarization modes. The strain power spectral density is simply related to the fractional energy density,

$$S_{h}(f)=\frac{3H_{0}^{2}}{2\pi^{2}}\frac{\Omega_{\mbox{gw}}(f)}{f^{3}}.$$

(5)

The fractional energy density is defined as

$$\Omega_{\mbox{gw}}(f)=\frac{1}{\rho_{c}}\frac{d\rho_{\mbox{gw}}}{d\ln f}$$

(6)

where $\rho_{c}\equiv 3c^{2}H_{0}^{2}/8\pi G$ is the critical energy density required for closed universe, and $H_{0}$, $G$ are the Hubble constant and gravitational constant respectively.

The metric perturbation is weak, so the single pulsar timing signal can be written as

$$h(t)=\int_{-\infty}^{\infty}df\int d^{2}\Omega_{\hat{n}}R^{ab}(f,\hat{n})h_{ab}(f,\hat{n})e^{i2\pi ft}.$$

(7)

Here the response function for Doppler frequency measurements is 2017Review_detection_GWB

$$R^{ab}(f,\hat{n})=\frac{1}{2}\frac{u^{a}u^{b}}{1+\hat{n}\cdot\hat{u}}\left(1-e^{-\frac{i2\pi fL}{c}(1+\hat{n}\cdot\hat{u})}\right)e^{i2\pi f\hat{n}\cdot\vec{r}_{2}/c},$$

(8)

where $\hat{u}$ is the unit vector of the pulsar pointing to the earth, $L$ is the distance from the pulsar to the earth, and $\vec{r}_{2}$ is the position vector of the earth. In the frequency domain, the signal is written in terms of polarized basis as

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

(9)

where $R^{A}(f,\hat{n})=R^{ab}(f,\hat{n})e^{A}_{ab}(\hat{n})$.

The cross-correlation of two pulsars signal $I$ and $J$ is

$$\left\langle\tilde{h}_{I}(f)\tilde{h}_{J}^{*}(f^{\prime})\right\rangle=\frac{1}{2}\delta(f-f^{\prime})\sum_{A}\Gamma^{A}_{IJ}(f)S^{A}_{h}(f),$$

(10)

where $A=\{T,V,B,L\}$. The ORF for tensor mode, vector mode, scalar-breathing mode and scalar-longitudinal is

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

The ORF of the tensor mode is

	$\displaystyle\Gamma$	$\displaystyle{}^{T}_{IJ}(f)=\Gamma^{T}_{IJ}(\beta_{u},\beta_{v},\gamma)$		(13)
		$\displaystyle=\frac{1}{32\pi}\int d^{2}\Omega_{\hat{n}}\frac{(1-e^{-i\beta_{u}(1+\sin\theta\cos\underset{\sim}{\phi})})}{(1+\sin\theta\cos\underset{\sim}{\phi})}$
		$\displaystyle\times\frac{(1-e^{i\beta_{v}(1+\sin\theta\cos\tilde{\phi})})}{(1+\sin\theta\cos\tilde{\phi})}$
		$\displaystyle\times((1-\sin^{2}\theta\cos^{2}\underset{\sim}{\phi})(1-\sin^{2}\theta\cos^{2}\tilde{\phi})-2\cos^{2}\theta\sin^{2}\gamma),$

where $\beta_{u}=2\pi fL_{u}/c$, $\beta_{v}=2\pi fL_{v}/c$, and $L_{u},L_{v}$ are the distances from the two pulsars to the Earth, respectively. Change the integral variables $\{\theta,\phi\}$ by $\{x,y\}$, so that the integral region is transformed from the unit sphere to the elliptical area on the $xy$ plane. Then,

	$\displaystyle\Gamma$	$\displaystyle{}^{T}_{IJ}(f)=\frac{1}{16\pi}\int_{-1}^{1}dx\int_{y_{l}}^{y_{u}}dy(1-e^{-i\beta_{u}(1+x)})$		(14)
		$\displaystyle\times(1-e^{i\beta_{v}(1+y)})\left(\frac{(1-x)(1-y)}{\eta(x,y)}-\frac{2\eta(x,y)}{(1+x)(1+y)}\right),$

where $y_{l},y_{u}=x\cos\gamma\mp\sin\gamma\sqrt{1-x^{2}}$ and

$$\eta(x,y)=\sqrt{\sin^{2}\gamma-(x^{2}+y^{2}-2xy\cos\gamma)}$$

is the boundary of the integration area. The integral can be done analytically (the details of the calculation are provided in Appendix), and the result is

	$\displaystyle\Gamma^{T}_{IJ}(f)$	$\displaystyle=\frac{1}{3}\gamma^{T}_{IJ}(f)=\frac{1}{16}\bigg{\{}2+\frac{2\cos\gamma}{3}+16\sin^{2}\frac{\gamma}{2}\mbox{ln}(\sin\frac{\gamma}{2})-2i(1-\cos\gamma)\left(\frac{e^{-2i\beta_{u}}}{\beta_{u}}-\frac{e^{2i\beta_{v}}}{\beta_{v}}\right)$		(15)
		$\displaystyle-\frac{(3+e^{-2i\beta_{u}})\cos\gamma+1-e^{-2i\beta_{u}}}{\beta_{u}^{2}}-\frac{(3+e^{2i\beta_{v}})\cos\gamma+1-e^{2i\beta_{v}}}{\beta_{v}^{2}}-2i\cos\gamma\left(\frac{1-e^{-2i\beta_{u}}}{\beta_{u}^{3}}-\frac{1-e^{2i\beta_{v}}}{\beta_{v}^{3}}\right)$
		$\displaystyle+8\sin^{2}\frac{\gamma}{2}\left(\mbox{Ei}(-2i\beta_{u})+\mbox{Ei}(2i\beta_{v})-\mbox{Ei}(-i(\beta_{s}+\beta_{u}-\beta_{v}))-\mbox{Ei}(i(\beta_{s}-\beta_{u}+\beta_{v}))\right)$
		$\displaystyle+\frac{e^{-i(\beta_{u}-\beta_{v})}}{\beta_{u}\beta_{v}}\left[\cos\beta_{s}\left(-\frac{1}{2}-3i(\beta_{u}-\beta_{v})-\frac{\beta^{2}}{\beta_{s}^{2}}+\frac{3(\beta_{u}^{2}-\beta_{v}^{2})^{2}}{2\beta_{s}^{4}}-\frac{i(\beta_{u}-\beta_{v})(\beta_{u}+\beta_{v})^{2}}{\beta_{s}^{2}}\right)\right.$
		$\displaystyle+\sin\beta_{s}\left(-\frac{3(\beta_{u}^{2}-\beta_{v}^{2})^{2}}{2\beta_{s}^{5}}+\frac{2\beta^{2}+(\beta_{u}^{2}-\beta_{v}^{2})^{2}}{2\beta_{s}^{3}}\right.\left.\left.+\frac{1+4\beta_{u}\beta_{v}-2i(\beta_{u}-\beta_{v})}{2\beta_{s}}+\frac{i(\beta_{u}-\beta_{v})(\beta_{u}+\beta_{v})^{2}}{\beta_{s}^{3}}+\frac{7\beta_{s}}{2}\right)\right]\bigg{\}},$

where

$$\beta=\sqrt{\beta_{u}^{2}+\beta_{v}^{2}},$$

$$\beta_{s}=\sqrt{\beta_{u}^{2}+\beta_{v}^{2}-2\beta_{u}\beta_{v}\cos\gamma},$$

and $\text{Ei}(z)=-\int_{-z}^{\infty}e^{-t}/tdt$ is the exponential integral function. $\gamma^{T}_{IJ}(f)$ is the normalized ORF, so that $\gamma^{T}_{IJ}(f)=1$ for a single pulsar. The normalization holds only for the short-wave approximation. Using the short-wave approximation $\beta_{u},\beta_{v}\gg 1$, Eq. (15) will reduce to the HD function 1983apj_HD_curve ,

$$\gamma_{HD}(\gamma)=\frac{3}{2}\left(\frac{1-\cos\gamma}{2}\right)\ln\frac{1-\cos\gamma}{2}-\frac{1-\cos\gamma}{8}+\frac{1}{2}.$$

(16)

In addition, it can also accurately converge to 1 when the two pulsars tend to coincide, because the exponential terms is included.

Since no approximations are used, the analytic expressions agree exactly with the values obtained by direct numerical integration. When the distances of two pulsars are the same ($\beta_{u}=\beta_{v}$), the ORF is a real function. As shown in Fig. 2, our result is consistent with the HD curve as $\beta_{u}$ increase, and the curve approaches 1 smoothly when the two pulsars are very close. For the long-wave wavelengths, $\lim_{\gamma\rightarrow 0}\gamma^{T}_{IJ}$ is a value less than 0.5, which depends on $\beta_{u}$. However, the two pulsars are usually at different distances in reality, which implies that the ORF is a complex function and satisfies that $\Gamma_{IJ}(-f)=\Gamma^{*}_{IJ}(f)$. As shown in Fig. 3, for different distances of pulsars ($\beta_{u}\neq\beta_{v}$), the imaginary part of the ORF depends on the ratio of $\beta_{u}$ and $\beta_{v}$. The greater the difference between them, the greater the value of the imaginary part. But when both $\beta_{u}$ and $\beta_{v}$ are very large, the imaginary part is close to 0, no matter how large their ratio is. And the real part is tend to HD function in this case.

We show the behavior in Fig. LABEL:fig3 when the angle between the two pulsars is small. If $\beta_{u}=\beta_{v}$, $\gamma^{T}_{IJ}$ decays from 1 to 0.5 faster and faster as $\beta_{u}$ increases. From their similar curves, we can conclude that the ORF for two nearby pulsars depends on their distance from each other. For pulsars in the same direction $\gamma=0$,

	$\displaystyle\lim_{\gamma\rightarrow 0}$	$\displaystyle\gamma^{T}_{IJ}(f)=\frac{1}{2}-\frac{3}{4\beta_{u}^{2}}-\frac{3}{4\beta_{v}^{2}}+\frac{3}{4(\beta_{u}-\beta_{v})^{2}}$		(17)
		$\displaystyle-i\frac{3}{8}\left(\frac{1-e^{-2i\beta_{u}}}{\beta_{u}^{3}}-\frac{1-e^{2i\beta_{v}}}{\beta_{v}^{3}}\right)$
		$\displaystyle-e^{-i(\beta_{u}-\beta_{v})}\frac{3\sin(\beta_{u}-\beta_{v})}{4(\beta_{u}-\beta_{v})^{3}}-\frac{3i(\beta_{u}^{2}-\beta_{u}\beta_{v}+\beta_{v}^{2})}{4\beta_{u}\beta_{v}(\beta_{u}-\beta_{v})}.$

In addition, we noticed that $\lim_{\gamma\rightarrow 0}\gamma^{T}_{IJ}\neq 1$ when $\beta_{u}$ is not very large. Actually, the auto-correlated response for the tensor mode of a single detector is

$$\gamma^{T}_{II}(f)=\lim_{\begin{subarray}{c}\gamma\rightarrow 0\\ \beta_{v}\rightarrow\beta_{u}\end{subarray}}\gamma^{T}_{IJ}(f)=1-\frac{3}{2\beta_{u}^{2}}+\frac{3\sin 2\beta_{u}}{4\beta_{u}^{3}}.$$

(18)

It dose converge to 1 for $\beta_{u}\gg 1$, but gradually decreases to 0 as $\beta_{u}$ decreases. The normalization expression $\gamma^{T}_{IJ}(f)$ is only valid for short-wave approximation.

The ORF of the vector mode isotropic is

$\displaystyle\Gamma_{IJ}^{V}(f)=\frac{1}{32\pi}\int d^{2}\Omega_{\hat{n}}\frac{(1-e^{-i\beta_{u}(1+\sin\theta\cos\underset{\sim}{\phi})})}{(1+\sin\theta\cos\underset{\sim}{\phi})}\frac{(1-e^{i\beta_{v}(1+\sin\theta\cos\tilde{\phi})})}{(1+\sin\theta\cos\tilde{\phi})}4\left(\sin^{2}\theta\cos\underset{\sim}{\phi}\cos\tilde{\phi}(\cos\gamma-\sin^{2}\theta\cos\underset{\sim}{\phi}\cos\tilde{\phi})\right).$

(19)

Changing the variables, then we get

$\displaystyle\Gamma$

$\displaystyle{}^{V}_{IJ}(f)=\frac{1}{16\pi}\int_{-1}^{1}dx\int_{y_{l}}^{y_{u}}dy(1-e^{-i\beta_{u}(1+x)})(1-e^{i\beta_{v}(1+y)})\frac{4xy(\cos\gamma-xy)}{\eta(x,y)(1+x)(1+y)},$

(20)

Finishing the integral, and we find that

	$\displaystyle\Gamma_{IJ}^{V}(f)$	$\displaystyle=-\frac{1}{2}-\frac{2\cos\gamma}{3}-\mbox{ln}\sin\frac{\gamma}{2}+\frac{i\left(2\cos\gamma+1-e^{2i\beta_{v}}\right)}{4\beta_{v}}-\frac{i\left(2\cos\gamma+1-e^{-2i\beta_{u}}\right)}{4\beta_{u}}$		(21)
		$\displaystyle+\frac{\left(3+e^{-2i\beta_{u}}\right)\cos\gamma-e^{-2i\beta_{u}}+1}{4\beta_{u}^{2}}+\frac{\left(3+e^{2i{\beta_{v}}}\right)\cos\gamma-e^{2i\beta_{v}}+1}{4\beta_{v}^{2}}+\frac{i\cos\gamma}{2}\left(\frac{1-e^{-2i\beta_{u}}}{\beta_{u}^{3}}-\frac{1-e^{2i\beta_{v}}}{\beta_{v}^{3}}\right)$
		$\displaystyle+\frac{1}{2}e^{-i(\beta_{u}-\beta_{v})}\left\{\cos\beta_{s}\left(\frac{3\beta_{u}\beta_{v}\sin^{2}\gamma}{\beta_{s}^{4}}+\frac{-2\cos\gamma+i(\beta_{u}-\beta_{v})(\cos\gamma+1)}{\beta_{s}^{2}}-i\left(\frac{1}{\beta_{u}}-\frac{1}{\beta_{v}}\right)\right)\right.$
		$\displaystyle\left.+\sin\beta_{s}\left(-\frac{3\beta_{u}\beta_{v}\sin^{2}\gamma}{\beta_{s}^{5}}+\frac{\beta_{u}\beta_{v}\sin^{2}\gamma-i(\beta_{u}-\beta_{v})(\cos\gamma+1)+2\cos\gamma}{\beta_{s}^{3}}-\frac{2\cos\gamma+1}{\beta_{s}}-\frac{\beta_{s}}{\beta_{u}\beta_{v}}\right)\right\}$
		$\displaystyle+\frac{1}{2}\left(\mbox{Ei}(-i(\beta_{s}+\beta_{u}-{\beta_{v}}))+\mbox{Ei}(i(\beta_{s}-\beta_{u}+\beta_{v}))-\mbox{Ei}(-2i\beta_{u})-\mbox{Ei}(2i\beta_{v})\right).$

The ORF for vector mode cannot be normalized, because $\Gamma_{IJ}^{V}(\gamma=0)$ increase with $\beta_{u}$. The expression with using short-wave approximation is 2008apj_Lee_ORF_nonGR

$$\Gamma^{V}_{sw}(\gamma)=-\frac{1}{2}-\frac{2}{3}\cos\gamma-\ln(\sin\frac{\gamma}{2}),$$

(22)

where the angular separation $\gamma$ is assumed to be not small. In the limit $\gamma\rightarrow 0$, Eq. (21) will give the exact value while the short-wave expression $\Gamma^{V}_{sw}$ diverges. The auto-correlated response for vector mode of a single detector is

	$\displaystyle\Gamma^{V}_{II}(f)=$	$\displaystyle\lim_{\begin{subarray}{c}\gamma\rightarrow 0\\ \beta_{v}\rightarrow\beta_{u}\end{subarray}}\Gamma^{V}_{IJ}(f)=-\frac{7}{3}+\gamma_{E}+\text{ln}(2\beta_{u})-\text{Ci}(2\beta_{u})$		(23)
		$\displaystyle+\frac{2}{\beta^{2}_{u}}+\frac{(\beta^{2}_{u}-2)\cos\beta_{u}\sin\beta_{u}}{\beta^{3}_{u}},$

where $\gamma_{E}$ is the Euler constant and

$$\text{Ci}(z)=-\int_{z}^{\infty}\cos t/tdt$$

is the cosine-integral function. In the limit $\gamma\rightarrow 0$, $\Gamma^{V}$ grows logarithmically with $2\beta_{u}$. It is consistent with Eq.(A42) in 2008apj_Lee_ORF_nonGR .

The behaviors of ORF for vector mode are shown in Fig. 2 and Fig. 3 in the case $\beta_{u}=\beta_{v}$ and $\beta_{u}\neq\beta_{v}$ respectively. For two pulsars at the same distance, the ORF for vector mode gradually approaches the short-wave expression as $\beta_{u}$ increases. From the short-wave expression (22), we know that $\Gamma^{V}_{sw}(\gamma)$ diverges at $\gamma=0$. However, at $\gamma=0$, the full analytic expression tends to have finite values $\Gamma^{V}_{II}(\beta_{u})$, which grows logarithmically with $2\beta_{u}$. For two pulsars at different distances, the ORF for vector mode is complex. The real part will also gradually approaches the short-wave expression for not too small of $\gamma$. However, when $\gamma$ is very small, the imaginary part is not negligible compared to the real part even $\beta_{u}$ and $\beta_{v}$ are pretty large. The small angle behavior of vector mode is quite different from the tensor mode. As shown in Fig. LABEL:fig3, the ORF for vector mode converges to a finite value, which depends on $\beta_{u}$ and increases with $\beta_{u}$ logarithmically. There is also no fast halving decay, and the transition of the curve is smooth.

The ORF of scalar-breathing mode is

	$\displaystyle\Gamma_{IJ}^{B}(f)=$	$\displaystyle\frac{1}{16\pi}\int d^{2}\Omega_{\hat{n}}\frac{(1-e^{-i\beta_{u}(1+\sin\theta\cos\underset{\sim}{\phi})})}{(1+\sin\theta\cos\underset{\sim}{\phi})}\frac{(1-e^{i\beta_{v}(1+\sin\theta\cos\tilde{\phi})})}{(1+\sin\theta\cos\tilde{\phi})}(1-\sin^{2}\theta\cos^{2}\underset{\sim}{\phi})(1-\sin^{2}\theta\cos^{2}\tilde{\phi})$		(24)
	$\displaystyle=$	$\displaystyle\frac{1}{8\pi}\int_{-1}^{1}dx\int_{y_{l}}^{y_{u}}dy(1-e^{-i\beta_{u}(1+x)})(1-e^{i\beta_{v}(1+y)})\frac{(1-x)(1-y)}{\eta(x,y)}.$

The normalized analytic expression of ORF for scalar-breathing mode is

	$\displaystyle\gamma_{IJ}^{B}(f)$	$\displaystyle=\frac{3}{2}\Gamma_{IJ}^{B}(f)=\frac{3}{16}\bigg{\{}2+\frac{2\cos\gamma}{3}+2i(1+\cos\gamma)\left(\frac{1}{\beta_{u}}-\frac{1}{\beta_{v}}\right)-2i\cos\gamma\left(\frac{1-e^{-2i\beta_{u}}}{\beta_{u}^{3}}-\frac{1-e^{2i\beta_{v}}}{\beta_{v}^{3}}\right)$		(25)
		$\displaystyle+\frac{e^{-2i\beta_{u}}-1-\left(3+e^{-2i\beta_{u}}\right)\cos\gamma}{\beta_{u}^{2}}+\frac{e^{2i\beta_{v}}-1-\left(3+e^{2i\beta_{v}}\right)\cos\gamma}{\beta_{v}^{2}}$
		$\displaystyle+\frac{e^{-i(\beta_{u}-\beta_{v})}}{\beta_{u}\beta_{v}}\left[\sin\beta_{s}\left(\frac{2\beta^{2}+\left(\beta_{u}^{2}-\beta_{v}^{2}\right)^{2}}{2\beta_{s}^{3}}-\frac{3\left(\beta_{u}^{2}-\beta_{v}^{2}\right)^{2}}{2\beta_{s}^{5}}\right.\right.$
		$\displaystyle\left.+\frac{i(\beta_{u}-\beta_{v})(\beta_{u}+\beta_{v})^{2}}{\beta_{s}^{3}}+\frac{-2i({\beta_{u}}-\beta_{v})+4\beta_{u}\beta_{v}+1}{2\beta_{s}}-\frac{\beta_{s}}{2}\right)$
		$\displaystyle\left.+\cos\beta_{s}\left(-\frac{\beta^{2}}{\beta_{s}^{2}}+\frac{3\left(\beta_{u}^{2}-\beta_{v}^{2}\right)^{2}}{2\beta_{s}^{4}}-\frac{i(\beta_{u}-\beta_{v})(\beta_{u}+\beta_{v})^{2}}{\beta_{s}^{2}}+i(\beta_{u}-\beta_{v})-\frac{1}{2}\right)\right]\bigg{\}}.$

The short-wave expression is 2008apj_Lee_ORF_nonGR

$$\gamma^{B}_{sw}(\gamma)=\frac{3}{8}+\frac{1}{8}\cos\gamma.$$

(26)

The auto-correlated response for the scalar-breathing mode of a single detector is the same as that of the tensor mode,

$$\gamma^{B}_{II}(f)=\lim_{\begin{subarray}{c}\gamma\rightarrow 0\\ \beta_{v}\rightarrow\beta_{u}\end{subarray}}\gamma^{B}_{IJ}(f)=1-\frac{3}{2\beta_{u}^{2}}+\frac{3\sin 2\beta_{u}}{4\beta_{u}^{3}}.$$

(27)

But the actual responses satisfy that $\Gamma^{B}_{II}(f)=1/2\Gamma^{T}_{II}(f)$, because their normalization factors are different. Overall, the ORF for the scalar-breathing mode is very similar to that for the tensor mode, probably because they are both transverse modes. In particular, it should be noticed that $\Gamma_{IJ}^{B}(f)\geq 0$ for any configuration of pulsar pairs.

The ORF of scalar-longitudinal mode is

	$\displaystyle\Gamma_{IJ}^{L}(f)=$	$\displaystyle\frac{1}{16\pi}\int d^{2}\Omega_{\hat{n}}\frac{(1-e^{-i\beta_{u}(1+\sin\theta\cos\underset{\sim}{\phi})})}{(1+\sin\theta\cos\underset{\sim}{\phi})}\frac{(1-e^{i\beta_{v}(1+\sin\theta\cos\tilde{\phi})})}{(1+\sin\theta\cos\tilde{\phi})}2(\sin^{4}\theta\cos^{2}\underset{\sim}{\phi}\cos^{2}\tilde{\phi})$		(28)
	$\displaystyle=$	$\displaystyle\frac{1}{8\pi}\int_{-1}^{1}dx\int_{y_{l}}^{y_{u}}dy(1-e^{-i\beta_{u}(1+x)})(1-e^{i\beta_{v}(1+y)})\frac{2x^{2}y^{2}}{\eta(x,y)(1+x)(1+y)}.$

The analytic expression of ORF for scalar-longitudinal mode is

	$\displaystyle\Gamma_{IJ}^{L}(f)$	$\displaystyle=\frac{1}{8}\biggl{\{}4+\frac{28\cos\gamma}{3}+\csc^{2}\frac{\gamma}{2}\left(4\ln(\sin\frac{\gamma}{2})+2\gamma_{E}\cos^{2}\gamma+\cos^{2}\gamma\ln(4\beta_{u}\beta_{v})\right)$		(29)
		$\displaystyle+2i\left(\frac{1}{\beta_{u}}-\frac{1}{\beta_{v}}\right)\left(\frac{1}{2}\cos\beta_{s}e^{-i(\beta_{u}-\beta_{v})}\left(4-\csc^{2}\frac{\gamma}{2}-\sec^{2}\frac{\gamma}{2}\right)+3\cos\gamma+1\right)$
		$\displaystyle-2i\left(\frac{e^{-2i\beta_{u}}}{\beta_{u}}-\frac{e^{2i\beta_{v}}}{\beta_{v}}\right)(\cos\gamma+1)-4i\cos\gamma\left(\frac{1-e^{-2i\beta_{u}}}{\beta_{u}^{3}}-\frac{1-e^{2i\beta_{v}}}{\beta_{v}^{3}}\right)$
		$\displaystyle+\frac{-2\left(3+e^{-2i\beta_{u}}\right)\cos\gamma+2e^{-2i\beta_{u}}-2}{\beta_{u}^{2}}+\frac{-2\left(3+e^{2i\beta_{v}}\right)\cos\gamma+2e^{2i\beta_{v}}-2}{\beta_{v}^{2}}$
		$\displaystyle-\frac{1}{2}(\cos 2\gamma-3)\csc^{2}\frac{\gamma}{2}\left(\mbox{Ei}(-2i\beta_{u})+\mbox{Ei}(2i\beta_{v})\right)-2\csc^{2}\frac{\gamma}{2}\left(\mbox{Ei}(-i(\beta_{s}+\beta_{u}-\beta_{v}))+\mbox{Ei}(i(\beta_{s}-\beta_{u}+\beta_{v}))\right)$
		$\displaystyle-\cos^{2}\gamma\csc^{2}\frac{\gamma}{2}\left[e^{-2i\beta_{u}\sin^{2}\frac{\gamma}{2}}\left(\mbox{Ei}(-i\beta_{u}(\cos\gamma-1))+\mbox{Ei}(-i\beta_{u}(\cos\gamma+1))\right.\right.$
		$\displaystyle\left.-\mbox{Ei}(i(\beta_{s}+\beta_{v}-\beta_{u}\cos\gamma))-\mbox{Ei}(-i(\beta_{s}-\beta_{v}+\beta_{u}\cos\gamma))\right)$
		$\displaystyle\left.+e^{2i\beta_{v}\sin^{2}\frac{\gamma}{2}}\left(\mbox{Ei}(i\beta_{v}(\cos\gamma-1))+\mbox{Ei}(i\beta_{v}(\cos\gamma+1))\right.\left.-\mbox{Ei}(-i(\beta_{s}+\beta_{u}-\beta_{v}\cos\gamma))-\mbox{Ei}(i(\beta_{s}-\beta_{u}+\beta_{v}\cos\gamma))\right)\right]$
		$\displaystyle+e^{-i(\beta_{u}-\beta_{v})}\left[\cos\beta_{s}\left(-\frac{12\beta_{u}\beta_{v}\sin^{2}\gamma}{\beta_{s}^{4}}-\frac{-8\cos\gamma+4i(\beta_{u}-\beta_{v})(\cos\gamma+1)}{\beta_{s}^{2}}-\frac{4i(\beta_{u}-\beta_{v})\csc^{2}\gamma}{\beta_{u}\beta_{v}}\right)\right.$
		$\displaystyle\left.+\sin\beta_{s}\left(\frac{12\beta_{u}\beta_{v}\sin^{2}\gamma}{\beta_{s}^{5}}-\frac{4\beta_{u}\beta_{v}\sin^{2}\gamma-4i(\beta_{u}-\beta_{v})(\cos\gamma+1)+8\cos\gamma}{\beta_{s}^{3}}+\frac{4\beta_{s}}{\beta_{u}\beta_{v}}+\frac{12\cos\gamma+4}{\beta_{s}}\right)\right]\biggl{\}}.$

Like the vector mode, the ORF of the scalar-longitudinal mode cannot be normalized. The auto-correlated response for the scalar-longitudinal mode of a single detector is

	$\displaystyle\Gamma^{L}_{II}(f)=$	$\displaystyle\lim_{\begin{subarray}{c}\gamma\rightarrow 0\\ \beta_{v}\rightarrow\beta_{u}\end{subarray}}\Gamma^{L}_{IJ}(f)=\frac{37}{12}-2\gamma_{E}-2\ln(2\beta_{u})$		(30)
		$\displaystyle+\frac{1}{2}\beta_{u}\mbox{Si}(2\beta_{u})+\frac{1}{4}\cos(2\beta_{u})$
		$\displaystyle+2\mbox{Ci}(2\beta_{u})-\frac{\sin(2\beta_{u})}{\beta_{u}}-\frac{2}{\beta^{2}_{u}}+\frac{\sin(2\beta_{u})}{\beta^{3}_{u}}.$

For short wavelengths, $\Gamma^{L}_{II}\propto\frac{\pi}{4}\beta_{u}$ and it grows faster than the logarithmically growing vector mode as $\beta_{u}$ increases. This is consistent with Eq.(A36) in 2008apj_Lee_ORF_nonGR and Eq.(41) in 2012prd_ORF_nonGR .