# Scalar Stochastic Gravitational-Wave Background in Brans-Dicke Theory of
Gravity

Song Ming Du<EMAIL_ADDRESS>TAPIR, Walter Burke Institute for Theoretical
Physics, California Institute of Technology, Pasadena, California 91125, USA

###### Abstract

We study the scalar stochastic gravitational-wave background (SGWB) from
astrophysical sources, including compact binary mergers and stellar collapses,
in the Bras-Dicke theory of gravity. By contrast to tensor waves, we found the
scalar SGWB to predominantly arise from stellar collapses. These collapses not
only take place at higher astrophysical rates, but emit more energy. This is
because, unlike tensor radiation, which mainly starts from quadrupole order,
the scalar perturbation can be excited by changes in the monopole moment. In
particular, in the case of stellar collapse into a neutron star or a black
hole, the monopole radiation, at frequencies below 100 Hz, is dominated by the
memory effect. At low frequencies, the scalar SGWB spectrum follows a power
law of $\Omega_{\text{S}}\propto f^{\alpha}$, with $\alpha=1$. We predict that
$\Omega_{\text{S}}$ is inversely proportional to the square of $\omega_{\rm
BD}+2$, with $\quad(\omega_{\rm BD}+2)^{2}\Omega_{S}(f=25\,{\rm Hz})=2.8\times
10^{-6}$. We also estimate the detectability of the scalar SGWB for current
and third-generation detector networks, and the bound on $\omega_{\rm BD}$
that can be imposed from these observations.

##  I Introduction

The first direct detection of gravitational waves (GWs) from the merger of
binary black holes (BBHs) by the LIGO-Virgo collaboration Abbott  _et al._
(2016a) marks the beginning of gravitational-wave astronomy and opens up a new
window to the Universe. Since then, more GW events, both from BBH mergers and
from binary neutron star (BNS) mergers, are detected by the Advanced
LIGO/Virgo network Abbott  _et al._ (2016b, 2017a, 2017b, 2017c). Besides
these resolvable, individual GW sources we have discovered so far, a
Stochastic Gravitational-Wave Background (SGWB), which arises from the
population of unresolved GW events at larger distances, is anticipated to be
detectable in the upcoming years Abbott  _et al._ (2016c, 2018a).

Gravitational wave signals provide us with unprecedented opportunities to test
general relativity (GR) and study modified theories of gravity Abbott  _et
al._ (2016a, 2017a, d, 2017d). One significant prediction of general
relativity is that gravitational waves only contain two tensor polarization
modes ($+$ mode and $\times$ mode). On the other hand, additional polarization
modes do exist in modified theories of gravity; if directly detected, they
would become a strong evidence for extensions to Einstein’s original theory
Eardley  _et al._ (1973a, b). For example, the Brans-Dicke (BD) theory Brans
and Dicke (1961), which minimally extends Einstein’s gravity by incorporating
a scalar field (Brans-Dicke field) coupled to the metric tensor, predicts the
existence of a transverse scalar polarization mode (also referred to as the
breathing, or the “$\circ$ mode”). Previously, several detection strategies
for the non-tensorial SGWB have been proposed Maggiore and Nicolis (2000);
Nishizawa  _et al._ (2009); Callister  _et al._ (2017). A recent study Abbott
_et al._ (2018b), based on the method in Callister  _et al._ (2017) and the
data from LIGO’s O1 observing run, has placed the first constraints on the
contributions from non-tensorial polarizations to the SGWB.

All the works so far have assumed general, power-law models for the energy
spectra of the non-tensorial SGWB — without considering its specific origins.
However, in order to theoretically estimate the plausible magnitudes of the
non-tensorial SGWB, and to experimentally make statistical inferences on
parameters of modified gravity models from detector data, it is necessary to
consider the astrophysical origins of the non-tensorial SGWB. Furthermore,
obtaining astrophysicically motivated energy spectra may allow us to more
efficiently search for the non-tensorial SGWB using a more optimal matched
filtering technique Allen and Romano (1999) — than simply assuming a power-law
spectrum.

In this paper, we focus on the SGWB in the transverse scalar mode of the
Brans-Dicke (BD) theory: identifying its astrophysical origins, and obtaining
its energy spectrum (as a function of the BD coupling constant). Candidates
for sources of the SGWB include gravitational stellar collapses and compact
binary mergers. As we will see, the existence of monopole scalar radiations
makes stellar collapses by far the major contributor to this SGWB. This
differs significantly from the tensorial SGWB in GR, which is dominated by BBH
and BNS mergers.

This paper is organized as follows: In Sec. II, we will give an overview of
the scalar GW in Brans-Dicke theory and its relation to scalar SGWB. In Secs.
III and IV, we will calculate the contributions to the scalar SGWB from
compact binary coalescences, including BBH and BNS mergers, as well as the
contribution from gravitational stellar collapses. In Sec. V, we will explore
how the scalar SGWB depends on variations in the underlying population models.
In Sec. VI, we will discuss the detectability and possible constraints on the
BD coupling constant $\omega_{\rm BD}$ from current and future observations.
Finally, in Sec. VII, we will draw conclusions and suggest prospective
research directions.

##  II Scalar GW in Brans-Dicke Theory and Relation to SGWB

In the Brans-Dicke (BD) theory, the Lagrangian density of the gravity sector
in the original conformal frame is given by

$\displaystyle\mathcal{L}_{\text{BD}}=\sqrt{-g}\left[\phi
R-\omega_{\text{BD}}\
\frac{\partial^{\mu}\phi\partial_{\mu}\phi}{\phi}\right],$ (1)

where the Ricci scalar $R$ is associated with the spacetime metric tensor
$g_{\mu\nu}$. The scalar field $\phi$ is related to the gravitational constant
$G$ via the relation

$G\phi_{0}=\frac{2\omega_{\text{BD}}+4}{2\omega_{\text{BD}}+3}$ (2)

where $\phi_{0}$ is the value of the scalar field at null infinity. The matter
sector Lagrangian density remains the same as in GR, which means the scalar
field does not couple to the matter fields directly. When the model parameter
$\omega_{\text{BD}}$ approaches to infinity, Brans-Dicke theory recovers to
GR. In this rest of this section, we shall discuss the polarization content of
GWs in the BD theory, and describe the energy content of the scalar SGWB.
Details can be found in Refs. Brans and Dicke (1961); Maggiore and Nicolis
(2000).

###  II.1 GWs in the BD Theory

To study GWs in Brans-Dicke theory, we perform a perturbation of the metric
tensor and the scalar field around the Minkowski spacetime and the null
infinity value, respectively:

$\displaystyle
g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu},\qquad\phi=\phi_{0}+\delta\phi,$ (3)

where components of the metric tensor in Minkowski spacetime is chosen to be
$\eta_{\mu\nu}=\text{diag}(-1,1,1,1)$. The perturbative Lagrangian contains a
quadratic cross term:
$h_{\mu\nu}(\partial^{\mu}\partial^{\nu}-\eta^{\mu\nu}\partial^{2})\delta\phi$.
To eliminate this term we redefine the following physical degrees of freedom:

$\displaystyle
H_{\mu\nu}=h_{\mu\nu}+\eta_{\mu\nu}\frac{\delta\phi}{\phi_{0}}\,.$ (4)

Under these treatments, the perturbative Lagrangian is expressed as:

$\displaystyle\mathcal{L}_{\text{BD}}=\mathcal{L}_{\text{BD}}^{\text{kin}}+\mathcal{L}_{\text{BD}}^{\text{S}}+\mathcal{L}_{\text{BD}}^{\text{other}}\,,$
$\displaystyle\mathcal{L}_{\text{BD}}^{\text{kin}}=\frac{\phi_{0}}{2}H^{\mu\nu}V_{\mu\nu\rho\sigma}H^{\rho\sigma}+\frac{\omega_{\text{BD}}+3/2}{\phi_{0}}\
\eta^{\mu\nu}\delta\phi\ \partial_{\mu}\partial_{\nu}\delta\phi\,,$
$\displaystyle\mathcal{L}_{\text{BD}}^{\text{S}}=\frac{\omega_{\text{BD}}+3/2}{\phi_{0}}\left(H^{\mu\nu}-\frac{1}{2}\eta^{\mu\nu}H\right)\partial_{\mu}\delta\phi\partial_{\nu}\delta\phi\
.$ (5)

Here $\mathcal{L}_{\text{BD}}^{\text{kin}}$ represents the kinetic terms for
the tensor field $H_{\mu\nu}$ and the scalar field $\delta\phi$, where the
operator $V_{\mu\nu\rho\sigma}$ is defined as:

$\displaystyle V_{\mu\nu\rho\sigma}=\frac{1}{2}$
$\displaystyle\Big{[}(\eta_{\mu\rho}\eta_{\nu\sigma}-\eta_{\mu\nu}\eta_{\rho\sigma})\partial^{2}+\eta_{\mu\nu}\partial_{\rho}\partial_{\sigma}$
$\displaystyle+\eta_{\rho\sigma}\partial_{\mu}\partial_{\nu}-\eta_{\mu\rho}\partial_{\nu}\partial_{\sigma}-\eta_{\nu\sigma}\partial_{\mu}\partial_{\rho}\Big{]}\
.$

The Lagrangian $\mathcal{L}_{\text{BD}}^{\text{S}}$ contains the leading
interaction terms between the scalar and the tensor fields. Later we shall
show that it relates to the scalar stress-energy tensor. The third term
$\mathcal{L}_{\text{BD}}^{\text{other}}$ contains other higher order
interaction terms. Notice that the Lagrangians in Eq. (5) is invariant under
the infinitesimal diffeomorphism transformation $x^{\mu}\rightarrow
x^{\prime\mu}=x^{\mu}+\xi^{\mu}(x)$:

$\displaystyle\begin{cases}&H_{\mu\nu}\rightarrow
H^{\prime}_{\mu\nu}=H_{\mu\nu}-\partial_{\mu}\xi_{\nu}-\partial_{\nu}\xi_{\mu}\\\
\\\ &\delta\phi\rightarrow\delta\phi^{\prime}=\delta\phi\ .\end{cases}$ (6)

The vacuum field equation for $H_{\mu\nu}$ is obtained from
$\delta\mathcal{L}_{\text{BD}}^{\text{kin}}/\delta H_{\mu\nu}=0$, which gives
$V^{\mu\nu\rho\sigma}H_{\rho\sigma}=0$. If we choose the harmonic coordinate
condition, this equation is reduced to

$\displaystyle\begin{cases}\partial^{2}H_{\mu\nu}=0\\\
\partial^{\mu}H_{\mu\nu}-\frac{1}{2}\partial_{\nu}H=0\ .\end{cases}$ (7)

