# High-energy Neutrinos from Outflows Powered by Kicked Remnants of Binary
Black Hole Mergers in AGN Accretion Disks

Zhi-Peng Ma Department of Astronomy, School of Physics, Huazhong University of
Science and Technology, Wuhan, Hubei 430074, China  Kai Wang Department of
Astronomy, School of Physics, Huazhong University of Science and Technology,
Wuhan, Hubei 430074, China Kai Wang<EMAIL_ADDRESS>
###### Abstract

Merging of stellar-mass binary black holes (BBH) could take place within the
accretion disk of active galactic nuclei (AGN). The resulting BH remnant is
likely to accrete the disk gas at a super-Eddington rate, launching a fast,
quasi-spherical outflow (wind). Particles will be accelerated by shocks driven
by the wind, subsequently interacting with the shocked disk gas or radiation
field through hadronic processes and resulting in the production of high-
energy neutrinos and potential electromagnetic (EM) emissions. This study
delves into the intricate evolution of the shock driven by the remnant BH wind
within AGN disks. Subsequently, we calculated the production of neutrinos and
the expected detection numbers for a single event, along with their
contributions to the overall diffuse neutrino background. Our analysis,
considering various scenarios, reveals considerable neutrino production and
possible detection by IceCube for nearby events. The contribution of the
remnant BH winds on the diffuse neutrino background is minor due to the low
event rate density, but it can be improved to some extent for some optimistic
parameters. We also propose that there could be two neutrino/EM bursts, one
originating from the premerger BBH wind and the other from the remnant BH
wind, with the latter typically having a time gap to the GW event of around
tens of days. When combined with the anticipated gravitational waves (GW)
emitted during the BBH merger, such a system emerges as a promising candidate
for joint observations involving neutrinos, GWs, and EM signals.

Neutrino astronomy — High energy astrophysics — Active galactic nuclei —
Gravitational waves — Black holes

##  1 Introduction

It is believed that many stars and stellar remnants can be embedded within the
accretion disk of AGN (Artymowicz et al., 1993; Cantiello et al., 2021;
Dittmann et al., 2021; Collin & Zahn, 1999; Fabj et al., 2020; Wang et al.,
2011). The abundant stellar remnants would migrate inward the disk and be
trapped at the specific orbit (Bellovary et al., 2016; Tagawa et al., 2020),
leading to the high possibility of compact object mergers (Cheng & Wang, 1999;
McKernan et al., 2020), such as binary neutron star (BNS) mergers (McKernan et
al., 2020), binary black hole (BBH) mergers, and NS-BH (NSBH) mergers (Bartos
et al., 2017; Leigh et al., 2018; Kimura et al., 2021; Tagawa et al., 2022).
Especially, among these compact object mergers in the AGN disk, BBH mergers
have attracted the most attention since they have been prime targets for
Earth-based gravitational wave (GW) detectors. Another important point is that
BBH mergers embedded in AGN disks can potentially power their electromagnetic
counterparts by accreting the abundant disk gas (Wang et al., 2021; Kimura et
al., 2021; Tagawa et al., 2022; Chen & Dai, 2024), while for the normal
environment, the accompanying EM radiations with BBH mergers are not expected
because there are no accretion materials to power a jet. Recently, seven AGN
flares have been suggested to have a statistical association with LIGO/Virgo
BBH merging events (Graham et al., 2023). In addition, the Zwicky Transient
Facility (ZTF) has reported an optical counterpart candidate to the GW event
GW190521 (Graham et al., 2020), which comes from the most massive stellar-mass
BBH system to date (up to $150M_{\odot}$) (Abbott et al., 2020).

Besides the GW and EM signals, astrophysical activities within AGN disks,
e.g., supernovae (Zhu et al., 2021; Zhou et al., 2023), gamma-ray bursts (Zhu
et al., 2021), stellar-mass BH jets (Tagawa et al., 2023) and the remnant BH
jets of BBH mergers (Zhou & Wang, 2023), can also be the potential high-energy
neutrino source. For the BBH system embedded in AGN disks, the merged BH
remnant will be kicked into an intact AGN disk environment. Then, the
surrounding disk gas will be accreted with a super-Eddington accretion rate,
and the jet will be launched by Blandford–Znajek (BZ) mechanism. The
interaction between the postmerger BZ jet and the AGN disk materials will
create diverse shocks, accelerate particles, and produce neutrinos by hadronic
processes (Zhou & Wang, 2023). In addition to the jet, the remnant BH of the
BBH post-merger will also power the radiation-driven wind (Jiang et al., 2014;
Ohsuga et al., 2005; Skadowski et al., 2013), and this wind interacting with
the disk medium may also play an important role in multimessenger signals of
BBH mergers in AGN disks.

Figure 1: Schematic picture for the system. Left: The BBH wind will flow out
the ambient gas and create a cavity around the merged BH. Right: When the
merged BH is kicked out of the cavity into a new disk environment, it will
accrete the disk gas and power a new wind, and this wind will produce a
forward shock sweeping across the disk gas and a reverse shock sweeping across
the wind, dividing the system into four regions. See texts for more details.

BBH mergers are the potential GW-neutrino association source, and Zhou & Wang
(2023) estimate the joint detection of GW from BBH merger and high-energy
neutrinos from the jet of the postmerger remnant BH within AGN disks. Compared
with the jet of remnant BH, the sub-relativistic wind driven by the merged BH
is quasi-isotropic without the beaming effect, leading to a more promising
detection of GW-neutrino association for the current small number of detected
BBH mergers. In this paper, we study the neutrino emission from the
interaction between the wind driven by the remnant BH and the disk medium. In
our scenario, the premerger BBH system will produce a radiation-driven wind,
and create a cavity around the merged BH in the AGN disk. When the merged BH
is kicked into a new, gas-abundant disk environment, it will continue to
accrete the disk gas and power a new material wind (Kimura et al., 2021). Such
a wind will generate a forward shock in the disk gas, hence particles can be
accelerated at the shock front. Then, the accelerated protons will interact
with the disk materials and thermal radiations of the shocked gas to produce
high-energy neutrinos. We show a schematic picture of our scenario in Fig. 1.

This paper is organized as follows. We show the equations and solutions for
shock dynamics in Section 2. The method for estimating neutrino emission is
shown in Section 3. The neutrino energy fluence, expected numbers, and
contributions to diffuse neutrino background are shown in Section 4. We
summarize our work in Section 5.

##  2 Dynamical model

In this section, we describe the properties of the disk, including density
profile and temperature distribution. Then, we set up the dynamic equations
for the forward shock and illustrate the solutions. Notice that we only focus
on the outflow driven by the merged BH. However, such equations for shock
dynamics and neutrino productions are also valid for the premerger outflow and
the massive BH formed through other ways (such as Core-collapsed).

###  2.1 Disk Model and Outflow Luminosity

Since BBH and the final merged BH would migrate within the disk of the
supermassive black hole (SMBH), we set the distance $r$ from the merged BH to
SMBH as a free parameter, in the range of $10^{3}R_{\rm sch}\leq r\leq
10^{4}R_{\rm sch}$, where $R_{\rm sch}=2GM_{\rm SMBH}/c^{2}$ is the
Schwarzschild radius of SMBH, G is the gravitational constant, $M_{\rm SMBH}$
is the SMBH mass, and c is the speed of light. We adopt a modified
$\alpha$-disk model (Cantiello et al., 2021) to calculate the gas temperature
$T_{\rm mid}$ and gas density profile $\rho_{\rm mid}$ for the disk midplane,
also the scale height H and surface density $\Sigma$ for the disk. By
considering gravitational stability, all these physical quantities have
different expressions in the inner and outer regions of the disk. The detailed
formulas are given in the Appendix A. Once the expressions for the midplane
are confirmed, one can get the vertical density profile in Gaussian form (Zhu
et al., 2021) (with fixed SMBH mass),

$\rho_{\rm disk}(r,h)=\rho_{\rm mid}(r){\rm exp}(-h^{2}/2H^{2}),$ (1)

where $h=r/H$ is the aspect ratio. For the gas temperature distribution in the
vertical direction, we simply assume that $T_{\rm disk}(r,h)=T_{\rm mid}(r)$.
The pressure can be calculated by $P_{\rm disk}(r,h)=\rho_{\rm
disk}/m_{p}k_{b}T_{\rm disk}$ under the assumption of the ideal gas, where
$k_{b}$ is the Boltzmann constant. We note that the pressure of the disk gas
may not follow the ideal gas relation, but such an estimation will not make a
significant deviation in calculations, because the deceleration of the shocked
gas is mainly caused by the gravity of the merged black hole rather than the
pressure of the disk gas.

When the merged BH is kicked out of the cavity into a gas-abundant region on
the disk, it will accrete the disk gas in the way of Bondi-Hoyle-Littleton
(BHL) accretion and power a fast wind. The accretion rate can be estimated
through (Bondi, 1952; Hoyle & Lyttleton, 1939; Edgar, 2004)

$\dot{M}_{\rm bh}(r)=\frac{4\pi G^{2}M^{2}_{\rm bh}\rho_{\rm mid}}{v^{3}_{\rm
kick}},$ (2)

where $M_{\rm bh}$ is the mass of the merged black hole, which is fixed as
$150M_{\odot}$. $v_{\rm kick}\simeq 1\times 10^{7}\rm{cm/s}$ is the kick
velocity (Campanelli et al., 2007; Herrmann et al., 2007; Gonzalez et al.,
2007). The kinetic power of the wind can be estimated through

$L_{\rm w}=\eta_{w}\dot{M}_{\rm bh}v^{2}_{w},$ (3)

where $\eta_{w}$ is the energy distribution factor for the BHL accretion rate
and $v_{w}$ is the velocity of the wind.

###  2.2 Dynamic Equations

We adopt a 1D shock model whose physical picture is similar to that of stellar
wind or AGN wind. We assume a steady spherical wind launched from BH with a
velocity $v_{\rm w}=10^{10}\rm cm/s$. The disk gas around BH is also assumed
to be spherically symmetric, based on the direction vertical to the midplane.
When the fast wind is injected into the disk, a forward shock and a reverse
shock will be formed at the radii of $R_{\rm s}$ and $R_{\rm rs}$
respectively. Two shock radii together with the contact discontinuity at
radius $R_{\rm cd}$ divide the system into four regions: (a) the unshocked
wind moving with $v_{w}$, (b) the shocked winds, (c) the shocked disk gas, and
(d) the unshocked disk gas. The schematic picture can be seen in Fig. 1.

The forward shock will sweep the disk gas in Region (d) and compress it into a
shell. For simplicity, we assume that $R_{s}\simeq R_{cd}$ (thin shell
approximation). Under such an assumption (Weaver et al., 1977; Faucher-Giguère
& Quataert, 2012; Wang & Loeb, 2015; Liu et al., 2018), the forward shock and
the gases in Region (c) have a unified velocity, notated $v_{\rm s}$. For the
shocked wind in Region (b), we follow the treatment proposed by Liu et al.
(2018), which considers a steady shocked wind with homogeneous number density
$n_{\rm sw}$, temperature $T_{\rm sw}$ and gas pressure $P_{\rm sw}$ in the
region at any time111Such an assumption is validated when the local sound
cross time $\sqrt{P_{\rm sw}/\rho_{\rm sw}}$ is shorter than the system
dynamical timescale, which is generally satisfied in our scenario.. The number
density is $n_{\rm sw}=4n_{w}=\dot{M}_{w}/\pi R^{2}_{\rm rs}m_{\rm p}v_{w}$,
where the wind injection rate is $\dot{M}_{w}=2L_{w}/v^{2}_{w}$. Under the
condition of mass conservation, $R^{2}v_{\rm sw}(R)$ is a constant in Region
(b), where $R$ is the distance to the merged BH and $v_{\rm sw}(R)$ is the
velocity of shocked wind at radius $R$. There are two boundary conditions for
$v_{\rm sw}(R)$ at $R_{\rm s}$ and at $R_{\rm rs}$, which are $v_{\rm
sw}(R_{\rm s})=v_{\rm s}$ and $v_{\rm sw}(R_{\rm rs})=v^{\prime}_{\rm sw}$,
where $v^{\prime}_{\rm sw}$ is the velocity of the shocked wind just behind
$R_{\rm rs}$. By combining two boundary conditions with mass conservation, one
can find the relation,

