#  Rheomergy: Collective behavior mediated by active flow-based recruitment

S Ganga Prasath  School of Engineering and Applied Sciences, Harvard
University, Cambridge MA 02138.  L Mahadevan<EMAIL_ADDRESS>School of
Engineering and Applied Sciences, Harvard University, Cambridge MA 02138.
Department of Physics, Harvard University, Cambridge MA 02138.  Department of
Organismic and Evolutionary Biology, Harvard University, Cambridge 02138

###### Abstract

The physics of signal propagation in a collection of organisms that
communicate with each other both enables and limits how active excitations at
the individual level reach, recruit and lead to collective patterning.
Inspired by the patterns in a planar swarm of bees that release pheromones,
and use fanning flows to recruit additional bees, we develop a theoretical
framework for patterning via active flow-based recruitment. Our model
generalizes the well-known Patlak-Keller-Segel model of diffusion-dominated
aggregation and leads to more complex phase space of patterns spanned by two
dimensionless parameters that measure the scaled stimulus/activity and the
scaled chemotactic response. Together these determine the efficacy of signal
communication that leads to a variety of migration and aggregation patterns
consistent with observations.

In groups of individual agents in both natural and artificial settings,
communication is central to achieving a coordinated collective response.
Indeed, the survival of the colony pivots around robust communication in the
context of foraging Couvillon (2012); Jackson and Ratnieks (2006); Hölldobler
(1978), protection against invasion by an intruder Nouvian  _et al._ (2016);
Trhlin  _et al._ (2011), brood care in eusocial clusters Gordon (1996);
Hölldobler  _et al._ (1990) and other tasks inside the nest Greene and Gordon
(2003); Free (1987). The efficacy of communication is limited by the dynamics
of signal propagation among agents through the environment Wingreen and Levin
(2006); Seeley (2009); Gernat  _et al._ (2018); Prasath  _et al._ (2021). Over
the last century, stigmergy has emerged as one of the most well studied
mechanisms to achieve robust communication Theraulaz and Bonabeau (1999);
Perna and Theraulaz (2017); Camazine  _et al._ (2003): here animals modify the
environment by using localized signals that serve as recruitment cues in the
context of functional behavior. However, stigmergy is just one way of
environmentally modulated communication. Organisms use slow diffusive signals
Patlak (1953); Keller and Segel (1970); Hillen and Painter (2009) that
underlie the principle of quorum sensing Waters and Bassler (2005), passive
and active fluid flows Ocko and Mahadevan (2015); Jost  _et al._ (2007);
Mathijssen  _et al._ (2019), and even elastic deformations Peleg  _et al._
(2018) to transmit information via physical communication channels. Indeed,
social insects such as honey bees (e.g. the species apis mellifera) use all
these different modes of communication and actively modify their local
mechanical, thermal and hydrodynamic environment to transmit information and
coordinate their activities and functions.

Figure 1: $(a)$ Sequence of images showing recruitment of worker bees by
collective fanning to the nest where the queen resides (Image: Peleg Lab at UC
Boulder). We see a zoomed in view close to the queen’s cage (see also Nguyen
_et al._ (2021) and Ranch (2013)). $(b)$ Interaction rules between the
variables in our model in Eqs. 1-5: agent density $\varrho(\mathbf{x}_{||},t)$
which is confined to a 2D surface $\mathbf{x}_{||}=(x,y)$ generates a
pheromone field of concentration $c(\mathbf{x},t)$ that gets fanned through
the environment by the actively generated bee velocity field
$\mathbf{u}_{f}(\mathbf{x},t)$. The advected pheromone causes either
additional agents (bees) to be recruitment or causes them to migrate. $(c)$
Phase diagram of rheomergy represented using two non-dimensional ratios:
Active Peclet number which is the ratio of pheromone decay rate to active
fluid shear rate, $\mathsf{R}_{c}/\mathsf{V}_{c}$ and recruitment strength to
chemotaxis strength, $\mathsf{R}_{\varrho}/\mathsf{V}_{\varrho}$ (see text for
details). Four phases characterizing the collective behavior of agents exist
depending on the region in phase-space: polarized recruitment, migration,
local aggregation and isotropic recruitment (similar to that seen in the
Patlak-Keller-Segel model).

A particularly striking example of active, long-range and fast communication
is seen in bees that settle on a substrate and flap their wings to generate a
fluid flow parallel to their body axis while releasing pheromones into it
Peters  _et al._ (2019); Nguyen  _et al._ (2021); Reinhard and Srinivasan
(2009); Peters  _et al._ (2017, 2017); Nguyen  _et al._ (2021); Ranch (2013).
In Fig. 1$(a)$ and SI Fig. S1, we show how fluid flow mediated communication
leads to a characteristic triangular swarm shape that migrates towards the
source of pheromones Nguyen  _et al._ (2021). These patterns arise on time
scales of minutes spontaneously from a disordered state and can extend to
include hundreds of bees.

To understand the nature and efficacy of flow-mediated patterning, we need to
couple the three observable spatiotemporal fields that are shown schematically
in Fig. 1$(b)$: the bee density $\varrho(\mathbf{x}_{||},t)$ (and orientation
$\hat{\mathbf{p}}(\mathbf{x}_{||},t)$ denotes head forward) as a function of
location on the plane $z=0$ denoted by $\mathbf{x}_{||}=(x,y)$ and time, the
scalar pheromone concentration field $c(\mathbf{x},t)$ and the active fluid
flow generated by the bees $\mathbf{u}_{f}(\mathbf{x},t)$. Bees disperse
pheromones from their Nasonov glands into the active flow that they generate
by fanning. This leads to the active recruitment of bees from the environment
and their migration which then leads to an enhanced rate of pheromone
dispersal. The positive feedback leads to spatio-temporal patterns of growing,
self-focusing and migrating swarms. While there are many qualitative
observations of the process Reinhard and Srinivasan (2009); Winston (1991),
the quantitative collective dynamics of recruitment, agent-environment
feedback, and inter-agent interaction mediated by volatile pheromones is
largely unexplored. Here, we develop a minimal model for fluid flow mediated
communication and recruitment which we dub rheomergy. As we will see, our
model broadens the scope of the classical Patlak-Keller-Segel (PKS) model
Patlak (1953); Keller and Segel (1970); Hillen and Painter (2009) for
chemotaxis by accounting for dynamics of active fluid flow and the behavioral
response of the agents to flow-mediated communication leading to a more
complex phase space of patterns.

Continuum model of rheomergy: The evolution equations for the agent density
$\varrho(\mathbf{x}_{||},t)$, the pheromone field $c(\mathbf{x},t)$ and the
fluid flow generated by the agents, $\mathbf{u}_{f}(\mathbf{x},t)$ can be
written as,

$\displaystyle\partial_{t}\varrho+\nabla_{||}\cdot(\mathbf{u}_{a}\varrho)=$
$\displaystyle\
D_{a}\nabla_{||}^{2}\varrho+(1-\varrho/\varrho_{\text{max}})[k_{b}\delta\varrho$
$\displaystyle\ +\alpha_{p}\Theta(\delta c)],$ (1)
$\displaystyle\mathbf{u}_{a}(\mathbf{x}_{||},t)=$ $\displaystyle\
\chi\nabla_{||}c,$ (2)
$\displaystyle\partial_{t}c+\nabla\cdot(\mathbf{u}_{f}c)=$ $\displaystyle\
D_{p}\nabla^{2}c+k_{+}\varrho-k_{-}c,$ (3)
$\displaystyle\partial_{t}\mathbf{u}_{f}=$ $\displaystyle\ -\frac{\nabla
p}{\rho_{f}}+\nu\nabla^{2}\mathbf{u}_{f}-\nu_{s}\mathbf{u}_{f},$ (4)
$\displaystyle\nabla\cdot\mathbf{u}_{f}=$ $\displaystyle\ 0.$ (5)

Equation 1 describes the migration, relaxation and recruitment of bees as a
consequence of the agents relaxing to a nominal ambient density; here the
chemotactic velocity defined by Eq. 2 is $\mathbf{u}_{a}(\mathbf{x}_{||},t)$
with a chemotactic gain $\chi$ and $\nabla_{||}$ denotes in-plane gradients,
defined along the plane $z=0$, the first term on the right side captures
diffusion of agents with a diffusivity $D_{a}$, the second term denotes
relaxation to a nominal density $\varrho_{o}$ with
$\delta\varrho=(\varrho_{o}-\varrho)$ and $k_{b}$ denotes the relaxation rate,
and the last term characterizes recruitment with $\alpha_{p}$ being the
sensitivity to pheromone above a threshold $c^{*}$ where $\delta c=(c-c^{*})$,
and $\Theta(\bullet)$ is the Heaviside function. The agent density is limited
to a maximum $\varrho_{\text{max}}$ through the term
$(1-\varrho/\varrho_{\text{max}})$. Equation 3 characterizes the dynamics of
the pheromone that is advected by the active flow velocity $\mathbf{u}_{f}$
generated by the swarm, diffuses into the environment with a diffusivity
$D_{p}$, is produced by bees at a rate $k_{+}$ and decays at a rate $k_{-}$.
Equations 4-5 follow from the dynamics of (incompressible) flow of fluid with
viscosity $\nu$ and density $\rho_{f}$, and an effective friction factor
$\nu_{s}$ associated with fluid motion relative to the substrate. Since this
flow is actively generated by bees, we model this using a shear rate at the
boundary proportional to their local density:
$\partial_{z}\mathbf{u}_{f}(z=0)=\beta\varrho\hat{\mathbf{p}}$, where
$\hat{\mathbf{p}}$ is the local polarity of the bees which generates a flow
velocity opposite to their orientation and parallel to the surface (see Fig.
1$(a)$ and Peters  _et al._ (2017)), with $\beta>0$. Equations 1-5 form a
closed system for the variables $\\{\varrho,c,\mathbf{u}_{f},p\\}$ which are
complete when we specify initial and boundary conditions for all the fields
(see SI sec. III for details).

We make Eqs. 1-5 dimensionless using the natural length-scale
$\zeta^{-1}=(\nu/\nu_{s})^{1/2}$ and time-scale $k_{b}^{-1}$ for Eq. 1, and
$k_{-}^{-1}$ for Eq. 3 (see SI tab. I for other length and time-scales in the
system), the homogeneous relaxation density $\varrho_{o}$ as the density-
scale, the maximum pheromone concentration $c_{o}$ as the concentration-scale
(see SI sec. I for details), resulting in 7 non-dimensional numbers associated
with scaled magnitudes of chemotactic strength, $\mathsf{V}_{\varrho}=(\chi
c_{o}\zeta^{2})/k_{b}$; agent diffusivity,
$\mathsf{D}_{\varrho}=(D_{a}\zeta^{2}/k_{b})$; recruitment strength,
$\mathsf{R}_{\varrho}=\alpha_{p}c_{o}/(k_{b}\varrho_{o})$; shear,
$\mathsf{V}_{c}=(\beta\varrho_{o}/k_{-})$; pheromone decay,
$\mathsf{R}_{c}=k_{+}\varrho_{o}/(k_{-}c_{o})$; pheromone diffusivity,
$\mathsf{D}_{c}=(D_{p}\zeta^{2}/k_{-})$; friction factor,
$\mathsf{D}_{u}=\nu_{s}/(\beta\varrho_{o})$. This allows us to rewrite the
variables after scaling as:
$\tilde{\varrho}\rightarrow\varrho/\varrho_{o},\tilde{c}\rightarrow
c/c_{o},\tilde{\mathbf{u}}_{f}\rightarrow\mathbf{u}_{f}(\zeta/\beta\varrho_{o}),\tilde{p}\rightarrow
p\ \zeta^{2}/(\beta\varrho_{o})^{2}\rho_{f}$.

In the case of active flow-driven communication, relevant for social insects,
we assume that the agents operate in the interesting neighborhood of the
threshold pheromone concentration $c\sim c^{*}$, we therefore replace
$\Theta(\delta c)\rightarrow\mathsf{R}_{\varrho}\delta c$. In the limit of
rapid equilibration of the pheromone concentration and a rapid relaxation to
steady velocity relative to agent recruitment rates and pheromone decay rates
relevant for social insects, we assume that $\partial_{t}c\approx 0$ and
$\partial_{t}\mathbf{u}_{f}\approx 0$, so that Eqs. 1-5 in scaled units reduce
to

$\displaystyle\partial_{t}\varrho+\mathsf{V}_{\varrho}\nabla_{||}\cdot(\varrho\nabla_{||}c)=$
$\displaystyle\
\mathsf{D}_{\varrho}\nabla_{||}^{2}\varrho+(1-{\varrho}/{\varrho_{\text{max}}})[(1-\varrho)$
$\displaystyle\ +\mathsf{R}_{\varrho}\delta c],$ (6)
$\displaystyle\mathsf{D}_{c}\nabla^{2}c+\mathsf{V}_{c}\nabla\cdot(\mathbf{u}_{f}c)=$
$\displaystyle\ \mathsf{R}_{c}\varrho-c,$ (7)
$\displaystyle\mathsf{D}_{u}(\nabla^{2}\mathbf{u}_{f}-\mathbf{u}_{f})=$
$\displaystyle\ \nabla p$ (8) $\displaystyle\nabla\cdot\mathbf{u}_{f}=$
$\displaystyle\ 0.$ (9)