Notice that the vacuum field equation Eq. (7) and the gauge transformation Eq.
(6) have the same form as in GR, hence we can similarly gauge away redundant
degrees of freedom which leave only two physical ones.

However, GW detectors respond directly to the change in the spacetime metric,
i.e. $h_{\mu\nu}$, which depends both on $H_{\mu\nu}$ and $\delta\phi$. As a
result, three physical degrees of freedom remain for $h_{\mu\nu}$ Maggiore and
Nicolis (2000). More specifically, within a spatial slice, let
$\Sigma_{\Omega}$ be the 2-D plane perpendicular to the plain wave propagation
direction $\hat{\mathbb{\Omega}}$, and let $\hat{\mathbb{m}}$,
$\hat{\mathbb{n}}$ be two the orthogonal unit vectors in $\Sigma_{\Omega}$,
then we can find a gauge in which the plain wave can be expanded as,

$\displaystyle
h_{ij}(x)=h_{+}(x)e_{ij}^{+}+h_{\times}(x)e_{ij}^{+}+h_{\text{S}}(x)e_{ij}^{\text{S}}.$
(8)

Here, the amplitudes are related to $H_{ij}$ and $\delta\phi$ via

$\displaystyle h_{+}=e_{+}^{ij}H^{\text{TT}}_{ij}/2,\quad
h_{\times}=e_{\times}^{ij}H^{\text{TT}}_{ij}/2,\quad
h_{\text{S}}=-\delta\phi/\phi_{0},$ (9)

here $H^{\text{TT}}_{ij}$ is the transverse-traceless part of $H_{ij}$ Thorne
(1989). The polarization tensors are expressed as

$\displaystyle e^{+}_{ij}=\hat{m}_{i}\hat{m}_{j}-\hat{n}_{i}\hat{n}_{j},\qquad
e^{\times}_{ij}=\hat{m}_{i}\hat{n}_{j}+\hat{n}_{i}\hat{m}_{j},$ $\displaystyle
e^{S}_{ij}=\hat{m}_{i}\hat{m}_{j}+\hat{n}_{i}\hat{n}_{j}.$ (10)

Here $\mathbb{e}^{+,\times}$ and $\mathbb{e}^{\text{S}}$ represent tensor and
scalar polarizations of GW respectively because under an $\mathrm{SO}(2)$
rotation in $\Sigma_{\Omega}$ plane:
$\hat{\mathbb{m}}^{\prime}+i\hat{\mathbb{n}}^{\prime}=\text{exp}(i\theta)(\hat{\mathbb{m}}+i\hat{\mathbb{n}})$,
they behave as
${\mathbb{e}^{+}}^{\prime}+i{\mathbb{e}^{\times}}^{\prime}=\text{exp}(2i\theta)(\mathbb{e}^{+}+i\mathbb{e}^{\times})$
and ${\mathbb{e}^{S}}^{\prime}=\mathbb{e}^{S}$.

###  II.2 The scalar SGWB

We expect the presence of the scalar GW would give rise to a stochastic
background with scalar polarization, which is described by the dimensionless
energy density spectrum:

$\displaystyle\tilde{\Omega}_{\text{S}}(f)=\frac{1}{\rho_{c}}\frac{d\rho_{\text{S}}}{d\ln{f}}.$
(11)

In this equation, $\rho_{c}=3H_{0}^{2}/8\pi G$ is the critical density to
close the Universe with $H_{0}$ representing the Hubble constant. The energy
density of the scalar GW $\rho_{\text{S}}$ relates to the scalar stress-energy
tensor $T_{\text{S}}^{\mu\nu}$ via $\rho_{\text{S}}=T_{\text{S}}^{00}$, with

$\displaystyle
T_{\text{S}}^{\mu\nu}=\frac{1}{8\pi}\frac{\partial\mathcal{L}_{\text{BD}}^{\text{S}}}{\partial
H_{\mu\nu}}=\frac{\omega_{\text{BD}}+2}{8\pi
G}\left(\partial^{\mu}h_{\text{S}}\partial^{\nu}h_{\text{S}}-\frac{\eta^{\mu\nu}}{2}\partial^{\rho}h_{\text{S}}\partial_{\rho}h_{\text{S}}\right)\
.$ (12)

Combining with the field equation $\partial^{2}h_{\text{S}}=0$ and averaging
over several wave length, we obtain Maggiore and Nicolis (2000)

$\displaystyle\rho_{\text{S}}=\frac{\omega_{\text{BD}}+2}{8\pi
G}\left<\dot{h}_{\text{S}}^{2}(x)\right>$ (13)

Under the assumption that the stochastic background is stationary, isotropic
and Gaussian, the ensemble average of the Fourier transformed amplitude
$\tilde{h}_{\text{S}}(f,\hat{\Omega})$ satisfies

$\displaystyle\left<\tilde{h}^{*}_{\text{S}}(f,\hat{\Omega})\tilde{h}_{\text{S}}(f^{\prime},\hat{\Omega}^{\prime})\right>=\frac{5}{8\pi}\delta(\hat{\Omega}-\hat{\Omega}^{\prime})\delta(f-f^{\prime})H_{\text{S}}(f).$
(14)

Here $\tilde{h}_{\text{S}}(f,\hat{\Omega})$ relates to $h(t,\mathbf{x})$ via

$h_{S}(t,\mathbf{x})=\int\frac{d^{3}\mathbf{k}}{(2\pi)^{3}}e^{-i\omega(\mathbf{k})t+i\mathbf{k}\cdot\mathbf{x}}\tilde{h}_{S}+\mbox{c.c.}$
(15)

with $\omega(\mathbf{k})=|\mathbf{k}|/c$, and

$\mathbf{k}=2\pi f\hat{\Omega}/c$ (16)

where $\hat{\Omega}$ is the unit vector along the direction of $\mathbf{k}$.
The quantity $H_{\text{S}}$ is defined as the spectral density for scalar GW.
The factor of $5/8\pi$ follows the same convention in Abbott  _et al._
(2018b); Callister  _et al._ (2017).

Under this definition, $\tilde{H}_{S}$ is related to the scalar spectral
density $\tilde{\Omega}_{\text{GW}}^{\text{S}}$ [defined in Eq. (11)–(13)],
via

$\displaystyle\tilde{\Omega}_{\text{S}}(f)=(\omega_{BD}+2)\frac{20\pi^{2}}{3H_{0}^{2}}f^{3}H_{\text{S}}(f)\
.$ (17)

As we shall see later in Sec. V, the quantity $H_{S}$ is directly related to
the detectability of the scalar SGWB [see Eq. (33)]. In this way, even though
$\tilde{\Omega}_{S}$ is directly proportional to the energy density of the
scalar wave, detectability of the background, given the same
$\tilde{\Omega}_{S}$, still depends on the BD coupling constant $\omega_{\rm
BD}$ . This is related to the violation of the Isaacson formula in BD theory
Allen and Romano (1999); Isaacson (1968). Instead, following the same
convention as Ref. Abbott  _et al._ (2018b), we define a new quantity

$\displaystyle\Omega_{\text{S}}(f)=\frac{\tilde{\Omega}_{\text{S}}(f)}{\omega_{BD}+2}=\frac{20\pi^{2}}{3H_{0}^{2}}f^{3}H_{\text{S}}(f)\,,$
(18)

In the following discussions, we will keep using this redefined energy density
spectrum to describe the scalar SGWB.

##  III Scalar and Tensor SGWB from Mergers of Compact Binary System

###  III.1 Tensor SGWB from Compact Binary Mergers in BD Theory

In GR, the SGWB has only tensor polarization and the major contribution within
the bandwidth of ground based GW detectors is from the mergers of BBH, with
$\Omega_{\text{T}}(f=25\text{Hz})\simeq 1.1\times 10^{-9}$ Abbott  _et al._
(2016c). Besides BBH, the mergers of BNS has a comparable contribution to the
SGWB, with $\Omega_{\text{T}}(f=25\text{Hz})\simeq 0.7\times 10^{-9}$ Abbott
_et al._ (2018a). In BD, we expect the tensor SGWB takes approximately the
same value as in GR, which is predicted from the relation Brunetti  _et al._
(1999); O’Connell and Salmona (1967):

$\displaystyle P_{\text{T}}^{\text{(BD)}}=\frac{2\omega_{\rm
BD}+3}{2\omega_{\rm BD}+4}P_{\text{T}}^{\text{(GR)}},$ (19)

where $P_{\text{T}}^{\text{(BD)}}$ and $P_{\text{T}}^{\text{(GR)}}$ denote the
power emitted in GW with tensor polarization from a system of binary stars in
BD and in GR respectively, at the same orbital frequency. For a large
$\omega_{\text{BD}}$, we expect the ratio between the two powers is
approximately equal to one. As will be shown later in the next section, in BD
most of gravitational radiation by binary stars is from tensor GW — with
scalar radiation suppressed by $\omega_{\rm BD}$. Consequently, the coalescing
trajectory of the compact binary system which is mainly a result from GW
radiation, as well as the spectrum of tensor GW radiation, is nearly unchanged
as in GR.

###  III.2 Scalar Radiation from a Compact Binary in BD Theory

As for the scalar part, the story is completely different: the contribution to
the scalar SGWB from the mergers of BBH is exactly zero. This is a direct
implication from Hawking’s no scalar-hair theorem of black holes in BD theory
of gravity Hawking (1972). The theorem states that for black holes in BD the
exterior spacetine geometry is the same as in GR and the scalar field $\phi$
takes a constant value. Since $h_{\text{S}}=0$ everywhere, there is no scalar
GW radiation from the merger of BBH.

On the other hand, the no scalar-hair theorem does not forbid scalar GWs
emitted from mergers of BNS. Within the bandwidth of ground based GW
detectors, the background is mainly from the inspiral stage, since BNS merger
frequency is above 2 kHz Abbott  _et al._ (2018a). The power of scalar GW
emission from inspiraling binary systems has been studied in Brunetti  _et
al._ (1999); O’Connell and Salmona (1967). Contrary to the tensor case, the
scalar GW has monopole and dipole radiations in addition to quadrupole
radiation. In the limit of vanishing eccentricity $e\rightarrow 0$ (this
assumption should be valid since the orbital angular momentum should have been
radiated away from GW emission for coalescing binary systems as they enter the
band of ground-based detectors) the scalar energy spectrum for the monopole
radiation ($j=0$) is negligible (binary systems with circular orbit have no
monopole moment), while the dipole ($j=1$) and quadruple ($j=2$) are given by

