††thanks: email<EMAIL_ADDRESS>
# Probing Formation of Double Neutron Star Binaries around 1mHz with LISA

Lucy O. McNeill  Naoki Seto  Department of Physics, Kyoto University, Kyoto
606-8502, Japan

###### Abstract

We propose a novel method to examine whether Galactic double neutron star
binaries are formed in the LISA band. In our method, we assign an effective
time fraction $\tau$ to each double neutron star binary detected by LISA. This
fraction is given as a function of the observed orbital period and
eccentricity and should be uniformly distributed in the absence of in-band
binary formation. Applying statistical techniques such as the Kolmogorov-
Smirnov test to the actual list of $\tau$, we can inspect the signature of the
in-band binary formation. We discuss the prospects of this method, paying
close attention to the available sample number of Galactic double neutron star
binaries around 1mHz.

††preprint: APS/123-QED

##  I Introduction

Double neutron star binaries (hereafter DNSBs) are bountiful astrophysical
targets. They have been detected as radio pulsars [1, 2], and the current
sample has orbital periods $P$ from 1.9 hours [3] to six weeks [4]. This
sample contains around 20 DNSBs, and searches might be incomplete at the faint
end of the luminosity function (also limited by the beaming fraction). Shorter
period ($P<$ 1 hour) DNSBs ought to exist in the Galaxy as well. However, due
to Doppler smearing and shorter merger timescales, their detection would be
more difficult than the longer period ones [5].

The Laser Interferometer Space Antenna (LISA) is planned in the 2030’s [6] and
is sensitive to gravitational waves (GWs) around 0.1-100mHz. It will detect
all Galactic DNSBs in the frequency range $f\gtrsim$1.5mHz (corresponding to
the orbital period $P=2/f\lesssim 20$ min), unlike the longer-period radio
sample. Observationally motivated estimates [7] suggest that at least dozens
of DNSBs exist in the Galaxy at $f\gtrsim$1.5mHz. Numerical galactic modelling
[8, 9] predicts that LISA will detect altogether from a few to upto hundreds
of DNSBs.

DNSB formation depends on the complex interplay between many astrophysical
processes. In the general picture, the binary must survive two supernova
explosions. Preceding these supernovae, various mechanisms have been theorized
with respect to mass loss/exchange in the binary [10, 11, 12, 13] after
hydrogen burning. These processes all play a role in determining the
separation at formation, if the binary survives.

However, the related efficiencies and rates in populations are not well
established (see [14] for a review). In particular, it’s unclear whether there
is a critical minimum orbital period (or maximum orbital frequency) for
isolated DNSB formation.

Dynamical encounters in star clusters is an alternative pathway for short
period DNSB formation, though their contribution to the LISA sample is
estimated to be small [15, 16]. Also, for the dynamical scenario it will be
difficult to solidly estimate the distribution function for the orbital
periods of the generated DNSBs.

In this paper, we are interested in the possibility of DNSB formation
specifically at $f\gtrsim 1$mHz. We hereafter call this channel as “in-band”
(mHz) DNSB formation, or simply “injections”. Considering the aforementioned
theoretical uncertainties, it will be fruitful to observationally examine the
in-band formation in a model independent manner.

We thus develop a statistical method to examine the in-band formation with
LISA (see also [17] for formation between the LISA band and the lower
frequencies already probed by the radio sample). Recently, several studies
have proposed to statistically deal with multiple LISA sources in the Galaxy.
Among others, a large number ($\sim 10^{4}$) of white dwarf binaries (WDBs)
will be a powerful data set for various astronomical analyses (e.g., probing
the Galactic structure [18, 19]). In this context, one of the authors
suggested to measure the fluxes of the Galactic WDBs in frequency space [20].
He pointed out that the measurement will enable us to follow the collective
evolution of the WDBs, resulting in mergers or stable mass transfers.

One might imagine that we can get some information about the in-band formation
by similarly measuring the DNSB flux with LISA at various frequencies.
Unfortunately, LISA will detect much fewer DNSBs than WDBs, and the small
number statistics will severely limit the flux approach for DNSBs. On the
other hand, unlike WDBs, DNSBs can be well regarded as point particle systems
in the LISA band, and their long-term orbital evolution from GW emission can
be predicted quite accurately [1]. Considering these pros and cons of DNSBs,
we newly introduce the effective time fraction $\tau$, corresponding to the
fraction of total time that each binary has spent in the mHz band. Without in-
band formation, the fraction $\tau$ should be uniformly distributed. By
analysing the observed list of $\tau$, we can examine potential in-band
formation through its deviation from a uniform distribution.

The basic assumption of our study is that the mHz DNSB population is in the
“steady state” [21]. This is a reasonable assumption, since a DNSB passes
through the mHz band in the timescale of Myr but the Galactic DNSB merger rate
will change in the Hubble timescale [22, 23].

This paper is organized as follows. In section II, we roughly estimate how
many Galactic DNSBs are likely to be detected by LISA. In section III we study
the long-term orbital evolution of DNSBs and define the effective time
fraction $\tau$, as a function of the GW frequency and the orbital
eccentricity. Then, in section IV, using statistical tools such as the
Kolmogorov-Smirnov test, we discuss how well we can probe in-band formation
with LISA. In section V, we mention potential extensions of this study. We
summarize our findings in section VI.

##  II Galactic double neutron star binaries

###  II.1 Expected number in the LISA band

We first estimate the merger rate $R_{\mathrm{MW}}$ of DNSBs in our Galaxy
(Milky Way). As a basic observational input, we use the comoving merger rate
${\mathcal{R}}=660^{+1040}_{-530}{\rm Gpc^{-3}yr^{-1}}$ from a recent report
by the LVK collaboration [24] (Multi source model).

To relate the two rates $\cal R$ and $R_{\rm MW}$, we apply the traditional
argument based on the effective number density of Milky Way equivalent
galaxies [15, 25], and put $R_{\rm MW}=L_{B,\mathrm{MW}}{\cal
R}/\mathcal{L}_{B}$. Here $\mathcal{L}_{B}$ is the B-band luminosity per
comoving volume and $L_{B,\mathrm{MW}}$ is the B-band luminosity of our
Galaxy. Using their typical values, we obtain

$\displaystyle R_{\rm MW}=$ $\displaystyle 6.0\times
10^{-5}\mathrm{yr}^{-1}\left(\frac{\mathcal{R}}{660\mathrm{Gpc}^{-3}\mathrm{yr}^{-1}}\right)$
$\displaystyle\times\left(\frac{L_{B,\mathrm{MW}}}{9\times
10^{9}L_{\odot}}\right)\left(\frac{\mathcal{L}_{B}}{10^{17}L_{\odot}\mathrm{Gpc}^{-3}}\right)^{-1}.$
(1)

We should notice that the Galactic merger rate $R_{\mathrm{MW}}$ still has
large uncertainties (at least a factor of three).

Next, we roughly estimate the total number of Galactic DNSBs in the LISA band.
In this paper, we use the notation $f$ specifically for the second harmonic
frequency (given by $f=2f_{\mathrm{orb}}$ with the orbital frequency
$f_{\mathrm{orb}}$). Due to radiation reaction, the GW frequency $f$ evolves
as [26]

$\frac{\mathrm{d}f}{\mathrm{d}t}=\frac{96\pi^{8/3}G^{5/3}f^{11/3}\mathcal{M}^{5/3}}{5c^{5}(1-e^{2})^{7/2}}\left(1+\frac{73}{24}e^{2}+\frac{37}{96}e^{4}\right)$
(2)

with the orbital eccentricity $e$ and the chirp mass $\cal M$ of the binary
[26]. At the same time, eccentricity evolves according to

$\frac{\mathrm{d}e}{\mathrm{d}t}=-\frac{304e\pi^{8/3}G^{5/3}f^{8/3}\mathcal{M}^{5/3}}{15c^{5}(1-e^{2})^{5/2}}\left(1+\frac{121}{304}e^{2}\right).$
(3)

These equations are a good approximation for DNSBs in the LISA band since the
relativistic corrections are small (except for $1-e\ll 1$).

Simply assuming (i) the steady state condition for the Galactic DNSB
population at $f\gtrsim 1$mHz, and (ii) no binary formation there, we have the
frequency distribution

$\frac{\mathrm{d}N}{\mathrm{d}f}=R_{\rm
MW}\left(\frac{\mathrm{d}f}{\mathrm{d}t}\right)^{-1}\propto R_{\rm
MW}{f^{-11/3}}{\mathcal{M}^{-5/3}}.$ (4)

Here we ignored the eccentricity dependence of $\dot{f}$. After the frequency
integral, we obtain the cumulative number

$\begin{split}N(>f)&=30\left(\frac{\mathcal{M}}{1.2M_{\odot}}\right)^{-5/3}\left(\frac{f}{1.5\mathrm{mHz}}\right)^{-8/3}\\\
&\times\left(\frac{R_{\rm MW}}{6.0\times
10^{-5}\mathrm{yr}^{-1}}\right).\end{split}$ (5)

