# Effects of spin-orbit coupling on gravitational waveforms from a triaxial
non-aligned neutron star in a binary system

Wen-Fan Feng  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  Tan Liu<EMAIL_ADDRESS>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  School of Physics, Hubei University, Wuhan 430062, China  Jie-
Wen Chen  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  National Time Service Center, Chinese Academy
of Sciences, Xi’an 710600, China  Key Laboratory of Time and Frequency Primary
Standards, Chinese Academy of Sciences, Xi’an 710600, China  Yan Wang
<EMAIL_ADDRESS>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  Soumya D. Mohanty  Department of Physics
and Astronomy, University of Texas Rio Grande Valley, Brownsville, Texas
78520, USA  Department of Physics, IIT Hyderabad, Kandai, Telangana-502284,
India

###### Abstract

Spinning neutron stars (NSs) can emit continuous gravitational waves (GWs)
that carry a wealth of information about the compact object. If such a signal
is detected, it will provide us with new insight into the physical properties
of matter under extreme conditions. According to binary population synthesis
simulations, future space-based GW detectors, such as LISA and TianQin, can
potentially detect some double NSs in tight binaries with orbital periods
shorter than 10 minutes. The possibility of a successful directed search for
continuous GWs from the spinning NS in such a binary system identified by
LISA/TianQin will be significantly increased with the proposed next-generation
ground-based GW observatories, such as Cosmic Explorer and Einstein Telescope.
Searching for continuous GWs from such a tight binary system requires highly
accurate waveform templates that account for the interaction of the NS with
its companion. In this spirit, we derive analytic approximations that describe
the GWs emitted by a triaxial non-aligned NS in a binary system in which the
effects of spin-orbit coupling have been incorporated. The difference with the
widely used waveform for the isolated NS is estimated and the parameter
estimation accuracy of an example signal using Cosmic Explorer is calculated.
For a typical tight double NS system with a 6 min orbital period, the angular
frequency correction of the spinning NS in this binary due to spin precession
is $\sim 10^{-6}~{}{\rm Hz}$, which is in the same order of magnitude as the
angular frequency of orbital precession. The fitting factor between the
waveforms with and without spin precession will drop to less than 0.97 after a
few days ($\sim 10^{5}~{}{\rm s}$). We find that spin-orbit coupling has the
potential to improve the accuracy of parameter estimation, especially for the
binary inclination angle and spin precession cone opening angle, by up to 3
orders of magnitude.

##  I Introduction

Rapidly spinning neutron stars (NSs) are promising sources of a long-lasting
form of gravitational waves (GWs), namely continuous waves (CWs) [1, 2, 3, 4,
5]. Detecting these potential CWs would help us solve some of the mysteries in
NS physics, such as NS’s equation of state, deformability, and magnetic field
[6, 7, 8].

There are two types of simplified waveforms that are commonly used in current
searches for CWs emitted by NSs (modeled as Newtonian rigid bodies) with
Advanced LIGO [9] and Advanced Virgo [10]. One is the mass quadrupole mode
with a frequency at twice the rotation frequency of the pulsar, which comes
from a triaxial rigid body rotating about one of its principal axes with
assumed principal moments of inertia $I_{1}<I_{2}<I_{3}$ (hereafter referred
to as the triaxial aligned waveforms, e.g., [11, 12, 13, 14]); the other is
the mode with frequencies at both once and twice the rotation frequency, which
comes from an axisymmetric freely precessing rigid body with assumed
$I_{1}=I_{2}\neq I_{3}$ (hereafter referred to as the biaxial waveforms, e.g.,
[15, 16, 17]). Similar two-frequency mode searches are also performed in [18,
19, 20, 21, 22, 23, 24]. The GW emission due to $r$-modes [25, 26] in a
rotating perfect fluid star is not the subject of our work. So far, no
credible detection has been reported in these searches [5].

The more general waveforms that come from a freely precessing triaxial rigid
body (hereafter referred to as the triaxial non-aligned waveforms) were first
calculated by Zimmermann [27]. The dominant waveform components are obtained
by expanding the quadrupole moment formula in terms of small parameters, such
as wobble angle, oblateness, and non-axisymmetry parameters. These waveforms
have been extended to include higher than first-order expansion terms of the
wobble angle and non-axisymmetry [28, 29] in order to extract more physical
information.

In addition to the waveform modeling of isolated NS discussed above, there are
also considerations about the NS located in a binary system since the
electromagnetic observations show that nearly half of the known pulsars within
the most sensitive band of the ground-based GW detectors belong to binary
systems [30, 31, 32]. Some search schemes are proposed for this type of CWs
[33, 34, 30, 35, 36]. However, the waveform model used in these searches is
obtained by simply incorporating the Doppler frequency modulation into the
phase of the triaxial aligned waveforms emitted by the isolated NS.

For future space-borne GW detectors, the detectability and parameter
estimation accuracy of double NS systems that will merge within the next 10
Myr have been studied for LISA [37] and TianQin [38]. Based on the merger rate
density ($920~{}\mathrm{Gpc^{-3}~{}yr^{-1}}$) inferred from GWTC-1 [39], about
300 double NS systems are expected to be detected in the mHz band during the
4-year observation period, including binaries with orbital periods shorter
than 10 minutes. Proposed next-generation ground-based GW observatories, such
as Cosmic Explorer [40] and Einstein Telescope [41], are expected to operate
concurrently with LISA and TianQin in the 2030s. Searching for CWs from the
spinning NS in such a tight system identified by LISA and/or TianQin requires
consideration of the interaction of the rapidly spinning NS with its
companion.

In this paper, we incorporate the effects of spin-orbit coupling to the GWs
emitted by the spinning NS in a circular orbital binary and extend the
triaxial aligned NS to a general triaxial non-aligned NS. Other effects, such
as magnetic dipole field [42], tidal interaction [43], and radiation reaction
[44] are neglected in the current work. Spin-orbit coupling causes spin
precession and orbital precession around the total angular momentum [44].
Similar to the treatment in previous works, such as that of LIGO [33], the GWs
emitted by the spinning NS in a binary are obtained by simply incorporating
the Doppler frequency modulation (including the effects of orbital precession)
into the phase of the triaxial non-aligned waveforms of the NS with spin
precession. In contrast to the isolated case (neglecting the electromagnetic
and the gravitational radiation-reaction torque as in [27]), the spin angular
momentum of the NS is no longer constant in the binary. Instead, it will be
precessed due to the spin-orbit coupling. We analytically solve the spin
precession equation for the NS using the perturbation method to obtain the
spin angular frequency evolution and calculate the waveforms based on the
quadrupole moment formula. Next, the waveforms are expanded into some simple
components in the small parameter case for the subsequent analysis of CW
detection. Finally, using these easy-to-use waveform components, we
investigate the impact of spin-orbit coupling on the parameter estimation
accuracy of the spinning NS. Calculations along these lines yield the
following results: (i) The waveforms of the NS undergoing spin precession will
deviate from the isolated ones after a few days ($\sim 10^{5}~{}{\rm s}$) when
the fitting factor between the two waveforms drops to less than 0.97. (ii)
Spin-orbit coupling has the potential to improve the parameter estimation
accuracy, specifically for the cosine of the binary inclination $\cos\iota$
and spin precession cone opening angle $\theta_{S}$, by up to 3 orders of
magnitude.

The rest of this paper is organized as follows. In Sec. II, we briefly review
the mathematical formalism for GWs from an isolated spinning NS, which will be
used for subsequent calculations for NS in a binary system. Analytical
approximations for the GWs from a spinning NS in a binary system, taking into
account spin-orbit coupling effects, are given in Sec. III. The comparison of
results derived using waveforms with and without spin precession is given in
Sec. IV. The parameter estimation accuracy of the waveforms with and without
spin-orbit coupling using Cosmic Explorer [40] are given in Sec. V. Our
conclusions are discussed in Sec. VI. Some details of our calculation have
been relegated to the appendix in order to keep the main ideas of the paper as
clear as possible.

##  II Gravitational waveforms from isolated NS

Since the waveforms emitted by spinning NS undergoing spin-orbit coupling are
based on the waveforms emitted by the isolated NS, we will first discuss the
case for the isolated NS. Following the conventions of Landau and Lifshitz
[45] and Zimmermann [27], in Fig. 1, the inertial coordinate system is denoted
as $(X,Y,Z)$ with basis vectors
$(\boldsymbol{e}_{x},\boldsymbol{e}_{y},\boldsymbol{e}_{z})$ and
$\boldsymbol{e}_{z}$ along the body’s angular momentum, and the body
coordinate system $(x_{1},x_{2},x_{3})$ with basis vectors
$(\boldsymbol{e}_{1},\boldsymbol{e}_{2},\boldsymbol{e}_{3})$ parallel to the
eigenvectors of the body’s moment of inertia tensor and satisfying
$I_{3}>I_{2}\geq I_{1}$. The origins of the two systems are placed at the
center of mass of the NS. The Euler angles ($\theta,\phi,\psi$) describe the
orientation of the body coordinate system with respect to the inertial
coordinate system. We use the Latin subscripts (e.g., $x,y,z$) for components
evaluated in the inertial coordinate system, and the Greek ones (e.g.,
$\mu,\nu$) in the body coordinate system.

Figure 1: The inertial coordinate system $(X,Y,Z)$ and the body coordinate
system $(x_{1},x_{2},x_{3})$ are both centered on the center of mass ($O_{1}$)
of the isolated NS. The $x_{1}-x_{2}$ plane intersects the $X-Y$ plane at the
line of nodes $O_{1}N$. ($\theta,\phi,\psi$) are Euler angles. The distant
observer with position vector $\boldsymbol{D}$ is assumed in the $Y-Z$ plane
at colatitude $i$ from the $Z$ axis.

The metric perturbation under the transverse-traceless gauge can be written in
terms of two GW polarizations,
$h_{jk}^{\mathrm{TT}}=h_{+}\left(\hat{e}_{+}\right)_{jk}+h_{\times}\left(\hat{e}_{\times}\right)_{jk}$,
with the polarization tensors defined as

$\hat{e}_{+}\equiv\hat{v}\otimes\hat{v}-\hat{w}\otimes\hat{w}\,,\quad\hat{e}_{\times}\equiv\hat{v}\otimes\hat{w}+\hat{w}\otimes\hat{v}\,,$
(1)

where $\hat{v}$ and $\hat{w}$ are the transverse basis vectors perpendicular
to the wave’s propagation direction, and $\otimes$ denotes the tensor product.
Without loss of generality, we assume that the observer is located in the
$Y-Z$ plane with colatitude $i$ from the $Z$ axis and distance
$D=|\boldsymbol{D}|$. In this configuration,

$\hat{v}\equiv\hat{e}_{y}\cos i-\hat{e}_{z}\sin
i\,,\quad\hat{w}\equiv-\hat{e}_{x}\,,$ (2)

and the two GW polarizations can be written as [27]

$\displaystyle h_{+}$ $\displaystyle=-\frac{G}{c^{4}D}[(R_{y\mu}\cos
i-R_{z\mu}\sin i)$ (3a) $\displaystyle\times(R_{y\nu}\cos i-R_{z\nu}\sin
i)-R_{x\mu}R_{x\nu}]A_{\mu\nu}\,,$ $\displaystyle h_{\times}$
$\displaystyle=\frac{2G}{c^{4}D}(R_{y\mu}\cos i-R_{z\mu}\sin
i)R_{x\nu}A_{\mu\nu}\,,$ (3b)

where Einstein summation is performed for $\mu$ and $\nu$, both of which take
the values $\\{1,2,3\\}$. The components of the symmetric matrix $A_{\mu\nu}$
are given by [27]

$\displaystyle A_{11}$
$\displaystyle=2\left(\Delta_{2}\Omega_{2}^{2}-\Delta_{3}\Omega_{3}^{2}\right)\,,$
(4a) $\displaystyle A_{22}$
$\displaystyle=2\left(\Delta_{3}\Omega_{3}^{2}-\Delta_{1}\Omega_{1}^{2}\right)\,,$
(4b) $\displaystyle A_{33}$
$\displaystyle=2\left(\Delta_{1}\Omega_{1}^{2}-\Delta_{2}\Omega_{2}^{2}\right)\,,$
(4c) $\displaystyle A_{12}$
$\displaystyle=A_{21}=\left(\Delta_{1}-\Delta_{2}\right)\Omega_{1}\Omega_{2}+\Delta_{3}\dot{\Omega}_{3}\,,$
(4d) $\displaystyle A_{23}$
$\displaystyle=A_{32}=\left(\Delta_{2}-\Delta_{3}\right)\Omega_{2}\Omega_{3}+\Delta_{1}\dot{\Omega}_{1}\,,$
(4e) $\displaystyle A_{31}$
$\displaystyle=A_{13}=\left(\Delta_{3}-\Delta_{1}\right)\Omega_{3}\Omega_{1}+\Delta_{2}\dot{\Omega}_{2}\,,$
(4f)

with $\Delta_{1}\equiv I_{2}-I_{3}$, $\Delta_{2}\equiv I_{3}-I_{1}$,
$\Delta_{3}\equiv I_{1}-I_{2}$. Here, $(\Omega_{1},\Omega_{2},\Omega_{3})$
denotes the angular frequency of NS in the body coordinate system. The
rotation matrix that transforms from the body coordinate system to the
inertial coordinate system (e.g., $R_{x2}$ denotes the entry in row $x=1$ and
column $\mu=2$) is

$\displaystyle
R=\left(\begin{array}[]{ccc}\cos\psi\cos\phi-\cos\theta\sin\psi\sin\phi&-\cos\theta\cos\psi\sin\phi-\sin\psi\cos\phi&\sin\theta\sin\phi\\\
\cos\theta\sin\psi\cos\phi+\cos\psi\sin\phi&\cos\theta\cos\psi\cos\phi-\sin\psi\sin\phi&-\sin\theta\cos\phi\\\
\sin\theta\sin\psi&\sin\theta\cos\psi&\cos\theta\\\ \end{array}\right)\,.$ (8)

According to Euler’s equations of free rotation of a rigid body and the
initial conditions ($\Omega_{1}(0)=a,\Omega_{2}(0)=0,\Omega_{3}(0)=b$), the
angular frequencies in the body coordinate system are [27]

