# Perturbative Frequency Expansion for Nearly Monochromatic Binary Black Holes
Detectable with LISA

Naoki Seto  Department of Physics, Kyoto University, Kyoto 606-8502, Japan

###### Abstract

The proposed space gravitational wave (GW) detector LISA has potential to
detect stellar-mass black hole binaries (BBHs). The majority of the detected
BBHs are expected to emit nearly monochromatic GWs, whose frequency evolution
will be efficiently described by Taylor expansions. We study the measurability
of the associated time derivative coefficients of the frequencies, by
extending a recent work based on a simplified Fisher matrix analysis.
Additionally, we provide qualitative discussions on how to extract
astrophysical information, such as orbital eccentricity and tertiary
perturbation, from the observed derivative coefficients.

###### pacs:

PACS number(s): 95.55.Ym 98.80.Es,95.85.Sz

##  I introduction

The LIGO-Virgo-Kagra network had detected $O(100)$ binary black holes (BBHs)
by the end of its O3 run KAGRA:2021vkt . The proposed space GW detectors such
as LISA LISA:2017pwj , Taiji Ruan:2018tsw  and TianQin TianQin:2015yph  can
observe similar BBHs in the lower frequency regime 2022arXiv220306016A ;
Sesana:2016ljz .

According to a recent statistical analysis employing the LVK GWTC-3 catalog,
LISA is expected to detect $N_{\rm tot}\sim 2.3(T/{\rm 4yr})^{3/2}(\rho_{\rm
thr}/10)^{-3}$ BBHs similar to those cataloged Seto:2022xmh  (see also
Gerosa:2019dbe ; Klein:2022rbf  for multiband observations). Here $\rho_{\rm
thr}$ is the detection threshold for the signal-to-noise ratio, and $T$ is the
observational period. The estimated number $N_{\rm tot}$ contains an
uncertainty of a factor of two, partly due to the unclear redshift dependence
of the merger rate. In the statistical analysis, the BBHs are simply assumed
to be nearly circular still in the LISA band.

In the total number $N_{\rm tot}$, the fraction of BBHs that will merge in the
observational period $T=4$yr was estimated to be $\sim 5\%$ ($\sim 10\%$ for
$T=10$yr) Seto:2022xmh . In fact, the majority of the detected BBHs will be
nearly monochromatic GW sources, with the peak number density (per logarithmic
frequency interval) around $f\sim 5$mHz. Here the characteristic frequency
$\sim 5$mHz is given as the tangential point of the noise spectrum $S_{n}(f)$
with the slope $f^{-2/9}$ which is derived from the radiation reaction. In
Seto:2022xmh , the median chirp mass of the detected BBHs was also estimated
to be ${\cal M}\sim 30M_{\odot}$ with a relatively large uncertainty for the
distribution function at ${\cal M}\gtrsim 50M_{\odot}$.

For the nearly monochromatic GWs from these BBHs, we will be able to
efficiently describe the frequency evolution, using the time derivative
coefficients ${\dot{f}},{\ddot{f}},\cdots$ (also denoting
$f^{(n)}=\partial_{t}^{n}f$), as conventionally applied in similar occasions
(e.g. spin-down of a pulsar). This approach condenses the evolutionary
information to a finite number of parameters and will be mathematically more
tractable, rather than directly dealing with the functional degrees of
freedom. Meanwhile, stimulated by the recent reports on the long-term orbital
evolution of HM Cancri (an interacting double white dwarf binary at
$f=6.2$mHz) 2021ApJ…912L…8S ; 2023MNRAS.518.5123M , the author developed a
simplified Fisher matrix analysis on the estimation errors of the frequency
derivatives up to the second order ${\ddot{f}}$ Seto:2023skl  (see also
Ebadi:2024oaq ). In this paper, extending the analytical approach to include
higher derivatives such as ${\dddot{f}}$ and beyond, we apply it to the nearly
monochromatic BBHs in the LISA band.

Importantly, compared with the BBHs in the LVK band, those in the LISA band
could have much larger eccentricities Peters:1964zz ; Seto:2016wom ;
Nishizawa:2016eza . Likewise, in the LISA band, the environmental effects
(e.g. tertiary perturbation) could be more strongly imprinted in the emitted
GWs, because of the longer gravitational radiation timescale Seto:2008di ;
Seto:2016zfy ; Tamanini:2018cqb ; Wong:2019hsq ; Kang:2021bmp ; Robson:2018svj
; Tamanini:2019usx ; Xuan, Peng, & Chen (2021); Zhang:2021pwe ; Maeda:2023tao
. In this paper, we also discuss the possibilities of examining the orbital
eccentricities and the potential tertiary perturbations from the observables
$f,{\dot{f}},{\ddot{f}},{\dddot{f}},\cdots$ (see also Takahashi:2023flk ). To
this end, we introduce the non-dimensional parameters $C_{n}\propto
f^{(n)}f^{n-1}/({\dot{f}})^{n}$ ($n=2,3,\cdots$) which are adjusted to be
$C_{n}=1$ for an isolated circular BBH. We apply Fisher matrix results to
evaluate the measurement errors for the parameters $C_{n}$.

To provide a concrete picture for our study, we consider an equal mass BBH
with the characteristic values mentioned earlier: a chirp mass ${\cal
M}=30M_{\odot}$, a GW frequency $f=$5mHz and a signal strength $\rho=10$ (for
$T=4$yr with LISA). We set this BBH as our fiducial observational target.

This paper is organized as follows. In section II, we evaluate the frequency
derivatives $f^{(n)}$ for circular and eccentric BBHs. We introduce the non-
dimensional measures $C_{n}$ and qualitatively argue the corrections due to
tertiary perturbations. In section III, based on the Fisher matrix analysis,
we present analytical expressions for the estimation errors of the time
derivatives $f^{(n)}$. In section IV, we discuss the measurability of these
derivative coefficients. We relatedly discuss the prospects for studying
orbital eccentricity and tertiary effects by using the quantities $C_{n}$. In
section V, we mention potential extensions of this study. Section VI is a
brief summary of this paper.

##  II evolution of binary black holes

###  II.1 Circular Systems

Let us first discuss the orbital evolution of a circular BBH. In the
quadrupole approximation, the chirp rate of the GW frequency $f$ (twice the
orbital frequency) is given by

$\frac{df}{dt}=\frac{96\pi^{8/3}G^{5/3}{\cal M}^{5/3}f^{11/3}}{5c^{5}}$ (1)

with the chirp mass ${\cal M}$, the gravitational constant $G$ and the speed
of light $c$ Peters:1964zz . For our fiducial model parameters (summarized in
Table 1) we have

$\displaystyle{\dot{f}}=6.2\times 10^{-13}\left(\frac{{\cal
M}}{30M_{\odot}}\right)^{5/3}\left(\frac{f}{\rm 5mHz}\right)^{11/3}{\rm
Hz\,s^{-1}}.$ (2)

Correspondingly, in the observational period $T$, the GW frequency $f$
increases only by $1.5(T/{\rm 4yr})$ percent.

Table 1: Fiducial binary parameters and the associated frequency derivatives due to the radiation reaction for a circular orbit.  frequency $f$ | 5 mHz
---|---
chirp mass ${\cal M}$ | 30 $M_{\odot}$
SNR $\rho$ (4yr) | 10
$\dot{f}$ | $6.2\times 10^{-13}\rm Hz~{}s^{-1}$
$\ddot{f}$ | $2.8\times 10^{-22}\rm Hz~{}s^{-2}$
$\dddot{f}$ | $2.2\times 10^{-31}\rm Hz~{}s^{-3}$
$\dddot{f}$ | $2.4\times 10^{-40}\rm Hz~{}s^{-4}$

Recursively taking the time derivatives of Eq. (1), we can readily obtain the
higher time derivatives such as

$\displaystyle{\ddot{f}}=2.8\times 10^{-22}\left(\frac{{\cal
M}}{30M_{\odot}}\right)^{10/3}\left(\frac{f}{\rm 5mHz}\right)^{19/3}{\rm
Hz\,s^{-2}},$ (3) $\displaystyle{\dddot{f}}=2.2\times
10^{-31}\left(\frac{{\cal M}}{30M_{\odot}}\right)^{5}\left(\frac{f}{\rm
5mHz}\right)^{9}{\rm Hz\,s^{-3}},$ (4) $\displaystyle{\ddddot{f}}=2.4\times
10^{-40}\left(\frac{{\cal M}}{30M_{\odot}}\right)^{20/3}\left(\frac{f}{\rm
5mHz}\right)^{35/3}{\rm Hz\,s^{-4}}.$ (5)