Note that the chirp mass distribution of known DNSBs is narrow and centred
around $\mathcal{M}=1.2M_{\odot}$ [27].

In fact, we will later relax the assumption (ii). But the above result will
serve as a rough guide, except for extreme model settings.

###  II.2 Gravitational wave observations

####  II.2.1 Identification of binary neutron stars

Next, we consider a DNSB at distance $D$, with chirp mass $\mathcal{M}$,
eccentricity $e$ and gravitational wave frequency $f$. When the binary is
approximated as monochromatic ($f=$ constant), the dimensionless gravitational
wave strain amplitude from the binary in the second orbital harmonic is

$h_{2}=\frac{8G^{5/3}\mathcal{M}^{5/3}\pi^{2/3}f^{2/3}}{5^{1/2}Dc^{4}}\left[1-\frac{5}{2}e^{2}+\frac{35}{24}e^{4}+O(e^{6})\right]$
(6)

[28, 29, 30]. This expression is obtained using the strain amplitude in the
$n$th orbital harmonic $h_{n}\propto g(n,e)^{1/2}/n$, with $g(n,e)$ given by
Equation (20) in [31]. For binaries upto $e\leq 0.2$, the correction for the
eccentricity is less than $\sim$10 percent.

In a gravitational wave detector with a sensitivity curve $S_{n}(f)$, the
angle averaged signal-to-noise ratio $\bar{\rho}_{2}$ over an observing time
$T$ is given by

$\bar{\rho}_{2}=\frac{h_{2}}{S_{n}(f)^{1/2}}\sqrt{T}.$ (7)

We put the noise curve of LISA by $S_{n}(f)=S_{d}(f)+S_{c}(f)$ with the
detector noise $S_{d}(f)$ and astrophysical foreground confusion noise
$S_{c}(f)$ [28], where the latter is a function of $T$ [32].

Now we will calculate some estimates related to gravitational wave detection
of Galactic DNSBs. Taking a circular binary with the chirp mass
$1.2M_{\odot}$, conservatively located at 20kpc, Equations (6) and (7) are
used to obtain $\bar{\rho}_{2}=7.8$ and 452 for $f=$ 1.5 mHz and 50mHz
respectively for an observation time $T=$ 4 years. If the observation time is
increased to $T=$ 10 years, these increase to $\bar{\rho}_{2}=$ 17 and 715,
respectively.

The gravitational wave frequency derivative, $\dot{f}$, characterises small
frequency drifts due to GW emission in Equation (2). Even though we used the
approximation that $\dot{f}=0$ in Equations (6) and (7), measuring this small
quantity is of paramount importance. Specifically, $\dot{f}$ can be used in
conjunction with Equation (6) to constrain the chirp mass
$\mathcal{M}\propto{\dot{f}}^{3/5}$ (also accounting for the measured
eccentricity discussed shortly).

In terms of the signal–to–noise ratio $\bar{\rho}_{2}$ and observation time
$T$ (longer than 2yr), the frequency derivative can be measured with
resolution [33]:

$\Delta\dot{f}\simeq 0.43\left(\frac{\bar{\rho}_{2}}{10}\right)^{-1}T^{-2}.$
(8)

Using this and Equation (2), we roughly estimate the fractional frequency
resolution ${\dot{f}/}{\Delta\dot{f}}$ for low eccentricity DNSBs located at
20kpc. If we consider binaries with $f=$ 1.5 and 2mHz, then over a 4 year
observational period, we have ${\dot{f}/}{\Delta\dot{f}}=$ 1.0 and 7.4,
respectively. These results are summarized in the upper part of Table 1. When
the observing time is increased to 10 years, we obtain
${\dot{f}/}{\Delta\dot{f}}=$ 14 and 110.

Therefore, for a 10 year observational period, the majority of low
eccentricity DNSBs will have a measurable $\dot{f}$ accurate to within 10%.
Over 4 years, when $f>$2mHz, $\dot{f}$ is accurate within 15%. The fractional
chirp mass resolution follows

$\frac{\Delta\mathcal{M}}{\mathcal{M}}\simeq\frac{3}{5}\frac{\Delta\dot{f}}{\dot{f}}.$
(9)

It will be possible to select Galactic DNSBs based on the chirp masses
(expected to be narrowly distributed around $\sim 1.2M_{\odot}$), distinct
from the much more numerous WDBs in most cases. This is due to the expected
rarity of mHz WDBs with mass $\mathcal{M}>1M_{\odot}$ (see e.g., Fig. 1 in
[23]). We will later discuss potential issues related to binaries including
high mass white dwarfs.

####  II.2.2 Eccentricity

Known Galactic DNSBs could have residual eccentricities $e\sim 0.1$ by the
time they enter the LISA band [34], even despite the tendency to rapidly
circularize by Equations (2) and (3) as they evolve towards mHz frequencies.

Table 1: Summary of detection properties of the Hulse-Taylor like DNSB. The binary is conservatively located at 20kpc, and has a chirp mass of $\mathcal{M}=1.2M_{\odot}$.  Property | Hulse–Taylor like pulsar
---|---
$\mathcal{M}$ | 1.2 $M_{\odot}$
${D}$ | 20 kpc
SNR(1.5mHz, 4yr) | 7.8
SNR(50mHz, 4yr) | 452
SNR(1.5mHz, 10yr) | 17
SNR(50mHz, 10yr) | 715
${\dot{f}}/{\Delta\dot{f}}$(1.5mHz, 4yr) | 1.0
${\dot{f}}/{\Delta\dot{f}}$(2mHz, 4yr) | 7.4
${\dot{f}}/{\Delta\dot{f}}$(1.5mHz, 10yr) | 14
${\dot{f}}/{\Delta\dot{f}}$(2mHz, 10yr) | 110
$e$(1.5mHz)=$e_{\mathrm{i}}$ | 0.057
$e$(50mHz)=$e_{\mathrm{f}}$ | 0.0014
$T_{\mathrm{merge,\ i}}$ | 486 kyr
$T_{\mathrm{merge,\ f}}$ | 42 yr
$T_{\mathrm{band}}$ | 486 kyr

The Hulse–Taylor binary pulsar [1] (HT) has an eccentric ($e=0.6$) 7.8 hour
orbit ($f=$ 0.07mHz), and a chirp mass of $1.2M_{\odot}$. If we consider an
HT–like pulsar with today’s properties and evolve it from gravitational wave
emission through Equations (2) and (3), it will have eccentricity $e=0.057$
upon entering the LISA band at 1.5mHz (see the lower part of Table 1).

For such binaries with small eccentricities $e\ll 1$, whose strain in the
second orbital harmonic $f_{\mathrm{GW}}=2f_{\mathrm{orb}}(=f)$ is given by
Equation (6), the sub-leading strains ($\propto e$) will be present in the
first and third orbital harmonics (at $f_{\mathrm{GW}}=1f_{\mathrm{orb}}$ and
$3f_{\mathrm{orb}}$). Using the same technique to obtain Equation (6), the
strain in the third harmonic at $f_{\mathrm{GW}}=3f_{\mathrm{orb}}$ is
approximately

$h_{3}=\frac{9e}{4}h_{2}.$ (10)

Ignoring changes in $S_{n}(f)$ at $f_{\mathrm{GW}}=3f_{\mathrm{orb}}$ compared
to $f_{\mathrm{GW}}=2f_{\mathrm{orb}}$, the signal–to–noise ratio in the third
orbital harmonic $\bar{\rho}_{3}$ is then:

$\bar{\rho}_{3}\approx\frac{9e}{4}\bar{\rho}_{2}.$ (11)

For $(f,e)=(1.5\mathrm{mHz},0.057)$, the strain in the third orbital harmonic
($h_{3}$) is 13 percent of the dominant harmonic ($h_{2}$) via Equation (10).
If conservatively located at 20kpc, using Equations (7) and (11) we get
$\bar{\rho}_{3}=1$ for a 4yr observation. Therefore, this eccentric binary
will not produce a detectable $\bar{\rho}_{3}$. At $50$mHz, the eccentricity
has decreased by a factor of 40 ($e=0.0014$), but the signal-to-noise has
increased by a factor of 60 (see Table 1). This gives $\bar{\rho}_{3}={1.4}$,
which is also not detectable. However, at intermediate frequencies between the
frequency boundaries, the HT binary’s eccentricity is marginally detectable.
For example, at $f=15$mHz, $\bar{\rho}_{3}=5$.

On the $f$-$e$ plane, we define “the detectable region” as the area where we
can make a complete survey for Galactic DNSBs. So far, we have mainly
discussed DNSBs with relatively small residual eccentricities in the LISA
band. The detection for these binaries will be complete at $f\gtrsim 1.5$mHz
for an observational period of $\sim 10$ yr. Meanwhile, we are particularly
interested in the DNSBs which formed in the LISA band. If they have large
eccentricities, their higher harmonic strains will be essential for detecting
them.