$\displaystyle\Omega_{1}$ $\displaystyle=a\operatorname{cn}(\tau,m)\,,$ (9a)
$\displaystyle\Omega_{2}$
$\displaystyle=a\left[\frac{I_{1}\left(I_{3}-I_{1}\right)}{I_{2}\left(I_{3}-I_{2}\right)}\right]^{1/2}\operatorname{sn}(\tau,m)\,,$
(9b) $\displaystyle\Omega_{3}$ $\displaystyle=b\operatorname{dn}(\tau,m)\,,$
(9c)

where $\rm cn$, $\rm sn$ and $\rm dn$ are Jacobian elliptic functions [46]
with the parameters

$\displaystyle\tau$
$\displaystyle=bt\left[\frac{\left(I_{3}-I_{2}\right)\left(I_{3}-I_{1}\right)}{I_{1}I_{2}}\right]^{1/2},$
(10) $\displaystyle m$
$\displaystyle=\frac{\left(I_{2}-I_{1}\right)I_{1}a^{2}}{\left(I_{3}-I_{2}\right)I_{3}b^{2}}\,.$
(11)

The Euler angles can be expressed in terms of Jacobian elliptic functions and
the fourth theta functions [46] ($\vartheta_{4}$ and their derivatives
$\vartheta_{4}^{\prime}$):

$\displaystyle\cos\theta$
$\displaystyle=\frac{I_{3}b}{S}\operatorname{dn}(\tau,m)\,,$ (12a)
$\displaystyle\tan\psi$
$\displaystyle=\left[\frac{I_{1}\left(I_{3}-I_{2}\right)}{I_{2}\left(I_{3}-I_{1}\right)}\right]^{1/2}\frac{\operatorname{cn}(\tau,m)}{\operatorname{sn}(\tau,m)}\,,$
(12b) $\displaystyle\phi$ $\displaystyle=\phi_{1}+\phi_{2}\,,$ (12c)

with $S$ being the magnitude of the spin angular momentum of the NS and
$\phi_{1,2}$ given by [29]

$\displaystyle\exp\left[2\mathrm{i}\phi_{1}(t)\right]=\frac{\vartheta_{4}\left(\frac{2\pi
t}{T}+\mathrm{i}\pi\alpha,q\right)}{\vartheta_{4}\left(\frac{2\pi
t}{T}-\mathrm{i}\pi\alpha,q\right)}\,,$ (13a)
$\displaystyle\phi_{2}=\frac{2\pi
t}{T^{\prime}}=\left[\frac{S}{I_{1}}+\frac{2\pi\mathrm{i}}{T}\frac{\vartheta_{4}^{\prime}(\mathrm{i}\pi\alpha,q)}{\vartheta_{4}(\mathrm{i}\pi\alpha,q)}\right]t\,.$
(13b)

Here, $T$ is the period of the angular frequency in the body coordinate
system,

$T=\frac{4K(m)}{b}\left[\frac{I_{1}I_{2}}{\left(I_{3}-I_{2}\right)\left(I_{3}-I_{1}\right)}\right]^{1/2}\,.$
(14)

$\alpha$ satisfies ${\rm sn}[2\mathrm{i}\alpha
K(m)]=\mathrm{i}I_{3}b/(I_{1}a)$. $\mathrm{i}$ is the imaginary unit.
$q=\exp[-\pi K(1-m)/K(m)]$, where $K(m)$ is the complete elliptic integral of
the first kind [46]. Since $\cos\phi_{2}$ has a period $T^{\prime}$ which is
generally not commensurate with $T$, the motion of the NS is usually
nonperiodic. When the NS becomes axisymmetric, $T^{\prime}\to 2\pi I_{1}/S$
[27].

Explicit waveforms are usually expressed in a series expansion of some small
parameters [27, 28, 29]. To facilitate the following calculation, we define
the spinning NS’s free precession angular frequency and rotation angular
frequency

$\displaystyle\Omega_{\rm p}\equiv\frac{2\pi}{T}\,,\quad\Omega_{\rm
r}\equiv\frac{2\pi}{T^{\prime}}-\frac{2\pi}{T}\,,$ (15)

and three parameters that characterize NS’s properties

$\displaystyle\epsilon\equiv\frac{I_{3}-I_{1}}{I_{3}}\,,\quad\kappa\equiv\frac{1}{16}\frac{I_{3}}{I_{1}}\frac{I_{2}-I_{1}}{I_{3}-I_{2}}\,,\quad\gamma\equiv\frac{aI_{1}}{bI_{3}}\,,$
(16)

where $\epsilon$ is called the oblateness (or poloidal ellipticity [5]) of the
NS, $\kappa$ describes the ($I_{2}-I_{1}$) with respect to the axisymmetric
non-sphericity ($I_{3}-I_{2}$), while $\gamma$ is called the wobble angle.
Their characteristic values are discussed in [28]. For the small quantities
above, the expansions of the sines and cosines of the Euler angles (cf. Eqs.
(12)) in $S$-aligned coordinate system up to terms of $O(\gamma^{2})$ and
$O(\kappa^{2})$ are [29]

$\displaystyle\cos{\theta}$
$\displaystyle=1-{\gamma^{2}}/{2}\,,\quad\sin{\theta}=\gamma+8\gamma\kappa\sin^{2}(t\Omega_{\rm
p})\,,$ (17a) $\displaystyle\cos{\phi}$ $\displaystyle=\cos[t(\Omega_{\rm
r}+\Omega_{\rm p})]\,,\,\,\sin{\phi}=\sin[t(\Omega_{\rm r}+\Omega_{\rm
p})]\,,$ (17b) $\displaystyle\cos{\psi}$ $\displaystyle=\sin\left(t\Omega_{\rm
p}\right)+8\kappa\sin\left(t\Omega_{\rm p}\right)\cos^{2}\left(t\Omega_{\rm
p}\right)$ (17c) $\displaystyle+8\kappa^{2}\left(3\sin\left(3t\Omega_{\rm
p}\right)-13\sin\left(t\Omega_{\rm p}\right)\right)\cos^{2}\left(t\Omega_{\rm
p}\right)\,,$ $\displaystyle\sin{\psi}$ $\displaystyle=\cos\left(t\Omega_{\rm
p}\right)-8\kappa\sin^{2}\left(t\Omega_{\rm p}\right)\cos\left(t\Omega_{\rm
p}\right)$ (17d) $\displaystyle+96\kappa^{2}\sin^{4}\left(t\Omega_{\rm
p}\right)\cos\left(t\Omega_{\rm p}\right)\,.$

##  III Gravitational waveforms from spinning NS in a binary

Following the treatment in previous works, such as [33], the GWs emitted by
the spinning NS in a binary can be obtained by incorporating the Doppler
frequency modulation (modulated by orbital precession) into the phase of the
triaxial non-aligned waveforms of the NS with spin precession.

First, we calculate the GWs emitted by a spinning NS undergoing spin
precession. Consider a binary system consisting of a spinning NS with spin
angular momentum $\boldsymbol{S}$ and a nonspinning NS (or a slowly spinning
NS of which the spin effects can be ignored). This is consistent with the
standard evolution scenario of the double NS formed in an isolated system, in
which one of the NSs is a rapidly spinning millisecond pulsar and the other is
a normal pulsar [47, 48]. If the companion star has spin, then spin-spin
coupling will also cause the rapidly spinning NS to precess. The ratio of the
precessional angular frequency due to the spin-spin (SS) coupling to that due
to the spin-orbit (SO) coupling satisfies $\Omega_{\rm
pre}({\rm{SS}})/\Omega_{\rm
pre}({\rm{SO}})<4M/(3M+m_{1})(m_{2}/M)^{1/2}(R_{2}/r)^{1/2}$, where
$M=m_{1}+m_{2}$ is the total mass of the binary, $m_{1}$ is the mass of the
rapidly spinning NS, $m_{2}$ and $R_{2}$ are the mass and the radius of the
companion, and $r$ is the orbital separation (cf. Eq. (10.179) in [49]). For a
double NS system with an orbital period of 10 min, $\Omega_{\rm
pre}({\rm{SS}})/\Omega_{\rm pre}({\rm{SO}})<0.01$ for a companion with maximum
spin, i.e., the spin-spin precession angular frequency is at least two orders
of magnitude smaller than the spin-orbit precession angular frequency,
therefore we can ignore the spin-spin coupling. Furthermore, in this work we
assume that the orbits of binary stars are circular. A recent population
synthesis simulation in [50] shows that the eccentricities satisfy $e\lesssim
0.01$ for the double NS systems with orbital periods shorter than 10 min. On
the other hand, in the extreme cases where eccentricity is important, we need
to generalize our current work to incorporate the effects of eccentricity.
This can be a subject of our future work.

For a binary system in which only one of the bodies has spin, the precession
equations [44] for the spin of the body and the orbit of the binary show that
to a reasonable approximation the total angular momentum $\boldsymbol{J}$
maintains its direction, $\boldsymbol{S}$ keeps its magnitude constant and
precesses around $\boldsymbol{J}$ with

$\displaystyle\frac{d\boldsymbol{S}}{dt}=\boldsymbol{\Omega}_{\rm
pre}\times\boldsymbol{S}\,,$ (18)

where

$\boldsymbol{\Omega}_{\rm
pre}=\frac{G}{2c^{2}r^{3}}\left(1+\frac{3M}{m_{1}}\right)\boldsymbol{J}$ (19)

is the angular frequency of spin precession and orbital precession induced by
spin-orbit coupling. Although the decreasing $r$ and the magnitude of
$\boldsymbol{J}$ due to the radiation reaction cause the magnitude
${\Omega}_{\rm pre}$ of $\boldsymbol{\Omega}_{\rm pre}$ to vary with time, for
a typical double NS system with a merger time $\sim O(10^{4}{\rm yr})$ and a
half-year observation time (see Sec. IV), the relative variation of
${\Omega}_{\rm pre}$ is $\lesssim 10^{-5}$, so we can assume that
${\Omega}_{\rm pre}$ remains approximately constant in the case considered.

Figure 2: $J$-aligned coordinate system $(X_{J},Y_{J},Z_{J})$ with $Z_{J}$
axis parallel to the total angular momentum $\boldsymbol{J}$ and $S$-aligned
coordinate system $(X_{S},Y_{S},Z_{S})$ with $Z_{S}$ axis aligned with spin
angular momentum $\boldsymbol{S}$ are both centered at the center of mass
($O_{1}$) of the spinning NS, where spin-orbit coupling causes
$\boldsymbol{S}$ to precess about $\boldsymbol{J}$ at an angular frequency
${\Omega}_{\rm pre}$. The opening angle of $\boldsymbol{S}$ precession cone is
$\theta_{S}$. Suppose that at initial time, $X_{S}$ axis coincides with
$X_{J}$ axis, and $\boldsymbol{S}$ is in the $Y_{J}-Z_{J}$ plane. The
projection of the spin evolution $\boldsymbol{S}(t)$ coincides with $O_{1}N$
in the $X_{J}-Y_{J}$ plane, its precession angle $\alpha=\Omega_{\rm pre}t$.
The distant observer with position vector $\boldsymbol{D}$ is also assumed in
the $Y_{J}-Z_{J}$ plane at colatitude $\iota$ from $\boldsymbol{J}$.

As shown in Fig. 2, $(X_{J},Y_{J},Z_{J})$ is the coordinate system with the
origin placed at the center of mass ($O_{1}$) of the spinning NS and $Z_{J}$
axis parallel to $\boldsymbol{J}$ (hereafter referred to as the $J$-aligned
coordinate system), in which the distant observer is assumed to be in the
$Y_{J}-Z_{J}$ plane with the inclination $\iota$ and the position vector
$\boldsymbol{D}$. The opening angle of $\boldsymbol{S}$ precession cone is
$\theta_{S}$. The coordinate system $(X_{S},Y_{S},Z_{S})$ constructed with
$Z_{S}$ axis aligned with $\boldsymbol{S}$ is referred to as the $S$-aligned
coordinate system. Without loss of generality, we assume that at initial time,
the $X_{S}$ axis coincides with the $X_{J}$ axis, and $\boldsymbol{S}$ is in
the $Y_{J}-Z_{J}$ plane. The evolution of $\boldsymbol{S}$ over time is
represented by $\boldsymbol{S}(t)$ with the precession angle
$\alpha=\Omega_{\rm pre}t$ measured in the $X_{J}-Y_{J}$ plane.

To simplify the calculation of waveforms, we use a similar convention for two
polarization tensors (cf. Eq. (1) and Eq. (2)) as in the calculation of the
isolated NS in Sec. II, so that the waveforms in Eqs. (3) also apply to the
triaxial non-aligned NS in a binary, but with ($i,R,A$) in Eqs. (3) replaced
by ($\iota,\mathcal{R},\mathcal{A}$), which are the quantities calculated in a
coordinate system at rest with respect to the center of mass of the spinning
NS, i.e., the $J$-aligned coordinate system shown in Fig 2.

$\displaystyle i\rightarrow\iota,\quad R\rightarrow\mathcal{R},\quad
A\rightarrow\mathcal{A}\,.$ (20)

According to the discussion in [28], the small parameters in Eq. (16) satisfy
$\epsilon\ll\kappa\ll\gamma_{\rm max}$ and $\kappa\sim O(\gamma^{2})$. Our
goal is to expand the CWs of the spinning NS in a tight binary up to order
$O(\gamma^{2})$ and $O(\kappa)$. We first calculate $\mathcal{R}$ which
represents the rotation matrix from the body coordinate system of the spinning
NS to $J$-aligned coordinate system. It can be obtained by the following
rotation transformations: first, from the body coordinate system to the
$S$-aligned coordinate system by $R$ (cf. Eq. (8)), and then from the
$S$-aligned coordinate system to the $J$-aligned coordinate system by
$\mathcal{T}_{S\rightarrow J}$. Thus,

$\mathcal{R}=\mathcal{T}_{S\rightarrow J}\cdot R\,,$ (21)

where