When the diffusion of pheromone dominates over advection i.e.
$\mathsf{D}_{c}\sim\mathsf{R}_{c}\gg 1,\mathsf{V}_{c}$ and in the absence of
recruitment and relaxation ($\alpha_{p},k_{b}=0$) or when the agent density
has reached the maximum value ($\varrho=\varrho_{\text{max}}$), Eqs. 6-7
reduce to the PKS equations (see SI sec. IV B for details). In this article,
however, we are interested in the limit relevant for social insects where the
dynamics of the agents is governed primarily by chemotaxis and recruitment,
i.e. $\mathsf{D}_{\varrho}\ll\mathsf{V}_{\varrho}\sim\mathsf{R}_{\varrho}$ and
that the pheromone rapidly equilibrates due to advection, i.e.
$\mathsf{D}_{c}\ll\mathsf{V}_{c}\sim\mathsf{R}_{c}$. From now on, we will
therefore ignore the effects of agent and pheromone diffusivity.

We explore the different limits of the phase-space defined by the pheromone
dynamics (governed by $\mathsf{V}_{c},\mathsf{R}_{c}$) and by agent dynamics
on the other (governed by $\mathsf{V}_{\varrho},\mathsf{R}_{\varrho}$) to
understand the different modes of collective patterning determined by our
model. There are four different regimes characterized by relative magnitudes
of $\mathsf{R}_{\varrho},\mathsf{V}_{\varrho}$ and
$\mathsf{R}_{c},\mathsf{V}_{c}$ (see Fig. 1$(c)$) which we address in a
1-dimensional setting before turning briefly to the 2-dimensional case (see SI
sec. IV for details).

Figure 2: Dynamics of agent density, $\varrho(x,t)$ and pheromone
concentration, $c(x,t)$ in the $(a)$ flow (and recruitment) dominated regime
(Ia in the text) where agents polarize in response to autogenic flow and
recruitment, leading to an increase in the agent density (for
$\mathsf{V}_{c}=0.5,\mathsf{R}_{c}=0.5,\mathsf{V}_{\varrho}=0.0,\mathsf{R}_{\varrho}=300.0$).
The plots on the right show the self-similar evolution of agent density and
the pheromone concentration, consistent with a the collapse of the different
curves when scaled by $xt,x$ before saturation (see SI sec. IV A for details).
$(b)$ Flow (and chemotaxis) dominated regime phase (Ib in the text) where the
agents migrate upstream in the direction of flow through (for
$\mathsf{V}_{c}=5,\mathsf{R}_{c}=0.9,\mathsf{V}_{\varrho}=280.0,\mathsf{R}_{\varrho}=0.0$).
The results were obtained by numerically solving Eqs. 10, 11 (see SI sec. V
for details).

(I) Flow dominated regimes ($\mathsf{V}_{c},\mathsf{R}_{c}\gg 1$) : In the
flow-dominated regime, we explore the collective behavior of the agents as we
vary $\mathsf{R}_{\varrho}/\mathsf{V}_{\varrho}$ \- see Fig. 1$(c)$. Assuming
that the colony is uniformly active and generating flow everywhere, i.e.
$\varrho=1$, the velocity field is given by the solution of Eqs. 8, 9 subject
to the shear boundary condition at the wall:
$\partial_{z}(\mathbf{u}_{f}\cdot\hat{\mathbf{x}}_{||})(z=0)=1$, which yields
a polarized flow field in the in-plane direction
$u_{||}(z)\equiv\mathbf{u}_{f}\cdot\hat{\mathbf{x}}_{||}=-\exp(-z)$. Then the
pheromone concentration (obtained by substituting the above result in Eq. 7)
yields $c(x)~{}=~{}\mathsf{R}_{c}[1-\exp{(x/\mathsf{V}_{c})}]$. This
inhomogeneous pheromone concentration profile driven by a homogeneous agent
density $\varrho=1$ causes the agents to migrate when the chemotactic response
$\mathsf{V}_{\varrho}>0$ (see from Eq. 6), indicating the unstable nature of
this density profile.

Substituting the exponential flow profile into Eq. 7 for the pheromone
dynamics (see SI sec. III for details) in the limit of advection-dominated
pheromone transport leads to the shear becoming proportional to agent density
so that $\partial_{z}u_{||}\sim 1$, and Eqs. 6-8 simplify to the coupled
system for the agent density and the pheromone concentration,

$\displaystyle\partial_{t}\varrho+\mathsf{V}_{\varrho}\partial_{x}(\varrho\partial_{x}c)=$
$\displaystyle\
(1-\varrho/\varrho_{\text{max}})[(1-\varrho)+\mathsf{R}_{\varrho}c],$ (10)
$\displaystyle-\mathsf{V}_{c}\partial_{x}(\varrho c)=$ $\displaystyle\
\mathsf{R}_{c}\varrho-c.$ (11)

To explore the behavior of the collective agent dynamics in the different
limits of phase-space we solve Eqs. 10, 11 (see SI sec. V for details).

(Ia) Recruitment dominated regime ($\mathsf{R}_{\varrho}\gg
1,\mathsf{V}_{\varrho}\ll 1$): In Fig. 2$(a)$ we see that an initially
localized gaussian profile of the agent density $\varrho(x,0)$ centered around
$x=0$ grows and spreads along the direction of pheromone advection due to
flow-driven agent recruitment. Eventually, the agent density saturates to
maximal packing density $\varrho=\varrho_{\text{max}}$. Simultaneously, the
pheromone concentration $c(x,t)$ increases with increasing recruitment and is
also advected by flow (see Fig. 2$(a)$). To understand the scaling behavior of
the spatio-temporal profiles of the agent and pheromone density in the growth
phase, we note that when the pheromone is produced by agents and transported
by the active flow with decay being insignificant, Eq. 11 reduces to
$-\mathsf{V}_{c}\partial_{x}(\varrho c)\sim\mathsf{R}_{c}\varrho$.
Furthermore, when the effects of chemotaxis are small (i.e.
$\mathsf{V}_{\varrho}\rightarrow 0$), the agent density scales as $\varrho\sim
xt$ while the concentration of the pheromone scales as $c\sim x$ (see SI sec.
IV A for a detailed analysis). Normalizing $\varrho(x,t)$ and $c(x,t)$ using
this scaling causes the curves to collapse on top of each other (see rightmost
figures in Fig. 2$(a)$).

(Ib) Chemotaxis dominated regime ($\mathsf{V}_{\varrho}\gg
1,\mathsf{R}_{\varrho}\ll 1$): In Fig. 2$(b)$ we see that an initially
localized agent density collectively migrates up the pheromone gradient.
Simultaneously the concentration of pheromone decreases in magnitude and
migrates with the agents. To understand these observations we note that in the
flow-dominated regime, Eq. 11 implies that
$-\mathsf{V}_{c}\partial_{x}(\varrho c)\sim\mathsf{R}_{c}\varrho$, which when
substituted into Eq. 10 yields:
$\partial_{t}\varrho-(\mathsf{V}_{\varrho}\mathsf{R}_{c}/\mathsf{V}_{c})\partial_{x}\varrho=\mathsf{V}_{\varrho}\partial_{x}(c\partial_{x}\varrho)$.
We see that agents migrate with a speed
$(\mathsf{V}_{\varrho}\mathsf{R}_{c}/\mathsf{V}_{c})$ and diffuse
inhomogeneously due to the pheromone concentration that acts as a spatially
varying diffusivity. Here the emergence of collectively polarized migration
even in the absence of recruitment ($\mathsf{V}_{\varrho}\gg
1,\mathsf{R}_{\varrho}\rightarrow 0$) is due to autonomous fluid flow coupled
with pheromone production.

(II) Pheromone decay dominated regimes ($\mathsf{V}_{c}\ll\mathsf{R}_{c}\sim
1$): In the pheromone dominated regime, we explore the collective behavior of
the agents as we vary $\mathsf{R}_{c}/\mathsf{V}_{c}$ \- see Fig. 1$(c)$. In
the limit of a homogeneous agent density $\varrho=1$, when the effects of
decay dominate ($\mathsf{R}_{c}\gtrsim 1,\mathsf{V}_{c}\rightarrow 0$), Eq. 11
yields $c^{\text{ss}}=\mathsf{R}_{c}$. Perturbing about this steady state
($\varrho^{\text{ss}}=1,c^{\text{ss}}=\mathsf{R}_{c}$) using the ansatz
$\\{\varrho(x,t)-\varrho^{\text{ss}},c(x,t)-c^{\text{ss}}\\}=\\{\hat{\varrho}(q),\hat{c}(q)\\}\exp(iqx+\Omega
t)$ where $\Omega(q)$ is the frequency of oscillation and $q$ the scalar
wavenumber, $\hat{\varrho}(q),\hat{c}(q)$ are the complex amplitudes of agent
density and pheromone concentration and substituting into Eqs. 6, 7 yields the
dispersion relation (derived in SI sec. II),

$\displaystyle\Omega(q)=$ $\displaystyle\
\frac{1}{2}\bigg{\\{}-(1+\lambda)\pm[(\lambda-1)^{2}+4\lambda\mathsf{R}_{c}(\mathsf{R}_{\varrho}+4q^{2}\mathsf{V}_{\varrho})]^{1/2}\bigg{\\}},$
(12)

where we have $\lambda=(k_{-}/k_{b})$ and we use $k_{b}^{-1}$ as the time-
scale to non-dimensionalize the frequency. When $\Re(\Omega(q))>0$ the agent
density becomes unstable; for long wavelengths ($q\rightarrow 0$) Eq. 12
becomes
$\Omega(q)\approx[-(1+\lambda)+\sqrt{(\lambda-1)^{2}+4\lambda\mathsf{R}_{c}\mathsf{R}_{\varrho}}]/2$
and we see that the dynamics is governed by competition between recruitment
sensitivity (from $\mathsf{R}_{\varrho}$) and pheromone decay (from
$\mathsf{R}_{c}$). When $\mathsf{R}_{\varrho}\gg 1$ corresponding to large
recruitment sensitivity this triggers long wavelength perhaps suggestive of
such a mechanism at play when the colony in under threat from an intruder. For
smaller wavelengths ($q\gg 1$) the instability is due to competition between
pheromone decay (contributions from $\mathsf{R}_{c}$) and chemotaxis (from
$\mathsf{V}_{\varrho}$), as $\Omega(q)\approx
2\sqrt{\lambda\mathsf{R}_{c}\mathsf{V}_{\varrho}}q$. This linearly growing
dependence of the unstable frequency on wavenumber is saturated by the effects
of diffusion at large wavenumbers as $\Omega(q)\approx-\mathsf{D}_{c}\lambda
q^{2}$.

Figure 3: Polarized recruitment of active agents and corresponding pheromone
concentration in two-dimensions obtained by solving Eqs. 6, 11 and plotted at
$t=1.5$. Agent density along $x=0$ is set to
$\varrho(0,y,t)=\exp(-y^{2}/2w^{2})$ with $w=2$ and these agents generate flow
$u_{||}$ along positive $x$-axis that results in recruitment of other agents.
The triangular shape is a consequence of chemotaxis driven by gradients in
pheromone along $y$-axis. We choose Neumann boundary condition for
$\varrho(x,y,t),c(x,y,t)$ along other boundaries. The non-dimensional numbers
are chosen to be
$\mathsf{V}_{c}=0.5,\mathsf{R}_{c}=0.5,\mathsf{V}_{\varrho}=25.0,\mathsf{R}_{\varrho}=300.0$
(see SI sec. V for further details).

To understand the nonlinear evolution of the agents in the pheromone-decay
dominated regime, the concentration of the pheromone is proportional to the
density everywhere, $c\approx\mathsf{R}_{c}\varrho$. Substituting this into
Eq. 11 we get,

$\displaystyle\partial_{t}\varrho+(\mathsf{V}_{\varrho}\mathsf{R}_{c}/2)\partial^{2}_{y}\varrho^{2}=$
$\displaystyle\ \mathsf{R}_{\varrho}\mathsf{R}_{c}\varrho+(1-\varrho),$ (13)

The agent density can exhibit two distinct behaviors, determined by the
relative magnitude of chemotaxis and recruitment.

(IIa) Chemotaxis dominated regime ($\mathsf{V}_{\varrho}\gg
1,\mathsf{R}_{\varrho}\ll 1$): In this limit, Eq. 13 reduces to a nonlinear
heat equation with a negative diffusivity. As the agents sense the pheromone
generated by their immediate neighbors (as $\mathsf{V}_{\varrho}\gg 1$) they
migrate towards each other which results in a pile-up of the density (see SI
sec. IV A for details of the similarity solution associated with aggregation).
This aggregation behavior is qualitatively similar to what is seen in the PKS
model, however the dynamics follows a different scaling law.

(IIb) Recruitment dominated regime ($\mathsf{R}_{\varrho}\gg
1,\mathsf{V}_{\varrho}\ll 1$): In this limit, Eq. 13 reduces to a linear
homogeneous equation and the agent density grows exponentially as
$\varrho\sim\exp(\mathsf{R}_{\varrho}\mathsf{R}_{c}t)$. This leads to
homogeneous growth (see SI Fig. S2$(d)$) unlike the flow-dominated scenario
where advection of the pheromone leads to polarized recruitment and
anisotropic patterning.