$v^{\prime}_{\rm sw}=(\frac{R_{\rm s}}{R_{\rm rs}})^{2}v_{s}.$ (4)

Together with the Rankine–Hugoniot jump relation $v_{\rm w}-v_{\rm
rs}=4(v^{\prime}_{\rm sw}-v_{\rm rs})$, where $v_{\rm rs}$ is the velocity of
the reverse shock, one can get the useful expression for $v_{\rm rs}$, that
is,

$v_{\rm rs}=\frac{4}{3}(\frac{R_{\rm s}}{R_{\rm rs}})^{2}v_{\rm
s}-\frac{1}{3}v_{\rm w}.$ (5)

Once the evolution of the forward shock is known, one can estimate the motion
of reverse shock through Eq.(5) directly.

Next, we need to describe the evolution of the forward shock. For the shocked
disk gas in Region (c), we have the motion equation,

$\frac{d(M_{\rm s}v_{\rm s})}{dt}=4\pi R^{2}_{\rm s}(P_{\rm tot}-P_{\rm
disk})-\frac{GM_{\rm tot}M_{\rm s}}{R^{2}_{\rm s}},$ (6)

where $P_{\rm tot}=P_{\rm sw}+P_{\rm ph}$ is the pressure provided by the gas
and photons in Region (b), which pushes the shell forward. The reason for
considering the photon pressure is that the radiation from electrons and
protons in Region (b) hardly escapes until the shock breaks out, and most of
the photons move with the shocked wind and provide pressure together with the
matter. $P_{\rm disk}$ is the pressure of the disk gas, which resists
expansion of the shell. The second term on the right-hand side represents the
gravity exerted by the total gravitational mass within $R_{\rm s}$, which
mainly causes the deceleration of the shell, where $M_{\rm tot}=M_{\rm
bh}+M_{\rm s}/2$. $M_{\rm s}$ is the total mass of the shocked disk gas, which
can be obtained by integrating the mass swept by the forward shock: $M_{\rm
s}(R_{\rm s})=\int_{0}^{R_{\rm s}}{4\pi R^{2}\rho_{\rm disk}(r,R)dR}$. We
ignore the effect of gravity of the SMBH, since the Hill radius $R_{\rm
Hill}=(M_{\rm bh}/3M_{\rm SMBH})^{1/3}r$ is comparable to the scale height
$H$, suggesting that gravity of the merged stellar-mass BH is dominated.

Figure 2: The evolution of the forward shock velocity $v_{\rm s}$ (top-left),
the total pressure of the shocked wind $P_{\rm tot}$ (top-right), and the
thermal energy injection rate of shocked gas $L_{\rm th}$ (bottom), in the
case of $M_{\rm SMBH}=10^{7}M_{\odot},r=10^{3}R_{\rm sch}$. Solid lines in
different colors represent different $\eta_{\rm w}$ selections. The vertical
red dashed line represents the scale height of the AGN disk. The vertical
dark-blue, blue, and red dot-dashed line represents the breakout radius
$R_{\rm br}$ in each $\eta_{\rm w}$ selection. Note that the red dot-dashed
line is overlapped by the dark-blue line and becomes invisible.

To calculate the pressure $P_{\rm tot}$, we have the equation for the total
internal energy rate of Region (b),

$\frac{dE_{\rm tot}}{dt}=\frac{1}{2}\dot{M}_{w}(1-\frac{v_{\rm rs}}{v_{\rm
w}})(v^{2}_{\rm w}-v^{\prime}_{\rm sw})-4\pi R^{2}_{\rm s}P_{\rm tot}v_{\rm
s}.$ (7)

The first term on the right-hand side of the above equation represents the
injection rate of matter thermal energy in Region (b) (the difference between
the wind kinetic power and shocked wind kinetic power). The second term is the
energy loss due to the expansion. We also have the energy-pressure relation
for the ideal gas,

$E_{\rm tot}=\frac{3}{2}P_{\rm tot}V_{\rm sw}=2\pi P_{\rm tot}(R^{3}_{\rm
s}-R^{3}_{\rm rs}),$ (8)

with the volume of $V_{\rm sw}=4\pi/3(R^{3}_{\rm s}-R^{3}_{\rm rs})$ . Notice
that Eq.(7) and Eq.(8) are the simplified equations, the original version is a
set of differential equations describing the energy evolution of protons,
electrons, and photons respectively. The original energy-pressure relation
also includes that of photons. For simplification, We adopt this form in
calculations, but the final result will not deviate so much. One can see the
detailed argument in Appendix B. When the evolution of $v_{\rm s}$ is
obtained, we can also estimate the thermal energy injection rate of the
forward shock by $L_{\rm th}=4\pi R^{2}_{\rm s}u_{\rm s}v_{\rm s}$, where
$u_{\rm s}=\frac{3}{2}n_{\rm s}k_{\rm b}T_{\rm s}=\frac{9}{32}\rho_{\rm
s}v^{2}_{\rm s}$ is the thermal energy density of the shocked disk gas,
$T_{\rm s}=3m_{\rm p}v^{2}_{\rm s}/16k_{\rm b}$ is the temperature of the
shocked disk gas, $\rho_{\rm s}=4\rho_{\rm disk}$ is the density of the
shocked disk gas.

By equating the radiation-diffusive velocity and the shock velocity, one can
obtain the characteristic optical depth $\tau_{\rm c}=c/v_{\rm s}$. The shock
will break out at the radius $R_{\rm br}$ where $\tau_{\rm c}=\tau_{\rm\rm
m}$, $\tau_{\rm m}\simeq(\rho_{\rm disk}/m_{\rm p})\sigma_{\rm T}R$ is the
material optical depth ahead of the shock (Wang et al., 2019). When $R_{\rm
s}>R_{\rm br}$, the shock becomes collisionless and the particle acceleration
starts to occur.

###  2.3 Solutions for Dynamic Equations

Figure 3: The timescales for various processes of protons when the forward
shock is at $R_{\rm s}=1\times 10^{14}\rm cm$ (left) and $R_{\rm s}=3\times
10^{14}\rm cm$ (right) under the case of $M_{\rm
SMBH}=10^{7}M_{\odot},r=10^{3}R_{\rm sch}$, where $R_{\rm br}\sim 2\times
10^{14}\rm cm$. The maximum proton energy can be obtained by equating the
timescale of the acceleration and the most effective cooling process. Note
that the particle acceleration in the left panel (prior to shock breakout) is
not allowed, and it is solely used for comparing the timescales of various
processes at different radii.

We set the wind initially launched from the location $R=H/100$ to the merged
BH in the midplane. Based on Eq.(6)-Eq.(8), we can obtain the evolution of
$v_{\rm s}$, $P_{\rm tot}$ and $L_{\rm th}$, we show our results in Fig 2 with
$M_{\rm SMBH}=10^{7}M_{\odot},r=10^{3}R_{\rm sch}$, and diverse $\eta_{\rm
w}$. We can see that the forward shock decelerates at the beginning and, when
it passes through the scale height of the disk, it will continue to accelerate
even beyond the breakout radius. The injected thermal luminosity is almost
unchanged before the scale height, and then it becomes decreasing. When the
shock reaches the breakout radius, $L_{\rm th}$ is relatively low compared to
the initial value. We introduce the factor $\epsilon=0.1$ to represent the
fraction of thermal luminosity converted to cosmic-ray luminosity (Caprioli &
Spitkovsky, 2014), that is, $L_{\rm cr}=\epsilon L_{\rm th}$.

In our solutions, the shock breakout radius is always larger than the scale
height of the disk, i.e. $R_{\rm br}>H$. Kimura et al. (2021) proposed an
analytic formula to estimate the breakout time:
$t^{\prime}_{br}\simeq\left(\frac{\rho_{\rm disk}H^{5}}{0.53L_{\rm
w}}\right)^{1/3}$. In the case of $\eta_{\rm w}=0.1$, this formula gives
$t^{\prime}_{\rm br}\simeq 80h$, which is comparable to our result $t_{\rm
br}=44h$.

##  3 Neutrino Production

After the shock breakout, particles are accelerated at the forward shock front
and produce high-energy neutrinos through hadronic processes. We consider that
the accelerated protons follow the power-law distribution, namely, $dN_{\rm
p}/dE_{\rm p}\propto E^{-s}_{\rm p}$, where $E_{\rm p}$ is the proton energy,
$s=2$ is the proton spectrum index (Blandford & Eichler, 1987; Malkov & Drury,
2001). The accelerated protons will interact with the shocked gas ($pp$
process) or the thermal photons produced by the shocked disk gas ($p\gamma$
process) and thus produce high-energy neutrino. Since the energy of secondary
neutrinos comes from the parent protons, we can calculate the neutrino flux by
(Murase et al., 2016)

$\epsilon^{2}_{\nu}F_{\nu}=\frac{3K}{4(1+K)}\frac{\zeta_{{\rm p}}L_{{\rm
cr}}}{4\pi D^{2}_{{\rm L}}{\rm ln}(E_{\rm p,max}/E_{\rm p,min})},$ (9)

where $\epsilon_{\nu}=0.05E_{\rm p}$ is the neutrino energy, $D_{\rm L}$ is
the the luminosity distance, ${\rm ln}(E_{\rm p,max}/E_{\rm p,min})$ is the
normalized factor. $K=2$ for $t^{-1}_{\rm pp}>t^{-1}_{\rm p\gamma}$, $K=1$ for
$t^{-1}_{\rm pp}<t^{-1}_{\rm p\gamma}$, where $t_{\rm pp}$ and $t_{\rm
p\gamma}$ are timescales for $pp$ process and $p\gamma$ process. The time
scale for $pp$ process is

$t^{-1}_{\rm pp}=0.5\sigma_{\rm pp}n_{\rm s}c,$ (10)

where $\sigma_{\rm pp}\sim 5\times 10^{-26}cm^{-2}$ is the cross section for
$pp$ process, $n_{s}=4n_{\rm disk}=4\rho_{\rm disk}/m_{\rm p}$ is the number
density of shocked disk gas (Stecker, 1968; Murase, 2007).The time scale for
$p\gamma$ process is

$t^{-1}_{\rm p\gamma}=\frac{c}{2\gamma^{2}_{\rm p}}\int_{\tilde{E}_{\rm
th}}^{\infty}d\tilde{E}\sigma_{\rm p\gamma}(\tilde{E})\kappa_{\rm
p\gamma}(\tilde{E})\tilde{E}\int_{\tilde{E}/2\gamma_{\rm
p}}^{\infty}dE_{\gamma}E_{\gamma}^{-2}\frac{dN_{\gamma}}{dE_{\gamma}},$ (11)

