# Gravitational Waveform and Polarization from Binary Black Hole Inspiral in
Dynamical Chern-Simons Gravity: From Generation to Propagation

Zhao Li<EMAIL_ADDRESS>Department of Astronomy, University of
Science and Technology of China, Hefei, Anhui 230026, China,
and School of Astronomy and Space Science, University of Science and
Technology of China, Hefei 230026, China  Jin Qiao  Purple Mountain
Observatory, Chinese Academy of Sciences, Nanjing, 210023, China,
and School of Astronomy and Space Science, University of Science and
Technology of China, Hefei 230026, China  Tan Liu  MOE Key Laboratory of
Fundamental Physical Quantities Measurements, Hubei Key Laboratory of
Gravitation and Quantum Physics, PGMF, Department of Astronomy, and School of
Physics, Huazhong University of Science and Technology, Wuhan 430074, China
Tao Zhu  Institute for theoretical physics and cosmology, Zhejiang University
of Technology, Hangzhou, 310032, China,
United Center for Gravitational Wave Physics (UCGWP), Zhejiang University of
Technology, Hangzhou, 310032, China  Wen Zhao<EMAIL_ADDRESS>Department
of Astronomy, University of Science and Technology of China, Hefei, Anhui
230026, China,
and School of Astronomy and Space Science, University of Science and
Technology of China, Hefei 230026, China

###### Abstract

We calculate the gravitational waveform radiated from spinning black holes
(BHs) binary in dynamical Chern-Simons (dCS) gravity. The equation of motion
(EOM) of the spinning binary BHs is derived based on the modified Mathisson-
Papapetrou-Dixon equation for the spin-aligned circular orbits. The leading-
order effects induced by the dCS theory contain spin-spin interaction and
monopole-quadrupole interaction, which influences both the EOM of the binary
system and corresponding gravitational waveform at the second post-Newtonian
(PN) order (i.e., 2PN order). After reporting the waveforms, we investigate
the polarization modes of gravitational waves (GWs) in dCS theory. None of the
extra modes appears in this theory up to the considered PN order. Moreover,
since the time scale of the binary merger is much smaller than that of the
cosmological expansion, the parity-violating effect of the dCS theory does not
appear in the process of GW generation. However, during the process of GW
propagation, amplitude birefringence, a typical parity-violating effect, makes
plus and cross modes convert to each other, which modifies the gravitational
waveform at 1.5PN order.

##  I Introduction

General relativity (GR) is always considered as the most successful theory of
gravity. However, various difficulties of this theory are also well known. For
instance, in the theoretical side, GR has the singularity and quantization
problems Kiefer (2007). In the experimental side, all the observations in
cosmological scale indicate the existence of so-called dark matter and dark
energy Sahni (2005), which might mean that GR is invalid at this scale. For
these reasons, since Einstein proposed the theory of GR, many experiments have
tested the validity of this theory of gravity on various scales from sub-
millimeter-scale tests in the laboratory to tests at solar system and
cosmological scales Weinberg (1972); Poisson and Will (2014); Will (2014);
Hoyle  _et al._ (2001); Everitt  _et al._ (2011); Sabulsky  _et al._ (2019).
However, most of these experiments focus on the weak-field effects.

Gravitational wave (GW) is one of the most important predictions of GR, which
can only be generated in strong gravitational fields and hardly interacts with
matter, there is an excellent opportunity to test the theory of gravity in the
strong-field regime. In recent years, with the detection of the first GW event
from a binary black holes (BBH) system GW150914 Abbott  _et al._ (2016a) and
other GW events Abbott  _et al._ (2017, 2020, 2019a, 2021a); Abbott and other
(2021), the theory and detection of GWs have attracted much attention.
Numerous works have used GWs as a new probe for gravitational testing Abbott
_et al._ (2016b, 2019b, 2021b, 2021c).

The BBH systems radiate GWs during the process of inspiraling to merging
Poisson and Will (2014); Misner  _et al._ (1971); Maggiore (2008). Because the
binary has a very strong gravitational field in the pre-merger phase, the
radiated GW may encode the distinction between the real gravitational theory
and GR. Therefore, the detection of GWs generated from BBH is an essential
aspect of the experiment of gravitational theory Yunes and Pretorius (2009a);
Li  _et al._ (2012); Agathos  _et al._ (2014). With the improvement of
sensitivities of GW detectors, future GW tests of gravitational theory impose
higher requirements on the calculation of waveform templates Jaranowski and
Królak (2012). Constructing such accurate templates has motivated significant
research on post-Newtonian (PN) approximation method, which is applied to
compute the GW waveforms in the regime where the two bodies have large
separation. In the PN framework, the separation remains large with respect to
the radii of both objects, thus the bodies can be regarded effectively as
point particles.

The testing of GR by GW observations entails comparing the predictions of GW
signals in GR and those in the alternative theories and constraining their
differences by observations. Therefore, the choice of typical alternative
gravitational theory Brans and Dicke (1961); Jackiw and Pi (2003); Kanti  _et
al._ (1996); Horndeski (1974); Eling  _et al._ (2004); Moffat (2015); Gao
(2014); Gao  _et al._ (2019) and the calculation of GW waveforms in the theory
have crucial roles. For example, the gravitational waveforms have been derived
in Brans-Dicke (BD) gravity Zhang  _et al._ (2017a); Liu  _et al._ (2020);
Zhao  _et al._ (2019), Einstein-aether theory Zhang  _et al._ (2020); Lin  _et
al._ (2019); Zhao  _et al._ (2019), screen modified gravity Zhang  _et al._
(2016, 2017b); Liu  _et al._ (2018); Zhang  _et al._ (2019a), and Horndeski
theory Chowdhuri and Bhattacharyya (2022). Current and future GW observations
can constrain these gravitational theories with high accuracy Zhang  _et al._
(2019b); Niu and Zhao (2019); Zhang  _et al._ (2019c); Niu  _et al._ (2020,
2021); Shi  _et al._ (2022); Carson and Yagi (2020); Chamberlain and Yunes
(2017); Barausse  _et al._ (2016); Perkins  _et al._ (2021). Some kinds of
modified theories, including parity-violating effects, have been studied in
recent years Smith  _et al._ (2008); Yunes and Spergel (2009); Yunes  _et al._
(2010); Mirshekari  _et al._ (2012); Zhao  _et al._ (2020a); Wu  _et al._
(2022); Gong  _et al._ (2022); Wang  _et al._ (2013); Zhu  _et al._ (2013);
Qiao  _et al._ (2020); Zhu  _et al._ (2023); Qiao  _et al._ (2023); Boudet
_et al._ (2023); Bombacigno  _et al._ (2023), like Chern-Simons (CS) gravity
Jackiw and Pi (2003); Alexander and Yunes (2009) (consisting of dynamical and
non-dynamical cases), ghost-free parity-violating theory Qiao  _et al._
(2019); Zhao  _et al._ (2020b); Wang  _et al._ (2021, 2022); Zhao  _et al._
(2022), Nieh-Yan modified teleparallel gravity Li  _et al._ (2020, 2021),
parity-violating symmetric teleparallel gravities Li  _et al._ (2022a) and the
general spatial covariant gravities Takahashi and Soda (2009); Wang  _et al._
(2013); Wang (2017); Zhu  _et al._ (2013); Gao (2014); Gao  _et al._ (2019);
Zhu  _et al._ (2022, 2023). As a typical example, CS theory introduces a
coupling term between the pseudo-scalar and Pontryagin density in the
Einstein-Hilbert action, which causes the non-conservation of the CS
topological current. One of the most important predictions of CS theory is
that the amplitude of the left-hand circular polarization mode of GWs
increases (or decreases) during the propagation while the amplitude for the
right-hand mode decreases (or increases). This phenomenon is always called
amplitude birefringence of GWs Alexander and Yunes (2009).

blackThis work aims to derive the gravitational waveform in dynamical CS (dCS)
gravity. The GW polarization modes are investigated in both generation and
propagation processes. Although the binary equation of motion (EOM) and PN
waveform in the dCS theory are obtained by Refs. Yagi  _et al._ (2012a, b,
2016), the polarization of GWs is not studied, and the energy flux carried by
tensor radiation is absent. The latter can result in a problematic
parameterized post-Einsteinian (ppE) parameters. In this work, we will
recalculate the gravitational waveforms generated by the inspirals of BBH,
which are different from Refs. Yagi  _et al._ (2012a, b, 2016) in the
following three aspects. First, we derive the EOM of BBH from the modified
Mathisson-Papapetrou-Dixon (MPD) equation Bini and Geralico (2014), a powerful
tool to describe the motion of spinning point particle, which is based on the
effective field theory (EFT) approach Steinhoff (2011); Puetzfeld  _et al._
(2015) and is extended by Loutrel  _et al._ (2018); Loutrel and Yunes (2022)
to study the motion and precession of spinning particles in dCS theory.
Second, the investigation of GW polarization modes and improved energy flux
are included. Third, the propagation effect, i.e., amplitude birefringence, is
discussed. In conclusion, the improved ppE parameters and waveform including
propagation effects are obtained in this work.

This paper is organized as follows. In Section II, we briefly review the dCS
theory. In Section III, we focus on the spin-aligned case and obtain the EOM
of BBHs up to 2PN approximation. Section IV calculates the scalar and tensor
radiation fields by multipole moment formulae. The polarization modes of
generated GW are obtained through Newman-Penrose (NP) tetrad Will (2014);
Wagle  _et al._ (2019) in Section V. Furthermore, the energy flux carried by
radiation and the orbital evolution are given in Section VI. In Section VII,
we obtain the ready-used frequency-domain waveform. The birefringence effect
of GWs when propagating in the cosmological background is discussed in Section
VIII. Finally, we make a summary and discussion in Section IX.

Throughout the paper, we adopt the following conventions: We work in four
dimensions with metric signature $(-,+,+,+)$, Latin indices
$(a,b,c,\cdots,j,k,...)$ in the index list represent spatial indices, whereas
Greek indices $(\alpha,\beta,\cdots)$ represent spacetime indices, round
brackets around indices represent symmetrization, square brackets represent
anti symmetrization, $\partial_{\mu}$ represents a partial derivative,
$\nabla_{\mu}$ represents a covariant derivative,
$\Box^{2}_{\eta}\equiv\partial_{\mu}\partial^{\mu}$,
$\Box^{2}_{g}\equiv\nabla_{\mu}\nabla^{\mu}$, whereas
$\bm{\nabla}^{2}\equiv\partial_{j}\partial_{j}$, the Einstein summation
convention is employed, and we work in geometric units in which $c=G=1$, where
$c$ is the speed of light in the vacuum and $G$ is the gravitational constant.

##  II Dynamical Chern-Simons Gravity

###  II.1 Action and Field Equations

The full action of the dCS theory is Alexander and Yunes (2009)

$\displaystyle S=\int
d^{4}x\sqrt{-g}\left[\frac{1}{16\pi}R+\frac{\alpha}{4}\vartheta
R\hat{R}-\frac{\beta}{2}(\nabla_{\mu}\vartheta)(\nabla^{\mu}\vartheta)+\mathcal{L}_{m}\right],$
(1)

where the gravity is described by a pseudo scalar field $\vartheta$ and the
metric $g_{\mu\nu}$. As in redRef. Yagi  _et al._ (2012a), we do not consider
the potential of the scalar field. In Eq. (1), $g$ is the determinant of the
metric $g_{\mu\nu}$ and $R$ is the Ricci scalar. $\mathcal{L}_{m}$ is the
Lagrangian density of the matter field. $\alpha$ and $\beta$ are the coupling
parameters.
$R\hat{R}\equiv(1/2)\varepsilon^{\rho\sigma\alpha\beta}R_{\nu\mu\rho\sigma}R^{\mu\nu}_{\
\ \alpha\beta}$ is the Potryagin density, with $R_{\nu\mu\rho\sigma}$ being
the Riemann tensor and $\varepsilon^{\rho\sigma\alpha\beta}$ being the Levi-
Civitá tensor defined in terms of the antisymmetric symbol
$\epsilon^{\rho\sigma\alpha\beta}$ as
$\varepsilon^{\rho\sigma\alpha\beta}=(1/\sqrt{-g})\epsilon^{\rho\sigma\alpha\beta}$.

The variation of the total action (1) with respect to the metric $g_{\mu\nu}$
yields the modified field equations Jackiw and Pi (2003); Alexander and Yunes
(2009),

$R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R+16\pi\alpha
C_{\mu\nu}=8\pi\left[T_{\mu\nu}^{(m)}+T_{\mu\nu}^{(\vartheta)}\right],$ (2)

where $R_{\mu\nu}$ is Ricci tensor and $C_{\mu\nu}$ is Cotton tensor defined
as

$C^{\mu\nu}=-\varepsilon^{\rho(\mu|\alpha\beta}\left[\nabla_{\alpha}R^{\nu)}_{\
\beta}\right](\nabla_{\rho}\vartheta)-\hat{R}^{\kappa(\mu|\rho|\nu)}(\nabla_{\kappa}\nabla_{\rho}\vartheta).$
(3)

Note that the Cotton tensor $C_{\mu\nu}$ is traceless,
$g^{\mu\nu}C_{\mu\nu}=0$, and satisfies the Bianchi identity,
$\nabla^{\mu}C_{\mu\nu}=0$. $T_{\mu\nu}^{(m)}$ denotes the energy-momentum
tensor (EMT) of the matter fields, and $T_{\mu\nu}^{(\vartheta)}$ represents
the EMT of the coupled CS scalar field,

$T_{\mu\nu}^{(\vartheta)}=\beta\left[(\nabla_{\mu}\vartheta)(\nabla_{\nu}\vartheta)-\frac{1}{2}g_{\mu\nu}(\nabla_{\alpha}\vartheta)(\nabla^{\alpha}\vartheta)\right].$
(4)

The equation of the scalar field can also be derived by variation the action
(1) with respect to the scalar field $\vartheta$, which leads to

$\beta\Box^{2}_{g}\vartheta=-\frac{\alpha}{4}R\hat{R}.$ (5)

We would like to mention here that when the coupling constant $\beta=0$, the
total action (1) reduces to that of the non-dynamical CS gravity. In the non-
dynamical case, the scalar field equation (5) becomes an additional
differential constraint, i.e, the Pontryagin constraint on the space of
allowed solutions,

$R\hat{R}=0.$ (6)

In this paper we will not consider this case and only focus on the dCS gravity
in which the parameter $\beta\neq 0$.

###  II.2 Slowly-Rotating Black Hole

One of important features of the dCS theory is the parity violation induced by
the coupling between the scalar field $\vartheta$ and the Pontryagin density
$R\hat{R}$. Thus, the Pontryagin density in general disappears in spherically
symmetric spacetime. For this reason, the Schwarzschild black hole (BH) is
still an exact solution to the dCS theory Grumiller and Yunes (2008). The GWs
radiated from binary Schwarzschild BHs are the same as those in GR Alexander
and Yunes (2009). Since the Pontryagin density normally is not vanishing for
non-spherically symmetric gravitational systems, it is expected that the dCS
theory can induce significant contributions to the GW waveform generated by
binary spinning BHs. In this work, we focus on the gravitational radiation
from binary spinning BHs. For this purpose, one has to first derive the
spinning BH solution in dCS theory. Currently, the slowly rotating solution
have been found within the small-coupling and small-spin approximation Yunes
and Pretorius (2009b); Yagi  _et al._ (2012c), which describes a slowly-
rotating BH in the vacuum. The full metric $g_{\mu\nu}$ of the slowly-rotating
BH solution and the scalar field profile $\vartheta(r)$ are expanded in terms
of coupling constant $\alpha$ and dimensionless spin $\chi$ as follows,

$g_{\mu\nu}=g_{\mu\nu}^{\rm(K)}+\alpha^{\prime
2}\chi^{\prime}g^{(2,1)}_{\mu\nu}+\alpha^{\prime 2}\chi^{\prime
2}g^{(2,2)}_{\mu\nu}+\cdots$ (7)
$\vartheta=\alpha^{\prime}\chi^{\prime}\vartheta^{(1,1)}+\alpha^{\prime}\chi^{\prime
2}\vartheta^{(1,2)}+\cdots.$ (8)

The bookkeeping parameters $\alpha^{\prime}$ and $\chi^{\prime}$ denote the
expansion orders. In Boyer-Lindquist coordinate Boyer and Lindquist (1967),
$(\tilde{t},\tilde{r},\tilde{\theta},\tilde{\varphi})$, $g^{(2,1)}_{\mu\nu}$,
$g^{(2,2)}_{\mu\nu}$, $\vartheta^{(1,1)}$ and $\vartheta^{(1,2)}$ are
undetermined functions of $\tilde{r}$ and $\tilde{\theta}$. The zero-order
approximation $g_{\mu\nu}^{\rm(K)}$ is the Kerr metric Kerr (1963). The
superscripts $(m,n)$ represent the expansion order of $\alpha^{\prime}$ and
$\chi^{\prime}$. These undetermined functions can be solved by substituting
these expansions, (7) and (8), into the field equations, (2) and (5) Yagi  _et
al._ (2012c). Up to the required PN order, the scalar field profile at each
order are given by

$\vartheta^{(1,1)}=-\frac{5}{8}\chi\frac{\alpha}{\beta}\frac{\cos\tilde{\theta}}{\tilde{r}^{2}},\quad\text{and}\quad\vartheta^{(1,2)}=0.$
(9)

And the metric field are

$g_{\tilde{t}\tilde{\varphi}}^{(2,1)}=10\frac{\pi\alpha^{2}}{\beta
m_{0}^{4}}\left(\frac{m_{0}}{\tilde{r}}\right)^{4}m_{0}\chi\sin^{2}\tilde{\theta},$
(10)

and

$g^{(2,2)}_{\tilde{t}\tilde{t}}=g^{(2,2)}_{\tilde{r}\tilde{r}}=\frac{1}{\tilde{r}^{2}}g^{(2,2)}_{\tilde{\theta}\tilde{\theta}}=\frac{1}{\tilde{r}^{2}\sin^{2}\tilde{\theta}}g^{(2,2)}_{\tilde{\varphi}\tilde{\varphi}}=\frac{201}{112}\frac{\pi\alpha^{2}}{\beta
m_{0}^{4}}\left(\frac{m_{0}}{\tilde{r}}\right)^{3}\chi^{2}(3\cos^{2}\tilde{\theta}-1).$
(11)

The scalar and metric solutions (9), (10), and (11) are of $\mathcal{O}(\chi)$
order and $\mathcal{O}(\chi^{2})$ order, respectively. $m_{0}$ is the mass
parameter of the slowly-rotating BH. Other components in the metric function
are zero up to $\mathcal{O}(\alpha^{2}\chi^{2})$ order.

###  II.3 Scalar and Tensor Perturbation

The full metric is decomposed into a flat metric background $\eta_{\mu\nu}$
and a full metric perturbation $H_{\mu\nu}$,

$g_{\mu\nu}=\eta_{\mu\nu}+H_{\mu\nu},$ (12)

where the perturbation $H_{\mu\nu}$ is further decomposed into the GR part
$h_{\mu\nu}$ and a metric deformation away from GR $k_{\mu\nu}$ Yagi  _et al._
(2012a), i.e.,

$H_{\mu\nu}=h_{\mu\nu}+k_{\mu\nu}.$ (13)

Accordingly, the inverse metric can be approximately written in the form of
Misner  _et al._ (1971)

$g^{\mu\nu}\approx\eta^{\mu\nu}-H^{\mu\nu}+H^{\mu\alpha}H^{\nu}_{\alpha},$
(14)

up to the second-order perturbation. Similar to the calculation of GW
radiation from binary system in GR with the relaxed Einstein equation, it is
convenient to define the standard trace-reversed metric perturbation
$\bar{H}^{\mu\nu}$ as Landau and Lifchitz (2004); Maggiore (2008),

$\bar{H}^{\mu\nu}\equiv\eta^{\mu\nu}-\sqrt{-g}g^{\mu\nu},$ (15)

and here one requires $\bar{H}^{\mu\nu}$ satisfying the Lorenz gauge,

$\partial_{\mu}\bar{H}^{\mu\nu}=0.$ (16)

At the same time, the trace-reversed metrics $\bar{h}_{\mu\nu}$ and
$\bar{k}_{\mu\nu}$ are defined as

$\bar{h}^{\mu\nu}=h^{\mu\nu}-\frac{1}{2}\eta^{\mu\nu}h+\frac{1}{2}hh^{\mu\nu}+\frac{1}{4}\eta^{\mu\nu}h_{\alpha\beta}h^{\alpha\beta}-\frac{1}{8}\eta^{\mu\nu}h^{2}-h^{\mu\alpha}h_{\alpha}^{\nu},$
(17)

and

$\bar{k}^{\mu\nu}=k^{\mu\nu}-\frac{1}{2}\eta^{\mu\nu}k,$ (18)

respectively. The traces of GR perturbation then given by
$h\equiv\eta^{\mu\nu}h_{\mu\nu}$ and
$\bar{h}\equiv\eta^{\mu\nu}\bar{h}_{\mu\nu}$, and the traces of metric
deformation are defined as $k\equiv\eta^{\mu\nu}k_{\mu\nu}$ and
$\bar{k}\equiv\eta^{\mu\nu}\bar{k}_{\mu\nu}$.