Analysis in 2D: Our 1D setting and analysis translates naturally to 2D
scenarios. When recruitment is driven by fluid flow
($\mathsf{V}_{c}\sim\mathsf{R}_{c}$), we expect accumulation along the
transverse direction ($x$-axis) when $\mathsf{R}_{\varrho}\gg 1$ as shown in
Fig. 3, with a concomitant narrowing of the accumulated agent density along
$y$-axis seen in Fig. 3. This results in the evolution of the density and the
pheromone concentration towards a triangular shape seen in recent experiments
with bees Nguyen  _et al._ (2021) (shown in Fig. 2$(a)$ inset).

Our study of the spatio-temporal organization of a collective of agents that
actively generate and advect pheromones using autonomous flows to recruit
additional agents generalizes the classic Patlak-Keller-Segel diffusive
aggregation framework and shows the emergence of spatio-temporal patterns of
agents that collectively communicate, migrate and aggregate. We delineate this
in a phase-space defined by the ratio of recruitment strength to chemotaxis
i.e. $\mathsf{R}_{\varrho}/\mathsf{V}_{\varrho}$ and the ratio of shear
strength to pheromone decay strength i.e. $\mathsf{V}_{c}/\mathsf{R}_{c}$. For
large shear strength and recruitment strength
($\mathsf{V}_{c}/\mathsf{R}_{c},\mathsf{R}_{\varrho}/\mathsf{V}_{\varrho}\gg
1$) the agent density exhibits successful polarized recruitment driven by the
autogenic fluid flow. On the other hand, when the chemotaxis strength
dominates
($\mathsf{V}_{c}/\mathsf{R}_{c},\mathsf{V}_{\varrho}/\mathsf{R}_{\varrho}\gg
1$) the agents are able to collectively migrate along the direction of the
fluid flow. When the pheromone decay effects dominate
($\mathsf{R}_{c}/\mathsf{V}_{c},\mathsf{V}_{\varrho}/\mathsf{R}_{\varrho}\gg
1$) agents collectively aggregate which eventually results in the formation of
shocks. Lastly, when the chemotactic velocity is weak but the recruitment
strength is large
($\mathsf{R}_{c}/\mathsf{V}_{c},\mathsf{R}_{\varrho}/\mathsf{V}_{\varrho}\gg
1$) the agents can perform isotropic recruitment. Generalizing our framework
to include the orientation dynamics of agents while accounting for the role of
noise in sensing and actuation are natural next steps in characterizing
decision making and collective response in active, embodied agents that use
rheomergy.

###### Acknowledgements.

We thank Orit Peleg for Fig. 1$(a)$, and the NSF Simons Center for
Mathematical and Statistical Analysis of Biology, the Simons Foundation and
the Seydoux Fund for partial support.

## References

  * Couvillon (2012) M. Couvillon, Insectes sociaux 59, 297 (2012).
  * Jackson and Ratnieks (2006) D. E. Jackson and F. L. Ratnieks, Current biology 16, R570 (2006).
  * Hölldobler (1978) B. Hölldobler, in  _Advances in the Study of Behavior_ , Vol. 8 (Elsevier, 1978) pp. 75–115.
  * Nouvian  _et al._ (2016) M. Nouvian, J. Reinhard, and M. Giurfa, Journal of Experimental Biology 219, 3505 (2016).
  * Trhlin  _et al._ (2011) M. Trhlin, J. Rajchard, _et al._ , Vet. Med 56, 265 (2011).
  * Gordon (1996) D. M. Gordon, Nature 380, 121 (1996).
  * Hölldobler  _et al._ (1990) B. Hölldobler, E. O. Wilson, _et al._ ,  _The ants_ (Harvard University Press, 1990).
  * Greene and Gordon (2003) M. J. Greene and D. M. Gordon, Nature 423, 32 (2003).
  * Free (1987) J. B. Free,  _Pheromones of social bees._ (Chapman and Hall, 1987).
  * Wingreen and Levin (2006) N. S. Wingreen and S. A. Levin, PLoS biology 4, e299 (2006).
  * Seeley (2009) T. D. Seeley,  _The wisdom of the hive: the social physiology of honey bee colonies_ (Harvard University Press, 2009).
  * Gernat  _et al._ (2018) T. Gernat, V. D. Rao, M. Middendorf, H. Dankowicz, N. Goldenfeld, and G. E. Robinson, Proceedings of the National Academy of Sciences 115, 1433 (2018).
  * Prasath  _et al._ (2021) S. G. Prasath, S. Mandal, F. Giardina, J. Kennedy, V. N. Murthy, and L. Mahadevan, bioRxiv (2021).
  * Theraulaz and Bonabeau (1999) G. Theraulaz and E. Bonabeau, Artificial life 5, 97 (1999).
  * Perna and Theraulaz (2017) A. Perna and G. Theraulaz, Journal of Experimental Biology 220, 83 (2017).
  * Camazine  _et al._ (2003) S. Camazine, J.-L. Deneubourg, N. R. Franks, J. Sneyd, E. Bonabeau, and G. Theraulaz,  _Self-organization in biological systems_ (Princeton university press, 2003).
  * Patlak (1953) C. S. Patlak, The bulletin of mathematical biophysics 15, 311 (1953).
  * Keller and Segel (1970) E. F. Keller and L. A. Segel, Journal of theoretical biology 26, 399 (1970).
  * Hillen and Painter (2009) T. Hillen and K. J. Painter, Journal of mathematical biology 58, 183 (2009).
  * Waters and Bassler (2005) C. M. Waters and B. L. Bassler, Annu. Rev. Cell Dev. Biol. 21, 319 (2005).
  * Ocko and Mahadevan (2015) S. A. Ocko and L. Mahadevan, Physical review letters 114, 134501 (2015).
  * Jost  _et al._ (2007) C. Jost, J. Verret, E. Casellas, J. Gautrais, M. Challet, J. Lluc, S. Blanco, M. J. Clifton, and G. Theraulaz, Journal of the Royal Society Interface 4, 107 (2007).
  * Mathijssen  _et al._ (2019) A. J. Mathijssen, J. Culver, M. S. Bhamla, and M. Prakash, Nature 571, 560 (2019).
  * Peleg  _et al._ (2018) O. Peleg, J. M. Peters, M. K. Salcedo, and L. Mahadevan, Nature Physics 14, 1193 (2018).
  * Nguyen  _et al._ (2021) D. M. T. Nguyen, M. L. Iuzzolino, A. Mankel, K. Bozek, G. J. Stephens, and O. Peleg, Proceedings of the National Academy of Sciences 118 (2021).
  * Ranch (2013) G. B. Ranch, “Bees fanning at hive entrance,”  (2013).
  * Peters  _et al._ (2019) J. M. Peters, O. Peleg, and L. Mahadevan, Journal of the Royal Society Interface 16, 20180561 (2019).
  * Reinhard and Srinivasan (2009) J. Reinhard and M. V. Srinivasan, Food exploitation by social insects: ecological, behavioral, and theoretical approaches 1, 165 (2009).
  * Peters  _et al._ (2017) J. M. Peters, N. Gravish, and S. A. Combes, Journal of Experimental Biology 220, 2203 (2017).
  * Winston (1991) M. L. Winston,  _The biology of the honey bee_ (harvard university press, 1991).


11institutetext: 1 National Astronomical Observatory of Japan, National
Institutes of Natural Sciences, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
2 School of Physics and Astronomy, Faculty of Science, Monash University,
Clayton, Victoria 3800, Australia
3 School of Physics and Astronomy, Cardiff University, The Parade, Cardiff,
CF24 3AA, UK
4 Astronomical Institute, Tohoku University, Sendai 980-8578, Japan
5 Division for the Establishment of Frontier Sciences, Organization for
Advanced Studies, Tohoku University, Sendai 980-8577, Japan
6 INAF-Osservatorio Astronomico di Padova, Via dell’Osservatorio 5, I-35122
Padova, Italy
7 INAF-Osservatorio Astronomico di Capodimonte, Via Moiariello 16, I-80131
Napoli, Italy
8 INFN section of Naples, Via Cinthia 6, I-80126, Napoli, Italy
9 International Center for Relativistic Astrophysics, Piazzale della Republica
2, I-65122 Pescara, Italy
10 Department of Physics, Lancaster University, Lancaster, LA1 4YB, UK
11 Department of Physics and Astronomy, Vesilinnantie 5, 20014 University of
Turku, Finland
12 Department of Physics ”E. Pancini”, University Federico II, Via Cinthia 6,
I-80126, Napoli, Italy
13 INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125,
Firenze, Italy
14 Aix-Marseille Univ, CNRS/IN2P3, CPPM, Marseille, France
15 ESAC/ESA, Camino Bajo del Castillo, s/n., Urb. Villafranca del Castillo,
28692 Villanueva de la Cañada, Madrid, Spain
16 Institute of Cosmology and Gravitation, University of Portsmouth,
Portsmouth PO1 3FX, UK
17 INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via
Piero Gobetti 93/3, I-40129 Bologna, Italy
18 INAF-Osservatorio Astrofisico di Torino, Via Osservatorio 20, I-10025 Pino
Torinese (TO), Italy
19 INFN-Sezione di Roma Tre, Via della Vasca Navale 84, I-00146, Roma, Italy
20 Dipartimento di Fisica, Universitá degli studi di Genova, and INFN-Sezione
di Genova, via Dodecaneso 33, I-16146, Genova, Italy
21 Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP,
Rua das Estrelas, PT4150-762 Porto, Portugal
22 Dipartimento di Fisica, Universitá degli Studi di Torino, Via P. Giuria 1,
I-10125 Torino, Italy
23 INFN-Sezione di Torino, Via P. Giuria 1, I-10125 Torino, Italy
24 INAF-IASF Milano, Via Alfonso Corti 12, I-20133 Milano, Italy
25 Port d’Informació Científica, Campus UAB, C. Albareda s/n, 08193 Bellaterra
(Barcelona), Spain
26 Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of
Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona), Spain
27 INAF-Osservatorio Astronomico di Roma, Via Frascati 33, I-00078 Monteporzio
Catone, Italy
28 Dipartimento di Fisica e Astronomia ”Augusto Righi” - Alma Mater Studiorum
Università di Bologna, via Piero Gobetti 93/2, I-40129 Bologna, Italy
29 Institut national de physique nucléaire et de physique des particules, 3
rue Michel-Ange, 75794 Paris Cédex 16, France
30 Centre National d’Etudes Spatiales, Toulouse, France
31 Institute for Astronomy, University of Edinburgh, Royal Observatory,
Blackford Hill, Edinburgh EH9 3HJ, UK
32 Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy,
University of Manchester, Oxford Road, Manchester M13 9PL, UK
33 European Space Agency/ESRIN, Largo Galileo Galilei 1, 00044 Frascati, Roma,
Italy
34 Univ Lyon, Univ Claude Bernard Lyon 1, CNRS/IN2P3, IP2I Lyon, UMR 5822,
F-69622, Villeurbanne, France
35 Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique
Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix,
Switzerland
36 Mullard Space Science Laboratory, University College London, Holmbury St
Mary, Dorking, Surrey RH5 6NT, UK
37 Departamento de Física, Faculdade de Ciências, Universidade de Lisboa,
Edifício C8, Campo Grande, PT1749-016 Lisboa, Portugal
38 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências,
Universidade de Lisboa, Campo Grande, PT-1749-016 Lisboa, Portugal
39 Department of Astronomy, University of Geneva, ch. dÉcogia 16, CH-1290
Versoix, Switzerland
40 Université Paris-Saclay, CNRS, Institut d’astrophysique spatiale, 91405,
Orsay, France
41 Department of Physics, Oxford University, Keble Road, Oxford OX1 3RH, UK
42 INFN-Padova, Via Marzolo 8, I-35131 Padova, Italy
43 Université Paris-Saclay, Université Paris Cité, CEA, CNRS, Astrophysique,
Instrumentation et Modélisation Paris-Saclay, 91191 Gif-sur-Yvette, France
44 INAF-Osservatorio Astronomico di Trieste, Via G. B. Tiepolo 11, I-34143
Trieste, Italy
45 Istituto Nazionale di Astrofisica (INAF) - Osservatorio di Astrofisica e
Scienza dello Spazio (OAS), Via Gobetti 93/3, I-40127 Bologna, Italy
46 Istituto Nazionale di Fisica Nucleare, Sezione di Bologna, Via Irnerio 46,
I-40126 Bologna, Italy
47 Max Planck Institute for Extraterrestrial Physics, Giessenbachstr. 1,
D-85748 Garching, Germany
48 Universitäts-Sternwarte München, Fakultät für Physik, Ludwig-Maximilians-
Universität München, Scheinerstrasse 1, 81679 München, Germany
49 Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029
Blindern, N-0315 Oslo, Norway
50 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak
Grove Drive, Pasadena, CA, 91109, USA
51 von Hoerner & Sulger GmbH, SchloßPlatz 8, D-68723 Schwetzingen, Germany
52 Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark
53 Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg,
Germany
54 Department of Physics and Helsinki Institute of Physics, Gustaf Hällströmin
katu 2, 00014 University of Helsinki, Finland
55 NOVA optical infrared instrumentation group at ASTRON, Oude Hoogeveensedijk
4, 7991PD, Dwingeloo, The Netherlands
56 Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71,
53121 Bonn, Germany
57 INFN-Sezione di Bologna, Viale Berti Pichat 6/2, I-40127 Bologna, Italy
58 Department of Physics, Institute for Computational Cosmology, Durham
University, South Road, DH1 3LE, UK
59 Institut d’Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université,
98 bis boulevard Arago, 75014 Paris, France
60 University of Applied Sciences and Arts of Northwestern Switzerland, School
of Engineering, 5210 Windisch, Switzerland
61 European Space Agency/ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The
Netherlands
62 Department of Physics and Astronomy, University of Aarhus, Ny Munkegade
120, DK-8000 Aarhus C, Denmark
63 Institute of Space Science, Bucharest, Ro-077125, Romania
64 IFPU, Institute for Fundamental Physics of the Universe, via Beirut 2,
34151 Trieste, Italy
65 Dipartimento di Fisica e Astronomia ”G.Galilei”, Universitá di Padova, Via
Marzolo 8, I-35131 Padova, Italy
66 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas
(CIEMAT), Avenida Complutense 40, 28040 Madrid, Spain
67 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências,
Universidade de Lisboa, Tapada da Ajuda, PT-1349-018 Lisboa, Portugal
68 Universidad Politécnica de Cartagena, Departamento de Electrónica y
Tecnología de Computadoras, 30202 Cartagena, Spain
69 Infrared Processing and Analysis Center, California Institute of
Technology, Pasadena, CA 91125, USA
70 INAF-Osservatorio Astronomico di Brera, Via Brera 28, I-20122 Milano, Italy
71 Junia, EPA department, F 59000 Lille, France
72 California institute of Technology, 1200 E California Blvd, Pasadena, CA
91125, USA