$\displaystyle\frac{dE_{\text{S}}^{j=1}}{df}=\frac{1}{\omega_{\text{BD}}+2}\frac{5}{48}\left(\frac{\text{BE}_{1}}{m_{1}}-\frac{\text{BE}_{2}}{m_{2}}\right)\frac{m_{1}m_{2}}{m_{1}+m_{2}}f^{-1}$
$\displaystyle\frac{dE_{\text{S}}^{j=2}}{df}(f)=\frac{1}{\omega_{\text{BD}}+2}\frac{(\pi
G)^{\frac{2}{3}}}{36}\frac{m_{1}m_{2}}{(m_{1}+m_{2})^{\frac{1}{3}}}f^{-\frac{1}{3}}\
.$ (20)

Here $f$ represents the frequency of GW, $m_{1}$ and $m_{2}$ the masses of the
two neutron stars in the binary system. The energy spectrum is derived from
the relation to the power: $dE_{\text{S}}^{j}/df=P_{\text{S}}^{j}/\dot{f}$,
where we adopt the power of scalar GW emission $P_{\text{S}}^{j}$ calculated
by Brunetti et al. Brunetti  _et al._ (1999). In the limit of $e\rightarrow
0$, the orbital frequency $F$ and the GW frequency $f$ are related by $f=jF$
for $j=1,2$. The rate of change of the orbital frequency due to GW emission is
the same as in GR Misner  _et al._ (2017):

$\displaystyle\dot{F}=\frac{48\pi^{\frac{8}{3}}G^{\frac{5}{3}}}{5}\frac{m_{1}m_{2}}{(m_{1}+m_{2})^{\frac{1}{3}}}(2F)^{\frac{11}{3}}$
(21)

In Eq. (20), BE represents the binding energy of neutron stars and we adopt
the model by Lattimer and Prakash Lattimer and Prakash (2001), which reads

$\displaystyle\frac{\text{BE}}{m}\simeq\frac{0.6\ \beta}{1-0.5\ \beta},$ (22)

where $\beta=Gm/R$ with $R$ denoting the radius of the neutron star.

###  III.3 Scalar SGWB from Compact Binaries in BD Theory

The energy density spectrum of the produced SGWB can be obtained from the
emission spectrum of a single BNS merger event via Phinney (2001); Abbott  _et
al._ (2016c)

$\displaystyle\Omega_{\text{S}}^{j}(f)=\frac{1}{\omega_{\text{BD}}+2}\frac{f}{\rho_{c}}\int_{0}^{z_{\text{max}}}dz\
\frac{R_{m}(z)\ \frac{dE^{j}_{\text{S}}}{df}(f_{z})}{(1+z)H(z)}\ ,$ (23)

where $f_{z}=(1+z)f$ is the frequency at emission. Note that the factor of
$1/(\omega_{\text{BD}}+2)$ is from the definition of Eq. (18). Here we adopt
the $\Lambda\text{CDM}$ cosmological model, with

$H(z)=H_{0}[\Omega_{M}(1+z)^{3}+\Omega_{\Lambda}]^{1/2},$ (24)

where the Hubble constant $H_{0}=70\text{km/s Mpc}$, $\Omega_{M}=0.3$ and
$\Omega_{\Lambda}=0.7$. The redshift cutoff is chosen as $z_{\text{max}}=10$.
In Eq. (23), $R_{m}(z)$ is the BNS merger rate per comoving volume at redshift
$z$. We adopt the same merger rate as in Abbott  _et al._ (2018a), which is
expressed as

$\displaystyle
R_{m}(z)=R_{m}(0)\frac{\int_{t_{\text{min}}}^{t_{\text{max}}}dt_{d}\
R_{f}[z_{f}(z,t_{d})]p(t_{d})}{\int_{t_{\text{min}}}^{t_{\text{max}}}dt_{d}\
R_{f}[z_{f}(0,t_{d})]p(t_{d})}\ .$ (25)

Here, $t_{d}$ denotes the time delay between formation and merger of BNS and
$p(t_{d})$ is its probability distribution function. We assume
$p(t_{d})\propto 1/t_{d}$ for $t_{\text{min}}<t_{d}<t_{\text{max}}$, with
$t_{\text{min}}=20\ \text{Mpc}$ and $t_{\text{max}}$ equal to the Hubble time
$H_{0}^{-1}$. The BNS formation rate $R_{f}(z)$ is assumed to be proportional
to the star formation rate (SFR): $R_{f}(z)\propto\dot{\rho}_{*}(z)$. As in
Abbott  _et al._ (2016c, 2018a) we adopt the GRB-based SFR model given in
Vangioni  _et al._ (2015), which is inferred from observed gamma-ray burst
data at high redshift Kistler  _et al._ (2013). The local BNS merger rate is
inferred from GW170817 Abbott  _et al._ (2017c) with $R_{m}(0)=1540\
\text{Gpc}^{-3}\text{yr}^{-1}$, and $z_{f}(z,t_{d})$ is the redshift at the
binary formation time $t_{f}=t(z)-t_{d}$, with $t(z)$ the age of the Universe
at merger.

In Eq. (23), the energy spectrum is given by Eq. (20) with the observed GW
frequency $f$ replaced by the frequency at emission $f_{z}$. The frequency
cutoff is at the innermost stable circular orbit (ISCO) Zhu  _et al._ (2013):
$f_{\text{max}}=f_{\text{ISCO}}\simeq 4400/(m_{1}+m_{2})\ \text{Hz}$, with the
mass in the unit of $M_{\odot}$. As in Abbott  _et al._ (2018a), the neutron
star masses $m_{1}$ and $m_{2}$ in the binary are assumed to follow a uniform
distribution ranging from $1$ to $2\ M_{\odot}$. We adopt the neutron star
mass-radius relation from the baseline model of Steiner et al. Steiner  _et
al._ (2013). Within our range of $m$, the radius is around $R\approx 12\
\text{km}$.

Figure 1: The scalar SGWB from compact binary systems. The yellow curve is the
contribution from mergers of BNS with $j=1$, at low frequencies it follows a
power law of $f^{0}$. The blue curve is the contribution from mergers of BNS
with $j=2$, at low frequencies it follows a power law of $f^{2/3}$. The green
curve is the total BNS scalar SGWB. The mergers of BBH has no contribution to
the scalar SGWB, which is a direct consequence of Hawking’s no scalar-hair
theorem Hawking (1972).

We show the resulting scalar SGWB in Fig 1. Note that the energy density
spectrum $\Omega_{\rm S}$ we have chosen to use scales with the BD parameter
as $\Omega_{\text{S}}\propto(\omega_{\text{BD}}+2)^{-2}$. For BNS mergers, we
predict $(\omega_{\text{BD}}+2)^{2}\
\Omega_{\text{S}}^{j=1}(f=25\text{Hz})=1.1\times 10^{-11}$ with a power law of
$f^{0}$ at low frequencies and $(\omega_{\text{BD}}+2)^{2}\
\Omega_{\text{S}}^{j=2}(f=25\text{Hz})=6.0\times 10^{-11}$ with a power law of
$f^{2/3}$ at low frequencies. For $f>10\text{Hz}$, the dipole $(j=1)$
contribution to the scalar SGWB is much less than the quadrupole $(j=2)$,
which is a consequence from the small asymmetry between the two neutron stars
in the binary system. Also, as discussed earlier BBH has no contribution to
the scalar SGWB.

##  IV Scalar SGWB from Stellar Gravitational Core Collapse

It is well known that massive stars end their lives through gravitational core
collapse. In GR, stellar core collapses only contribute a minor fraction to
the total SGWB. For example, Crocker et al. Crocker  _et al._ (2015) predict
an SGWB from the black hole ringdown following the collapse with
$\Omega_{\text{T}}(f=25\text{Hz})\simeq 2\times 10^{-12}$, (Fig. 6 of Crocker
_et al._ (2015), model 2 & 3). In Buonanno  _et al._ (2005) Buonanno et al.
predict the background from the neutrino burst associated with the core
collapse with $\Omega_{\text{T}}(f=25\text{Hz})\simeq 1\times 10^{-13}$, (Fig.
1 of Buonanno  _et al._ (2005), the optimistic model). The small contribution
to the SGWB given the greater event rate of stellar collapses compared to
binary mergers can be explained by the fact that in GR, the tensor GWs are
emitted through secondary effects of stellar collapse: only the small
asymmetry in the collapse gives rise to a non-zero quadrupole moment.

###  IV.1 Scalar Emission from Gravitational Core Collapse

However, we expect a different picture in BD: the scalar GW emission starts
from the monopole order, which indicates even the perfectly spherical
collapses are able to emit scalar GW. Further, the progenitors of collapse are
sources of the scalar field, with a monopole moment proportional to
$m/(2\omega_{\text{BD}}+4)$ Shibata  _et al._ (1994). As the progenitor
collapses, this scalar monopole is radiated away. In this way, the scalar GW
is dominated by the memory effect at low frequencies Du and Nishizawa (2016):
the scalar GW slumps from a nonzero initial value $h_{\text{S}}^{\text{ini}}$
to a zero final value if the collapse remnant is a black hole or a different
nonzero final value if the remnant is a neutron star. The change in the
amplitude of the scalar field is expressed as Du and Nishizawa (2016):

$\displaystyle\Delta
h_{ij}=\begin{cases}\displaystyle-\frac{1}{\omega_{\text{BD}}+2}\frac{G\
m}{r}\ e^{\text{S}}_{ij}&\quad\text{black hole remnant}\\\ \\\
\displaystyle-\frac{1}{\omega_{\text{BD}}+2}\frac{G(m-m_{\text{NS}})}{r}\
e^{\text{S}}_{ij}&\quad\text{neutron star remnant}\end{cases}.$ (26)

Here, $m$ and $m_{\text{NS}}$ represent the mass of the progenitor and the
mass of the remnant neutron star, respectively.

As discussed in Du and Nishizawa (2016), for ground-based GW detectors, most
of the detection band is dominated over by the memory as the “zero-frequency
limit”. The resulting scalar energy spectrum from the memory effect is

$\displaystyle\frac{dE_{\text{S}}}{df}=\frac{G\left[m-m_{\text{NS}}\Theta(M_{\text{BH}}-m)\right]^{2}}{\omega_{\text{BD}}+2}\Theta(m-M_{C})\Theta(f_{\text{cut}}-f),$
(27)