However, on the $f$-$e$ plane, we simply put the detectable region as the
rectangular area which is bounded by the inequalities $f_{\mathrm{i}}\equiv
1.5{\rm mHz}\leq f\leq f_{\mathrm{f}}\equiv{50}{\rm mHz}$ and $0\leq e<1$,
ignoring the eccentricity dependence for the boundary frequencies $f_{\rm i}$
and $f_{\rm f}$. Our choice $f_{\rm f}=50$mHz is somewhat arbitrary, and most
of our results below are almost independent of it.

In Fig. 1, the detectable region is shown as the colorful rectangular area.
While it will not be difficult to more precisely include the eccentricity
dependence of the boundaries, this task is beyond the scope of our conceptual
study. We hereafter put $N$ as the total number of DNSBs in the detectable
region.

Figure 1: Evolution of Galactic DNSBs on the $f$-$e$ plane. The detectable
region is shown in pastel colors, bounded by $f_{\rm i}=1.5$mHz and $f_{\rm
f}={50}$mHz. The black curves are the flow lines from Eq. (12), moving
rightward in time. For each point, using the associated flow line, we can
define the effective time fraction $\tau(f,e)$, which is shown as a contour
map. We have $\tau=0$ and 1 at the two boundary frequencies. The red curve
corresponds to $\tau=0.5$ (halfway in this dimensionless time unit). The
purple point on the second flow line corresponds to a binary with $(f,e)=$
(5mHz, 0.5).

##  III flow lines in the detectable region

As discussed in the previous section, LISA will provide us with Galactic DNSBs
scattered in the detectable region on the $f$-$e$ plane. Our primary interest
in this paper is to detect the signature of DNSB formation in the LISA band.
If all of the detected DNSBs were formed at low frequencies (much lower than
mHz) and resultantly have low eccentricities $e\ll 1$, the frequency
distribution $\mathrm{d}N/\mathrm{d}f$ would be the appropriate data to be
analyzed. Indeed, in such a case, we will have the profile
$\mathrm{d}N/\mathrm{d}f\propto f^{-11/3}$ as shown in Eq. (4).

However, we might detect DNSBs with non-negligible eccentricities. In
particular, those formed around the mHz band could have $e=O(1)$. For our
study, we thus need to deal with the two dimensional distribution
$\mathrm{d}^{2}N/\mathrm{d}f\mathrm{d}e$. It is nevertheless advantageous to
compress the two dimensional data into much simpler one dimensional data. In
fact, there are many useful tools to statistically examine one-dimensional
patterns, as we see later.

The basic question here is whether we have a data compression scheme suitable
for studying the potential mHz formation for DNSBs. In this section, after
discussing the long-term evolution of DNSBs in the mHz band, we formally
define the effective time fraction $\tau(f,e)$ as a function of $f$ and $e$,
for the data compression.

###  III.1 Binary evolution in the mHz band

First, we consider the long-term evolution of DNSBs in the mHz band due to the
GW emission. From Eqs. (2) and (3), we have

$\frac{\mathrm{d}\ln
f}{\mathrm{d}e}=-\frac{18}{19}\frac{(1+73/24e^{2}+37/96e^{4})}{e(1-e^{2})(1+121/304e^{2})}.$
(12)

This equation is independent of the mass parameters. As mentioned earlier,
except for $1-e\ll 1$, the relativistic correction is small for DNSBs in the
LISA band, and this equation is an excellent approximation. We can easily
integrate Eq. (12) and obtain the flow lines as shown in Fig. 1 with the black
lines. Due to the structure of Eq. (12), this family satisfies a self-similar
relation.

For a given point in the detectable region, we can identify an associated flow
line. For example, with respect to the point $(f,e)=$ (5mHz, 0.5) shown in
Fig. 1, we have the flow line that enters the detectable region at $(f_{\rm
i},e_{\rm i})=(1.5{\rm mHz},0.74)$ and escapes it at $(f_{\rm f},e_{\rm
f})=(50{\rm mHz},0.073)$. Similarly, as presented in Table 1, the HT-like
binary would have $e_{\rm i}=0.057$ and $e_{\rm f}=0.0014$.

###  III.2 Effective time fraction on a flow line

From Eqs. (2) and (3), we can also estimate the remaining time before the
merger [26]

$\begin{split}&T_{\mathrm{merge}}(f,e)=\frac{15}{304}\frac{c^{5}}{\pi^{8/3}G^{5/3}}\mathcal{M}^{-5/3}\\\
&\times
f^{-8/3}\left(\frac{1-e^{2}}{e^{12/19}\left(1+\frac{121}{304}e^{2}\right)^{870/2299}}\right)^{4}\\\
&\times\int_{0}^{e}\frac{e^{\prime 29/19}\left(1+\frac{121}{304}e^{\prime
2}\right)^{1181/2299}}{\left(1-e^{\prime 2}\right)^{3/2}}\
\mathrm{d}e^{\prime}\,\end{split}$ (13)

as a function of $f$ and $e$. Then, for a given DNSB at $(f,e)$ in the
detectable region, we can evaluate the total time $T_{\mathrm{total}}(f,e)$
that the binary spends on its flow line in the detectable region (i.e. from
the minimum frequency $f_{\rm i}$ to the maximum one $f_{\rm f}$). More
specifically, we can put

$T_{\mathrm{total}}(f,e)=T_{\mathrm{merge}}(f_{\rm
i},e_{\mathrm{i}}(e,f))-T_{\mathrm{merge}}(f_{\rm f},e_{\mathrm{f}}(e,f)),$
(14)

where the eccentricities $e_{\rm i}$ and $e_{\rm f}$ at the two boundary
frequencies should be regarded as functions of $f$ and $e$ (through the
corresponding flow line). For $f_{\rm f}=50$mHz, the second term in Eq. (14)
is virtually ignorable. For the HT-like binary in Table 1, we obtain
$T_{\mathrm{total}}=$ 486 kyr. Similarly, the time elapsed after entering into
the detectable region at $f=f_{\rm i}$ is given by

$T_{\mathrm{band}}(f,e)=T_{\mathrm{merge}}(f_{\rm
i},e_{\mathrm{i}}(f,e))-T_{\mathrm{merge}}(f,e).$ (15)

Then, we define the effective time fraction

$\tau(f,e)=\frac{T_{\mathrm{band}}(f,e)}{T_{\mathrm{total}}(f,e)}$ (16)

for characterizing the position of a DNSB on its flow line in the detectable
region. We note that the effective time fraction $\tau$ is independent of the
chirp mass $\mathcal{M}$.

In Fig. 1, we show the contour levels for the effective time fraction $\tau$.
We have $\tau=0$ at the lower bound $f=f_{\rm i}$ and $\tau=1$ at the upper
end $f=f_{\rm f}$. The DNSB at $(f,e)=$ (5mHz, 0.5) has $\tau=0.77$, shown by
the purple point. At $e\lesssim 0.05$, the function $\tau(f,e)$ depends very
weakly on the eccentricity $e$, as understood from the weak correction
$O(e^{2})$ in Eq. (2). We can approximately put

${\tau\approx\frac{f_{\mathrm{i}}^{-8/3}-f^{-8/3}}{f_{\mathrm{i}}^{-8/3}-f_{\mathrm{f}}^{-8/3}}\approx
1-\left(\frac{f}{1.5\mathrm{mHz}}\right)^{-8/3}.}$ (17)

at $e\ll 1$. In the high eccentricity limit $e\to 1$, we have

${\tau\approx{1.45}\left[1-\left(\frac{f}{1.5\mathrm{mHz}}\right)^{-1/3}\right].}$
(18)

For the detectable region in Fig. 1, we have simply ignored the eccentricity
dependence of the upper and lower frequency boundaries. Even if the detectable
region is deformed on the $f$-$e$ plane, it is straightforward to evaluate the
effective time fraction $\tau(f,e)$ by using a family of the flow lines.

##  IV statistical tests for in-band binary formation

In this section, we discuss how to examine DNSB injections by using the
cumulative distribution function (CDF) of the observed time fractions $\tau$.
First, we derive theoretical expressions for the CDFs without injections and
with injections (based on a simple model). In reality, due to the finiteness
of the sample size $N$, the observed CDF will have fluctuations around the
theoretical curves. We discuss statistical tests to determine the potential
injections, under the presence of these fluctuations.

###  IV.1 Without injections

Without injections, each binary goes through the detectable region on a flow
line at a constant speed $\mathrm{d}\tau/\mathrm{d}t=$ const. As long as the
DNSBs are in a steady state, after counting contributions from all flow lines,
we have a uniform probability distribution function (PDF) for $\tau$

${\rm Pr}(\tau)={\rm const}.$ (19)