As shown in Section II.2, the scalar and gravitational fields should be of
order $\mathcal{O}(\alpha)$ and $\mathcal{O}(\alpha^{2})$, respectively.
Therefore, to derive the scalar equation, the terms who is of order
$\mathcal{O}(\vartheta h)$ and $\mathcal{O}(h^{2})$ are retained, and that of
order $\mathcal{O}(hk)$ are removed. And to derive the tensor equation, the
perturbation order is of $\mathcal{O}(\vartheta^{2})$ and $\mathcal{O}(k)$.

With the above setups, one is able to get the second-order perturbated
equations by substituting decomposition (12) into Eqs. (2) and (5). For the
scalar field equation, it becomes Yagi  _et al._ (2016)

$\Box^{2}_{\eta}\vartheta=-16\pi\sigma,$ (19)

where the source term $\sigma$ is given by

$16\pi\sigma=\frac{1}{2}(\partial_{\mu}h)(\partial^{\mu}\vartheta)-(\partial_{\mu}h^{\mu\nu})(\partial_{\nu}\vartheta)-h^{\mu\nu}(\partial_{\mu}\partial_{\nu}\vartheta)+\frac{\alpha}{4\beta}\epsilon^{\rho\sigma\alpha\beta}(\partial_{\rho}\partial^{\lambda}h_{\mu\sigma})(\partial_{\alpha}\partial^{\mu}h_{\lambda\beta}-\partial_{\alpha}\partial_{\lambda}h^{\mu}_{\beta}).$
(20)

As we have mentined in the above, the full metric perturbation $H_{\mu\nu}$ is
divided into two parts: the GR part $h_{\mu\nu}$ and metric deformation
$k_{\mu\nu}$ away from GR. The former satisfies the quadratic-linearized
relaxed Einstein equation Landau and Lifchitz (2004); Maggiore (2008), while
the latter one satisfies Yagi  _et al._ (2016)

$\Box^{2}_{\eta}\bar{k}_{\mu\nu}=-16\pi K_{\mu\nu},$ (21)

where the source term $K_{\mu\nu}$ is

$K_{\mu\nu}=kT_{\mu\nu}^{(m)}+\left(1+h\right)\delta
T_{\mu\nu}^{(m)}+T_{\mu\nu}^{(\vartheta)}+\frac{1}{16\pi}\tilde{\Lambda}_{\mu\nu}-2\alpha\tilde{C}_{\mu\nu}.$
(22)

Here the quantity $\delta T_{\mu\nu}^{(m)}$ stands for the perturbation to the
EMT of the matter fields, and the tensor $\tilde{\Lambda}^{(2)}_{\mu\nu}$ is
defined as

$\tilde{\Lambda}_{\mu\nu}\approx 16\pi(-g)\tilde{t}_{\mu\nu}^{\rm
LL}-(\partial^{\alpha}\partial^{\beta}\bar{h}_{\mu\nu})\bar{k}_{\alpha\beta}-(\partial^{\alpha}\partial^{\beta}\bar{k}_{\mu\nu})\bar{h}_{\alpha\beta}+(\partial^{\beta}\bar{h}_{\alpha\mu})(\partial^{\alpha}\bar{k}_{\beta\nu})+(\partial^{\beta}\bar{k}_{\alpha\mu})(\partial^{\alpha}\bar{h}_{\beta\nu}),$
(23)

up to the linear order of metric deformation. The quantity $\tilde{t}^{\rm
LL}_{\mu\nu}$ in the expression of $\tilde{\Lambda}^{(2)}_{\mu\nu}$ denotes
the linearized Landau-Lifshitz energy-momentum pseudotensor, which is given by
Landau and Lifchitz (2004)

$\displaystyle 16\pi(-g)\tilde{t}^{\rm LL}_{\mu\nu}$
$\displaystyle=(\partial_{\alpha}\bar{h}_{\mu}^{\lambda})(\partial^{\alpha}\bar{k}_{\lambda\nu})+(\partial_{\alpha}\bar{k}_{\mu}^{\lambda})(\partial^{\alpha}\bar{h}_{\lambda\nu})+\eta_{\mu\nu}(\partial_{\beta}\bar{h}^{\alpha\lambda})(\partial_{\alpha}\bar{k}_{\lambda}^{\beta})$
(24)
$\displaystyle-(\partial_{\mu}\bar{h}^{\beta\lambda})(\partial_{\beta}\bar{k}_{\lambda\nu})-(\partial_{\mu}\bar{k}^{\beta\lambda})(\partial_{\beta}\bar{h}_{\lambda\nu})-(\partial_{\nu}\bar{h}^{\beta\lambda})(\partial_{\beta}\bar{k}_{\lambda\mu})-(\partial_{\nu}\bar{k}^{\beta\lambda})(\partial_{\beta}\bar{h}_{\lambda\mu})$
$\displaystyle+(1/2)(\partial_{\mu}\bar{h}^{\lambda\sigma})(\partial_{\nu}\bar{k}_{\lambda\sigma})+(1/2)(\partial_{\mu}\bar{k}^{\lambda\sigma})(\partial_{\nu}\bar{h}_{\lambda\sigma})-(1/4)(\partial_{\mu}\bar{h})(\partial_{\nu}\bar{k})$
$\displaystyle-(1/4)(\partial_{\mu}\bar{k})(\partial_{\nu}\bar{h})-(1/2)\eta_{\mu\nu}(\partial_{\alpha}\bar{h}^{\lambda\sigma})(\partial^{\alpha}\bar{k}_{\lambda\sigma})+(1/4)\eta_{\mu\nu}(\partial_{\alpha}\bar{h})(\partial^{\alpha}\bar{k}).$

Up to the linear order of the scalar field $\vartheta$ and the GR perturbation
$h_{\mu\nu}$, the C-tensor becomes

$\tilde{C}_{\mu\nu}=\frac{1}{2}\eta_{(\mu|\lambda}\epsilon^{\rho\lambda\alpha\beta}\partial_{\alpha}(\Box^{2}_{\eta}h_{\nu)\beta})(\partial_{\rho}\vartheta)-\frac{1}{2}\eta_{(\nu|\lambda}\epsilon^{\rho\lambda\alpha\beta}\left(\partial_{\alpha}\partial_{\mu)}h^{\kappa}_{\beta}-\partial_{\alpha}\partial^{\kappa}h_{\beta\mu}\right)(\partial_{\kappa}\partial_{\rho}\vartheta).$
(25)

The source $\delta T_{\mu\nu}^{(m)}$ in Eq. (22) contains contributions
proportional to the density, pressure and velocity of matter, which has
_compact support_ , that is, is nonzero only in finite spatial regions where
the matter resides. Inversely, other sources in Eq. (22) contain contributions
proportional to products of fields and their derivatives. These contributions
extend over all spacetime, although they generally decrease with distance from
the sources Will (2014). Differently, the source of scalar field (20) only
involves non-compact support terms, while the compact-support terms will be
derived from the EFT approach (see subsection II.4). The near-zone solution
and the motion of the binary system are mainly determined by the compact
support terms, while the gravitational radiation is sourced from both the
compact and non-compact support terms.

###  II.4 Modified MPD Equations

In this work, the BHs are modelled as a spinning point mass using the EFT
method extended by Ref. Loutrel  _et al._ (2018) to dCS gravity. The motion of
BBH is influenced by the mass and the spin angular momentum (SAM) of the point
particles. By considering a set of symmetries associated with theory,
specifically (1) re-parameterization invariance of proper time $\tau$, (2)
Lorentz invariance, and (3) diffeomorphism invariance, one can write the
quadrupole-approximated Lagrangian of a spinning point mass as Steinhoff
(2011); Loutrel  _et al._ (2018)

$L_{m}(u^{\mu},\Omega^{\mu\nu},R_{\alpha\beta\gamma\delta})=p_{\mu}u^{\mu}+\frac{1}{2}S_{\alpha\beta}\Omega^{\alpha\beta}-\frac{1}{6}J^{\alpha\beta\gamma\delta}R_{\alpha\beta\gamma\delta}.$
(26)

In Eq. (26), $u_{\mu}$, $p_{\mu}$, $\Omega^{\alpha\beta}$ and
$S^{\alpha\beta}$ are the 4-velocity, 4-momentum, angular velocity tensor and
SAM tensor, respectively. The residual degree of freedom (DOF) of the SAM
tensor is removed by the Tulczyjew spin supplementary condition Blanchet
(2014); Loutrel  _et al._ (2018), $S^{\mu\nu}p_{\mu}=0$. The Dixon’s
quadrupole moment is defined as Steinhoff (2011); Bini and Geralico (2014);
Puetzfeld  _et al._ (2015); Loutrel  _et al._ (2018)

$J^{\alpha\beta\gamma\delta}\equiv-\frac{3}{m_{0}}u^{[\alpha}S^{\beta]\lambda}S_{\
\
\lambda}^{[\gamma}u^{\delta]}=\frac{3}{m_{0}}\left(u^{[\alpha}S^{\beta]}S^{[\gamma}u^{\delta]}-S^{2}u^{[\alpha}g^{\beta][\gamma}u^{\delta]}\right),$
(27)

where $S^{2}\equiv S^{\mu}S_{\mu}$ and $S^{\mu}$ is the spin 4-vector defined
as
$S^{\mu}\equiv-(1/2)\epsilon^{\alpha\beta\gamma\mu}u_{\gamma}S_{\alpha\beta}$,
automatically satisfying the Tulczyjew condition. Then, the residual DOF is
removed by gauge condition, $S^{\mu}u_{\mu}=0$.

However, dCS theory brings additional modification to the Lagrangian Loutrel
_et al._ (2018). The scalar field modifies the Lagrangian (26), and the
slowly-rotating BH metric described by Eqs. (9-11) violates the standard Kerr
metric, an additional contribution to the Newtonian gravitational potential.
By requiring the Lagrangian to be invariant under (1) a shift of the scalar
field and (2) a parity transformation, the modified Lagrangian function is
finally given by

$L_{m}(u^{\mu},\Omega^{\mu\nu},R_{\alpha\beta\gamma\delta},\nabla^{\alpha}\vartheta)=p_{\mu}u^{\mu}+\frac{1}{2}S_{\alpha\beta}\Omega^{\alpha\beta}\\\
-\frac{1}{6}(1+\delta
C_{Q})J^{\alpha\beta\gamma\delta}R_{\alpha\beta\gamma\delta}+\frac{1}{m_{0}^{2}}\delta
C_{\vartheta}\hat{S}_{\mu}^{\ \alpha}u^{\mu}(\nabla_{\alpha}\vartheta),$ (28)

where the correction constants, $\delta C_{\vartheta}$ and $\delta C_{Q}$, are
determined by comparing the near-zone and the slowly-rotating BH solution. The
dual tensor of SAM is $\hat{S}_{\mu\nu}\equiv(1/2)\epsilon_{\mu\nu}^{\ \
\lambda\rho}S_{\lambda\rho}$.

From the modified Lagrangian (28), the modified MPD equations are given by
Loutrel  _et al._ (2018)

$\displaystyle u^{\alpha}\nabla_{\alpha}p_{\mu}$
$\displaystyle=\frac{1}{2}S_{\alpha\beta}u^{\lambda}R^{\alpha\beta}_{\ \
\lambda\mu}+\frac{1}{m_{0}^{2}}\delta C_{\vartheta}\hat{S}_{\lambda}^{\
\alpha}u^{\lambda}(\nabla_{\mu}\nabla_{\alpha}\vartheta)-\frac{1}{6}(1+\delta
C_{Q})J^{\alpha\beta\gamma\delta}\nabla_{\mu}R_{\alpha\beta\gamma\delta},$
(29a) $\displaystyle u^{\alpha}\nabla_{\alpha}S_{\mu\nu}$
$\displaystyle=2p_{[\mu}u_{\nu]}+\frac{2}{m_{0}^{2}}\delta
C_{\vartheta}u^{\lambda}\hat{S}_{\lambda[\mu}(\nabla_{\nu]}\vartheta)+\frac{4}{3}(1+\delta
C_{Q})J^{\rho\lambda\alpha}_{\ \ \ \ [\mu}R_{\nu]\alpha\rho\lambda}.$ (29b)

These two equations reduce to standard MPD equations by setting $\delta
C_{\vartheta}$ and $\delta C_{Q}$ to zeros Bini and Geralico (2014). In a BBH
system, Eq. (29a) determines the acceleration of the particle in dCS gravity
and Eq. (29b) is the spin-precession equation for the particle.

From the modified Lagrangian (28), the modified EMT of a spinning point mass
is written as Blanchet (2014); Loutrel  _et al._ (2018)

$T_{(m)}^{\mu\nu}=\int\frac{d\tau}{\sqrt{-g}}\left[m_{0}u^{\mu}u^{\nu}\delta^{(4)}-\nabla_{\lambda}(S^{\lambda(\mu}u^{\nu)}\delta^{(4)})\\\
-\frac{2}{3}(1+\delta
C_{Q})\nabla_{\alpha}\nabla_{\beta}(J^{\mu\alpha\beta\nu}\delta^{(4)})\right](x^{\mu}-z^{\mu}(\tau)).$
(30)

The three terms in the above EMT tensor act as mass, spin and quadrupole
source terms, respectively. $x^{\mu}$ is the coordinate of a point and
$z^{\mu}(\tau)$ is the trajectory of the point mass. Some terms describing
scalar-spin coupling and quadrupole-spacetime coupling are omitted because
they do not contribute to the gravitational field after regularization Loutrel
_et al._ (2018).

EFT shows that the point-particle approximation brings an additional compact
support term to the scalar field equation (5) Yagi  _et al._ (2016); Loutrel
_et al._ (2018),

$\rho_{\vartheta}=-\frac{\delta C_{\vartheta}}{4\pi\beta m_{0}^{2}}\int
d\tau\hat{S}_{\mu}^{\
\alpha}u^{\mu}\nabla_{\alpha}\left[\frac{\delta^{(4)}(x^{\mu}-z^{\mu})}{\sqrt{-g}}\right].$
(31)

This quantity acts as an effective dipolar source for the scalar field.
$\delta^{(4)}(x^{\mu}-z^{\mu})$ is four-dimensional Dirac delta function. Then
the field equation (19) is rewritten as

$\Box^{2}_{\eta}\vartheta=-16\pi\sigma+4\pi\rho_{\vartheta}.$ (32)

This compact support term will play an important role in the near-zone
solution and radiated waveform.

##  III Equation of Motion

In this Section, The EOM of the BBH system with quasi-circular orbit is
derived up to 2PN approximation. The SAMs of BHs are assumed to be aligned to
the orbital angular momentum of the system. At this time, the precession
equation becomes $u^{\alpha}\nabla_{\alpha}S_{\mu\nu}=0$, indicating that the
SAM of BHs is conserved. First, the approximate solution of the scalar field
and the PN metric in the near zone are derived for preparation. Second, the
2PN EOM of BBH is presented. Third, the modified relative EOM, Kepler’s third
law, and conserved energy are given in subsequent subsections.

###  III.1 Zones

This paper will discuss the generation and propagation of GWs on multiple
scales. For the convenience of later discussions, we divide the whole space
into four zones: _inner zone_ , _near zone_ , black _wave zone_ , and
_propagation zone_. The inner zones are centered on each object of the binary
system. Their radiuses are much smaller than the typical orbital radius of the
binary system. The gravitational fields in these zones are too strong to apply
the PN approximation and linear GW theory. The near zone is centered on the
binary’s center of mass, whose radius is comparable with the typical
wavelength of emitted GWs. The wave zone is also centered on the binary, whose
radius is about several times the wavelength but far less than the distance
between galaxies. The background can be considered flat, and GW propagates
freely in this zone. The propagation zone includes all other regions far away
from the source. The cosmological background must be considered when GWs
propagate in this zone. In dCS theory, the homogeneous background scalar field
induces the amplitude birefringence effect of GWs.

###  III.2 Near-Zone Solution: Scalar

The scalar field generated by a slowly-rotating BH acts like a magnetic
dipole. For simplification, we define the scalar dipole moment (SDM) of one BH
in a binary system (the $A$-th BH) as

$\mu^{\alpha}_{A}\equiv\frac{\delta C_{\vartheta}}{4\pi\beta
m_{A}^{2}}S_{A}^{\alpha},$ (33)

where $m_{A}$ and $S_{A}^{\alpha}$ are the mass and SAM of the $A$-th BH. The
SDM vector is constrained by gauge condition
$\mu_{A}^{\alpha}u^{A}_{\alpha}=0$. In the Newtonian limit, the proper time of
a particle $\tau$ is just the coordinate time $t$. Thus the time component of
4-velocity is $u^{A}_{0}=-1$, and the spatial components are replaced by
$v_{A}^{i}\equiv u_{A}^{i}/u_{A}^{0}$. Such that the time component of SDM is
$\mu_{A}^{0}=\mu_{A}^{i}v^{A}_{i}$. Then the compact-support source
$\rho_{\vartheta}$ (88) is written as Yagi  _et al._ (2016)

$\rho_{\vartheta}=\rho_{1}+\rho_{2}+\rho_{3}.$ (34)

The three terms in Eq. (34) are

$\displaystyle\rho_{1}$
$\displaystyle=\sum_{A}\mu_{A}^{i}\partial_{i}\delta^{(3)}(\mathbf{x}-\mathbf{z}_{A}(t))\sim\mathcal{O}(1/c^{4}),$
(35a) $\displaystyle\rho_{2}$
$\displaystyle=\sum_{A}\mu_{A}^{i}\dot{v}_{A}^{i}\delta^{(3)}(\mathbf{x}-\mathbf{z}_{A}(t))\sim\mathcal{O}(1/c^{6}),$
(35b) $\displaystyle\rho_{3}$
$\displaystyle=-\sum_{A}\mu_{A}^{i}v_{A}^{i}v_{A}^{j}\partial_{j}\delta^{(3)}(\mathbf{x}-\mathbf{z}_{A}(t))\sim\mathcal{O}(1/c^{6}),$
(35c)

where $\mathbf{x}$ and $\mathbf{z}_{A}$ are the spatial components of
4-vectors $x^{\mu}$ and $z_{A}^{\mu}$. In above equations, $c$ is the light
speed in vacuum, indicating the magnitude of these terms.
$\delta^{(3)}(\mathbf{x}-\mathbf{z})$ is three-dimensional Dirac delta
function. With the help of the decomposition (35), the near-zone scalar
solution of BBH is mainly sourced by the first term $\rho_{1}$ (35a). We
finally obtain

$\vartheta^{\rm(B)}\approx\int\rho_{1}(t,\mathbf{x}^{\prime})\frac{1}{|\mathbf{x}-\mathbf{x}^{\prime}|}d^{3}\mathbf{x}^{\prime}=\sum_{A}\frac{\hat{\bm{n}}_{A}\cdot\bm{\mu}_{A}}{r_{A}^{2}},$
(36)

where the distance between point $\mathbf{x}$ and $A$-th BH is defined as
$r_{A}\equiv|\mathbf{x}-\mathbf{z}_{A}|$, and the unit directional vector is
defined as $\hat{\bm{n}}_{A}\equiv(\mathbf{x}-\mathbf{z}_{A})/r_{A}$.

On the other hand, the scalar field for the binary system is now just the sum
of the individual contributions, specifically

$\vartheta^{\rm(B)}=-\frac{5}{8}\sum_{A}\frac{\alpha}{\beta
m_{A}^{2}}\frac{\hat{\bm{n}}_{A}\cdot\hat{\bm{S}}_{A}}{r_{A}^{2}}.$ (37)

Comparing Eqs. (36) and (37), we can determine the value of the correction
constant,

$\delta C_{\vartheta}=-\frac{5}{2}\pi\alpha.$ (38)

The spatial components of SDM is re-expressed as

$\bm{\mu}_{A}=-\frac{5}{8}\frac{\alpha}{\beta m_{A}^{2}}\bm{S}_{A}.$ (39)

The SDM of a slowly-rotating BH is similar to the magnetic dipole moment
produced by electron spin Yagi  _et al._ (2012a), providing the spin-spin
coupling between two BHs. This completes the computation of the near-zone
scalar field of the binary system.

###  III.3 Near-Zone Solution: Tensor

We turn to investigate the near-zone solution of the metric. Up to 2PN order,
the PN metric is given by Blanchet  _et al._ (1995, 1998); Blanchet and Faye
(2001); Blanchet  _et al._ (2002); Blanchet (2014)

$\displaystyle h_{00}^{\rm(B)}$
$\displaystyle=2\bar{U}-2\bar{U}^{2}+8\left[\bar{X}-\bar{V}_{i}\bar{V}_{i}+(1/6)\bar{U}^{3}\right],$
(40a) $\displaystyle h_{0i}^{\rm(B)}$
$\displaystyle=-4\bar{V}_{i}-8\bar{R}_{i},$ (40b) $\displaystyle
h_{ij}^{\rm(B)}$
$\displaystyle=2\delta_{ij}\bar{U}+\left(2\delta_{ij}\bar{U}^{2}+4\bar{W}_{ij}\right).$
(40c)