#  Euclid: Searching for pair-instability supernovae with the Deep
Survey††thanks: This paper is published on behalf of the Euclid Consortium.

T.J. Moriya1,2<EMAIL_ADDRESS>C. Inserra3 M. Tanaka4,5 E.
Cappellaro6 M. Della Valle7,8,9 I. Hook10 R. Kotak11 G. Longo12,8 F.
Mannucci13 S. Mattila11 C. Tao14 B. Altieri15 A. Amara16 N. Auricchio17 D.
Bonino18 E. Branchini19,20 M. Brescia7 J. Brinchmann21 S. Camera22,23,18 V.
Capobianco18 C. Carbone24 J. Carretero25,26 M. Castellano27 S. Cavuoti7,8,12
A. Cimatti28,13 R. Cledassou29,30 G. Congedo31 C.J. Conselice32 L.
Conversi33,15 Y. Copin34 L. Corcione18 F. Courbin35 M. Cropper36 A. Da
Silva37,38 H. Degaudenzi39 M. Douspis40 F. Dubath39 C.A.J. Duncan41,32 X.
Dupac15 S. Dusini42 A. Ealet34 S. Farrens43 S. Ferriol34 M. Frailis44 E.
Franceschi17 M. Fumana24 B. Garilli24 W. Gillard14 B. Gillis31 C. Giocoli45,46
A. Grazian6 F. Grupp47,48 S.V.H. Haugan49 W. Holmes50 F. Hormuth51 A.
Hornstrup52 K. Jahnke53 S. Kermiche14 A. Kiessling50 M. Kilbinger43 T.
Kitching36 H. Kurki-Suonio54 S. Ligori18 P. B. Lilje49 I. Lloro55 E.
Maiorano17 O. Mansutti44 O. Marggraf56 K. Markovic50 F. Marulli28,17,57 R.
Massey58 H.J. McCracken59 M. Melchior60 M. Meneghetti17,57 G. Meylan35 M.
Moresco28,17 L. Moscardini28,17,57 E. Munari44 S.M. Niemi61 C. Padilla26 S.
Paltani39 F. Pasian44 K. Pedersen62 V. Pettorino43 M. Poncet30 L. Popa63 F.
Raison47 J. Rhodes50 G. Riccio7 E. Rossetti28 R. Saglia47,48 B. Sartoris44,64
P. Schneider56 A. Secroun14 G. Seidel53 C. Sirignano65,42 G. Sirri57 L.
Stanco42 P. Tallada-Crespí66,25 A.N. Taylor31 I. Tereno37,67 R. Toledo-Moreo68
F. Torradeflot66,25 Y. Wang69 G. Zamorani17 J. Zoubian14 S. Andreon70 V.
Scottez59,71 P.W. Morris72

(Received April 19, 2022; accepted August 24, 2022)

Pair-instability supernovae are theorized supernovae that have not yet been
observationally confirmed. They are predicted to exist in low-metallicity
environments. Because overall metallicity becomes lower at higher redshifts,
deep near-infrared transient surveys probing high-redshift supernovae are
suitable to discover pair-instability supernovae. The Euclid satellite, which
is planned to be launched in 2023, has a near-infrared wide-field instrument
that is suitable for a high-redshift supernova survey. The Euclid Deep Survey
is planned to make regular observations of three Euclid Deep Fields
($40\leavevmode\nobreak\ \mathrm{deg^{2}}$ in total) spanning the Euclid’s 6
year primary mission period. While the observations of the Euclid Deep Fields
are not frequent, we show that the predicted long duration of pair-instability
supernovae would allow us to search for high-redshift pair-instability
supernovae with the Euclid Deep Survey. Based on the current observational
plan of the Euclid mission, we conduct survey simulations in order to estimate
the expected numbers of pair-instability supernova discoveries. We find that
up to several hundred pair-instability supernovae at $z\lesssim 3.5$ can be
discovered within the Euclid Deep Survey. We also show that pair-instability
supernova candidates can be efficiently identified by their duration and color
that can be determined with the current Euclid Deep Survey plan. We conclude
that the Euclid mission can lead to the first confirmation of pair-instability
supernovae if their event rates are as high as those predicted by recent
theoretical studies. We also update the expected numbers of superluminous
supernova discoveries in the Euclid Deep Survey based on the latest
observational plan.

###### Key Words.:

stars: massive – supernovae: general – surveys

##  1 Introduction

Optical transient surveys have revealed the existence of huge diversity in
supernovae (SNe) and many new kinds of SNe that were not recognized before
have been discovered (Inserra, 2019). For example, SNe that are often more
than 10 times more luminous than other SNe, known as superluminous SNe (SLSNe,
Gal-Yam 2019 for a recent review), were discovered through unbiased optical
transient surveys (Quimby et al., 2011). Just after the existence of SLSNe was
recognized, they were often theoretically linked to pair-instability SNe
(PISNe, e.g., Smith et al. 2007; Gal-Yam et al. 2009). PISNe are theorized
thermonuclear explosions of very massive stars (Barkat et al., 1967; Rakavy &
Shaviv, 1967) having helium core mass between $\simeq 65\leavevmode\nobreak\
\mathrm{M}_{\odot}$ and $\simeq 135\leavevmode\nobreak\ \mathrm{M}_{\odot}$
(e.g., Heger & Woosley, 2002). Some PISNe are predicted to produce a large
amount of radioactive ${}^{56}\mathrm{Ni}$, up to $\simeq
60\leavevmode\nobreak\ \mathrm{M}_{\odot}$, that heats SN ejecta and they can
be as bright as SLSNe (Kasen et al., 2011; Dessart et al., 2013; Gilmer et
al., 2017). However, the light-curve and spectral properties of SLSNe are now
found to be different from those predicted for PISNe (e.g., Dessart et al.,
2012; Nicholl et al., 2013; Jerkstrand et al., 2016; Inserra et al., 2017;
Mazzali et al., 2019). Therefore, it is generally believed that most SLSNe are
not PISNe and they are powered by a different mechanism (see Moriya et al.
2018 for a review on SLSN powering mechanisms). No SNe have been confirmed to
be PISNe to date, although there are several possible PISN candidates (e.g.,
Lunnan et al., 2016; Terreran et al., 2017) and some possible indirect
signatures of PISNe (e.g., Aoki et al., 2014; Isobe et al., 2022). Confirming
the existence of PISNe is a critical test of stellar evolution theory (Langer,
2012). PISNe can play an important role in chemical evolution (e.g., Nomoto et
al., 2013). PISNe are also essential components for the proper understanding
of black hole mass distributions observed by gravitational waves (e.g.,
Marchant et al., 2016; Farmer et al., 2019).

Table 1: List of SN templates used in this study.

Name | Explosion | Progenitor | Mass
---|---|---|---
R250aa$a$Kasen et al. (2011) | PISN | RSG | 250 $\mathrm{M}_{\odot}$
R225aa$a$Kasen et al. (2011) | PISN | RSG | 225 $\mathrm{M}_{\odot}$
R200aa$a$Kasen et al. (2011) | PISN | RSG | 200 $\mathrm{M}_{\odot}$
R175aa$a$Kasen et al. (2011) | PISN | RSG | 175 $\mathrm{M}_{\odot}$
R150aa$a$Kasen et al. (2011) | PISN | RSG | 150 $\mathrm{M}_{\odot}$
He130aa$a$Kasen et al. (2011) | PISN | WR | 130 $\mathrm{M}_{\odot}$
He120aa$a$Kasen et al. (2011) | PISN | WR | 120 $\mathrm{M}_{\odot}$
He110aa$a$Kasen et al. (2011) | PISN | WR | 110 $\mathrm{M}_{\odot}$
He100aa$a$Kasen et al. (2011) | PISN | WR | 100 $\mathrm{M}_{\odot}$
He090aa$a$Kasen et al. (2011) | PISN | WR | 90 $\mathrm{M}_{\odot}$
SLSNbb$b$Moriya et al. (2022) | SLSN | - | -
111

Figure 1: Normalized Euclid filter transmission (Euclid Collaboration:
Scaramella et al. 2022; Euclid Collaboration: Schirmer et al. 2022). Figure 2:
Rise times and peak absolute magnitudes of the PISN and SLSN models in the
rest frame in the $I_{\scriptscriptstyle\rm E}$ band. The PISN models (R250,
R225, R200, R175, R150, He130, He120, He110, He100, and He090) are from Kasen
et al. (2011) and the SLSN template is from Moriya et al. (2022). We also show
rise times and peak absolute magnitudes of the SN Ia template from Hsiao et
al. (2007) and the SN II Nugent template for comparison (see Sect. 5 for
details). The error bars in the SLSN, SN Ia, and SN II templates show a
typical dispersion (e.g., Richardson et al., 2014; Miller et al., 2020;
González-Gaitán et al., 2015).  Figure 3: $H_{\scriptscriptstyle\rm E}$-band
light curves of the R250 PISN model, the He130 PISN model, and the SLSN
template at $z=2$ and 3. The horizontal dashed line at
$H_{\scriptscriptstyle\rm E}=24.0\leavevmode\nobreak\ \mathrm{mag}$ shows the
$5\sigma$ limiting magnitude of each visit planned in the EDS observations.

Figure 4:  Peak magnitudes of the PISN models. The point-source $5\sigma$
limiting magnitudes of each visit adopted in our survey simulations (25.5 mag
for the $I_{\scriptscriptstyle\rm E}$ band and 24.0 mag for the
$Y_{\scriptscriptstyle\rm E}$, $J_{\scriptscriptstyle\rm E}$, and
$H_{\scriptscriptstyle\rm E}$ bands in AB magnitudes) are shown with the
dashed lines. As discussed in Section 3.1, Euclid may actually have slightly
better limiting magnitudes.

PISNe are believed to appear in low-metallicity environments below
$\mathrm{Z_{\odot}}/3$ (Langer et al. 2007, but see also Georgy et al. 2017).
Thus, they may preferentially appear at high redshifts where metallicity is
low (e.g., Chruslinska & Nelemans, 2019). However, it is difficult to discover
high-redshift PISNe with optical transient surveys as they become faint in the
optical wavelength range because of redshift. For example, a deep optical
transient survey with Subaru/Hyper Suprime-Cam (HSC) carried out a search for
high-redshift PISN candidates, but no candidates have been found so far
(Moriya et al., 2021). The future Vera C. Rubin Observatory/Legacy Survey of
Space and Time (LSST) may be able to discover PISNe depending on its survey
strategy that is not determined yet (e.g., du Buisson et al., 2020). A deep
transient survey in the near-infrared (NIR) is observationally better suited
to discover PISNe appearing at high redshifts (Scannapieco et al., 2005; Pan
et al., 2012; Whalen et al., 2013; Moriya et al., 2019, 2022; Wong et al.,
2019; Regős et al., 2020).

There are several wide-field NIR instruments planned in 2020s such as the Wide
Field Instrument on Nancy Grace Roman Space Telescope from 2026 (Spergel et
al., 2015) and the ULTIMATE Wide Field Instrument on Subaru from
2027222https://ultimate.naoj.org/. Among them, Euclid is planned to be
launched in 2023 (Laureijs et al., 2011). It has wide-field
($0.53\leavevmode\nobreak\ \mathrm{deg^{2}}$) imagers, one in the visual range
(with a passband that we denote $I_{\scriptscriptstyle\rm E}$) and the other
in the NIR range with three filters ($Y_{\scriptscriptstyle\rm E}$,
$J_{\scriptscriptstyle\rm E}$, and $H_{\scriptscriptstyle\rm E}$). The Euclid
mission plans wide ($15\,000\leavevmode\nobreak\ \mathrm{deg^{2}}$) and deep
($40\leavevmode\nobreak\ \mathrm{deg^{2}}$) surveys in the four filters
(Euclid Collaboration: Scaramella et al., 2022). As introduced below, the
Euclid Deep Survey (EDS) covering $40\leavevmode\nobreak\ \mathrm{deg^{2}}$ in
total is planned to visit the survey fields regularly during the six-year
primary survey period of the Euclid mission (Laureijs et al., 2011; Euclid
Collaboration: Scaramella et al., 2022). Although the current EDS plan is not
suitable for searching for short-duration Type Ia SNe (SNe Ia, Astier et al.
2014), Inserra et al. (2018a) previously showed that the regular observations
planned in EDS allow us to search for high-redshift SLSNe up to $z\simeq 3.5$.