From the normalization condition $\int_{0}^{1}\mathrm{d}\tau{\rm Pr}(\tau)=1$,
we obtain

${\rm Pr}(\tau)=1.$ (20)

The corresponding CDF is given as

$F(\tau)=\int_{0}^{\tau}\mathrm{d}\tau^{\prime}{\rm Pr}(\tau^{\prime})=\tau.$
(21)

Figure 2:  The CDFs for the effective time fraction $\tau$, for the cases of
$a=0.8$ and $a=1.8$ are shown in the upper and lower panels, respectively. The
blue curves show $F_{0}(\tau)=\tau$ without injections (uniform distribution).
The dotted red curves show $F_{a}(\tau)$ (see Eq. (24)) with the injections
parameters $a=0.8$ and 1.8 (or the injected fractions $B=0.29$ and 0.47). The
solid red curves are the corresponding numerical realizations $F(\tau)_{R1}$
and $F(\tau)_{R2}$ with the sample size $N=200$. The maximum deviations $d$,
which are used in the Kolmogorov–Smirnov (KS) test, are given by the black
arrows. The respective KS $p$-value are 0.054 and $4.7\times 10^{-4}$.
Generally speaking, increasing $a$ leads to a larger deviation $d$ and hence
smaller $p_{\rm KS}$, while a larger $N$ gives a curve more similar to the
theoretical CDF $F_{a}(\tau)$. The $p$-values are also provided for the
Anderson-Darling (AD) test.

###  IV.2 With injections

Now we discuss the PDFs and CDFs for populations with injections. As a simple
model, we assume a constant injection rate over $\tau=$ [0,1] described as

$\frac{\mathrm{d}{\rm Pr}(\tau)}{\mathrm{d}\tau}={\rm const}.$ (22)

After integrating this equation and normalizing the result, we obtain

$\mathrm{Pr}(\tau)=\frac{1+a\tau}{1+a/2},$ (23)

where the positive parameter $a$ is related to the injection rate. In the
numerator of Eq. (23), the first and second terms respectively show the
component formed below $f_{\rm i}$ and the component injected in the band
$[f_{\rm i},f_{\rm f}]$. The associated CDF is given as

$F_{a}(\tau)=\frac{\tau+a\tau^{2}/2}{1+a/2},$ (24)

and is convex downward. Here we added the subscript $a$ to explicitly show the
parameter dependence. For $a=0$, we recover Eq. (21) without injections.

Note that we do not need to directly deal with detailed models for the
eccentricity dependence of the injections. This is another advantage of using
the compressed variable $\tau$.

From Eq. (24) we can easily confirm that, in the detectable region, the
fraction of injected DNSBs in the whole population is given by

${B}=\frac{a/2}{1+a/2}.$ (25)

Since this parameter is more comprehensive than the original one $a$, we use
them in parallel. Note that for our simple model (23), we have

${B}=\frac{1}{2}\frac{\mathrm{d}{\rm Pr}(\tau)}{\mathrm{d}\tau}.$ (26)

In Figure 2, with the red dotted curves, we show the two CDFs $F_{0.8}(\tau)$
(upper) and $F_{1.8}(\tau)$ (lower). We have the injected fractions $B=0.29$
and 0.47 respectively. We also present the CDF $F_{0}(\tau)=\tau$ without
injections (blue curves).

###  IV.3 Statistical tests

For a nonzero injection parameter $a$, the analytical function $F_{a}(\tau)$
is distinct from $F_{0}(\tau)$ without injections. However, we should recall
that the expected number $N$ of the DNSBs in the detectable region is
10’s–100’s. If we make the CDF for the observed time fractions $\tau$ of this
small sample size, we will have significant scatter due to the finiteness of
$N$.

To be concrete, we perform a numerical experiment for $a=0.8$ with the sample
number $N=200$. We employ the scheme known as inverse transform sampling, and
generate a list $\\{\tau_{\rm i}\\}$ (${\rm i}=1,\cdots,N$) whose individual
elements are drawn from the analytic CDF $F_{0.8}(\tau)$ without measurement
errors. After sorting, we obtain the CDF $F(\tau)_{R1}$, as shown by the solid
red line in the upper panel of Fig. 2. Similarly, we generate another
realization $F(\tau)_{R2}$ for $N=200$ and $a=1.8$ (given in the lower panel).
As expected, we can observe fluctuations around the original CDFs
$F_{0.8}(\tau)$ and $F_{1.8}(\tau)$ (Eq. 24) presented with the dotted red
curves.

Now, let us assume that LISA provides us with the sample corresponding to
$F(\tau)_{R1}$. Our central task here is to statistically determine how likely
this realization could have been drawn from the model CDF $F_{0}(\tau)$
without injections (the null hypothesis). If unlikely, it would be reasonable
to claim that we detect a signature of injections.

There are many sophisticated statistical techniques to check the consistency
of a data sample with respect to a reference CDF. In the Kolmogorov-Smirnov
(KS) test [35], we evaluate the maximum deviation between the two CDFs. In our
case, it is expressed as

$d=\max_{0\leq\tau\leq 1}|F(\tau)_{R1}-F_{0}(\tau)|.$ (27)

For an observed sample with size $N$, the probability of obtaining a maximum
deviation $\geq d$ is approximately given by [35]

$\begin{split}&p_{\mathrm{KS}}\approx\\\
&2\sum_{j=1}^{\infty}(-1)^{j-1}\mathrm{Exp}\left\\{-2j^{2}\left[\left(N_{\mathrm{}}^{1/2}+0.12+\frac{0.11}{N_{\mathrm{}}^{1/2}}\right)d\right]^{2}\right\\}.\end{split}$
(28)

We use the probability $p_{\rm KS}$ as the $p$-value for the sample
$F(\tau)_{R1}$ to be drawn from the CDF $F_{0}(\tau)$ without injections.

As shown in the upper panel of Fig. 2, for the realization $F(\tau)_{R1}$, we
have $d=0.094$ and $p_{\rm KS}=0.054$ (using Mathematica). Similarly, for the
other realization $F(\tau)_{R2}$ in the lower panel, we have $d=0.14$ and
$p_{\rm KS}=4.7\times 10^{-4}$. Note that, in each panel, due to the
statistical fluctuations, the black deviation $d$ is larger than the maximum
deviation between the two theoretical curves at $\tau=0.5$. These are
respectively given by ${a}/{(4a+8)}=B/4=0.071$ and 0.12.

The KS test is a simple method, and it is most sensitive to the data around
$\tau=0.5$, where we typically have the maximum deviation $d$. In contrast,
the Anderson-Darling (AD) test uses the whole range of a CDF [35], and often
provides us with a more stringent $p$-value $p_{\rm AD}<p_{\rm KS}$. For the
two realizations in Fig. 2, we have $p_{\rm AD}=0.018$ and $7.4\times 10^{-5}$
respectively (again using Mathematica). Interestingly, for the $a=0.8$ case
the upper sample is now below the standard threshold $p=0.05$.

###  IV.4 Systematic Study

We now systematically explore the parameter space $(N,B)$. Given the current
uncertainties of these parameters, we consider the range $0<N<500$ and
$0<B<0.5$, as shown in Fig. 3, and divide this space with a 20$\times$20 grid.

At each grid point, we simulate 1000 realizations for the $\tau$ distribution,
evaluate the $p$-values $p_{\rm KS}$ individually for the realizations, and
take their median value $\langle p\rangle_{\mathrm{KS}}$. The numerical
results $\langle p\rangle_{\mathrm{KS}}$ are shown in Figure 3 by the pastel
levels.

Figure 3: Contour plot of the median $p$-value $\log_{10}[\langle
p\rangle_{\mathrm{KS}}]$ for the Kolmogorov-Smirnov (KS) test. This is
calculated for sample size $N$ and injected fraction $B$ with a 20$\times$20
resolution. The black curve shows the boundary for $\langle
p\rangle_{\mathrm{KS}}=0.05$ (or $\log_{10}[\langle
p\rangle_{\mathrm{KS}}]=-1.3$). The red curve shows a similar boundary for the
Anderson-Darling (AD) test $\langle p\rangle_{\mathrm{AD}}=0.05$. Above these
curves, it is likely that the no injection scenario will be rejected due to a
$p$-value less than 0.05.

To determine the characteristic curve corresponding to $\langle
p\rangle_{\mathrm{KS}}=0.05$ in Fig. 3, we made iterative calculations and
obtained an approximate expression

$B=0.37\left(\frac{N}{100}\right)^{-1/2}.$ (29)

We show this expression with the black curve. For example, at $B=0.37$, we
need $N\sim 100$, which is 3.3 times larger than the reference value for the
expected number of Galactic DNSBs in Eq. (5).

In the same manner, we obtain an expression for the AD test ($\langle
p\rangle_{\mathrm{AD}}=0.05$) as

$B=0.35\left(\frac{N}{100}\right)^{-1/2}$ (30)