By integrating Eq. (1), we can also estimate the time before the merger as

$\displaystyle t_{c}=97\left(\frac{{\cal
M}}{30M_{\odot}}\right)^{-5/3}\left(\frac{f}{\rm 5mHz}\right)^{-8/3}{\rm yr}.$
(6)

In terms of $t_{c}$, the time derivatives are roughly given by

$\partial^{n}_{t}f=f^{(n)}\sim\frac{f}{t_{c}^{n}}.$ (7)

More precisely, from Eq. (1), we can show

$\displaystyle C_{2}\equiv\frac{3{\ddot{f}}f}{11{\dot{f}}^{2}}$
$\displaystyle=$ $\displaystyle 1,$ (8) $\displaystyle
C_{3}\equiv\frac{9{\dddot{f}}f^{2}}{209{\dot{f}}^{3}}$ $\displaystyle=$
$\displaystyle 1,$ (9) $\displaystyle
C_{4}\equiv\frac{{\ddddot{f}}f^{3}}{209{\dot{f}}^{4}}$ $\displaystyle=$
$\displaystyle 1.$ (10)

Up to now, we ignored the post-Newtonian effects. The leading order (1PN)
correction to the product $C_{2}$ takes a negative value Blanchet:2013haa ;
Cutler:1994ys

$\displaystyle-\frac{2}{11}\left(\frac{743}{336}+\frac{11\eta}{4}\right)\left(\frac{GM_{t}\pi
f}{c^{3}}\right)^{2/3}$ (11)

with the symmetric mass ratio $\eta$ and the total mass
$M_{t}=\eta^{-3/5}{\cal M}$. In the LISA era, from the preceding ground based
detectors including the Einstein Telescope Punturo:2010zz , we will be able to
use the detailed prior information on the ratio $\eta$ for stellar mass BBHs
and better model the PN corrections. For an equal mass binary (with
$\eta=1/4$), Eq. (11) is expressed as

$\displaystyle-1.6\times 10^{-4}\left(\frac{{\cal
M}}{30M_{\odot}}\right)^{2/3}\left(\frac{f}{\rm 5mHz}\right)^{2/3}.$ (12)

###  II.2 Eccentric Systems

Next we discuss the time evolution of an eccentric BBH (see also Tanay:2016zog
). We continue to use $f$ as twice the orbital frequency. At the quadrupole
order, using the orbital eccentricity $e$, we can write down the evolution
equations Peters:1964zz

$\displaystyle\frac{df}{dt}$ $\displaystyle=$
$\displaystyle\frac{96\pi^{8/3}G^{5/3}{\cal
M}^{5/3}f^{11/3}}{5c^{5}\left(1-e^{2}\right)^{7/2}}{\left(1+\frac{73e^{2}}{24}+\frac{37e^{4}}{96}\right)}$
(13) $\displaystyle\equiv$ $\displaystyle F_{1}(f,e),$ (14)
$\displaystyle\frac{de}{dt}$ $\displaystyle=$
$\displaystyle-\frac{304\pi^{8/3}G^{5/3}f^{8/3}{\cal
M}^{5/3}e}{15c^{5}(1-e^{2})^{5/2}}\left(1+\frac{121e^{2}}{304}\right)$ (15)
$\displaystyle\equiv$ $\displaystyle E_{1}(f,e).$ (16)

Applying the chain rule for derivatives, we can evaluate the second derivative
$\ddot{f}$ as follows

$\displaystyle{\ddot{f}}=\frac{\partial F_{1}}{\partial f}F_{1}+\frac{\partial
F_{1}}{\partial e}E_{1}\equiv F_{2}(f,e).$ (17)

The product $C_{2}$ in Eq. (8) is now given as

$\displaystyle C_{2}$ $\displaystyle=$
$\displaystyle\frac{3{\ddot{f}}f}{11{\dot{f}}^{2}}$ $\displaystyle=$
$\displaystyle\frac{4\left(407e^{8}+2344e^{6}+93810e^{4}+58720e^{2}+25344\right)}{11\left(37e^{4}+292e^{2}+96\right)^{2}}$
$\displaystyle\equiv$ $\displaystyle{{K}}_{2}(e)$ (19)

with a perturbative expression

$\displaystyle{{K}}_{2}(e)=1-\frac{2983e^{2}}{792}+O(e^{4}).$ (20)

In contrast to Eqs. (13) and (15), the function ${{K}}_{2}(e)$ is regular at
$e=1$ with  ${{K}}_{2}(1)=4/11$.

In the same manner as Eq. (17), we can evaluate the third derivative
$\dddot{f}$ for an eccentric orbit as

$\displaystyle{\dddot{f}}=\frac{\partial F_{2}}{\partial
f}F_{1}+\frac{\partial F_{2}}{\partial e}E_{1}\equiv F_{3}(f,e).$ (21)

We then evaluate the product $C_{3}$ and define the function ${K}_{3}(e)$

$\displaystyle{C}_{3}$ $\displaystyle=$
$\displaystyle\frac{9{\dddot{f}}f^{2}}{209{\dot{f}}^{3}}$ $\displaystyle=$
$\displaystyle\frac{4}{209\left(37e^{4}+292e^{2}+96\right)^{3}}(-9361e^{12}+500532e^{10}$
$\displaystyle-42206966e^{8}+145890160e^{6}+129278336e^{4}$
$\displaystyle+104147968e^{2}+46227456)$ $\displaystyle\equiv$
$\displaystyle{{K}}_{3}(e)$ (23)

with ${{K}}_{3}(0)=1$, ${{K}}_{3}(1)=20/209$ and

$\displaystyle{K}_{3}(e)=1-\frac{2041e^{2}}{297}+O(e^{4}).$ (24)

Likewise, we can evaluate the fourth derivative $\ddddot{f}(\equiv f^{(4)})$
and the product

$\displaystyle C_{4}=\frac{{\ddddot{f}}f^{3}}{209{\dot{f}}^{4}}.$ (25)

The analytic expression for the corresponding function ${K}_{4}(e)$ is given
as

$\displaystyle{K}_{4}(e)$ $\displaystyle=$
$\displaystyle\frac{8}{5643(37e^{4}+292e^{2}+96)^{4}}(421245e^{16}$ (26)
$\displaystyle-100387440e^{14}+12951727630e^{12}$
$\displaystyle-122186028060e^{10}+117746520816e^{8}$
$\displaystyle+132600265216e^{6}+150472180736e^{4}$
$\displaystyle+137985376256e^{2}+59910782976)$

with ${K}_{4}(0)=1$, ${K}_{4}(1)=40/1881$ and

$\displaystyle K_{4}(e)=1-\frac{1898287e^{2}}{192456}+O(e^{4}).$ (27)

In Fig. 1, we present the functions ${K}_{2}(e)$, ${K}_{3}(e)$ and
${K}_{4}(e)$. All of them decrease monotonically with $e$. By evaluating the
products $C_{n}$ from the observed frequency derivatives, we might inversely
estimate the orbital eccentricity $e$.

Figure 1: The three functions ${{K}}_{2}(e)$, ${K}_{3}(e)$ and ${K}_{4}(e)$
defined in Eqs. (19) (23) and (26). For a circular binary ($e=0$), we have
${K}_{2}={K}_{3}={K}_{4}=1$. At $e=0.1$, we have ${K}_{2}=0.96$,
${K}_{3}=0.93$ and ${K}_{4}=0.91$.

The PN corrections to the products $C_{n}$ depend on the eccentricity $e$. In
this paper, we approximately use the results (12) also for mildly eccentric
BBHs.

###  II.3 Tertiary Effects

Here, we briefly discuss the GW phase modulation induced by a tertiary (of
mass $m_{s}$) around a BBH (see, e.g. Seto:2023rfh ; Bassa:2016fiy ;
Kaplan:2016ymq  for the effects on the derivative coefficients $f^{(n)}$), as
a representative environmental perturbation. In this subsection, we
temporarily attach the subscript $s$ to the tertiary related quantities and
$r$ to the radiational effects (under the approximation ignoring the cross
terms Seto:2023rfh ). For simplicity, we assume that the outer orbits is
hierarchical and circular with the period $P_{s}$ and the semimajor axis
$R_{s}$. We have