where $\sigma_{\rm p\gamma}$ is the cross section for $p\gamma$ process
(Patrignani et al., 2016), $\tilde{E}$ is the energy of the photon in the rest
frame of the proton, $\tilde{E}_{\rm th}\simeq 145\rm MeV$ is the threshold
energy, $\gamma_{\rm p}=E_{\rm p}/m_{\rm p}c^{2}$ is the Lorentz factor of the
proton, and $dN_{\gamma}/dE_{\gamma}$ is the photon spectrum.

$\zeta_{\rm p}$ in Eq.(9) represents the efficiency of $pp$ and $p\gamma$
reactions (Murase, 2008), which can be estimated by $\zeta_{\rm
p}=(t^{-1}_{\rm pp}+t^{-1}_{\rm p\gamma})/(t^{-1}_{\rm cool})$, where
$t^{-1}_{\rm cool}=t^{-1}_{\rm pp}+t^{-1}_{\rm p\gamma}+t^{-1}_{\rm
BH}+t^{-1}_{\rm syn}+t^{-1}_{\rm IC}+t^{-1}_{\rm adi}$ is the total cooling
timescale for protons. The accelerated protons will lose energy via different
channels besides $pp$ or $p\gamma$ processes, including adiabatic cooling with
timescale $t_{\rm adi}=R_{\rm s}/v_{\rm s}$, synchrotron radiation cooling
with timescale

$t_{\rm syn}=\frac{6\pi m^{4}_{\rm p}c^{3}}{\sigma_{\rm T}m_{\rm e}B^{2}E_{\rm
p}},$ (12)

where the magnetic field strength can be estimated by
$B=\sqrt{8\pi\epsilon_{\rm B}u_{\rm s}}=\sqrt{12\pi\epsilon_{\rm B}n_{\rm
s}k_{\rm b}T_{\rm s}}$, $\epsilon_{\rm B}=0.1$. The timescale of inverse-
Compton scattering is (Denton & Tamborra, 2018):

$t_{\rm
IC}=\begin{cases}\frac{3m_{p}^{4}c^{3}}{4\sigma_{T}m_{e}^{2}n_{\gamma}E_{\gamma}E_{\rm
p}},&E_{\gamma}E_{p}<m_{p}^{2}c^{4},\\\ \\\
\frac{3E_{\gamma}E_{p}}{4\sigma_{T}m_{e}^{2}n_{\gamma}c^{5}}\,\,,&E_{\gamma}E_{p}>m_{p}^{2}c^{4}.\\\
\end{cases}$ (13)

where $n_{\gamma}$ is the number density of the thermal photon produced by
shocked disk gas, it can be estimated by $n_{\gamma}=\epsilon_{\rm e}u_{\rm
s}/E_{\gamma}$, where $E_{\gamma}=2.7k_{\rm b}T_{\rm s}$ is the typical
thermal photon energy, $\epsilon_{e}=0.1$. Also, protons may take the
Bethe–Heitler process with thermal photons, which timescale is (Denton &
Tamborra, 2018),

$t_{\rm BH}=\frac{E_{\rm p}\sqrt{m^{2}_{p}c^{4}+2E_{\rm
p}E_{\rm\gamma}}}{2n_{\gamma}\sigma_{\rm BH}m_{\rm e}c^{3}(E_{\rm
p}+E_{\gamma})}.$ (14)

The cross-section of the Bethe-Heitler process can be written as

$\sigma_{\rm BH}=\alpha r^{2}_{\rm e}[\frac{28}{9}{\rm ln}(\frac{2E_{\rm
p}E_{\gamma}}{m_{\rm p}m_{\rm e}c^{4}})-\frac{106}{9}],$ (15)

where $\alpha$ is the fine structure constant, $r_{\rm e}$ is the classical
electron radius.

Figure 4: The all flavor neutrino fluence for an individual event occurring at
$D_{\rm L}=10~{}\rm Mpc$ for the diverse SMBH masses ($M_{\rm SMBH}$) and
locations of merged BH from the SMBH (r). $\eta_{\rm w}$ is fixed as 0.1.

For the normalized factor in Eq.(9), we set $E_{\rm p,min}=1\,\rm TeV$ and the
$E_{\rm p,max}$ can be obtained by equating the accelerated timescale
$t_{acc}=20E_{\rm p}c/3eBv_{\rm s}$ with the minimum cooling timescale, where
$e$ is the charge of the electron. In Fig. 3, we show the timescales for
various processes when the forward shock expands at $R_{\rm s}=1\times
10^{14}\rm cm$ and $R_{\rm s}=3\times 10^{14}\rm cm$ (before and beyond
$R_{\rm br}$ ) in the case of $M_{\rm SMBH}=10^{7}M_{\odot},r=10^{3}R_{\rm
sch}$. We notice that the $pp$ process is generally more efficient than the
$p\gamma$ process. For all considered cases, we summarize that the maximum
proton energy after the shock breakout is in the range of 10 TeV to 20 PeV.

##  4 Neutrino fluence, light curves, and Diffuse Neutrino Emission

###  4.1 Individual Fluences, Light curves, and Expected Numbers

The neutrino fluence for an individual event can be obtained by

$\epsilon^{2}_{\nu}\phi_{\nu}\simeq\int_{t_{\rm
br}}^{\infty}\epsilon^{2}_{\nu}F_{\nu}dt,$ (16)

where $\epsilon^{2}_{\nu}F_{\nu}$ can be obtained through Eq.(9). In Fig. 4,
we show the all-flavor neutrino fluence under various conditions, where the
luminosity distance of the SMBH is $D_{\rm L}=10~{}\rm Mpc$ and $\eta_{\rm
w}=0.1$. Based on the effective area $A_{\rm eff}(\epsilon_{\nu})$ of the
IceCube detector, we also estimate the expected neutrino detection number
(Aartsen et al., 2020),

$N_{\nu}=\int A_{\rm eff}(\epsilon_{\nu})\phi_{\nu}d\epsilon_{\nu},$ (17)

and present them in Table 1.

Table 1: Number of detected $\nu_{\mu}$ from individual event at 10 Mpc by IceCube. Models | Detected Number
---|---
$10^{8}M_{\odot},10^{3}R_{sch}$ | $0.137$
$10^{8}M_{\odot},10^{4}R_{sch}$ | $0.687$
$10^{7}M_{\odot},10^{3}R_{sch}$ | $0.01$
$10^{7}M_{\odot},10^{4}R_{sch}$ | $0.163$
$10^{6}M_{\odot},10^{3}R_{sch}$ | $0.0007$
$10^{6}M_{\odot},10^{4}R_{sch}$ | $0.017$
Figure 5: The contribution to the diffuse neutrino background by merged BH
winds with different BBH merging rate densities, where $r=10^{3}R_{\rm sch}$,
$\eta_{\rm w}=0.1$. The black data points are observed astrophysical diffuse
neutrino flux measured by the IceCube detector (Aartsen et al., 2015).

From Fig.4, we can see that neutrino fluence roughly increases with the mass
of SMBH and radial location $r$. This is because the forward shock can exist
longer after shock breakout when $M_{\rm SMBH}$ or $r$ are larger. The shock
can carry more energy, and neutrino production can be more efficient, leading
to higher neutrino fluence. This feature is also illustrated by the light
curves shown in Fig. 6. One can see from the light curves that the neutrino
burst persists for several hours when $M_{\rm SMBH}$ or $r$ is small, while it
can persist for even more than 400 hours when $M_{\rm SMBH}$ or $r$ is large.
We can also find that the neutrino cutoff energy is located at the lower
energy when the fluence is higher. In the case of $M_{\rm
SMBH}=10^{7}M_{\odot},r=10^{3}R_{\rm sch}$, the neutrino cutoff energy is
located at several hundred TeV, while for $M_{\rm
SMBH}=10^{8}M_{\odot},r=10^{4}R_{\rm sch}$ it is at several ten GeV. The
maximum proton energy will be suppressed for a more massive SMBH or a further
location to SMBH because of more significant cooling processes for protons,
inducing a smaller maximum proton energy and a smaller neutrino cutoff energy.

###  4.2  Diffuse Neutrino Emission

The diffuse neutrino emission can be obtained through (Razzaque et al., 2004;
Xiao et al., 2016)

$E^{2}_{\nu}\Phi_{\nu}(E_{\nu})=\frac{c}{4\pi H_{0}}\int_{0}^{z_{{\rm
max}}}{\frac{\rho_{0}f(z)\epsilon^{2}_{\nu}\phi_{\nu}(\epsilon_{\nu})}{(1+z)\sqrt{\Omega_{{\rm
m}}(1+z)^{3}+\Omega_{\Lambda}}}}dz,$ (18)

where $z$ is the redshift, $E_{\nu}=\epsilon_{\nu}/(1+z)$ is the observed
neutrino energy, $H_{0}=67.8~{}\rm km/s$ is the Hubble constant, $\Omega_{\rm
m}=0.308$, $\Omega_{\Lambda}=0.692$ are the cosmology constants (Ade et al.,
2016). $f(z)$ is the redshift evolution factor and is roughly taken as the
same as for short gamma-ray bursts because both are binary compact object
mergers, that is,
$f(z)=[(1+z)^{5.7\eta}+(\frac{1+z}{0.36})^{1.3\eta}+(\frac{1+z}{3.3})^{-9.5\eta}+(\frac{1+z}{3.3})^{-24.5\eta}]^{1/\eta}$
with $\eta=-2$ (Sun et al., 2015). The density of BBH merging rate within AGN
disks $\rho_{0}$ is still quite uncertain, and we adopt
$\rho_{0}\simeq[0.002,20]~{}\rm Gpc^{-3}yr^{-1}$ proposed by Gröbner et al.
(2020). The diffuse neutrino emissions produced by the merged BH wind are
shown in Fig. 5. We find that neutrino production from the remnant BH wind has
little contribution to the diffuse neutrino background due to the low event
rate density, but the contribution can be improved if the larger $\eta_{\rm
w}$ and location $r$ from SMBH are involved. In addition, the neutrino from
the remnant BH wind can still be identified from the near individual source as
shown in Table 1.

Figure 6: Top panel: The light curves for neutrino bursts under various
conditions. The vertical dotted line represents the shock breakout time (as
the initial time). Bottom panel: Schematic light curve for the neutrino/EM and
GW joint observation. There might be twice neutrino/EM bursts in the premerger
phase and postmerger phase, respectively. The shaded region represents the
gravitational wave duration time.

###  4.3 Multimessenger Observations

Although the diffuse neutrino contribution from merged BH winds is relatively
low, the neutrino fluence for an individual nearby event is still
considerable, even compared to that of other neutrino-produced activities on
the disk, such as supernova explosions (Zhu et al., 2021; Zhou et al., 2023),
choked gammar-ray bursts (Zhu et al., 2021), and jet powered by the merged BH
(Zhou & Wang, 2023).