presented with the red curve in Fig. 3. For a given injected fraction $B$, the
AD test requires $\sim 10$% smaller sample size $N$ than the KS test.

###  IV.5 Parameter Estimation

So far, we have studied the statistical tests to check potential DNSB
injections. Here we briefly discuss how well we can estimate the injected
fraction $B$ (or equivalently $a$) from the CDF of the observed effective time
fractions $\tau$. The standard approach for the estimation is to compare the
overall profile of the observed CDF with the theoretical expression (24) and
fit the parameter $B$ (see e.g., the solid and dotted red curves in Fig. 2).
Roughly speaking, this method is similar to the concept of the AD test.
Instead, below, we examine a simple method only using the maximum deviation
$d$, as in the case of the KS test.

As mentioned earlier, the maximum deviation $d$ has statistical fluctuations,
due to the finiteness of the sample size $N$. For a given $B$ and $N$, we put
the mean value of $d$ by ${\bar{d}}(B,N)$ and the root-mean-square (rms)
scatter by $\sigma_{d}(B,N)$. If we try to estimate the parameter $B$ from the
observed maximum deviation $d$ (inversely using the relation
$d={\bar{d}}(B,N)$), the estimated parameter $B$ contains the rms error
$\sigma_{B}(B,N)$ roughly given by

${\sigma_{B}}(B,N)={\sigma_{d}}(B,N){\left\\{\frac{\partial{\bar{d}}(B,N)}{\partial
B}\right\\}}^{-1}$ (31)

with the Jacobian ${\partial{\bar{d}}(B,N)}/{\partial B}$.

At various points in Fig. 3, we numerically evaluated the two factors
$\sigma_{d}$ and ${\partial{\bar{d}}(B,N)}/{\partial B}$ by generating a large
number of realizations. For the scatter $\sigma_{d}$, we found an approximate
relation

$\sigma_{d}(B,N)\simeq 0.04\left(\frac{N}{100}\right)^{-1/2},$ (32)

which is independent of $B$. By taking finite differences instead of the
derivatives, we also found

$\frac{\partial{\bar{d}}(B,N)}{\partial B}\simeq 0.25.$ (33)

This simple relation seems reasonable, given the maximum deviation $B/4$
between the two theoretical curves in Fig. 2 (as mentioned earlier). Then, for
the estimation error of $B$, we have

${\sigma_{B}}(B,N)\simeq 0.16\left(\frac{N}{100}\right)^{-1/2}.$ (34)

Comparing this result with the critical curves given by Eqs. (29) and (30), we
can see that these curves roughly correspond to the condition $B\sim
2\sigma_{B}$.

##  V Discussion

In this paper we have discussed the basic idea of probing DNSB injections in
the detectable region on the $(f,e)$ plane. Here we mention potential
extensions of our study.

###  V.1 Other projections

For simplicity, we proposed the projection to the single variable $\tau$.
Considering the dimensionality of the original data, we might develop other
projection methods, as illustrated in the following example.

Using the variable $\tau$ in Fig. 1, we first divide the DNSBs into two
groups: $G^{-}$ (with $0\leq\tau\leq 0.5$) and $G^{+}$ (with $0.5<\tau\leq
1.0$). Then, along the flow lines, we move the DNSBs in the group $G^{-}$, to
the boundary curve $\tau=0.5$ (the red line in Fig. 1). We put
$P^{-}_{0.5}(e)$ as the resultant eccentricity distribution at $\tau=0.5$. In
the same manner, we can obtain $P_{0.5}^{+}(e)$ for the group $G^{+}$.

Now, without the DNSB injections, the two profiles $P_{0.5}^{-}(e)$ and
$P_{0.5}^{+}(e)$ should be the same, except for some statistical fluctuations.
With injections, this is generally not the case. Therefore, we can apply
various techniques, such as the two sample KS-test [35], for probing potential
injections.

###  V.2 Lower frequency regime

We have dealt with DNSBs only in the detectable region at $f\geq 1.5$mHz, by
evaluating the signal-to-noise ratios at a conservative distance of 20kpc.
However, it would be highly desirable to probe the injections in the lower
frequency regime. In the extended regime, LISA’s Galactic DNSB survey is not
complete. For applying our method, we thus need to correct for the selection
bias, by using appropriate models for the distribution of Galactic DNSBs.

###  V.3 Massive white dwarfs

PSR B2303+46 [36, 37] and J0453+1559 [38] are Galactic neutron star white
dwarf (NS-WD) binary candidates with relatively large chirp masses of
$\mathcal{M}=1.05M_{\odot}$ and $\mathcal{M}=$1.17$M_{\odot}$ respectively.
While they have orbital periods of days (merger timescales larger than Hubble
time), we might actually detect similar NS-WD systems with DNSB–like chirp
masses in the mHz band. This is because NS-WD are likely to be the next most
abundant compact binary after WDBs [39]. A binary’s white dwarf component
could be confirmed with targeted optical/IR followup, but these searches may
not be successful. Observational challenges include interstellar extinction,
and limitations in the sky localisation with LISA. In relation to the
electromagnetic wave observation, radio followup will be also important for
estimating the age of a DNSB [7]. Here, we should note that a dynamically
formed DNSB may have an age largely different from the spin-down age of the
component pulsars.

In any case, we can rather decrease the chirp mass threshold to add NS-WDs and
WDBs to our binary sample, and thereby increase the total number $N$. However,
in Fig. 1, the flow lines and the associated time fraction $\tau$ are based on
the point particle approximation in Eqs. (2) and (3). We thus need to
carefully examine the possible modifications induced by the finiteness of
white dwarfs [40, 41].

##  VI Summary

Due to the complex interplay between various physical processes, it is still
difficult to solidly predict the orbital parameters of compact binaries at
their formation. New observational results will thus help us to refine
theoretical modelings. In the near future, the space interferometer LISA will
explore GWs around 0.1-100mHz and will make a complete survey for Galactic
DNSBs above $\sim 1.5$mHz after $\sim$10yr operation.

In this paper, we discussed how well we can detect a signature of DNSB
formation (injections) around 1mHz with LISA. We first introduced the
effective time fraction $\tau$ for each DNSB to compress the original two-
dimensional data $(f,e)$, as shown in Fig. 1. The probability distribution
function ${\rm Pr(\tau)}$ for the measured fractions $\tau$ plays the central
role in our method. There was a flat profile $\mathrm{d}{\rm
Pr}(\tau)/\mathrm{d}\tau=0$ without injections, but $\mathrm{d}{\rm
Pr}(\tau)/\mathrm{d}\tau>0$ with injections.

In reality, we need to discriminate the differences between the profiles under
the existence of scatter due to the finite sample size $N$. To be concrete, we
made a simple model for the DNSB injections with $\mathrm{d}{\rm
Pr}(\tau)/\mathrm{d}\tau=2B(={\rm const})$, characterized by the injection
fraction $B$. Then, we examined the prospects for discriminating the profile
differences with the Kolmogorov–Smirnov and Anderson–Darling tests. Our main
results are presented in Fig. 3 with the characteristic relations (29) and
(30) for the $p$-value of 0.05. For example, with the sample size $N=100$, we
need an injection fraction $B\gtrsim 0.35$ to detect injections.

In this paper, we have discussed the very basic idea of studying potential
DNSB formation around 1mHz. Our approach can be extended in various
directions, including those mentioned in section V.

## Acknowledgements

LM acknowledges financial support from the Japan Society for the Promotion of
Science (JSPS) International Research Fellow program (Graduate School of
Science, Kyoto University JSPS P21017). NS is supported by JSPS Kakenhi Grant-
in-Aid for Scientific Research (Nos. 17H06358 and 19K03870). Mathematica was
used to perform calculations and generate figures in this work.