In Eq. (40), several potentials
$\bar{\mathbf{X}}=\\{\bar{U},\bar{V}_{i},\bar{X},\bar{R}_{i},\bar{W}_{ij}\\}$
are defined as the integrations of source terms (22). The potential $\bar{U}$
reduces to the Newtonian potential at leading order. Because the matter EMT
(30) contains mass, spin, and quadrupole parts, these potentials are also
correspondingly divided into these three parts, i.e.
$\bar{\mathbf{X}}=\bar{\mathbf{X}}^{\rm(M)}+\bar{\mathbf{X}}^{\rm(S)}+\bar{\mathbf{X}}^{\rm(Q)}$,
which are not shown here. The readers can find the explicit expressions of
$\bar{\mathbf{X}}^{\rm(M)}$ in Ref. Blanchet  _et al._ (1998),
$\bar{\mathbf{X}}^{\rm(S)}$ in Ref. Faye  _et al._ (2006), and
$\bar{\mathbf{X}}^{\rm(Q)}$ in Refs. Buonanno  _et al._ (2013); Bohé  _et al._
(2015). For instance, the monopole, dipole, and quadrupole sectors of
Newtonian potential are

$\bar{U}^{\rm(M)}=\sum_{A}\frac{m_{A}}{r_{A}}+\text{(higher PN order without
spin effects)}\sim\mathcal{O}(1/c^{2}),$ (41)
$\bar{U}^{\rm(S)}=-2\varepsilon_{ijk}\sum_{A}v_{A}^{i}S_{A}^{j}\partial_{k}\left(\frac{1}{r_{A}}\right)\sim\mathcal{O}(1/c^{5}),$
(42)

and

$\bar{U}^{\rm(Q)}=-\frac{3}{2}\sum_{A}\left(\frac{m_{A}}{r_{A}}\right)^{3}\frac{1}{m_{A}^{4}}\left[(\hat{\bm{n}}_{A}\cdot\bm{S}_{A})^{2}-\frac{1}{3}(\bm{S}_{A}\cdot\bm{S}_{A})\right]\sim\mathcal{O}(1/c^{6}).$
(43)

As we have shown, the quadratic-spin sector of 2PN Newtonian potential in GR
is proportional to $S_{A}^{2}/r_{A}^{3}$.

The lowest-order dCS modification to the metric perturbation (40a) is just

$\bar{U}\rightarrow\bar{U}+\delta U.$ (44)

$\delta U$ is the retarded potential integration of the non-GR sector of
matter EMT in Eq. (30). The $00$ component of the leading-order contribution
is

$\delta T_{00}^{(m)}=\frac{2}{3}\sum_{A}\delta
C_{Q}J_{A}^{i0k0}\partial_{i}\partial_{k}\delta^{(3)}(\mathbf{x}-\mathbf{z}_{A}(t)).$
(45)

Thus, by integrating the EMT (45), we eventually find

$\displaystyle\delta U\equiv\int\frac{\delta
T^{(m)}_{00}}{|\mathbf{x}-\mathbf{x}^{\prime}|}d^{3}\mathbf{x}^{\prime}$
$\displaystyle=-\frac{3}{2}\sum_{A}\frac{\delta
C_{Q}}{m_{A}}\left(S_{A}^{i}S_{A}^{j}-\frac{1}{3}\delta_{ij}S_{A}^{k}S_{A}^{k}\right)\frac{\hat{n}_{A}^{i}\hat{n}_{A}^{j}}{r_{A}^{3}}$
(46) $\displaystyle=-\frac{3}{2}\sum_{A}\frac{\delta
C_{Q}}{m_{A}^{4}}\left(\frac{m_{A}}{r_{A}}\right)^{3}\left[(\hat{\bm{n}}_{A}\cdot\bm{S}_{A})^{2}-\frac{1}{3}(\bm{S}_{A}\cdot\bm{S}_{A})\right]\sim\mathcal{O}(1/c^{6}).$

The modified potential function $\delta U$ in Eq. (46) is similar to
quadrupole gravitational potential $\bar{U}^{\rm(Q)}$, which contains
$\mathcal{O}(S^{2})$ terms and enters 2PN correction. The near-zone metric
deformation tensors are

$k^{\rm(B)}_{00}=2\delta U,\quad\text{and}\quad
k^{\rm(B)}_{0i}=k^{\rm(B)}_{ij}=0.$ (47)

at the leading PN order. We find that the $\delta U$ is proportional to
$\delta C_{Q}\cdot S_{A}^{2}/r_{A}^{3}$, which is similar to
$\bar{U}^{\rm(Q)}$ except a constant $\delta C_{Q}$, which is different for
each object. Such that, the quadratic-spin effect and the dCS modification are
both 2PN higher than the Newtonian potential, $m_{A}/r_{A}$. Other potential
in the PN metric, $\\{\bar{V}_{i},\bar{X},\bar{R}_{i},\bar{W}_{ij}\\}$, are
higher PN-order effects relative to the Newtonian potential. The dCS
modification to these potential functions,
$\\{\delta{V}_{i},\delta{X},\delta{R}_{i},\delta{W}_{ij}\\}$, are 2PN effects
relative to leading terms of
$\\{\bar{V}_{i},\bar{X},\bar{R}_{i},\bar{W}_{ij}\\}$. Therefore, these
corrections enter higher order correction. Accordingly, only the correction to
the potential $\bar{U}$ is necessary to be considered when deriving the near-
zone solution and EOM.

At the same time, the modified potential $\delta U$ also can be defined from
the slowly-rotating BH solution as $g_{\tilde{t}\tilde{t}}=-1+2(\bar{U}+\delta
U)$, with $g_{\tilde{t}\tilde{t}}$ being shown in Eq. (11). In the harmonic
coordinates, which are identical to the Boyer-Lindquist coordinates at leading
PN order Loutrel  _et al._ (2018), the potential is written as

$\delta U=-\frac{201}{224}\sum_{A}\frac{\pi\alpha^{2}}{\beta
m_{A}^{4}}\chi_{A}^{2}\frac{m_{A}^{3}}{\tilde{r}_{A}^{3}}\left(3\cos^{2}\tilde{\theta}_{A}-1\right)=3\sum_{A}\delta
Q_{A}^{ij}\frac{\hat{n}_{A}^{i}\hat{n}_{A}^{j}}{\tilde{r}_{A}^{3}},$ (48)

where $\delta Q_{A}^{ij}$ is the quadrupole moment tensor of the $A$-th BH
induced by dCS modification, defined by

$\delta Q_{A}^{ij}=\frac{201}{224}\frac{\pi\alpha^{2}}{\beta
m_{A}^{4}}m_{A}^{3}\chi_{A}^{2}\left(\begin{array}[]{ccc}-1/3&0&0\\\
0&-1/3&0\\\ 0&0&2/3\end{array}\right).$ (49)

To find the value of constant $\delta C_{Q}$, we match the modified potential
(48) and Eq. (46) in the binary reference system, taking the center of mass of
the binary system as the origin and the vertical direction of the orbital
plane as the $z$-axis. Transforming the quadrupole moment tensors $\delta
Q_{A}^{ij}$ (49) from the Boyer-Lindquist coordinates to the binary reference
system, the modified potential is now written as

$\delta
U=\frac{603}{3584}\sum_{A}\frac{\zeta_{A}}{m_{A}^{4}}\left(\frac{m_{A}}{r_{A}}\right)^{3}\left[(\hat{\bm{n}}_{A}\cdot\bm{S}_{A})^{2}-\frac{1}{3}(\bm{S}_{A}\cdot\bm{S}_{A})\right].$
(50)

The redefined dimensionless coupling constant is defined as
$\zeta_{A}=16\pi\alpha^{2}/\beta m_{A}^{4}$. Comparing Eq. (46) with (50), we
find that the exact result of modification constant $\delta C_{Q}$ is

$\delta C_{Q}=-\frac{201}{1792}\zeta_{A},$ (51)

which depends on the typical mass of each object.

So far, the near-zone metric and scalar field are completely determined. On
the one hand, the near-zone fields are the sum of the contributions from the
individual objects. On the other hand, the near-zone solutions provide a way
to find the undetermined constants in EFT and modified MPD equations.

###  III.4 2PN EOM

The modified MPD equations (29) describe the motion and spin precession of BBH
in dCS gravity. As we have mentioned, the SAM of each object is conserved for
spin-aligned orbits. Thus, we only focus on the first equation (29a), which
describes how the SAMs, scalar fields, and Dixon’s quadrupole moments modify
the accelerations. Four kinds of forces appear in the EOM, connection
($\Gamma$) force, spin-orbit (SO) force, dipole-dipole (DD) force, and
monopole-quadrupole (MQ) force, defined as

$\displaystyle Y_{(\Gamma)}^{\mu}$
$\displaystyle\equiv-m\Gamma_{\lambda\rho}^{\mu}u^{\lambda}u^{\rho},$ (52a)
$\displaystyle Y_{\rm(SO)}^{\mu}$
$\displaystyle\equiv-\frac{1}{2}\epsilon_{\alpha\beta\gamma\delta}u^{\gamma}u^{\lambda}S^{\delta}R^{\alpha\beta\mu}_{\
\ \ \ \lambda},$ (52b) $\displaystyle Y_{\rm(DD)}^{\mu}$
$\displaystyle\equiv\frac{1}{m^{2}}\delta
C_{\vartheta}S^{\alpha}(\nabla^{\mu}\nabla_{\alpha}\vartheta),$ (52c)
$\displaystyle Y_{\rm(MQ)}^{\mu}$ $\displaystyle\equiv-\frac{1}{6}(1+\delta
C_{Q})J^{\alpha\beta\gamma\delta}\nabla^{\mu}R_{\alpha\beta\gamma\delta},$
(52d)

respectively. The Christoffel connection in Eq. (52a) is calculated from
metric, including the SO, SS, and MQ effects. The constants $\delta
C_{\vartheta}$ and $\delta C_{Q}$ are determined in the last two subsections
III.2 and III.3. Calculating the above geometric quantities up to the required
order, we can write down the explicit expression of the EOM at 2PN
approximation. It should be noted that, when deriving the EOM, the self field
is removed via Hadamard regularization Yagi  _et al._ (2012a); Blanchet
(2014), which is briefly introduced in Appendix A. We labeled these two BHs
with $1$ and $2$. At the Newtonian order and leading dCS order, the motion of
the $1$-st BH depends only on the mass, spin, and relative position of the
$2$-nd BH, and vice versa. In terms of the 3-velocity, $\bm{v}_{1}$, the
spatial components of 4-velocity, we can write down the acceleration of $1$-st
BH Loutrel  _et al._ (2018),

$\frac{dv_{1}^{k}}{dt}=a^{1,k}_{(\Gamma)}+a^{1,k}_{\rm(SO)}+a^{1,k}_{\rm(DD)}+a^{1,k}_{\rm(MQ)}.$
(53)

Each term in Eq. (53) can be divided into two parts: the GR part,
$\bar{a}^{1,k}_{(\cdots)}$, and the part from dCS theory, $\delta
a^{1,k}_{(\cdots)}$. Explicitly, they are expressed as

$\bar{a}^{1,k}_{(\Gamma)}=-(\bar{\Gamma}_{1,00}^{k}-2v_{1}^{j}\bar{\Gamma}_{1,0j}^{k}+v_{1}^{i}v_{1}^{j}\bar{\Gamma}_{1,ij}^{k})+v_{1}^{k}(\bar{\Gamma}_{1,00}^{0}-2v_{1}^{j}\bar{\Gamma}_{1,0j}^{0}+v_{1}^{i}v_{1}^{j}\bar{\Gamma}_{1,ij}^{0}),$
(54) $\delta a^{1,k}_{(\Gamma)}=\partial_{k}\delta U_{1},$ (55)
$\bar{a}^{1,k}_{\rm(SO)}=\frac{1}{2m_{1}}\epsilon_{ijl}S_{1}^{l}(2v_{1}^{j}\bar{R}_{1}^{i0k0}+\bar{R}_{1}^{ijk0}-v_{1}^{q}\bar{R}_{1}^{ijkq}),$
(56) $\delta a^{1,k}_{\rm(SO)}=0,$ (57) $\bar{a}^{1,k}_{\rm(DD)}=0,$ (58)
$\delta a^{1,k}_{\rm(DD)}=\frac{1}{m_{1}^{3}}\delta
C_{\vartheta}S_{1}^{j}(\partial_{k}\partial_{j}\vartheta_{1}),$ (59)
$\bar{a}^{1,k}_{\rm(MQ)}=-\frac{1}{2m_{1}^{2}}(S_{1}^{i}S_{1}^{j}-S_{1}^{2}\delta_{ij})(\partial_{k}\partial_{i}\partial_{j}\bar{U}_{1}),$
(60) $\delta a^{1,k}_{\rm(MQ)}=-\frac{1}{2m_{1}^{2}}\delta
C_{Q}(S_{1}^{i}S_{1}^{j}-S_{1}^{2}\delta_{ij})(\partial_{k}\partial_{i}\partial_{j}\bar{U}_{1}),$
(61)

where $\bar{\Gamma}_{1,\mu\nu}^{\lambda}$, $\bar{R}_{\mu\nu\rho\lambda}^{1}$,
and $\bar{U}_{1}$ are the GR parts of the Christoffel connection, Riemann
tensor, and Newtonian potential, Hadamard-regularized at the location of the
$1$-st BH. In the above expressions, $\vartheta_{1}$ and $\delta U_{1}$ are
Hadamard-regularized near-zone scalar field and modified Newtonian potential
around the $1$-st BH.

For the GR case, all of the non-spinning effects (0PN, 1PN, and 2PN) are
included in Eq. (56). The SO coupling, at linear spin order (54), is a 1.5PN
effect. The SS and MQ couplings, at quadratic-spin order, are both 2PN
effects. The former is included in (58), while the latter is included in Eqs.
(54) and (60). In particular, the dCS theory does not contribute to the SO
coupling in the BBH system. The DD interaction in Eq. (59) is proportional to
$S_{1}S_{2}$, which is the modification of the SS coupling. In conclusion, the
lowest-order modification induced by the dCS theory is at the 2PN order.

From the leading-order modified Newtonian potential (50) and near-zone scalar
field (37), the Hadamard-regularized values of $\bar{U}_{1}$, $\vartheta_{1}$,
and $\delta U_{1}$ at the position of the $1$-st BH are given by Appendix A,

$\bar{U}_{1}=\frac{m_{2}}{r},\quad\vartheta_{1}=\frac{\delta
C_{\vartheta}}{4\pi\beta}\frac{S_{2}^{k}}{m_{2}^{2}}\frac{\hat{n}^{k}}{r^{2}},\quad\text{and}\quad\delta
U_{1}=\frac{603}{3584}\zeta\frac{m}{m_{2}}\gamma^{3}\frac{1}{m_{2}^{4}}\left[(\hat{\bm{n}}\cdot\bm{S}_{2})^{2}-\frac{1}{3}S_{2}^{2}\right].$
(62)

$r\equiv|\mathbf{z}_{1}-\mathbf{z}_{2}|$ is the relative distance (orbital
radius) of the BBH, and $\hat{\bm{n}}\equiv\hat{\bm{n}}_{1}-\hat{\bm{n}}_{2}$
is the unit relative position vector. Additionally, we define the total mass
of the binary system $m\equiv m_{1}+m_{2}$ and $\gamma\equiv m/r$. The
dimensionless coupling normalized by the total mass is

$\zeta\equiv 16\pi\frac{\alpha^{2}}{\beta m^{4}}.$ (63)

Substituting Eq. (62) into (53)-(61), the dCS modification to the acceleration
is thus calculated explicitly as

$\displaystyle\delta\bm{a}^{1}$
$\displaystyle=\delta\bm{a}^{1}_{(\Gamma)}+\delta\bm{a}^{1}_{\rm(DD)}+\delta\bm{a}^{1}_{\rm(MQ)}$
(64)
$\displaystyle=\gamma^{4}\zeta\left\\{\frac{603}{3584}\left[\frac{1}{m_{2}}\frac{1}{m_{2}^{4}}\left[2(\hat{\bm{n}}\cdot\bm{S}_{2})\bm{S}_{2}-5(\hat{\bm{n}}\cdot\bm{S}_{2})^{2}\hat{\bm{n}}+S_{2}^{2}\hat{\bm{n}}\right]+\frac{m_{2}}{m_{1}^{2}}\frac{1}{m_{1}^{4}}\left[2(\hat{\bm{n}}\cdot\bm{S}_{1})\bm{S}_{1}-5(\hat{\bm{n}}\cdot\bm{S}_{1})^{2}\hat{\bm{n}}+S_{1}^{2}\hat{\bm{n}}\right]\right]\right.$
$\displaystyle\qquad\qquad\left.-\frac{75}{256}\frac{1}{m_{1}}\frac{1}{m_{1}^{2}}\frac{1}{m_{2}^{2}}\left[(\bm{S}_{1}\cdot\bm{S}_{2})\hat{\bm{n}}+(\hat{\bm{n}}\cdot\bm{S}_{1})\bm{S}_{2}+(\hat{\bm{n}}\cdot\bm{S}_{2})\bm{S}_{1}-5(\hat{\bm{n}}\cdot\bm{S}_{1})(\hat{\bm{n}}\cdot\bm{S}_{2})\hat{\bm{n}}\right]\right\\}.$

This modification (64) is very similar to the SS and MQ terms in GR, which are
shown in Refs. Thorne and Hartle (1985); Kidder (1995). Under Newton’s
gravity, SAM does not influence the motion of objects. However, in GR and its
modification, SAM will produce the components of gravitational force along the
SAM directions.

###  III.5 Relative EOM for Spin-Aligned Quasi-Circular Orbits

Now, let us focus on the spin-aligned quasi-circular orbits, in which the SAM
is perpendicular to the orbital plane, and the SAM of each object is
conserved. The spin-aligned assumption means
$\hat{\bm{n}}\cdot\bm{S}_{1}=\hat{\bm{n}}\cdot\bm{S}_{2}=0$. At the same time,
when the energy dissipation is not taken into consideration, the relative
orbit of the binary BHs is assumed to be circular, indicating that
$\hat{\bm{n}}\cdot\bm{v}=0$, where $\bm{v}$ is the relative velocity, defined
as $\bm{v}\equiv\bm{v}_{1}-\bm{v}_{2}$. The total acceleration of $1$-st BH
(64) thus becomes

$\delta\bm{a}^{1}=\zeta\gamma^{4}\Bigg{\\{}-\frac{75}{256}\frac{1}{m_{1}}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\cdot\frac{\bm{S}_{2}}{m_{2}^{2}}\right)\\\
+\frac{603}{3584}\left[\frac{m_{2}}{m_{1}^{2}}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\right)^{2}+\frac{1}{m_{2}}\left(\frac{\bm{S}_{2}}{m_{2}^{2}}\right)^{2}\right]\Bigg{\\}}\hat{\bm{n}}.$
(65)

Exchanging the label $1$ and $2$ in Eq. (65), we obtain the EOM of $2$-nd BH,
governed by the gravitational field of $1$-st BH, $\delta\bm{a}^{2}$. The
relative acceleration $\delta\bm{a}$ is the difference between
$\delta\bm{a}^{1}$ and $\delta\bm{a}^{2}$,

$\delta\bm{a}\equiv\delta\bm{a}^{1}-\delta\bm{a}^{2}=\zeta\frac{\gamma^{3}}{m}\Bigg{\\{}-\frac{75}{256}\frac{1}{\nu}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\cdot\frac{\bm{S}_{2}}{m_{2}^{2}}\right)+\frac{603}{3584}\left[\frac{m^{2}}{m_{1}^{2}}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\right)^{2}+\frac{m^{2}}{m_{2}^{2}}\left(\frac{\bm{S}_{2}}{m_{2}^{2}}\right)^{2}\right]\Bigg{\\}}\hat{\bm{n}},$
(66)

where $\nu\equiv m_{1}m_{2}/m^{2}$ is the symmetric mass ratio. Combined with
the GR part, the total acceleration is written as

$\frac{d\bm{v}}{dt}=-\frac{\gamma^{2}}{m}\left[1+\bar{\mathcal{A}}_{1}+\bar{\mathcal{A}}_{1.5}+(\bar{\mathcal{A}}_{2}+\delta\mathcal{A})\right]\hat{\bm{n}},$
(67)

up to 2PN order. For spin-aligned quasi-circular orbits, the relative
acceleration (66) is completely along the line between the two BHs. The
coefficients $\bar{\mathcal{A}}_{1}$, $\bar{\mathcal{A}}_{1.5}$, and
$\bar{\mathcal{A}}_{2}$, are 1PN, 1.5PN, and 2PN coefficients of the GR part.
These coefficients are Blanchet  _et al._ (1995); Kidder (1995); Faye  _et
al._ (2006); Buonanno  _et al._ (2013)

$\bar{\mathcal{A}}_{1}=\left(1+3\nu\right)v^{2}-2\left(2+\nu\right)\gamma,$
(68) $\bar{\mathcal{A}}_{1.5}=-\gamma\left(5\frac{S}{m^{2}}+3\frac{\delta
m}{m}\frac{\Sigma}{m^{2}}\right)v,$ (69)

and

$\bar{\mathcal{A}}_{2}=\nu(3-4\nu)v^{4}-\gamma\nu\left(\frac{13}{2}-2\nu\right)v^{2}+\gamma^{2}\left(9+\frac{87}{4}\nu\right)+3\gamma^{2}\nu\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\cdot\frac{\bm{S}_{2}}{m_{2}^{2}}\right)+\frac{3}{2}\gamma^{2}\left[\frac{m_{1}^{2}}{m^{2}}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\right)^{2}+\frac{m_{2}^{2}}{m^{2}}\left(\frac{\bm{S}_{2}}{m_{2}^{2}}\right)^{2}\right].$
(70)