Besides, there would be a detectable electromagnetic signal (thermal emission
of the shocked gas) after the shock breakout, which luminosity can be
estimated by $L_{\rm br}\simeq 7.1\times 10^{45}~{}{\rm erg/s}(\frac{\rho_{\rm
disk}}{\rm 10^{-9}g~{}cm^{-}{3}})(\frac{H}{10^{13.5}{\rm
cm}})^{2}(\frac{v_{\rm s}}{10^{9}{\rm cm~{}s^{-1}}})^{3}$ (Kimura et al.,
2021). Such a signal may be distinguished from the background AGN radiation
becoming an EM transient with the neutrino signal simultaneously. The BBH
merger would produce detectable GW signals for Ligo/Virgo, thus we can expect
a joint light curve for GW, EM, and neutrino, as illustrated in Fig. 6. Notice
that there might be twice neutrino/EM bursts in the premerger and postmerger
phases, respectively. The former is produced by the BBH wind, which can be
modeled by the same method as in Sections 2 and 3. After the BBH merger, the
kicked newborn remnant BH can spend a time of $t_{\rm cross}\simeq R_{\rm
cav}/v_{\rm kick}$ to get out of the cavity formed by premerger BBH winds,
where $R_{\rm cav}\sim R_{\rm Hill}$ is the radius of the cavity. Therefore,
there will be a time delay between the GW signal from the BBH merger and the
second EM/neutrino transient, i.e., $t_{\rm delay}=t_{\rm cross}+t_{\rm br}$.
In the case of $M_{\rm SMBH}=10^{7}M_{\odot},r=10^{3}R_{\rm sch}$, $v_{\rm
kick}=10^{7}{\rm cm/s}$, the gap time is $\sim 60$ days.

##  5 Summary

The BBH within the AGN disks would merge and leave the remnant BH. The
premerger outflows can create a cavity around the merged BH. When the newborn
merged BH is kicked out of the cavity into a gas-abundant disk environment, it
will accrete the disk gas in the way of the BHL accretion and power a fast
wind. Such fast outflow will drive a forward shock in the disk gas and thus
accelerate the protons and electrons. The accelerated protons will
subsequently interact with photons/matters of the shocked gas to produce high-
energy neutrinos. In this paper, we investigate the detailed dynamic evolution
of the interaction between the remnant BH wind and the disk medium, calculate
the neutrino production and expected detection numbers for an individual
event, and estimate the contribution of these events to the diffuse neutrino
background.

The stellar-mass BBH merger events have become the most important detection
targets for ground-based GW detectors including Ligo, Virgo and KAGRA (Abbott
et al., 2019). In addition to GW bursts, BBH mergers embedded in AGN disks are
expected to produce potential EM and neutrino signals, and thus the
multimessenger study for BBH mergers in AGN disks is very desirable.
Especially, high-energy neutrinos will be an ideal and unique messenger for
studying BBH mergers within AGN disks, as they can easily escape from the
complicated environments of AGNs. Recently, Abbasi et al. (2023) performed the
high-energy neutrino searches for GW events and no significant neutrino
emission was found, which is consistent with our results for their far
luminosity distances and low expectations of neutrino fluxes. The neutrino
follow-up search on the candidate optical counterpart of the BBH merger event
GW190521 was carried out as well. However, the relatively far distance of
GW190521, i.e., $3931\pm 953\,\rm Mpc$ (Graham et al., 2020), makes the
expected neutrino fluence low and can not provide the useful constraints
currently. Future GW detectors will further significantly increase the
detection number of BBH mergers, which is expected to reach
$10^{5}-10^{6}\,\rm yr^{-1}$ (Maggiore et al., 2020), making them one of the
most important astrophysical sources. The association between GWs and high-
energy neutrinos can be expected.

## Acknowledgments

We acknowledge support from the National Natural Science Foundation of China
under grant No.12003007 and the Fundamental Research Funds for the Central
Universities (No. 2020kfyXJJS039).

##  Appendix A The modified $\alpha$-disk model

To describe the midplane of the disk, a modified $\alpha$-disk model is
adopted (Cantiello et al., 2021) and we assume that all processes occur in the
gas-pressure-dominated region. By considering the gravitational stability, the
disk is separated into two regions by a characteristic radius,

$R_{Q}\simeq
0.02\left(\frac{\kappa^{3}\alpha^{7}}{\mu^{12}}\right)^{2/27}\frac{M^{1/27}_{8}}{f^{8/27}_{m}}\
{\rm pc},$ (A1)

where $\kappa\sim 1$ is the opacity, $\alpha\sim 0.1$ is the viscosity factor
and $\mu\sim 0.5$ is the molecular weight. $M_{8}=M_{\rm
SMBH}/10^{8}M_{\odot}$ with the mass of central SMBH $M_{\rm SMBH}$ and the
mass of the Sun $M_{\odot}$. At a distance $r\leq R_{Q}$ from the SMBH, the
Toomre parameter Q satisfies $Q\gtrsim 1$, suggesting that the disk is
gravitationally stable. In such region, the midplane gas temperature and
density can be described by

$T_{\rm mid}=3.7\times
10^{3}\left(\frac{\kappa\mu}{\alpha}\right)^{1/5}\frac{f^{2/5}_{m}M^{7/10}_{8}}{r^{9/10}_{\rm
pc}}\ {\rm K},$ (A2) $\rho_{\rm mid}=1.5\times
10^{-13}\frac{\mu^{6/5}f^{2/5}_{m}M^{19/20}_{8}}{\kappa^{3/10}\alpha^{7/10}r^{33/20}_{\rm
pc}}\ {\rm g/cm^{3}},$ (A3)

respectively, where $r_{\rm pc}=r/1\ {\rm pc}$, $f_{m}\sim 1$ is the
dispersion factor. We also have the scale height for the disk,

$H\simeq\sqrt{\frac{2R_{g}T_{\rm mid}}{\mu\Omega}},$ (A4)

where $R_{g}$ is the special gas constant defined as $R_{g}=k_{b}/m_{\rm p}$,
$\Omega=\sqrt{GM_{\rm SMBH}/r^{3}}$ is the angular frequency at a distance r
from the SMBH and surface density is

$\Sigma=2\rho_{\rm mid}H,$ (A5)

At the region $r>R_{Q}$, the conventional $\alpha$-disk model should be
modified by considering the gravity normal to the disk plane, momentum
transfer, and possible thermal and momentum feedback from newly born/captured
stars within the disk. The gravitational stability will be maintained with
$Q\sim 1$ by assuming a self-regulated star formation rate, which leads to the
modified formulas at $r>R_{Q}$,

$T^{\prime}_{\rm mid}=1.3\times 10^{4}M^{2/3}_{8}\left(\frac{f_{\rm
m}Q}{\alpha}\right)^{2/3}\ {\rm K},$ (A6) $\rho^{\prime}_{\rm mid}=8.3\times
10^{-16}\frac{M_{8}}{r^{3}_{\rm pc}Q}\ {\rm g/cm^{3}},$ (A7)

Also, the scale height and surface density are

$H^{\prime}=\frac{0.025r^{1/2}_{\rm pc}f^{1/3}_{\rm
m}Q^{1/3}r}{\alpha^{1/3}M^{1/6}_{8}}\ {\rm cm},$ (A8)
$\Sigma^{\prime}=\frac{180M^{5/6}_{8}f^{1/3}_{\rm m}}{r^{2/3}_{\rm
pc}Q^{2/3}\alpha^{1/3}}\ {\rm g/cm^{3}},$ (A9)

When the location of the merged black hole is in the region of $r\leq R_{Q}$,
we use the formulas Eq.(A2)–Eq.(A4) to describe the midplane properties of the
disk; otherwise, the formulas Eq.(A6)–Eq.(A9) are adopted.

##  Appendix B simplication for dynamic equations

In this section, we give the original dynamical equations for the shocked wind
region (Region (b)), our model is similar to that of Liu et al. (2018),Wang &
Loeb (2015) and Wang & Loeb (2016). The electrons, protons, and photons in
Region (b) will interact with each other through various channels. First,
protons and electrons will lose energy through free–free emission and
synchrotron emission. In addition, protons and electrons may collide with each
other to exchange energy, namely the Coulomb collision. Finally, they will
change energy with photons through Compton scattering. By considering these
particle interactions, the internal energy rates for protons, electrons, and
photons are described by

$\frac{dE_{\rm p}}{dt}=\frac{1}{2}\dot{M}_{\rm w}(1-\frac{v_{\rm rs}}{v_{\rm
w}})(v^{2}_{\rm w}-v^{\prime}_{\rm sw})-4\pi R^{2}_{\rm s}P_{\rm psw}v_{\rm
s}-L_{\rm cou}-L_{\rm rad,p}-L_{\rm com,p}$ (B1) $\frac{dE_{\rm
e}}{dt}=\frac{1}{2}\dot{M}_{\rm w}(\frac{m_{\rm e}}{m_{\rm p}})(1-\frac{v_{\rm
rs}}{v_{\rm w}})(v^{2}_{\rm w}-v^{\prime}_{\rm sw})-4\pi R^{2}_{\rm s}P_{\rm
esw}v_{\rm s}+L_{\rm cou}-L_{\rm rad,e}-L_{\rm com,e}$ (B2) $\frac{dE_{\rm
ph}}{dt}=L_{\rm rad,p}+L_{\rm rad,e}+L_{\rm com,p}+L_{\rm com,e}-4\pi
R^{2}_{\rm s}P_{\rm ph}v_{\rm s},$ (B3)

where the first terms on the right hand of Eq.(B1) and Eq.(B2) are the
injected thermal energy rate for protons and electrons. $L_{\mathrm{rad},i}$
and $L_{\mathrm{com},i}$ ($i=e,p$) are the radiation cooling rate and Compton
cooling/heating rate for electrons and protons (Sazonov & Sunyaev, 2001;
Faucher-Giguère & Quataert, 2012), while they are energy injection term for
photons. $L_{\rm cou}$ is the energy changing rate through Coulomb collision.
$4\pi R^{2}_{\rm s}P_{\rm j}v_{\rm s}$ (j=psw,esw,ph) is the expanding cooling
term for protons, electrons, and photons. We also have the energy-pressure
relations,

$E_{\rm p/e}=\frac{3}{2}P_{\rm p/e}V_{\rm sw}$ (B4)

for gas, which comes from the ideal gas relation $P_{g}=nk_{\rm
b}T=\frac{2}{3}u_{\rm g}$ with the energy density of the gas $u_{\rm g}$.
Also, we have

$E_{\rm ph}=3P_{\rm ph}V_{\rm sw}$ (B5)

for photons, which comes from the relation $P_{\rm ph}=\frac{1}{3}u_{\rm ph}$
for uniform photon fields.

In principle, one can solve the above equations together with Eq.(6) to obtain
the thermal energy/pressure evolution of three species of particles. However,
we only care about the total pressure $P_{\rm tot}=P_{\rm psw}+P_{\rm
esw}+P_{\rm ph}$. Therefore, we can simplify the equations Eq.(B1)–Eq.(B3) by
adding them together. Then all the energy interaction terms will vanish and
the injected power of protons will be much larger than that of electrons
($m_{\rm e}/m_{\rm p}\simeq 1/1000$), which gives the following

$\frac{dE_{\rm tot}}{dt}=\frac{1}{2}\dot{M}_{w}(1-\frac{v_{\rm rs}}{v_{\rm
w}})(v^{2}_{\rm w}-v^{\prime}_{\rm sw})-4\pi R^{2}_{\rm s}P_{\rm tot}v_{\rm
s}.$ (B6)

It is just Eq.(7). For the energy-pressure relationship, we notice that the
right hand of Eq.(B4) and Eq.(B5) is nearly comparable with only twice the
difference. Namely, the pressure provided by gas/photon is comparable when
they have the same energy. Therefore, we uniformly use Eq.(B4) for
simplification.