$\displaystyle\mathcal{T}_{S\rightarrow J}=$ (22)
$\displaystyle\left(\begin{array}[]{ccc}\cos(\Omega_{\rm
pre}t)&-\sin(\Omega_{\rm pre}t)&0\\\ \sin(\Omega_{\rm pre}t)&\cos(\Omega_{\rm
pre}t)&0\\\ 0&0&1\\\ \end{array}\right)\left(\begin{array}[]{ccc}1&0&0\\\
0&\cos\theta_{S}&\sin\theta_{S}\\\ 0&-\sin\theta_{S}&\cos\theta_{S}\\\
\end{array}\right)\,.$ (29)

After inserting Eqs. (17) into Eq. (21), $\mathcal{R}$ can be expanded as Eqs.
(55) given in Appendix A. Although $\mathcal{R}$ contains $O(\kappa^{2})$
terms as in [29], we keep them in Eqs. (55) in order to facilitate the
generalization of the waveforms to higher orders in the future.

Next, we need to calculate $\mathcal{A}$ which depends on the angular
frequencies $\omega_{i}$ of the spinning NS under spin-orbit coupling. The
$(1,1)$ and $(1,2)$ components of $\mathcal{A}$ read

$\displaystyle\mathcal{A}_{11}$
$\displaystyle=2\left(\Delta_{2}\omega_{2}^{2}-\Delta_{3}\omega_{3}^{2}\right)\,,$
(30a) $\displaystyle\mathcal{A}_{12}$
$\displaystyle=\mathcal{A}_{21}=\left(\Delta_{1}-\Delta_{2}\right)\omega_{1}\omega_{2}+\Delta_{3}\dot{\omega}_{3}\,.$
(30b)

The other components of $\mathcal{A}$ are the same as the corresponding ones
of $A$ in Eqs. (4) but with $\Omega_{i}$ replaced by $\omega_{i}$. Since in
the body coordinate system
$\boldsymbol{S}=S_{\mu}\boldsymbol{e}_{\mu}=I_{1}\omega_{1}\boldsymbol{e}_{1}+I_{2}\omega_{2}\boldsymbol{e}_{2}+I_{3}\omega_{3}\boldsymbol{e}_{3}$
and $\boldsymbol{\Omega}_{\rm pre}={\Omega}_{\rm
pre}\boldsymbol{e}_{z}={\Omega}_{\rm
pre}\mathcal{R}_{z\mu}\boldsymbol{e}_{\mu}$, then Eq. (18) can be expressed as
follows

$\frac{d{S_{\mu}}}{dt}\boldsymbol{e}_{\mu}+S_{\mu}{\omega}_{\nu}\boldsymbol{e}_{\nu}\times\boldsymbol{e}_{\mu}={\Omega}_{\rm
pre}\mathcal{R}_{z\mu}S_{\nu}\boldsymbol{e}_{\mu}\times\boldsymbol{e}_{\nu}\,,$
(31)

with components

$\displaystyle\frac{\mathrm{d}\omega_{1}}{\mathrm{d}t}$
$\displaystyle=\frac{\Delta_{1}}{I_{1}}\omega_{2}\omega_{3}+\frac{\Omega_{\rm
pre}}{I_{1}}\left(\mathcal{R}_{z2}I_{3}\omega_{3}-\mathcal{R}_{z3}I_{2}\omega_{2}\right)\,,$
(32a) $\displaystyle\frac{\mathrm{d}\omega_{2}}{\mathrm{d}t}$
$\displaystyle=\frac{\Delta_{2}}{I_{2}}\omega_{3}\omega_{1}+\frac{\Omega_{\rm
pre}}{I_{2}}\left(\mathcal{R}_{z3}I_{1}\omega_{1}-\mathcal{R}_{z1}I_{3}\omega_{3}\right)\,,$
(32b) $\displaystyle\frac{\mathrm{d}\omega_{3}}{\mathrm{d}t}$
$\displaystyle=\frac{\Delta_{3}}{I_{3}}\omega_{1}\omega_{2}+\frac{\Omega_{\rm
pre}}{I_{3}}\left(\mathcal{R}_{z1}I_{2}\omega_{2}-\mathcal{R}_{z2}I_{1}\omega_{1}\right)\,.$
(32c)

The second terms on the right hand sides are the spin-orbit coupling terms.

Although the general solution of the above precession Eqs. (32) can be
obtained numerically, the analytic solution is more favorable for GW detection
and parameter estimation because it can be incorporated directly into search
algorithms and is more manageable and efficient in data analysis. We expect
the difference between the angular frequency of the spinning NS in an isolated
case and that in a binary system to be small (see Fig. 3), since the spin-
orbit coupling is of 1.5 post-Newtonian order [49]. In this sense, we use the
perturbation method to solve Eqs. (32) analytically by assuming

$\displaystyle\omega_{i}\approx\Omega_{i}+\delta\Omega_{i}\quad(i=1,2,3)\,.$
(33)

Inserting the above expression into Eqs. (32) yields the linearized evolution
equations

$\left(\begin{array}[]{c}\frac{d\delta\Omega_{1}}{dt}\\\\[5.69054pt]
\frac{d\delta\Omega_{2}}{dt}\\\\[5.69054pt]
\frac{d\delta\Omega_{3}}{dt}\end{array}\right)=\left(\begin{array}[]{ccc}0&\frac{\Delta_{1}\Omega_{3}-\Omega_{\rm
pre}\mathcal{R}_{z3}I_{2}}{I_{1}}&\frac{\Delta_{1}\Omega_{2}+\Omega_{\rm
pre}\mathcal{R}_{z2}I_{3}}{I_{1}}\\\\[5.69054pt]
\frac{\Delta_{2}\Omega_{3}+\Omega_{\rm
pre}\mathcal{R}_{z3}I_{1}}{I_{2}}&0&\frac{\Delta_{2}\Omega_{1}-\Omega_{\rm
pre}\mathcal{R}_{z1}I_{3}}{I_{2}}\\\\[5.69054pt]
\frac{\Delta_{3}\Omega_{2}-\Omega_{\rm
pre}\mathcal{R}_{z2}I_{1}}{I_{3}}&\frac{\Delta_{3}\Omega_{1}+\Omega_{\rm
pre}\mathcal{R}_{z1}I_{2}}{I_{3}}&0\end{array}\right)\left(\begin{array}[]{c}\delta\Omega_{1}\\\\[2.84526pt]
\delta\Omega_{2}\\\\[2.84526pt]
\delta\Omega_{3}\end{array}\right)+\left(\begin{array}[]{c}\frac{\Omega_{\rm
pre}}{I_{1}}\left(\mathcal{R}_{z2}I_{3}\Omega_{3}-\mathcal{R}_{z3}I_{2}\Omega_{2}\right)\\\\[5.69054pt]
\frac{\Omega_{\rm
pre}}{I_{2}}\left(\mathcal{R}_{z3}I_{1}\Omega_{1}-\mathcal{R}_{z1}I_{3}\Omega_{3}\right)\\\\[5.69054pt]
\frac{\Omega_{\rm
pre}}{I_{3}}\left(\mathcal{R}_{z1}I_{2}\Omega_{2}-\mathcal{R}_{z2}I_{1}\Omega_{1}\right)\end{array}\right)\,.$
(34)

According to Appendix C, $\delta\Omega_{3}\ll\delta\Omega_{1,2}$, so the
contribution of $\delta\Omega_{3}$ can be ignored when solving for
$\delta\Omega_{1,2}$. For a typical double NS system with 6 min orbital period
discussed below, $\Omega_{\rm pre}/{\rm
Hz}\lesssim\epsilon\ll\kappa\lesssim\gamma^{2}$. We first expand $\Omega_{i}$
(cf. Eqs. (9)) to leading order. The combination of Eq. (10), Eq. (14), and
Eq. (15) gives $\tau=(2K(m)/{\pi})t\Omega_{\rm p}=(1+O(m))t\Omega_{\rm p}$.
From Eq. (11) and Eq. (16), we have $m=16\gamma^{2}\kappa$. Therefore,
$\tau\simeq t\Omega_{\rm p}$ to leading order. In a similar way, to leading
order, the Jacobian elliptic functions become ${\rm
cn}(\tau,m)\simeq\cos(t\Omega_{\rm p})$, ${\rm
sn}(\tau,m)\simeq\sin(t\Omega_{\rm p})$, and ${\rm dn}(\tau,m)\simeq 1$, and
the factor
$\left[{I_{1}\left(I_{3}-I_{1}\right)}/{I_{2}/\left(I_{3}-I_{2}\right)}\right]^{1/2}\simeq
1$. Consequently, we have the leading order result

$\displaystyle\Omega_{1}$ $\displaystyle=a\cos(t\Omega_{\rm p})\,,$ (35a)
$\displaystyle\Omega_{2}$ $\displaystyle=a\sin(t\Omega_{\rm p})\,,$ (35b)
$\displaystyle\Omega_{3}$ $\displaystyle=b\,.$ (35c)

By inserting Eqs. (35) and Eqs. (55) into Eq. (34), we expand the linearized
evolution equations up to order $O(\gamma^{2})$ and $O(\kappa)$, and ignore
higher order terms such as $O(\gamma^{2}\kappa)$, $O(\gamma\Omega_{\rm pre})$
and $O(\kappa\Omega_{\rm pre})$. To give analytic approximations for
$\delta\Omega_{i}$, we further simplify the equations by setting
$1+\epsilon\simeq 1$ and discard $O(\Omega_{\rm pre})$ in $O(\Omega_{\rm
pre})+O(b\epsilon)$ due to $O(b)\sim{2\pi}\times 100~{}\rm{Hz}$. Finally, we
obtain

$\displaystyle\frac{d\delta\Omega_{1}}{dt}$ $\displaystyle\simeq
b\epsilon(-1+16\kappa)\delta\Omega_{2}+b\Omega_{\rm pre}\sin\theta_{\rm
S}\sin(t\Omega_{\rm r})\,,$ (36a) $\displaystyle\frac{d\delta\Omega_{2}}{dt}$
$\displaystyle\simeq b\epsilon\delta\Omega_{1}+b\Omega_{\rm
pre}\sin\theta_{\rm S}\cos(t\Omega_{\rm r})\,.$ (36b)

Setting the initial conditions
$(\delta\Omega_{1}(0),\delta\Omega_{2}(0),\delta\Omega_{3}(0))=(0,0,0)$, and
ignoring $O(b\epsilon)$ in $O(b\epsilon)+O(\Omega_{\rm r})$ due to
$b\simeq\Omega_{\rm r}$, the integration of Eqs. (36) yields

$\displaystyle\delta\Omega_{1}$ $\displaystyle\simeq{\Omega_{\rm
pre}\sin\theta_{\rm
S}\left[\cos(tb\epsilon\sqrt{1-16\kappa})-\cos(t\Omega_{\rm r})\right]}\,,$
(37a) $\displaystyle\delta\Omega_{2}$ $\displaystyle\simeq{\Omega_{\rm
pre}\sin\theta_{\rm
S}\left[\sin(tb\epsilon\sqrt{1-16\kappa})+\sin(t\Omega_{\rm r})\right]}\,.$
(37b)

Using Eqs. (37), the integration of the third row in Eq. (34) leads to

$\delta\Omega_{3}\simeq{\gamma\Omega_{\rm
pre}\sin{\theta_{S}}[\cos(t(\Omega_{\rm p}+\Omega_{\rm r}))-1]}\,.$ (38)

The error of this analytic approximation is shown and discussed in Appendix C.

Now, $\mathcal{A}$ can be calculated by inserting $\Omega_{i}$ (cf. Eqs. (35))
and $\delta\Omega_{i}$ (cf. Eqs. (37) and Eq. (38)) into Eqs. (30). The result
is given in Eqs. (54) of Appendix A.

With $\mathcal{R}$ and $\mathcal{A}$, the waveforms emitted by the spinning NS
undergoing spin precession can be expressed in terms of the series expansion
as follows

$\displaystyle h_{+}$
$\displaystyle=h_{+}^{(1)}+h_{+}^{(2)}+h_{+}^{(3)}+\cdots$ (39a)
$\displaystyle h_{\times}$
$\displaystyle=h_{\times}^{(1)}+h_{\times}^{(2)}+h_{\times}^{(3)}+\cdots$
(39b)

where