The total spin is $\bm{S}\equiv\bm{S}_{1}+\bm{S}_{2}$ and $S\equiv|\bm{S}|$.
The mass difference is $\delta m\equiv m_{1}-m_{2}$ and the spin difference is
$\bm{\Sigma}\equiv m(\bm{S}_{2}/m_{2}-\bm{S}_{1}/m_{1})$, with
$\Sigma\equiv|\bm{\Sigma}|$. As we have mentioned, the SO coupling is included
in the 1.5PN correction (69), and the SS and MQ effects are included in the
2PN correction (70). The dCS modification

$\displaystyle\delta\mathcal{A}$
$\displaystyle=\gamma^{2}\left\\{\frac{75}{256}\nu\cdot\zeta_{12}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\cdot\frac{\bm{S}_{2}}{m_{2}^{2}}\right)-\frac{603}{3584}\left[\zeta_{1}\frac{m_{1}^{2}}{m^{2}}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\right)^{2}+\zeta_{2}\frac{m_{2}^{2}}{m^{2}}\left(\frac{\bm{S}_{2}}{m_{2}^{2}}\right)^{2}\right]\right\\}$
(71)
$\displaystyle=\zeta\gamma^{2}\left\\{\frac{75}{256}\frac{1}{\nu}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\cdot\frac{\bm{S}_{2}}{m_{2}^{2}}\right)-\frac{603}{3584}\left[\frac{m^{2}}{m_{1}^{2}}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\right)^{2}+\frac{m^{2}}{m_{2}^{2}}\left(\frac{\bm{S}_{2}}{m_{2}^{2}}\right)^{2}\right]\right\\}$

is also 2PN approximation. As shown in above equations, the dCS modification
$\delta\mathcal{A}$ (71) has similar dependence on spin with the quadratic-
spin terms of $\bar{\mathcal{A}}_{2}$ (70). To obtain the dCS modification, we
just need to take the replacements
$\bm{S}_{1}\cdot\bm{S}_{2}\rightarrow(25/256)\zeta_{12}(\bm{S}_{1}\cdot\bm{S}_{2})$
for DD coupling and
$\bm{S}_{A}^{2}\rightarrow-(201/1792)\zeta_{A}\bm{S}_{A}^{2}$ for MQ coupling.
The DD effect of dCS modification in Eq. (71) may provide attraction or
repulsion. The scalar fields of two BHs with parallel (antiparallel) SAM will
attract (repulse) each other, just like two parallel (antiparallel) magnetic
moments, and provide negative (positive) potential energy. The modified MQ
effect in Eq. (71) generates repulsive force rather than gravity because
$\delta C_{Q}<0$, as shown in Eq. (51).

###  III.6 Modified Kepler’s Third Law

In Newtonian gravity, the orbital frequency $\omega$ is given by EOM (67),
i.e., $|a|=\omega^{2}r$. Following the standard textbook, we define
dimensionless frequency $x\equiv(m\omega)^{2/3}$. The modified Kepler’s third
law is the relationship between dimensionless frequency and distance, $x$ and
$\gamma$, which are given by

$x=\gamma\left[1+\frac{1}{3}\gamma\bar{\varpi}_{1}+\frac{1}{3}\gamma^{1.5}\bar{\varpi}_{1.5}+\frac{1}{3}\gamma^{2}\left(-\frac{1}{3}\bar{\varpi}_{1}^{2}+\bar{\varpi}_{2}+\delta\varpi\right)\right],$
(72)
$\gamma=x\left[1-\frac{1}{3}x\bar{\varpi}_{1}-\frac{1}{3}x^{1.5}\bar{\varpi}_{1.5}+\frac{1}{3}x^{2}\left(\bar{\varpi}_{1}^{2}-\bar{\varpi}_{2}-\delta\varpi\right)\right].$
(73)

A set of higher PN coefficients are

$\bar{\varpi}_{1}=-3+\nu,$ (74)
$\bar{\varpi}_{1.5}=-5\frac{S}{m^{2}}-3\frac{\delta
m}{m}\frac{\Sigma}{m^{2}},$ (75)
$\bar{\varpi}_{2}=6+\frac{41}{4}\nu+\nu^{2}+3\nu\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\cdot\frac{\bm{S}_{2}}{m_{2}^{2}}\right)+\frac{3}{2}\left[\frac{m_{1}^{2}}{m^{2}}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\right)^{2}+\frac{m_{2}^{2}}{m^{2}}\left(\frac{\bm{S}_{2}}{m_{2}^{2}}\right)^{2}\right],$
(76)

and

$\delta\varpi=\zeta\Bigg{\\{}\frac{75}{256}\frac{1}{\nu}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\cdot\frac{\bm{S}_{2}}{m_{2}^{2}}\right)-\frac{603}{3584}\left[\frac{m^{2}}{m_{1}^{2}}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\right)^{2}+\frac{m^{2}}{m_{2}^{2}}\left(\frac{\bm{S}_{2}}{m_{2}^{2}}\right)^{2}\right]\Bigg{\\}}.$
(77)

Combined with frequency evolution, the modified Kepler’s third law can be used
to calculate the parameterized post-Einsteinian (ppE) parameters without
calculating the full waveform. In terms of $\gamma$ and Eqs. (74)-(77), the
relative EOM of the binary system (67) is rewritten as

$\frac{d\bm{v}}{dt}=-\frac{\gamma^{2}}{m}\left[1+\gamma\bar{\varpi}_{1}+\gamma^{1.5}\bar{\varpi}_{1.5}+\gamma^{2}(\bar{\varpi}_{2}+\delta\varpi)\right]\hat{\bm{n}},$
(78)

up to the 2PN order. The EOM is also can re-expressed in terms of symbol
dimensionless frequency $x$.

###  III.7 Conserved Energy

blackFrom the EOM (78), the conserved energy of the binary system can be given
by “guess-work” Faye  _et al._ (2006); Poisson and Will (2014). We know that
at the Newtonian limit, the conserved energy is $E=(1/2)\mu v^{2}-\mu\gamma$.
All the possible contributions of conserved energy from dCS modification must
be proportional to $\bm{S}_{1}\cdot\bm{S}_{2}$,
$(\hat{\bm{n}}\cdot\bm{S}_{1})(\hat{\bm{n}}\cdot\bm{S}_{2})$, $S_{1}^{2}$,
$S_{2}^{2}$, $(\hat{\bm{n}}\cdot\bm{S}_{1})^{2}$, and
$(\hat{\bm{n}}\cdot\bm{S}_{2})^{2}$. Then the correct combination of such
terms can be identified by including them all (with unknown coefficients) in a
trial expression for $E$ and demanding that $dE/dt=0$ by virtue of the full
2PN EOM. The binding energy of the BBH system with spin-aligned quasi-circular
orbits is given by

$E(x)=-\frac{\mu
x}{2}\left[1+x\bar{\varepsilon}_{1}+x^{1.5}\bar{\varepsilon}_{1.5}+x^{2}(\bar{\varepsilon}_{2}+\delta\varepsilon)\right],$
(79)

with reduced mass being defined as $\mu\equiv m_{1}m_{2}/m=\nu m$. The higher-
order coefficients in Eq. (79) are

$\bar{\varepsilon}_{1}=-\frac{3}{4}-\frac{1}{12}\nu,$ (80)
$\bar{\varepsilon}_{1.5}=\frac{14}{3}\frac{S}{m^{2}}+2\cdot\frac{\delta
m}{m}\frac{\Sigma}{m^{2}},$ (81)
$\bar{\varepsilon}_{2}=-\frac{27}{8}+\frac{19}{8}\nu-\frac{1}{24}\nu^{2}-2\nu\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\cdot\frac{\bm{S}_{2}}{m_{2}^{2}}\right)-\left[\frac{m_{1}^{2}}{m^{2}}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\right)^{2}+\frac{m_{2}^{2}}{m^{2}}\left(\frac{\bm{S}_{2}}{m_{2}^{2}}\right)^{2}\right],$
(82)

and

$\delta\varepsilon=-\zeta\left\\{\frac{25}{128}\frac{1}{\nu}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\cdot\frac{\bm{S}_{2}}{m_{2}^{2}}\right)-\frac{201}{1792}\left[\frac{m^{2}}{m_{1}^{2}}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\right)^{2}+\frac{m^{2}}{m_{2}^{2}}\left(\frac{\bm{S}_{2}}{m_{2}^{2}}\right)^{2}\right]\right\\}.$
(83)

As we have derived, the DD interaction potential is positive (negative) for
antiparallel (parallel) SAM binary BHs and the MQ potential is positive for
both parallel and antiparallel cases. The conserved energy provides the energy
of gravitational radiation. This result is consistent with that in Refs. Yagi
_et al._ (2012b, 2016), in which the scalar fields are seen as magnetic
dipoles, and the modified Dixon tensor is seen as the effective density of BH.
In the simplest case, where the $1$-st and $2$-nd BHs have roughly the same
mass and SAM, the ratio of binding energy caused by the SS and MQ effect is
thus about $0.87$. In conclusion, the correction of SS and MQ effect on
binding energy is comparable.

blackThe acceleration (71), modified Kepler’s third law (72) and binding
energy (79), are consistent with that in Refs. Yagi  _et al._ (2012b); Tahura
and Yagi (2018), although completely different methods are adopted. Our method
is based on the modified MPD equations in dCS theory and guess-work provided
by Refs. Poisson and Will (2014); Faye  _et al._ (2006), which generally
describe the motion and precession of a spinning particle. In Ref. Yagi  _et
al._ (2012b), the modified Newtonian potential satisfies the Poisson equation,
$\bm{\nabla}^{2}\delta U=4\pi(\delta\rho_{1}+\delta\rho_{2})$, with
$\delta\rho_{A}$ being the effective density induced by the dCS modification
of the $A$-th BH. Then, the potential energy density from modified MQ coupling
is determined by $\delta\rho_{1}U_{2}+\delta\rho_{2}U_{1}+\rho_{1}\delta
U_{2}+\rho_{2}\delta U_{1}$, with $\rho_{A}$ being the mass density of the
$A$-th BH. At the same time, the binding energy between the scalar fields of
each BH is calculated by regarding them as two magnetic dipoles. The
interaction force is defined as the derivative of effective potential energy
with respect to the distance between BBH. The approach in Ref. Yagi  _et al._
(2012b) is reasonable up to the lowest-order PN correction. But our method is
more general, especially for the higher accuracy template calculation in the
future.

##  IV Radiation Field

###  IV.1 Scalar Radiation

The scalar radiation is obtained from the scalar perturbation equation (32),
$\Box^{2}_{\eta}\vartheta=-16\pi\sigma+4\pi\rho_{\vartheta}$, where
$\sigma\equiv\sigma_{1}+\sigma_{2}$ and $\rho_{\vartheta}$ are the non-compact
and compact sources of the scalar field. In terms of the near-zone solution of
the scalar field (37) and the PN metric (40), the source terms $\sigma_{1}$
and $\sigma_{2}$ are written as

$16\pi\sigma_{1}=-2\bar{U}(\partial_{i}\partial_{i}\vartheta^{\rm(B)})\sim\mathcal{O}(1/c^{6}),\quad\text{and}\quad
16\pi\sigma_{2}=\frac{8\alpha}{\beta}\epsilon^{ijk}(\partial_{i}\partial_{m}\bar{U})(\partial_{j}\partial_{m}\bar{V}_{k})\sim\mathcal{O}(1/c^{9}).$
(84)

respectively. The monopole, dipole, and quadrupole radiation Maggiore (2008)
generated by non-compact sources are

$\displaystyle\vartheta^{(\sigma)}_{\rm mon}$
$\displaystyle=\frac{4}{R}\int\sigma(t_{r},\mathbf{x}^{\prime})d^{3}\mathbf{x}^{\prime},$
(85a) $\displaystyle\vartheta^{(\sigma)}_{\rm dip}$
$\displaystyle=\frac{4}{R}\frac{\partial}{\partial
t}\int\sigma(t_{r},\mathbf{x}^{\prime})(\hat{\mathbf{N}}\cdot\mathbf{x}^{\prime})d^{3}\mathbf{x}^{\prime},$
(85b) $\displaystyle\vartheta^{(\sigma)}_{\rm quad}$
$\displaystyle=\frac{4}{R}\frac{1}{2}\frac{\partial^{2}}{\partial
t^{2}}\int\sigma(t_{r},\mathbf{x}^{\prime})(\hat{\mathbf{N}}\cdot\mathbf{x}^{\prime})^{2}d^{3}\mathbf{x}^{\prime}.$
(85c)

In Eq. (85), $R$ is the distance between the observer and the GW source system
111Please distinguish the distance $R$ and the Ricci scalar $R$ according to
the situation., $t_{r}\equiv t-R$ is retarded time, and $\hat{\mathbf{N}}$ is
the GW propagation direction. The source $\sigma$ (84) is the product of two
near-zone solutions, e.g., $\bar{U}\vartheta^{\rm(B)}$ and
$\bar{U}\bar{V}_{k}$, etc. Every solution is written in the form, e.g.,
$\bar{U}=\bar{U}_{1}+\bar{U}_{2}$,
$\bar{V}^{k}=\bar{V}_{1}^{k}+\bar{V}_{2}^{k}$, etc. Thus, the retarded
integration consists of the self-interaction terms like
$\bar{U}_{1}\bar{V}_{1}^{k}$ and cross-interaction terms like
$\bar{U}_{1}\bar{V}_{2}^{k}$ (85). The full volume integration of
$1/r_{1}^{2}$ is divergent, without physical meaning. The subsequent
calculation will ignore the self-interaction term. The integration of cross-
interaction terms like $1/(r_{1}r_{2})$ can be managed appropriately by
Hadamard regularization Yagi  _et al._ (2012a); Blanchet (2014),

$\int\frac{d^{3}\mathbf{x}}{r_{1}r_{2}}=-2\pi r.$ (86)

The dominant contribution comes from the monopole term (85a), given by

$\vartheta^{(\sigma)}_{\rm
mon}=-\frac{2}{R}\frac{m}{r^{2}}\left[\frac{m_{2}}{m}(\hat{\bm{n}}\cdot\bm{\mu}_{1})-\frac{m_{1}}{m}(\hat{\bm{n}}\cdot\bm{\mu}_{2})\right].$
(87)

The dipole (85b) and quadrupole radiation (85c) enter higher-order PN
approximation Yagi  _et al._ (2012a, 2016). The estimations are shown in
Appendix B.

Similarly, the multipole moment formula gives the radiation by the source with
compact support,

$\displaystyle\vartheta^{(\rho)}_{\rm mon}$
$\displaystyle=-\frac{1}{R}\int\rho_{\vartheta}(t_{r},\mathbf{x}^{\prime})d^{3}\mathbf{x}^{\prime},$
(88a) $\displaystyle\vartheta^{(\rho)}_{\rm dip}$
$\displaystyle=-\frac{1}{R}\frac{\partial}{\partial
t}\int\rho_{\vartheta}(t_{r},\mathbf{x}^{\prime})(\hat{\mathbf{N}}\cdot\mathbf{x}^{\prime})d^{3}\mathbf{x}^{\prime},$
(88b) $\displaystyle\vartheta^{(\rho)}_{\rm quad}$
$\displaystyle=-\frac{1}{R}\frac{1}{2}\frac{\partial^{2}}{\partial
t^{2}}\int\rho_{\vartheta}(t_{r},\mathbf{x}^{\prime})(\hat{\mathbf{N}}\cdot\mathbf{x}^{\prime})^{2}d^{3}\mathbf{x}^{\prime}.$
(88c)

After calculating Eq. (88), the dominant terms are

$\vartheta^{(\rho)}_{\rm
mon}=\frac{1}{R}\frac{m}{r^{2}}\left[\frac{m_{2}}{m}(\hat{\bm{n}}\cdot\bm{\mu}_{1})-\frac{m_{1}}{m}(\hat{\bm{n}}\cdot\bm{\mu}_{2})\right],$
(89)

and

$\vartheta^{(\rho)}_{\rm
quad}=-\frac{1}{R}\frac{m}{r^{2}}\left[\frac{m_{2}}{m}(\hat{\mathbf{N}}\cdot\bm{\mu}_{1})-\frac{m_{1}}{m}(\hat{\mathbf{N}}\cdot\bm{\mu}_{2})\right](\hat{\mathbf{N}}\cdot\hat{\bm{n}}).$
(90)

The dipole radiation (88b) enters higher PN order. The estimations are also
listed in Appendix B.

As shown in Eqs. (35, 84), the magnitudes of the scalar source terms are
$\sigma_{1}\sim\mathcal{O}(1/c^{6})$, $\sigma_{2}\sim\mathcal{O}(1/c^{9})$,
$\rho_{1}\sim\mathcal{O}(1/c^{4})$, $\rho_{2}\sim\mathcal{O}(1/c^{6})$, and
$\rho_{3}\sim\mathcal{O}(1/c^{6})$, respectively. Generally, the dominant
contribution is sourced by $\rho_{1}$. However, the monopole and dipole
radiations of $\rho_{1}$ vanish because the spin-aligned assumption and that
dCS theory requires the scalar field to be dipole (see Eqs. (185a, 185b)) Yagi
_et al._ (2012c). Such that, the monopole radiation of $\sigma_{1}$, the
quadrupole radiation of $\rho_{1}$, and the monopole radiation of $\rho_{2}$
are of order $1/c^{6}$, while the monopole radiation of $\rho_{3}$ vanishes
due to dipole scalar field. These terms contribute to the dominant scalar
radiation, which are $1/c^{2}$ higher than leading-order gravitational
radiation. The leading-order dCS scalar radiation enters 2PN modification in
energy flux relative to GW radiation (of order $\mathcal{O}(1/c^{4})$).

Let’s make a comparison between dCS theory and BD theory. In the BD gravity,
the scalar field $\varphi$ satisfies equation
$\Box_{\eta}^{2}\varphi=-16\pi\rho(1-2s)/(4\omega_{0}+6)$, with $s$ being the
sensitivity of objects and $\omega_{0}$ being the coupling constant. The
matter energy density is of order $1/c^{2}$. The monopole scalar radiation
vanishes as shown in Zhang  _et al._ (2017a); Liu  _et al._ (2020). The
leading-order (dipole) scalar radiation is of order $1/c^{3}$, which is $1/c$
lower than leading-order gravitational radiation. Finally, the leading-order
BD scalar radiation enters -1PN modification in energy flux relative to GW
radiation.

Combined Eqs. (87), (89), and (90), the total quadrupole scalar radiation is
written as

$\vartheta=\frac{5}{8}\frac{\mu}{R}\frac{\alpha}{\beta
m^{2}}\frac{\gamma^{2}}{\nu}(\hat{\mathbf{N}}\cdot\bm{\Delta})(\hat{\mathbf{N}}\cdot\hat{\bm{n}})\sim\mathcal{O}(1/c^{6}),\quad\text{with}\quad\bm{\Delta}\equiv\frac{m_{2}}{m}\frac{\bm{S}_{1}}{m_{1}^{2}}-\frac{m_{1}}{m}\frac{\bm{S}_{2}}{m_{2}^{2}}.$
(91)

This result (91) will be used to calculate the energy flux carried by
gravitational radiation. Comparing Eq. (91) with the previous results given in
Refs. Yagi  _et al._ (2012a, 2016), some terms are omitted because of the
spin-aligned assumption.

###  IV.2 Tensor Radiation

In the PN framework, the gravitational waveform contains an “instantaneous”
term, depending on the state of the binary at the retarded time $t_{r}$ only,
and a “tail” term, which is sensitive to the wave field at all previous time
$t_{r}-\tau$. We write the waveform as two pieces,

$\bar{H}_{ij}=(\bar{H}_{ij})_{\rm inst}+(\bar{H}_{ij})_{\rm tail}.$ (92)

The instantaneous and tail terms are given by Blanchet  _et al._ (1995);
Blanchet (2014)

$\displaystyle(\bar{H}_{ij})_{\rm inst}$
$\displaystyle=\frac{2}{R}\Bigg{\\{}\bar{\mathcal{I}}^{(2)}_{ij}+\left[\frac{1}{3}\bar{\mathcal{I}}^{(3)}_{ija}\hat{N}^{a}+\frac{4}{3}\epsilon_{ab(i}\bar{\mathcal{J}}^{(2)}_{j)a}\hat{N}^{b}\right]+\left[\frac{1}{12}\bar{\mathcal{I}}^{(4)}_{ijab}\hat{N}^{ab}+\frac{1}{2}\epsilon_{ab(i}\bar{\mathcal{J}}^{(3)}_{j)ac}\hat{N}^{bc}\right]$
(93)
$\displaystyle+\left[\frac{1}{60}\bar{\mathcal{I}}^{(5)}_{ijabc}\hat{N}^{abc}+\frac{2}{15}\epsilon_{ab(i}\bar{\mathcal{J}}_{j)acd}\hat{N}^{bcd}\right]+\left[\frac{1}{360}\bar{\mathcal{I}}^{(6)}_{ijabcd}\hat{N}^{abcd}+\frac{1}{36}\epsilon_{ab(i}\bar{\mathcal{J}}^{(5)}_{j)acde}\hat{N}^{bcde}\right]+\delta{\mathcal{I}}^{(2)}_{ij}\Bigg{\\}},$