## References

  * Aartsen et al. (2015) Aartsen, M., Ackermann, M., Adams, J., et al. 2015, The Astrophysical Journal Letters, 805, L5
  * Aartsen et al. (2020) —. 2020, Physical review letters, 124, 051103
  * Abbasi et al. (2023) Abbasi, R., Ackermann, M., Adams, J., et al. 2023, The Astrophysical Journal, 944, 80, doi: 10.3847/1538-4357/aca5fc
  * Abbott et al. (2019) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, Phys. Rev. X, 9, 031040, doi: 10.1103/PhysRevX.9.031040
  * Abbott et al. (2020) Abbott, R., Abbott, T., Abraham, S., et al. 2020, Physical review letters, 125, 101102
  * Ade et al. (2016) Ade, P. A., Aghanim, N., Arnaud, M., et al. 2016, Astronomy & Astrophysics, 594, A13
  * Artymowicz et al. (1993) Artymowicz, P., Lin, D., & Wampler, E. 1993, Astrophysical Journal, Part 1 (ISSN 0004-637X), vol. 409, no. 2, p. 592-603., 409, 592
  * Bartos et al. (2017) Bartos, I., Kocsis, B., Haiman, Z., & Márka, S. 2017, The Astrophysical Journal, 835, 165
  * Bellovary et al. (2016) Bellovary, J. M., Mac Low, M.-M., McKernan, B., & Ford, K. S. 2016, The Astrophysical Journal Letters, 819, L17
  * Blandford & Eichler (1987) Blandford, R., & Eichler, D. 1987, Physics Reports, 154, 1
  * Bondi (1952) Bondi, H. 1952, Monthly Notices of the Royal Astronomical Society, 112, 195
  * Campanelli et al. (2007) Campanelli, M., Lousto, C. O., Zlochower, Y., & Merritt, D. 2007, Physical Review Letters, 98, 231102
  * Cantiello et al. (2021) Cantiello, M., Jermyn, A. S., & Lin, D. N. 2021, The Astrophysical Journal, 910, 94
  * Caprioli & Spitkovsky (2014) Caprioli, D., & Spitkovsky, A. 2014, The Astrophysical Journal, 783, 91
  * Chen & Dai (2024) Chen, K., & Dai, Z.-G. 2024, ApJ, 961, 206, doi: 10.3847/1538-4357/ad0dfd
  * Cheng & Wang (1999) Cheng, K., & Wang, J.-M. 1999, The Astrophysical Journal, 521, 502
  * Collin & Zahn (1999) Collin, S., & Zahn, J.-P. 1999, Ap&SS, 265, 501, doi: 10.1023/A:1002191506811
  * Denton & Tamborra (2018) Denton, P. B., & Tamborra, I. 2018, The Astrophysical Journal, 855, 37
  * Dittmann et al. (2021) Dittmann, A. J., Cantiello, M., & Jermyn, A. S. 2021, The Astrophysical Journal, 916, 48
  * Edgar (2004) Edgar, R. 2004, New Astronomy Reviews, 48, 843
  * Fabj et al. (2020) Fabj, G., Nasim, S. S., Caban, F., et al. 2020, Monthly Notices of the Royal Astronomical Society, 499, 2608
  * Faucher-Giguère & Quataert (2012) Faucher-Giguère, C.-A., & Quataert, E. 2012, Monthly Notices of the Royal Astronomical Society, 425, 605
  * Gonzalez et al. (2007) Gonzalez, J. A., Sperhake, U., Bruegmann, B., Hannam, M., & Husa, S. 2007, Physical Review Letters, 98, 091101
  * Graham et al. (2020) Graham, M., Ford, K., McKernan, B., et al. 2020, Physical review letters, 124, 251102
  * Graham et al. (2023) Graham, M. J., McKernan, B., Ford, K. S., et al. 2023, The Astrophysical Journal, 942, 99
  * Gröbner et al. (2020) Gröbner, M., Ishibashi, W., Tiwari, S., Haney, M., & Jetzer, P. 2020, Astronomy & Astrophysics, 638, A119
  * Herrmann et al. (2007) Herrmann, F., Hinder, I., Shoemaker, D., & Laguna, P. 2007, Classical and Quantum Gravity, 24, S33
  * Hoyle & Lyttleton (1939) Hoyle, F., & Lyttleton, R. A. 1939in , Cambridge University Press, 405–415
  * Jiang et al. (2014) Jiang, Y.-F., Stone, J. M., & Davis, S. W. 2014, The Astrophysical Journal, 796, 106
  * Kimura et al. (2021) Kimura, S. S., Murase, K., & Bartos, I. 2021, The Astrophysical Journal, 916, 111
  * Leigh et al. (2018) Leigh, N. W., Geller, A. M., McKernan, B., et al. 2018, Monthly Notices of the Royal Astronomical Society, 474, 5672
  * Liu et al. (2018) Liu, R.-Y., Murase, K., Inoue, S., Ge, C., & Wang, X.-Y. 2018, The Astrophysical Journal, 858, 9
  * Maggiore et al. (2020) Maggiore, M., Van Den Broeck, C., Bartolo, N., et al. 2020, J. Cosmology Astropart. Phys, 2020, 050, doi: 10.1088/1475-7516/2020/03/050
  * Malkov & Drury (2001) Malkov, M., & Drury, L. O. 2001, Reports on Progress in Physics, 64, 429
  * McKernan et al. (2020) McKernan, B., Ford, K., & O’Shaughnessy, R. 2020, Monthly Notices of the Royal Astronomical Society, 498, 4088
  * Murase (2007) Murase, K. 2007, Physical Review D, 76, 123001
  * Murase (2008) —. 2008, Physical Review D, 78, 101302
  * Murase et al. (2016) Murase, K., Guetta, D., & Ahlers, M. 2016, Physical Review Letters, 116, 071101
  * Ohsuga et al. (2005) Ohsuga, K., Mori, M., Nakamoto, T., & Mineshige, S. 2005, The Astrophysical Journal, 628, 368
  * Patrignani et al. (2016) Patrignani, C., Agashe, K., Aielli, G., et al. 2016
  * Razzaque et al. (2004) Razzaque, S., Mészáros, P., & Waxman, E. 2004, Physical Review Letters, 93, 181101
  * Sazonov & Sunyaev (2001) Sazonov, S. Y., & Sunyaev, R. 2001, Astronomy Letters, 27, 481
  * Skadowski et al. (2013) Skadowski, A., Narayan, R., Penna, R., & Zhu, Y. 2013, Monthly Notices of the Royal Astronomical Society, 436, 3856
  * Stecker (1968) Stecker, F. 1968, Physical Review Letters, 21, 1016
  * Sun et al. (2015) Sun, H., Zhang, B., & Li, Z. 2015, The Astrophysical Journal, 812, 33
  * Tagawa et al. (2020) Tagawa, H., Haiman, Z., & Kocsis, B. 2020, The Astrophysical Journal, 898, 25
  * Tagawa et al. (2023) Tagawa, H., Kimura, S. S., & Haiman, Z. 2023, The Astrophysical Journal, 955, 23, doi: 10.3847/1538-4357/ace71d
  * Tagawa et al. (2022) Tagawa, H., Kimura, S. S., Haiman, Z., et al. 2022, The Astrophysical Journal, 927, 41
  * Wang et al. (2021) Wang, J.-M., Liu, J.-R., Ho, L. C., Li, Y.-R., & Du, P. 2021, The Astrophysical Journal Letters, 916, L17
  * Wang et al. (2011) Wang, J.-M., Ge, J.-Q., Hu, C., et al. 2011, The Astrophysical Journal, 739, 3
  * Wang et al. (2019) Wang, K., Huang, T.-Q., & Li, Z. 2019, The Astrophysical Journal, 872, 157
  * Wang & Loeb (2015) Wang, X., & Loeb, A. 2015, Monthly Notices of the Royal Astronomical Society, 453, 837
  * Wang & Loeb (2016) —. 2016, Journal of Cosmology and Astroparticle Physics, 2016, 012
  * Weaver et al. (1977) Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, Astrophysical Journal, Part 1, vol. 218, Dec. 1, 1977, p. 377-395., 218, 377
  * Xiao et al. (2016) Xiao, D., Mészáros, P., Murase, K., & Dai, Z.-G. 2016, The Astrophysical Journal, 832, 20
  * Zhou & Wang (2023) Zhou, Z.-H., & Wang, K. 2023, ApJ, 958, L12, doi: 10.3847/2041-8213/ad096f
  * Zhou et al. (2023) Zhou, Z.-H., Zhu, J.-P., & Wang, K. 2023, The Astrophysical Journal, 951, 74
  * Zhu et al. (2021) Zhu, J.-P., Wang, K., & Zhang, B. 2021, ApJ, 917, L28
  * Zhu et al. (2021) Zhu, J.-P., Wang, K., Zhang, B., et al. 2021, The Astrophysical Journal Letters, 911, L19


Rule 60 cellular automaton, Mersenne numbers, and the Newman-Moore spin
lattice

Jonás Carmona-Pírez1, Adrian J. Peguero2, Vanja Dunjko2, Maxim Olshanii2,
Joanna Ruhl2*

1 Technical Advisory Subdirectorate of Information Management (STAGI),
Andalusian Health Service (SAS), 41071 Seville, Spain

2 Department of Physics, University of Massachusetts Boston, Boston
Massachusetts 02125, USA

*<EMAIL_ADDRESS>
## Abstract

The goal of this paper is to review the properties of a Rule 60 cellular
automaton on a ring with a Mersenne number circumference and to use this
knowledge to explicitly construct all the ground state configurations of the
classical Newman-Moore model (a particular two-dimensional spin lattice model
with a specific three-spin interaction) on a square lattice of the same size.
In this particular case, the number of ground states is equal to half of the
available spin configurations in any given row of the lattice.

###### Contents

  1. 1 Introduction
  2. 2 Rule 60 cellular automaton
    1. 2.1 Periodic trajectories
    2. 2.2 Configurations and their characterization
    3. 2.3 Useful lemmas
  3. 3 Fundamental periods and a propagator for the Rule 60 on a ring of a Mersenne number circumference
    1. 3.1 Fundamental periods
    2. 3.2 Propagator
  4. 4 Application to the Newman-Moore model
    1. 4.1 Connection between the Rule 60 automaton and the Newman-Moore model
    2. 4.2 Newman-Moore results inferred from the properties of the Rule 60 automaton
  5. 5 Summary of results

##  1 Introduction

Both the classical Newman-Moore model (a particular two-dimensional spin
lattice model with a specific three-spin interaction) [1] and it’s quantum
counterpart [2] play an important role in studies of systems with glassy
dynamics (see, for example,[3], [4], [5] and references therein). In this
paper, we review the properties of a Rule 60 cellular automaton on a ring with
Mersenne number circumference and use this knowledge to explicitly construct
all ground state configurations of the classical Newman-Moore model on a
square lattice of the same size.

##  2 Rule 60 cellular automaton

For the one-dimensional cellular automata [6], the state of the system at the
instant of time $j$ is given by a length $L$ string:

$\displaystyle\vec{x}_{j}=[x_{0,j},\,x_{1,j},\,\ldots\,\,x_{i,j},\,\ldots,\,x_{L-1,j}]\,\,,$
(1)

where $x_{i,j}$ is the state of the $i$’th bit at the instance of time $j$.
Each $x_{i,j}$ is either $0$, if the site is unoccupied, or $1$, if the site
is occupied. Consider a Rule 60 automaton [6] on a ring of a circumference
$L$. According to the rule, the state of the automaton at the instance of time
$j+1$ is related to its state at time $j$ by