$\displaystyle h_{+}^{(1)}$
$\displaystyle=\frac{{G{b^{2}}{I_{3}}\epsilon\gamma}}{{4{c^{4}}D}}\Big{[}\cos(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))[\sin(2{\theta_{S}})(6{\sin^{2}}\iota-(3+\cos(2\iota))\cos(2t{\Omega_{{\rm{pre}}}}))+4\cos(2{\theta_{S}})\cos(t{\Omega_{{\rm{pre}}}})\sin(2\iota)]$
(40a)
$\displaystyle+2[-2\cos{\theta_{S}}\sin(2\iota)\sin(t{\Omega_{{\rm{pre}}}})+(3+\cos(2\iota))\sin{\theta_{S}}\sin(2t{\Omega_{{\rm{pre}}}})]\sin(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))\Big{]}\,,$
$\displaystyle h_{\times}^{(1)}$
$\displaystyle=\frac{{G{b^{2}}{I_{3}}\epsilon\gamma}}{{{c^{4}}D}}\Big{[}\cos(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))[2\sin\iota\cos(2{\theta_{S}})\sin(t{\Omega_{{\rm{pre}}}})-\cos\iota\sin(2{\theta_{S}})\sin(2t{\Omega_{{\rm{pre}}}})]$
(40b)
$\displaystyle+2[\sin\iota\cos{\theta_{S}}\cos(t{\Omega_{{\rm{pre}}}})-\cos\iota\cos(2t{\Omega_{{\rm{pre}}}})\sin{\theta_{S}}]\sin(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))\Big{]}\,,$
$\displaystyle h_{+}^{(2)}$
$\displaystyle=-\frac{{16G{b^{2}}{I_{3}}\epsilon\kappa}}{{{c^{4}}D}}\Big{[}(3+\cos(2\iota))[\cos^{4}\left({\frac{{{\theta_{S}}}}{2}}\right)\cos(2t({\Omega_{{\rm{pre}}}}+{\Omega_{\rm{r}}}))+\cos(2t({\Omega_{{\rm{pre}}}}-{\Omega_{\rm{r}}})){\sin^{4}}\left({\frac{{{\theta_{S}}}}{2}}\right)]$
(40c)
$\displaystyle+\cos(2t{\Omega_{\rm{r}}})[3{\sin^{2}}{\theta_{S}}{\sin^{2}}\iota+\cos(t{\Omega_{{\rm{pre}}}})\sin(2{\theta_{S}})\sin(2\iota)]-2\sin{\theta_{S}}\sin(2\iota)\sin(t{\Omega_{{\rm{pre}}}})\sin(2t{\Omega_{\rm{r}}})\Big{]}\,,$
$\displaystyle h_{\times}^{(2)}$
$\displaystyle=\frac{{16G{b^{2}}{I_{3}}\epsilon\kappa}}{{{c^{4}}D}}\Big{[}\cos(2t{\Omega_{\rm{r}}})[-2\sin\iota\sin(2{\theta_{S}})\sin(t{\Omega_{{\rm{pre}}}})-\cos\iota(3+\cos(2{\theta_{S}}))\sin(2t{\Omega_{{\rm{pre}}}})]$
(40d)
$\displaystyle-4[\cos\iota\cos{\theta_{S}}\cos(2t{\Omega_{{\rm{pre}}}})+\sin\iota\cos(t{\Omega_{{\rm{pre}}}})\sin{\theta_{S}}]\sin(2t{\Omega_{\rm{r}}})\Big{]}\,,$
$\displaystyle h_{+}^{(3)}$
$\displaystyle=\frac{Gb^{2}{I_{3}}\epsilon\gamma^{2}}{c^{4}D}\Big{[}(3+\cos(2\iota))[\cos^{4}\left({\frac{{{\theta_{S}}}}{2}}\right)\cos(2t({\Omega_{{\rm{pre}}}}+{\Omega_{\rm{r}}}+\Omega_{\rm{p}}))+\cos(2t({\Omega_{{\rm{pre}}}}-{\Omega_{\rm{r}}}-\Omega_{\rm{p}})){\sin^{4}}\left({\frac{{{\theta_{S}}}}{2}}\right)]$
(40e)
$\displaystyle+\cos(2t(\Omega_{\rm{r}}+\Omega_{\rm{p}}))[3{\sin^{2}}{\theta_{S}}{\sin^{2}}\iota+\cos(t{\Omega_{{\rm{pre}}}})\sin(2{\theta_{S}})\sin(2\iota)]-2\sin{\theta_{S}}\sin(2\iota)\sin(t{\Omega_{{\rm{pre}}}})\sin(2t(\Omega_{\rm{r}}+\Omega_{\rm{p}}))\Big{]}\,,$
$\displaystyle h_{\times}^{(3)}$
$\displaystyle=\frac{Gb^{2}{I_{3}}\epsilon\gamma^{2}}{{{c^{4}}D}}\Big{[}\cos(2t(\Omega_{\rm{p}}+\Omega_{\rm{r}}))[2\sin(2{\theta_{S}})\sin\iota\sin(t{\Omega_{{\rm{pre}}}})+(3+\cos(2\theta_{S}))\cos\iota\sin(2t{\Omega_{{\rm{pre}}}})]$
(40f)
$\displaystyle+4[\cos{\theta_{S}}\cos(2t{\Omega_{{\rm{pre}}}})\cos\iota+\cos(t\Omega_{\rm{pre}})\sin{\theta_{S}}\sin{\iota}]\sin(2t(\Omega_{\rm{p}}+\Omega_{\rm{r}}))\Big{]}\,.$

Note that $\iota$ in the above expressions is the angle between
$\boldsymbol{J}$ and the line of sight, not to be confused with the
inclination defined by the angle between $\boldsymbol{S}$ and the line of
sight for the isolated case.

Here, $h_{+,\times}^{(1)}$ are of order $O(\gamma)$, $h_{+,\times}^{(2)}$ are
of order $O(\kappa)$, $h_{+,\times}^{(3)}$ are of order $O(\gamma^{2})$. Note
that only the waveform components up to terms of $O(\gamma)$, $O(\kappa)$ and
$O(\gamma^{2})$ are shown here, according to the discussion of the
characteristic value of $\gamma$ and $\kappa$ in [28]. The $O\left(\Omega_{\rm
pre}\right)$ components (see Appendix B), and higher-order $O(\gamma\kappa)$,
and $O(\kappa^{2})$ components can be ignored for the parameter values we
adopt in the following analysis.

Then, the waveforms of the spinning NS in a binary system with spin-orbit
coupling effects considered are completed by incorporating the Doppler
frequency modulation of this NS around the binary barycenter (BB) into the
phases of waveforms Eqs. (40), which is done by the second term on the right-
hand side of Eq. (49) in Sec. V. We leave this for further discussion in Sec.
V, where this Doppler modulation and the Doppler modulation due to the motion
of the GW detector around the solar system barycenter (SSB) are combined, as
in [33].

##  IV Comparison with waveforms from isolated NS

As a limiting case, when the spinning NS is isolated and spin $\boldsymbol{S}$
is along the $z$-axis of the coordinate system, i.e., $\Omega_{\rm
pre}=\theta_{S}=0$, the waveforms in Eqs. (40) will reduce to those given in
[29]:

$\displaystyle h_{+}^{(1)}$
$\displaystyle=\frac{Gb^{2}I_{3}\epsilon\gamma\sin(2\iota)\cos(t(\Omega_{\rm
p}+\Omega_{\rm r}))}{c^{4}D}\,,$ (41a) $\displaystyle h_{\times}^{(1)}$
$\displaystyle=\frac{2Gb^{2}I_{3}\epsilon\gamma\sin\iota\sin(t(\Omega_{\rm
p}+\Omega_{\rm r}))}{c^{4}D}\,,$ (41b) $\displaystyle h_{+}^{(2)}$
$\displaystyle=-\frac{32Gb^{2}I_{3}\epsilon\kappa\left(1+\cos^{2}\iota\right)\cos(2t\Omega_{\rm
r})}{c^{4}D}\,,$ (41c) $\displaystyle h_{\times}^{(2)}$
$\displaystyle=-\frac{64Gb^{2}I_{3}\epsilon\kappa\cos\iota\sin(2t\Omega_{\rm
r})}{c^{4}D}\,,$ (41d) $\displaystyle h_{+}^{(3)}$
$\displaystyle=\frac{2Gb^{2}I_{3}\epsilon\gamma^{2}\left(1+\cos^{2}\iota\right)\cos(2t(\Omega_{\rm
p}+\Omega_{\rm r}))}{c^{4}D}\,,$ (41e) $\displaystyle h_{\times}^{(3)}$
$\displaystyle=\frac{4Gb^{2}I_{3}\epsilon\gamma^{2}\cos\iota\sin(2t(\Omega_{\rm
p}+\Omega_{\rm r}))}{c^{4}D}\,.$ (41f)

Note that some of the signs in the above equations are reversed because the
convention for the inclination ($\iota$) we use here [27] is equivalent to
($\pi-\iota$) in [29].

For small equatorial ellipticity $\varepsilon\equiv|I_{1}-I_{2}|/I_{3}\simeq
16\epsilon\kappa$ and $\gamma=0$, $h_{+,\times}^{(1)}=h_{+,\times}^{(3)}=0$
and $h_{+,\times}^{(2)}$ will reduce to the triaxial aligned waveforms (e.g.,
Eq. (4.223) of [51]). If $\kappa=0$ and $\gamma\neq 0$,
$h_{+,\times}^{(2)}=0$, $h_{+,\times}^{(1)}$ and $h_{+,\times}^{(3)}$ will
reduce to the biaxial waveforms (cf. Eq. (1) of [52]).

The values of the parameters used in the following analysis are intended to
make the signal amplitude as large as possible under current observational and
theoretical constraints. The NSs, measured by pulsar timing [53] or GW
observation [54], can have masses up to about two solar masses. The widely
accepted range of the moment of inertia for NSs resides in $1-3\times
10^{38}~{}{\rm kg~{}m^{2}}$ (see [55] and references therein). The parameters
that characterize the properties of the spinning NS are in accordance with the
discussion in [28], in which $\epsilon\ll\kappa\ll\gamma_{\rm max}$ and
$\kappa\sim O(\gamma^{2})$. Recent observations from Advanced LIGO and
Advanced Virgo constrain two recycled pulsars (PSR J0437-4715 and PSR
J0711-6830) to have equatorial ellipticities
($\varepsilon\equiv|I_{1}-I_{2}|/I_{3}\simeq 16\epsilon\kappa$) of less than
$10^{-8}$ [56]. An ellipticity of $10^{-8}$ is also a typical value used in
searching for CWs from small-ellipticity sources [57]. Population syntheses of
Galactic disk double NS systems detectable by LISA [58, 37, 50] and TianQin
[38] suggest that the existence of double NSs with orbital periods as low as 6
minutes. We assume that the orbital period under consideration is 6 minutes,
which maximizes the strength of the orbital precession. Based on the current
observations of 22 double NSs [31, 32], we select $10~{}{\rm ms}$ as the
typical spin period of the rapidly spinning NS [38], and $1~{}{\rm kpc}$ as
the typical distance since about half of the known double NSs are located near
that value [31, 32].

As an example, we assume the component masses of a double NS system
$m_{1}=m_{2}=2.0~{}M_{\odot}$, the orbital period $P_{\rm b}=6~{}{\rm min}$,
and the opening angle of $\boldsymbol{S}$ precession cone is $\theta_{\rm
S}=5\pi/12$. The spinning NS’s characteristic parameters are $I_{3}=2.0\times
10^{38}~{}{\rm kg~{}m^{2}}$, $\epsilon=3.6\times 10^{-6}$, $\kappa=1.75\times
10^{-4}$, $\gamma=5.0\times 10^{-2}$ (equatorial ellipticity
$\varepsilon\equiv|I_{1}-I_{2}|/I_{3}=1.0\times 10^{-8}$). From Eq. (19) and
Eq. (15), one can obtain $\Omega_{\rm r}=628.32~{}{\rm Hz}$, $\Omega_{\rm
p}=2.26\times 10^{-3}~{}{\rm Hz}$, and $\Omega_{\rm pre}=7.50\times
10^{-7}~{}{\rm Hz}$. During a free precession period ($T=2\pi/\Omega_{\rm
p}=2785~{}{\rm s}$), Fig. 3 shows the approximate solution for
$\delta\Omega_{i}$ in Eqs. (37) and (38), which can reach up to $1.45\times
10^{-6}~{}{\rm Hz}$ for $\delta\Omega_{1,2}$ and $3.62\times 10^{-8}~{}{\rm
Hz}$ for $\delta\Omega_{3}$. Note that the jagged profiles in this figure are
due to a reduced sampling rate over a long spin precession period, as in the
following figures.

Figure 3: The analytic approximate solution for $\delta\Omega_{i}$ in Eqs.
(37) and (38) during a free precession period ($T=2\pi/\Omega_{\rm
p}=2785~{}{\rm s}$). Assume that the binary mass
$m_{1}=m_{2}=2.0~{}M_{\odot}$, the orbital period $P_{\rm b}=6~{}{\rm min}$,
and the opening angle of the spin precession cone is $\theta_{\rm S}=5\pi/12$.
The spinning NS’s characteristic parameters are $I_{3}=2.0\times
10^{38}~{}{\rm kg~{}m^{2}}$, $\epsilon=3.6\times 10^{-6}$, $\kappa=1.75\times
10^{-4}$, $\gamma=5.0\times 10^{-2}$ (equatorial ellipticity
$\varepsilon\equiv|I_{1}-I_{2}|/I_{3}=1.0\times 10^{-8}$).
$\delta\Omega_{1,2}$ ($\delta\Omega_{3}$) can reach up to $1.45\times
10^{-6}~{}{\rm Hz}$ ($3.62\times 10^{-8}~{}{\rm Hz}$) for $\Omega_{\rm
r}=628.32~{}{\rm Hz}$, $\Omega_{\rm p}=2.26\times 10^{-3}~{}{\rm Hz}$,
$\Omega_{\rm pre}=7.50\times 10^{-7}~{}{\rm Hz}$ obtained with the above
values of the parameters.

Fig. 4 shows the different waveform components with spin precession
incorporated during a spin precession period ($T_{\rm pre}=2\pi/\Omega_{\rm
pre}=8.382\times 10^{6}~{}{\rm s}$). The values of the parameters used here
are the same as in Fig. 3, and we set the spin period of the NS $P_{\rm
s}=10~{}{\rm ms}$, the inclination angle of $\boldsymbol{J}$ with respect to
the line of sight $\iota=\pi/4$, and the distance to the observer $D=1~{}{\rm
kpc}$. The modulated amplitude profiles for a binary depend on $\iota$,
$\theta_{S}$, and $\Omega_{\rm pre}$, which can be several times larger or
smaller than the isolated case for different waveform components. The
amplitude modulations shown here are only due to the NS’s spin precession
caused by spin-orbit coupling, the Doppler modulation will be included in Sec.
V (see the discussion at the end of Sec. III).

Figure 4: The waveform components from the triaxial non-aligned NS with spin
precession incorporated (blue) during a spin precession period and in
comparison with those from the isolated NS (red). The values of the parameters
used here are the same as in Fig. 3, and set the spin period of the NS $P_{\rm
s}=10~{}{\rm ms}$, the inclination angle of $\boldsymbol{J}$ with respect to
the line of sight $\iota=\pi/4$, and the distance to the observer $D=1~{}{\rm
kpc}$. The modulated amplitude profiles in a binary depend on $\iota$,
$\theta_{S}$, and $\Omega_{\rm pre}$, which can be several times larger or
smaller than the isolated case for different waveform components. The inset in
each subplot shows the waveform during the first spin period of the NS.

According to Eqs. (41), the GW angular frequency components of an isolated NS
are $\Omega_{+,\times}^{\rm iso}=\Omega_{\rm r}+\Omega_{\rm p}$, $2\Omega_{\rm
r}$, and $2(\Omega_{\rm r}+\Omega_{\rm p})$ for both $+$ and $\times$
polarizations. Spin-orbit coupling in the binary can split these frequencies
in the following way