where $M_{C}$ is the minimum mass for the progenitor to end its life via core
collapse and $M_{\text{BH}}$ is the mass threshold above which the final
product from collapse is a black hole instead of a neutron star. As suggested
in Smartt (2009), we choose $M_{C}=8M_{\odot}$ and
$M_{\text{BH}}=25M_{\odot}$. The NS mass is chosen as
$M_{\text{NS}}=1.4M_{\odot}$. The cutoff frequency of the memory effect is
$f_{\text{cut}}\simeq 1/\tau_{c}$, where the collapsing time is approximated
from the Oppenheimer-Snyder model Oppenheimer and Snyder (1939):

$\displaystyle\tau_{c}\simeq Gm\frac{\pi}{\sqrt{8\beta^{3}(1-2\beta)}},$ (28)

where $\beta$ is the same as in Eq. (22). Here we choose $\beta=0.1$ for the
progenitor as in Du and Nishizawa (2016); Shibata  _et al._ (1994).

###  IV.2 Scalar SGWB from Core Collapse

From the individual energy spectrum, the total scalar SGWB energy density
spectrum can be obtained using knowledge of collapse rates throughout the age
of the universe Crocker  _et al._ (2015),

$\displaystyle\Omega_{\text{S}}(f)=\frac{1}{\omega_{\text{BD}}+2}\frac{f}{\rho_{c}}\int_{0}^{z_{\text{max}}}dz\int_{M_{C}}^{M_{\text{max}}}dm\
\frac{\frac{dR_{c}}{dm}(z,m)\ \frac{dE_{\text{S}}}{df}(f_{z},m)}{(1+z)H(z)}\
.$ (29)

In this equation, $M_{\text{max}}$ is the upper limit of the mass of massive
stars and here we choose $M_{\text{max}}=100M_{\odot}$ as reference. The
energy spectrum is from Eq. (27) and the other parameters are the same as in
Eq. (23). The collapse rate density $dR_{c}/dm$ (the number of collapses per
unit proper time, per unit co-moving volume and per progenitor mass) is
related to the Star Formation Rate (SFR) and initial mass function $\xi$ via
Ferrari  _et al._ (1999); Zhu  _et al._ (2010)

$\displaystyle\frac{dR_{c}}{dm}(z,m)=\frac{\dot{\rho}_{*}(z)\xi(m)}{\int_{M_{\text{min}}}^{M_{\text{max}}}d\mu\
\mu\ \xi(\mu)}.$ (30)

Here we use the same SFR as in Section III, and choose the Salpeter IMF:
$\xi(m)\propto m^{-2.35}$, with $M_{\text{min}}=0.1M_{\odot}$ and
$M_{\text{max}}=100M_{\odot}$ Vangioni  _et al._ (2015). The total merger
$R_{c}(z)$ rate between $M_{C}$ and $M_{\text{max}}$ together with the BNS
merger rate $R_{m}(z)$ are shown in Fig 2.

Figure 2: Comparison between the BNS merger rate $R_{m}(z)$ and the core
collapse rate $R_{c}(z)$.

In Fig. 3, we show the resulting scalar SGWB from core collapse. Same as the
BNS scalar SGWB, the predicted energy density spectrum scales with the BD
parameter as $\Omega_{\text{S}}\propto(\omega_{\text{BD}}+2)^{-2}$. At the
reference frequency, $(\omega_{\text{BD}}+2)^{2}\
\Omega_{\text{S}}(f=25\text{Hz})=2.8\times 10^{-6}$. At frequencies below
$\sim 40\,$Hz, $\Omega_{S}$ follows a power law of $f^{\alpha}$, with
$\alpha=1$. Note that the core collapse scalar SGWB is around four orders of
magnitude greater than BNS mergers. This difference can be accounted for by
considering two factors. First, the collapse rate is much larger than the
merger rate: at their peak values $R_{c}\simeq 10^{6}\
\text{Gpc}^{-}3\text{yr}^{-1}$ and $R_{m}\simeq 5\times 10^{3}\
\text{Gpc}^{-}3\text{yr}^{-1}$. Second, the energy emitted from scalar GW
radiation for a single collapse event is much larger than a merger event:
notice that $E_{\text{S}}\propto m^{2}$, for mergers $m\sim 1M_{\odot}$ and
for collapses $m\sim 10M_{\odot}$.

Figure 3: Red curve: The scalar SGWB from core collapse. The model parameters
are given in Section IV. Green curve: The scalar SGWB from BNS merger, which
is a sum of $j=1$ and $j=2$ radiation. The model parameters are given in
Section III.

In Table 1 we summarize the energy densities of SGWB with scalar polarization
from varied sources, compared with the tensor SGWB.

| $\Omega_{T}(f=25\,{\rm Hz})\quad$ | $\quad(\omega_{\rm BD}+2)^{2}\Omega_{S}(f=25\,{\rm Hz})$
---|---|---
BBH |  $1.1\times 10^{-9}$ Abbott  _et al._ (2016c) | 0
BNS |  $0.7\times 10^{-9}$ Abbott  _et al._ (2018a) | $7.1\times 10^{-11}$
Collapse |  $2\times 10^{-12}$ Crocker  _et al._ (2015) | $2.8\times 10^{-6}$
Table 1: Energy density of tensor and scalar SGWB at 25 Hz, from various
origins.

###  IV.3 Model dependence of the Core Collapse Scalar SGWB

In this section we want to explore the influence to the core collapse scalar
SGWB from alternative models. From now on we refer the model described in
Section IV as the Baseline model. More specifically, we consider four
alternative models that follow.

(i) The TimeDelay model. In this model we take into account the time delay
between the formation of a massive star and its core collapse. In this case,
the collapse rate is modified as

$\displaystyle\frac{dR_{c}}{dm}(z,m)=\frac{\int_{t_{\text{min}}}^{t_{\text{max}}}dt_{d}\
\dot{\rho}_{*}[z_{f}(z,t_{d})]\xi(m)p(t_{d})}{\int_{M_{\text{min}}}^{M_{\text{max}}}d\mu\
\mu\ \xi(\mu)}.$ (31)

With the other parameters the same as in Eq. (30), we assume the distribution
as $p(t_{d})=\delta(t_{d}-T(m))$, with $T(m)$ the lifetime of a star with mass
$m$. In addition, we use the relation $T(m)=T_{\odot}(m/M_{\odot})^{-2.5}$ for
main sequence stars, with $M_{\odot}$ representing the solar mass and
$T_{\odot}=10^{4}\ \text{Myr}$.

(ii) The AltSFR model. In the baseline model we adopt an SFR model which is
based on the GRB rate. In the AltSFR model we consider an alternative SFR
model Vangioni  _et al._ (2015), which based on the luminosity of star-forming
galaxies Behroozi  _et al._ (2013). This model is more conservative than the
GRB-based SFR at high redshifts. We compare the two SFR models in Fig. 4.

Figure 4: Comparison between the GRB-based SFR and the Luminosity-based SFR.

(iii) The BHonly model. In this model we only consider the scalar SGWB from
core collapses into BHs. In this case the scalar energy spectrum is given by

$\displaystyle\frac{dE_{\text{S}}}{df}=\frac{G\
m^{2}}{\omega_{\text{BD}}+2}\Theta(m-M_{\text{BH}})\Theta(f_{\text{cut}}-f),$
(32)

where the BH mass threshold $M_{\text{BH}}$ and the cutoff frequency
$f_{\text{cut}}$ are the same as in the Baseline model.

(iv) The HighMass model. To reflect the recent observations of massive stars
with $M\sim 200-300M_{\odot}$ Crowther  _et al._ (2010), we replace the mass
upper limit $M_{\text{max}}$ to $200M_{\odot}$, with other parameters
remaining the same.

Figure 5: The scalar overlap reduction function $\gamma_{\text{S}}$ for
LIGO/Voyager and Einstein Telescope.

We show the scalar SGWB from the alternative models in Fig. 5. We can see that
TimeDelay and AltSFR have negligible influence on the background inside the
detection band of ground-based detectors. The BHonly and HighMass alter the
background in low and high frequencies respectively. At the reference
frequency $f=25$ Hz, HighMass predicts a value of $\Omega_{\text{S}}$ that is
$1.1$ times the Baseline value, while BHonly predicts 0.7 the baseline value.
At this frequency, the impact from TimeDelay and AltSFR to the scalar SGWB
spectra is less than $1\%$. At frequencies below $\sim 10\,$Hz, the HighMass
model predicts somewhat higher $\Omega_{S}$, due to contributions from
collapses of higher-mass objects. At frequency $f>100~{}\text{Hz}$, the only
non-negligible change to the spectrum is from the BHonly model. This is
because the stars which collapse into NSs have lower mass compared to those
collapse into BHs. The cut-off frequency in Eq. 32 is related to the
collapsing time in Eq. 28 which is shorter for collapsing stars with lower
mass. As a result, the high frequency part of the spectrum is suppressed from
the missing low mass progenitors.

##  V Detectability

Since the dominant contribution to the scalar SGWB in BD theory is from the
core collapses, in this section we will focus on the scalar background
predicted by the baseline core collapse model as described in Section IV.

A resent analysis Abbott  _et al._ (2018b) based on Advanced LIGO’s first
observing run (O1) has put the first upper limit on the scalar SGWB, with
$\Omega_{\text{S}}(f=25\text{Hz})<1.1\times 10^{-7}$. Compared with our
prediction $(\omega_{\text{BD}}+2)^{2}\
\Omega_{\text{S}}(f=25\text{Hz})=2.8\times 10^{-6}$, it is straightforward to
obtain $\omega_{\text{BD}}>3$. Much better upper limits are expected since the
O1 data only includes an observation time of four months and the detectors are
running below the designed sensitivity.

Next we want to explore the detectability from Advanced LIGO at its designed
sensitivity and the planed future ground based GW detectors. The optimal
signal-to-noise ratio (SNR) for the scalar SGWB between a pair of detectors is
given by Nishizawa  _et al._ (2009); Callister  _et al._ (2017)

$\displaystyle\text{SNR}=\frac{3H_{0}^{2}}{10\pi^{2}}\sqrt{2T}\left(\int_{0}^{\infty}df\
\frac{\gamma_{\text{S}}(f)^{2}\Omega_{\text{S}}(f)^{2}}{f^{6}P_{1}(f)P_{2}(f)}\right)^{1/2},$
(33)

where $P_{1,2}(f)$ are the detectors’ noise spectral density and
$\gamma_{\text{S}}(f)$ is the scalar overlap reduction function between the
detectors Maggiore and Nicolis (2000). Here we recall that it was the choice
we had made in Eq. (18) for $\Omega_{S}$ that would lead to this expression
for the SNR, which is similar to that for a tensor gravitational wave
background.