and

$\displaystyle(\bar{H}_{ij})_{\rm
tail}=\frac{4m}{R}\int_{0}^{+\infty}d\tau\left[\ln\left(\frac{\tau}{\tau_{1}}\right)\left(\bar{\mathcal{I}}_{ij}+\delta\mathcal{I}_{ij}\right)^{(4)}+\frac{1}{3}\ln\left(\frac{\tau}{\tau_{2}}\right)\hat{N}_{a}\bar{\mathcal{I}}_{ija}^{(5)}+\frac{4}{3}\ln\left(\frac{\tau}{\tau_{3}}\right)\epsilon_{ab(i}\hat{N}_{b}\bar{\mathcal{J}}_{j)a}^{(4)}\right](t_{r}-\tau),$
(94)

with $\hat{N}^{i_{1}i_{2}\cdots
i_{l}}\equiv\hat{N}^{i_{1}}\hat{N}^{i_{2}}\cdots\hat{N}^{i_{l}}$.
$\bar{\mathcal{I}}_{L}$ and $\bar{\mathcal{J}}_{iL}$ are mass and current
multipole moments in GR. $\tau_{1},\tau_{2},\tau_{3}$ are arbitrary constants
Blanchet  _et al._ (1995). All of these moments can be divided into mass (M),
spin (S), and quadratic-spin (Q) parts,

$\bar{\mathcal{I}}_{L}=\bar{\mathcal{I}}_{L}^{\rm(M)}+\bar{\mathcal{I}}_{L}^{\rm(S)}+\bar{\mathcal{I}}_{L}^{\rm(Q)},\quad\bar{\mathcal{J}}_{iL}=\bar{\mathcal{J}}_{iL}^{\rm(M)}+\bar{\mathcal{J}}_{iL}^{\rm(S)}+\bar{\mathcal{J}}_{iL}^{\rm(Q)}.$
(95)

For example, the mass quadrupole moment in GR is Blanchet  _et al._ (1995);
Kidder (1995); Faye  _et al._ (2006); Bohé  _et al._ (2015); Buonanno  _et
al._ (2013)

$\displaystyle\bar{\mathcal{I}}_{ij}^{\rm(M)}$ $\displaystyle=\mu
r^{2}\Bigg{\\{}\left[1+\gamma\left(-\frac{1}{42}-\frac{13}{14}\nu\right)+\gamma^{2}\left(-\frac{461}{1512}-\frac{18395}{1512}\nu-\frac{241}{1512}\nu^{2}\right)\right](\hat{n}^{i}\hat{n}^{j})^{\rm
STF}$ (96)
$\displaystyle\qquad+\left[\left(\frac{11}{21}-\frac{11}{7}\nu\right)+\gamma\left(\frac{1607}{378}-\frac{1681}{378}\nu+\frac{229}{378}\nu^{2}\right)\right](v^{i}v^{j})^{\rm
STF}\Bigg{\\}},$ $\bar{\mathcal{I}}_{ij}^{\rm(S)}=\nu
r\left\\{\frac{8}{3}[\hat{n}_{i}(\bm{v}\times\bm{S})_{j}]^{\rm
STF}-\frac{4}{3}[v_{i}(\hat{\bm{n}}\times\bm{S})_{j}]^{\rm
STF}+\frac{8}{3}\frac{\delta
m}{m}[\hat{n}_{i}(\bm{v}\times\bm{\Sigma})_{j}]^{\rm
STF}-\frac{4}{3}\frac{\delta
m}{m}[v_{i}(\hat{\bm{n}}\times\bm{\Sigma})_{j}]^{\rm STF}\right\\},$ (97)

and

$\bar{\mathcal{I}}_{ij}^{\rm(Q)}=-\left[\frac{1}{m_{1}}(S_{1}^{i}S_{1}^{j})^{\rm
STF}+\frac{1}{m_{2}}(S_{2}^{i}S_{2}^{j})^{\rm STF}\right].$ (98)

STF means symmetric tracefree tensor, e.g., $(\hat{n}_{i}\hat{n}_{j})^{\rm
STF}=\hat{n}_{i}\hat{n}_{j}-(1/3)\delta_{ij}$ Maggiore (2008). The higher-
order mass and current moments can be found in Refs. Blanchet  _et al._ (1995,
2006); Buonanno  _et al._ (2013); Bohé  _et al._ (2015). For spin-aligned
quasi-circular orbits, the quadratic-spin part (98) is conserved. Thus they
will not contribute to gravitational radiation. Similarly, the dCS
modification of mass quadrupole moment,

$\delta\mathcal{I}_{ij}=\frac{201}{1792}\zeta\left[\frac{m^{4}}{m_{1}^{5}}(S_{1}^{i}S_{1}^{j})^{\rm
STF}+\frac{m^{4}}{m_{2}^{5}}(S_{2}^{i}S_{2}^{j})^{\rm STF}\right],$ (99)

is also conserved, i.e., $\delta\bar{\mathcal{I}}_{ij}^{(2)}=0$ and
$\delta\bar{\mathcal{I}}_{ij}^{(4)}=0$. The dCS modification to the fourth-
order time derivative of mass quadrupole moment enters higher order. Thus the
tail term (94) has no dCS modification. The dCS modification only comes from
the second-order time derivative of the mass quadrupole moment. The total
gravitational radiation waveform is Blanchet  _et al._ (1995); Arun  _et al._
(2009)

$\displaystyle\bar{H}_{ij}$
$\displaystyle=\frac{2\mu}{R}\cdot\Big{[}2(v_{i}v_{j}-\gamma\hat{n}_{i}\hat{n}_{j})+\\{\text{Higher-
Order PN Waveforms in GR}\\}\Big{]}+\bar{k}_{ij}.$ (100)

We list the leading-order and dCS waveform in Eq. (100) for simplification.
The 0.5PN, $\cdots$, 2PN waveform in GR can be found in Refs. Blanchet  _et
al._ (1995); Kidder (1995); Buonanno  _et al._ (2013); Bohé  _et al._ (2015);
Blanchet  _et al._ (2006), including the SO effect at 1.5PN order, SS and MQ
effects at 2PN order.

The metric-deformation radiation $\bar{k}_{ij}$ is explicitly expressed as

$\bar{k}_{ij}=-\frac{4\mu\gamma^{3}}{R}\delta\varpi\cdot\hat{n}^{i}\hat{n}^{j}\approx-\frac{4\mu
x^{3}}{R}\delta\varpi\cdot\hat{n}^{i}\hat{n}^{j}.$ (101)

The coefficient $\delta\varpi$ is proportional to the dimensionless coupling
constant $\zeta$, containing $\bm{S}_{1}^{2}$ -, $\bm{S}_{2}^{2}$ -, and
$\bm{S}_{1}\cdot\bm{S}_{2}$ -terms, defined in Eqs. (72) and (77). The
waveform (101) is similar to the quadratic-spin waveform in GR except for a
coupling constant.

The complete metric-deformation waveform is not presented in Ref. Yagi  _et
al._ (2012a) because some source terms in the tensor perturbation equation are
missed. The erratum Yagi  _et al._ (2016) corrected the scalar and tensor
perturbation equation, including all possible source terms, concluded that the
energy flux carried by metric-deformation radiation is at 3PN order, far less
than that of the scalar radiation (2PN). Such that the metric-deformation
waveform is not given. However, we find that the scalar and tensor flux are
both at 2PN order. The reasons are as follows. The lowest-order radiation is
given by retarded-potential integral of the source and can be rewritten via
conservation law, $\partial_{\mu}K^{\mu\nu}=0$,

$\bar{k}_{ij}\approx\frac{4}{R}\int
K^{ij}d^{3}\mathbf{x}^{\prime}=\frac{2}{R}\left[2\int\partial_{k}(x^{\prime
i}K^{kj})d^{3}\mathbf{x}^{\prime}+\partial_{0}\int\partial_{k}\left(x^{\prime
i}x^{\prime
j}K^{0k}\right)d^{3}\mathbf{x}^{\prime}+\partial_{0}\partial_{0}\int x^{\prime
i}x^{\prime j}K^{00}d^{3}\mathbf{x}^{\prime}\right].$ (102)

The source $K_{ij}$ contains compact-support terms (the first two terms in Eq.
(22)) and non-compact terms (the last three terms in Eq. (22)). The compact
sources tend to $0$, strictly, while the non-compact sources tend to $1/R^{2}$
at far regions. For the compact case, the two complete-derivative integrations
in the square bracket of the Eq. (102) can be transformed into two infinity-
sphere integrations and then vanish. However, for the non-compact case, these
two volume integrations should be calculated as follow.

$\displaystyle\int\partial_{k}(x^{\prime i}K^{kj})d^{3}\mathbf{x}^{\prime}$
$\displaystyle=\oint(x^{\prime i}K^{kj})dS_{k}=\oint\hat{N}^{k}(x^{\prime
i}K^{kj})R^{2}d\Omega\propto R\frac{1}{R^{4}}R^{2}=\frac{1}{R}\rightarrow 0,$
(103) $\displaystyle\int\partial_{k}(x^{\prime i}x^{\prime
j}K^{0k})d^{3}\mathbf{x}^{\prime}$ $\displaystyle=\oint(x^{\prime i}x^{\prime
j}K^{0k})dS_{k}=\oint\hat{N}^{k}(x^{\prime i}x^{\prime
j}K^{0k})R^{2}d\Omega\propto R^{2}\frac{1}{R^{4}}R^{2}\rightarrow\text{finite
value}.$

According to the above discussion about gravitational radiation, we conclude
that

$\frac{4}{R}\int
K^{ij}d^{3}\mathbf{x}^{\prime}\neq\frac{2}{R}\partial_{0}\partial_{0}\int
x^{\prime i}x^{\prime j}K^{00}d^{3}\mathbf{x}^{\prime}.$ (104)

Such that, the estimations to the PN order in the Eqs. (11, 16, 17) in Ref.
Yagi  _et al._ (2016) are correct, which are 2PN corrections, but Eqs. (19-21)
(in Ref. Yagi  _et al._ (2016)) are incorrect, because they assumed that the
left-hand side equals to the right-hand side of Eq. (104) for the non-compact
source. Therefore, we recalculate the full metric-deformation waveform in this
section for spin-aligned case, the standard PN method Blanchet (2014) is
adopted. As shown in our result (101), the gravitational radiation contributes
to the energy flux at 2PN order rather than 3PN as predicted by Ref. Yagi  _et
al._ (2016). Accordingly, the calculation of total radiant energy flux and ppE
parameters in Refs. Yagi  _et al._ (2016, 2012b) are also inaccurate.

Moreover, another equivalent derivation is shown in Appendix C. The radiation
field of the metric deformation is calculated through the monopole formula of
$K_{ij}$ rather than the quadrupole formula of $K_{00}$ because of Eq. (104).
Thus, the EOM of the BBH is not required for calculating the radiation in the
wave zone. On the other hand, the time derivative of $K_{ij}$ is of higher PN
order than we required, so it will not be included during calculation by the
multipole moment formulas. However, the time derivative of the lowest-order
moment will bring the modification of the EOM of the BBH. In the required
order, the two methods are equivalent. Because the conservation of the matter
EMT and the gravitational field determines the motion of the matter field
Roshan (2013), The influence caused by the modified EOM on the lowest-order
moment is equivalent to that caused by the modified source term on the higher-
order moments.

###  IV.3 Propagation Frame

Following Ref. Maggiore (2008), we set the propagation direction of GWs as
$\hat{\mathbf{N}}=(0,\sin\iota,\cos\iota)$, with the inclination angle $\iota$
between the orbital plane and propagation direction. The relative position and
velocity of BBH are $\hat{\bm{n}}=(\cos\phi,\sin\phi,0)$ and $\bm{v}=\omega
r(-\sin\phi,\cos\phi,0)$, where $\phi$ is the orbital phase evaluated at
retarded time, defined as $\phi(t)\equiv\omega t_{r}=\omega(t-R)$. We then
transform the metric deformation $\bar{k}_{ij}$ (101) from binary frame to
propagation frame, in which propagation direction is regarded as the $z$-axis.
Through the rotation matrix,

$\bm{\mathcal{R}}(\hat{\mathbf{N}})=\left(\begin{array}[]{ccc}1&0&0\\\
0&\cos\iota&-\sin\iota\\\ 0&\sin\iota&\cos\iota\end{array}\right),$ (105)

the metric deformation $\bar{k}_{ij}$ (101) becomes

$\bar{k}_{ij}(\hat{\mathbf{N}})=-\frac{2\mu
x^{3}}{R}\delta\varpi\left(\begin{array}[]{ccc}\cos 2\phi&\cos\iota\sin
2\phi&\sin\iota\sin 2\phi\\\ \cos\iota\sin 2\phi&-\cos^{2}\iota\cos
2\phi&-\sin\iota\cos\iota\cos 2\phi\\\ \sin\iota\sin
2\phi&-\sin\iota\cos\iota\cos 2\phi&-\sin^{2}\iota\cos
2\phi\end{array}\right).$ (106)

Several constants have been omitted because they do not contribute to energy
flux carried by radiation. This result (106) will be used to study the
polarization modes and the energy flux carried by the tensor radiation.

##  V Polarization Modes

The polarization modes of GWs from BBH system with spin-aligned quasi-circular
orbit in dCS gravity are studied through NP formalism Newman and Penrose
(1962). The NP null tetrads are reviewed, and the polarization modes of metric
deformation are obtained.

###  V.1 Null Tetrad

The NP null tetrad is $\bm{e}_{(a)}\equiv(\bm{l},\bm{q},\bm{m},\bar{\bm{m}})$
Chandrasekhar (1983); Wagle  _et al._ (2019), where $\bar{\bm{m}}$ is the
complex conjugation of vector, $\bm{m}$ and the definitions of $\bm{l}$,
$\bm{q}$, and $\bm{m}$ are

$\bm{l}=\frac{1}{\sqrt{2}}(1,0,0,1),\quad\bm{q}=\frac{1}{\sqrt{2}}(1,0,0,-1),\quad\text{and}\quad\bm{m}=\frac{1}{\sqrt{2}}(0,1,i,0).$
(107)

These basis vectors are lightlike and satisfy the orthogonality,
$\bm{l}\cdot\bm{l}=\bm{q}\cdot\bm{q}=\bm{m}\cdot\bm{m}=\bar{\bm{m}}\cdot\bar{\bm{m}}=0$
and $\bm{l}\cdot\bm{q}=-1$, $\bm{m}\cdot\bar{\bm{m}}=1$. The flat metric
relates the null tetrads by
$\eta^{\mu\nu}=-l^{\mu}q^{\nu}-q^{\mu}l^{\nu}+m^{\mu}\bar{m}^{\nu}+\bar{m}^{\mu}m^{\nu}$.
Projecting this equation onto null tetrad $\bm{e}_{(a)}$, we obtain the flat
metric in NP formalism Chandrasekhar (1983),

$\eta_{(a)(b)}\equiv
e_{(a)}^{\mu}e_{(b)}^{\nu}\eta_{\mu\nu}=\left(\begin{array}[]{cccc}0&-1&0&0\\\
-1&0&0&0\\\ 0&0&0&1\\\ 0&0&1&0\\\ \end{array}\right).$ (108)

$\eta_{\mu\nu}$ and $\eta^{\mu\nu}$ are used to raise and fall the coordinate
indices, and $\eta_{(a)(b)}$ and $\eta^{(a)(b)}$ (108) to raise and fall the
tetrad indices. All geometric quantities can be projected onto null tetrad
$\bm{e}_{(a)}$, e.g., Riemann tensor

$R_{(a)(b)(c)(d)}=R_{\mu\nu\rho\lambda}e_{(a)}^{\mu}e_{(b)}^{\nu}e_{(c)}^{\rho}e_{(d)}^{\lambda}.$
(109)

In NP formalism, $10$ Ricci scalars ($\Phi_{ij}$ ($i,j=0,1,2$) and $\Lambda$)
and $5$ complex Weyl scalars ($\Psi_{0},\cdots,\Psi_{4}$) are introduced to
describe spacetime. Four of them,

$\Psi_{2}=\frac{1}{6}R_{(2)(1)(2)(1)},\quad\Psi_{3}=\frac{1}{2}R_{(2)(1)(2)(4)},\quad\Psi_{4}=R_{(2)(4)(2)(4)},\quad\text{and}\quad\Phi_{22}=R_{(2)(4)(2)(3)},$
(110)

describe all possible polarization modes of gravitational radiation.

On the other hand, the geodesic deviation equation,

$\frac{d^{2}\xi_{i}}{dt^{2}}=-R_{0i0j}\xi^{j},$ (111)

indicates that some components of the Riemann tensor describe the
observational effects of the GWs, where $\xi^{j}$ is the deviation vector. In
order to obtain the polarization modes of GW radiation, it is convenient to
adopt a specific gauge for the spacetime metric. For GR theory, we know that
there are just two DOFs in the GW metric through the gauge covariance, such
that the transverse traceless (TT) gauge is adopted to obtain the polarization
modes. However, for modified gravity, the number of polarization modes may
exceed two, and the TT gauge can not be adopted arbitrarily. The NP quantities
(110) are a set of gauge invariants. They relate to the observational effects
of GWs, such that NP tetrad becomes a powerful tool to study the GW
polarization modes in modified gravity. The components $R_{0i0j}$ (111) can be
written in terms of Ricci scalar $\Phi_{22}$ and Weyl scalars
$\\{\Psi_{2},\Psi_{3},\Psi_{4}\\}$ (110) as

$\displaystyle R_{0101}$ $\displaystyle=\frac{1}{2}(\Phi_{22}+{\rm
Re}\Psi_{4}),\quad R_{0102}=-\frac{1}{2}{\rm Im}\Psi_{4},\quad R_{0103}=2{\rm
Re}\Psi_{3},$ (112) $\displaystyle R_{0202}$
$\displaystyle=\frac{1}{2}(\Phi_{22}-{\rm Re}\Psi_{4}),\quad R_{0203}=-2{\rm
Im}\Psi_{3},\quad R_{0303}=6\Psi_{2}.$

“Re” and “Im” represent the real and imaginary parts of a quantity. On the
other hand, the polarization modes are defined in terms of the components of
the Riemann tensor (111) in the following form,

$R_{0i0j}=\left(\begin{array}[]{ccc}\ddot{H}_{b}+\ddot{H}_{+}&\ddot{H}_{\times}&\ddot{H}_{x}\\\
\ddot{H}_{\times}&\ddot{H}_{b}-\ddot{H}_{+}&\ddot{H}_{y}\\\
\ddot{H}_{x}&\ddot{H}_{y}&\ddot{H}_{L}\end{array}\right)=\left(\begin{array}[]{ccc}\frac{1}{2}(\Phi_{22}+{\rm
Re}\Psi_{4})&-\frac{1}{2}{\rm Im}\Psi_{4}&2{\rm Re}\Psi_{3}\\\
-\frac{1}{2}{\rm Im}\Psi_{4}&\frac{1}{2}(\Phi_{22}-{\rm Re}\Psi_{4})&-2{\rm
Im}\Psi_{3}\\\ 2{\rm Re}\Psi_{3}&-2{\rm
Im}\Psi_{3}&6\Psi_{2}\end{array}\right).$ (113)

In Eq. (113), $H_{+}$ and $H_{\times}$ are the plus and cross modes predicted
by GR. The other four independent components include two scalar modes,
breathing mode $H_{b}$ and longitudinal mode $H_{L}$, and two vector modes,
$x$-mode $H_{x}$ and $y$-mode $H_{y}$ Will (2014).

###  V.2 Polarization Modes

In subsection IV.3, we have given the spatial components of trace-reverse
metric deformation in the propagation frame, (106). Accordingly, the time
component and time-spatial components are obtained by the Lorenz gauge (16),

$\bar{k}_{0i}=\int
dt(\partial_{j}\bar{k}_{ji})=-\bar{k}_{zi},\quad\text{and}\quad\bar{k}_{00}=\int
dt(\partial_{j}\bar{k}_{j0})=-\bar{k}_{z0}=\bar{k}_{zz}.$ (114)

Thus the complete metric deformation $k_{\mu\nu}$ is

$\displaystyle k_{\mu\nu}(\hat{\mathbf{N}})=-\frac{4\mu x^{3}}{R}\delta\varpi$
(115)
$\displaystyle\times\left(\begin{array}[]{cccc}\frac{1}{2}\sin^{2}\iota\sin^{2}\phi&-\frac{1}{2}\sin\iota\sin
2\phi&-\sin\iota\cos\iota\sin^{2}\phi&-\sin^{2}\iota\sin^{2}\phi\\\
-\frac{1}{2}\sin\iota\sin
2\phi&\frac{1}{2}(\cos^{2}\phi-\cos^{2}\iota\sin^{2}\phi)&\frac{1}{2}\cos\iota\sin
2\phi&\frac{1}{2}\sin\iota\sin 2\phi\\\
-\sin\iota\cos\iota\sin^{2}\phi&\frac{1}{2}\cos\iota\sin
2\phi&\frac{1}{2}(\cos^{2}\iota\sin^{2}\phi-\cos^{2}\phi)&\sin\iota\cos\iota\sin^{2}\phi\\\
-\sin^{2}\iota\sin^{2}\phi&\frac{1}{2}\sin\iota\sin
2\phi&\sin\iota\cos\iota\sin^{2}\phi&\left(1-\frac{3}{2}\cos^{2}\iota\right)\sin^{2}\phi-\frac{1}{2}\cos^{2}\phi\\\
\end{array}\right).$