$\displaystyle\Omega_{+}^{\rm iso}$
$\displaystyle\longrightarrow\Omega_{+}^{\rm iso}\pm n\Omega_{\rm
pre}~{}(n=0,1,2)\,,$ (42a) $\displaystyle\Omega_{\times}^{\rm iso}$
$\displaystyle\longrightarrow\Omega_{\times}^{\rm iso}\pm n\Omega_{\rm
pre}~{}(n=1,2)\,.$ (42b)

Similar to the analysis of the GW spectrum of isolated systems [28], the
spectral analysis of the above frequency components can be used to infer the
orbital period of the binary and the characteristic parameters of the NS.

In order to quantitatively measure the degree of matching between these two
types of waveforms, one can calculate the fitting factor (FF) [59] between the
genuine GW waveforms generated by the spinning NS in a binary (denoted as
$h_{\rm b}$) and the ones by an isolated NS (denoted as $h_{\rm i}$), given
the latter has been used in the matched filtering of CW data analysis

${\rm{FF}}\equiv\max_{\boldsymbol{\lambda}}\frac{\left(h_{\rm b},h_{\rm
i}\right)}{\sqrt{\left(h_{\rm b},h_{\rm b}\right)(h_{\rm i},h_{\rm i})}}\,,$
(43)

where ${\boldsymbol{\lambda}}$ is the set of the parameters that characterize
the waveforms. For a quasi-monochromatic signal, the inner product
$\left(h_{\rm b},h_{\rm i}\right)$ can be simplified as $\left(h_{\rm
b},h_{\rm i}\right)\equiv\int_{0}^{T_{\rm obs}}h_{\rm b}(t)h_{\rm i}(t)dt\,,$
where $T_{\rm obs}$ is the observation time. Fig. 5 shows the FFs of different
waveform components in $h_{\rm b}$ and $h_{\rm i}$ as a function of $T_{\rm
obs}$. As we can see, the two waveforms only match well (${\rm FF}>0.97$)
within roughly a few days ($\sim O(10^{5}~{}{\rm s})$) and then start to
diverge rapidly.

Figure 5: The fitting factors between the isolated case and the one with spin
precession for different waveform components shown in Fig. 4.

##  V Detecting GW from spinning NS in a binary

###  V.1 Signal model in detector coordinate system

The GW strain signal from a spinning triaxial non-aligned NS in a binary can
be written as a sum of different waveform components $h_{n}(t)~{}(n=1,2,3)$ as
follows [17]

$\displaystyle
h(t)=\sum_{n=1}^{3}h_{n}(t)=\sum_{n=1}^{3}(F_{+}(t)H_{+}^{(n)}(t)+F_{\times}(t)H_{\times}^{(n)}(t))\,,$
(44)

where $H_{+,\times}^{(n)}(t)$ are the Doppler-modulated waveforms (cf. Eqs.
(51) below), and $F_{+,\times}(t)$ are the antenna pattern functions of GW
detector, of which the explicit expressions can be found in [17].
$F_{+,\times}(t)$ depend on a set of parameters listed below: $\gamma_{\rm o}$
characterizes the orientation of the detector with respect to the local
geographical directions, $\zeta$ denotes the angle between the interferometer
arms, $\lambda$ is the geographical latitude of the detector’s site,
$(\alpha,\delta)$ are the right ascension and declination of the source,
$\psi_{\rm p}$ is the GW polarization angle, $\Omega_{\rm Er}$ is the
rotational angular frequency of the Earth, and $\phi_{\rm r}$ is the initial
phase of the Earth’s diurnal motion.

Below we will consider how to incorporate Doppler modulation into the phases
of the waveforms from the spinning NS undergoing spin precession (cf. Eqs.
(40)). Fig. 6 shows the binary coordinate system $(X_{b},Y_{b},Z_{b})$ with
$Z_{b}$ axis aligned with the binary’s total angular momentum $\boldsymbol{J}$
and the motion of spinning NS within it. $(X_{L},Y_{L},Z_{L})$ is the
coordinate system with $Z_{L}$ axis aligned with the binary’s orbital angular
momentum $\boldsymbol{L}$ (hereafter referred to as $L$-aligned coordinate
system), in which $X_{L}-Y_{L}$ represents the orbital plane. The origins of
these two frames are both placed at the BB ($O$). The opening angle of
$\boldsymbol{L}$ precession cone is $\theta_{L}$. Suppose that at initial
time, $\boldsymbol{S}$ and $\boldsymbol{L}$ are both in the $Y_{b}-Z_{b}$
plane and the spinning NS sits on the $X_{b}$ axis. After a period of time
$t$, the orbital plane precessed by $\phi_{\rm N}$, the spinning NS’s position
vector and the orbital longitude are $\boldsymbol{r}_{1}$ and $\psi_{1}$, and
its spin is $\boldsymbol{S}(t)$. (The $Y$-axes not drawn in Fig. 6 are
determined by the right-hand rule.)

In the $L$-aligned coordinate system, the position vector of the spinning NS
$\boldsymbol{r}_{1L}=(r_{1}\cos(\omega_{\rm b}t),r_{1}\sin(\omega_{\rm
b}t),0)$. By two rotation transformations (first rotates by $\theta_{L}$
clockwise about the $X_{L}$ axis, and then rotates by $\phi_{N}$ clockwise
about the $Z_{b}$ axis), the position vector $\boldsymbol{r}_{1}$ in the
binary coordinate system can be given by

$\displaystyle\boldsymbol{r}_{1}=r_{1}\left(\begin{array}[]{c}\cos\psi_{1}\cos\phi_{N}-\cos\theta_{L}\sin\psi_{1}\sin\phi_{N}\\\
\cos\theta_{L}\sin\psi_{1}\cos\phi_{N}+\cos\psi_{1}\sin\phi_{N}\\\
\sin\theta_{L}\sin\psi_{1}\\\ \end{array}\right)$ (48)

with $\psi_{1}=\omega_{\rm b}t$ and $\phi_{N}=\Omega_{\rm pre}t$.

Figure 6: The binary coordinate system $(X_{b},Y_{b},Z_{b})$ and the motion of
the spinning NS within it. $(X_{L},Y_{L},Z_{L})$ is the coordinate system with
$Z_{L}$ axis aligned with the binary’s orbital angular momentum
$\boldsymbol{L}$, $X_{L}-Y_{L}$ represents the orbital plane. The opening
angle of $\boldsymbol{L}$ precession cone is $\theta_{L}$. The distant
detector with position vector $\boldsymbol{D}$ is assumed in the $Y_{b}-Z_{b}$
plane at colatitude $\iota$ from $\boldsymbol{J}$.

In the detector coordinate system, the Doppler shift of the GW frequency
$f_{0}$ from the spinning NS in a binary system can be expressed as the
combination of detector Doppler shift around the SSB and source Doppler shift
around the BB

$\Delta f_{\rm{D}}\approx f_{0}\left(\frac{{\boldsymbol{n}}\cdot{\rm
d}\boldsymbol{r}_{\rm d}/{\rm d}t}{c}+\frac{{\boldsymbol{n}}_{\rm b}\cdot{\rm
d}\boldsymbol{r}_{1}/{\rm d}t}{c}\right)\,,$ (49)

with [17]

$\displaystyle\boldsymbol{n}\cdot\boldsymbol{r}_{\rm d}$
$\displaystyle=R_{\rm{E}}[\cos\lambda\cos\delta\cos(\alpha-\phi_{r}-\Omega_{\rm
Er}t)$ (50)
$\displaystyle+\sin\lambda\sin\delta]+R_{\rm{ES}}[\cos\alpha\cos\delta\cos(\phi_{\rm
o}+\Omega_{\rm Eo}t)$ $\displaystyle+(\cos\varepsilon_{\rm
e}\sin\alpha\cos\delta+\sin\varepsilon_{\rm e}\sin\delta)\sin(\phi_{\rm
o}+\Omega_{\rm Eo}t)]$

being the projection of the detector’s position vector $\boldsymbol{r}_{\rm
d}$ along the spinning NS’s line of sight $\boldsymbol{n}$ in SSB coordinate
system, where $R_{\rm E}$ and $R_{\rm ES}$ are the mean radius of the Earth
and the mean distance from the Earth’s center to the SSB, $\Omega_{\rm Eo}$ is
the mean orbital angular frequency of the Earth, $\phi_{\rm o}$ is the initial
phase of the Earth’s annual motion, and $\varepsilon_{\rm e}$ is the ecliptic
obliquity. $-{\boldsymbol{n}}_{\rm b}=(0,\sin{\iota},\cos{\iota})$ is the
SSB’s location in the binary coordinate system, orbital radius
$r_{1}=G^{1/3}m_{2}/(\omega_{\rm b}M)^{2/3}$ and angular frequency
$\omega_{\rm b}=2\pi/P_{\rm b}$ for a circular orbit. If there is no orbital
precession, i.e., $\Omega_{\rm pre}=0$, then the Doppler shift in Eq. (49)
reduces to the simple case as Eq. (6) in [33].

Since the waveform components in Eqs. (40) can be decomposed into a series of
sine and cosine functions in which the frequencies are linear combinations of
$\Omega_{\rm r}$, $\Omega_{\rm p}$, and $\Omega_{\rm pre}$, the Doppler-
modulated waveforms $H_{+,\times}^{(n)}(t)$ in Eq. (44) can be obtained by

$\displaystyle H_{+}^{(n)}(t)$
$\displaystyle=h_{+}^{(n)}(\Omega_{l}\rightarrow\Omega_{l}+2\pi\Delta
f_{\rm{D}})\,,$ (51a) $\displaystyle H_{\times}^{(n)}(t)$
$\displaystyle=h_{\times}^{(n)}(\Omega_{l}\rightarrow\Omega_{l}+2\pi\Delta
f_{\rm{D}})\,,$ (51b)

with $l\in\\{\rm r,p,pre\\}$.

The binary is assumed to follow an invariant circular orbit when calculating
the Doppler shift of the spinning NS. In fact, the orbit is constantly
shrinking due to gravitational radiation. According to the orbital velocity
evolution equation ($d(v/c)/dt\propto(v/c)^{9}$, cf. [60]), its relative
variation is $\lesssim 10^{-5}$ for a half-year observation, resulting in a
relative variation for the Doppler shift of $\lesssim 10^{-5}$. Thus, we can
ignore the reaction of gravitational radiation on the orbit for the
observation time under consideration.

###  V.2 Effects of spin-orbit coupling on parameter estimation

In addition to spin precession, spin-orbit coupling also causes orbital plane
precession, which is expected to introduce additional information into the
Doppler-modulated waveforms $H_{+,\times}^{(n)}(t)$ (cf. Eqs. (51)). We use
the Fisher information matrix (FIM) to obtain a quantitative assessment of the
parameter estimation accuracy for GW detection (e.g., see [61]). For GW signal
$h(t)$ (cf. Eq. (44)) with parameter set $\boldsymbol{\lambda}$, FIM is
defined as

${{\Gamma}^{ij}}\equiv\left(\frac{\partial
h}{\partial{{\lambda}_{i}}},\frac{\partial
h}{\partial{{\lambda}_{j}}}\right)\,.$ (52)

For a monochromatic signal of frequency $f$, the noise-weighted inner product
$\left(a,b\right)\simeq\frac{2}{S_{n}\left(f\right)}\int_{0}^{T_{\rm
obs}}a(t)b(t)dt$ [62], where $S_{n}(f)$ is the power spectral density of the
instrumental noise at frequency $f$, and $T_{\rm obs}$ is the observation
time. The optimal signal-to-noise ratio (SNR) for signal detection is defined
as ${\rm SNR}\equiv(h,h)^{1/2}$. The root-mean-square (RMS) error of parameter
${\lambda_{i}}$ is estimated as $\Delta{\lambda_{i}}=\sqrt{{\Sigma_{ii}}}$,
where the covariance matrix $\Sigma$ is the inverse of the FIM, i.e.,
$\Sigma={\Gamma^{-1}}$. We extend the parameter set in [61], as our work
considers the spin-orbit coupling for a triaxial non-aligned NS, which results
in a parameter set
$\boldsymbol{\lambda}=(\ln{h_{10}},\ln{h_{20}},\ln{h_{30}},\alpha,\sin{\delta},\ln{P_{\rm
b}},\cos{\iota},\theta_{S},\psi_{\rm p},\ln{\Omega_{\rm r}},\ln{\Omega_{\rm
p}})$, in which the amplitudes for the different waveform components are
defined as

$h_{10}=\frac{{2G{b^{2}}{I_{3}}\epsilon\gamma}}{{{c^{4}}D}},h_{20}=\frac{{64G{b^{2}}{I_{3}}\epsilon\kappa}}{{{c^{4}}D}},h_{30}=\frac{4Gb^{2}{I_{3}}\epsilon\gamma^{2}}{c^{4}D}\,.$
(53)

Logarithms are taken for some parameters because the relative errors for them
are more meaningful than the absolute errors. For example,
$\Delta\ln{h_{10}}=\Delta h_{10}/h_{10}$ is simply the relative error in
$h_{10}$. Since $\boldsymbol{J}=\boldsymbol{L}+\boldsymbol{S}$, then
$\theta_{L}$ is related to $\theta_{S}$ by
$S/\sin(\theta_{L})=L/\sin(\theta_{S})$. According to Eq. (19) and Kepler’s
third law, $\Omega_{\rm pre}\propto P_{\rm b}^{-5/3}$ and $r_{1}\propto P_{\rm
b}^{2/3}$. Therefore, the parameters $(\theta_{L},\Omega_{\rm pre},r_{1})$ are
related to $\theta_{S}$ and $P_{\rm b}$ used in the FIM, so they are not
included in the parameter set $\boldsymbol{\lambda}$. Similar to the sky
localization error defined in [63], we define the corresponding one for a
source located at $(\alpha,\delta)$,
$\Delta\Omega=2\pi\sqrt{\Sigma_{\alpha\alpha}\Sigma_{\sin\delta\sin\delta}-\Sigma_{\alpha\sin\delta}^{2}}$.