Here we consider the design sensitivity of Advanced LIGO Ajith  _et al._
(2008) and the planed sensitivities of LIGO Voyager voy (2015) and Einstein
Telescope (ET) Sathyaprakash and Schutz (2009). The scalar overlap reduction
function between the detectors at Hanford and Livingston is calculated in
Callister  _et al._ (2017); Nishizawa  _et al._ (2009), here we adopt the
normalization convention as Callister  _et al._ (2017). Voyager has the same
overlap reduction function as LIGO. The co-located ET detectors have a
constant $\gamma_{\text{S}}=-1/16$ for $f<1000\text{Hz}$ (see Appendix for
more details). These overlap reduction functions are shown in Fig. 6. Note
that our $\gamma_{\text{S}}$ for LIGO is one half of Callister  _et al._
(2017), which is due to the fact that, as explained in Section II, only the
breathing (traverse) and no longitudinal part of the scalar polarization exist
in BD theory.

Figure 6: The scalar overlap reduction function $\gamma_{\text{S}}$ for
LIGO/Voyager and Einstein Telescope.

We show the maximal detectable $\omega_{\text{BD}}$ for LIGO, Voyager and ET
to reach an SNR threshold of 3 in Table 2 with observation times of 1 year and
5 years — and in Fig. 7 for a range of observation times. With 5-year
integration, to reach $\text{SNR}>3$ at ET, the BD parameter should be no less
than $264$. On the other hand, the current cosmological constraints on BD set
$\omega_{\text{BD}}>692$ Avilez and Skordis (2014) and the solar system data
from the Cassini mission put a stronger constraint that
$\omega_{\text{BD}}>40000$ Bertotti  _et al._ (2003); Will (2006).

$T$ | LIGO | Voyager | ET
---|---|---|---
1 yr | $10.8$ | $54.1$ | $175.8$
5 yrs | $17.1$ | $81.8$ | $263.8$
Table 2: Maximal detectable BD parameter $\omega_{\text{BD}}$ to reach an SNR
threshold of $3$ from the scalar SGWB with observation times of 1 year and 5
years. Figure 7: Maximal detectable BD parameter $\omega_{\text{BD}}$ to reach
an SNR threshold of $3$ from the scalar SGWB as a function of observation of
time.

##  VI Conclusions and Discussions

In this paper, we studied the scalar SGWB in BD theory, from astrophysical
sources, in particular compact binary mergers and stellar collapses. Unlike
the tensor SGWB in GR, we found that the scalar SGWB in the BD theory is
dominated by stellar collapses, by roughly 4 orders of magnitude, over compact
binary mergers. We have attributed this dominance to the higher rate of
gravitational collapses than binary mergers, as well as the fact that scalar
radiation does not require asymmetry.

Furthermore, scalar radiation from stellar collapses, in the LIGO band, is
mainly dominated by the memory wave — as pointed out in an earlier paper Du
and Nishizawa (2016). Since the memory wave has a simple frequency dependence
of $h(f)\sim 1/f$, this has lead to a $\Omega_{S}(f)\propto f$, which differs
from the tenor SGWB, which as $\Omega_{\rm T}\propto f^{2/3}$.

For the dominant stellar-collapse scalar SGWB, we have studied a range of
models, which had lead to consistent predictions, with the most significant
uncertainty lying at low frequencies: up to within 30% at $f=25~{}\text{Hz}$,
mainly due to possible existence of heavier stars and the exclusion of the
collapses whose remnant are NSs.

Upon obtaining the SGWB spectrum, we have estimated the detectability for
current and future detector networks. It is estimated that 3rd-generation
ground-based detectors can pose upper limit for $\omega_{\rm BD}$ around $\sim
300$.

The potential bound for $\omega_{\rm BD}$ from our calculation is low compare
with solar-system bounds, and somewhat lower than cosmological bounds, this
nevertheless provides an independent test. More importantly, having
established that the scalar SGWB mainly arise from stellar collapses, we can
further investigate other models that lead to scalar radiations, e.g., scalar-
tensor theories in which $\omega(\phi)$ depends on the value of $\phi$ instead
of being a constant. As we had shown in Ref. Du and Nishizawa (2016), in such
models the scalar memory, which dominates scalar radiation during collapse,
can be significantly enhanced by such dependencies through scalariation Damour
and Esposito-Farese (1993), therefore might lead to much stronger SGWB
enhanced by several orders of magnitude Novak (1998); Harada (1998). In that
case, we expect a considerable increase in the detectability from the current
and the next generation of detectors. We leave these for further studies.

##  VII Acknowledgements

The author thanks Yanbei Chen for valuable discussions and the comments on the
manuscript. The author is also grateful to Xi-Long Fan, Xiang-Cheng Ma and
Atsushi Nishizawa for discussions. We acknowledge support from the Brinson
Foundation and the Simons Foundation. Our research has also been supported by
the National Science Foundation, through grants PHY-1404569, PHY-1708212 and
PHY-1708213.

##  VIII Appendix: Scalar Overlap Reduction Function For Einstein Telescope

Figure 8: Configuration of Einstein Telescope.

In this appendix we calculate the scalar overlap reduction function
$\gamma_{\text{S}}(f)$ for Einstein Telescope (ET). The configuration of ET is
shown in Fig. 8. The coordinate system for the detectors are

$\displaystyle\begin{cases}\hat{\mathbb{x}}=(1,0,0)\\\
\hat{\mathbb{y}}=(0,1,0)\\\ \hat{\mathbb{z}}=(0,0,1)\end{cases}$ (34)

In this coordinate system, the unit vectors for the ET detector arms are

$\displaystyle\hat{\mathbb{l}}_{11}=\hat{\mathbb{x}},\qquad\hat{\mathbb{l}}_{12}=\frac{1}{2}\hat{\mathbb{x}}+\frac{\sqrt{3}}{2}\hat{\mathbb{y}},$
$\displaystyle\hat{\mathbb{l}}_{21}=-\frac{1}{2}\hat{\mathbb{x}}+\frac{\sqrt{3}}{2}\hat{\mathbb{y}},\qquad\hat{\mathbb{l}}_{22}=-\hat{\mathbb{l}}_{11},$
$\displaystyle\hat{\mathbb{l}}_{31}=-\hat{\mathbb{l}}_{12},\qquad\hat{\mathbb{l}}_{32}=-\hat{\mathbb{l}}_{21},$
(35)

The detector tensors for ET are express by

$\displaystyle
D^{ij}_{a}=\frac{1}{2}(\hat{l}_{a1}^{i}\hat{l}_{a1}^{j}-\hat{l}_{a2}^{i}\hat{l}_{a2}^{j}),$
(36)

where $a=1,2,3$.

Suppose the GW is propagating along the angle $(\theta,\phi)$, the GW
coordinate system can be constructed as

$\displaystyle\begin{cases}&\hat{\mathbb{m}}=\cos\theta\cos\phi\
\hat{\mathbb{x}}+\cos\theta\sin\phi\ \hat{\mathbb{y}}-\sin\theta\
\hat{\mathbb{z}}\\\ &\hat{\mathbb{n}}=-\sin\phi\ \hat{\mathbb{x}}+\cos\phi\
\hat{\mathbb{y}}\\\ &\hat{\mathbb{\Omega}}=\sin\theta\cos\phi\
\hat{\mathbb{x}}+\sin\theta\sin\phi\ \hat{\mathbb{y}}+\cos\theta\
\hat{\mathbb{z}}\end{cases}$ (37)

Then the angular pattern functions for scalar polarization are Allen and
Romano (1999); Nishizawa  _et al._ (2009)

$\displaystyle
F^{\text{S}}_{a}(\hat{\mathbb{\Omega}})=\sum_{ij}D_{a}^{ij}e_{\text{S}}^{ij},$
(38)

where the scalar polarization tensor $e_{\text{S}}^{ij}$ are given in Eq.
(10). It is straightforward to find that

$\displaystyle
F^{\text{S}}_{1}(\hat{\mathbb{\Omega}})=\frac{1}{8}\sin^{2}{\theta}(-3\cos{2\phi}+\sqrt{3}\sin{2\phi})$
$\displaystyle
F^{\text{S}}_{2}(\hat{\mathbb{\Omega}})=\frac{1}{8}\sin^{2}{\theta}(3\cos{2\phi}+\sqrt{3}\sin{2\phi})$
$\displaystyle
F^{\text{S}}_{3}(\hat{\mathbb{\Omega}})=-\frac{\sqrt{3}}{4}\sin^{2}{\theta}\sin{2\phi}$
(39)

The scalar overlap reduction function is defined as Nishizawa  _et al._
(2009); Callister  _et al._ (2017)

$\displaystyle\gamma_{ab}^{\text{S}}(f)=\frac{5}{8\pi}\int
d\hat{\mathbb{\Omega}}\ e^{2\pi
if\mathbb{\hat{\Omega}}\cdot\Delta\mathbb{x}}F^{\text{S}}_{a}(\hat{\mathbb{\Omega}})F^{\text{S}}_{b}(\hat{\mathbb{\Omega}}).$
(40)

Here we use the same normalization as Callister  _et al._ (2017). For ET, the
separation $|\Delta\mathbb{x}|$ is equal to the arm length $d=10\ \text{km}$.
Hence, for $f<10^{3}\ \text{Hz}$, the exponential function $e^{2\pi
if\mathbb{\hat{\Omega}}\cdot\Delta\mathbb{x}}\simeq 1$. In this case,

$\displaystyle\gamma^{\text{S}}_{12}=\frac{5}{8\pi}\int_{0}^{\pi}d\theta\
\int_{0}^{2\pi}d\phi\left[-\frac{3}{64}(1+2\cos{4\phi})\sin^{5}{\theta}\right]=-\frac{1}{16}.$
(41)

Similarly, we can show $\gamma^{\text{S}}_{23}=\gamma^{\text{S}}_{31}=-1/16$
for $f<10^{3}\ \text{Hz}$.