High-redshift PISNe are expected to be bright in the NIR and have long
duration (Sect. 2.1). Therefore, it is possible that high-redshift PISNe can
be identified even with the long cadence observations planned in EDS (Sect.
3.1). In this paper, we perform transient survey simulations in order to
quantify the expected numbers of PISN discoveries in EDS. We also investigate
the expected numbers of SLSN discoveries in EDS, because the EDS observational
plan has been updated since the publication of Inserra et al. (2018a).

The rest of this paper is organized as follows. We first introduce the light
curve templates of PISNe and SLSNe we use in our survey simulations in Sect.
2. We then introduce the assumptions in our survey simulations in Sect. 3. The
survey simulation results are presented in Sect. 4. We discuss how to identify
high-redshift PISN and SLSN candidates efficiently with photometry in Sect. 5.
Section 6 concludes this paper. All magnitudes presented in this paper are in
the AB system. We adopt the $\Lambda$CDM cosmology with
$H_{0}=70\leavevmode\nobreak\ \mathrm{km\leavevmode\nobreak\
s^{-1}\leavevmode\nobreak\ Mpc^{-1}}$, $\Omega_{\Lambda}=0.7$, and
$\Omega_{\mathrm{m}}=0.3$.

##  2 SN templates

SN templates are required to estimate the observational properties at high
redshifts. The SN templates we use in this study is introduced in this
section. They are summarized in Table 1.

###  2.1 PISNe

We use the synthetic spectral evolution of PISNe obtained by Kasen et al.
(2011), who took the PISN explosion models from Heger & Woosley (2002) and
calculated their spectral evolution by performing radiative transfer
calculations. Kasen et al. (2011) calculated PISN properties for red
supergiant (RSG), blue supergiant (BSG), and Wolf-Rayet (WR) progenitors. The
RSG and BSG PISN progenitors have hydrogen-rich envelopes, while WR PISN
progenitors are hydrogen-free. The BSG PISN progenitors are expected to exist
mainly at zero metallicity (Heger & Woosley, 2002). The highest redshift PISNe
that can be discovered by the Euclid survey, which is determined by the
maximum expected luminosity of PISNe, is $z\simeq 3.5$ as we show later and
most of them are not likely to have zero metallicity. Thus, a significant
fraction of PISNe discovered by Euclid would be from RSG and WR progenitors
and we adopt RSG and WR PISN models in this study. The RSG PISN models have a
metallicity of $10^{-4}\leavevmode\nobreak\ \mathrm{Z_{\odot}}$. Although the
metallicity of the host galaxies in the redshift range we can probe with EDS
is higher than $10^{-4}\leavevmode\nobreak\ \mathrm{Z_{\odot}}$ (e.g., Steidel
et al., 2014; Chruslinska & Nelemans, 2019), their light-curve properties are
not significantly different from those with high metallicity (e.g., Kozyreva
et al., 2014). We also take PISNe with WR progenitors into account in this
study. The WR PISN models in Kasen et al. (2011) are calculated with the zero-
metallicity condition (Heger & Woosley, 2002), but WR PISNe can appear with
finite metallicity (e.g., du Buisson et al., 2020). The light-curve properties
of zero-metallicity WR PISN models are not significantly different from those
of higher metallicity WR PISNe up to the proposed PISN metallicity threshold
of around $\mathrm{Z_{\odot}}/3$ (T.J. Moriya et al. in prep.). Thus, we adopt
the WR PISN properties from the zero metallicity progenitors in this study. We
note that the adopted progenitor models do not take stellar rotation into
account (e.g., Chatzopoulos & Wheeler, 2012). In addition, only the wind mass
loss is considered.

Kasen et al. (2011) provide the RSG PISN explosion models for 250
$\mathrm{M}_{\odot}$ (R250), 225 $\mathrm{M}_{\odot}$ (R225), 200
$\mathrm{M}_{\odot}$ (R200), 175 $\mathrm{M}_{\odot}$ (R175), and 150
$\mathrm{M}_{\odot}$ (R150) stars at the zero-age main sequence and we
investigate them all. We adopt the WR PISN models whose progenitor evolution
is calculated from the helium zero-age main sequence with a helium core mass
of 130 $\mathrm{M}_{\odot}$ (He130), 120 $\mathrm{M}_{\odot}$ (He120), 110
$\mathrm{M}_{\odot}$ (He110), 100 $\mathrm{M}_{\odot}$ (He100), and 90
$\mathrm{M}_{\odot}$ (He090). They are the five brightest WR PISNe presented
in Kasen et al. (2011) whereas the WR PISN models below 90
$\mathrm{M}_{\odot}$ are too faint to be observed in our redshift range of
interest (see Sect. 3.1).

In order to estimate the PISN properties at high redshifts, we take the rest-
frame synthetic spectra, shift them to an assumed redshift, and convolve the
shifted observer-frame spectra with the Euclid filter transmission (Fig. 1).
The rise times and peak absolute magnitudes of the PISN models in the rest
frame are summarized in Fig. 2. Figure 3 shows the $H_{\scriptscriptstyle\rm
E}$-band light curves of the brightest RSG PISN model (R250) and the brightest
WR PISN model (He130) at $z=2$ and 3. Figure 4 shows the peak magnitudes of
the PISN models at $0.5<z<4$ within the four Euclid filters.

We do not consider Ly$\alpha$ absorption by the intergalactic medium (IGM) and
circumgalactic medium (CGM) in this study. The IGM and CGM absorption affects
the brightness in the $I_{\scriptscriptstyle\rm E}$ band for SNe mainly at
$z\gtrsim 3$. PISNe at $z\gtrsim 3$ are fainter than the
$I_{\scriptscriptstyle\rm E}$-band detection limit (25.5 mag) except for the
R250 PISN model (Fig. 4). In addition, the maximum redshift up to which we can
discover R250 PISNe with the assumed cadence is estimated to be $z\simeq 2.5$,
as shown in the following sections. Thus, the Ly$\alpha$ absorption does not
affect our conclusions in this study.

Figure 5: Peak magnitudes of the SLSN template. The limiting magnitudes of
each visit by the EDS observations (25.5 mag for the $I_{\scriptscriptstyle\rm
E}$ band and 24.0 mag for the $Y_{\scriptscriptstyle\rm E}$,
$J_{\scriptscriptstyle\rm E}$, and $H_{\scriptscriptstyle\rm E}$ bands in AB
magnitudes) are shown with the dashed lines. The template magnitudes are
randomly scattered by the standard deviation of 0.7 mag (both to brighter and
fainter magnitudes) assuming the Gaussian probability distribution in the
survey simulations to match the observed luminosity function as discussed in
Sect. 3.

###  2.2 SLSNe

In this study, we only consider hydrogen-poor (Type I) SLSNe as in our
previous study (Inserra et al., 2018a) and we refer to hydrogen-poor SLSNe
simply as SLSNe. This is because the observational properties of nearby
hydrogen-poor SLSNe have been well investigated and we can reliably estimate
their high-redshift properties. Contrarily, hydrogen-rich (Type II) SLSNe are
known to exist (e.g., Smith et al., 2007; Inserra et al., 2018b), but their
rest-frame ultraviolet properties have not been studied well enough to
estimate their high-redshift properties.

We use the SLSN template presented in Moriya et al. (2022). The SLSN spectral
template is constructed by assuming a broken blackbody function as originally
introduced by Prajs et al. (2017). The template spectra above 300 nm follow
the blackbody function as observed in SLSNe (e.g., Quimby et al., 2018). Below
300 nm, the blackbody function is truncated to match the spectral shape of the
well-observed SLSN iPTF13ajg (Vreeswijk et al., 2014). iPTF13ajg has the
typical shape of SLSN ultraviolet spectra at around peak luminosity, although
some SLSNe are known to be brighter in the ultraviolet (Nicholl et al., 2017).
We only take the spectral shape of iPTF13ajg and the template peak luminosity
is adjusted to match an average value. The template has a peak u-band absolute
magnitude of $-21.1\leavevmode\nobreak\ \mathrm{mag}$, which is an average
SLSN peak brightness estimated by Lunnan et al. (2018). The u-band light-curve
rise time of the template is 30 days.

The rise time and peak absolute magnitude of the SLSN template in the rest
frame are compared to those of the PISN models in Fig. 2. The peak observer-
frame magnitudes of the SLSN template at $0.5<z<4$ obtained by convolving the
Euclid filters are presented in Fig. 5. In the following survey simulations,
the peak SLSN magnitudes are assumed to randomly scatter from the average
value by a standard deviation of 0.7 mag with a Gaussian probability
distribution in order to account for the SLSN luminosity diversity (Lunnan et
al., 2018). While the actual luminosity function of SLSNe may not have a
Gaussian distribution function (e.g., De Cia et al., 2018; Perley et al.,
2020), a Gaussian distribution with 0.7 mag standard deviation covers the wide
peak luminosity range observed in SLSNe. Similarly, the rise times are
independently set to have a standard deviation of 10 days with a Gaussian
probability distribution, but the rise times are not allowed to be shorter
than 10 days. As mentioned in Sect. 2.1, no Ly$\alpha$ absorption in the IGM
and CGM is taken into account because no SLSNe at $z\gtrsim 3$ are brighter
than the detection limit.

Figure 6:  EDF observation plan of the Euclid mission (Euclid Collaboration:
Scaramella et al., 2022). The date is from the beginning of Euclid operation.
The observational epochs at each EDF are shown with squares. The numbers above
the squares show the number of visits within each epoch. For example, 5 visits
for the entire $10\leavevmode\nobreak\ \mathrm{deg^{2}}$ of EDFF will be made
at every epoch shown with squares. The visits within one epoch are made within
15 days. We assume that the images taken in the first 1.25 years of the
mission are used to construct the reference images and the SN search starts
after 1.25 years.

##  3 Survey simulation setups

In this section, we summarize the assumptions in our survey simulations.

###  3.1 Euclid Deep Fields, cadence, and sensitivity

The EDS will cover three Euclid Deep Fields
(EDFs)333https://sci.esa.int/web/euclid/-/61407-the-euclid-wide-and-deep-
surveys. Two fields are planned at the north and south ecliptic poles (EDFN
and EDFS, respectively). The third field is in the Chandra Deep Field South
located in Fornax (EDFF). We introduce the current observational plan at each
field in this section which is also available in Euclid Collaboration:
Scaramella et al. (2022).

EDFN will cover $10\leavevmode\nobreak\ \mathrm{deg^{2}}$ and it will be
visited 10 times at the beginning of the Euclid mission for calibration
purposes. Then, it will be observed regularly roughly every 6 months during
the six-year primary mission period. EDFS will cover $20\leavevmode\nobreak\
\mathrm{deg^{2}}$ and the whole field will be visited 45 times. The first 10
visits will be made at the beginning of the mission and then the subsequent 35
visits will span 5 years. Finally, the $10\leavevmode\nobreak\
\mathrm{deg^{2}}$ area in EDFF will be observed 56 times in total over 6
years.

Figure 6 summarizes the current observational plan for the EDFs. The
observations at each EDF will be conducted at the epochs indicated by squares.
The numbers above show the number of visits planned in each epoch. For
example, EDFF is planned to have 5 visits for the entire
$10\leavevmode\nobreak\ \mathrm{deg^{2}}$ area in every epoch shown with the
squares. The visits within each epoch are made within 15 days. The brightness
of high-redshift PISNe does not change significantly in 15 days (Fig. 3), and
thus for PISNe we regard each square containing several visits as one epoch.
At each visit, images in the $I_{\scriptscriptstyle\rm E}$,
$Y_{\scriptscriptstyle\rm E}$, $J_{\scriptscriptstyle\rm E}$, and
$H_{\scriptscriptstyle\rm E}$ filters will be taken covering the entire survey
field. The $5\sigma$ limiting magnitudes of each visit are estimated to be
$I_{\scriptscriptstyle\rm E}=25.5\leavevmode\nobreak\ \mathrm{mag}$ and
$Y_{\scriptscriptstyle\rm E}=J_{\scriptscriptstyle\rm
E}=H_{\scriptscriptstyle\rm E}=24.0\leavevmode\nobreak\ \mathrm{mag}$ in
Inserra et al. (2018a) and we adopted these limiting magnitudes in this study.
For example, EDFF is planned to have 5 visits with the above $5\sigma$
limiting magnitudes in each epoch. We note that the recently updated estimates
for the limiting magnitudes are better than those we use in this study (Euclid
Collaboration: Scaramella et al., 2022; Guy et al., 2022). In addition, the
limiting magnitudes are for each visit. Because we have multiple visits in one
epoch, we can obtain slightly deeper images by stacking all visits in one
epoch. Meanwhile, when we perform image subtraction for transient
identification, noise will increase and the detection limits can be reduced.
The noise due to image subtraction may not be significant thanks to stable
Euclid point spread functions, although the diffraction-limited NIR point
spread functions are under-sampled. The actual limiting magnitudes will be
revealed when the Euclid observations start.