$R_{s}=\left[\frac{G(M_{t}+m_{s})P_{s}^{2}}{4\pi^{2}}\right]^{1/3}.$ (28)

from the third Kepler law. In this simplified case, the BBH moves
approximately on a circular orbit around the barycenter of the triple system.

Following Seto:2023rfh , we briefly discuss the tertiary corrections
(equivalent to the Rømer effects) to the apparent GW frequency and its
derivative coefficients. The positional variation induces the phase shift as

$\Phi_{s}(t)=-2\pi f_{r}D(t)/c$ (29)

where $D(t)$ is the line-of-sight distance variation and expressed as
($\phi_{s}$: the initial outer orbital phase)

$D(t)=D_{s}\cos[2\pi t/P_{s}+\phi_{s}]$ (30)

with the projected distance

$D_{s}=\frac{m_{s}}{M_{t}+m_{s}}R_{s}\sin I_{s}$ (31)

($I_{s}$: the outer inclination angle). The tertiary corrections to the
apparent GW frequency and its time derivatives are given as

$f_{s}=\frac{1}{2\pi}\frac{d\Phi_{s}}{dt},~{}~{}{\dot{f}}_{s}=\frac{1}{2\pi}\frac{d^{2}\Phi_{s}}{dt^{2}},~{}~{}{\ddot{f}}_{s}=\frac{1}{2\pi}\frac{d^{3}\Phi_{s}}{dt^{3}}$
(32)

and so on.

We put $v_{s}\sim 2\pi D_{s}/P_{s}$ as the line-of-sight component of the
outer circular velocity of the binary. Then, for an observational time $T(\ll
P_{s})$, the tertiary effects $f_{s}^{(n)}$ have the characteristic magnitudes

$\displaystyle|f_{s}^{(n)}|\sim
f\left(\frac{v_{s}}{c}\right)\left(\frac{2\pi}{P_{s}}\right)^{n}.$ (33)

A typical field system will have $v_{s}/c\ll 1$. Depending on the outer
orbital phase, the tertiary corrections $f_{s}^{(n)}$ can take both signs.

Now we compare the tertiary and radiational effects, using Eqs. (7) and (33).
We can roughly evaluate the ratios

$\displaystyle\frac{|f_{s}^{(n)}|}{f_{r}^{(n)}}\sim\left(\frac{v_{s}}{c}\right)\left(\frac{2\pi
t_{c}}{P_{s}}\right)^{n}$ (34)

with the time $t_{c}$ before the merger. For $2\pi t_{c}/P_{s}\gg 1$, the
tertiary effects could be relatively more important at larger $n$. However, in
contrast to Galactic white dwarf binaries ($t_{c}\gtrsim 10^{5}$yr in the LISA
band), our BBHs have much shorter merger time $t_{c}\sim 100$yr, and the
amplification factor $\left({2\pi t_{c}}/{P_{s}}\right)$ will work less
efficiently in Eq. (34).

As a concrete reference model, we set the tertiary mass $M_{s}\sim
25M_{\odot}$ and an edge-on outer orbit ($I_{s}=\pi/2$) of $P_{s}=60$yr ($\ll
2\pi t_{c}$). Then, using Eqs. (29)-(32), we obtain

$\displaystyle\frac{{f}_{s}}{{f}_{r}}$ $\displaystyle=$ $\displaystyle
1.1\times 10^{-4}\sin\psi_{s},$ (35)
$\displaystyle\frac{{\dot{f}}_{s}}{{\dot{f}}_{r}}$ $\displaystyle=$
$\displaystyle 8.2\times 10^{-4}\cos\psi_{s},$ (36)
$\displaystyle\frac{{\ddot{f}}_{s}}{{\ddot{f}}_{r}}$ $\displaystyle=$
$\displaystyle-6.0\times 10^{-3}\sin\psi_{s},$ (37)
$\displaystyle\frac{{\dddot{f}}_{s}}{{\dddot{f}}_{r}}$ $\displaystyle=$
$\displaystyle-2.6\times 10^{-2}\cos\psi_{s}$ (38)

with the outer orbital phase parameter $\psi_{s}=2\pi t/P_{s}+\phi_{s}$.

In this paper, we ignore the dynamical interaction between the inner and outer
orbits (see, e.g. Samsing:2024syt ). This would be a reasonable approximation
for our reference tertiary system, at least for the Kozai-Lidov mechanism.

If the tertiary effects are much smaller than the radiational effects (as for
the reference model above), we can perturbatively deal with the former. For
$n=2$, we have

$C_{2}=\frac{3({\ddot{f}}_{r}+{\ddot{f}}_{s})(f_{r}+f_{s})}{11({\dot{f}}_{r}+{\dot{f}}_{s})^{2}}\simeq
K_{2}(e)\left(1+\frac{f_{s}}{f_{r}}-2\frac{\dot{f}_{s}}{\dot{f}_{r}}+\frac{\ddot{f}_{s}}{\ddot{f}_{r}}\right)$
(39)

and similarly

$C_{n}\simeq{K}_{n}(e)\left(1+(n-1)\frac{f_{s}}{f_{r}}-n\frac{\dot{f}_{s}}{\dot{f}_{r}}+\frac{f_{s}^{(n)}}{f_{r}^{(n)}}\right)$
(40)

with the inequalities ${K}_{n}(e)\leq K_{n}(0)=1$ from the definitions of
${K}_{n}(e)$. As shown in Eq. (40) The eccentricity and tertiary effects
degenerate in the product $C_{n}$.

For a nearly circular binary $e\sim 0$, we could have a chance to observe the
excess $C_{n}-1>0$ due to the second factor in Eq. (40), depending on the
outer orbital phase $\psi_{s}$. This inequality cannot be realized for an
isolated eccentric BBH and suggests the existence of a third body, as one of
its causes.

The situation will be more complicated for an eccentric binary, for which the
tertiary effect might be actually more interesting in view of dynamical
processes. In some cases, we will be able to measure the higher order product
$C_{3}$ in addition to the lowest one $C_{2}$. Then the orbital eccentricity
$e$ can be estimated independently for $n=2$ and 3 by $e=K^{-1}_{n}(C_{n})$.
The tertiary effects could induce a mismatch between the inverted values. Note
that we have $K^{\prime}_{n}(0)=0$, and the inversions work less efficiently
for nearly circular binary $e\simeq 0$.

Below, in particular for the expected values for the frequency derivatives
$f^{(n)}$, we mainly deal with the radiational effects, dropping the
subscripts $r$ and $s$ (if unnecessarily). We will return to the tertiary
effect later in Sec. IV D.

##  III fisher matrix analysis

In this section, we apply the Fisher matrix analysis to analytically evaluate
how well we can measure the frequency derivatives in a matched filtering
analysis for a nearly monochromatic BBH.

###  III.1 Simplified Waveform Model

We basically extend the earlier work by the author Seto:2023skl  (see also
Seto:2002dz ; Robson:2018svj ) and apply the simplified waveform model

$\displaystyle h(t)=A\cos[\Phi_{n}(t)]$ (41)

with the constant amplitude $A$ and the phase function

$\displaystyle\Phi_{n}(t)=2\pi\left(ft+\frac{{\dot{f}}t^{2}}{2!}+\frac{{\ddot{f}}t^{3}}{3!}+\cdots+\frac{f^{(n)}t^{n+1}}{(n+1)!}\right)+\varphi.$
(42)

This phase function is perturbatively expanded at the time origin $t=0$. Here
the index $n$ represents the order of the highest frequency derivative
$f^{(n)}$ (different from the associated power index of $t^{n+1}$).

In this paper, we take the time origin $t=0$ (for the expansion (42)) at the
midpoint of the observational period $T$. More specifically, we set our
observational period by $t\in[-T/2,T/2]$ not by $t\in[0,T]$. As we see below,
this choice is for reducing some of the Fisher matrix elements by a block
diagonalization. Applying the Parseval’s theorem for a nearly monochromatic
signal, we obtain the expected signal-to-noise ratios $\rho$ as

$\displaystyle\rho^{2}$ $\displaystyle=$ $\displaystyle\frac{2}{S_{\rm
n}(f)}\int_{-T/2}^{T/2}h(t)h(t)dt$ (43) $\displaystyle=$
$\displaystyle\frac{A^{2}T}{S_{\rm n}(f)}$ (44)