## References

  * Abbott  _et al._ (2016a) B. P. Abbott _et al._ , Physical review letters 116, 061102 (2016a).
  * Abbott  _et al._ (2016b) B. P. Abbott _et al._ , Phys Rev Lett 116, 241103 (2016b).
  * Abbott  _et al._ (2017a) B. P. Abbott _et al._ , Physical Review Letters 118, 221101 (2017a).
  * Abbott  _et al._ (2017b) B. P. Abbott _et al._ , Physical review letters 119, 141101 (2017b).
  * Abbott  _et al._ (2017c) B. P. Abbott _et al._ , Physical Review Letters 119, 161101 (2017c).
  * Abbott  _et al._ (2016c) B. Abbott _et al._ , Physical review letters 116, 131102 (2016c).
  * Abbott  _et al._ (2018a) B. P. Abbott _et al._ , Physical review letters 120, 091101 (2018a).
  * Abbott  _et al._ (2016d) B. P. Abbott _et al._ , Physical Review Letters 116, 221101 (2016d).
  * Abbott  _et al._ (2017d) B. P. Abbott _et al._ , The Astrophysical Journal Letters 848, L13 (2017d).
  * Eardley  _et al._ (1973a) D. M. Eardley, D. L. Lee, A. P. Lightman, R. V. Wagoner, and C. M. Will, Physical Review Letters 30, 884 (1973a).
  * Eardley  _et al._ (1973b) D. M. Eardley, D. L. Lee, and A. P. Lightman, Physical Review D 8, 3308 (1973b).
  * Brans and Dicke (1961) C. Brans and R. H. Dicke, Physical review 124, 925 (1961).
  * Maggiore and Nicolis (2000) M. Maggiore and A. Nicolis, Physical Review D 62, 024004 (2000).
  * Nishizawa  _et al._ (2009) A. Nishizawa, A. Taruya, K. Hayama, S. Kawamura, and M.-a. Sakagami, Physical Review D 79, 082002 (2009).
  * Callister  _et al._ (2017) T. Callister, A. S. Biscoveanu, N. Christensen, M. Isi, A. Matas, O. Minazzoli, T. Regimbau, M. Sakellariadou, J. Tasson, and E. Thrane, Physical Review X 7, 041058 (2017).
  * Abbott  _et al._ (2018b) B. Abbott _et al._ , arXiv preprint arXiv:1802.10194 (2018b).
  * Allen and Romano (1999) B. Allen and J. D. Romano, Physical Review D 59, 102001 (1999).
  * Thorne (1989) K. Thorne,  _Three hundred years of gravitation_ (Cambridge University Press, 1989).
  * Isaacson (1968) R. A. Isaacson, Physical Review 166, 1263 (1968).
  * Brunetti  _et al._ (1999) M. Brunetti, E. Coccia, V. Fafone, and F. Fucito, Physical Review D 59, 044027 (1999).
  * O’Connell and Salmona (1967) R. O’Connell and A. Salmona, Physical Review 160, 1108 (1967).
  * Hawking (1972) S. Hawking, Communications in Mathematical Physics 25, 167 (1972).
  * Misner  _et al._ (2017) C. W. Misner, K. S. Thorne, and J. A. Wheeler,  _Gravitation_ (Princeton University Press, 2017).
  * Lattimer and Prakash (2001) J. Lattimer and M. Prakash, The Astrophysical Journal 550, 426 (2001).
  * Phinney (2001) E. Phinney, arXiv preprint astro-ph/0108028 (2001).
  * Vangioni  _et al._ (2015) E. Vangioni, K. A. Olive, T. Prestegard, J. Silk, P. Petitjean, and V. Mandic, Monthly Notices of the Royal Astronomical Society 447, 2575 (2015).
  * Kistler  _et al._ (2013) M. D. Kistler, H. Yuksel, and A. M. Hopkins, arXiv preprint arXiv:1305.1630 (2013).
  * Zhu  _et al._ (2013) X.-J. Zhu, E. J. Howell, D. G. Blair, and Z.-H. Zhu, Monthly Notices of the Royal astronomical Society 431, 882 (2013).
  * Steiner  _et al._ (2013) A. W. Steiner, J. M. Lattimer, and E. F. Brown, The Astrophysical Journal Letters 765, L5 (2013).
  * Crocker  _et al._ (2015) K. Crocker, V. Mandic, T. Regimbau, K. Belczynski, W. Gladysz, K. Olive, T. Prestegard, and E. Vangioni, Physical Review D 92, 063005 (2015).
  * Buonanno  _et al._ (2005) A. Buonanno, G. Sigl, G. G. Raffelt, H.-T. Janka, and E. Müller, Physical Review D 72, 084001 (2005).
  * Shibata  _et al._ (1994) M. Shibata, K. Nakao, and T. Nakamura, Physical Review D 50, 7304 (1994).
  * Du and Nishizawa (2016) S. M. Du and A. Nishizawa, Physical Review D 94, 104063 (2016).
  * Smartt (2009) S. J. Smartt, Annual Review of Astronomy and Astrophysics 47, 63 (2009).
  * Oppenheimer and Snyder (1939) J. R. Oppenheimer and H. Snyder, Physical Review 56, 455 (1939).
  * Ferrari  _et al._ (1999) V. Ferrari, S. Matarrese, and R. Schneider, Monthly Notices of the Royal Astronomical Society 303, 247 (1999).
  * Zhu  _et al._ (2010) X.-J. Zhu, E. Howell, and D. Blair, Monthly Notices of the Royal Astronomical Society: Letters 409, L132 (2010).
  * Behroozi  _et al._ (2013) P. S. Behroozi, R. H. Wechsler, and C. Conroy, The Astrophysical Journal 770, 57 (2013).
  * Crowther  _et al._ (2010) P. Crowther, O. Schurr, R. Mirschi, N. Yusof, R. Parker, S. Goodwin, and H. Kassim, MNRAS 408, 731 (2010).
  * Ajith  _et al._ (2008) P. Ajith, S. Babak, Y. Chen, M. Hewitson, B. Krishnan, A. Sintes, J. T. Whelan, B. Brügmann, P. Diener, N. Dorband, _et al._ , Physical Review D 77, 104017 (2008).
  * voy (2015) LIGO Document T1400316 (2015).
  * Sathyaprakash and Schutz (2009) B. S. Sathyaprakash and B. F. Schutz, Living Reviews in Relativity 12, 2 (2009).
  * Avilez and Skordis (2014) A. Avilez and C. Skordis, Physical review letters 113, 011101 (2014).
  * Bertotti  _et al._ (2003) B. Bertotti, L. Iess, and P. Tortora, Nature 425, 374 (2003).
  * Will (2006) C. M. Will, Living reviews in relativity 9, 3 (2006).
  * Damour and Esposito-Farese (1993) T. Damour and G. Esposito-Farese, Physical Review Letters 70, 2220 (1993).
  * Novak (1998) J. Novak, Physical Review D 58, 064019 (1998).
  * Harada (1998) T. Harada, Physical Review D 57, 4802 (1998).


Research

#  Discovering a change point and piecewise linear structure in a time series
of organoid networks
via the iso-mirror

TCTianyi Chen  YPYoungser Park  ASEAli Saad-Eldin  ZLZachary Lubberts
AAAvanti Athreya  BDPBenjamin D. Pedigo  JTVJoshua T. Vogelstein  FPFrancesca
Puppo  GASGabriel A. Silva  ARMAlysson R. Muotri  WYWeiwei Yang
CMWChristopher M. White  CEPCarey E. Priebe  Department of Applied Mathematics
and Statistics, Johns Hopkins University, Baltimore, USA  Department of
Pediatrics, University of California, San Diego, San Diego, USA  Department of
Neurosciences, University of California, San Diego, San Diego, USA  Microsoft
Research, Seattle, USA

###### Abstract

Recent advancements have been made in the development of cell-based in-vitro
neuronal networks, or organoids. In order to better understand the network
structure of these organoids, a super-selective algorithm has been proposed
for inferring the effective connectivity networks from multi-electrode array
data. In this paper, we apply a novel statistical method called spectral
mirror estimation to the time series of inferred effective connectivity
organoid networks. This method produces a one-dimensional iso-mirror
representation of the dynamics of the time series of the networks which
exhibits a piecewise linear structure. A classical change point algorithm is
then applied to this representation, which successfully detects a change point
coinciding with the neuroscientifically significant time inhibitory neurons
start appearing and the percentage of astrocytes increases dramatically. This
finding demonstrates the potential utility of applying the iso-mirror dynamic
structure discovery method to inferred effective connectivity time series of
organoid networks.

organoid,

time series network analysis,

multidimensional scaling,

manifold learning,

change point detection,

###### keywords:

## Introduction

Detecting structural changes in time series of networks is central to many
modern network science applications. However, due to the complexity of
temporal network data and the myriad possible aspects for potential structural
change, this problem can be daunting. For discovering underlying dynamics in
time series of networks, [1] proposes theory and methods for representing
temporal network structure with a curve, or ‘mirror’, in low dimensional
Euclidean space, enabling the use of classical change point detection
algorithms. In this paper, we estimate the mirror for a time series of brain
organoid connectivity networks and subsequently identify change points.
Because the mirror estimation method requires a 1-1 vertex correspondence for
the networks across time, we first demonstrate that the putative 1-1
correspondence obtained directly from data collection is sufficiently accurate
by comparing it to the vertex correspondence obtained via graph matching [2].
Thence, mirror estimation and manifold learning recovers a 1-dimensional
piecewise linear ‘iso-mirror’ representation, with an evident slope change. By
using the change point detection algorithm from [3] and break point estimation
for piecewise linear models from [4], we identify a change points
neuroscientific significance coinciding with development stages.

We organize this paper as follows. In Section Organoids we introduce our brain
organoids data and the extraction of effective connectivity networks based on
extracellular electrophysiology recordings. In Section Graph Matching we
define the graph matching problem and present the fast approximate assignment
algorithm. In Section Discovering underlying dynamics in times series of
networks we introduce the mirror estimation method and relevant model
assumptions. In the Results section we apply these methods to the time series
of organoid networks and present the graph matching results and the change
point detection results. We conclude the paper with a discussion.

## Organoids

The brain organoids we consider are self-organizing structures composed of
roughly 2.5 million neural cells. They are generated from human induced
pluripotent stem cells (hiPSCs) [5]. After growing for 6 weeks, they are
plated in 8 wells of a multi-electrode array (MEA) plate (Axion Biosystem,
Atlanta,GA, USA). MEA contains 64 low-impedance (0.04 MU) platinum
microelectrodes with 30 µm of diameter spaced by 200 µm. Each well contains
two or three organoids. Then, to characterize the functional development of
the organoids, extracellular spontaneous electrical activity is recorded
weekly using Maestro MEA system and AxIS Software Spontaneous Neural
Configuration (Axion Biosystems) with a 0.1-Hz to 5-kHz band-pass filter. Then
spikes are detected with an adaptive threshold crossing set to 5.5 times the
standard deviation of the estimated noise for each electrode. Each time series
consists of five minutes of recorded neural activity across 10 months and the
data are recorded irregularly – not exactly once a week. Cortical organoids
show low and sparse activity during the first 2 months with an average firing
frequency of 0.5-Hz, then they start exhibiting highly synchronized and
stereotypical network activity which transitions into 2-Hz and 3-Hz rhythmic
activity by 4-6 months. At later stages (6 to 9 months), the activity includes
high-rate spiking with peak of activity reaching a 20-Hz pace and highly
complex bursting behaviors with cross-frequency coupling [6]. A combination of
principal component analysis and k-means clustering is used to spike sort the
multi-unit activity from the 64 electrodes of each well. The average number of
neurons detected in each well increases with the maturation of the organoids
onto the electrodes, sometime reaching saturation after 6 months. Similar
behavior is observed in all MEA wells. For example, the number of spike-sorted
neurons in one well is detected as 122, 160, 189, 174, 190 at 2, 4, 6, 8 and
10 months. The other well has 81, 160, 174, 171, 170 neurons detected at 2, 4,
6, 8, and 10 months.