As shown in Fig. 6, the EDFs will be repeatedly observed with about a half-
year separation. The high-redshift PISNe are expected to be brighter than 24.0
mag in the NIR bands for a year or more (Fig. 3). Thus, the planned
observations every half a year are suitable for searching for high-redshift
PISNe. On the other hand, the typical duration of high-redshift SLSNe is not
long enough for them to be captured by the half-year cadence. However, they
can be still identified because there are multiple (3–5) visits made in a
single epoch. We note that the EDS observation plan has been updated since the
previous study of Inserra et al. (2018a). Some EDFs were planned to have more
frequent visits in the previous plan, which was therefore more suitable for
SLSN discovery.

In order to discover transients on the sky, deep reference images are
required. We assume that deep reference images are created using images taken
during the first 1.25 years of the Euclid mission as indicated in Fig. 6. Many
calibration visits are planned in the first year covering the three EDFs. The
EDFs will be observed 7–13 times, depending on the field, in the first 1.25
years (Euclid Collaboration: Scaramella et al., 2022) and this number of
observations is enough to create reference images of sufficient depth for
subtraction.

Figure 7: Event rates used in the survey simulations. The SLSN event rate is
estimated by setting the SLSN event rate at $z=0.2$ to be
$30\leavevmode\nobreak\ \mathrm{Gpc^{-3}\leavevmode\nobreak\ yr^{-1}}$ (Quimby
et al., 2013; Frohmaier et al., 2021) and extrapolating it following the
cosmic SFR density of Madau & Dickinson (2014). For reference, the SLSN rates
constrained by observations at $z\sim 0.2$ (Frohmaier et al., 2021) and $z\sim
1.1$ (Prajs et al., 2017) are shown. We also use the 10% and 1% SLSN rates for
PISNe. Four theoretical PISN rate estimates by Briel et al. (2022, dotted
lines) are shown for comparison. The four models have different assumptions on
the cosmic SFR history.  Table 2: Expected numbers of PISN and SLSN
discoveries with EDS. The total discovery numbers in the three EDFs are shown.

Model | Rateaa$a$The event rate relative to the SLSN rate presented in Fig. 7. | Detectionbb$b$The minimum number of detected epochs required to be regarded as a discovery. | Total | Redshift Range |
---|---|---|---|---|---
|  | (epochs) |  | 0–0.5 | 0.5–1.0 | 1.0–1.5 | 1.5–2.0 | 2.0–2.5 | 2.5–3.0 | 3.0–3.5 | 3.5–4.0
R250 | 1 | 2 | 145 | 4.5 | 25 | 47 | 49 | 19 | 0.030 | 0 | 0
R250 | 0.1 | 2 | 14 | 0.46 | 2.5 | 4.7 | 4.9 | 1.9 | 0.0030 | 0 | 0
R250 | 0.01 | 2 | 1.5 | 0.045 | 0.27 | 0.52 | 0.49 | 0.19 | 0.0010 | 0 | 0
R225 | 1 | 2 | 158 | 4.4 | 25 | 47 | 53 | 28 | 0 | 0 | 0
R225 | 0.1 | 2 | 16 | 0.46 | 2.6 | 4.8 | 5.2 | 2.7 | 0 | 0 | 0
R225 | 0.01 | 2 | 1.6 | 0.036 | 0.25 | 0.43 | 0.54 | 0.29 | 0 | 0 | 0
R200 | 1 | 2 | 49 | 4.4 | 23 | 21 | 0 | 0 | 0 | 0 | 0
R200 | 0.1 | 2 | 4.9 | 0.47 | 2.4 | 2.1 | 0 | 0 | 0 | 0 | 0
R200 | 0.01 | 2 | 0.52 | 0.056 | 0.24 | 0.23 | 0 | 0 | 0 | 0 | 0
R175 | 1 | 2 | 7.5 | 2.9 | 4.6 | 0 | 0 | 0 | 0 | 0 | 0
R175 | 0.1 | 2 | 0.79 | 0.29 | 0.49 | 0 | 0 | 0 | 0 | 0 | 0
R175 | 0.01 | 2 | 0.081 | 0.038 | 0.044 | 0 | 0 | 0 | 0 | 0 | 0
R150 | 1 | 2 | 3.2 | 1.5 | 1.7 | 0 | 0 | 0 | 0 | 0 | 0
R150 | 0.1 | 2 | 0.32 | 0.13 | 0.19 | 0 | 0 | 0 | 0 | 0 | 0
R150 | 0.01 | 2 | 0.033 | 0.018 | 0.015 | 0 | 0 | 0 | 0 | 0 | 0
He130 | 1 | 2 | 210 | 4.2 | 24 | 47 | 53 | 43 | 28 | 12 | 0.0060
He130 | 0.1 | 2 | 21 | 0.41 | 2.3 | 4.7 | 5.3 | 4.3 | 2.7 | 1.2 | 0
He130 | 0.01 | 2 | 2.1 | 0.040 | 0.26 | 0.46 | 0.51 | 0.45 | 0.26 | 0.12 | 0.00050
He120 | 1 | 2 | 164 | 4.3 | 21 | 47 | 51 | 36 | 3.8 | 0 | 0
He120 | 0.1 | 2 | 16 | 0.46 | 2.1 | 4.6 | 5.1 | 3.6 | 0.37 | 0 | 0
He120 | 0.01 | 2 | 1.6 | 0.048 | 0.20 | 0.43 | 0.47 | 0.38 | 0.037 | 0 | 0
He110 | 1 | 2 | 103 | 4.3 | 20 | 45 | 32 | 2.2 | 0 | 0 | 0
He110 | 0.1 | 2 | 10 | 0.43 | 2.1 | 4.5 | 3.2 | 0.24 | 0 | 0 | 0
He110 | 0.01 | 2 | 1.1 | 0.051 | 0.20 | 0.47 | 0.32 | 0.027 | 0 | 0 | 0
He100 | 1 | 2 | 25 | 4.1 | 16 | 4.7 | 0 | 0 | 0 | 0 | 0
He100 | 0.1 | 2 | 2.4 | 0.40 | 1.6 | 0.45 | 0 | 0 | 0 | 0 | 0
He100 | 0.01 | 2 | 0.24 | 0.049 | 0.16 | 0.038 | 0 | 0 | 0 | 0 | 0
He090 | 1 | 2 | 2.0 | 1.9 | 0.16 | 0 | 0 | 0 | 0 | 0 | 0
He090 | 0.1 | 2 | 0.20 | 0.18 | 0.024 | 0 | 0 | 0 | 0 | 0 | 0
He090 | 0.01 | 2 | 0.027 | 0.027 | 0 | 0 | 0 | 0 | 0 | 0 | 0
SLSN | 1 | 2 | 28 | 2.6 | 8.5 | 8.9 | 5.1 | 2.1 | 0.75 | 0.25 | 0.096
SLSN | 1 | 1 | 121 | 4.9 | 23 | 37 | 31 | 16 | 6.4 | 2.1 | 0.42

444

Figure 8: Cumulative redshift distributions of PISNe discovered in the survey
simulations. The total discovery numbers in the three EDFs are presented. The
top, middle, and bottom panels show the results for the event rates of 100%,
10%, and 1% of the SLSNe rate for each PISN model.

Figure 9:  Redshift distributions of PISNe (top panel) and SLSNe (bottom
panel) discovered in the survey simulations. The total discovery numbers in
the three EDFs are presented. The event rates are set to be the same as the
SLSN rate in all the models in this figure, but the shape of the redshift
distribution is not affected by the assumed rate.

###  3.2 Event rates

Event rates are an essential assumption in the transient survey simulations.
Because no confident PISNe have been discovered so far, the PISN event rate is
unknown. Observationally, the PISN event rate is constrained to be less than
the SLSN event rate (e.g., Nicholl et al., 2013; Moriya et al., 2021). The
SLSN event rate and its redshift evolution have been constrained by optical
transient surveys (Quimby et al., 2013; Prajs et al., 2017; Frohmaier et al.,
2021). The SLSN event rate at $z\simeq 0.2$ is estimated to be
$30\leavevmode\nobreak\ \mathrm{Gpc^{-3}\leavevmode\nobreak\ yr^{-1}}$ (Quimby
et al., 2013). The SLSN event rate up to $z\simeq 1$ is estimated to be
consistent with that obtained by extrapolating the $z\simeq 0.2$ rate
following the cosmic star-formation rate (SFR) density of Madau & Dickinson
(2014, Fig. 7). We currently have a few possible PISN candidates as discussed
in Section 1. The total number of SLSN discoveries is currently around
$100-200$ (e.g., Nicholl, 2021; Chen et al., 2022a, b). Although it is likely
inappropriate to compare the numbers of the PISN candidates and the SLSN
discoveries because of their different properties especially in spectral
energy distributions (Section 2), the current fraction of PISN candidates to
SLSNe should be around a few per cent.

Theoretical predictions for the PISN rates are also very different (see Langer
et al. 2007; Magg et al. 2016; Takahashi 2018; Farmer et al. 2019; du Buisson
et al. 2020; Tanikawa et al. 2022). If we simply assume a Salpeter initial
mass function and a PISN mass range 150–300 $\mathrm{M}_{\odot}$ (Heger &
Woosley, 2002), 1% of core-collapse SNe (8–25 $\mathrm{M}_{\odot}$) would be
PISNe. Given core-collapse SN rates of $10^{5}$–$10^{6}\leavevmode\nobreak\
\mathrm{Gpc^{-3}\leavevmode\nobreak\ yr^{-1}}$ at $z=0$–2.5 (Strolger et al.,
2015), the corresponding PISN rate is $10^{3}$–$10^{4}\leavevmode\nobreak\
\mathrm{Gpc^{-3}\leavevmode\nobreak\ yr^{-1}}$, which is much higher than the
aforementioned SLSN rate. Recently, Briel et al. (2022) provided theoretical
predictions for the PISN rates. By using several cosmic SFR densities from
both observations and cosmological simulations, they predicted PISN rates
based on stellar evolution models from Binary Population and Spectral
Synthesis (BPASS, Eldridge et al. 2017; Stanway & Eldridge 2018). The PISN
rate predictions at $z\lesssim 3.5$ are shown to be 10–100% of the SLSN rate
(Fig. 7).

Although the PISN event rate is unknown, we need to assume some PISN event
rates for our survey simulations. We think it is reasonable to assume PISN
rates that are scaled with the SLSN event rate for our purpose. The SLSN event
rate is constrained observationally and we do not need to rely on uncertain
theoretical rate predictions in this way. A fraction of SLSNe may be PISNe
even though we do not have a conclusive case. In addition, the latest
theoretical study by Briel et al. (2022) predicts that the PISN rate at
$z\lesssim 3.5$ is 10–100% of the SLSN rate (Fig. 7). In light of the above
and given the large uncertainties involved in this type of estimates, below we
make three assumptions for the PISN event rate: (i) the same as the SLSN event
rate, (ii) 10% of the SLSN event rate, and (iii) 1% of the SLSN event rate.

NIR and radio surveys for core-collapse SNe in the local Universe have shown
that a significant fraction of the events are lost from current optical
surveys as they explode in very dusty environments (e.g., Mannucci et al.,
2007; Pérez-Torres et al., 2009; Mattila et al., 2012; Kool et al., 2018;
Jencson et al., 2019; Kankare et al., 2021; Fox et al., 2021). In particular,
the fraction of star-formation activity that takes place in luminous
starbursting and dusty galaxies increases sharply with redshift and becomes
the dominant component of star formation at $z\simeq 1$. Consequently, an
increasing fraction of SNe is expected to be lost from high-redshift optical
searches. Up to 60% of core-collapse SNe are likely to be lost at $z\simeq 2$
(Mannucci et al., 2007). Since the SLSN rates discussed in this section are
based on optical surveys, the rates extrapolated to $z\simeq 4$ should be
taken with some caution. They can probably be considered as lower limits to
the actual events rate.

Figure 10: Cumulative redshift distributions of SLSNe discovered in the survey
simulations. The total discovery numbers in the three EDFs are presented.
Figure 11: Distribution of the number of epochs for which luminous PISNe
(R250, R225, R200, He130, He120, and He110) and SLSNe are detected with the
EDFF cadence (Fig. 6). The distribution is obtained by using the whole
detected sample at all redshifts.

###  3.3 Discovery criteria

PISNe and SLSNe are assumed to appear randomly following the assumed event
rates in the survey simulations. Once a PISN or SLSN appears, the expected
magnitudes in each epoch are estimated for the event. If it satisfies the
discovery criteria presented below, the event is regarded as a discovery.

High-redshift PISNe remain bright for a year or more (Fig. 3), while other
common high-redshift SNe such as SNe Ia do not remain bright for such a long
time (Hounsell et al., 2018), as can be seen in Fig. 2. Thus, the transient
duration is important information to identify rare high-redshift PISNe among
many other kinds of high-redshift SNe. We, therefore, count the numbers of
PISNe that become brighter than the limiting magnitude for 2 epochs or more in
any band in our survey simulations. We also investigate the PISN candidate
selection through their color in Sect. 5. Even if we have a detection only in
a single band, non-detection in the other bands can constrain the color of SN
candidates.