with the effective noise spectrum $S_{n}(f)$. We have ignored the time
dependence of the confusion noise, which will be unimportant at $f\gtrsim
5$mHz and $T\gtrsim 4$yr Robson:2018ifk .

Next, we apply the Fisher matrix analysis to evaluate the measurement errors
for the fitting parameters
${\mbox{\boldmath${\theta}$}}_{n}=\\{\varphi,f,{\dot{f}},\cdots,f^{(n)}\\}$ in
the phase function (42). Note that the error for the amplitude $A$ has no
correlation with those for the parameters ${\mbox{\boldmath${\theta}$}}_{n}$.
The Fisher matrix elements are formally given as

$\displaystyle F_{\theta_{i}\theta_{j}}=\frac{2}{S_{\rm
n}(f)}\int_{-T/2}^{T/2}\partial_{\theta_{i}}h(t)\partial_{\theta_{j}}h(t)dt.$
(45)

We can analytically evaluate these matrix elements and replace the common
factor $A^{2}/S_{n}(f)$ by $\rho^{2}/T$. The covariance matrix
$\left\langle\delta{\theta_{i}}\delta{\theta_{j}}\right\rangle$ for the
measurement errors $\delta\theta_{i}$ can be estimated by the inverse of the
Fisher matrix as

$\displaystyle\left\langle\delta{\theta_{i}}\delta{\theta_{j}}\right\rangle=(F^{-1})_{\theta_{i}\theta_{j}}.$
(46)

Below, for notational simplicity, we put the rms value of the measurement
error by

$\displaystyle\Delta{\theta_{i}}\equiv\left\langle\delta{\theta_{i}}\delta{\theta_{i}}\right\rangle^{1/2}.$
(47)

Here, importantly, the Fisher matrix and its inverse matrix have block
diagonal structures. We can divide the parameters
${\mbox{\boldmath${\theta}$}}_{n}$ into two subsets, on the basis of the
evenness and oddness of the time derivatives (effectively regarding $\varphi$
as the minus first derivative). The odd ones are
$\\{\varphi,{\dot{f}},{\dddot{f}},\cdots\\}$ and the even ones are
$\\{f,{\ddot{f}},\cdots\\}$. The two groups have no error correlation, owing
to the symmetric cancellations in the time integrals (45). For example, we
have

$\displaystyle F_{\varphi{\ddot{f}}}\propto\int_{-T/2}^{T/2}t^{3}dt=0.$ (48)

Here our midpoint expansion plays the critical role.

Our simplified waveform model (41) does not contain other parameters such as
the source direction and orientation. In particular, the phase function (42)
lacks the Doppler phase modulation induced by the annual rotation of the
detector around the Sun. Therefore, one might be skeptical to the validity of
our approach. These concerns were closely examined in Seto:2023skl  by
comparison of our analytical approach with the full Fisher matrix analysis
including other fitting parameters and the Doppler and amplitude modulations
Cutler:1997ta ; Takahashi:2002ky . It was shown that, for observational period
$T\gtrsim 2$yr, our simplified evaluations work quite efficiently for the
phase related parameters
${\mbox{\boldmath${\theta}$}}_{2}=\\{\varphi,f,{\dot{f}},{\ddot{f}}\\}$. This
is because, at $T\gtrsim 2$yr, we can distinguish the periodic Doppler
modulation from the perturbative frequency drift (42) truncated at a finite
order. As a result, the mutual correlation between the source geometric
parameters and those in Eq. (42) becomes weak, and our simple approach will be
valid even for $n=4$. Note that, For $T\gtrsim 2$yr, the cancellations similar
to Eq. (48) effectively work, even if we include the annual amplitude
modulation Seto:2023skl .

###  III.2 Analytical Results

In this subsection, for various orders $n$ of the phase expansion, we present
the concrete expressions for the estimation errors of the fitting parameters
${\mbox{\boldmath${\theta}$}}_{n}=\\{\varphi,f,{\dot{f}},\cdots,f^{(n)}\\}$.
As mentioned in the previous subsection, these expressions will be valid for
$T\gtrsim 2$yr.

If the apparent frequency evolution is well described by the second-order
expansion $\Phi_{2}(t)$ with the four parameters
${\mbox{\boldmath${\theta}$}}_{2}=\\{\varphi,f,{\dot{f}},{\ddot{f}}\\}$, we
can obtain the Fisher matrix results

$\displaystyle\Delta\varphi=\frac{3}{2\rho},~{}\Delta
f=\frac{5\sqrt{3}}{2\pi\rho T},~{}\Delta{\dot{f}}=\frac{6\sqrt{5}}{\pi\rho
T^{2}},~{}\Delta{\ddot{f}}=\frac{60\sqrt{7}}{\pi\rho T^{3}},$ (49)

as presented in Seto:2023skl . For the third-order expansion $\Phi_{3}(t)$
with the five fitting parameters
${\mbox{\boldmath${\theta}$}}_{3}=\\{\varphi,f,{\dot{f}},{\ddot{f}},{\dddot{f}}\\}$
(newly including $\dddot{f}$), we obtain

$\displaystyle\Delta\varphi$ $\displaystyle=$
$\displaystyle\frac{15}{8\rho},~{}\Delta f=\frac{5\sqrt{3}}{2\pi\rho
T},~{}\Delta{\dot{f}}=\frac{21\sqrt{5}}{\pi\rho
T^{2}},~{}\Delta{\ddot{f}}=\frac{60\sqrt{7}}{\pi\rho T^{3}}$
$\displaystyle\Delta{\dddot{f}}$ $\displaystyle=$
$\displaystyle\frac{2520}{\pi\rho T^{4}}.$ (50)

Because of the block diagonalization, the results $\\{\Delta
f,\Delta{\ddot{f}}\\}$ are the same as those in Eq. (49). In contrast, the
remaining ones $\\{\Delta\varphi,\Delta{\dot{f}}\\}$ become larger, reflecting
their correlation with the additional error $\Delta{\dddot{f}}$.

Next we move to the fourth-order function $\Phi_{4}(t)$ with the six fitting
parameters
${\mbox{\boldmath${\theta}$}}_{4}=\\{\varphi,f,{\dot{f}},{\ddot{f}},{\dddot{f}},{\ddddot{f}}\\}$.
We obtain

$\displaystyle\Delta\varphi$ $\displaystyle=$
$\displaystyle\frac{15}{8\rho},~{}\Delta f=\frac{35\sqrt{3}}{8\pi\rho
T},~{}\Delta{\dot{f}}=\frac{21\sqrt{5}}{\pi\rho
T^{2}},~{}\Delta{\ddot{f}}=\frac{270\sqrt{7}}{\pi\rho T^{3}}$
$\displaystyle\Delta{\dddot{f}}$ $\displaystyle=$
$\displaystyle\frac{2520}{\pi\rho
T^{4}},~{}\Delta{\ddddot{f}}=\frac{15120\sqrt{11}}{\pi\rho T^{5}}.$ (51)

Here the results $\\{\Delta\varphi,\Delta{\dot{f}},\Delta{\dddot{f}}\\}$ are
the same as those in Eq.(50), but others $\\{\Delta f,\Delta{\ddot{f}}\\}$
become larger. In this manner, under the midpoint expansion, we can append the
highest derivative coefficients $f^{(n+1)}$ to the fitting parameters, without
sacrificing the accuracy of the second highest one $f^{(n)}$. As we see in the
next section, this will be preferable at actual data analysis for nearly
monochromatic BBHs.

In Fig. 2, to show the dependence on the order of the phase expansion ($n=2,3$
and 4), we compare the numerical coefficients for the errors presented in Eqs.
(49)-(51) normalized by those presented in Eq. (51). We can see that the lower
order parameters such as $\varphi$ are less affected by the order of the
derivative expansion. This seems reasonable, given the differences of the
power indexes in the phase function (42).

Figure 2: Comparison of the numerical factors in Eqs. (49)-(51) normalized by
those in Eq. (51).

So far, we have applied the GW phase expansion (42) at the midpoint $t=0$ of
the observation period $[-T/2,T/2]$. For comparison, we evaluated the Fisher
matrix results, expanding the GW phase at the initial point $t=0$ of the
integration period $[0,T]$. For the six fitting parameters
${\mbox{\boldmath${\theta}$}}_{4}=\\{\varphi,f,{\dot{f}},{\ddot{f}},{\dddot{f}},{\ddddot{f}}\\}$,
we obtained