We denote the dCS part of Ricci and Weyl scalars as $\delta\Psi_{2}$,
$\delta\Psi_{3}$, $\delta\Psi_{4}$, and $\delta\Phi_{22}$. The only non-zero
quantity is

$k_{+}-ik_{\times}=\frac{1}{2}\int dt\int
dt^{\prime}\delta\Psi_{4}=-\frac{2\mu
x^{3}}{R}\delta\varpi\left(\frac{1+\cos^{2}\iota}{2}\cos 2\phi-i\cos\iota\sin
2\phi\right).$ (116)

Other quantities in Eq. (110) vanish, i.e.,
$\delta\Psi_{2}=\delta\Psi_{3}=\delta\Phi_{22}=0$. This result indicates that
there is no extra polarization mode in dCS theory, which is consistent with
the conclusion in Ref. Wagle  _et al._ (2019). Finally, the polarization modes
of tensor radiation are

$\displaystyle H_{+,\times}(t)=h_{+,\times}(t)+k_{+,\times}(t)$
$\displaystyle=\frac{4\mu
x}{R}\Bigg{\\{}h_{+,\times}^{(0)}+x^{1/2}h_{+,\times}^{(1/2)}+x\left[h_{+,\times}^{(1)}+h_{+,\times}^{\rm(1,SO)}\right]$
(117)
$\displaystyle+x^{3/2}\left[h_{+,\times}^{(3/2)}+h_{+,\times}^{\rm(3/2,SO)}\right]+x^{2}\left[h_{+,\times}^{(2)}+h_{+,\times}^{\rm(2,SO)}+h_{+,\times}^{\rm(2,SS)}+h_{+,\times}^{\rm(2,MQ)}+k_{+,\times}^{(2)}\right]\Bigg{\\}}.$

For example, the leading-order coefficients are

$h_{+}^{(0)}=-\frac{1+\cos^{2}\iota}{2}\cos 2\psi(t)\quad\text{and}\quad
h_{\times}^{(0)}=-\cos\iota\sin 2\psi(t),$ (118)

and the subleading-order (0.5PN) coefficients are

$\displaystyle h_{+}^{(1/2)}$ $\displaystyle=-\frac{\sin\iota}{16}\frac{\delta
m}{m}\left[(5+\cos^{2}\iota)\cos\psi(t)-9(1+\cos^{2}\iota)\cos
3\psi(t)\right],$ (119) $\displaystyle h_{\times}^{(1/2)}$
$\displaystyle=-\frac{3}{8}\sin\iota\cos\iota\frac{\delta
m}{m}[\sin\psi(t)-3\sin 3\psi(t)],$

where $\psi(t)$ is a phase related to $\phi(t)$ by

$\psi(t)=\phi(t)-2m\omega\ln\left(\frac{\omega(t)}{\omega_{0}}\right).$ (120)

The constant frequency $\omega_{0}$ can be conveniently chosen as the entry
frequency of an interferometric detector. The use of the $\psi$ instead of the
actual phase $\phi$ of the source is convenient because it allows collecting
the logarithmic terms that come out of the computation of the tail effects
Maggiore (2008).

Additionally, the SS and MQ effects in GR are

$h_{+,\times}^{\rm(2,SS)}=2\nu\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\cdot\frac{\bm{S}_{2}}{m_{2}^{2}}\right)h_{+,\times}^{(0)},\quad\text{and}\quad
h_{+,\times}^{\rm(2,MQ)}=\frac{1}{2}\left[\frac{m_{1}^{2}}{m^{2}}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\right)^{2}+\frac{m_{2}^{2}}{m^{2}}\left(\frac{\bm{S}_{2}}{m_{2}^{2}}\right)^{2}\right]h_{+,\times}^{(0)}.$
(121)

The dCS modification is

$k_{+,\times}^{(2)}=\frac{2}{3}\cdot\delta\varpi\cdot h_{+,\times}^{(0)}.$
(122)

Other coefficients in Eq. (117) can be found in Refs. Blanchet (2014); Kidder
(1995); Buonanno  _et al._ (2013). Notice here that we have replaced $\gamma$
and $v^{2}$ in Eq. (100) with $x$ by the modified Kepler’s third law (73).

From Eqs. (116) and (122), we find that there is no scalar polarization mode
in the dCS gravity, although the dCS theory includes a scalar DOF. The primary
reason is that the dCS theory is a quadratic theory. In the dCS action (1),
the scalar field couples with the square of the Riemann tensor. Therefore, the
first-order scalar field and metric field are independent in the perturbation
equation. Such that the metric deformation does not be affected by the scalar
DOF, and then the breathing mode does not appear. At the same time, the Cotton
tensor (3) is traceless, such that the dCS gravity is a massless theory.
Therefore, the longitudinal mode does not appear as well. Our analysis is
based on the quadratic-order perturbation equation rather than PN waveforms.
Such that, the above conclusion is generic for linear GWs and independent of
the PN orders. Additionally, this statement is consistent with the conclusion
of Ref. Wagle  _et al._ (2019), who adopted the Newman-Penrose analysis, which
is independent of the PN expansion.

As we have derived in Eq. (122), the two polarization modes are generated
symmetrically. Thus there is no parity-violating effect in the modified GW
waveform. This is because we mainly consider the orbital evolution and
gravitational radiation of binary systems in the wave zone and in a short time
scale before merging. This time scale is very short compared with the
expansion time scale of the Universe, so that the homogeneous background
scalar field in the dCS theory does not need to be considered. The amplitude
birefringence, a kind of important parity-violating effect, will be
investigated in the propagation zone (see Section VIII).

##  VI Energy Flux and Orbital Evolution

###  VI.1 Scalar Field

From the EMT of the scalar field (4), the time-spatial component is Yagi  _et
al._ (2012a)

$T_{0i}^{(\vartheta)}=\beta\dot{\vartheta}(\partial_{i}\vartheta),$ (123)

and then the scalar energy flux is

$\mathcal{F}_{S}=-\oint_{\partial\Omega}\left\langle\hat{N}_{i}T_{0i}^{(\vartheta)}\right\rangle
R^{2}d\Omega=\beta
R^{2}\oint_{\partial\Omega}\langle\dot{\vartheta}^{2}\rangle d\Omega,$ (124)

where $\left\langle\cdots\right\rangle$ represents averaging over a period.
Then, substituting the scalar waveform (91) into Eq. (124), using the integral
formula,

$\int\hat{N}_{i}\hat{N}_{j}\hat{N}_{k}\hat{N}_{m}d\Omega=\frac{4\pi}{15}\left(\delta_{ij}\delta_{km}+\delta_{ik}\delta_{jm}+\delta_{im}\delta_{jk}\right),$
(125)

and averaging formula, $\left\langle
v^{2}\right\rangle=\left\langle(\hat{\bm{v}}\cdot\bm{\Delta})^{2}\right\rangle=0$,
we get the final scalar energy flux,

$\mathcal{F}_{S}=\frac{32}{5}\nu^{2}x^{7}\cdot\frac{25}{24576}\frac{\zeta}{\nu^{2}}\bm{\Delta}^{2}.$
(126)

The symbol $\bm{\Delta}$ has been defined in Eq. (91). This result (126) is
consistent with that given by Ref. Yagi  _et al._ (2016), while the term
$27\left\langle(\bm{\Delta}\cdot\bm{v})^{2}\right\rangle$ is omitted because
of spin-aligned assumption.

###  VI.2 Tensor Field

In terms of metric tensor (12), the time-spatial component of the EMT of the
gravitational field is

$T_{0i}^{(H)}=\frac{1}{32\pi}\left\langle(\partial_{0}\bar{H}^{\rm
TT}_{jk})(\partial_{i}\bar{H}^{\rm TT}_{jk})\right\rangle.$ (127)

The radiation energy flux is written as

$\mathcal{F}_{H}=-\oint_{\partial\Omega}\left\langle\hat{N}_{i}T_{0i}^{(H)}\right\rangle
R^{2}d\Omega\approx\frac{R^{2}}{32\pi}\oint_{\partial\Omega}\left[\langle\dot{h}^{\rm
TT}_{jk}\dot{h}^{\rm TT}_{jk}\rangle+2\langle\dot{h}^{\rm TT}_{jk}\dot{k}^{\rm
TT}_{jk}\rangle\right]d\Omega.$ (128)

The first term in the square bracket is the GR part, while the second
represents the contributions from the dCS theory. We rewrite them as

$\bar{\mathcal{F}}_{H}=\frac{1}{16\pi}R^{2}\oint_{\partial\Omega}\left[\langle\dot{h}_{+}^{2}\rangle+\langle\dot{h}_{\times}^{2}\rangle\right]d\Omega\quad\text{and}\quad\delta\mathcal{F}_{H}=\frac{1}{8\pi}R^{2}\oint_{\partial\Omega}\left[\langle\dot{h}_{+}\dot{k}_{+}\rangle+\langle\dot{h}_{\times}\dot{k}_{\times}\rangle\right]d\Omega$
(129)

in terms of polarization modes (100). The GR part has been fully calculated in
previous work Kidder (1995); Blanchet  _et al._ (2002, 2006); Blanchet (2014);
Bohé  _et al._ (2015); Cho  _et al._ (2021, 2022). Finally, we obtain the dCS
modification of the energy flux carried by tensor radiation,

$\delta\mathcal{F}_{H}=\frac{32}{5}\nu^{2}x^{7}\cdot\frac{4}{3}\delta\varpi.$
(130)

We have used
$\langle\sin^{2}[2\phi(t)]\rangle=\langle\cos^{2}[2\phi(t)]\rangle=1/2$.
Combined Eqs. (126) and (130), the total energy flux is written as

$\mathcal{F}=\frac{32}{5}\nu^{2}x^{5}\left[1+x\bar{\mathscr{F}}_{1}+x^{3/2}\bar{\mathscr{F}}_{1.5}+x^{2}(\bar{\mathscr{F}}_{2}+\delta\mathscr{F})\right].$
(131)

The 1PN, 1.5PN, and 2PN coefficients are given by

$\bar{\mathscr{F}}_{1}=-\frac{1247}{336}-\frac{35}{12}\nu,$ (132)
$\bar{\mathscr{F}}_{1.5}=4\pi-\frac{1}{m^{2}}\left[4S+\frac{5}{4}\frac{\delta
m}{m}\Sigma\right],$ (133)

and

$\mathscr{F}_{2}=-\frac{44711}{9072}+\frac{9271}{504}\nu+\frac{65}{18}\nu^{2}+\frac{31}{8}\nu\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\cdot\frac{\bm{S}_{2}}{m_{2}^{2}}\right)+\frac{33}{16}\left[\frac{m_{1}^{2}}{m^{2}}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\right)^{2}+\frac{m_{2}^{2}}{m^{2}}\left(\frac{\bm{S}_{2}}{m_{2}^{2}}\right)^{2}\right],$
(134)

respectively. The factor $4\pi$ in (133) comes from tail terms of GW
radiation. The dCS modification of energy flux is

$\displaystyle\delta\mathscr{F}=\frac{25}{24576}\frac{\zeta}{\nu^{2}}\bm{\Delta}^{2}+\frac{8}{3}\delta\varpi=\zeta\left\\{\frac{9575}{12288}\frac{1}{\nu}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\cdot\frac{\bm{S}_{2}}{m_{2}^{2}}\right)-\frac{77009}{172032}\left[\frac{m^{2}}{m_{1}^{2}}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\right)^{2}+\frac{m^{2}}{m_{2}^{2}}\left(\frac{\bm{S}_{2}}{m_{2}^{2}}\right)^{2}\right]\right\\}.$
(135)

Previous work has also given the dCS flux. However, the scalar flux given by
Ref. Yagi  _et al._ (2012a) is not complete, and the metric deformation flux
is absent in Ref. Yagi  _et al._ (2016) because of unreasonable magnitude
estimation. Similar to the discussion after Eq. (83), for equal-mass and equal
spin systems, the ratio of radiation energy flux caused by SS and MQ effect is
about $0.14$. That is, the correction of the MQ effect on radiant energy flux
and orbital evolution is much more obvious than that of the SS effect.

###  VI.3 Orbital Evolution

The energy flux carried by radiation comes from the binding energy of the
binary system. Using the conservation of energy, $\mathcal{F}=-\dot{E}$, the
differential equation of frequency evolution with time can be established in
terms of the dimensionless frequency $x$,

$\dot{x}=\frac{64}{5}\frac{\nu}{m}x^{5}\left[1+x\bar{\Omega}_{1}+x^{3/2}\bar{\Omega}_{1.5}+x^{2}\left(\bar{\Omega}_{2}+\delta\Omega\right)\right].$
(136)

The PN coefficients are

$\bar{\Omega}_{1}=-\frac{743}{336}-\frac{11}{4}\nu,$ (137)
$\bar{\Omega}_{1.5}=4\pi-\frac{1}{m^{2}}\left[\frac{47}{3}S+\frac{25}{4}\frac{\delta
m}{m}\Sigma\right],$ (138)
$\bar{\Omega}_{2}=\frac{34103}{18144}+\frac{13661}{2016}\nu+\frac{59}{18}\nu^{2}-\frac{79}{8}\nu\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\cdot\frac{\bm{S}_{2}}{m_{2}^{2}}\right)+\frac{81}{16}\left[\frac{m_{1}^{2}}{m^{2}}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\right)^{2}+\frac{m_{2}^{2}}{m^{2}}\left(\frac{\bm{S}_{2}}{m_{2}^{2}}\right)^{2}\right].$
(139)

The dCS modification takes the form of

$\delta{\Omega}=\zeta\left\\{\frac{16775}{12288}\frac{1}{\nu}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\cdot\frac{\bm{S}_{2}}{m_{2}^{2}}\right)-\frac{19271}{24576}\left[\frac{m^{2}}{m_{1}^{2}}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\right)^{2}+\frac{m^{2}}{m_{2}^{2}}\left(\frac{\bm{S}_{2}}{m_{2}^{2}}\right)^{2}\right]\right\\}.$
(140)

With the energy dissipation, the binding energy of the binary system
decreases. Combined with modified Kepler’s third law, one can calculate the
ppE parameters directly without using the waveform. The orbital radius of the
binary becomes smaller and smaller, and thus the angular frequency, ($\omega$
or $x$), becomes larger and larger. When the orbital radius is close to the
gravitational radius of the BH, the near-zone gravitational field is so strong
that the PN approximation is no longer valid. Let the merging time of the
binary BHs be $t_{c}$, and define the time from merging time as
$\tau=t_{c}-t$. In terms of $\tau$, above differential equation (136) is
solved perturbatively as

$\tau=\int_{\tau}^{0}d\tau=\frac{5}{256}\frac{m}{\nu}\frac{1}{x^{4}}\\\
\left[1-\frac{4}{3}x\bar{\Omega}_{1}+\frac{8}{5}x^{3/2}\bar{\Omega}_{1.5}+2x^{2}\left(\bar{\Omega}_{1}^{2}-\bar{\Omega}_{2}-\delta\Omega\right)\right].$
(141)

Inversely, the solution to $x(\tau)$ is

$x=\frac{1}{4}\Theta^{-1/4}\left\\{1-\frac{1}{12}\bar{\Omega}_{1}\Theta^{-1/4}-\frac{1}{20}\bar{\Omega}_{1.5}\Theta^{-3/8}+\frac{1}{288}\left[8\bar{\Omega}_{1}^{2}-9(\bar{\Omega}_{2}+\delta\Omega)\right]\Theta^{-1/2}\right\\},$
(142)

where the dimensionless time is defined as $\Theta(\tau)\equiv\nu\tau/5m$,
relating to the orbital phase by

$x^{3/2}=-\frac{\nu}{5}\frac{d\phi}{d\Theta}.$ (143)

Integrating Eq. (143) gives

$\displaystyle\phi(\Theta)=\phi_{c}-\frac{5}{\nu}\int^{\Theta}_{0}x^{3/2}(\Theta)d\Theta=\phi_{c}-\frac{1}{\nu}\Theta^{5/8}\left[1+\bar{\kappa}_{1}\Theta^{-1/4}+\bar{\kappa}_{1.5}\Theta^{-3/8}+(\bar{\kappa}_{2}+\delta\kappa)\Theta^{-1/2}\right],$
(144)

where $\phi_{c}$ is the initial orbital phase. The higher-order PN
coefficients are given by

$\bar{\kappa}_{1}=\frac{3715}{8064}+\frac{55}{96}\nu,$ (145)
$\bar{\kappa}_{1.5}=-\frac{3\pi}{4}+\frac{1}{m^{2}}\left[\frac{47}{16}S+\frac{75}{64}\frac{\delta
m}{m}\Sigma\right],$ (146)

and

$\bar{\kappa}_{2}=\frac{9275495}{14450688}+\frac{284875}{258048}\nu+\frac{1855}{2048}\nu^{2}-\frac{1185}{512}\nu\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\cdot\frac{\bm{S}_{2}}{m_{2}^{2}}\right)-\frac{1215}{1024}\left[\frac{m_{1}^{2}}{m^{2}}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\right)^{2}+\frac{m_{2}^{2}}{m^{2}}\left(\frac{\bm{S}_{2}}{m_{2}^{2}}\right)^{2}\right],$
(147)

The dCS correction is

$\delta\kappa=\zeta\left\\{-\frac{83875}{262144}\frac{1}{\nu}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\cdot\frac{\bm{S}_{2}}{m_{2}^{2}}\right)+\frac{96355}{524288}\left[\frac{m^{2}}{m_{1}^{2}}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\right)^{2}+\frac{m^{2}}{m_{2}^{2}}\left(\frac{\bm{S}_{2}}{m_{2}^{2}}\right)^{2}\right]\right\\}.$
(148)

Similar to the discussions around Eqs. (83) and (135), for the equal-mass and
equal spin binary BH system, the ratio of correction of orbital phase caused
by SS and MQ effect is about $0.87$. The contribution from the SS effect is
comparable with the MQ effect for orbital phase evolution.

##  VII Frequency-Domain Waveform and ppE parameters

In extracting GWs from noisy data, it is customary to employ the frequency-
domain gravitational waveform in Fourier space, which can be computed
analytically via the stationary phase approximation Maggiore (2008); Arun  _et
al._ (2009); Zhang  _et al._ (2017b, 2019a); Liu  _et al._ (2020). For
preparation, let us first write the detected GW signal in the following form,

$\displaystyle
H(t)=H_{+}F_{+}+H_{\times}F_{\times}=\frac{4\mu}{R}x(t_{r})\sum_{n=1}^{6}\sum_{k=0}^{4}x^{k/2}(t_{r})\times\left\\{c_{n,k}\cos\left[n\phi(t)\right]+s_{n,k}\sin\left[n\phi(t)\right]\right\\},$
(149)

where $F_{+}$ and $F_{\times}$ are antenna pattern functions of GW detectors
Maggiore (2008), the coefficients $c_{n,k}$ and $s_{n,k}$ are a set of
constants depending on the mass and spin of binary BHs. The index $k$ denotes
the order or PN expansion, and $n$ is the multiple of GW frequency relative to
the orbital frequency $\omega$. These constants can be read out from the
gravitational waveforms Blanchet  _et al._ (1995); Kidder (1995); Buonanno
_et al._ (2013); Bohé  _et al._ (2015). The Fourier transformation of the
signal (149) is

$\displaystyle\tilde{H}(f)=\int_{-\infty}^{\infty}H(t)e^{i2\pi
ft}dt=\frac{2\mu}{R}e^{i2\pi
fR}\sum_{n=1}^{6}\sum_{k=0}^{4}\int_{-\infty}^{\infty}dtx^{k/2+1}e^{i2\pi
ft}\left[(c_{n,k}-is_{n,k})e^{in\phi}+(c_{n,k}+is_{n,k})e^{-in\phi}\right],$
(150)

where $f$ is the oscillation frequency of GW signals. Eq. (150) is an
oscillation integration, mainly contributed by the value of the integrated
function at the stationary point. The stationary points are determined by
equations $n\dot{\phi}=2\pi f$ and $n\dot{\phi}=-2\pi f$. The former equation
gives

$\displaystyle\Theta^{*}_{n}=\frac{1}{256}u_{n}^{-8/3}\left\\{1+\frac{32}{5}\kappa_{1}u_{n}^{2/3}+\frac{128}{15}\kappa_{1.5}u_{n}+\frac{64}{75}\left[3\kappa_{1}^{2}+10(\kappa_{2}+\delta\kappa)\right]u_{n}^{4/3}\right\\},$
(151)

while the latter will not give a stationary point. Eq. (151) introduces
another dimensionless frequency as

$u_{n}\equiv\frac{2}{n}\pi mf$ (152)

for convenience. The frequency-domain waveform is thus obtained,

$\displaystyle\tilde{H}(f)$
$\displaystyle=\sqrt{\frac{5}{12}}\frac{m^{2}}{R}\sum_{n=1}^{6}\sum_{k=0}^{4}\sqrt{\frac{\pi\nu}{n}}(c_{n,k}+is_{n,k})\times
u_{n}^{(2k-7)/6}e^{i\Psi_{n}(f)}$ (153)
$\displaystyle\qquad\times\left\\{1+\frac{12}{5}\kappa_{1}u_{n}^{2/3}+\frac{8}{3}\kappa_{1.5}u_{n}+\frac{32}{15}\left[(\kappa_{2}+\delta\kappa)-\frac{21}{20}\kappa_{1}^{2}\right]u_{n}^{4/3}\right\\}.$