In the following analysis, since the sensitivity of Einstein Telescope [41] in
the frequency band of interest is not as good as that of Cosmic Explorer, we
use Cosmic Explorer, which consists of two facilities (one 40 km on a side and
one 20 km on a side), each with a single L-shaped detector [40]. Various
angular parameters are taken as $\zeta=\pi/2,\lambda=0.764,\gamma_{\rm
o}=1.5,\phi_{\rm r}=\phi_{\rm o}=0,\alpha=1.209,\delta=1.475,\psi_{\rm
p}=1.0$. Assuming a 40 km arm length with a low-frequency optimized
sensitivity is used, the amplitude spectral density of the instrumental noise
$\sqrt{S_{n}(f_{0})}=1.36\times 10^{-25}~{}{\rm Hz}^{-1/2}$ and
$\sqrt{S_{n}(2f_{0})}=1.66\times 10^{-25}~{}{\rm Hz}^{-1/2}$ for
$f_{0}=100~{}{\rm Hz}$, the observation time $T_{\rm obs}=2T_{\rm
pre}=1.676\times 10^{7}~{}{\rm s}$.

After calculating the FIM numerically with Mathematica for Eq. (52), the RMS
errors of the estimated parameters for a typical spinning NS in a binary (the
parameters used here are the same as in Fig. 3 and Fig. 4) are shown in Fig. 7
as red downward-pointing triangles. For comparison, the results without spin-
orbit coupling are represented by blue open squares and the results for the
triaxial aligned case are shown as black circles. For the first three Doppler-
modulated signal components $h_{n}(t)~{}(n=1,2,3)$ in Eq. (44), their SNRs are
104, 8.5, and 7.6, respectively. In comparison, the corresponding SNRs are
74.3, 12.6, and 11.3 for the signals without spin-orbit coupling. The total
SNR for the triaxial aligned case is also 12.6 as that for $h_{2}$ without
spin-orbit coupling, since the SNRs for $h_{1}$ and $h_{3}$ become zero when
$\gamma=0$ (cf. Sec. IV). The spin and orbital precession modulated signal
$h_{1}(t)$ increases its SNR by 40% compared with the case without spin-orbit
coupling, while the other two signals $h_{2}(t)$ and $h_{3}(t)$ both decrease
by 33%. We can get some clues to understand the changes in SNRs from the
limiting case of the waveforms in Sec. IV. From Eqs. (41) ($\iota$ corresponds
to $\iota-\theta_{S}$ under spin-orbit coupling), we can see that
$h_{+,\times}^{(1)}$ is proportional to $\sin(2\iota)$ or $\sin{\iota}$, while
$h_{+,\times}^{(2)}$ and $h_{+,\times}^{(3)}$ are proportional to
$(1+\cos^{2}{\iota})$ or $\cos{\iota}$. The spin-orbit coupling makes $\iota$
larger during one precession period (cf. $\iota-\theta_{S}$ in Fig. 2),
therefore the SNR of $h_{1}(t)$ becomes larger and that of $h_{2}(t)$ and
$h_{3}(t)$ becomes smaller.

The fractional estimation errors for three amplitudes are inversely
proportional to their SNRs. So the parameter estimations for the triaxial non-
aligned cases are more accurate than those for the triaxial aligned case. For
the triaxial non-aligned case, the precession improves the sky localization by
a factor of two and slightly improves the estimation of the orbital period.
The improvement in sky localization is small because for double NS systems,
whose orbital angular momentum is very close to the total angular momentum,
the modulation of the orbital angular momentum is relatively small. The most
significant improvement comes from the estimation of the angles $\cos{\iota}$
and $\theta_{S}$, both of which are improved by about 3 orders of magnitude.
This is because these two angles are encoded in the amplitudes of the
waveforms (cf. Eqs. (40)), and they modulate the profiles of the waveforms in
Fig. 4.

Figure 7: The RMS errors (red downward-pointing triangles) of the signal
parameters for a typical spinning NS in a binary using Cosmic Explorer
detector. The effects of spin-orbit coupling are considered in the waveforms.
For comparison, the blue open squares represent the result without spin-orbit
coupling and the black circles represent the triaxial aligned case (the
parameters $(h_{10},h_{30},\Omega_{\rm p})$ do not exist in this simplest
model, cf. Sec. IV). The error of the sky location $(\alpha,\delta)$ is
combined into the elliptical area $\Delta\Omega$. The parameter estimations
for the triaxial non-aligned cases are more accurate than those for the
triaxial aligned case. After taking into account spin-orbit coupling, the
estimation accuracy of the parameters $\cos{\iota}$ and $\theta_{S}$ (in green
shadow) are improved by about 3 orders of magnitude.

##  VI Conclusions

In this work, we calculate the gravitational waveforms of a triaxial non-
aligned NS in a compact binary system in which the effects of spin-orbit
coupling have been incorporated. Then, we compare our waveforms with the ones
commonly used in current CWs searches. Finally, we evaluated the parameter
estimation accuracy for the signal detected by the proposed next-generation GW
detector using the Fisher information matrix method.

For a tight double NS system with a 6-min orbital period, by solving the
precession equation with the perturbation method, we find that spin
precession-induced correction to the spin angular frequencies of NS is in the
same order of magnitude as the angular frequency of orbital precession. The
fitting factor between the waveforms with and without spin precession will
drop to less than 0.97 after a few days. The analytic waveforms show that
spin-orbit coupling introduces additional modulation information that will
help in improving the accuracy of parameter estimation in CW detection.

The double NS system (consisting of a rapidly spinning NS and a nonspinning
NS) considered in this work can be seen as a dual-line GW source, in which the
orbital motion of the binary will emit low-frequency GWs in the mHz band in
addition to the high-frequency GWs from the spinning NS. This dual-line GW
source is of astrophysical interest, such as constraining the NS’s moment of
inertia and ellipticity using the ratio of the strain amplitudes of the low-
and high-frequency GWs [64] or combining the angular momentum loss of the NS
[65]. Since the angular frequency of the orbital precession contains
information about the orbital period and mass of the binary, we can use it to
infer binary parameters by combining information from the emitted GWs of the
dual-line sources, and such studies are currently under our investigation.

###### Acknowledgements.

Y.W. gratefully acknowledges support from the National Key Research and
Development Program of China (No. 2022YFC2205201 and No. 2020YFC2201400), the
National Natural Science Foundation of China (NSFC) under Grants No. 11973024,
Major Science and Technology Program of Xinjiang Uygur Autonomous Region (No.
2022A03013-4), and Guangdong Major Project of Basic and Applied Basic Research
(Grant No. 2019B030302001). T. L. is supported by NSFC Grant No. 12003008 and
the China Postdoctoral Science Foundation Grant No. 2020M682393. J.-W. C.
acknowledges the support from China Postdoctoral Science Foundation under
Grant No. 2021M691146. S.D.M is supported by U.S. National Science Foundation
(NSF) grant PHY-2207935. We thank the anonymous referee for helpful comments
and suggestions.

##  Appendix A Expressions for $\mathcal{R}$ and $\mathcal{A}$

The symmetric matrix $\mathcal{A}$ used in calculating the waveforms (cf. Eqs.
(40)) can be explicitly expressed as follows

$\displaystyle\mathcal{A}_{11}$
$\displaystyle=2bI_{3}\epsilon[16b\kappa+\gamma\sin(t\Omega_{\rm
p})[b\gamma\sin(t\Omega_{\rm p})+2\Omega_{\rm
pre}\sin\theta_{S}\sin(t\Omega_{\rm r})]]\,,$ (54a)
$\displaystyle\mathcal{A}_{22}$
$\displaystyle=2bI_{3}\epsilon[-16b\kappa+\gamma\cos(t\Omega_{\rm
p})[b\gamma\cos(t\Omega_{\rm p})+4\Omega_{\rm
pre}\sin\theta_{S}\sin^{2}(t\Omega_{\rm r}/2)]]\,,$ (54b)
$\displaystyle\mathcal{A}_{33}$ $\displaystyle=-2b\gamma
I_{3}\epsilon[b\gamma+4\Omega_{\rm pre}\sin\theta_{S}\sin(t\Omega_{\rm
r}/2)\sin(t(2\Omega_{\rm p}+\Omega_{\rm r})/2)]\,,$ (54c)
$\displaystyle\mathcal{A}_{12}$ $\displaystyle=-b\gamma
I_{3}\epsilon[b\gamma\sin(2t\Omega_{\rm p})+2\Omega_{\rm
pre}\sin\theta_{S}[\sin(t\Omega_{\rm p})-\sin(t(\Omega_{\rm p}-\Omega_{\rm
r}))]]\,,$ (54d) $\displaystyle\mathcal{A}_{23}$
$\displaystyle=b^{2}I_{3}\epsilon\gamma\sin(t\Omega_{\rm p})\,,$ (54e)
$\displaystyle\mathcal{A}_{31}$
$\displaystyle=bI_{3}\epsilon[b\gamma\cos(t\Omega_{\rm p})+\Omega_{\rm
pre}\sin\theta_{S}]\,.$ (54f)

The transformation matrix $\mathcal{R}$ used in calculating the waveforms (cf.
Eqs. (40)) can be explicitly expressed as follows

$\displaystyle{\mathcal{R}_{x1}}$
$\displaystyle=\sin(t{\Omega_{\rm{p}}})\Big{[}\cos(t{\Omega_{{\rm{pre}}}})\cos(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))-\cos{\theta_{S}}\sin(t{\Omega_{{\rm{pre}}}})\sin(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))\Big{]}+\frac{1}{2}\cos(t{\Omega_{\rm{p}}})$
(55a)
$\displaystyle\times\Big{[}\cos(t{\Omega_{{\rm{pre}}}})[16\kappa\sin(t{\Omega_{\rm{p}}})\big{(}\cos(t{\Omega_{\rm{r}}})+\kappa[3\cos(t(2{\Omega_{\rm{p}}}-{\Omega_{\rm{r}}}))-8\cos(t{\Omega_{\rm{r}}})+\cos(t(2{\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))]\big{)}$
$\displaystyle+({\gamma^{2}}-2)\sin(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))]+\sin(t{\Omega_{{\rm{pre}}}})\big{[}-2\gamma\sin{\theta_{S}}+\cos{\theta_{S}}\Big{(}({\gamma^{2}}-2)\cos(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))$
$\displaystyle+16\kappa\sin(t{\Omega_{\rm{p}}})\big{(}3\kappa\sin(t(2{\Omega_{\rm{p}}}-{\Omega_{\rm{r}}}))+(8\kappa-1)\sin(t{\Omega_{r}})-\kappa\sin(t(2{\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))\big{)}\Big{)}\big{]}\Big{]}\,,$
$\displaystyle{\mathcal{R}_{x2}}$
$\displaystyle=\frac{1}{2}\sin(t{\Omega_{\rm{p}}})\Big{[}\sin(t{\Omega_{{\rm{pre}}}})(({\gamma^{2}}-2)\cos{\theta_{S}}\cos(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))-2\gamma(1+8\kappa)\sin{\theta_{S}})$
(55b)
$\displaystyle+({\gamma^{2}}-2)\cos(t{\Omega_{{\rm{pre}}}})\sin(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))\Big{]}-\cos(t{\Omega_{\rm{p}}})\Big{[}\cos{\theta_{S}}\sin(t{\Omega_{{\rm{pre}}}})(8\kappa[\cos(t{\Omega_{\rm{r}}})+\kappa(3\cos(t(2{\Omega_{\rm{p}}}-{\Omega_{\rm{r}}}))$
$\displaystyle-8\cos(t{\Omega_{\rm{r}}})+\cos(t(2{\Omega_{\rm{p}}}+{\Omega_{\rm{r}}})))]\sin(t{\Omega_{\rm{p}}})-\sin(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}})))+\cos(t{\Omega_{{\rm{pre}}}})\Big{(}\cos(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))$
$\displaystyle+8\kappa\sin(t{\Omega_{\rm{p}}})[\sin(t{\Omega_{\rm{r}}})+\kappa(-3\sin(t(2{\Omega_{\rm{p}}}-{\Omega_{\rm{r}}}))-8\sin(t{\Omega_{\rm{r}}})+\sin(t(2{\Omega_{\rm{p}}}+{\Omega_{\rm{r}}})))]\Big{)}\Big{]}\,,$
$\displaystyle{\mathcal{R}_{x3}}$
$\displaystyle=\gamma\cos{\theta_{S}}(1+4\kappa-4\kappa\cos(2t{\Omega_{\rm{p}}}))\cos(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))\sin(t{\Omega_{{\rm{pre}}}})+\frac{1}{2}({\gamma^{2}}-2)\sin{\theta_{S}}\sin(t{\Omega_{{\rm{pre}}}})$
(55c)
$\displaystyle+\gamma\cos(t{\Omega_{{\rm{pre}}}})(1+8\kappa{\sin^{2}}(t{\Omega_{\rm{p}}}))\sin(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))\,,$
$\displaystyle{\mathcal{R}_{y1}}$
$\displaystyle=\frac{1}{2}\Big{[}2(\cos(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))\sin(t{\Omega_{\rm{p}}})+4\kappa[\cos(t{\Omega_{\rm{r}}})+\kappa(3\cos(t(2{\Omega_{\rm{p}}}-{\Omega_{\rm{r}}}))-8\cos(t{\Omega_{\rm{r}}})$
(55d)
$\displaystyle+\cos(t(2{\Omega_{\rm{p}}}+{\Omega_{\rm{r}}})))]\sin(2t{\Omega_{\rm{p}}}))\sin(t{\Omega_{{\rm{pre}}}})+\cos(t{\Omega_{\rm{p}}})(2\gamma\cos(t{\Omega_{{\rm{pre}}}})\sin{\theta_{S}}+({\gamma^{2}}-2)\sin(t{\Omega_{{\rm{pre}}}})\sin(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}})))$
$\displaystyle+\cos{\theta_{S}}\cos(t{\Omega_{\rm{pre}}})\Big{(}-({\gamma^{2}}-2)\cos(t{\Omega_{\rm{p}}})\cos(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))+2\sin(t{\Omega_{\rm{p}}})\sin(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))+8\kappa\sin(2t{\Omega_{\rm{p}}})[\sin(t{\Omega_{\rm{r}}})$
$\displaystyle+\kappa(-3\sin(t(2{\Omega_{\rm{p}}}-{\Omega_{\rm{r}}}))-8\sin(t{\Omega_{\rm{r}}})+\sin(t(2{\Omega_{\rm{p}}}+{\Omega_{\rm{r}}})))]\Big{)}\Big{]}\,,$
$\displaystyle{\mathcal{R}_{y2}}$
$\displaystyle=\frac{1}{2}\Big{[}2\gamma(1+8\kappa)\cos(t{\Omega_{{\rm{pre}}}})\sin{\theta_{S}}\sin(t{\Omega_{\rm{p}}})+\cos{\theta_{S}}\cos(t{\Omega_{{\rm{pre}}}})(-({\gamma^{2}}-2)\sin(t{\Omega_{\rm{p}}})\cos(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))$
(55e)
$\displaystyle+8\kappa[\cos(t{\Omega_{\rm{r}}})+\kappa(3\cos(t(2{\Omega_{\rm{p}}}-{\Omega_{\rm{r}}}))-8\cos(t{\Omega_{\rm{r}}})+\cos(t(2{\Omega_{\rm{p}}}+{\Omega_{\rm{r}}})))]\sin(2t{\Omega_{\rm{p}}})-2\cos(t{\Omega_{\rm{p}}})\sin(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}})))$
$\displaystyle+\sin(t{\Omega_{{\rm{pre}}}})\Big{(}-2\cos(t{\Omega_{\rm{p}}})\cos(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))+({\gamma^{2}}-2)\sin(t{\Omega_{\rm{p}}})\sin(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))+8\kappa\sin(2t{\Omega_{\rm{p}}})[3\kappa\sin(t(2{\Omega_{\rm{p}}}-{\Omega_{\rm{r}}}))$
$\displaystyle+(8\kappa-1)\sin(t{\Omega_{\rm{r}}})-\kappa\sin(t(2{\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))]\Big{)}\Big{]}\,,$
$\displaystyle{\mathcal{R}_{y3}}$
$\displaystyle=\gamma\cos{\theta_{S}}(-1-4\kappa+4\kappa\cos(2t{\Omega_{\rm{p}}}))\cos(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))\cos(t{\Omega_{{\rm{pre}}}})-\frac{1}{2}({\gamma^{2}}-2)\sin{\theta_{S}}\cos(t{\Omega_{{\rm{pre}}}})$
(55f)
$\displaystyle+\gamma\sin(t{\Omega_{{\rm{pre}}}})(1+8\kappa{\sin^{2}}(t{\Omega_{\rm{p}}}))\sin(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))\,,$
$\displaystyle{\mathcal{R}_{z1}}$
$\displaystyle=\gamma\cos{\theta_{S}}\cos(t{\Omega_{\rm{p}}})+\frac{1}{2}\sin{\theta_{S}}[({\gamma^{2}}-2)\cos(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))\cos(t{\Omega_{\rm{p}}})-2\sin(t{\Omega_{\rm{p}}})\sin(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))$
(55g)
$\displaystyle+8\kappa\sin(2t{\Omega_{\rm{p}}})(3\kappa\sin(t(2{\Omega_{\rm{p}}}-{\Omega_{\rm{r}}}))+(8\kappa-1)\sin(t{\Omega_{\rm{r}}})-\kappa\sin(t(2{\Omega_{\rm{p}}}+{\Omega_{\rm{r}}})))]\,,$
$\displaystyle{\mathcal{R}_{z2}}$
$\displaystyle=\gamma(1+8\kappa)\cos{\theta_{S}}\sin(t{\Omega_{\rm{p}}})+\frac{1}{2}\sin{\theta_{S}}\Big{[}({\gamma^{2}}-2)\cos(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))\sin(t{\Omega_{\rm{p}}})+2\cos(t{\Omega_{\rm{p}}})\sin(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))$
(55h)
$\displaystyle-8\kappa[\cos(t{\Omega_{\rm{r}}})+\kappa(3\cos(t(2{\Omega_{\rm{p}}}-{\Omega_{\rm{r}}}))-8\cos(t{\Omega_{\rm{r}}})+\cos(t(2{\Omega_{\rm{p}}}+{\Omega_{\rm{r}}})))]\sin(2t{\Omega_{\rm{p}}})\Big{]}\,,$
$\displaystyle{\mathcal{R}_{z3}}$
$\displaystyle=\frac{1}{2}(2-{\gamma^{2}})\cos{\theta_{S}}+\gamma\cos(t({\Omega_{\rm{p}}}+{\Omega_{\rm{r}}}))\sin{\theta_{S}}(1+8\kappa\sin^{2}(t{\Omega_{\rm{p}}}))\,.$
(55i)