$\displaystyle\Delta\varphi$ $\displaystyle=$
$\displaystyle\frac{6}{\rho},~{}\Delta f=\frac{35\sqrt{3}}{\pi\rho
T},~{}\Delta{\dot{f}}=\frac{336\sqrt{5}}{\pi\rho
T^{2}},~{}\Delta{\ddot{f}}=\frac{2160\sqrt{7}}{\pi\rho T^{3}}$
$\displaystyle\Delta{\dddot{f}}$ $\displaystyle=$
$\displaystyle\frac{25200}{\pi\rho
T^{4}},~{}\Delta{\ddddot{f}}=\frac{15120\sqrt{11}}{\pi\rho T^{5}}.$ (52)

The magnitude $\Delta{\ddddot{f}}$ for the highest derivative is identical to
Eq. (51) for the midpoint expansion (as explained below), but other parameters
have much larger estimation errors. Note that, in the case of the initial
point expansion, the diagonal elements $F_{f^{(n)}f^{(n)}}$ of the Fisher
matrix (45) are $2^{2n+2}$ times larger than those for the midpoint expansion
(again regarding $\varphi\propto f^{(-1)}$). However, the diagonal elements of
the inverted matrix (see Eqs. (46) and (47)) become generally larger for the
initial point expansion, as shown in Eqs. (51)-(52). This is due to the off-
diagonal elements of the Fisher matrix (related to correlations of the
parameters), and the block diagonal structure of the midpoint expansion would
play an important role.

Eqs. (51) and (52) are obtained under the assumption that the phase function
$\Phi(t)$ is well approximated by the fifth order polynomial expansion of the
time variable $t$. Using the difference of the time origin $t=0$, we can write
down the third derivative coefficient of the initial point expansion

${\dddot{f}}\simeq{\dddot{f}}_{\rm m}-\frac{T}{2}{\ddddot{f}}_{\rm m}$ (53)

in terms of these evaluated under the midpoint expansion ${\dddot{f}}_{\rm m}$
and ${\ddddot{f}}_{\rm m}$ (temporarily adding the subscript m). Then, we have

$\Delta{\dddot{f}}\simeq\left(\Delta{\dddot{f}}_{\rm
m}^{2}+\frac{T^{2}}{4}\Delta{\ddddot{f}}_{\rm m}^{2}\right)^{1/2}$ (54)

due to the block diagonalization for the midpoint expansion. We can easily
confirm the expression $\Delta{\dddot{f}}$ in (52) from Eq. (51). We also have
${\ddddot{f}}={\ddddot{f}_{\rm m}}$ and thus
$\Delta{\ddddot{f}}=\Delta{\ddddot{f}_{\rm m}}$, as commented earlier.

###  III.3 Order of the Expansion

In any case, we should carefully adjust the order $n$ of the phase expansion
at the matched filtering analysis. An unnecessarily large $n$ results in
avoidable degradation of the important parameters such as $\ddot{f}$ (see Fig.
2). At the same time, they will be biased, if we take a too small $n$. In
reality, we should adjust the order $n$, also taking into account the observed
properties of the actual BBH samples detected by space detectors. Below, as a
working hypothesis, we stop the expansion at the lowest order $n$ satisfying
the resolution condition

$f^{(n)}\lesssim\Delta f^{(n)}.$ (55)

Here, as the reference values, we use the derivatives $f^{(n)}$ given in Eqs.
(2)-(5) for radiation reaction. At the next order expansion $n+1$, for our
fiducial BBHs, we will have $f^{(n+1)}\ll\Delta f^{(n+1)}$ from the scaling
relations $f^{(n)}\propto t_{c}^{-n}$ and $\Delta f^{(n)}\propto T^{-n+1/2}$,
together with the increasingly large numerical factors of the latter (see Eqs.
(49)-(51)).

In Fig. 3, assuming circular BBHs and using expressions in Secs. II A and III
B, we present the integration time $T$ satisfying the condition $\Delta
f^{(n)}=f^{(n)}$ for various chirp masses ${\cal M}$ at $f=5$mHz. For each
index $n$, we use the error $\Delta f^{(n)}$ obtained with the corresponding
expansion $\Phi_{n}(t)$ (e.g. $\Delta{\ddot{f}}$ from Eq. (49) not from Eq.
(51)). We can roughly understand the required order $n$ of the expansion as a
function of the integration time $T$.

In this figure, considering the applicability of our simplified phase model
(41), the results will be valid at $T\gtrsim 2$ yr. For the Fisher matrix, we
have also used the nearly monochromatic approximation for evaluating the time
integrals (45). Therefore, in Fig. 3, we show the time $T$ during which the GW
frequency $f$ changes by 1% and 5% (namely ${\dot{f}}T/f=0.01$ and $0.05$). In
Fig. 4, we preset a similar plot, now fixing the chirp mass at ${\cal
M}=30M_{\odot}$.

Figure 3: The integration time $T$ satisfying the condition $\Delta
f^{(n)}=f^{(n)}$ at $f=5$mHz. The results are given for circular BBHs which
have the signal-to-noise ratio $\rho=10$ in $T=4$yr. The GW frequency $f$
changes by 1% and 5% during the time $T$ shown by the black lines.  Figure 4:
The integration time $T$ satisfying the condition $\Delta f^{(n)}=f^{(n)}$ at
various GW frequency $f$. The results are given for circular BBHs with the
fixed chirp mass ${\cal M}=30M_{\odot}$ and $\rho=10$ in $T=4$yr. The results
at $f=5$mHz are identical to those for ${\cal M}=30M_{\odot}$ in Fig. 3.
Figure 5: The integration time $T$ satisfying the condition (57) for the phase
mismatch. The results are given for circular BBHs which have the GW frequency
$f=5$mHz and various chirp masses ${\cal M}$.  Figure 6: The integration time
$T$ satisfying the condition (57) for the phase mismatch. The results are
given for circular BBHs which have the fixed chirp mass ${\cal M}=30M_{\odot}$
various GW frequencies. The results at $f=5$mHz are identical to those for
${\cal M}=30M_{\odot}$ in Fig. 5.

We have assumed to determine the order $n$ of the Taylor expansion, following
the condition (55). This choice is equivalent to set
$f^{(n+1)}=f^{(n+2)}=\cdots=0$ for the higher derivative coefficients. The
associated phase mismatch is roughly given by

$\displaystyle\frac{f^{(n+1)}}{(n+2)!}t^{n+2},$ (56)

which takes the maximum magnitude at the initial and end points $t=\pm T/2$.
From the viewpoint of matched filtering analysis, the phase mismatch due to
the incomplete template should be much less than unity (see e.g.
Creighton:2011zz ).

Here we evaluate the integration time $T$ with which the phase mismatch
becomes unity as

$\displaystyle\frac{f^{(n+1)}}{(n+2)!}\left(\frac{T}{2}\right)^{n+2}=1.$ (57)

The results are presented in Figs. 5 and 6. Comparing Figs. 3 and 5, we can
see that, for a given $n$ (e.g. $n=3$), the integration time $T$ in Fig. 3 is
generally smaller than that in Fig. 5. Therefore, following the condition
(55), the associated phase mismatch will be much less than unity. We can
confirm similar results in Figs. 4 and 6. We should also note that the left-
hand side of Eq. (57) is proportional to $f^{(n+1)}T^{(n+2)}$, in a similar
way to the ratio $f^{(n+1)}/(\Delta f^{(n+1)})$.

##  IV observational prospects

We now discuss the prospects for measuring the frequency derivatives for our
fiducial BBH. For a circular orbit, from Eqs. (2)-(5) (also Table I), we have

$\displaystyle{\dot{f}}=6.2\times 10^{-13}{\rm Hz\,s^{-1}},$ (58)
$\displaystyle{\ddot{f}}=2.8\times 10^{-22}{\rm Hz\,s^{-2}},$ (59)
$\displaystyle{\dddot{f}}=2.2\times 10^{-31}{\rm Hz\,s^{-3}},$ (60)
$\displaystyle{\ddddot{f}}=2.4\times 10^{-40}{\rm Hz\,s^{-4}}.$ (61)