$\displaystyle\begin{array}[]{lcc}x_{i,j+1}=\text{Xor}(x_{i-1,j},\,x_{i,j})&\text{
for }&i=1,\,2,\,\ldots,\,L-1\\\ x_{0,j+1}=\text{Xor}(x_{L-1,j},\,x_{0,j})\\\
j=0,\,1,\,\ldots\quad.\end{array}$ (5)

###  2.1 Periodic trajectories

In what follows, we will be interested in periodic trajectories generated by
the rule (5).

######  Definition 2.1.

An initial condition $\vec{x}_{0}$ initiates a $\bm{M}$-periodic trajectory
$\vec{x}_{j}$ if

$\displaystyle\vec{x}_{M}=\vec{x}_{0}\,\,.$

######  Definition 2.2.

An initial condition $\vec{x}_{0}$ initiates a $\bm{M}$-fundamental-periodic
trajectory $\vec{x}_{j}$. if

$\displaystyle\vec{x}_{M}=\vec{x}_{0}$ and $\displaystyle\text{for any
}0<j<M,\quad\vec{x}_{j}\neq\vec{x}_{0}\,\,.$

###  2.2 Configurations and their characterization

The state of the system at a given instant of time is called a configuration.
E.g. for the trajectory

$\displaystyle\begin{array}[]{lll}j=0&&[100]\\\ j=1&&[110]\\\ j=2&&[011]\\\
j=3&&[101]\\\ &\ldots&\end{array}\,\,,$

$[011]$ is the state of the size $L=3$ lattice at the instance of time $j=2$,
given that the initial condition at $j=0$ was $[100]$.

######  Definition 2.3.

Let ${}_{L}N_{\text{occupied}}$-type state indicate a configuration that
features $N_{\text{occupied}}$ occupied sites total on a length $L$ lattice.

E.g. $[110]$ will be characterised as a configuration of a ${}_{3}2$-type.

###  2.3 Useful lemmas

######  Lemma 2.1.

For any $L$, there are as many ${}_{L}(\text{\rm even})$-type as there are
${}_{L}(\text{\rm odd})$-type states.

###### Proof.

One can establish a one-to-one correspondence between ${}_{L}(\text{\rm
even})$ and ${}_{L}(\text{\rm odd})$ states simply by controlling the state of
the first bit, leaving other bits intact. ∎

######  Lemma 2.2.

For any $L$, the state $[00\ldots 0]$ is a stationary state (a $1$-periodic
trajectory).

###### Proof.

This directly follows from the rule (5). ∎

######  Lemma 2.3.

For any ring of circumference $L$, both ${}_{L}(\text{\rm odd})$-type and
${}_{L}(\text{\rm even})$-type states become an ${}_{L}(\text{\rm even})$-type
state after one time step:

${}_{L}(\text{\rm odd})\to_{L}(\text{\rm even})$ ${}_{L}(\text{\rm
even})\to_{L}(\text{\rm even})\,\,.$

###### Proof.

One proof is presented in [6] as Lemma 3.1. We present an alternative proof
here. Assume that the state of the system at an instance $j$ is

$[x_{0,j},\,x_{1,j},\,\ldots\,\,x_{i,j},\,\ldots,\,x_{L-1,j}]\,\,.$

According to the rule (5), at the next step, $1$’s will be marking points
where the site state changes from $0$ to $1$ or vice versa, along the lattice.
For example, in

$\displaystyle\begin{array}[]{lll}\vec{x}_{j}&=&[011\bm{01}1001]\\\
\vec{x}_{j+1}&=&[1101\bm{1}0101]\end{array}\,\,,$

the boldface “$1$” at $j+1$ marks a lateral switch from “$0$” to “$1$”, in
bold, at $j$. Likewise in

$\displaystyle\begin{array}[]{lll}\vec{x}_{j}&=&[\bm{0}1101100\bm{1}]\\\
\vec{x}_{j+1}&=&[\bm{1}10110101]\end{array}\,\,,$

the boldface “$1$” at $j+1$ marks a lateral switch from “$1$” to “$0$” (in
boldface), as the ring closes.

Observe now that on a ring, the site state can only change laterally an _even_
number of times. This proves the Lemma. ∎

######  Corollary 2.3.1.

Type ${}_{L}(\text{\rm odd})$ can not initiate a periodic motion.

###### Proof.

For an $M$-periodic motion, the state must return to itself after $M$ steps.
But according to the Lemma 2.3, an ${}_{L}(\text{\rm odd})$-type state can
never become an ${}_{L}(\text{\rm odd})$-type state. ∎

######  Lemma 2.4.

Let $M_{\text{\rm max}}$ be the maximum fundamental period for a given $L$.
Let $M_{\mathcal{P}}$ be the period which the ${}_{L}1$ system converges to.
Then

$M_{\text{\rm max}}=M_{\mathcal{P}}\,\,.$

Other periods divide $M_{\text{\rm max}}$.

###### Proof.

The Lemma 2.4 appears in [7] as Theorem 1. An alternative formulation related
to cyclotomic polynomials is given in [8], Proposition 3.6. ∎

An explicit form of a propagator, for $L=2^{n}-1$ ($n=2,\,3,\,\ldots$) is
given by equation 19 in Theorem 3.6.

######  Lemma 2.5.

Let k and n be positive integers. Then, the binomial coefficient
$\binom{2^{n}}{k}$ is even for all $0<k<2^{n}$.

Several proofs for this property are known, the following is perhaps the
simplest one, adapted from [9] with minor changes.

###### Proof.

According to the binomial theorem,

$(x+1)^{2^{n}}=\sum_{k=0}^{2^{n}}\binom{2^{n}}{k}x^{k}\,\,.$

Therefore the lemma will be established if we prove that for all positive
integer $n$ there is some integer $a_{k}$ with $0<k<2^{n}$ such that

$\displaystyle(x+1)^{2^{n}}=x^{2^{n}}+1+\sum_{k=1}^{2^{n}-1}2a_{k}\,x^{k}\,\,.$
(6)

We will prove this by induction on $n$. First we prove the base case where
$n=1$, so that

$(x+1)^{2^{n}}=(x+1)^{2}=x^{2}+1+2x\,\,.$

which has the form of 6 with $a_{1}=1$.

The induction hypothesis is that this property holds for $n=m$, for a given
$m\geq 1$, so that

$(x+1)^{2^{m}}=x^{2^{m}}+1+\sum_{k=1}^{2^{m}-1}2a_{k}\,x^{k}\,\,.$

We now prove the property holds for $m+1$:

$\displaystyle(x+1)^{2^{m+1}}$ $\displaystyle=\left((x+1)^{2^{m}}\right)^{2}$

By the induction hypothesis this is

$\displaystyle\left((x+1)^{2^{m}}\right)^{2}$
$\displaystyle=\left(x^{2^{m}}+1+\sum_{k=1}^{2^{m}-1}2a_{k}\,x^{k}\right)^{2}$
$\displaystyle=\left\\{\left(x^{2^{m}}+1\right)^{2}\right\\}+\left\\{2\left(x^{2^{m}}+1\right)\left(2\sum_{k=1}^{2^{m}-1}a_{k}\,x^{k}\right)\right\\}+\left\\{\left(2\sum_{k=1}^{2^{m}-1}a_{k}\,x^{k}\right)^{2}\right\\}$
$\displaystyle=x^{2^{m+1}}+1+2\,x^{2^{m}}$
$\displaystyle\quad+4\left\\{\left(x^{2^{m}}+1\right)\sum_{k=1}^{2^{m}-1}a_{k}\,x^{k}+\left(\sum_{k=1}^{2^{m}-1}a_{k}\,x^{k}\right)^{2}\right\\}$

The first two terms to the right of the final equality are the first two terms
of 6 with $n=m+1$. The binomial theorem applied to $(x+1)^{2^{m+1}}$ produces
these same two terms, and then a polynomial in $x$ with powers running from
$1$ to $2^{m+1}-1$, inclusive. Because all terms in this particular remainder
are multiplied either by two or by four, they are all explicitly even, and the
resulting polynomial will be sums of even numbers, which must therefore also
be even. So it is possible to find numbers $b_{k}$ such that

$(x+1)^{2^{m+1}}=x^{2^{m+1}}+1+\sum_{k=1}^{2^{m+1}-1}2b_{k}\,x^{k}\,\,.$

In particular,

$b_{k}=\begin{cases}2a_{1},&\text{for $k=1$;}\\\
2a_{k}+\sum_{m=1}^{k-1}a_{k-m}a_{m},&\text{for $k=2,\,\ldots,\,2^{m}-1$;}\\\
1+2\sum_{m=1}^{2^{m}-1}a_{2^{m}-m}a_{m},&\text{for $k=2^{m}$;}\\\
2a_{k-2^{m}}+\sum_{m=k-2^{m}+1}^{2^{m}-1}a_{m}a_{k-m},&\text{for
$k=2^{m}+1,\,\ldots,\,2^{m+1}-2$; and}\\\ 2a_{2^{m}-1},&\text{for
$k=2^{m+1}-1$.}\end{cases}$

∎

######  Lemma 2.6.

For the Rule 60 automaton on a ring of circumference $L=2^{n}$
($n=2,\,3,\,\ldots$), the only allowed fundamental period is $M=1$.

This result is well known [6, 5, 7]. However, in what follows, we will be
proving a statement concerning temporal periods on a $L=2^{n}-1$ ring, and we
are going to use the proof of the Lemma above as a model.

###### Proof.

[7] Lemma 1, Corollary 1; [5] Eq. 13. ∎

##  3 Fundamental periods and a propagator for the Rule 60 on a ring of a
Mersenne number circumference

###  3.1 Fundamental periods

######  Definition 3.1.

A domain is the maximal contiguous set of sites having the same occupation
state bounded by a right domain boundary and a left domain boundary.

######  Definition 3.2.

A site $x_{i,j}$ is said to be a right domain boundary (RDB) if the site
$x_{i+1,j}$ has a different occupation than $x_{i,j}$. Similarly, a site
$x_{i,j}$ is said to be a left domain boundary (LDB) if the site $x_{i-1,j}$
has a different occupation. If $x_{i,j}$ is both a RDB and a LDB, then it is a
single-site domain. If $x_{i,j}$ is neither a RDB or a LDB, then it is an
interior site of a domain.

######  Definition 3.3.

A bit flip on $x_{i,j}$ is an operation which changes the occupation of that
site, either from $1$ to $0$ or from $0$ to $1$, and leaves all other sites in
$\vec{x}_{j}$ unchanged.

######  Definition 3.4.

A maximal state is a state of $L$ with the maximum number of domains.

######  Lemma 3.1.

Every ring of odd circumference $L$ has an even number of domains for all
occupation states.

###### Proof.

Since $L$ is odd, it is not possible to have only single-site domains on a
ring. For all states of odd $L$, at least one domain must contain at least two
sites. For the maximal state, only one domain will contain two sites, all
others will be a single-state domain. Thus with periodic boundary conditions,
when $L$ is odd, the maximum number of domains is $L-1$, which is even.

Any other state can be made from a maximal state by a series of bit flips on
successive sites, up to $L$ sites. On each site there are four possible cases:

Case 1: The site $x_{i,j}$ is a single-site domain. By definition, this means
the occupation states of $x_{i-1,j}$ and $x_{i+1,j}$ are the same. Therefore
when the $x_{i,j}$ site is bit flipped, its single-site domain will merge with
the two domains to the right and to the left, reducing the total number of
domains by two (three original domains become one). In a maximal state, there
are an even number of domains, and bit flipping any of the single site domains
results in a net loss of two domains, so the total number of domains remains
even.