The frequency-domain phase is

$\Psi_{n}(f)=2\pi
f(t_{c}+R)-\frac{\pi}{4}-n\phi_{c}+\frac{3}{256}\frac{n}{\nu}u_{n}^{-5/3}\left\\{1+\frac{32}{3}\kappa_{1}u_{n}^{2/3}+\frac{64}{3}\kappa_{1.5}u_{n}+\frac{64}{15}\left[3\kappa_{1}^{2}+10(\kappa_{2}+\delta\kappa)\right]u_{n}^{4/3}\right\\}.$
(154)

So far, we have completed the calculation of ready-to-use gravitational
waveform radiated from binary slowly-rotating BHs in dCS gravity.

We rewrite this frequency-domain waveform in terms of ppE parameters,

$\tilde{H}(f)=\tilde{h}(f)\left[1+\alpha_{\rm g}\cdot u^{a_{\rm
g}}\right]\exp\left[i\beta_{\rm g}\cdot u^{b_{\rm g}}\right],$ (155)

where $\tilde{h}(f)$ is the GR waveform. $u\equiv(\pi\mathcal{M}f)^{1/3}$ is
another dimensionless frequency, with chirp mass defined as
$\mathcal{M}=\nu^{3/5}m$. The ppE parameters are

$a_{\rm g}=4,\quad b_{\rm g}=-1,$ (156)

and

$\displaystyle\alpha_{\rm g}$
$\displaystyle=\nu^{-4/5}\cdot\zeta\left\\{-\frac{11975}{24576}\frac{1}{\nu}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\cdot\frac{\bm{S}_{2}}{m_{2}^{2}}\right)+\frac{96305}{344064}\left[\frac{m^{2}}{m_{1}^{2}}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\right)^{2}+\frac{m^{2}}{m_{2}^{2}}\left(\frac{\bm{S}_{2}}{m_{2}^{2}}\right)^{2}\right]\right\\},$
(157) $\displaystyle\beta_{\rm g}$
$\displaystyle=\nu^{-4/5}\cdot\zeta\left\\{-\frac{83875}{262144}\frac{1}{\nu}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\cdot\frac{\bm{S}_{2}}{m_{2}^{2}}\right)+\frac{96355}{524288}\left[\frac{m^{2}}{m_{1}^{2}}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\right)^{2}+\frac{m^{2}}{m_{2}^{2}}\left(\frac{\bm{S}_{2}}{m_{2}^{2}}\right)^{2}\right]\right\\}.$

The dCS theory not only corrects the phase of GWs predicted by GR but also
corrects its amplitude. This result (155) gives the gravitational waveform
radiated by the spin-aligned circular-orbits binary BHs in the wave zone,
where the effects of cosmological background evolution and parity-violating
amplitude birefringence effect are not considered. This waveform can be
directly used to extract GW signals, estimate the constrain ability of the
future GW detector on dCS theory and constrain the dCS theory through the
current observation data. The same results will be given by a simple formulas,
$\alpha_{\rm g}=(2/3)\delta\varpi-(1/2)\delta\Omega$ and $\beta_{\rm
g}=(15/64)\delta\Omega$, provided by Ref. Tahura and Yagi (2018), in which the
different result is given, because the tensor-deformation flux is not
included.

##  VIII Propagation Effect of GWs

As mentioned earlier, the inspiraling binary BH system generates only two
polarization modes of GWs. Due to the time scale of the binary merge being too
short, the parity-violating effect will not appear in the wave zone. However,
when the GW enters the propagation zone, the background scalar field will
induce the amplitude birefringence effect, and the two polarization modes will
transform each other.

###  VIII.1 Evolution Equation of Circular Polarization Modes

The background metric is assumed as flat Robertson-Walker geometry,

$ds^{2}=g_{\mu\nu}^{\rm(B)}dx^{\mu}dx^{\nu}=a^{2}(\eta)(-d\eta^{2}+d\mathbf{x}^{2}).$
(158)

$\eta$ denotes the conformal time, and $a(\eta)$ denotes the scale factor of
the Universe. Throughout this section, we set the present scale factor
$a_{0}\equiv a(\eta_{0})=1$. conformal time relates to the cosmic time $t$ by
$dt=a(\eta)d\eta$. The metric perturbation is also denoted by $H_{\mu\nu}$,
which satisfies the TT gauge

$H_{0\mu}=0,\quad\partial^{i}H_{ij}=0,\quad H=g^{ij}H_{ij}=0.$ (159)

From the total metric $g^{\rm(B)}_{\mu\nu}+H_{\mu\nu}$ and field equation (2),
the evolution equation of $H_{\mu\nu}$ is given by

$H^{\prime\prime}_{jk}+2\mathcal{H}H^{\prime}_{jk}-(\Box^{2}_{\eta}H_{jk})+\frac{16\pi\alpha}{a^{2}}\cdot\epsilon^{sm}_{\
\ \
(j}\partial_{m}\left[\vartheta^{\prime}_{0}H_{sk)}^{\prime\prime}+\vartheta^{\prime\prime}_{0}H_{sk)}^{\prime}-\vartheta^{\prime}_{0}(\Box^{2}_{\eta}H_{sk)})\right]=0.$
(160)

In above Eq. (160), $\mathcal{H}\equiv a^{\prime}/a$ is the Hubble parameter.
A prime $(\prime)$ denotes the derivative with respect to the conformal time
$\eta$. $\vartheta_{0}$ is the background scalar field, which homogeneously
distributes in the entire Universe. It is just a function of cosmic time and
changes with the expansion of the Universe. Because of the background scalar
field, GWs exhibit an amplitude birefringence effect during propagation, which
violates the prediction by GR. We decompose the metric perturbation according
to different wave vectors and divide it into left- and right-hand polarization
modes,

$H_{jk}=\sum_{A=L,R}\int\frac{d^{3}\bm{\kappa}}{(2\pi)^{3}}H_{A}(\eta,\kappa_{i})\exp\left(i\bm{\kappa}\cdot\mathbf{x}\right)e_{jk}^{A}.$
(161)

The subscripts $A=L,R$ represent two kinds of circular modes. $\bm{\kappa}$ is
the wavevector of propagating GWs, and we define
$\kappa^{2}\equiv\bm{\kappa}\cdot\bm{\kappa}$ and
$\hat{\bm{\kappa}}\equiv\bm{\kappa}/\kappa$. The circular polarization basis
tensors are defined as $e_{ij}^{L,R}=(1/\sqrt{2})(e_{ij}^{+}\mp
ie_{ij}^{\times})$. $e^{+,\times}_{ij}$ describe the plus and cross
polarization modes. Additionally, the basis tensors satisfy

$\hat{\kappa}_{m}\epsilon^{sm}_{\ \ \ j}e^{A}_{sk}=i\rho_{A}e_{jk}^{A}$ (162)

with $\rho_{L}=-1$ and $\rho_{R}=+1$. Decomposition (161) gives two
independent equations of circular polarizations,

$H^{\prime\prime}_{A}+2\mathcal{H}H^{\prime}_{A}+\kappa^{2}H_{A}-\frac{16\pi\alpha\kappa}{a^{2}}\rho_{A}\left[\vartheta^{\prime}_{0}H_{A}^{\prime\prime}+\vartheta^{\prime\prime}_{0}H_{A}^{\prime}+\kappa^{2}\vartheta^{\prime}_{0}H_{A}\right]=0,$
(163)

or equivalently,

$H_{A}^{\prime\prime}+\left(2+\nu_{A}\right)\mathcal{H}H_{A}^{\prime}+\kappa^{2}H_{A}=0,\quad\text{with}\quad\nu_{A}=\frac{(16\pi\alpha)\rho_{A}\kappa(2\mathcal{H}\vartheta^{\prime}_{0}-\vartheta^{\prime\prime}_{0})/(a^{2}\mathcal{H})}{1-(16\pi\alpha)\rho_{A}\kappa\vartheta^{\prime}_{0}/a^{2}}.$
(164)

$\nu_{A}$ depends on the modes and describes the amplitude birefringence
effect, that is, the amplitudes of left- and right-hand waves will gradually
increase (decrease) during propagation. It can be seen that the dCS theory
does not modify the wave number of GWs, so GWs still travel in vacuum at the
speed of light, whether left- or right-hand polarization mode.

###  VIII.2 Amplitude Birefringence

The amplitude corrections to the gravitational waveform arising from the
parameters $\nu_{A}$ are given from the evolution equation (164). We further
decompose circular polarizations as

$H_{A}=\mathcal{A}_{A}e^{-i\Phi_{A}(\eta)}e^{-i\theta_{A}(\eta)}.$ (165)

$\mathcal{A}_{A}$ and $\Phi_{A}$ represent the original amplitude and phase of
circular polarization waveforms. The phase factor $e^{-i\theta}$ denotes
amplitude correction induced by parameter $\nu_{A}$. Substituting Eq. (165)
into the perturbation equation, we give the equation of $\theta_{A}$,

$i\theta^{\prime\prime}+2\theta^{\prime}\Phi^{\prime}+\theta^{\prime
2}+i\left(2+\nu_{A}\right)\mathcal{H}\theta^{\prime}+i\mathcal{H}\nu_{A}\Phi^{\prime}=0,$
(166)

Because of $\theta_{A}\sim\nu_{A}$ and $d/d\eta\sim\mathcal{H}$, we have
$\theta^{\prime}_{A}\sim\mathcal{H}\theta$ and
$\theta^{\prime\prime}_{A}\sim\mathcal{H}^{2}\theta$. At the same time, the
wavelength is assumed to be much smaller than the scale of the Universe, such
that $\kappa\gg\mathcal{H}$. Finally, Eq. (166) is approximated as

$2\theta^{\prime}\Phi^{\prime}+i\mathcal{H}\nu_{A}\Phi^{\prime}=0.$ (167)

This is integrated from the conformal time of the GW source to the present
time. The solution to $\theta_{A}$ is obtained,

$\theta_{A}=-\frac{i}{2}\int_{\eta_{s}}^{\eta_{0}}\mathcal{H}\nu_{A}d\eta.$
(168)

Compared with the original one, $\mathcal{A}_{A}e^{-i\Phi_{A}(\eta)}$, the
detected waveform is written as

$H_{A}^{\rm(o)}=H_{A}(1+\delta H_{A}),$ (169)

where $\delta H_{A}$ is the amplitude correction for each polarization mode,
which is calculated as follows,

$1+\delta
H_{A}=e^{-i\theta_{A}}=1-\frac{1}{2}\int_{\eta_{s}}^{\eta_{0}}\mathcal{H}\nu_{A}d\eta=1-\ln\sqrt{\frac{1-(16\pi\alpha)\rho_{A}\kappa[\vartheta^{\prime}_{0}(\eta_{0})/a^{2}(\eta_{0})]}{1-(16\pi\alpha)\rho_{A}\kappa[\vartheta^{\prime}_{0}(\eta_{s})/a^{2}(\eta_{s})]}}\approx
1+\rho_{A}\xi u^{3},$ (170)

where

$\xi\equiv\frac{16\pi\alpha}{\mathcal{M}}\left[\vartheta^{\prime}_{0}(\eta_{0})-(1+z)^{2}\vartheta^{\prime}_{0}(\eta_{s})\right].$
(171)

The cosmological redshift of the wave source is defined as $z\equiv
1/a(\eta_{s})-1$. In terms of plus and cross modes, the detected polarizations
are written as

$H_{+}^{\rm(o)}=H_{+}-i\xi u^{3}H_{\times},\quad
H_{\times}^{\rm(o)}=H_{\times}+i\xi u^{3}H_{+}.$ (172)

Amplitude birefringence will enter at least 1.5PN order correction of the
waveform. Like the gravitational Faraday rotation of GWs Li  _et al._ (2022b),
the polarization tensor will rotate an angle on the polarization plane, and
the plus and cross modes then transform each other.

###  VIII.3 The Detected Signal

The detected GW signal is a linear combination of plus and cross modes,

$H^{\rm(o)}=H_{+}^{\rm(o)}F_{+}+H_{\times}^{\rm(o)}F_{\times}=H_{+}F_{+}+H_{\times}F_{\times}-i\xi
u^{3}\left(H_{\times}F_{+}-H_{+}F_{\times}\right).$ (173)

During the process of GW generation, the dCS theory influences the waveform at
2PN order. However, in the propagation zone, the leading-order modification by
amplitude birefringence enters the 1.5PN waveform, while the modified 0.5PN
waveform enters the 2PN order. In Fourier space, the waveform is modified as

$\tilde{H}^{\rm(o)}(f)=\tilde{H}(f)-i\xi\tilde{h}_{\xi}(f).$ (174)

The first term is given in Eq. (150). Up to 2PN order, the second term is

$\tilde{h}_{\xi}(f)=\sqrt{\frac{5\pi}{12}}\frac{m^{2}}{R}\nu^{6/5}u^{-1/2}\Bigg{\\{}\alpha_{02}e^{i\beta_{2}(f)}+2^{-5/6}\nu^{-1/5}u\left[\alpha_{11}e^{i\beta_{1}(f)}+\alpha_{13}e^{i\beta_{3}(f)}\right]\Bigg{\\}}e^{2\pi
f(t_{c}+R)-\pi/4},$ (175)

without dCS modification. The amplitude correction parameters are

$\displaystyle\alpha_{02}$
$\displaystyle=\frac{1}{\sqrt{2}}\left[\frac{1}{2}(1+\cos^{2}\iota)F_{\times}-i\cos\iota
F_{+}\right],$ (176a) $\displaystyle\alpha_{11}$
$\displaystyle=\frac{1}{4}\frac{\delta
m}{m}\sin\iota\left[(5+\cos^{2}\iota)F_{\times}-3i\cos\iota F_{+}\right],$
(176b) $\displaystyle\alpha_{13}$
$\displaystyle=-\frac{9}{4}\sqrt[3]{3}\frac{\delta
m}{m}\sin\iota\left[(1+\cos^{2}\iota)F_{\times}-i\cos\iota F_{+}\right].$
(176c)

The phase correction parameter is

$\beta_{n}(f)=-n\left[\Phi_{c}+\frac{3}{256}(n/2)^{5/3}u^{-5}\right].$ (177)

The dCS modification during GW generation and propagation can be put together.
The detected frequency-domain signal is then written as

$\tilde{H}^{\rm(o)}(f)=\tilde{h}(f)\cdot(1+\alpha_{\rm g}\cdot
u^{4})\cdot\exp(i\beta_{\rm g}\cdot u^{-1})-i\xi\tilde{h}_{\xi}(f).$ (178)

The parameters $\alpha_{\rm g}$ and $\beta_{\rm g}$ describe correction to the
generation process, which enters 2PN approximation, measured by coupling
$\alpha$ or $\zeta$. Ref. Silva  _et al._ (2021) constrains $\sqrt{\alpha}\leq
8.5\,{\rm km}$ by multimessenger neutron star observations. Another term in
Eq. (178) describes amplitude correction induced by parity-violating effect,
which enters 1.5PN and 2PN corrections. This effect is measured by the
coupling parameter and the time derivative of the background scalar field,
$\alpha\vartheta^{\prime}_{0}$. Ref. Smith  _et al._ (2008) used solar-system
measurements of frame-dragging from LAGEOS and Gravity Probe B to bound
$\ell_{0}\equiv\alpha\dot{\vartheta}_{0}/16\pi\leq 3000\,{\rm km}$. Ref.
Okounkova  _et al._ (2022) improved this constraint to $\ell_{0}\leq
1000\,{\rm km}$ using GWTC-2 catalog Abbott  _et al._ (2021a).

##  IX Conclusion and Discussion

DCS theory is a parity-violating gravitational theory that has recently
attracted more and more attention. In this article, we calculate the
gravitational waveform radiated by binary slowly-rotating BHs, where the
polarization modes of GWs are studied in the generation and propagation
process. The EOM is established from the modified MPD equation, which is a
powerful tool to describe the motion of spinning point mass based on the EFT
method. Compared with Ref. Yagi  _et al._ (2012b), we emphasize that the EOM
is derived for spin-aligned quasi-circular orbits case, in which the SAM of
each object, $\bm{S}_{A}$, is conserved, without spin precession. The complete
relative acceleration of the binary system, $\bm{a}$, is given in Eq. (67),
including non-spinning effect, SO, SS, and MQ coupling. The dCS modification,
$\delta\bm{a}$, is similar to SS and MQ effects predicted by GR, which enters
2PN corrections.

Based on the EOM of binary, the scalar and tensor-deformation radiations,
$\vartheta$ and $\bar{k}_{ij}$, are calculated through the multipole moment
formula. Through NP tetrad, the polarization modes are investigated in this
work. We find that DCS theory cannot produce any extra scalar and vector modes
since dCS is a quadratic and massless theory, although a scalar DOF is
introduced in the action. The plus and cross modes, $k_{+}$ and $k_{\times}$,
are presented in Eq. (116). In the wave zone, there is no parity-violating
effect for plus and cross modes because of the short time scale of the binary
merger, in which the expansion of the Universe and homogeneous background
scalar field are not taken into consideration.

The energy flux carried by radiation consists of scalar and tensor fluxes.
Compared with Refs. Yagi  _et al._ (2012b, 2016), the flux of metric
deformation, $\delta\mathcal{F}_{H}$, is also included, whose leading-order
modification is at 2PN order, same as the scalar flux, rather than 3PN.
Following the standard PN method, the orbital decay $\dot{x}$, frequency-
domain waveform $\tilde{H}(f)$ and improved ppE parameters $(\alpha_{\rm
p},\beta_{\rm p})$ are also given in Eqs. (136), (150) and (155). The GWs need
to travel for a long time before reaching the detector. The effects of
cosmological metric and background scalar field $\vartheta_{0}$ are
significant in this process. In dCS gravity, the amplitudes of the left-hand
circular polarization mode of GW, $H_{L}$, will increase (decrease) during its
propagation, while the right-hand mode, $H_{R}$, will decrease (increase),
which is called amplitude birefringence (172). Thus, more amplitude and phase
corrections are introduced into the detected frequency-domain waveform (178).
This effect enters 1.5PN correction, which is more obvious than that produced
in the generation process. This result may be helpful for the joint constraint
on the background scalar field $\theta_{0}$ and dCS coupling parameter
$\alpha$. The future third-generation ground-based detectors, including
Einstein Telescope and Cosmic Explorer, and the space-based GW detectors, like
LISA, TianQin and Taiji, will hopefully detect the higher PN waveform Shi  _et
al._ (2022). Our results can be used for GW signal extraction, detector
ability prediction, and gravitational test in future work.

###### Acknowledgements.

We would like to thank Kent Yagi, Alessandra Buonanno, Guillaume Faye,
Nicholas Loutrel, Sharaban Tahura, Shaoqi Hou, Rui Niu and Yupeng Zhang for
helpful discussions and comments. This work is supported by the National Key
R&D Program of China Grant No. 2022YFC2200100 and 2021YFC2203102, NSFC No.
12273035 and 11903030 the Fundamental Research Funds for the Central
Universities under Grant No. WK2030000036 and WK3440000004, and the science
research grants from the China Manned Space Project with No.CMS-CSST-2021-B01,
CMS-CSST-2021-B11. T. Z. is supported in part by the National Key Research and
Development Program of China under Grant No. 2020YFC2201503, the National
Natural Science Foundation of China under Grant No. 12275238 and No. 11675143,
the Zhejiang Provincial Natural Science Foundation of China under Grant No.
LR21A050001 and LY20A050002, and the Fundamental Research Funds for the
Provincial Universities of Zhejiang in China under Grant No. RF-A2019015. T.
L. is supported by the National Natural Science Foundation of China (NSFC)
Grant No. 12003008 and the China Postdoctoral Science Foundation Grant No.
2020M682393.

##  Appendix A Hadamard Regularization

Let us consider the class of functions $F$ depending on the field point
$\mathbf{x}$ as well as on two source points $\mathbf{z}_{1}$ and
$\mathbf{z}_{2}$, and admitting, when the field point approaches one of the
source points ($r_{1}\rightarrow|\mathbf{x}-\mathbf{z}_{1}|\rightarrow 0$ for
instance), an expansion of the type

$F(\mathbf{x})=\sum_{-k_{0}\leqslant k\leqslant
0}r_{1}^{k}f_{k}(\mathbf{x};\mathbf{z}_{1},\mathbf{z}_{2}).$ (179)

We define the value of the function $F$ at the location $1$ to be the so-
called Hadamard finite part, which is the average, with respect to the
direction $\hat{\bm{n}}_{1}$, of the approach to point $1$, of the term with
lowest-order coefficient in Eq. (179),

$F_{1}\equiv\frac{1}{4\pi}\int
f_{0}(\mathbf{z}_{1};\mathbf{z}_{1},\mathbf{z}_{2})d\Omega_{1},$ (180)

with $d\Omega_{1}$ a solid angle element centered on the location of the
$1$-st object.

To derive the EOM, we need to give the Hadamard-regularized Newtonian
potential $\bar{U}$, dCS scalar field $\vartheta$, and dCS-modified potential
$\delta U$. These function can be rewritten via Eq. (179) at the location of
$1$-st object,