Below, also for mildly eccentric BBHs ($e\lesssim 0.1$), we use these results
as the reference values.

In this section, we mainly set $T=4$ yr and $\rho=10$, given the nominal
operation period of LISA. In section IV D, we discuss the case for an extended
case $T=10$ yr and $\rho=(10/4)^{1/2}\times 10\simeq 16$.

###  IV.1 Estimation Errors

Using expressions (49) for the second order phase expansion $\Phi_{2}(t)$, we
obtain

$\displaystyle\Delta{\dot{f}}=2.7\times 10^{-17}{\rm Hz\,s^{-1}}$ (62)
$\displaystyle\Delta{\ddot{f}}=2.5\times 10^{-24}{\rm Hz\,s^{-2}}$ (63)

for $\rho=10$ and $T=4$yr. Given the relation $\Delta{\ddot{f}}\ll{\ddot{f}}$
for the highest derivative $\ddot{f}$ (presented in Eq. (59)), we examine the
next-order function $\Phi_{3}(t)$. From expression (50), we obtain

$\displaystyle\Delta{\dot{f}}=9.4\times 10^{-17}{\rm Hz\,s^{-1}}$ (64)
$\displaystyle\Delta{\ddot{f}}=2.5\times 10^{-24}{\rm Hz\,s^{-2}}$ (65)
$\displaystyle\Delta{\dddot{f}}=3.2\times 10^{-31}{\rm Hz\,s^{-3}}.$ (66)

Now we have a relation ${\dddot{f}}\lesssim\Delta{\dddot{f}}$ for the highest
derivative. As shown in Eqs. (63) and (65), due to the block diaginalization,
the error $\Delta{\ddot{f}}$ is invariant. In contrast, the error
$\Delta{\dot{f}}$ in Eq. (64) becomes 7/2 times larger than that in Eq. (62).
We still have $\Delta{\dot{f}}\ll{\dot{f}}$ and this degradation of the factor
7/2 might not be a practial problem from astrophysical point of views. This is
because we usually do not need to very accurately (sub-percent level) measure
the rate ${\dot{f}}$ (or the redshifted chirp mass) alone. When we extract
some astrophysical information by combining the rate $\dot{f}$ with other
fitting parameters (e.g. $\ddot{f}$ for the combination $C_{2}$ as discussed
in section II B), the error $\Delta{\dot{f}}$ will have a minor contribution
to the total error.

Using Eq. (51) for the fourth-order expansion $\Phi_{4}(t)$, we have the
estimation errors

$\displaystyle\Delta{\dot{f}}=9.4\times 10^{-17}{\rm Hz\,s^{-1}}$ (67)
$\displaystyle\Delta{\ddot{f}}=1.1\times 10^{-23}{\rm Hz\,s^{-2}}$ (68)
$\displaystyle\Delta{\dddot{f}}=3.2\times 10^{-31}{\rm Hz\,s^{-3}}$ (69)
$\displaystyle\Delta{\ddddot{f}}=5.0\times 10^{-38}{\rm Hz\,s^{-4}}.$ (70)

We now have ${\ddddot{f}}\ll\Delta{\ddddot{f}}$ for the highest derivative,
and this expansion will be excessive in view of the criteria (55), hampering
the resolution $\Delta{\ddot{f}}$ (increased by a factor of 9/2 from Eq.
(65)). The reasonable choice will be $n=3$, and we below use Eqs. (64)-(66)
for the estimation errors.

###  IV.2 Eccentricity Estimation

We next discuss the extraction of the astrophysical information from the
observed products $C_{n}$.

Given the hierarchy of the relative errors $\Delta f^{(n)}/f^{(n)}$, the
measurement error $\Delta C_{n}$ is estimated to be

$\displaystyle\Delta C_{n}\sim\frac{\Delta
C_{n}}{C_{n}}\sim\frac{\Delta{f^{(n)}}}{f^{(n)}}$ (71)

up to moderate eccentricity with $C_{n}\sim 1$. For our fiducial BBH, from
Eqs. (59) and (65), we have $\Delta C_{2}\sim 9\times 10^{-3}$. With the
relation $\Delta{\dddot{f}}\sim{\dddot{f}}$ for the next order one $C_{3}$, we
cannot measure $C_{3}$ in real terms. In addition, even for the reference
system in Sec. II C, the shift of $C_{2}$ due to the tertiary is comparable to
its error $\Delta C_{2}$. Therefore, below in this subsection, we concentrate
on the eccentricity measurement from the product $C_{2}$.

If the observed value $C_{2}$ is consistent with unity (namely
$|C_{2}-1|\lesssim\Delta C_{2}$), we can apply Eq. (20) and set an upper limit
to the eccentricity as

$\displaystyle e\lesssim\left\\{\frac{{\color[rgb]{0,0,0}792}}{2983}\Delta
C_{2}\right\\}^{1/2}.$ (72)

We have $e\lesssim 0.0{\color[rgb]{0,0,0}5}$ for our fiducial BBH.

In contrast, if the product $C_{2}$ satisfies the condition
$C_{2}-1{\color[rgb]{0,0,0}\gg}\Delta C_{2}$, we can solve Eq. (19) for the
eccentricity $e$ with the estimation error

$\displaystyle\Delta e\sim\frac{\Delta C_{2}}{|\partial_{e}{K}_{2}(e)|}.$ (73)

As an example, for the true value $e=0.1$, we have
$C_{2}={K}_{2}(0.1)={\color[rgb]{0,0,0}0.96}$, and the expected estimation
error becomes $\Delta e\sim 1.3\times 10^{-2}$.

###  IV.3 Extension to 10yr

In this subsection, we discuss the results for an extended observation of
$T=10$ yr. Compared with $T=4$ yr, the signal-to-noise $\rho$ increases by a
factor of $2.5^{1/2}\sim 1.6$ and the estimation errors $\Delta f^{(n)}$
formally shrink by $0.4^{n+3/2}$ (0.1 for $n=2$, 0.016 for $n=3$ and 0.006 for
$n=4$). However, for a longer $T$, we need to increase the order $n$ of the
phase expansion, following the criteria (55). Therefore, these scaling factors
cannot be simply applicable.

Table 2: The resolutions for the frequency derivatives at $T=4$ and 10 yr (respectively from Eqs. (64)-(66) and Eqs. (77)-(80)).  $T$ | 4yr | 10yr | ratio
---|---|---|---
$\Delta{\dot{f}}$ | $9.4\times 10^{-17}{\rm Hz\,s^{-1}}$ | $9.5\times 10^{-18}{\rm Hz\,s^{-1}}$ | 0.10
$\Delta{\ddot{f}}$ | $2.5\times 10^{-24}{\rm Hz\,s^{-2}}$ | $4.6\times 10^{-25}{\rm Hz\,s^{-2}}$ | 0.18
$\Delta{\dddot{f}}$ | $3.2\times 10^{-31}{\rm Hz\,s^{-3}}$ | $5.2\times 10^{-33}{\rm Hz\,s^{-3}}$ | 0.016
$\Delta{\ddddot{f}}$ |  | $3.3\times 10^{-40}{\rm Hz\,s^{-4}}$ |

For the third-order phase function $\Phi_{3}(t)$, from Eq. (50), we have the
estimation errors

$\displaystyle\Delta{\dot{f}}=9.5\times 10^{-18}{\rm Hz\,s^{-1}}$ (74)
$\displaystyle\Delta{\ddot{f}}=1.0\times 10^{-25}{\rm Hz\,s^{-2}}$ (75)
$\displaystyle\Delta{\dddot{f}}=5.2\times 10^{-33}{\rm Hz\,s^{-3}}$ (76)

with $\Delta{\dddot{f}}\ll{\dddot{f}}$ for the highest derivative. Similarly,
from Eq. (51) for the function $\Phi_{4}(t)$ we obtain

$\displaystyle\Delta{\dot{f}}=9.5\times 10^{-18}{\rm Hz\,s^{-1}}$ (77)
$\displaystyle\Delta{\ddot{f}}=4.6\times 10^{-25}{\rm Hz\,s^{-2}}$ (78)
$\displaystyle\Delta{\dddot{f}}=5.2\times 10^{-33}{\rm Hz\,s^{-3}}$ (79)
$\displaystyle\Delta{\ddddot{f}}=3.3\times 10^{-40}{\rm Hz\,s^{-4}}$ (80)