To infer the effective connections between neurons – i.e., the adjacency
matrices with neurons as vertices across time from the spike activity data –
the algorithm proposed in [6] is applied. This algorithm uses a super-
selection rule to individuate and discard correlation peaks corresponding to
apparent and indirect interactions, and reconstructs the effective
connectivity of the network considering the remaining correlation delays. In
the end, 45 effective connectivity networks on 127 vertices across 244 days
are obtained. This is our time series of organoid networks.

## Graph Matching

Given a time series of networks $G_{1},\cdots,G_{T}$ on the same set of
vertices $V$, a 1-1 vertex correspondence across all the networks facilitates
joint spectral embedding, which is a key step in the mirror estimation method
we will present in the next section, “Discovering underlying dynamics in time
series of networks”. Such a vertex correspondence may be available a priori in
labeled networks, or it can be inferred from unlabeled networks. This
inference problem – the so-called graph matching problem – is to find an
alignment of vertices between two graphs such that the corresponding edge
differences are minimized. We denote $A$, $B$ as two $n\times n$ adjacency
matrices for two graphs with $n$ vertices each. Then the graph matching
problem is to find the permutation matrix $P$ that maximizes the objective
function

$f(A,B;P)=-||A-P^{\top}BP||_{F}^{2}$

for $P\in\mathcal{P}$ where $\mathcal{P}=\\{P\in\\{0,1\\}^{n\times
n}:P^{\top}\mathbf{1}=P\mathbf{1}=\mathbf{1}\\}$,
$\mathbf{1}=(1,1,...1)^{\top}$, and $||\cdot||$ denotes the Frobenius norm.
This formulation is equivalent to maximizing

$f(A,B;P)=trace(APB^{T}P^{T}).$

Because solving this optimization problem is combinatorically difficult,
approximation algorithms have been proposed. We use the Fast Approximate
Quadratic (FAQ) assignment algorithm [2] to obtain an approximate solution.
FAQ iteratively finds a local solution to the relaxed problem – expanding the
constraint set to the convex hull of $\mathcal{P}$ – and then projects the
solution back to $\mathcal{P}$, and has been shown empirically to be
competitive with or superior to other state-of-art methods. In Section Results
we use FAQ to demonstrate that the given putative 1-1 correspondence for each
pair of networks is close to the solution to the corresponding graph matching
problem.

## Discovering underlying dynamics in time series of networks

To discover the underlying dynamics in time series of networks, we use the
model and method proposed in [1]. First we introduce the generative joint
model for time series of networks. We consider $T$ networks, each containing
$n$ vertices, with adjacency matrices $A_{t}$, $t\in\\{1,2,...,T\\}$. In the
model, each vertex is associated with a time varying $d$-dimensional latent
vector. These vectors $X_{t}$ are each one realization from a stochastic
process, called the latent position process (LPP) – each $X_{t}$ is a
$d$-dimensional random variable. For $n$ vertices, we generate $n$ i.i.d.
samples from the LPP, which collection then forms the latent position matrices
$\\{\mathbf{X}_{t}\\}$, where $\mathbf{X}_{t}\in\mathcal{R}^{n\times d}$ for
$t\in\\{1,2,...,T\\}$. The connection probability between vertex $i$ and
vertex $j$ at time $t$ is the inner product of the associated latent vectors
at time $t$. That is $\mathbb{E}(A_{t})=\mathbf{X}_{t}\mathbf{X}^{T}_{t}$.
Note that each network corresponds to a latent position random variable. Thus
we can capture the distance between graphs using the corresponding random
variables. We define the distance

$d_{MV}(X_{t},X_{t^{\prime}})=\min_{W\in\mathcal{O}^{d\times
d}}\|\mathbb{E}[(X_{t}-WX_{t^{\prime}})(X_{t}-WX_{t^{\prime}})^{\top}]\|^{1/2}_{2},$

where $\mathcal{O}^{d\times d}$ is the set of orthogonal transformation
matrices with dimension $d$. When $X_{t}$ and $X_{t^{\prime}}$ are centered,
the $d_{MV}$ distance can be interpreted as the maximum directional variation
for the random vector $X_{t}-WX_{t^{\prime}}$, where $W$ is an orthogonal
transformation used to align $X_{t}$ and $X_{t^{\prime}}$. We evaluate this
distance for every pair of random variables in the LPP and obtain a $T\times
T$ distance matrix $\mathcal{D}$. Then we apply classical multidimensional
scaling (MDS) [7] to $\mathcal{D}$ to get a low-dimensional Euclidean
representation of underlying network structure, called the mirror,
$\\{\psi(t)\\}$. In practice, the LPP is unknown and only network realizations
$\\{A_{t}\\}$ are observed. For the $n\times n$ symmetrized adjacency matrix
$A_{t}$, we use adjacency spectral embedding (ASE) [8] to obtain
$\hat{\mathbf{X}}_{t}=U\Sigma^{1/2}$, where the diagonal matrix $\Sigma$
contains the top $d$ eigenvalues of $A_{t}$ and $U$ contains the associated
eigenvectors. Then we use

$\hat{d}_{MV}(\hat{\boldsymbol{X}}_{t},\hat{\boldsymbol{X}}_{t^{\prime}})=\min_{W}\frac{1}{\sqrt{n}}\|\hat{\boldsymbol{X}}_{t}-\hat{\boldsymbol{X}}_{t^{\prime}}W\|^{1/2}_{2}$

to estimate the pairwise distance between networks, yielding
$\hat{\mathcal{D}}$. Applying MDS to $\hat{\mathcal{D}}$ yields the mirror
estimate $\\{\hat{\psi}(t)\\}$. When the mirror exhibits a manifold structure,
we can further simplify the change point detection problem by applying the
manifold learning method isometric mapping (ISOMAP) [9] to
$\\{\hat{\psi}(t)\\}$. This yields the _iso-mirror_ , which captures the
geodesic distance along the mirror and preserves it in lower dimensions with
Euclidean distance. Subsequent inference is then performed using the iso-
mirror representation. For convenience, the iso-mirror representation will
also be denoted as $\\{\hat{\psi}(t)\\}$.

## Results

### Time series of organoid networks data

All results in this section are based on data collected from well 8. For
analogous results from well 5, please refer to the Appendix. For well 8, the
time series of organoid networks consists of 45 time stamps
$\\{1,2,...,45\\}$; each time stamp corresponds an effective connectivity
graph $G_{t}$ with adjacency matrix $A_{t}$. Each graph is directed, weighted,
and hollow. We symmetrize the directed graphs, and use ranks in place of the
raw edge weights. All graphs have the same vertex set $V=\\{1,2,...n\\}$ with
$n=|V|=127$, although some of the graphs contain isolated vertices. See Figure
1.

Figure 1: The number of non-isolated vertices and the number of edges for the
graphs at each of the 45 time stamps. The number of edges are counted after
symmetrizing the directed graph.

### Putative 1-1 correspondence

For this time series of organoid networks, inferred neuron location gives a
putative 1-1 vertex correspondence across the graphs. We assess this
correspondence via graph matching using the objective function value (OFV)
$f(A,B;P)$. The OFV for the putative 1-1 correspondence is given by $f(A,B;I)$
where $I$ is the identity matrix. We denote the FAQ output for matching
adjacency matrices $A$ and $B$ initialized at $C$ as $P_{A,B;C}$. Typically we
choose the barycenter $b=\frac{\mathbf{1}\mathbf{1}^{T}}{n}$ as the initial
point. For all times $i\in\\{1,2,...44\\}$, we consider $A_{i}$, $A_{i+1}$ and
FAQ yields $P_{A_{i},A_{i+1};b}$ (denoted $P_{i,i+1;b}$ for short). In Figure
2 we compare $f(I)=f(A_{i},A_{i+1};I)$ and
$f(P_{i,i+1;b})=f(A_{i},A_{i+1};P_{i,i+1;b})$. Although $f(P_{i,i+1;b})$ is
always larger than $f(I)$, the two OFVs are close to each other for all time
stamps, indicating that the putative 1-1 correspondence is close to FAQ’s
solution.

(a) (b)

Figure 2: Comparison of OFV using $I$ and $P_{i,i+1;b}$ for 44 pairs of
graphs, demonstrating that FAQ increases the OFV only slightly. (a): $f(I)$
and $f(P_{i,i+1;b})$. (b): $\frac{f(P_{i,i+1;b})}{f(I)}$.

(a) (b)

Figure 3: Matching $A_{28}$ and $A_{29}$. Comparison of OFV for using $I$,$R$,
$P_{28,29;b}$, $P_{28,29;I}$ and $P_{28,29;R}$. (a): Histogram demonstrating
that OFV for $R$ is not nearly as good as for the others. (b): Enlargement of
the far right portion of the top figure, demonstrating that FAQ output at
different initial points $I$, $b$ and $R$ improves the OFV only slightly
compared to the putative 1-1, $I$.

To further asses the putative 1-1 correspondence, we consider a specific pair
of adjacency matrices: $A_{28}$, $A_{29}$. We uniformly generate 100,000
random permutation matrices $R$ and evaluate $f(R)=f(A_{28},A_{29};R)$. We
plot the histogram in Figure 3. We also indicate $f(A_{28},A_{29};I)$,
$f(A_{28},A_{29};P_{28,29;I})$, $f(A_{28},A_{29};P_{28,29;b})$ and
$f(A_{28},A_{29};P_{28,29;R})$, where $R$ are 100 randomly drawn permutation
matrices. We see that the putative 1-1 correspondence performs better than all
100,000 instantiations of $f(R)$ and is close to $P_{28,29;b}$, $P_{28,29;I}$
and $P_{28,29;R}$. Thus we conclude that the putative 1-1 correspondence is
sufficiently accurate, and we will proceed apace for mirror estimation and
change point detection.