High-redshift SLSNe usually have shorter durations than those of PISNe (Figs.
2 and 3) and it is difficult to distinguish them from other high-redshift SNe
by duration with the current EDS plan having a cadence of about a half year
(Fig. 6). However, high-redshift SLSNe can be identified among other short-
duration SNe by their color even with a single epoch when they are brighter
than 24 mag in the NIR bands (Sect. 5). We count the number of SLSNe that
become brighter than the detection limit at the first visit of each epoch in
any band, although many other short-duration SNe will be detected with this
criterion (Sect. 5). For comparison, we also show the number of SLSNe that
exceed the detection limits for two epochs.

##  4 Simulation results

Survey simulations with one set of assumptions are performed $10^{3}$ times
and we obtain the average numbers of discoveries. The survey simulations are
performed in the three EDFs independently with their areas and cadences as
discussed in Section 3.1. However, there are no significant differences in the
redshift distributions of the discovered PISNe and SLSNe among the three EDFs.
Therefore, we show and discuss the combined numbers of PISN and SLSN
discoveries in the three fields. The total numbers of PISN and SLSN
discoveries in the three EDFs are summarized in Table 2.

Figure 8 shows the cumulative redshift distributions of PISNe discovered in
our survey simulations. We show the distributions for the three different
assumptions for the event rates of the PISN models (100%, 10%, and 1% of the
SLSN rate), and the expected discovery numbers are roughly proportional to the
assumed rate. The brightest PISN models such as R250 and He130 can be
discovered with EDS if their rates are more than about 1% of the SLSN rate.
Meanwhile, the faintest PISN models in our survey simulations such as R150 and
He090 can only be discovered if their event rates are as high as that of
SLSNe. Figure 9 shows the redshift distributions of the discovered PISNe for
the same event rates as the SLSN event rate. The shapes of the redshift
distributions do not depend on the assumed rate. The brightest PISN model
He130 can be discovered up to $z\simeq 3.5$. The brightest RSG PISN models
(R250 and R225) can be discovered up to $z\simeq 2.5$. The faintest PISN
models (R150 and He090) can be discovered only up to $z\simeq 0.6$–0.8 and,
therefore, their expected discovery numbers are relatively small.

The expected cumulative number distributions for SLSNe are presented in Fig.
10 and their number density distributions are shown in Fig. 9. As we discussed
previously, the duration of high-redshift SLSNe are much shorter than that of
PISNe. Figure 11 shows the distribution of the number of detection epochs in
our survey simulations. Only about 20% of high-redshift SLSNe can be
identified if we require them to be detected at multiple epochs. 28 SLSNe are
expected to be detected at two epochs, while 121 SLSNe are discovered if SLSNe
with only a single epoch detection are included. High-redshift SLSNe could be
distinguished from high-redshift SNe Ia and II through color as we discuss in
Sect. 5. In addition, we have multiple visits at epoch and the luminosity
evolution in individual visits can also be used to distinguish high-redshift
SLSNe from high-redshift SNe Ia and II.

Although we lose a large fraction of high-redshift SLSNe if we take only those
detected at multiple epochs, we still have 28 SLSNe with multiple epoch
detection. Because of the unknown PISN rates, these SLSNe could be
contaminants. The color of PISNe and SLSNe can be similar and it is difficult
to distinguish high-redshift PISNe and SLSNe with color alone (Sect. 5).
Because almost no high-redshift SLSNe are expected to be detected at 3 epochs
and many high-redshift PISNe can be brighter than the detection limits for 3
epochs, transients with 3 epoch detection are likely to be genuine high-
redshift PISN candidates.

Finally, no extinction is assumed in our survey simulations. The EDFs have
little Galactic extinction, but the SNe could suffer from severe host galaxy
extinction. However, PISNe and SLSNe preferentially appear in low-metallicity
environments where little extinction is expected (Garn & Best, 2010). Thus,
overall, the effect of extinction may not be significant.

Figure 12: Color-magnitude diagrams in the $I_{\scriptscriptstyle\rm E}$ and
$H_{\scriptscriptstyle\rm E}$ bands for RSG PISNe (top), WR PISNe (middle),
and SLSNe (bottom). PISNe and SLSNe at $z>0.5$ at all phases are plotted with
a redshift interval of 0.1. The color evolution with time can be found in Fig.
13. The regions covered by SNe Ia and SNe II at $z>0.5$ are shown. No
extinction is assumed. The survey limiting magnitudes for the
$I_{\scriptscriptstyle\rm E}$ and $H_{\scriptscriptstyle\rm E}$ bands (25.5
mag and 24.0 mag in AB magnitudes, respectively) are indicated by the dashed
lines. Most SNe Ia and SNe II at $z>0.5$ are too faint to be detected in the
$H_{\scriptscriptstyle\rm E}$ band with the 24.0 mag limit. Figure 13:  Color
evolution of the PISN models and SLSN template presented in Fig. 3.

##  5 PISN and SLSN candidate identification with photometry

Efficient identification of rare PISN and SLSN candidates among many other
abundant SNe such as SNe Ia is critical for triggering spectroscopic follow-up
observations to confirm their discovery. The transient duration is essential
information to identify PISN candidates especially, as discussed in the
previous section. Figure 11 clearly shows that most luminous PISNe are
detected at more than 2 epochs, while other SNe are usually expected to be
bright for one epoch with the half-year cadence observations. Still, it is
important to have as much information as possible to screen the best PISN and
SLSN candidates to make the best use of follow-up resources. In this section,
we discuss the color of PISNe and SLSNe, which can be used to distinguish them
from other common SNe.

Figure 12 shows the color-magnitude diagrams in the $I_{\scriptscriptstyle\rm
E}$ and $H_{\scriptscriptstyle\rm E}$ bands for RSG PISNe at $z>0.5$ (top), WR
PISNe at $z>0.5$ (middle), and SLSNe at $z>0.5$ (bottom). Each point shows
color and magnitude for a certain epoch and all epochs in the models are
presented in the figure. We also present the expected color and magnitude
regions for SNe Ia and SNe II at $z>0.5$, assuming that those at $z<0.5$ can
be discovered by optical transient surveys with, e.g., Subaru/HSC, Blanco/Dark
Energy Camera, or the future Vera C. Rubin Observatory/LSST. A few SNe per
square degree at $z<0.5$ are likely to appear at each epoch (Yasuda et al.,
2019). The expected SN Ia brightness is obtained by using the Hsiao template
(Hsiao et al., 2007) with a peak B-band absolute magnitude of
$-19.26\leavevmode\nobreak\ \mathrm{mag}$ (Richardson et al., 2014), which is
an average absolute magnitude at peak. The SNe II includes SNe IIP, SNe IIL,
and SNe IIn from the Nugent
template555https://c3.lbl.gov/nugent/nugent_templates.html. The peak absolute
magnitudes for SNe II are set as $-18\leavevmode\nobreak\ \mathrm{mag}$ in the
V band in all the subtypes. The assumed SN II peak absolute magnitude is
brighter than the average (Richardson et al., 2014) but the bright SNe II are
preferentially discovered at high redshifts.

Figure 12 shows that high-redshift SNe Ia and SNe II can be as bright as high-
redshift PISNe and SLSNe in the $I_{\scriptscriptstyle\rm E}$ band, but high-
redshift SNe Ia and SNe II do not become as bright as high-redshift PISNe and
SLSNe in the $H_{\scriptscriptstyle\rm E}$ band. Thus, the color between the
$I_{\scriptscriptstyle\rm E}$ and $H_{\scriptscriptstyle\rm E}$ bands allows
us to efficiently distinguish high-redshift PISNe and SLSNe from other SNe.
While we do not assume any extinction in Fig. 12, the fact that SNe Ia and II
are detected in the $I_{\scriptscriptstyle\rm E}$ band but not in the
$H_{\scriptscriptstyle\rm E}$ band indicate that their detections.

Host galaxy information of transients is also an important clue to identify
their redshift information. The reference images taken at the beginning of EDS
as well as deep optical images available in the EDFs would allow us to
estimate host-galaxy photometric redshifts if the host galaxies are bright
enough. Some host-galaxy spectra needed to determine their redshifts could be
simultaneously obtained by the grisms on Euclid. The limiting magnitude of the
reference images in the $I_{\scriptscriptstyle\rm E}$ band are around 26 mag
or deeper. Therefore, host galaxies with absolute magnitudes brighter than
around $-19\leavevmode\nobreak\ \mathrm{mag}$ in the rest-frame ultraviolet
can be identified in the reference images. Many galaxies are fainter than
$-19\leavevmode\nobreak\ \mathrm{mag}$ in the ultraviolet (e.g., Bouwens et
al., 2021) and a large fraction of PISNe could be observed as hostless. Even
if we have the host galaxy photometry, the photometric redshift estimates at
high redshifts are challenging (e.g., Euclid Collaboration: Desprez et al.,
2020). Despite of the difficulties, host galaxy information provides
additional clues to identify the nature of transients. The NIR reference
images will be shallower than that of $I_{\scriptscriptstyle\rm E}$, but they
can trace fainter stellar masses and some host galaxies could only appear in
them.

Active galactic nuclei (AGN) can be variable over long timescales and such a
long-term AGN variability can contaminate our sample of high-redshift PISN and
SLSN candidates. The variable component of unobscured AGN spectra at
ultraviolet to optical wavelengths is estimated to be proportional to
$\nu^{1/3}$, where $\nu$ denotes the frequency (e.g., Kokubo et al., 2014).
Thus, the color of the AGN variability does not significantly change and it
remains at $I_{\scriptscriptstyle\rm E}-H_{\scriptscriptstyle\rm
E}\simeq-0.34\leavevmode\nobreak\ \mathrm{mag}$. Because the color of PISNe
and SLSNe changes much more with time (Fig. 13), the color evolution is
essential information to exclude AGN variability. In order to exclude AGN with
short-timescale variability, simultaneous short cadence optical transient
surveys are helpful. The location of transients in their host galaxies is also
important to exclude AGN and the high spatial resolution of Euclid can help to
constrain the nuclear vs non-nuclear location in host galaxies. Ultimately,
spectroscopic confirmation of the transient would be required to confirm the
classification as PISNe or SLSNe. It is possible that a part of the
spectroscopic follow-up observations will be realized by the grisms on Euclid.
James Webb Space Telescope will also be powerful for the spectroscopic follow-
up observations. Photometric observations in narrower optical bands than
$I_{\scriptscriptstyle\rm E}$ may also be useful to distinguish PISNe and
SLSNe because of the different rest-frame ultraviolet properties expected in
them (e.g., Dessart et al., 2012).

##  6 Conclusions

We have estimated the expected number of PISN and SLSN discoveries in EDS, in
which three EDFs ($40\leavevmode\nobreak\ \mathrm{deg^{2}}$ in total) are
planned to be observed regularly for six years during the Euclid mission.
Although the EDF observations are not frequent (Fig. 6), the intrinsically
long duration of PISNe and SLSNe allows us to discover them in the high-
redshift Universe with the current survey plan. With the currently planned
limiting magnitudes, it is possible to detect PISNe and SLSNe up to $z\simeq
3.5$ (Fig. 4). We find that we can discover up to several hundred PISNe with
EDS if the bright PISN rate is as high as the SLSN rate (Table 2). The
expected number of PISN discoveries exceeds unity as long as the bright PISN
rate is higher than 1% of the SLSN rate. Up to 30–120 SLSNe are expected to be
discovered with the same survey, depending on the discovery criteria (Table
2). We also show that high-redshift PISN and SLSN candidates can be
efficiently identified among abundant SNe such as SNe Ia using their duration,
color, and magnitude (Figs. 2, 11, and 12).

The predicted total PISN rates are around 10–100% of the SLSN rate at
$z\lesssim 3.5$ (Fig. 7). The bright PISNe that come from progenitor masses
above $\simeq 200\leavevmode\nobreak\ \mathrm{M}_{\odot}$ consist of about 50%
of the entire PISN population, assuming a Salpeter initial mass function and a
PISN progenitor mass range of 150–300 $\mathrm{M}_{\odot}$. Even if we take
50% of the estimated number of discoveries from our survey simulations, our
expected number of PISN discoveries exceed unity when the PISN rate is higher
than 10% of the SLSN rate. Therefore, our study indicates that the Euclid
mission provides an unprecedented opportunity to finally discover PISNe if the
latest PISN rate predictions are accurate. If the PISN rate is sufficiently
high, it is also possible that the Euclid mission will lead to the discovery
of many PISNe with which PISNe in different environments and cosmic time can
be studied. Identification and characterization of PISNe at $z\lesssim 3.5$ by
the Euclid mission will allow us to directly confirm the existence of such
objects that have never been identified with certainty before. The PISN
properties can be compared with the theoretical predictions aimed at testing
the theory of the evolution of massive stars in low-metallicity environments.
They can also constrain massive star formation in the early Universe by
providing the formation rate of the massive stars in the PISN mass range. The
PISN discovery by the Euclid mission will be an important base to discover and
characterize PISNe at higher redshifts ($z\gtrsim 3.5$) which can be reached
by future missions such as the Nancy Grace Roman Space Telescope (e.g., Moriya
et al., 2022; Lazar & Bromm, 2022). Even if we do not discover any PISNe, we
will be in a position to obtain a strong constraint on the PISN event rate at
$z\lesssim 3.5$ that will challenge the current understanding of stellar
evolution theory.