##  Appendix B Expressions for $O(\Omega_{\rm pre})$ waveform components

The waveform components to the order $O(\Omega_{\rm pre})$ can be expressed as
follows

$\displaystyle h_{+}^{({\rm pre})}$
$\displaystyle=\frac{{Gb{I_{3}}\epsilon{\Omega_{{\rm{pre}}}}\sin{\theta_{S}}}}{{4{c^{4}}D}}[\cos(t\Omega_{\rm{r}})(\sin(2{\theta_{S}})(-(3+\cos(2\iota))\cos(2t{\Omega_{{\rm{pre}}}})+6{\sin^{2}}\iota)$
(56a)
$\displaystyle+4\cos(2{\theta_{S}})\cos(t{\Omega_{\rm{pre}}})\sin(2\iota))+2(-2\cos{\theta_{S}}\sin(2\iota)\sin(t{\Omega_{{\rm{pre}}}})+(3+\cos(2\iota))\sin{\theta_{S}}\sin(2t{\Omega_{{\rm{pre}}}}))\sin(t{\Omega_{{\rm{r}}}})]\,,$
$\displaystyle h_{\times}^{({\rm{pre}})}$
$\displaystyle=\frac{{Gb{I_{3}}\epsilon{\Omega_{{\rm{pre}}}}\sin{\theta_{S}}}}{{{c^{4}}D}}[\cos(t\Omega_{\rm{r}})(2\cos(2{\theta_{S}})\sin\iota\sin(t{\Omega_{{\rm{pre}}}})-\cos\iota\sin(2{\theta_{S}})\sin(2t{\Omega_{{\rm{pre}}}}))$
(56b)
$\displaystyle+2(-\cos\iota\sin{\theta_{S}}\cos(2t{\Omega_{{\rm{pre}}}})+\cos{\theta_{S}}\cos(t{\Omega_{{\rm{pre}}}})\sin\iota)\sin(t\Omega_{\rm{r}})]\,.$

We can see that there are components with frequencies of $\Omega_{\rm r}$,
$\Omega_{\rm r}\pm\Omega_{\rm pre}$, and $\Omega_{\rm r}\pm 2\Omega_{\rm pre}$
in these waveforms. Since ${h_{+,\times}^{({\rm
pre})}}/{h_{+,\times}^{(1)}}\sim{\Omega_{\rm
pre}\sin\theta_{S}}/{(b\gamma)},~{}{h_{+,\times}^{({\rm
pre})}}/{h_{+,\times}^{(2)}}\sim{\Omega_{\rm
pre}\sin\theta_{S}}/{(16b\kappa)}$, and ${h_{+,\times}^{({\rm
pre})}}/{h_{+,\times}^{(3)}}\sim{\Omega_{\rm
pre}\sin\theta_{S}}/{(b\gamma^{2})}$, for typical parameters used as in Fig.
3, these $h_{+,\times}^{({\rm pre})}$ components are negligible compared to
$h_{+,\times}^{(1)}$, $h_{+,\times}^{(2)}$, and $h_{+,\times}^{(3)}$.

##  Appendix C The residual of two solutions

To confirm the fidelity of the approximate calculation, we use Euler angles to
accurately calculate the angular frequency. Although this method can give the
exact solution $\omega_{i}^{E}$, it is too complicated to give a simple
analytical waveform like Eqs. (40). They can be calculated as (an overdot
represents $d/dt$) [45]

$\displaystyle\omega_{1}^{E}$ $\displaystyle=\dot{\phi}_{\rm b}\sin\theta_{\rm
b}\sin\psi_{\rm b}+\dot{\theta}_{\rm b}\cos\psi_{\rm b}\,,$ (57a)
$\displaystyle\omega_{2}^{E}$ $\displaystyle=\dot{\phi}_{\rm b}\sin\theta_{\rm
b}\cos\psi_{\rm b}-\dot{\theta}_{\rm b}\sin\psi_{\rm b}\,,$ (57b)
$\displaystyle\omega_{3}^{E}$ $\displaystyle=\dot{\phi}_{\rm b}\cos\theta_{\rm
b}+\dot{\psi}_{\rm b}\,,$ (57c)

with the Euler angles derived from the primitive (no approximation) rotation
matrix (cf. Eq. (21)),

$\displaystyle\theta_{\rm b}$
$\displaystyle=\arccos\left(\mathcal{R}_{z3}\right)\,,$ (58a)
$\displaystyle\phi_{\rm b}$
$\displaystyle=\arctan\left(-\frac{\mathcal{R}_{x3}}{\mathcal{R}_{y3}}\right)\,,$
(58b) $\displaystyle\psi_{\rm b}$
$\displaystyle=\arctan\left(\frac{\mathcal{R}_{z1}}{\mathcal{R}_{z2}}\right)\,.$
(58c)

The absolute errors between the approximate angular frequency $\omega_{i}$ and
the exact angular frequency $\omega_{i}^{E}$ are shown in Fig. 8. During two
orbital precession periods, the relative deviation of the analytic
approximation from the exact solution is $\lesssim 10^{-9}$ for $\omega_{1,2}$
and $\lesssim 10^{-13}$ for $\omega_{3}$. Therefore, the solution of the
angular frequency is accurate enough for the calculation of the waveforms.

Figure 8: The absolute errors between the approximate angular frequency
$\omega_{i}$ and the exact angular frequency $\omega_{i}^{E}$ during two
orbital precession periods. The parameters used here are the same as in Fig.
3.

##  Appendix D Effects of orbital precession on Doppler shift

The Doppler shift correction is defined as

$\delta(\Delta\Omega_{n})=\Delta\Omega_{n}-\Delta\Omega_{n}(\Omega_{\rm
pre}=0)$ (59)

with $n=\\{\rm r,p,pre\\}$. It measures the effects of spin-orbit coupling on
the phase of GWs of the spinning NS in a binary. As seen in Fig. 9, the
deviation of the GW frequency can reach $\sim 0.5\%$ for $f_{0}$ and $\sim
1\%$ for $2f_{0}$ if we do not consider the orbital plane precession.

Figure 9: The Doppler shift corrections to the different frequency components
in the waveforms caused by orbital plane precession during one orbital
precession period. For visualisation purposes, $\delta(\Delta\Omega_{\rm p})$
and $\delta(\Delta\Omega_{\rm pre})$ are magnified by a factor of $10^{5}$ and
$10^{8}$, respectively. The parameters used here are the same as in Fig. 3.