Case 2: The site $x_{i,j}$ is an interior site of a domain. This is
essentially the reverse of Case 1; when this site is bit flipped, it will then
become a single site domain, with two new domains on either side. The original
single domain becomes three, for a net gain of two domains. Since the number
of domains in a maximal state is even, if a site interior to a domain is bit
flipped, two domains are added and the total remains even.

Case 3: The site $x_{i,j}$ is a LDB, but not a RDB. By definition, since
$x_{i,j}$ is a LDB but not a RDB it cannot be a single-site domain. Bit
flipping this site will therefore shift the domain boundary to the right one
site, merging the $x_{i,j}$ site with the domain that previously had the RDB
at $x_{i-1,j}$ and creating a new LDB at $x_{i+1,j}$. The total number of
domains remains the same. Since the maximal state has an even number of
domains, bit flipping a site which is a LDB but not a single-site domain
results in the same even number of domains.

Case 4: The site $x_{i,j}$ is a RDB, but not a LDB. Again by definition, this
site cannot be a single-site domain. Similarly to Case 3, bit flipping this
site will result in the domain boundary shifting left one site, making
$x_{i,j}$ the RDB previously at $x_{i+1,j}$. The total number of domains again
remains the same. Therefore all possible cases of bit flips from the maximal
state result in a state with an even number of domains. ∎

Remark: The above proof can be extended to even $L$ to formulate an alternate
proof of Lemma 2.3.

######  Definition 3.5.

A predecessor state of $\vec{x}_{j+1}$ is any state $\vec{x}_{j}$ that evolves
to $\vec{x}_{j+1}$ under the time evolution specified in (5)

######  Lemma 3.2.

A state on a ring of circumference $L$ has a predecessor state if, and only
if, it is even.

###### Proof.

We have already shown in Lemma 2.3 that all Rule 60 time evolution results in
even states, so no ${}_{L}(\text{odd})$ state has a predecessor. If a
predecessor does exist, it must be constructible by the following inverse of
Rule 60. Here we construct an $\vec{x}_{j}$ configuration from the knowledge
of the $\vec{x}_{j+1}$ configuration. This inverse rule is

$\displaystyle\begin{array}[]{lcc}x_{0,j}=1&\\\
x_{i,j}=\text{Xor}(x_{i-1,j},\,x_{i-1,j+1})&\text{ for
}&i=1,\,2,\,\ldots,\,L-1.\end{array}$ (9)

Because of the periodic boundary condition, there is also the self-consistency
condition

$x_{0,j}=\text{Xor}(x_{L-1,j},\,x_{L-1,j+1}).$ (10)

Note that there is a freedom in the initial choice of $x_{0,j}$, we have
chosen $1$ for this proof, but one can equally choose $0$ and the rest of the
construction goes through as above.

We will now prove that all even states have a predecessor by induction on $L$.

For the base case, we take $L=2$ at some time $j+1$. There are four possible
configurations:

$\displaystyle[00]$ $\displaystyle[01]$ $\displaystyle[10]$
$\displaystyle[11]$

Applying the inverse rule (9), we find:

$[00]$ has predecessor $[11]$

$[11]$ has predecessor $[10]$

The configurations $[01]$ and $[10]$ do not have predecessors, as application
of the inverse rule (9) can not satisfy the self-consistency condition in both
cases. Thus, we have proven by demonstration that all even states for $L=2$
have a predecessor state, and no odd states have a predecessor.

The induction hypothesis is that this holds for some $L=N$, so that given an
$\vec{x}_{j+1}$ of length $N>2$

$\displaystyle x_{i,j}$ $\displaystyle=\text{Xor}(x_{i-1,j},\,x_{i-1,j+1})$
$\displaystyle\text{ for }i=1,\,2,\,\ldots,\,N-1$
$\displaystyle\text{Xor}(x_{N-1,j},\,x_{N-1,j+1})=1$
$\displaystyle\text{for}_{N}(\text{even})$
$\displaystyle\text{Xor}(x_{N-1,j},\,x_{N-1,j+1})=0$
$\displaystyle\text{for}_{N}(\text{odd})$

We now prove the property holds for $N+1$. There are two cases:
Case 1: The ring of circumference $N$ had an even occupation, so by the
induction hypothesis

$\text{Xor}(x_{N-1,j},\,x_{N-1,j+1})=1.$

Extending the ring to $N+1$ means adding another site, $x_{N,j+1}$, which can
be either occupied or unoccupied. If the site is occupied then the total
occupation of $N+1$ is odd. By (9)

$x_{N,j}=\text{Xor}(x_{N-1,j},\,x_{N-1,j+1})=1$

but by the self-consistency condition (10)

$x_{0,j}=\text{Xor}(x_{N,j},\,x_{N,j+1})=0$

which is a contradiction, and so there can be no predecessor state.

If $x_{N,j+1}$ is unoccupied, then again

$x_{N,j}=\text{Xor}(x_{N-1,j},\,x_{N-1,j+1})=1$

but now

$x_{0,j}=\text{Xor}(x_{N,j},\,x_{N,j+1})=\text{Xor}(1,0)=1$

and the predecessor state has been successfully constructed using rule (9).

Case 2: The ring of circumference $N$ had an odd occupation, so by the
induction hypothesis

$\text{Xor}(x_{N-1,j},\,x_{N-1,j+1})=0.$

If the additional site is occupied, so $x_{N,j+1}=1$ then the occupation of
$N+1$ is now even, and

$x_{0,j}=\text{Xor}(x_{N,j},\,x_{N,j+1})=\text{Xor}(1,0)=1$

and the predecessor state has been successfully constructed.

If $x_{N,j+1}=0$, then the $N+1$ state remains odd and

$x_{0,j}=\text{Xor}(x_{N,j},\,x_{N,j+1})=\text{Xor}(0,0)=0$

in contradiction to (10) and a predecessor state can not be constructed.
Therefore, the ring of circumference $N+1$ has a predecessor state only if the
occupation is even. ∎

######  Lemma 3.3.

Every ${}_{L}(\text{even})$ state on a ring has two possible predecessor
states, which are related by bit flipping every site.

###### Proof.

In Lemma 3.2 we demonstrated that every ${}_{L}(\text{even})$ state has a
predecessor, but noted that there is a freedom in choosing the occupation of
the initial site of the predecessor state. Since the Xor operation only
considers relationships between pairs of sites, $x_{i,j}$ and $x_{i-1,j}$, bit
flipping every site preserves the outcome of the Xor operation and generates a
second predecessor. Since this is the only freedom in constructing the
predecessor states, these two states are the only predecessor states. ∎

######  Lemma 3.4.

On a ring of odd circumference $L$, every ${}_{L}(\text{even})$ state has one
${}_{L}(\text{odd})$ precursor.

###### Proof.

Construct an $\vec{x}_{j}$ state from a given $\vec{x}_{j+1}$ by applying (9).
If the occupation is odd, then the proof is complete. If the occupation is
even, then by Lemma 3.3 the second state can be found by bit flipping every
site of the even state. Since $L$ is odd, this flipped state will be odd. ∎

######  Theorem 3.5.

For a ring whose circumference is a Mersenne number, $L=2^{n}-1$
($n=2,\,3,\,\ldots$), every ${}_{(L=2^{n}-1)}(\text{even})$-type initial state
($2^{L-1}$ of them) initiates an $L$-periodic trajectory. These periodic
trajectories exhaust the set of periodic trajectories.

Remark: The proof that follows is not new. The most elegant one can be found
in [7] Lemma 1, Corollary 3. Our proof differs in several details from the one
in [7]. In particular, it does not rely on the properties of the propagators,
and as such is technically close to the one given in [5] for $L=2^{n}$.

###### Proof.

Regard the state of the system at an instance of time $j$ as a vector
$\vec{x}_{j}$. Observe that

$\displaystyle\vec{x}_{j}\stackrel{{\scriptstyle\mod
2}}{{=}}(\hat{I}+\hat{S}_{L})^{j}\vec{x}_{0}\,\,,$

where $\hat{I}$ is a unit matrix, $\hat{S}_{L}$ is a left shift operator:

$\displaystyle\begin{array}[]{lcc}(\hat{S}_{L}\vec{x})_{i}=(\vec{x})_{i-1}&\text{
for }&i=1,\,2,\,\ldots,\,L-1\\\
(\hat{S}_{L}\vec{x})_{0}=(\vec{x})_{L-1}&&\end{array}\,\,.$

Below we will use a particular obvious property of the shift operator:

$\displaystyle\left(\hat{S}_{L}\right)^{L}=\hat{I}\,\,.$ (11)

Now, let us consider $L=2^{n}-1$ and apply the Rule 60 $L+1=2^{n}$ times. We
get

$\displaystyle(\hat{I}+\hat{S}_{L})^{L+1}=\quad\,$
$\displaystyle\hat{I}+\sum_{r=1}^{L}\left(\begin{array}[]{c}L+1\\\
r\end{array}\right)(\hat{S}_{L})^{r}+(\hat{S}_{L})^{L+1}$
$\displaystyle=\quad\,$
$\displaystyle\hat{I}+\sum_{r=1}^{2^{n}-1}\left(\begin{array}[]{c}2^{n}\\\
r\end{array}\right)(\hat{S}_{L})^{r}+\hat{S}_{L}$
$\displaystyle\stackrel{{\scriptstyle\mod 2}}{{=}}$
$\displaystyle\hat{I}+\hat{S}_{L}\,\,.$

Above, we used the properties the property (11) and Lemma 2.5. In the other
words, for any initial state $\vec{x}_{0}$,

$\vec{x}_{L+1}=\vec{x}_{1}\,\,.$

Now, thanks to the Lemma 3.4, $\vec{x}_{1}$ can be represented by _any_
${}_{L}(\text{even})$ state. Shifting the time by one step back, we get that
for _any_ ${}_{(L=2^{n}-1)}(\text{even})$ initial state $\vec{x}_{0}$,

$\vec{x}_{L}=\vec{x}_{0}\,\,.$

This is exactly the property we need.

Finally, according to the Corollary 2.3.1, no other periodic motions are
possible. ∎

######  Corollary 3.5.1.

For a ring whose circumference $L=2^{n}-1=\text{\rm prime}$ is a Mersenne
prime, every ${}_{L}(\text{even})$-type initial state diferent from all zeros
($2^{L-1}-1$ of them) initiates a $L$-fundamental-periodic trajectory.

###### Proof.

According to the Theorem 3.5, $L$ is always one of the periods of _every_
periodic trajectory. Thus, it is the maximal fundamental period Lemma 2.4 is
devoted to:

$M_{\text{\rm max}}=L\,\,.$

The second part of the Lemma 2.4 says that all the smaller periods divide
$M_{\text{\rm max}}$. But in the case of $L$ being a Mersenne prime, $L$ is
prime and as such has no factors besides unity (and the “all zeros” initial
state is the only one that is capable of producing a unit period) and itself.
Hence $L$ is indeed the only available fundamental period.

Finally, while _any_ ${}_{(L=2^{n}-1)}(\text{even})$-type state initiates a
periodic trajectory, these exhaust the list. Indeed, according to the
Corollary 2.3.1, the ${}_{L}(\text{odd})$-type states cannot initiate periodic
trajectories. ∎

###  3.2 Propagator