###### Acknowledgements.

TJM thanks Mitsuru Kokubo for helpful discussions. TJM is supported by the
Grants-in-Aid for Scientific Research of the Japan Society for the Promotion
of Science (JP18K13585, JP20H00174, JP21K13966, JP21H04997). This work was
supported by JSPS Core-to-Core Program (JPJSCCA20210003). This research has
made use of the SVO Filter Profile Service (http://svo2.cab.inta-
csic.es/theory/fps/) supported from the Spanish MINECO through grant
AYA2017-84089. The Euclid Consortium acknowledges the European Space Agency
and a number of agencies and institutes that have supported the development of
Euclid, in particular the Academy of Finland, the Agenzia Spaziale Italiana,
the Belgian Science Policy, the Canadian Euclid Consortium, the French Centre
National d’Etudes Spatiales, the Deutsches Zentrum für Luft- und Raumfahrt,
the Danish Space Research Institute, the Fundação para a Ciência e a
Tecnologia, the Ministerio de Ciencia e Innovación, the National Aeronautics
and Space Administration, the National Astronomical Observatory of Japan, the
Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the
Romanian Space Agency, the State Secretariat for Education, Research and
Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom
Space Agency. A complete and detailed list is available on the Euclid web site
(http://www.euclid-ec.org).

## References

  * Aoki et al. (2014) Aoki, W., Tominaga, N., Beers, T. C., Honda, S., & Lee, Y. S. 2014, Science, 345, 912
  * Astier et al. (2014) Astier, P., Balland, C., Brescia, M., et al. 2014, A&A, 572, A80
  * Barkat et al. (1967) Barkat, Z., Rakavy, G., & Sack, N. 1967, Phys. Rev. Lett., 18, 379
  * Bouwens et al. (2021) Bouwens, R. J., Oesch, P. A., Stefanon, M., et al. 2021, AJ, 162, 47
  * Briel et al. (2022) Briel, M. M., Eldridge, J. J., Stanway, E. R., Stevance, H. F., & Chrimes, A. A. 2022, MNRAS, 514, 1315
  * Chatzopoulos & Wheeler (2012) Chatzopoulos, E. & Wheeler, J. C. 2012, ApJ, 748, 42
  * Chen et al. (2022a) Chen, Z. H., Yan, L., Kangas, T., et al. 2022a, arXiv e-prints, arXiv:2202.02059
  * Chen et al. (2022b) Chen, Z. H., Yan, L., Kangas, T., et al. 2022b, arXiv e-prints, arXiv:2202.02060
  * Chruslinska & Nelemans (2019) Chruslinska, M. & Nelemans, G. 2019, MNRAS, 488, 5300
  * De Cia et al. (2018) De Cia, A., Gal-Yam, A., Rubin, A., et al. 2018, ApJ, 860, 100
  * Dessart et al. (2012) Dessart, L., Hillier, D. J., Waldman, R., Livne, E., & Blondin, S. 2012, MNRAS, 426, L76
  * Dessart et al. (2013) Dessart, L., Waldman, R., Livne, E., Hillier, D. J., & Blondin, S. 2013, MNRAS, 428, 3227
  * du Buisson et al. (2020) du Buisson, L., Marchant, P., Podsiadlowski, P., et al. 2020, MNRAS, 499, 5941
  * Eldridge et al. (2017) Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34, e058
  * Euclid Collaboration: Desprez et al. (2020) Euclid Collaboration: Desprez, G., Paltani, S., Coupon, J., et al. 2020, A&A, 644, A31
  * Euclid Collaboration: Scaramella et al. (2022) Euclid Collaboration: Scaramella, R., Amiaux, J., Mellier, Y., et al. 2022, A&A, 662, A112
  * Euclid Collaboration: Schirmer et al. (2022) Euclid Collaboration: Schirmer, M., Jahnke, K., Seidel, G., et al. 2022, A&A, 662, A92
  * Farmer et al. (2019) Farmer, R., Renzo, M., de Mink, S. E., Marchant, P., & Justham, S. 2019, ApJ, 887, 53
  * Fox et al. (2021) Fox, O. D., Khandrika, H., Rubin, D., et al. 2021, MNRAS, 506, 4199
  * Frohmaier et al. (2021) Frohmaier, C., Angus, C. R., Vincenzi, M., et al. 2021, MNRAS, 500, 5142
  * Gal-Yam (2019) Gal-Yam, A. 2019, ARA&A, 57, 305
  * Gal-Yam et al. (2009) Gal-Yam, A., Mazzali, P., Ofek, E. O., et al. 2009, Nature, 462, 624
  * Garn & Best (2010) Garn, T. & Best, P. N. 2010, MNRAS, 409, 421
  * Georgy et al. (2017) Georgy, C., Meynet, G., Ekström, S., et al. 2017, A&A, 599, L5
  * Gilmer et al. (2017) Gilmer, M. S., Kozyreva, A., Hirschi, R., Fröhlich, C., & Yusof, N. 2017, ApJ, 846, 100
  * González-Gaitán et al. (2015) González-Gaitán, S., Tominaga, N., Molina, J., et al. 2015, MNRAS, 451, 2212
  * Guy et al. (2022) Guy, L. P., Cuillandre, J.-C., Bachelet, E., et al. 2022, arXiv e-prints, arXiv:2201.03862
  * Heger & Woosley (2002) Heger, A. & Woosley, S. E. 2002, ApJ, 567, 532
  * Hounsell et al. (2018) Hounsell, R., Scolnic, D., Foley, R. J., et al. 2018, ApJ, 867, 23
  * Hsiao et al. (2007) Hsiao, E. Y., Conley, A., Howell, D. A., et al. 2007, ApJ, 663, 1187
  * Inserra (2019) Inserra, C. 2019, Nature Astronomy, 3, 697
  * Inserra et al. (2018a) Inserra, C., Nichol, R. C., Scovacricchi, D., et al. 2018a, A&A, 609, A83
  * Inserra et al. (2017) Inserra, C., Nicholl, M., Chen, T. W., et al. 2017, MNRAS, 468, 4642
  * Inserra et al. (2018b) Inserra, C., Smartt, S. J., Gall, E. E. E., et al. 2018b, MNRAS, 475, 1046
  * Isobe et al. (2022) Isobe, Y., Ouchi, M., Suzuki, A., et al. 2022, ApJ, 925, 111
  * Jencson et al. (2019) Jencson, J. E., Kasliwal, M. M., Adams, S. M., et al. 2019, ApJ, 886, 40
  * Jerkstrand et al. (2016) Jerkstrand, A., Smartt, S. J., & Heger, A. 2016, MNRAS, 455, 3207
  * Kankare et al. (2021) Kankare, E., Efstathiou, A., Kotak, R., et al. 2021, A&A, 649, A134
  * Kasen et al. (2011) Kasen, D., Woosley, S. E., & Heger, A. 2011, ApJ, 734, 102
  * Kokubo et al. (2014) Kokubo, M., Morokuma, T., Minezaki, T., et al. 2014, ApJ, 783, 46
  * Kool et al. (2018) Kool, E. C., Ryder, S., Kankare, E., et al. 2018, MNRAS, 473, 5641
  * Kozyreva et al. (2014) Kozyreva, A., Blinnikov, S., Langer, N., & Yoon, S. C. 2014, A&A, 565, A70
  * Langer (2012) Langer, N. 2012, ARA&A, 50, 107
  * Langer et al. (2007) Langer, N., Norman, C. A., de Koter, A., et al. 2007, A&A, 475, L19
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints, arXiv:1110.3193
  * Lazar & Bromm (2022) Lazar, A. & Bromm, V. 2022, MNRAS, 511, 2505
  * Lunnan et al. (2018) Lunnan, R., Chornock, R., Berger, E., et al. 2018, ApJ, 852, 81
  * Lunnan et al. (2016) Lunnan, R., Chornock, R., Berger, E., et al. 2016, ApJ, 831, 144
  * Madau & Dickinson (2014) Madau, P. & Dickinson, M. 2014, ARA&A, 52, 415
  * Magg et al. (2016) Magg, M., Hartwig, T., Glover, S. C. O., Klessen, R. S., & Whalen, D. J. 2016, MNRAS, 462, 3591
  * Mannucci et al. (2007) Mannucci, F., Della Valle, M., & Panagia, N. 2007, MNRAS, 377, 1229
  * Marchant et al. (2016) Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., & Moriya, T. J. 2016, A&A, 588, A50
  * Mattila et al. (2012) Mattila, S., Dahlen, T., Efstathiou, A., et al. 2012, ApJ, 756, 111
  * Mazzali et al. (2019) Mazzali, P. A., Moriya, T. J., Tanaka, M., & Woosley, S. E. 2019, MNRAS, 484, 3451
  * Miller et al. (2020) Miller, A. A., Yao, Y., Bulla, M., et al. 2020, ApJ, 902, 47
  * Moriya et al. (2021) Moriya, T. J., Jiang, J.-a., Yasuda, N., et al. 2021, ApJ, 908, 249
  * Moriya et al. (2022) Moriya, T. J., Quimby, R. M., & Robertson, B. E. 2022, ApJ, 925, 211
  * Moriya et al. (2018) Moriya, T. J., Sorokina, E. I., & Chevalier, R. A. 2018, Space Sci. Rev., 214, 59
  * Moriya et al. (2019) Moriya, T. J., Wong, K. C., Koyama, Y., et al. 2019, PASJ, 71, 59
  * Nicholl (2021) Nicholl, M. 2021, Astronomy and Geophysics, 62, 5.34
  * Nicholl et al. (2017) Nicholl, M., Berger, E., Margutti, R., et al. 2017, ApJ, 835, L8
  * Nicholl et al. (2013) Nicholl, M., Smartt, S. J., Jerkstrand, A., et al. 2013, Nature, 502, 346
  * Nomoto et al. (2013) Nomoto, K., Kobayashi, C., & Tominaga, N. 2013, ARA&A, 51, 457
  * Pan et al. (2012) Pan, T., Kasen, D., & Loeb, A. 2012, MNRAS, 422, 2701
  * Pérez-Torres et al. (2009) Pérez-Torres, M. A., Romero-Cañizales, C., Alberdi, A., & Polatidis, A. 2009, A&A, 507, L17
  * Perley et al. (2020) Perley, D. A., Fremling, C., Sollerman, J., et al. 2020, ApJ, 904, 35
  * Prajs et al. (2017) Prajs, S., Sullivan, M., Smith, M., et al. 2017, MNRAS, 464, 3568
  * Quimby et al. (2018) Quimby, R. M., De Cia, A., Gal-Yam, A., et al. 2018, ApJ, 855, 2
  * Quimby et al. (2011) Quimby, R. M., Kulkarni, S. R., Kasliwal, M. M., et al. 2011, Nature, 474, 487
  * Quimby et al. (2013) Quimby, R. M., Yuan, F., Akerlof, C., & Wheeler, J. C. 2013, MNRAS, 431, 912
  * Rakavy & Shaviv (1967) Rakavy, G. & Shaviv, G. 1967, ApJ, 148, 803
  * Regős et al. (2020) Regős, E., Vinkó, J., & Ziegler, B. L. 2020, ApJ, 894, 94
  * Richardson et al. (2014) Richardson, D., Jenkins, Robert L., I., Wright, J., & Maddox, L. 2014, AJ, 147, 118
  * Scannapieco et al. (2005) Scannapieco, E., Madau, P., Woosley, S., Heger, A., & Ferrara, A. 2005, ApJ, 633, 1031
  * Smith et al. (2007) Smith, N., Li, W., Foley, R. J., et al. 2007, ApJ, 666, 1116
  * Spergel et al. (2015) Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv e-prints, arXiv:1503.03757
  * Stanway & Eldridge (2018) Stanway, E. R. & Eldridge, J. J. 2018, MNRAS, 479, 75
  * Steidel et al. (2014) Steidel, C. C., Rudie, G. C., Strom, A. L., et al. 2014, ApJ, 795, 165
  * Strolger et al. (2015) Strolger, L.-G., Dahlen, T., Rodney, S. A., et al. 2015, ApJ, 813, 93
  * Takahashi (2018) Takahashi, K. 2018, ApJ, 863, 153
  * Tanikawa et al. (2022) Tanikawa, A., Moriya, T. J., Tominaga, N., & Yoshida, N. 2022, arXiv e-prints, arXiv:2204.09402
  * Terreran et al. (2017) Terreran, G., Pumo, M. L., Chen, T. W., et al. 2017, Nature Astronomy, 1, 713
  * Vreeswijk et al. (2014) Vreeswijk, P. M., Savaglio, S., Gal-Yam, A., et al. 2014, ApJ, 797, 24
  * Whalen et al. (2013) Whalen, D. J., Even, W., Frey, L. H., et al. 2013, ApJ, 777, 110
  * Wong et al. (2019) Wong, K. C., Moriya, T. J., Oguri, M., et al. 2019, PASJ, 71, 60
  * Yasuda et al. (2019) Yasuda, N., Tanaka, M., Tominaga, N., et al. 2019, PASJ, 71, 74