now with the desired relation ${\ddddot{f}}\lesssim\Delta{\ddddot{f}}$. Below,
we use the estimation errors in Eq. (77)-(80). In Table II, we compare these
results with those estimated for $T=4$yr. Because of the block diagonal
structure, we can significantly improve the resolution $\Delta{\dddot{f}}$
with the full scaling relation $0.4^{9/2}$ mentioned before.

The error $\Delta{\ddot{f}}$ for $T=10$yr is 5 times smaller than that for
$T=4$yr. We then have

$\displaystyle\Delta C_{2}\sim\frac{\Delta{\ddot{f}}}{\ddot{f}}\sim 1.6\times
10^{-3},$ (81)

and the PN effects are likely to be still unresolvable. However, for a larger
chirp mass ${\cal M}$ and a higher frequency $f$, the PN effects should be
taken into account.

From Eqs. (60) and (79), we can evaluate the estimation error for $C_{3}$

$\displaystyle\Delta C_{3}\sim\frac{\Delta{\dddot{f}}}{\dddot{f}}\sim
2.4\times 10^{-2}.$ (82)

As expected, we have $\Delta C_{3}\gg\Delta C_{2}$.

###  IV.4 Tertiary Signatures

Now, based on the results for $T=10$yr, we discuss the search for the tertiary
effect. For a nearly circular BBH perturbed by the reference tertiary model in
Sec. II C, we have

$C_{2}-1\simeq\frac{{\ddot{f}}_{s}}{{\ddot{f}}_{r}}\simeq-6.0\times
10^{-3}\sin\psi_{s},$ (83)

which becomes incompatible with an isolated BBH for $\sin\psi_{s}<0$. The
excess can be measured with the resolution $\Delta C_{2}$ given in Eq. (81),
and we can detect it at $3.6\sigma$ for the optimal phase  $\sin\psi_{s}=-1$
(corresponding to ${\ddot{f}}_{s}\sim 2\times 10^{-24}{\rm Hz\,s^{-2}}$). For
$T=4$yr, the same excess will be totally masked by the measurement noise with
$|C_{2}-1|\ll\Delta C_{2}$.

Next, for an eccentric binary, we can examine the inconsistency between
$C_{2}$ and $C_{3}$ induced by the tertiary perturbation. In fact, given the
condition $\Delta C_{3}\gg\Delta C_{2}$, it will be more straightforward to
compare $C_{3}$ and the corresponding value $K_{3}[K_{2}^{-1}(C_{2})]$. The
uncertainty $\Delta C_{3}$ dominates the statistical error for the gap between
them. From Eqs. (36) and (38), the gap is approximately estimated to be

$C_{3}-K_{3}[K_{2}^{-1}(C_{2})]\simeq\frac{{\dddot{f}}_{s}}{{\dddot{f}}_{r}}\simeq-2.6\times
10^{-2}\cos\psi_{s}.$ (84)

Given the uncertainty (82), the third derivative ${\dddot{f}}_{s}$ (due to the
reference tertiary perturbation) can be detected only at $1\sigma$ level, even
for the optimal orbital phase $|\cos\psi_{s}|=1$.

##  V discussions

Here are some additional remarks regarding our current investigation.

###  V.1 Other Harmonic Modes

So far, we have mainly discussed observation for GW around a single frequency
$f$ (twice the orbital frequency). An eccentric binary emits GWs at integer
times the orbital frequency. Up to a moderate eccentricity $e\lesssim 0.3$,
the third harmonic mode at the frequency $3f/2$ will be the secondary target,
and its amplitude is proportional to $e$ with its SNR $\sim\rho e$ ($\rho$:
the SNR of the primary mode). As argued in Seto:2016wom ; Seto:2001pg , by
closely analyzing the signatures induced by the 1PN apsidal precession on the
third harmonics, we might determine the total mass $M_{T}$ of the BBH.

For our fiducial BBH ($\rho\sim 10$ and $f\sim 5$mHz), the third harmonics is
detectable for $e\sim 0.3$, and we can estimate the orbital eccentricity $e$
from the amplitude of the mode. Unfortunately, the resultant error $\Delta e$
will be much larger than that from $C_{2}$ as discussed in Eq. (73). However,
in the lower frequency regime e.g. $f\sim 1$ mHz, we cannot measure the second
derivative $\ddot{f}$ (see Fig. 4), and the third harmonics could be a
relatively better target.

###  V.2 Long-term Noise Variation

For the block diagonalization associated with the midpoint expansion, we have
assumed that the effective noise spectrum $S_{n}(f)$ does not depend on time
(more preciously, during the course of the matched filtering analysis). In
practice, detectors may degrade over time, and we might also operate a network
consisting of multiple detectors (e.g., LISA+Taiji Cai:2023ywp ) functioning
at different time intervals. In such scenarios, we cannot straightforwardly
apply the odd-even cancellation (48) based solely on the stationarity of the
noise $S_{n}(f)$. However, we can still adjust the time origin $t=0$ for the
perturbative phase expansion and achieve similar (but more restricted)
features to the structure of the Fisher matrix.

##  VI summary

LISA has a potential to detect stellar mass BBHs, which are similar to those
already observed by ground-based detectors. According to a recent statistical
study Seto:2022xmh , the majority of the detected BBHs will emit nearly
monochromatic GWs, exhibiting slight frequency variations within the nominal
operational period of $T\sim 4$ years. To characterize such signals, we can
efficiently employ a Taylor expansion for the evolution of the GW frequency
$f$. The associated derivative coefficients, such as ${\dot{f}}$ and
${\ddot{f}}$, may harbor intriguing astrophysical information, such as orbital
eccentricity and tertiary perturbation effects. To facilitate this analysis,
we introduced non-dimensional measures $C_{n}$ ($n=2,3,\cdots$) with $C_{n}=1$
for an isolated circular binary.

We provided simple analytical expressions for the parameter estimation errors
for the derivatives $f^{(n)}$, applying the Fisher matrix formulation to a
simplified phase model $\Phi_{n}(t)$. Then, we concretely discuss the
measurability of the aforementioned astrophysical information for our fiducial
BBH (chirp mass ${\cal M}\sim 30M_{\odot}$, the GW frequency $f\sim 5$mHz and
the signal-to-noise ratio $\rho=10$ for $T=4$yr). We found that LISA can make
a resolution $\Delta C_{2}\sim\Delta{\ddot{f}}/{\ddot{f}}\sim 10^{-2}$,
corresponding to the eccentricity error $\Delta e\sim 1.3\times 10^{-2}$ for
$e=0.1$. For an extended operation period of $T\sim 10$yr, we might also make
the next order combination $C_{3}$ and detect a tertiary perturbation as small
as ${\ddot{f}}_{s}\sim 2\times 10^{-24}{\rm Hz\,s^{-2}}$ for a nearly circular
BBH.

###### Acknowledgements.

Equation (15) plays an important role in this paper. The author is deeply
grateful to Tomoya Suzuguchi for pointing out errors in its previous version.