To aid in the aim of explicitly constructing all ground states of the Newman-
Moore model, let us first introduce a propagator, i.e. a solution of (5) that
has only one site occupied at $j=0$. Namely let us introduce

$\displaystyle\mathcal{P}_{(i,j)\leftarrow(i_{0},0)}$ (12)

that obeys

$\displaystyle\begin{array}[]{lcc}\mathcal{P}_{(i,0)\leftarrow(i_{0},0)}=\delta_{i,i_{0}}\\\
\mathcal{P}_{(i,j+1)\leftarrow(i_{0},0)}=\text{Xor}(\mathcal{P}_{(i-1,j)\leftarrow(i_{0},0)},\,\mathcal{P}_{(i,j)\leftarrow(i_{0},0)})&\text{
for }&i=1,\,2,\,\ldots,\,L-1\\\
\mathcal{P}_{(0,j+1)\leftarrow(i_{0},0)}=\text{Xor}(\mathcal{P}_{(L-1,j)\leftarrow(i_{0},0)},\,\mathcal{P}_{(0,j)\leftarrow(i_{0},0)})\\\
\quad j=0,\,1,\,\ldots\\\ \quad i_{0}=0,\,1,\,L-1\end{array}$ (18)

######  Theorem 3.6.

For $L=2^{n}-1$ ($n=2,\,3,\,\ldots$),

$\displaystyle\mathcal{P}_{(i,j)\leftarrow(0,0)}=\begin{cases}\binom{j}{i}_{L}\\!\\!\\!\\!\mod{2},\quad
0\leq j\leq L-1\\\ (1-\delta_{i,0}),\quad j=L\\\
\mathcal{P}_{(i,(j-1\\!\\!\\!\\!\mod L)+1)\leftarrow(0,0)},\quad
j>L\end{cases}$ (19)

where $\binom{j}{i}_{L}$ is the q-binomial coefficient defined as

$\binom{n}{m}_{q}=\prod_{k=0}^{m-1}\frac{1-q^{n-k}}{1-q^{k+1}}$ (20)

###### Proof.

For the initial state, note that

$\displaystyle\binom{0}{i}_{L}$ $\displaystyle=\begin{cases}1\quad i=0\\\
0\quad\text{otherwise}\end{cases}$ (21)

so the $j=0$ state has a single occupied site at $i=0$ as desired.

As stated in (5), , the time evolution is given by

$\displaystyle x_{i,j+1}=\text{Xor}(x_{i-1,j},\,x_{i,j}),$ $\displaystyle\quad
0<j\leq L-1.$ (22)

This is equivalent to addition mod 2, so

$x_{i,j+1}=\left[x_{i-1j}+x_{ij}\right]\mod{2}.$ (23)

We make use of the known recursion relation for the q-binomial coefficients

$\binom{j+1}{i}_{L}=\binom{j}{i-1}_{L}+L^{i}\binom{j}{i}_{L}.$ (24)

Note that for $L=2^{n}-1$

$L^{i}\mod{2}=1\quad\forall i$ (25)

since all terms in the expansion of $(2^{n}-1)^{i}$ will be powers of two and
their multiples, except the $\binom{i}{i}\left(2^{n}\right)^{0}(-1)^{i}$ term.
Thus

$\binom{j+1}{i}_{L}\mod{2}=\left[\binom{j}{i-1}_{L}+\binom{j}{i}_{L}\right]\mod{2},$
(26)

which is equivalent to (23). ∎

For other propagators, the translational invariance dictates

######  Theorem 3.7.

$\displaystyle\begin{split}&\mathcal{P}_{(i,j)\leftarrow(i_{0},0)}=\left\\{\begin{array}[]{lcc}\mathcal{P}_{(i+i_{0},j)\leftarrow(i_{0},0)}&\text{\rm
for }&i\leq 2(L-1)-i_{0}\\\
\mathcal{P}_{(i+i_{0}-(L-1),j)\leftarrow(i_{0},0)}&\text{\rm for
}&i>2(L-1)-i_{0}\end{array}\right.\\\ &i_{0}=0,\,1\,\ldots,\,L-1\\\
&i=0,\,1\,\ldots,\,L-1\\\ &j=0,\,1\,\ldots,\,L-1\end{split}$ (27)

###### Proof.

This follows directly from the definition of the propagator. ∎

The propagator (27) can be used to propagate the Rule 60 automaton with
$L=2^{n}-1$:

$\displaystyle
x_{i,j}=\sum_{i_{0}=0}^{L-1}\mathcal{P}_{(i,j)\leftarrow(i_{0},0)}x_{i_{0},j}$
$\displaystyle i=0,\,1\,\ldots,\,L-1$ $\displaystyle j=0,\,1\,\ldots,\,L-1$

Note that the formulas ((19), (27)) are not directly applicable for $L\neq
2^{n}-1$.

##  4 Application to the Newman-Moore model

###  4.1 Connection between the Rule 60 automaton and the Newman-Moore model

To connect the Rule 60 automaton and the two-dimensional spin lattices, one
considers the so-called Triangular Plaquette Model or Newman-Moore model [3,
1, 10]:

$\displaystyle
H=-J\sum_{i=0}^{L-1}\sum_{j=0}^{M-1}\sigma_{i,j}\sigma_{i+1,j}\sigma_{i+1,j+1}$
(28)

with periodic boundary conditions. Notice that according to the rule 60, it is
always the case that

$\displaystyle(1-2x_{i,j})(1-2x_{i+1,j})(1-2x_{i+1,j+1})=1\,\,.$

Hence, every $M$-periodic trajectory of (5) will produce a ground state of the
Hamiltonian (28), whose energy is

$\displaystyle E_{\text{ground}}=E_{\uparrow\ldots\uparrow}=-JLM\,\,,$

if one uses the following association:

$\sigma_{i,j}=(1-2x_{i,j})\,\,.$

Notice that “all spins up” is one of the ground states.

###  4.2 Newman-Moore results inferred from the properties of the Rule 60
automaton

We are now ready to address the principal application of our cellular automata
result: the Newman-Moore model. In what follows, we will be able to explicitly
list all the ground state configurations of the Newman-Moore model (28) with

$\displaystyle M=L=2^{n}-1$ $\displaystyle n=2,\,3,\,\ldots\,\,,$

i.e. for a square Newman-Moore lattice with the sides equal to one of the
Mersenne _numbers_.

We assert the following.

######  Theorem 4.1.

The ground state of the Newman-Moore model (28) with $M=L=2^{n}-1$ is
$2^{L-1}$ degenerate.

###### Proof.

The proof follows directly from the Theorem 3.5. ∎

We will now provide explicit expressions for the ground state configurations
of the Newman-Moore model (28) on an $L\times L$ square lattice with $L$ being
equal to one of the Mersenne numbers $2^{n}-1$.

The ground state configurations will be labeled by the state of the $0$-th
row:

$[\sigma_{0,0},\,\sigma_{1,0},\,\ldots,\,\sigma_{L-1,0}]\,\,,$

where number of spin-down sites, $\sigma=-1$ is even: there are $2^{L-1}$ such
configurations, consistent with the $2^{L-1}$-fold degeneracy.

######  Theorem 4.2.

The ground states of the Newman-Moore model (28) on an $L\times L$ square
lattice with $L=(2^{n}-1)$ ($n=2,\,3,\,\ldots$) read

$\displaystyle\begin{split}&\sigma_{i,j}=\sum_{i_{0}=0}^{L-1}\left(1-2\mathcal{P}_{(i,j)\leftarrow(i_{0},0)}\frac{1-\sigma_{i_{0},j}}{2}\right)\\\
&i=0,\,1\,\ldots,\,2^{n}-2\\\ &j=0,\,1\,\ldots,\,2^{n}-2\end{split}$ (29)

where the propagator $\mathcal{P}_{(i,j)\leftarrow(i_{0},0)}$ is given by (27)
and $[\sigma_{0,0},\,\sigma_{1,0},\,\ldots,\,\sigma_{2^{n}-2,0}]$ is any of
the $2^{L-1}$ spin configurations characterised by the number of the
$\sigma=-1$ sites being even.

##  5 Summary of results

The Theorem 3.5 plays the central role in our paper. It allows to identify all
the initial conditions that lead to periodic trajectories in a Rule 60
cellular automaton on a ring whose circumference $L$ is given by one of the
Mersenne numbers, $L=2^{n}-1$ . Furthermore, all periodic trajectories are
shown to share the period $L$. The Theorem 3.6 constructs the propagators for
this model.

The above results allow us to list all the ground states of the Newman-Moore
model explicitly. These results are summarised in Theorems 4.1 and 4.2

## Acknowledgements

The authors wish to thank Konstantinos Sfairopoulos for bringing this problem
to their attention, for providing invaluable guidance on references, and for
many useful discussions.

A significant portion of this work was produced during the thematic trimester
on “Quantum Many-Body Systems Out-of-Equilibrium”, at the Institut Henri
Poincaré (Paris): MO is immeasurably grateful to the organizers of the
trimester, Rosario Fazio, Thierry Giamarchi, Anna Minguzzi, and Patrizia
Vignolo, for an opportunity to be a part of it. MO and JC-P wish to express
their gratitude to the International High IQ Society for the networking
opportunities that it offers.

#### Funding information

MO was supported by the NSF Grant No. PHY-1912542. MO would like to thank the
Institut Henri Poincaré (UAR 839 CNRS-Sorbonne Université) and the LabEx
CARMIN (ANR-10-LABX-59-01) for their support.

## References

  * [1] M. E. J. Newman and C. Moore,  _Glassy dynamics and aging in an exactly solvable spin model_ ,  Phys. Rev. E 60, 5068 (1999),  10.1103/PhysRevE.60.5068.
  * [2] B. Yoshida and A. Kubica,  _Quantum criticality from ising model on fractal lattices_ (2014), 1404.6311.
  * [3] R. L. Jack, L. Berthier and J. P. Garrahan,  _Static and dynamic length scales in a simple glassy plaquette model_ ,  Phys. Rev. E 72, 016103 (2005),  10.1103/PhysRevE.72.016103.
  * [4] A. Lipowski, D. Johnston and D. Espriu,  _Slow dynamics of ising models with energy barriers_ ,  Phys. Rev. E 62, 3404 (2000),  10.1103/PhysRevE.62.3404.
  * [5] K. Sfairopoulos, L. Causer, J. F. Mair and J. P. Garrahan,  _Boundary conditions dependence of the phase transition in the quantum Newman-Moore model_ ,  Phys. Rev. B 108, 174107 (2023),  10.1103/PhysRevB.108.174107.
  * [6] O. Martin, A. M. Odlyzko and S. Wolfram,  _Algebraic properties of cellular automata_ ,  Commumun. Math. Phys. 93, 219 (1984).
  * [7] A. Ehrlich,  _Periods in Ducci’s $n$-number game of differences_,  The Fibonacci Quarterly 28, 302 (1990).
  * [8] F. Breuer, E. Lötter and B. van der Merwe,  _Ducci-sequences and cyclotomic polynomials_ ,  Finite Fields and Their Applications 13(2), 293 (2007),  https://doi.org/10.1016/j.ffa.2005.11.003.
  * [9] D. Schepler,  _How to prove binomial coefficient ${2^{n}\choose k}$ is even number?_,  Mathematics Stack Exchange,  URL:https://math.stackexchange.com/q/3087725 (version: 2019-01-25),  https://math.stackexchange.com/q/3087725.
  * [10] _Newman-Moore code_ ,  In V. V. Albert and P. Faist, eds., _The Error Correction Zoo_ (2023).