## References

  * Sieniawska and Bejger [2019] M. Sieniawska and M. Bejger, Continuous Gravitational Waves from Neutron Stars: Current Status and Prospects, Universe 5, 217 (2019), arXiv:1909.12600 [astro-ph.HE] .
  * Riles [2017] K. Riles, Recent searches for continuous gravitational waves, Modern Physics Letters A 32, 1730035-685 (2017), arXiv:1712.05897 [gr-qc] .
  * Lasky [2015] P. D. Lasky, Gravitational Waves from Neutron Stars: A Review, Publications of the Astronomical Society of Australia 32, e034 (2015), arXiv:1508.06643 [astro-ph.HE] .
  * Riles [2022] K. Riles, Searches for Continuous-Wave Gravitational Radiation, arXiv e-prints , arXiv:2206.06447 (2022), arXiv:2206.06447 [astro-ph.HE] .
  * Wette [2023] K. Wette, Searches for continuous gravitational waves from neutron stars: A twenty-year retrospective, arXiv e-prints , arXiv:2305.07106 (2023), arXiv:2305.07106 [gr-qc] .
  * Pitkin [2011] M. Pitkin, Prospects of observing continuous gravitational waves from known pulsars, Monthly Notices of the Royal Astronomical Society 415, 1849 (2011), arXiv:1103.5867 [astro-ph.HE] .
  * Soldateschi and Bucciantini [2021] J. Soldateschi and N. Bucciantini, Detectability of Continuous Gravitational Waves from Magnetically Deformed Neutron Stars, Galaxies 9, 101 (2021), arXiv:2110.06039 [astro-ph.HE] .
  * Lu  _et al._ [2022a] N. Lu, K. Wette, S. M. Scott, and A. Melatos, Inferring neutron star properties with continuous gravitational waves, arXiv e-prints , arXiv:2209.10981 (2022a), arXiv:2209.10981 [gr-qc] .
  * LIGO Scientific Collaboration [2015] LIGO Scientific Collaboration, Advanced LIGO, Classical and Quantum Gravity 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ [2015] F. Acernese, M. Agathos, _et al._ , Advanced Virgo: a second-generation interferometric gravitational wave detector, Classical and Quantum Gravity 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * The LIGO Scientific Collaboration  _et al._ [2022] The LIGO Scientific Collaboration, the Virgo Collaboration, and the KAGRA Collaboration, All-sky search for continuous gravitational waves from isolated neutron stars using Advanced LIGO and Advanced Virgo O3 data, arXiv e-prints , arXiv:2201.00697 (2022), arXiv:2201.00697 [gr-qc] .
  * Abbott  _et al._ [2021] R. Abbott, T. D. Abbott, and et al. (LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration), All-sky search for continuous gravitational waves from isolated neutron stars in the early o3 ligo data, Phys. Rev. D 104, 082004 (2021).
  * Abbott  _et al._ [2022a] R. Abbott, T. D. Abbott, and et al., Narrowband Searches for Continuous and Long-duration Transient Gravitational Waves from Known Pulsars in the LIGO-Virgo Third Observing Run, Astrophys. J. 932, 133 (2022a), arXiv:2112.10990 [gr-qc] .
  * Abbott  _et al._ [2022b] R. Abbott, T. D. Abbott, and et al., Search of the early O3 LIGO data for continuous gravitational waves from the Cassiopeia A and Vela Jr. supernova remnants, Phys. Rev. D 105, 082005 (2022b), arXiv:2111.15116 [gr-qc] .
  * Abbott  _et al._ [2021] R. Abbott, T. D. Abbott, S. Abraham, _et al._ , Searches for Continuous Gravitational Waves from Young Supernova Remnants in the Early Third Observing Run of Advanced LIGO and Virgo, Astrophys. J. 921, 80 (2021), arXiv:2105.11641 [astro-ph.HE] .
  * Sun  _et al._ [2019] L. Sun, A. Melatos, and P. D. Lasky, Tracking continuous gravitational waves from a neutron star at once and twice the spin frequency with a hidden Markov model, Phys. Rev. D 99, 123010 (2019), arXiv:1903.03866 [astro-ph.IM] .
  * Jaranowski  _et al._ [1998] P. Jaranowski, A. Krolak, and B. F. Schutz, Data analysis of gravitational-wave signals from spinning neutron stars: The signal and its detection, Phys. Rev. D 58, 063001 (1998), arXiv:gr-qc/9804014 [gr-qc] .
  * Abbott  _et al._ [2022c] R. Abbott, H. Abe, and et al., Searches for Gravitational Waves from Known Pulsars at Two Harmonics in the Second and Third LIGO-Virgo Observing Runs, Astrophys. J. 935, 1 (2022c), arXiv:2111.13106 [astro-ph.HE] .
  * Abbott  _et al._ [2019] B. P. Abbott, R. Abbott, _et al._ , Searches for Gravitational Waves from Known Pulsars at Two Harmonics in 2015-2017 LIGO Data, Astrophys. J. 879, 10 (2019), arXiv:1902.08507 [astro-ph.HE] .
  * Pitkin  _et al._ [2015] M. Pitkin, C. Gill, D. I. Jones, G. Woan, and G. S. Davies, First results and future prospects for dual-harmonic searches for gravitational waves from spinning neutron stars, Monthly Notices of the Royal Astronomical Society 453, 4399 (2015), arXiv:1508.00416 [astro-ph.HE] .
  * Jones [2015] D. I. Jones, Parameter choices and ranges for continuous gravitational wave searches for steadily spinning neutron stars, Monthly Notices of the Royal Astronomical Society 453, 53 (2015), arXiv:1501.05832 [gr-qc] .
  * Jones [2010] D. I. Jones, Gravitational wave emission from rotating superfluid neutron stars, Monthly Notices of the Royal Astronomical Society 402, 2503 (2010), arXiv:0909.4035 [astro-ph.SR] .
  * Gal’Tsov and Tsvetkov [1984] D. V. Gal’Tsov and V. P. Tsvetkov, On the gravitational radiation of an oblique rotator, Physics Letters A 103, 193 (1984).
  * Bonazzola and Gourgoulhon [1996] S. Bonazzola and E. Gourgoulhon, Gravitational waves from pulsars: emission by the magnetic-field-induced distortion., Astronomy and Astrophysics 312, 675 (1996), arXiv:astro-ph/9602107 [astro-ph] .
  * Lindblom  _et al._ [1998] L. Lindblom, B. J. Owen, and S. M. Morsink, Gravitational radiation instability in hot young neutron stars, Phys. Rev. Lett. 80, 4843 (1998).
  * Andersson [1998] N. Andersson, A New Class of Unstable Modes of Rotating Relativistic Stars, Astrophys. J. 502, 708 (1998), arXiv:gr-qc/9706075 [gr-qc] .
  * Zimmermann [1980] M. Zimmermann, Gravitational waves from rotating and precessing rigid bodies. II. General solutions and computationally useful formulas, Phys. Rev. D 21, 891 (1980).
  * Van Den Broeck [2005] C. Van Den Broeck, The gravitational wave spectrum of non-axisymmetric, freely precessing neutron stars, Classical and Quantum Gravity 22, 1825 (2005), arXiv:gr-qc/0411030 [gr-qc] .
  * Gao  _et al._ [2020] Y. Gao, L. Shao, R. Xu, L. Sun, C. Liu, and R.-X. Xu, Triaxially deformed freely precessing neutron stars: continuous electromagnetic and gravitational radiation, Monthly Notices of the Royal Astronomical Society 498, 1826 (2020), arXiv:2007.02528 [astro-ph.HE] .
  * Covas and Sintes [2019] P. B. Covas and A. M. Sintes, New method to search for continuous gravitational waves from unknown neutron stars in binary systems, Phys. Rev. D 99, 124019 (2019), arXiv:1904.04873 [astro-ph.IM] .
  * ATN [2021] ATNF Pulsar Catalogue, https://www.atnf.csiro.au/research/pulsar/psrcat/ (2021), [last accessed 1-Mar-2023].
  * Manchester  _et al._ [2005] R. N. Manchester, G. B. Hobbs, A. Teoh, and M. Hobbs, The Australia Telescope National Facility Pulsar Catalogue, The Astronomical Journal 129, 1993 (2005), arXiv:astro-ph/0412641 [astro-ph] .
  * Abbott  _et al._ [2021] R. Abbott, T. D. Abbott, and et al. (The LIGO Scientific Collaboration and the Virgo Collaboration), All-sky search in early o3 ligo data for continuous gravitational-wave signals from unknown neutron stars in binary systems, Phys. Rev. D 103, 064017 (2021).
  * Covas and Sintes [2020] P. B. Covas and A. M. Sintes, First All-Sky Search for Continuous Gravitational-Wave Signals from Unknown Neutron Stars in Binary Systems Using Advanced LIGO Data, Phys. Rev. Lett. 124, 191102 (2020), arXiv:2001.08411 [gr-qc] .
  * Zhang  _et al._ [2021] Y. Zhang, M. A. Papa, B. Krishnan, and A. L. Watts, Search for Continuous Gravitational Waves from Scorpius X-1 in LIGO O2 Data, The Astrophysical Journal Letters 906, L14 (2021), arXiv:2011.04414 [astro-ph.HE] .
  * Leaci and Prix [2015] P. Leaci and R. Prix, Directed searches for continuous gravitational waves from binary systems: Parameter-space metrics and optimal Scorpius X-1 sensitivity, Phys. Rev. D 91, 102003 (2015), arXiv:1502.00914 [gr-qc] .
  * Andrews  _et al._ [2020] J. J. Andrews, K. Breivik, C. Pankow, D. J. D’Orazio, and M. Safarzadeh, LISA and the Existence of a Fast-merging Double Neutron Star Formation Channel, The Astrophysical Journal Letters 892, L9 (2020), arXiv:1910.13436 [astro-ph.HE] .
  * Feng  _et al._ [2023] W.-F. Feng, J.-W. Chen, Y. Wang, S. D. Mohanty, and Y. Shao, Multimessenger observations of double neutron stars in the galactic disk with gravitational and radio waves, Phys. Rev. D 107, 103035 (2023).
  * Abbott  _et al._ [2019] B. Abbott, R. Abbott, _et al._ , GWTC-1: a gravitational-wave transient catalog of compact binary mergers observed by ligo and virgo during the first and second observing runs, Physical Review X 9, 031040 (2019).
  * Srivastava  _et al._ [2022] V. Srivastava, D. Davis, K. Kuns, P. Landry, S. Ballmer, M. Evans, E. D. Hall, J. Read, and B. S. Sathyaprakash, Science-driven Tunable Design of Cosmic Explorer Detectors, Astrophys. J. 931, 22 (2022), arXiv:2201.10668 [gr-qc] .
  * Punturo  _et al._ [2010] M. Punturo, M. Abernathy, and et al., The Einstein Telescope: a third-generation gravitational wave observatory, Classical and Quantum Gravity 27, 194002 (2010).
  * Davis and Goldstein [1970] L. Davis and M. Goldstein, Magnetic-Dipole Alignment in Pulsars, The Astrophysical Journal Letters 159, L81 (1970).
  * Bildsten and Cutler [1992] L. Bildsten and C. Cutler, Tidal Interactions of Inspiraling Compact Binaries, Astrophys. J. 400, 175 (1992).
  * Apostolatos  _et al._ [1994] T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, Spin-induced orbital precession and its modulation of the gravitational waveforms from merging binaries, Phys. Rev. D 49, 6274 (1994).
  * Landau and Lifshitz [1976] L. Landau and E. Lifshitz, Chapter vi - motion of a rigid body, in  _Mechanics (Third Edition)_, edited by L. Landau and E. Lifshitz (Butterworth-Heinemann, Oxford, 1976) third edition ed., pp. 96–130.
  * Abramowitz  _et al._ [1965] M. Abramowitz, I. A. Stegun, and D. M. Miller, Handbook of mathematical functions with formulas, graphs and mathematical tables (national bureau of standards applied mathematics series no. 55), Journal of Applied Mechanics 32, 239 (1965).
  * Tauris  _et al._ [2017] T. M. Tauris, M. Kramer, P. C. C. Freire, N. Wex, H. T. Janka, N. Langer, P. Podsiadlowski, E. Bozzo, S. Chaty, M. U. Kruckow, E. P. J. van den Heuvel, J. Antoniadis, R. P. Breton, and D. J. Champion, Formation of Double Neutron Star Systems, Astrophys. J. 846, 170 (2017), arXiv:1706.09438 [astro-ph.HE] .
  * Shao and Li [2018] Y. Shao and X.-D. Li, On the Role of Supernova Kicks in the Formation of Galactic Double Neutron Star Systems, Astrophys. J. 867, 124 (2018), arXiv:1810.03324 [astro-ph.HE] .
  * Poisson and Will [2014] E. Poisson and C. M. Will,  _Gravity: Newtonian, post-newtonian, relativistic_ (Cambridge University Press, 2014).
  * Wagg  _et al._ [2022] T. Wagg, F. S. Broekgaarden, S. E. de Mink, N. Frankel, L. A. C. van Son, and S. Justham, Gravitational Wave Sources in Our Galactic Backyard: Predictions for BHBH, BHNS, and NSNS Binaries Detectable with LISA, Astrophys. J. 937, 118 (2022), arXiv:2111.13704 [astro-ph.HE] .
  * Maggiore [2007] M. Maggiore,  _Gravitational waves: Volume 1: Theory and experiments_ (OUP Oxford, 2007).
  * Zimmermann and Szedenits [1979] M. Zimmermann and E. Szedenits, Gravitational waves from rotating and precessing rigid bodies: Simple models and applications to pulsars, Phys. Rev. D 20, 351 (1979).
  * Demorest  _et al._ [2010] P. B. Demorest, T. Pennucci, S. M. Ransom, M. S. E. Roberts, and J. W. T. Hessels, A two-solar-mass neutron star measured using Shapiro delay, Nature (London) 467, 1081 (2010), arXiv:1010.5788 [astro-ph.HE] .
  * Abbott  _et al._ [2020a] B. P. Abbott, R. Abbott, and et al., GW190425: Observation of a Compact Binary Coalescence with Total Mass $\sim$ 3.4 M⊙, The Astrophysical Journal Letters 892, L3 (2020a), arXiv:2001.01761 [astro-ph.HE] .
  * Lu  _et al._ [2022b] N. Lu, K. Wette, S. M. Scott, and A. Melatos, Inferring neutron star properties with continuous gravitational waves, arXiv e-prints , arXiv:2209.10981 (2022b), arXiv:2209.10981 [gr-qc] .
  * Abbott  _et al._ [2020b] R. Abbott, T. D. Abbott, _et al._ , Gravitational-wave Constraints on the Equatorial Ellipticity of Millisecond Pulsars, The Astrophysical Journal Letters 902, L21 (2020b), arXiv:2007.14251 [astro-ph.HE] .
  * Dergachev and Papa [2020] V. Dergachev and M. A. Papa, Results from the First All-Sky Search for Continuous Gravitational Waves from Small-Ellipticity Sources, Phys. Rev. Lett. 125, 171101 (2020), arXiv:2004.08334 [gr-qc] .
  * Lau  _et al._ [2020] M. Y. M. Lau, I. Mandel, A. Vigna-Gómez, C. J. Neijssel, S. Stevenson, and A. Sesana, Detecting double neutron stars with LISA, Monthly Notices of the Royal Astronomical Society 492, 3061 (2020), arXiv:1910.12422 [astro-ph.HE] .
  * Apostolatos [1995] T. A. Apostolatos, Search templates for gravitational waves from precessing, inspiraling binaries, Phys. Rev. D 52, 605 (1995).
  * Creighton and Anderson [2012] J. D. Creighton and W. G. Anderson,  _Gravitational-wave physics and astronomy: An introduction to theory, experiment and data analysis_ (John Wiley & Sons, 2012).
  * Jaranowski and Królak [1999] P. Jaranowski and A. Królak, Data analysis of gravitational-wave signals from spinning neutron stars. II. Accuracy of estimation of parameters, Phys. Rev. D 59, 063003 (1999), arXiv:gr-qc/9809046 [gr-qc] .
  * Shah  _et al._ [2012] S. Shah, M. van der Sluys, and G. Nelemans, Using electromagnetic observations to aid gravitational-wave parameter estimation of compact binaries observed with LISA, Astronomy and Astrophysics 544, A153 (2012), arXiv:1207.6770 [astro-ph.IM] .
  * Cutler and Vecchio [1998] C. Cutler and A. Vecchio, Lisa’s angular resolution for monochromatic sources, in  _AIP Conference Proceedings_ , Vol. 456 (American Institute of Physics, 1998) pp. 95–100.
  * Tauris [2018] T. M. Tauris, Disentangling Coalescing Neutron-Star-White-Dwarf Binaries for LISA, Phys. Rev. Lett. 121, 131105 (2018), arXiv:1809.03504 [astro-ph.SR] .
  * Chen [2021] W.-C. Chen, Ultracompact binary pulsars as continuous dual-line gravitational wave sources, Phys. Rev. D 103, 103004 (2021), arXiv:2104.05910 [astro-ph.HE] .