### Change point detection

For the 45 graphs, time stamps are from 1 to 244, in days. We choose time
stamps in [150,230] to avoid growth and death regimes.

For these graphs, we find the largest common connected component, which
contains 112 vertices. The average number of edges for the largest common
connected component is approximately 6130. We use the putative 1-1 vertex
correspondence across time. We apply our mirror estimation method to this time
series of networks, and ISOMAP manifold learning yields the 1-dimensional
representation of the dynamics $\\{\hat{\psi}_{t}\\}$. We choose 10 time
stamps shown in Figure 4. As we see, the representation is approximately
piecewise linear with an evident change of slope at $t=4$, day 188.

Figure 4: Discovering a change point in time series of inferred effective
connectivity organoid networks via the iso-mirror. The x-axis is time stamp,
in days. The y-axis is the ISOMAP representation of the estimated mirror
$\\{\hat{\psi}_{t}\\}$. This indicates a change in the network dynamics at day
188.

Assuming the true underlying 1-dimensional representation $\psi(t),t\in[0,T]$,
is piecewise linear and continuous, it is natural to define the change point
$t^{*}$ as the point when the slope changes. If we assume there is only one
change point, then we can write

$\psi(t)=\beta_{0}+\beta_{1}t+\beta(t-t^{*})I(t>t^{*}).$

Both the change point detection algorithm from [3] and the break point
estimation for piecewise linear models from [4] yield an estimated change
point $\hat{t^{*}}=4$, day 188, which coincides with neuroscientifically
significant developmental changes – inhibitory neurons start appearing and the
percentage of astrocytes increases dramatically – as described in [10]. Note
that the emergence of astrocytes does not happen at once but builds over time,
so there is no one precise date for the change point and detection of a time
coinciding with this change can be but suggestive.

## Conclusion

Reconstruction of effective connectivity networks of electrophysiologically
active brain organoids reflect their structural (increasing number of neurons
(nodes) and connections (edges)) and electrical development over time, as
previously demonstrated in [10].

By applying the spectral mirror estimation method to the time series of
organoid networks, we obtain a 1-dimensional iso-mirror representation of
dynamic inferred effective connectivity organoid networks. Two change point
detection algorithms successfully detect a change at day 188. At approximately
188 days (~6 months), cortical organoids start showing inhibitory neurons and
the percentage of astrocytes increases from 5% to 30-40% [10] resulting in
added complexity in the activity and network distribution of brain organoids.

There are several change point detection algorithms available for analyzing
time series of graphs, including the one proposed in [11]. However, it is
important to note that the spectral mirror estimation method used in our study
is not restricted to change point detection alone. Instead, it provides a low-
dimensional (in our case, one-dimensional Euclidean) representation of network
dynamics, enabling us to visualize network evolution. As illustrated in Figure
4, we observe piecewise linear structure and apply segment regression to
identify a significant increase in slope after day 188. This suggests that
organoid graphs drift continuously over time, but their rate of drift
accelerates significantly after day 188. Our method is preferred in this
regard as it provides us with more than just one change point. Additionally,
since the iso-mirror represents the underlying LPP, the detected change point
in the iso-mirror reflects a fundamental change in the underlying generative
LPP of the time series of networks, which may not necessarily correspond to
any specific network measure. Furthermore, the spectral mirror method can be
applied to other dynamic graphs that meet our model assumptions, namely that
the time series of networks are generated from a time series of latent
position random networks whose vertices have independent, identically
distributed latent positions given by an LPP.

Future work includes addressing two major theoretical issues of note, to make
the change point inference formally principled. First, the theory in [1]
requires a known 1-1 vertex correspondence across time. It remains to study
the effect of errors in this correspondence, such as those inherent in our
putative 1-1 correspondence deemed sufficiently accurate for practical
purposes. In addition, it remains to study the entry-wise behavior of the
error term $\epsilon(t)=\hat{\psi}(t)-\psi(t)$. For example: in [3] the proof
of consistency of the change point estimator requires the error process to be
stationary; in [4] the $\epsilon(t)$ are assumed to be i.i.d. normal to
construct a confidence interval for $\hat{t^{*}}$. For now, [1] has shown that
$\hat{\psi}(t)$ converges to $\psi(t)$ in Frobenius norm, that is
$\sum_{t=1}^{T}(\epsilon(t))^{2}\to 0$ with high probability, which is
insufficient to conclude normality or stationarity. What’s more, this
convergence result is for the mirror rather than the iso-mirror. As for
whether the same conclusion can be extended to the iso-mirror, further
investigation is needed.

## Appendix

### Time series of organoid networks data for well 5

For well 5, there are 40 graphs with time stamps [1,229] and all of them have
the same vertex set with $|V|=128$, although some of the graphs contain
isolated vertices. See figure 5. Each graph is directed, weighted, and hollow.
We symmetrize the directed graphs, and use ranks in place of the raw edge
weights.

Figure 5: The number of non-isolated vertices and the number of edges for the
graphs at each of the 40 time stamps for well 5. The number of edges are
counted after symmetrizing the directed graph.

### Putative 1-1 correspondence for well 5

We assess the putative 1-1 correspondence exactly the same way as in the paper
and Figure 6 indicates that the putative 1-1 correspondence is close to FAQ’s
solution.

(a) (b)

Figure 6: Comparison of OFV using $I$ and $P_{i,i+1;b}$ for 39 pairs of
graphs, demonstrating that FAQ increases the OFV only slightly. (a): $f(I)$
and $f(P_{i,i+1;b})$. (b): $\frac{f(P_{i,i+1;b})}{f(I)}$.

We also consider a specific pair of adjacency matrices: $A_{28}$, $A_{29}$ to
assess the 1-1 putative correspondence. See figure 7. The result is similar to
well 8 and we conclude the 1-1 putative correspondence is accurate enough.

Figure 7: Matching $A_{28}$ and $A_{29}$. Comparison of OFV for using $I$,$R$,
$P_{28,29;b}$, $P_{28,29;I}$ and $P_{28,29;R}$. Left: Histogram demonstrating
that OFV for $R$ is not nearly as good as for the others. Right: Enlargement
of the far right portion of the left figure, demonstrating that FAQ output at
different initial points $I$, $b$ and $R$ improves the OFV only slightly
compared to the putative 1-1, $I$.

### Change point detection for well 5

For the 40 graphs, time stamps are from 1 to 244, in days. We choose the same
time stamps in [150,230] as in the paper. We find the largest common connected
component for these graphs. We use the putative 1-1 vertex correspondence
across time. We apply our mirror estimation method to this time series of
networks, and ISOMAP manifold learning yields the 1-dimensional representation
of the dynamics $\\{\hat{\psi}_{t}\\}$. We choose the same 10 time stamps as
in the paper shown in Figure 8. As we see, the representation is approximately
piecewise linear. The break point estimation for piecewise linear models from
[4] yield an estimated change point $\hat{t^{*}}=4$, day 188.

Figure 8: For well 5, iso-mirror result.

## Ethical interest

not applicable

## Competing interests

The authors declare that they have no competing interests.

## Abbreviations

MEA: Multi-Electrode Array
PCA: Principal Component Analysis
FAQ: Fast Approximate Quadratic
LPP: Latent Position Process
MDS: Multidimensional Scaling
ISOMAP: Isometric Mapping
ASE: Adjacency Spectral Embedding
OFV: Objective Function Value

## Funding

TC’s work is partially supported by the Johns Hopkins Mathematical Institute
for Data Science (MINDS) Data Science Fellowship and Azure sponsorship credits
granted by Microsoft’s AI for Good Research Lab.

## Availability of data and materials

The organoid data and code are available at
https://github.com/youngser/organoid.

## Authors’ contributions

TC, ZL, AA, CEP developed the theory. ARM, FP, GAS provided the data. TC, YP,
AS, BPD designed & implemented the methods. TC, YP conducted the experiments.
TC, YP, CEP wrote the manuscript. WY, CW, JTV, CEP guided the whole process.
All authors read and approved the manuscript.

## References

  * [1] Athreya, A., Lubberts, Z., Park, Y., Priebe, C.E.: Discovering underlying dynamics in time series of networks. arXiv preprint arXiv:2205.06877 (2022)
  * [2] Vogelstein, J.T., Conroy, J.M., Lyzinski, V., Podrazik, L.J., Kratzer, S.G., Harley, E.T., Fishkind, D.E., Vogelstein, R.J., Priebe, C.E.: Fast approximate quadratic programming for graph matching. PLoS One 10(4), 0121002 (2015)
  * [3] Bücher, A., Dette, H., Heinrichs, F.: Are deviations in a gradually varying mean relevant? a testing approach based on sup-norm estimators. The Annals of Statistics 49(6), 3583–3617 (2021)
  * [4] Muggeo, V.M.: Interval estimation for the breakpoint in segmented regression: A smoothed score-based approach. Australian & New Zealand Journal of Statistics 59(3), 311–322 (2017)
  * [5] Muguruma, K., Nishiyama, A., Kawakami, H., Hashimoto, K., Sasai, Y.: Self-organization of polarized cerebellar tissue in 3d culture of human pluripotent stem cells. Cell reports 10(4), 537–550 (2015)
  * [6] Puppo, F., Pré, D., Bang, A.G., Silva, G.A.: Super-selective reconstruction of causal and direct connectivity with application to in vitro ipsc neuronal networks. Frontiers in Neuroscience 15, 647877 (2021)
  * [7] Torgerson, W.S.: Multidimensional scaling: I. theory and method. Psychometrika 17(4), 401–419 (1952)
  * [8] Athreya, A., Fishkind, D.E., Tang, M., Priebe, C.E., Park, Y., Vogelstein, J.T., Levin, K., Lyzinski, V., Qin, Y., Sussman, D.L.: Statistical Inference on Random Dot Product Graphs: a Survey. J. Mach. Learn. Res. 18(226), 1–92 (2018)
  * [9] Tenenbaum, J.B., Silva, V.d., Langford, J.C.: A global geometric framework for nonlinear dimensionality reduction. Science 290(5500), 2319–2323 (2000)
  * [10] Trujillo, C.A., Gao, R., Negraes, P.D., Gu, J., Buchanan, J., Preissl, S., Wang, A., Wu, W., Haddad, G.G., Chaim, I.A., et al.: Complex oscillatory waves emerging from cortical organoids model early human brain network development. Cell stem cell 25(4), 558–569 (2019)
  * [11] Wang, D., Yu, Y., Rinaldo, A.: Optimal change point detection and localization in sparse dynamic networks (2021)