## References

  * Hulse and Taylor [1975] R. A. Hulse and J. H. Taylor, ApJL 195, L51 (1975).
  * Manchester and Taylor [1977] R. N. Manchester and J. H. Taylor, _Pulsars_ (Freeman, San Francisco, 1977).
  * Stovall et al. [2018] K. Stovall, P. C. C. Freire, S. Chatterjee, P. B. Demorest, D. R. Lorimer, M. A. McLaughlin, N. Pol, J. van Leeuwen, R. S. Wharton, B. Allen, et al., The Astrophysical Journal Letters 854, L22 (2018), eprint 1802.01707.
  * Swiggum et al. [2015] J. K. Swiggum, R. Rosen, M. A. McLaughlin, D. R. Lorimer, S. Heatherly, R. Lynch, S. Scoles, T. Hockett, E. Filik, J. A. Marlowe, et al., Astrophys. J. 805, 156 (2015), eprint 1503.06276.
  * Bagchi et al. [2013] M. Bagchi, D. R. Lorimer, and S. Wolfe, MNRAS 432, 1303 (2013), eprint 1302.4914.
  * Amaro-Seoane et al. [2017] P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender, E. Berti, P. Binetruy, M. Born, D. Bortoluzzi, et al., arXiv e-prints arXiv:1702.00786 (2017), eprint 1702.00786.
  * Kyutoku et al. [2019] K. Kyutoku, Y. Nishino, and N. Seto, Monthly Notices of the Royal Astronomical Society 483, 2615 (2019), eprint 1812.02177.
  * Lau et al. [2020] M. Y. M. Lau, I. Mandel, A. Vigna-Gómez, C. J. Neijssel, S. Stevenson, and A. Sesana, MNRAS 492, 3061 (2020), eprint 1910.12422.
  * Wagg et al. [2021] T. Wagg, F. S. Broekgaarden, S. E. de Mink, L. A. C. van Son, N. Frankel, and S. Justham, arXiv e-prints arXiv:2111.13704 (2021), eprint 2111.13704.
  * Paczynski [1976] B. Paczynski, in _Structure and Evolution of Close Binary Systems_ , edited by P. Eggleton, S. Mitton, and J. Whelan (1976), vol. 73 of _IAU Symposium_ , p. 75.
  * Podsiadlowski et al. [1992] P. Podsiadlowski, P. C. Joss, and J. J. L. Hsu, Astrophys. J. 391, 246 (1992).
  * Vink et al. [2001] J. S. Vink, A. de Koter, and H. J. G. L. M. Lamers, AAP 369, 574 (2001), eprint astro-ph/0101509.
  * Dewi et al. [2005] J. D. M. Dewi, P. Podsiadlowski, and O. R. Pols, MNRAS 363, L71 (2005), eprint astro-ph/0507628.
  * Mandel and Broekgaarden [2022] I. Mandel and F. S. Broekgaarden, Living Reviews in Relativity 25, 1 (2022), eprint 2107.14239.
  * Phinney [1991] E. S. Phinney, The Astrophysical Journal Letters 380, L17 (1991).
  * Kremer et al. [2018] K. Kremer, S. Chatterjee, K. Breivik, C. L. Rodriguez, S. L. Larson, and F. A. Rasio, Physical Review Letters 120, 191103 (2018), eprint 1802.05661.
  * Andrews et al. [2020] J. J. Andrews, K. Breivik, C. Pankow, D. J. D’Orazio, and M. Safarzadeh, ApJL 892, L9 (2020), eprint 1910.13436.
  * Korol et al. [2019] V. Korol, E. M. Rossi, and E. Barausse, MNRAS 483, 5518 (2019), eprint 1806.03306.
  * Wilhelm et al. [2021] M. J. C. Wilhelm, V. Korol, E. M. Rossi, and E. D’Onghia, MNRAS 500, 4958 (2021), eprint 2003.11074.
  * Seto [2022] N. Seto, Phys. Rev. Lett. 128, 041101 (2022), eprint 2201.03685.
  * Kyutoku and Seto [2016] K. Kyutoku and N. Seto, MNRAS 462, 2177 (2016), eprint 1606.02298.
  * Nelemans et al. [2004] G. Nelemans, L. R. Yungelson, and S. F. Portegies Zwart, MNRAS 349, 181 (2004), eprint astro-ph/0312193.
  * Lamberts et al. [2019] A. Lamberts, S. Blunt, T. B. Littenberg, S. Garrison-Kimmel, T. Kupfer, and R. E. Sanderson, MNRAS p. 2426 (2019), eprint 1907.00014.
  * The LIGO Scientific Collaboration et al. [2021] The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, N. Adhikari, R. X. Adhikari, et al., arXiv e-prints arXiv:2111.03634 (2021), eprint 2111.03634.
  * Kalogera et al. [2001] V. Kalogera, R. Narayan, D. N. Spergel, and J. H. Taylor, ApJ 556, 340 (2001), eprint astro-ph/0012038.
  * Peters [1964] P. C. Peters, Phys. Rev. 136, B1224 (1964).
  * 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, et al., The Astrophysical Journal 846, 170 (2017), eprint 1706.09438.
  * Robson et al. [2019] T. Robson, N. J. Cornish, and C. Liu, Classical and Quantum Gravity 36, 105011 (2019), eprint 1803.01944.
  * Flanagan and Hughes [1998] É. É. Flanagan and S. A. Hughes, Phys. Rev. D. 57, 4535 (1998), eprint gr-qc/9701039.
  * Finn and Thorne [2000] L. S. Finn and K. S. Thorne, Phys. Rev. D. 62, 124021 (2000), eprint gr-qc/0007074.
  * Peters and Mathews [1963] P. C. Peters and J. Mathews, Physical Review 131, 435 (1963).
  * Seto [2019] N. Seto, Monthly Notices of the Royal Astronomical Society 489, 4513 (2019), eprint 1909.01471.
  * Takahashi and Seto [2002] R. Takahashi and N. Seto, The Astrophysical Journal 575, 1030 (2002).
  * Brown et al. [2001] G. E. Brown, C. H. Lee, S. F. Portegies Zwart, and H. A. Bethe, Astrophys. J. 547, 345 (2001), eprint astro-ph/9911130.
  * Press et al. [1986] W. H. Press, B. P. Flannery, and S. A. Teukolsky, _Numerical recipes. The art of scientific computing_ (Cambridge University Press, 1986).
  * van Kerkwijk and Kulkarni [1999] M. H. van Kerkwijk and S. R. Kulkarni, ApJL 516, L25 (1999), eprint astro-ph/9901149.
  * Tauris and Janka [2019] T. M. Tauris and H.-T. Janka, ApJL 886, L20 (2019), eprint 1909.12318.
  * Martinez et al. [2015] J. G. Martinez, K. Stovall, P. C. C. Freire, J. S. Deneva, F. A. Jenet, M. A. McLaughlin, M. Bagchi, S. D. Bates, and A. Ridolfi, ApJ 812, 143 (2015), eprint 1509.08805.
  * Nelemans et al. [2001] G. Nelemans, L. R. Yungelson, S. F. Portegies Zwart, and F. Verbunt, AAP 365, 491 (2001), eprint astro-ph/0010457.
  * McNeill et al. [2020] L. O. McNeill, R. A. Mardling, and B. Müller, MNRAS 491, 3000 (2020), eprint 1901.09045.
  * Wolz et al. [2021] A. Wolz, K. Yagi, N. Anderson, and A. J. Taylor, MNRAS 500, L52 (2021), eprint 2011.04722.


Automated interpretability research has recently attracted attention as a potential research direction that could scale explanations of neural network behavior to large models.
Existing automated circuit discovery work applies activation patching to identify subnetworks responsible for solving specific tasks (circuits).
In this work, we show that a simple method based on attribution patching outperforms all existing methods while requiring just two forward passes and a backward pass.
We apply a linear approximation to activation patching to estimate the importance of each edge in the computational subgraph. Using this approximation, we prune the least important edges of the network. We survey the performance and limitations of this method, finding that averaged over all tasks our method has greater AUC from circuit recovery than other methods.

§ INTRODUCTION

Mechanistic interpretability is a subfield of AI interpretability that focuses on attributing model behaviors to its components, thus reverse engineering the network <cit.>. This field aims to identify subnetworks (circuits) within the model which are responsible for solving specific tasks <cit.>. Prior attempts at finding circuits in language models have led to finding networks of attention heads and multi-layer perceptrons (MLPs) that partially or fully explain model behaviors at tasks such as indirect object identification, modular arithmetic, completion of docstrings, and predicting successive dates <cit.>. However, almost all previous work has been limited to relatively small models since manually applying mechanistic interpretability methods has not currently scaled to end-to-end circuits in larger models <cit.>.

It may be important to scale interpretability to large models as these are the neural networks most widely deployed and used by a wide range of people. Currently, we have little understanding into these models work and failure modes are not always found ahead of deployment. If successful, scaled interpretability could address a wide variety of concerns about the lack of transparency of language models <cit.>, in addition to speculative risks about the alignment of machine learning systems <cit.>.

Automated Circuit Discovery (ACDC; <cit.>) attempts to automate a large portion of the mechanistic interpretability workflow — the pruning of edges between attention heads and MLPs that do not affect the task being studied. ACDC begins with a computational graph, and recursively calculates the importance of an edge in the graph for a specific task. In our work, we use edges to refer to activations inside models between two components (<Ref> describes this motivation further). ACDC's pruning algorithm applies activation patching. (Note that activation patching is not attribution patching. Both are defined in full in <Ref>.) At a high level, activation patching edits a specific activation in a model forward pass and measures a model statistic (e.g loss) under this intervention. Activation patching is inefficient for circuit discovery because getting each statistic about model activations requires another forward pass. Our work uses attribution patching to recover circuits more efficiently (<Ref>).

Our main contributions are:

* Introducing a method for using attribution patching on all computational graph edges for automated circuit discovery (Edge Attribution Patching, <Ref>).