## References

  * (1) R. Abbott et al. [KAGRA, VIRGO and LIGO Scientific], Phys. Rev. X 13, no.4, 041039 (2023) doi:10.1103/PhysRevX.13.041039 [arXiv:2111.03606 [gr-qc]].
  * (2) P. Amaro-Seoane et al. [LISA], [arXiv:1702.00786 [astro-ph.IM]].
  * (3) W. H. Ruan, Z. K. Guo, R. G. Cai and Y. Z. Zhang, Int. J. Mod. Phys. A 35, no.17, 2050075 (2020) doi:10.1142/S0217751X2050075X [arXiv:1807.09495 [gr-qc]].
  * (4) J. Luo et al. [TianQin], Class. Quant. Grav. 33, no.3, 035010 (2016) doi:10.1088/0264-9381/33/3/035010 [arXiv:1512.02076 [astro-ph.IM]].
  * (5) P. Amaro-Seoane et al. [LISA], Living Rev. Rel. 26 (2023) no.1, 2 doi:10.1007/s41114-022-00041-y [arXiv:2203.06016 [gr-qc]].
  * (6) A. Sesana, Phys. Rev. Lett. 116, no.23, 231102 (2016) doi:10.1103/PhysRevLett.116.231102 [arXiv:1602.06951 [gr-qc]].
  * (7) N. Seto and K. Kyutoku, Mon. Not. Roy. Astron. Soc. 514, no.4, 4669-4675 (2022) doi:10.1093/mnras/stac1561 [arXiv:2201.02766 [astro-ph.HE]].
  * (8) D. Gerosa, S. Ma, K. W. K. Wong, E. Berti, R. O’Shaughnessy, Y. Chen and K. Belczynski, Phys. Rev. D 99, no.10, 103004 (2019) doi:10.1103/PhysRevD.99.103004 [arXiv:1902.00021 [astro-ph.HE]].
  * (9) A. Klein, G. Pratten, R. Buscicchio, P. Schmidt, C. J. Moore, E. Finch, A. Bonino, L. M. Thomas, N. Williams and D. Gerosa, et al. [arXiv:2204.03423 [astro-ph.HE]].
  * (10) T. E. Strohmayer 2021, ApJL, 912, L8. doi:10.3847/2041-8213/abf3cc
  * (11) J. Munday, T. R. Marsh, M. Hollands, I. Pelisoli, D. Steeghs, P. Hakala, E. Breedt, A. Brown, V. S. Dhillon and M. J. Dyer, et al. Mon. Not. Roy. Astron. Soc. 518 (2022) no.4, 5123-5139 doi:10.1093/mnras/stac3385 [arXiv:2211.09834 [astro-ph.SR]].
  * (12) N. Seto, Phys. Rev. D 108, no.2, 023019 (2023) doi:10.1103/PhysRevD.108.023019 [arXiv:2307.04453 [gr-qc]].
  * (13) R. Ebadi, V. Strokov, E. H. Tanin, E. Berti and R. L. Walsworth, [arXiv:2405.13109 [gr-qc]].
  * (14) P. C. Peters, Phys. Rev. 136, B1224-B1232 (1964) doi:10.1103/PhysRev.136.B1224
  * (15) N. Seto, Mon. Not. Roy. Astron. Soc. 460, no.1, L1-L4 (2016) doi:10.1093/mnrasl/slw060 [arXiv:1602.04715 [astro-ph.HE]].
  * (16) A. Nishizawa, A. Sesana, E. Berti and A. Klein, Mon. Not. Roy. Astron. Soc. 465, no.4, 4375-4380 (2017) doi:10.1093/mnras/stw2993 [arXiv:1606.09295 [astro-ph.HE]].
  * (17) N. Seto, Astrophys. J. Lett. 677, L55-L58 (2008) doi:10.1086/587785 [arXiv:0802.3411 [astro-ph]].
  * (18) N. Seto and K. Kyutoku, Phys. Rev. Lett. 118, no.15, 151101 (2017) doi:10.1103/PhysRevLett.118.151101 [arXiv:1611.05875 [astro-ph.HE]].
  * (19) N. Tamanini and C. Danielski, Nature Astron. 3, no.9, 858-866 (2019) doi:10.1038/s41550-019-0807-y [arXiv:1812.04330 [astro-ph.EP]].
  * (20) K. W. K. Wong, V. Baibhav and E. Berti, Mon. Not. Roy. Astron. Soc. 488, no.4, 5665-5670 (2019) doi:10.1093/mnras/stz2077 [arXiv:1902.01402 [astro-ph.HE]].
  * (21) Y. Kang, C. Liu and L. Shao, Astron. J. 162, no.6, 247 (2021) doi:10.3847/1538-3881/ac23d8 [arXiv:2108.01357 [astro-ph.EP]].
  * (22) T. Robson, N. J. Cornish, N. Tamanini and S. Toonen, Phys. Rev. D 98, no.6, 064012 (2018) doi:10.1103/PhysRevD.98.064012 [arXiv:1806.00500 [gr-qc]].
  * Xuan, Peng, & Chen (2021) Z. Xuan, P. Peng, X. Chen, 2021, MNRAS, 502, 4199. doi:10.1093/mnras/stab331
  * (24) F. Zhang, X. Chen, L. Shao and K. Inayoshi, Astrophys. J. 923, no.2, 139 (2021) doi:10.3847/1538-4357/ac2c07 [arXiv:2109.14842 [astro-ph.HE]].
  * (25) K. i. Maeda, P. Gupta and H. Okawa, Phys. Rev. D 107, no.12, 124039 (2023) doi:10.1103/PhysRevD.107.124039 [arXiv:2303.16553 [gr-qc]].
  * (26) N. Tamanini, A. Klein, C. Bonvin, E. Barausse and C. Caprini, Phys. Rev. D 101, no.6, 063002 (2020) doi:10.1103/PhysRevD.101.063002 [arXiv:1907.02018 [astro-ph.IM]].
  * (27) T. Takahashi, H. Omiya and T. Tanaka, Phys. Rev. D 107, no.10, 103020 (2023) doi:10.1103/PhysRevD.107.103020 [arXiv:2301.13213 [gr-qc]].
  * (28) L. Blanchet, Living Rev. Rel. 17, 2 (2014) doi:10.12942/lrr-2014-2 [arXiv:1310.1528 [gr-qc]].
  * (29) C. Cutler and E. E. Flanagan, Phys. Rev. D 49, 2658-2697 (1994) doi:10.1103/PhysRevD.49.2658 [arXiv:gr-qc/9402014 [gr-qc]].
  * (30) M. Punturo, M. Abernathy, F. Acernese, B. Allen, N. Andersson, K. Arun, F. Barone, B. Barr, M. Barsuglia and M. Beker, et al. Class. Quant. Grav. 27, 194002 (2010) doi:10.1088/0264-9381/27/19/194002
  * (31) S. Tanay, M. Haney and A. Gopakumar, Phys. Rev. D 93, no.6, 064031 (2016) doi:10.1103/PhysRevD.93.064031 [arXiv:1602.03081 [gr-qc]].
  * (32) N. Seto, Mon. Not. Roy. Astron. Soc. 524, no.4, 5442-5445 (2023) doi:10.1093/mnras/stad2187 [arXiv:2307.09739 [astro-ph.HE]].
  * (33) C. G. Bassa, G. H. Janssen, B. W. Stappers, T. M. Tauris, T. Wevers, P. G. Jonker, L. Lentati, J. P. W. Verbiest, G. Desvignes and E. Graikou, et al. Mon. Not. Roy. Astron. Soc. 460, no.2, 2207-2222 (2016) doi:10.1093/mnras/stw1134 [arXiv:1604.00129 [astro-ph.HE]].
  * (34) D. L. Kaplan, T. Kupfer, D. J. Nice, A. Irrgang, U. Heber, Z. Arzoumanian, E. Beklen, K. Crowter, M. E. DeCesar and P. B. Demorest, et al. Astrophys. J. 826, no.1, 86 (2016) doi:10.3847/0004-637X/826/1/86 [arXiv:1604.00131 [astro-ph.HE]].
  * (35) J. Samsing, K. Hendriks, L. Zwick, D. J. D’Orazio and B. Liu, [arXiv:2403.05625 [astro-ph.HE]].
  * (36) N. Seto, Mon. Not. Roy. Astron. Soc. 333, 469 (2002) doi:10.1046/j.1365-8711.2002.05432.x [arXiv:astro-ph/0202364 [astro-ph]].
  * (37) T. Robson, N. J. Cornish and C. Liu, Class. Quant. Grav. 36, no.10, 105011 (2019) doi:10.1088/1361-6382/ab1101 [arXiv:1803.01944 [astro-ph.HE]].
  * (38) R. Takahashi and N. Seto, Astrophys. J. 575, 1030-1036 (2002) doi:10.1086/341483 [arXiv:astro-ph/0204487 [astro-ph]].
  * (39) C. Cutler, Phys. Rev. D 57, 7089-7102 (1998) doi:10.1103/PhysRevD.57.7089 [arXiv:gr-qc/9703068 [gr-qc]].
  * (40) J. D. E. Creighton and W. G. Anderson, Gravitationalwave physics and astronomy: An introduction to theory, experiment and data analysis (2011).
  * (41) N. Seto, Phys. Rev. Lett. 87, 251101 (2001) [erratum: Phys. Rev. Lett. 101, 209901 (2008)] doi:10.1103/PhysRevLett.87.251101 [arXiv:astro-ph/0111107 [astro-ph]].
  * (42) R. G. Cai, Z. K. Guo, B. Hu, et al., [arXiv:2305.04551 [gr-qc]].