$\displaystyle\bar{U}=\sum_{-1\leqslant k\leqslant
0}r_{1}^{k}\bar{u}_{k}(\mathbf{x};\mathbf{z}_{1},\mathbf{z}_{2}),\quad$
$\displaystyle\text{with}\quad\bar{u}_{-1}=m_{1}\quad\text{and}\quad\bar{u}_{0}=\frac{m_{2}}{r_{2}}=\frac{m_{1}}{|\mathbf{x}-\mathbf{z}_{2}|},$
(181a) $\displaystyle\vartheta=\sum_{-2\leqslant k\leqslant
0}r_{1}^{k}\tilde{\vartheta}_{k}(\mathbf{x};\mathbf{z}_{1},\mathbf{z}_{2}),\quad$
$\displaystyle\text{with}\quad\tilde{\vartheta}_{0}=-\frac{\bm{n}\cdot\bm{\mu}_{2}}{r_{2}^{2}}=-\frac{\bm{n}\cdot\bm{\mu}_{2}}{|\mathbf{x}-\mathbf{z}_{2}|^{2}},$
(181b) $\displaystyle\delta U=\sum_{-3\leqslant k\leqslant 0}r_{1}^{k}\delta
u_{k}(\mathbf{x};\mathbf{z}_{1},\mathbf{z}_{2})\quad$
$\displaystyle\text{with}\quad\delta
u_{0}=\frac{603}{3584}\frac{\zeta_{2}}{m_{2}^{4}}\left(\frac{m_{2}}{r_{2}}\right)^{3}\left[(\hat{\bm{n}}_{2}\cdot\bm{S}_{2})^{2}-\frac{1}{3}(\bm{S}_{2}\cdot\bm{S}_{2})\right].$
(181c)

Such that the finite values of the Newtonian potential, scalar field, and dCS
potential at location $1$ are

$\displaystyle\bar{U}_{1}$
$\displaystyle=\frac{1}{4\pi}\int\bar{u}_{0}(\mathbf{z}_{1};\mathbf{z}_{1},\mathbf{z}_{2})d\Omega_{1}=\frac{m_{1}}{|\mathbf{z}_{1}-\mathbf{z}_{2}|}=\frac{m_{1}}{r},$
(182a) $\displaystyle\vartheta_{1}$
$\displaystyle=\frac{1}{4\pi}\int\tilde{\vartheta}_{0}(\mathbf{z}_{1};\mathbf{z}_{1},\mathbf{z}_{2})d\Omega_{1}=\frac{\bm{n}\cdot\bm{\mu}_{2}}{|\mathbf{z}_{1}-\mathbf{z}_{2}|^{2}}=\frac{\bm{n}\cdot\bm{\mu}_{2}}{r^{2}},$
(182b) $\displaystyle\delta{U}_{1}$
$\displaystyle=\frac{1}{4\pi}\int\delta{u}_{0}(\mathbf{z}_{1};\mathbf{z}_{1},\mathbf{z}_{2})d\Omega_{1}=\frac{603}{3584}\frac{\zeta_{2}}{m_{2}^{4}}\left(\frac{m_{2}}{r}\right)^{3}\left[(\hat{\bm{n}}\cdot\bm{S}_{2})^{2}-\frac{1}{3}(\bm{S}_{2}\cdot\bm{S}_{2})\right],$
(182c)

respectively. The results for the location $2$, $\bar{U}_{2}$,
$\vartheta_{2}$, and $\delta U_{2}$, are similar to Eq. (182). One can find
that the finite part of Newtonian potential at location $1$ is just the value
of Newtonian potential generated by the $2$-nd object at the $1$-st position,
which is equivalent to omitting the self field. These finite results are used
to derive the EOM of a black hole in the binary system (see Eq. (62)), which
is the same as that given in Ref. Yagi  _et al._ (2012b); Loutrel  _et al._
(2018), where the self-field contributions are omitted directly.

##  Appendix B Estimation to the PN Order of Scalar Radiation

In this appendix, we plan to analyze the PN order of scalar radiation.

The multipole scalar radiations from non-compact sources term are estimated as

$\displaystyle\vartheta^{(\sigma_{1})}_{\rm mon}$
$\displaystyle=-\frac{1}{R}\frac{1}{2\pi}\int
U(\partial_{i}\partial_{i}\vartheta^{\rm(B)})d^{3}\mathbf{x}\propto\frac{m}{R}\frac{\mu}{r^{2}}\sim\mathcal{O}(1/c^{6}),$
(183a) $\displaystyle\vartheta^{(\sigma_{1})}_{\rm dip}$
$\displaystyle=-\frac{1}{R}\frac{1}{2\pi}\frac{\partial}{\partial t}\int
U(\partial_{i}\partial_{i}\vartheta^{\rm(B)})(\hat{\mathbf{N}}\cdot\mathbf{x})d^{3}\mathbf{x}\propto\frac{m}{R}\frac{\mu}{r^{2}}v\sim\mathcal{O}(1/c^{7}),$
(183b) $\displaystyle\vartheta^{(\sigma_{1})}_{\rm quad}$
$\displaystyle=-\frac{1}{R}\frac{1}{4\pi}\frac{\partial^{2}}{\partial
t^{2}}\int
U(\partial_{i}\partial_{i}\vartheta^{\rm(B)})(\hat{\mathbf{N}}\cdot\mathbf{x})^{2}d^{3}\mathbf{x}\propto\frac{m}{R}\frac{\mu}{r^{2}}v^{2}\sim\mathcal{O}(1/c^{8}),$
(183c)

and

$\displaystyle\vartheta^{(\sigma_{2})}_{\rm
mon}=\frac{2}{R}\frac{1}{\pi}\frac{\alpha}{\beta}\int\epsilon^{ijk}(\partial_{i}\partial_{m}U)(\partial_{j}\partial_{m}V_{k})d^{3}\mathbf{x}=\frac{2}{R}\frac{1}{\pi}\epsilon^{ijk}\frac{\alpha}{\beta}\oint\hat{N}^{i}\left[(\partial_{m}U)(\partial_{j}\partial_{m}V_{k})\right]R^{2}d\Omega\propto\frac{1}{R^{3}}\rightarrow
0,$ (184a) $\displaystyle\begin{aligned} \vartheta^{(\sigma_{2})}_{\rm
dip}&=\frac{2}{R}\frac{1}{\pi}\frac{\alpha}{\beta}\epsilon^{ijk}\frac{\partial}{\partial
t}\int(\partial_{i}\partial_{m}U)(\partial_{j}\partial_{m}V_{k})(\hat{\mathbf{N}}\cdot\mathbf{x})d^{3}\mathbf{x}\\\
&=\frac{2}{R}\frac{1}{\pi}\frac{\alpha}{\beta}\epsilon^{ijk}\frac{\partial}{\partial
t}\left\\{\oint\hat{N}_{j}[(\partial_{i}\partial_{m}U)(\partial_{m}V_{k})(\hat{\mathbf{N}}\cdot\mathbf{x})]R^{2}d\Omega\right\\}\propto\frac{1}{R^{2}}\rightarrow
0,\end{aligned}$ (184b) $\displaystyle\vartheta^{(\sigma_{2})}_{\rm
quad}=\frac{1}{R}\frac{1}{\pi}\frac{\alpha}{\beta}\epsilon^{ijk}\frac{\partial^{2}}{\partial
t^{2}}\int(\partial_{i}\partial_{m}U)(\partial_{j}\partial_{m}V_{k})(\hat{\mathbf{N}}\cdot\mathbf{x})^{2}d^{3}\mathbf{x}\propto\frac{m}{R}\frac{m}{r^{3}}v^{3}\sim\mathcal{O}(1/c^{11}).$
(184c)

The multipole scalar radiations from compact sources $\rho_{1}$ are estimated
as

$\displaystyle\vartheta^{(\rho_{1})}_{\rm mon}$
$\displaystyle=-\frac{1}{R}\sum_{A}\mu_{A}^{i}\int\partial_{i}\delta^{(3)}(\mathbf{x}-\mathbf{z}_{A}(t_{r}))d^{3}\mathbf{x}=0,$
(185a) $\displaystyle\vartheta^{(\rho_{1})}_{\rm dip}$
$\displaystyle=-\frac{1}{R}\frac{\partial}{\partial
t}\sum_{A}\mu_{A}^{i}\int[\partial_{i}\delta^{(3)}(\mathbf{x}-\mathbf{z}_{A}(t_{r}))](\hat{\mathbf{N}}\cdot\mathbf{x})d^{3}\mathbf{x}=\frac{1}{R}\hat{N}^{j}\frac{\partial}{\partial
t}\sum_{A}\mu_{A}^{i}\delta_{ij}=0,$ (185b)
$\displaystyle\vartheta^{(\rho_{1})}_{\rm quad}$
$\displaystyle=-\frac{1}{R}\frac{1}{2}\frac{\partial^{2}}{\partial
t^{2}}\sum_{A}\mu_{A}^{i}\int[\partial_{i}\delta^{(3)}(\mathbf{x}-\mathbf{z}_{A}(t_{r}))](\hat{\mathbf{N}}\cdot\mathbf{x})^{2}d^{3}\mathbf{x}\propto\frac{m}{R}\frac{\mu}{r^{2}}\sim\mathcal{O}(1/c^{6}).$
(185c)

The multipole scalar radiations from $\rho_{2}$ are estimated as

$\displaystyle\vartheta^{(\rho_{2})}_{\rm mon}$
$\displaystyle=-\frac{1}{R}\sum_{A}\mu_{A}^{i}\dot{v}^{i}_{A}\int\delta^{(3)}(\mathbf{x}-\mathbf{z}_{A}(t_{r}))d^{3}\mathbf{x}\propto\frac{m}{R}\frac{\mu}{r^{2}}\sim\mathcal{O}(1/c^{6}),$
(186a) $\displaystyle\vartheta^{(\rho_{2})}_{\rm dip}$
$\displaystyle=-\frac{1}{R}\frac{\partial}{\partial
t}\sum_{A}\mu_{A}^{i}\dot{v}^{i}_{A}\int\delta^{(3)}(\mathbf{x}-\mathbf{z}_{A}(t_{r}))(\hat{\mathbf{N}}\cdot\mathbf{z})d^{3}\mathbf{x}\propto\frac{m}{R}\frac{\mu}{r^{2}}v\sim\mathcal{O}(1/c^{7}),$
(186b) $\displaystyle\vartheta^{(\rho_{2})}_{\rm quad}$
$\displaystyle=-\frac{1}{R}\frac{1}{2}\frac{\partial^{2}}{\partial
t^{2}}\sum_{A}\mu_{A}^{i}\dot{v}^{i}_{A}\int\delta^{(3)}(\mathbf{x}-\mathbf{z}_{A}(t_{r}))(\hat{\mathbf{N}}\cdot\mathbf{x})^{2}d^{3}\mathbf{x}\propto\frac{m}{R}\frac{\mu}{r^{2}}v^{2}\sim\mathcal{O}(1/c^{8}).$
(186c)

The multipole scalar radiations from $\rho_{3}$ are estimated as

$\displaystyle\vartheta^{(\rho_{3})}_{\rm mon}$
$\displaystyle=\frac{1}{R}\sum_{A}\mu_{A}^{i}v_{A}^{i}v_{A}^{j}\int\partial_{j}\delta^{(3)}(\mathbf{x}-\mathbf{z}_{A}(t_{r}))d^{3}\mathbf{x}=0,$
(187a) $\displaystyle\vartheta^{(\rho_{3})}_{\rm dip}$
$\displaystyle=\frac{1}{R}\frac{\partial}{\partial
t}\sum_{A}\mu_{A}^{i}v_{A}^{i}v_{A}^{j}\int\left[\partial_{j}\delta^{(3)}(\mathbf{x}-\mathbf{z}_{A}(t_{r}))\right](\hat{\mathbf{N}}\cdot\mathbf{x})d^{3}\mathbf{x}\propto\frac{m}{R}\frac{\mu}{r^{2}}v\sim\mathcal{O}(1/c^{7}),$
(187b) $\displaystyle\vartheta^{(\rho_{3})}_{\rm quad}$
$\displaystyle=\frac{1}{R}\frac{1}{2}\frac{\partial^{2}}{\partial
t^{2}}\sum_{A}\mu_{A}^{i}v_{A}^{i}v_{A}^{j}\int\left[\partial_{j}\delta^{(3)}(\mathbf{x}-\mathbf{z}_{A}(t_{r}))\right](\hat{\mathbf{N}}\cdot\mathbf{x})^{2}d^{3}\mathbf{x}\propto\frac{m}{R}\frac{\mu}{r^{2}}v^{2}\sim\mathcal{O}(1/c^{8}).$
(187c)

Combining all the above results, we find the leading-order scalar radiations
are given by $\vartheta_{\rm mon}^{(\sigma_{1})}$, $\vartheta_{\rm
mon}^{(\rho_{2})}$, and $\vartheta_{\rm quad}^{(\rho_{1})}$. All dipole
radiation is either strictly zero or higher than the order we consider.

##  Appendix C Waveform from Source Terms

The GW waveforms have given in subsection IV.2. In this appendix, we plan to
derive the same results through the method provided by Refs. Yagi  _et al._
(2012a, 2016). The monopole formula of metric-deformation radiation is

$\bar{k}_{ij}=\frac{4}{R}\int K_{ij}d^{3}\mathbf{x}.$ (188)

The source term $K_{ij}$ is given in Eq. (22). The spatial component of the
source of tensor perturbation is

$K_{ij}=kT_{ij}^{(m)}+(1+h)\delta
T_{ij}^{(m)}+T_{ij}^{(\vartheta)}+\frac{1}{16\pi}\tilde{\Lambda}_{ij}-2\alpha\tilde{C}_{ij}.$
(189)

The first and second terms in Eq. (189) are of order $v^{2}S^{2}$, where $v$
represents the velocity and $S$ represents the spin angular momentum of binary
BHs. Therefore, the contribution from these two terms will enter higher PN
order. The third term comes from the EMT of the scalar field,

$T_{ij}^{(\vartheta)}=\beta\left[(\partial_{i}\vartheta)(\partial_{j}\vartheta)-\frac{1}{2}\delta_{ij}(\partial_{k}\vartheta)(\partial_{k}\vartheta)\right],$
(190)

which will produce the part of SS (DD) coupling. The fourth term,

$\displaystyle\frac{1}{16\pi}\tilde{\Lambda}_{ij}=\frac{1}{4\pi}\left[-\delta_{ij}(\partial_{k}\bar{U})(\partial_{k}\delta
U)+(\partial_{i}\bar{U})(\partial_{j}\delta U)+(\partial_{i}\delta
U)(\partial_{j}\bar{U})\right],$ (191)

provides the MQ effect in gravitational waveforms. The last term is

$-2\alpha\tilde{C}_{ij}=-4\alpha\epsilon_{r(i|m}\partial_{m}(\Box^{2}_{\eta}\hat{V}_{j)})(\partial_{r}\vartheta)+8\alpha\epsilon_{r(j|m}(\partial_{m}\partial_{[i)}\hat{V}_{k]})(\partial_{k}\partial_{r}\vartheta),$
(192)

which implies the parity-violating effect of GWs in dCS theory. Later, we will
prove that the contribution of this term to GWs is zero. This shows again that
GWs will not show the parity-violating effect without considering the
cosmological expansion and background scalar field.

According to Eq. (188), the contribution from Eq. (190) is

$\displaystyle\frac{4}{R}\int T_{ij}^{(\vartheta)}d^{3}\mathbf{x}$
$\displaystyle=\frac{4\beta}{R}\int\left[(\partial_{i}\vartheta_{1})(\partial_{j}\vartheta_{2})-\frac{1}{2}\delta_{ij}(\partial_{k}\vartheta_{1})(\partial_{k}\vartheta_{2})\right]d^{3}\mathbf{x}+(1\leftrightarrow
2)$ (193)
$\displaystyle=\frac{4\beta}{R}\mu_{1}^{k}\mu_{2}^{l}\int\left[\partial_{i}\partial_{k}\left(\frac{1}{r_{1}}\right)\partial_{j}\partial_{l}\left(\frac{1}{r_{2}}\right)-\frac{1}{2}\delta_{ij}\partial_{q}\partial_{k}\left(\frac{1}{r_{1}}\right)\partial_{q}\partial_{l}\left(\frac{1}{r_{2}}\right)\right]d^{3}\mathbf{x}+(1\leftrightarrow
2)$
$\displaystyle=\frac{4\beta}{R}\mu_{1}^{k}\mu_{2}^{l}\Bigg{[}\partial^{(1)}_{ik}\partial^{(2)}_{jl}\int\frac{1}{r_{1}}\frac{1}{r_{2}}d^{3}\mathbf{x}-\frac{1}{2}\delta_{ij}\partial^{(1)}_{qk}\partial^{(2)}_{ql}\int\frac{1}{r_{1}}\frac{1}{r_{2}}d^{3}\mathbf{x}\Bigg{]}+(1\leftrightarrow
2)$
$\displaystyle=-\frac{8\pi\beta}{R}\mu_{1}^{k}\mu_{2}^{l}\left[\partial^{(1)}_{ik}\partial^{(2)}_{jl}r-\frac{1}{2}\delta_{ij}\partial^{(1)}_{qk}\partial^{(2)}_{ql}r\right]+(1\leftrightarrow
2)$
$\displaystyle=-\frac{4\mu}{R}\frac{75}{256}\zeta\gamma^{3}\frac{1}{\nu}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\cdot\frac{\bm{S}_{2}}{m_{2}^{2}}\right)\hat{n}^{i}\hat{n}^{j},$

and from Eq. (191) is

$\displaystyle\frac{4}{R}\int\frac{\tilde{\Lambda}_{ij}}{16\pi}d^{3}\mathbf{x}$
$\displaystyle=\frac{1}{R}\frac{1}{\pi}\int\left[(\partial_{i}\bar{U}_{1})(\partial_{j}\delta
U_{2})+(\partial_{j}\bar{U}_{1})(\partial_{i}\delta
U_{2})-\delta_{ij}(\partial_{k}\bar{U}_{1})(\partial_{k}\delta
U_{2})\right]d^{3}\mathbf{x}+(1\leftrightarrow 2)$ (194)
$\displaystyle=\frac{1}{R}\frac{1}{\pi}\delta
C_{Q}\frac{m_{1}}{m_{2}}S_{2}^{k}S_{2}^{l}\int\left[\partial_{k}\partial_{(i}\left(\frac{1}{r_{1}}\right)\partial_{j)}\partial_{l}\left(\frac{1}{r_{2}}\right)-\frac{1}{2}\delta_{ij}\partial_{q}\partial_{k}\left(\frac{1}{r_{1}}\right)\partial_{q}\partial_{l}\left(\frac{1}{r_{2}}\right)\right]d^{3}\mathbf{x}+(1\leftrightarrow
2)$ $\displaystyle=\frac{1}{R}\frac{1}{\pi}\delta
C_{Q}\frac{m_{1}}{m_{2}}S_{2}^{k}S_{2}^{l}\left[\partial^{(1)}_{ik}\partial^{(2)}_{jl}\int\frac{1}{r_{1}}\frac{1}{r_{2}}d^{3}\mathbf{x}-\frac{1}{2}\delta_{ij}\partial^{(1)}_{kq}\partial^{(2)}_{lq}\int\frac{1}{r_{1}}\frac{1}{r_{2}}d^{3}\mathbf{x}\right]+(1\leftrightarrow
2)$ $\displaystyle=-\frac{2}{R}\delta
C_{Q}\frac{m_{1}}{m_{2}}S_{2}^{k}S_{2}^{l}\left[\partial^{(1)}_{ik}\partial^{(2)}_{jl}r-\frac{1}{2}\delta_{ij}\partial^{(1)}_{qk}\partial^{(2)}_{ql}r\right]+(1\leftrightarrow
2)$
$\displaystyle=\frac{4\mu}{R}\gamma^{3}\frac{603}{3584}\zeta\left[\frac{m^{2}}{m_{1}^{2}}\left(\frac{\bm{S}_{1}}{m_{1}^{2}}\right)^{2}+\frac{m^{2}}{m_{2}^{2}}\left(\frac{\bm{S}_{2}}{m_{2}^{2}}\right)^{2}\right]\hat{n}^{i}\hat{n}^{j}.$

$(1\leftrightarrow 2)$ represents the exchange of labels $1$ and $2$. The
operators $\partial_{k}^{(1)}$ and $\partial_{k}^{(2)}$ are the derivative
with respect to the coordinates of objects, $\mathbf{z}_{1}$ and
$\mathbf{z}_{2}$. Here we only consider the contribution of cross-interaction
terms.

The GWs generated from Eq. (192) is

$\frac{4}{R}\int(-2\alpha\tilde{C}_{ij})d^{3}\mathbf{x}=\frac{8\alpha}{R}\int[\epsilon_{rjm}(\partial_{r}\partial_{k}\partial_{m}\partial_{i}\bar{V}_{k})+\epsilon_{rim}(\partial_{r}\partial_{k}\partial_{m}\partial_{j}\bar{V}_{k})]\vartheta^{\rm(B)}d^{3}\mathbf{x}=0.$
(195)

Combining Eq. (193) and (194), we obtain the gravitational waveforms of the
metric deformation,