* Benchmarking Edge Attribution Patching vs existing circuit discovery methods (<Ref>).

* Finding and explaining some limitations with Edge Attribution Patching (<Ref>).

§ RELATED WORK

Automated Circuit Discovery refers to finding the important subgraph of models' computational graphs for performance on particular tasks <cit.>. Existing algorithms include efficient heuristics <cit.> and gradient-descent based methods <cit.>.

Activation Patching is a technique for analyzing the role of individual components in a model. It involves targeted manipulations of activations during a forward pass (further explained in <Ref>). Previous works applied this technique under various names, such as Interchange Interventions <cit.>, Causal Mediation Analysis <cit.> and Causal Tracing <cit.>. We adapt the terminology used by Conmy et al. <cit.>.

Transformer Circuits. Our work builds upon the framework for understanding transformers for interpretability as introduced by <cit.>. The important details include how they formulate forward passes of transformer models. Individual attention heads and MLPs (collectively called nodes) read and write information to a central communication channel, also called the residual stream. In these terms we can examine dependencies of nodes with the output of earlier nodes, i.e we can measure the effect of attention heads in layer 0 on the attention heads in layer 2. In the following, we view these dependencies as edges between nodes, building on existing work using this perspective <cit.>.

§ EDGE ATTRIBUTION PATCHING

We present Edge Attribution Patching (EAP) as a technique to identify relevant model components for solving a specific task. In the following, we view language models as directed, acyclic graphs. In these terms, we aim to find small subgraphs that retain good performance on narrow tasks. We determine the importance of a specific edge through targeted manipulation of activations during a forward pass. We compare two approaches, Attribution Patching and Activation Patching, in order to motivate EAP.

§.§ Activation Patching

Activation patching refers to replacing the activations from one model forward pass with the activations from a different forward pass. This method is typically applied to measure the counterfactual importance of model components, i.e. to measure a statistic $L(x)$ from model outputs under the activation patching, where $x$ is an input prompt. For example, $L$ often represents loss or logit difference <cit.>.

Following existing work (<Ref>), we study the effect of activation patching on specific model edges by setting these equal to activations from different forward passes. Concretely, suppose that an edge $E$ in the computational graph has activation $e_\text{corr}$ on some corrupted prompt. In this work, we use the change in metric under activation patching

\begin{equation}
|L(\xclean |\ \text{do}(E = e_\text{corr})) - L(\xclean )|
\label{eqn:activation_patching}
\end{equation}

to measure the impact of edge $E$. We use Do-notation from causality <cit.> to emphasise that activation patching is a causal intervention.

§.§ Attribution Patching

Activation patching slows ACDC since each measurement (like <Ref>) requires another forward pass. Attribution patching <cit.> is a technique for estimating <Ref> for many different edges $E$ using only two forward passes and one backward pass.[Attribution patching (like activation patching) also applies to nodes and other model internal components that aren't edges, but we only use edges in this work.] It linearly approximates the metric difference after corrupting a single edge in the computational graph (<Ref>) by expanding $L$ as a function of the edge activation as a Taylor series with terms up to the first order:

\begin{equation}
L(\xclean |\ \text{do}(E=e_\text{corr})) \approx L(\xclean) + \underbrace{\left(e_{\text{corr}} - e_{\text{clean}}\right)^\top\frac{\partial}{\partial e_{\text{clean}}} L(\xclean |\ \text{do}(E=e_{\text{clean}}))}_{\text{Call this } \Delta_e L\text{, the \textbf{attribution score}.}}.
\label{eqn:attribution_patching}
\end{equation}

A simple rearrangement implies that <Ref> is approximately equal to eqn:abs_val$\left|\Delta_e L \right|$ eqn:abs_val(3) which we call the absolute attribution score for the rest of this paper.

In practice, all gradients needed to calculate the attribution scores come from intermediate terms computed in one ordinary backwards pass[In <Ref> we show how only one backwards pass is required.] in PyTorch <cit.>, hence attribution patching is extremely efficient.

§.§ Edge Attribution Patching

We can use the insights from <Ref> to build an automated circuit discovery algorithm. This takes two steps: i) use Equation eqn:abs_val(3) to obtain absolute attribution scores for the importance of all edges in the computational graph and then ii) sort these scores and keep the top $k$ edges in a circuit.
We use Edge Attribution Patching (EAP) to refer to this algorithm.
In the rest of the work we report results for all $k$ values when we evaluate EAP (similar to HISP in <cit.>).

axis lines=center,
zlabel=$z$: $L$,

3[surf, opacity=0.7, samples=20, domain=0:4, y domain=0:4] 4*x^2+4*y^2+5;

3[surf, opacity=0.7, color=gray, samples=2, domain=0.5:4.5, y domain=0.5:4.5] 16*x+16*y-27;

3[only marks, mark=*] coordinates (2,2,37) (3.5,3.5,80);

[above] at (axis cs:2.0,1.7,20) \(e_{\text{clean}}\);
[above] at (axis cs:4.0,3.9,57) \(e_{\text{corr}}\);
[above] at (axis cs:4.0,1.0,0) $(x, y)$: Activation;

[->, thick] (axis cs:2.01,2.01,37) – (axis cs:3.4,3.4,80);

Attribution Patching (<Ref>) approximates the difference in metric $L$ caused by corrupting edges.

Removing the least important edges.

Edge Attribution Patching (EAP)

Note that one limitation of attribution patching is that it will not work when the gradient of the metric is the zero vector. <cit.> recommended the use of KL divergence as a metric, which is i) equal to 0 when we run the model without ablations and ii) a non-negative metric. Therefore the zero point is a global minimum and hence all gradients are the zero vector at this point. In this work we use the `task-specific metrics' (not KL divergence) from <cit.> so avoid this issue.

§ RESULTS

§.§ Edge Attribution Patching vs Activation Patching vs ACDC

We compare Edge Attribution Patching (EAP) and ACDC on the Indirect Object Identification (IOI), Docstring, and Greater-Than tasks. For each of these tasks, previous studies identified a subgraph (circuit) relevant for solving the task. We use their results as a ground truth for benchmarking both methods. We also compare using ACDC with the task-specific metrics (e.g logit difference) and KL Divergence (which was originally recommended). For the docstring task, we also include repeated activation patching as another point of reference for performance comparisons. We applied repeated activation patching by running the same circuit discovery method described in <Ref> but using <Ref> rather than absolute attribution scores. Activation patching was not included in the other tasks as it was too computationally expensive to run on the GPT-2 small models used by IOI and Greater-Than. Subnetworks found using EAP for all three tasks are shown in <Ref>.

IOI task

Greater-Than task

Docstring task

ROC Curves comparing EAP, ACDC with task metric, and ACDC with KL Divergence. The Docstring plot also compares to Activation Patching.

The ROC curves in Figure <ref> suggest the performance of EAP is better than ACDC overall: it has the maximal AUC in <Ref>-<ref>, while ACDC used with the KL Divergence metric outperforms EAP in <Ref>.

ACDC outperformed the existing methods HISP and Subnetwork Probing methods <cit.>. We conclude EAP outperforms all previous methods for circuit discovery, since it is competitive with ACDC on recovering circuits while significantly reducing the computational demand: EAP only takes a constant number of forward and backwards passes while the number of forward passes required by ACDC is scaling exponentially with the number of nodes.

§.§ Validating EAP Attribution Scores

In this section, we look at the approximate metric change (attribution score) EAP assigns to each edge in the model. We aim to understand the relation between the attribution score and the function of the edge in the task being studied. First, we look at the distribution of scores for edges in the circuit compared to edges not in the circuit for each of the three tasks.

Distribution of Attribution Scores for the IOI Task (Logit Diff)

Figure <ref> shows the distribution of attribution scores for the IOI task. The distributions for the remaining tasks can be found in <Ref>. Qualitatively, attribution scores for edges in the circuit tend to be spread further from zero. Furthermore, there are only 6 edges outside of the interval $[-0.25, 0.25]$ that aren't part of the IOI circuit.
We further explore the attribution scores for the IOI circuit's classes of heads in <Ref>.

§ LIMITATIONS

We introduced edge activation patching as an approximation to activation patching. However, we found that edge activation patching outperformed ACDC, a technique based on activation patching (<Ref>). In this section, we investigate whether attribution patching's success is due to extremely accurate approximations (in <Ref> we find that the answer is no), and whether there is any further use for ACDC (in <Ref> we find that the answer is yes). We use the docstring task as a case study due to the small model size used.

§.§ How faithful are Attribution Patching's approximations?

To study how faithful the approximation <Ref> is, we plot the attribution patching scores (Equation eqn:abs_val(3)) against the activation patching scores (<Ref>) in <Ref>. Surprisingly, we find a fairly weak correlation between activation and attribution patching scores ($R^2 = 0.27$). Further, the line of best fit has gradient 0.531, suggesting that attribution patching estimates the effect of activation patching as twice as important as it really is.

Moreover, we can gain some sense for the discrepancy between activation and attribution patching by studying the continuous transition between clean ($\eclean$) and corrupted ($\ecorr$) activations in <Ref>, i.e studying the values $|L(\xclean |\ \text{do}(E = \lambda e_\text{corr} + (1 - \lambda) \eclean )) - L(\xclean )|$ for $0 \le \lambda \le 1$. We can compare this to the linear approximations of Attribution Patching $\lambda \Delta_e L$. <Ref> shows the result for one edge in the docstring circuit where the linear approximation to activation patching is not accurate.

Distribution of attribution scores for edges from activation patching and attribution patching. Circled: outlier EAP point studied in <Ref>.

Visualizing the rightmost point in <Ref>. Note that corrupting this edge (surprisingly) slightly increases the logit difference on the Docstring task (higher logit difference is better). However, EAP overestimates how large this increase is.

Visualizing Edge Attribution Patching.

We find that interpolating towards the corrupted input creates a concave curve (<Ref>) such that the linear approximation at $\lambda = 0$ overestimates the effect of activation patching this edge. In <Ref> we show that this also holds for the other outlier edges in the ellipse in <Ref>.

§.§ Is there any further use for ACDC?

In Section <ref> above, we found that EAP overestimates activation patching case the attribution score is concave. This suggests to refine the result by running ACDC on the pruned subgraph returned by EAP. We ran the combined methods on the Docstring task, varying both pruning thresholds for EAP and ACDC independently. Figure <ref> compares the TPR and FPR for the combined methods with the ROC curve of EAP only. The combined methods show increased performance compared to EAP only.

Comparing statistics of the combined EAP + ACDC methods with EAP only. The inset shows a zoom to the significant area of the statistics of the combined method.

Finally, one further limitation of this research is that the metrics used for interpretability do not precisely capture meaningful human understanding. Recovering a subgraph that humans previously recovered is limited because i) we can't evaluate this metric for interpretability tasks that we don't yet understand and ii) human-found circuits are imperfect which increases the noise in this measurement.

§ CONCLUSION

We provide evidence that Edge Attribution Patching (EAP) outperforms ACDC in identifying circuits while being substantially faster to compute. This result is surprising, as EAP is an approximation for activation patching, the method applied by ACDC. However, running ACDC on the prepruned subnetwork found by EAP can improve the identification of relevant edges. Therefore, we suggest future circuit discovery experiments to run EAP first and apply ACDC to the optimal subgraph found by EAP.

§ AUTHOR CONTRIBUTIONS

Aaquib Syed and Can Rager proposed combining ACDC with attribution patching methods and implemented initial prototypes. Arthur Conmy advised working on attributing edges rather than nodes and Aaquib made the first findings that this outperformed Automatic Circuit Discovery. All authors worked on the paper's figures, experiments and code.

§ ACKNOWLEDGEMENTS

We would like to thank Neel Nanda for a helpful discussion and János Kramár for suggestions based on a draft of this work.

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

titleMechanistic Interpretability, Variables, and the Importance of Interpretable Bases

titleZoom In: An Introduction to Circuits

booktitleThe Eleventh International Conference on Learning Representations
titleInterpretability in the Wild: a Circuit for Indirect Object Identification in GPT-2 Small

booktitleThe Eleventh International Conference on Learning Representations
titleProgress measures for grokking via mechanistic interpretability

titleA circuit for Python docstrings in a 4-layer attention-only transformer

titleHow does GPT-2 compute greater-than?: Interpreting mathematical abilities in a pre-trained language model

titleDoes Circuit Analysis Interpretability Scale? Evidence from Multiple Choice Capabilities in Chinchilla

titleCausal Mediation Analysis for Interpreting Neural NLP: The Case of Gender Bias

titleAn overview of 11 proposals for building safe advanced AI

titleTowards Automated Circuit Discovery for Mechanistic Interpretability

booktitleAdvances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada
titleAre Sixteen Heads Really Better than One?

booktitle6th International Conference on Learning Representations, ICLR 2018, Vancouver, BC, Canada, April 30 - May 3, 2018, Conference Track Proceedings
titleLearning Sparse Neural Networks through $L_0$ Regularization

Association for Computational Linguistics

booktitleProceedings of the 2021 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies
titleLow-Complexity Probing via Finding Subnetworks

titleCausal Abstractions of Neural Networks
ine Learning (cs.LG),FOS: Computer and information sciences,FOS: Computer and information sciences

booktitleAdvances in Neural Information Processing Systems
titleLocating and editing factual associations in GPT

journaltitleTransformer Circuits Thread
titleA Mathematical Framework for Transformer Circuits

Oxford University Press

titleCausal diagrams for empirical research

titleAttribution Patching: Activation Patching At Industrial Scale

Curran Associates, Inc.

booktitleAdvances in Neural Information Processing Systems 32
titlePyTorch: An Imperative Style, High-Performance Deep Learning Library

§ EAP SUBNETWORKS

IOI Subnetwork, Threshold=0.077

Docstring Subnetwork, Threshold=0.244

Greaterthan Subnetwork, Threshold=0.009

Resulting subnetworks after EAP at the given thresholds.

§ DISTRIBUTION OF EAP ATTRIBUTION SCORES

Distribution of Attribution Scores for the Docstring Task

Distribution of Attribution Scores for the Greater-Than Task

Distribution of Attribution Scores for the Docstring and Greater-Than tasks

§ FURTHER INVESTIGATION INTO COMBINING EAP WITH ACDC

Youdens-J statistic (maximum TPR minus FPR value) for combining EAP and ACDC methods on the docstring task. We applied ACDC to the pruned subgraph returned by EAP.

§ FURTHER FAILURES OF ATTRIBUTION PATCHING APPROXIMATION

In <Ref> we show further cases where in the docstring task attribution patching can be misleading. These cases all involve an edge that comes from the model's embeddings (positional and tokens). Our interpretation is that weighted averages of embeddings are anomalous inputs to the model and cause the concave change in docstring logit diff which doesn't occur when edges ae between non-embedding model components.

Visualizing Edge Attribution Patching in two further cases where the concave activation patching curve means the linear fit is poor.

§ EDGES ROLES IN IOI

We further explore the attribution scores for the IOI circuit. The IOI circuit is comprised of different attention head classes such as Induction heads, S-Inhibition heads, etc. <cit.>. Figure <ref> shows the distributions of scores stratified by the roles of the edges. The edge roles are defined according to the role of their origin node. While edge roles such as Previous Token, Duplicate Token, Induction, and S-Inhibition edges have attribution scores centered around zero, we see a bias in edge scores given to name mover and negative name mover edges. As the name mover edges are directly responsible for the model outputting the indirect object, the attribution scores are largely negative since ablating these edges removes the model's ability to output the indirect object, lowering the logit difference. Similarly, the negative name movers have attribution scores that are largely positive since ablating these edges improves the logit difference. This matches the intuitive function of the edges.

Distribution of Attribution Scores for each Edge Role in the IOI Task.

§ ONLY ONE BACKWARDS PASS REQUIRED FOR EAP

Note: it may be easier to understand our implementation <https://github.com/Aaquib111/acdcpp/blob/main/utils/prune_utils.py#L249> rather than read this explanation. Alternatively, this derivation uses essentially the same arguments as <cit.>[Specifically, this section: <https://www.neelnanda.io/mechanistic-interpretability/attribution-patching#how-to-think-about-activation-patching=: :text=axes%20of%20variation.-,Path%20patching,-The%20core%20intuition>] though with an updated codebase.

There are only two types of edges iterated over in ACDC: i) residual edges where the result is added at its endpoint, and ii) edges between the residual stream and the query, key and value calculations.[It may be worth looking at some ACDC outputs from <cit.>. See <https://colab.research.google.com/github/ArthurConmy/Automatic-Circuit-Discovery/blob/main/notebooks/colabs/ACDC_Implementation_Demo.ipynb> for an explanation of this design choice.] Clearly for all edges like ii) we can compute the gradient terms in <Ref> in one backwards pass.

Interestingly, for all $\Delta_e L$ terms where $e$ is a type i) edge (i.e added at the endpoint), we only need calculate the gradient with respect to the endpoint of the edge! For example, suppose we're calculating the effect of L0H0 on L1H0Q. If we represent the input to L1H0Q as a node $V$ in the computational graph then

\begin{equation}
\frac{\partial}{\partial e_{\text{clean}}} L(\xclean |\ \text{do}(E=e_{\text{clean}})) = \frac{\partial}{\partial v_{\text{clean}}} L(\xclean |\ \text{do}(V=v_{\text{clean}}))
\end{equation}

due to how $V$ is just the sum of all the edges entering $V$. This allows efficient calculation of all the $\Delta_e L$ values since gradients with respect to nodes in computational graphs are calculated by default in backwards passes.