# Supernova explosions in active galactic nuclear discs

Evgeni Grishin1,3,4, Alexey Bobrick2, Ryosuke Hirai3,4, Ilya Mandel3,4,5,
Hagai B. Perets 1
1Physics Department, Technion - Israel Institute of Technology, 3200003,
Haifa, Israel
2 Lund Observatory, Department of Astronomy and Theoretical Physics, Lund,
Sweden
3 Monash Centre for Astrophysics, School of Physics and Astronomy, Monash
University, Clayton, Victoria 3800, Australia
4 OzGrav: Australian Research Council Centre of Excellence for Gravitational
Wave Discovery, Clayton, VIC 3800, Australia
5 Institute of Gravitational Wave Astronomy and School of Physics and
Astronomy, University of Birmingham, Birmingham, B15 2TT, United Kingdom
E-mail<EMAIL_ADDRESS>(EG)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Active galactic nuclei (AGN) are prominent environments for stellar capture,
growth and formation. These environments may catalyze stellar mergers and
explosive transients, such as thermonuclear and core-collapse supernovae
(SNe). SN explosions in AGN discs generate strong shocks, leading to unique
observable signatures. We develop an analytical model which follows the
evolution of the shock propagating in the disc until it eventually breaks out.
We derive the peak luminosity, bolometric lightcurve, and breakout time. The
peak luminosities may exceed $10^{45}$ erg s-1 and last from hours to days.
The brightest explosions occur in regions of reduced density; either off-
plane, or in discs around low-mass central black holes ($\sim 10^{6}\
\textnormal{M}_{\odot}$), or in starved subluminous AGNs. Explosions in the
latter two sites are easier to observe due to a reduced AGN background
luminosity. We perform suites of 1D Lagrangian radiative hydrodynamics SNEC
code simulations to validate our results and obtain the luminosity in
different bands, and 2D axisymmetric Eulerian hydrodynamics code HORMONE
simulations to study the morphology of the ejecta and its deviation from
spherical symmetry. The observed signature is expected to be a bright blue, UV
or X-ray flare on top of the AGN luminosity from the initial shock breakout,
while the subsequent red part of the lightcurve will largely be unobservable.
We estimate the upper limit for the total event rate to be
$\mathcal{R}\lesssim 100\ \rm yr^{-1}\ Gpc^{-3}$ for optimal conditions and
discuss the large uncertainties in this estimate. Future high-cadence
transient searches may reveal these events. Some existing tidal disruption
event candidates may originate from AGN supernovae.

###### keywords:

galaxies: active – transients: supernovae – (stars:) supernovae: general –
hydrodynamics – shock waves – (stars:) circumstellar matter

††pubyear: 2021††pagerange: Supernova explosions in active galactic nuclear
discs–D

##  1 Introduction

The accretion discs of active galactic nuclei (AGN) are responsible for
feeding supermassive black holes (SMBHs) in the centres of galaxies (Lynden-
Bell, 1969). These AGN discs are often described as geometrically thin and
optically thick accretion discs, driven by the $\alpha$-viscosity (Shakura &
Sunyaev, 1973). The rapidly growing observational data on AGN discs has led to
more detailed models of their structure and evolution, either based on the
mass flow inferred from star formation (Thompson et al., 2005), or on the
optical/UV spectral energy distribution (SED) (Sirko & Goodman, 2003).

Galactic nuclei also possess a large number of stars and stellar remnants
(Miralda-Escudé & Gould, 2000; Bartko et al., 2010). Some regions of the discs
might become unstable and give rise to vigorous star-formation (Paczynski,
1978; Dittmann & Miller, 2020), and the interactions of stars and compact
objects with the discs could lead to their capture and embedding into the
discs (Artymowicz et al., 1993). Stars and compact objects in a gaseous
environment could accrete gas and grow over time (Davies & Lin, 2020;
Cantiello et al., 2020), as well as change their orbits due to the interaction
with the disc (Ostriker, 1983; Syer et al., 1991; Artymowicz et al., 1993).
The higher stellar density in the disc, the presence of gas, and the low
relative velocities may enhance the formation of binaries through tidal
effects (Cheng & Wang, 1999) and/or through gas-mediated capture (Goldreich et
al., 2002; Tagawa et al., 2020). The interactions of binary stars with the
gaseous environment could catalyze stellar mergers (Baruteau et al., 2011),
and the formation and gravitational-wave driven mergers of compact-object
binaries (McKernan et al., 2012).

Taken together, AGN disc environments could potentially increase the rate of
SNe, such as core-collapse SNe or type Ia SNe. The rate of core-collapse SNe
may be increased through potentially increased formation rate of massive stars
(Artymowicz et al., 1993). Similarly, the rate of type Ia SNe may be increased
through an increased number of white-dwarf (WD) mergers (McKernan et al.,
2020), or through accretion onto WDs, which can make them reach close-to-
Chandrasekhar masses and explode (Ostriker, 1983). Other explosive mergers of
compact objects with stars or planets could give rise to other transients such
as micro - tidal disruption events (when the stars are disrupted by such
compact objects, Perets et al., 2016).

SNe in AGNs were studied in the context of the feedback on the AGN structure
and evolution (Rozyczka et al., 1995; Ferrara & Tolstoy, 2000; Moranchel-
Basurto et al., 2021). Here we focus on the observational signatures of SN
explosions in AGN discs. The highly relativistic electromagnetic signatures of
$\gamma$-ray bursts (GRBs) and the accretion induced collapse of neutron stars
(NS) into BHs have been recently explored (Zhu et al., 2020; Zhu et al.,
2021a; Perna et al., 2021b; Perna et al., 2021a). Recently, Zhu et al. (2021b)
suggested that type Ia supernovae may potentially be observable in AGN discs.
However, a detailed model for the interaction of a non-relativistic SN
explosion with the AGN disc and its observational signatures is still lacking.
In this paper, we take the first steps in the detailed modelling of general SN
explosions in AGN environments.

Electromagnetic signatures from the interaction between SN ejecta and the
dense circumstellar medium (CSM) have been studied intensively in the context
of type IIn SNe both analytically and numerically (e.g. Woosley et al., 2007;
Chevalier & Irwin, 2011; Ginzburg & Balberg, 2012; Moriya et al., 2013;
Dessart et al., 2015; Tsuna et al., 2019; Takei & Shigeyama, 2020; Suzuki et
al., 2021). Most studies assume a wind-like dense CSM, which have density
profiles following $\rho\propto 1/r^{2}$, where $r$ is the distance from the
progenitor. Some studies have also dealt with asymmetrical CSM structures,
such as disc-like morphologies (Suzuki et al., 2019). In all cases, the main
feature is that the presence of the CSM allows the kinetic energy of the
ejecta to be efficiently dissipated and converted into radiation, increasing
the luminosity of these events. CSM interaction is considered to be the
dominant channel for producing SNe with narrow-line features (type IIn) and
the extremely bright SNe called superluminous SNe (SLSNe; see Gal-Yam, 2019,
for a review on SLSNe).

###  1.1 SNe in AGN discs

Motivated by the recent success of CSM interaction models in explaining SLSNe,
we develop a detailed model for the shock propagation and breakout of an SN
explosion inside an AGN disc. However, there are two main differences between
the wind-driven CSM SNe and our setting. While the origin of the CSM in CSM SN
is the wind itself, which is already expanding prior to the explosion (Nakar &
Piro, 2014; Morozova et al., 2017; Piro et al., 2020), in our case, the local
AGN gas may be assumed to be static at the time of the explosion. The second
difference is that the spatial extent of AGN discs overwhelmingly exceeds the
distance scales of wind-driven explosions. For this reason, as we show
further, explosions in an AGN disc can be analogous to explosions close to the
surface of a much more massive ball of gas (Yalinewich & Matzner, 2019).

The purpose of this paper is to describe the details of SN explosions in dense
environments, such as AGN discs. Here we provide an intuitive picture of the
underlying physical mechanisms and key takeaway messages.

After the explosion, a shock wave is generated and propagates outward.
Initially, the shock propagation velocity is much faster than the photon
diffusion speed in an optically thick environment. As the shock propagates
outward, the medium outside it becomes increasingly optically thin, until the
photons are able to overcome the shock velocity, break out before the shock
and diffuse to the photosphere. This is the breakout shell $z_{\rm bo}$,
formally defined as a solution to the equation $\tau(z_{\rm bo})=c/v(z_{\rm
bo})$, where $v(z_{\rm bo})$ is the shock velocity and $\tau(z_{\rm bo})$ is
the optical depth, all evaluated at $z_{\rm bo}$. At this point, there is a
surge of photons that escape out of the material, causing a sudden rise in the
observable luminosity. After the initial breakout, the photosphere starts to
expand and cool down, which will decrease the luminosity over time.

###  1.2 Overview and structure of the paper

The details of the shock propagation velocity and the breakout depend on the
structure of the medium. In sec. 2, we review the structure and properties of
AGN discs and derive the vertical density profile $\rho(z)$ and the optical
depth $\tau(z)$. In sec. 3, we discuss our assumptions in 3.1, derive and
compare the velocity profile $v(z)$ with numerical simulations in sec. 3.2.
Once $v(z)$ is established, we find the breakout shell $z_{\rm bo}$ and the
breakout time in sec. 3.3. The shock velocity depends on $\rho(z)$ and
$\tau(z)$, and also on the spatial location, and internal properties of the
explosion, namely the explosion energy $E_{0}$ and the ejecta mass $M_{\rm
ej}$.

Once the breakout shell $z_{\rm bo}$ and the photosphere $z_{\rm ph}$ (defined
by $\tau(z_{\rm ph})=1$) are established, the photon diffusion length scale is
$d=z_{\rm ph}-z_{\rm bo}$. Since the velocity $v(z_{\rm bo})$ is also known,
we can estimate the peak luminosity $L_{\rm peak}$ as the energy deposited
onto the breakout shell $e_{\rm bo}\sim\rho(z_{\rm bo})z_{\rm bo}^{2}dv_{\rm
bo}^{2}$ divided by the diffusion photon time $\sim d/v_{\rm bo}$. The
luminosity is then

$L_{\rm peak}=C_{L}\rho(z_{\rm bo})z_{\rm bo}^{2}v_{\rm bo}^{3}.$ (1)

The dimensionless prefactor $C_{L}$ depends on the geometry of the explosion.

The energy deposited in layers beneath the breakout slab cannot efficiently
diffuse out, and the material, therefore, cools mostly through adiabatic
expansion. Adiabatic cooling leads to a rapid decay of the lightcurve:
Initially, the volume is only linearly proportional with time $V\propto t$,
and for a pressure $P$, the change in the energy will be $\int dE=\int
PdV\propto\int V^{-\gamma}dV$. For a radiation dominated medium (which is the
case after it is shocked), we have $\gamma=4/3$, and $\int dE\propto
V^{-1/3}\propto t^{-1/3}$, where we used the connection $PV^{\gamma}=\rm
const$. The luminosity falls off as $L\propto E/t\propto t^{-4/3}$. The peak
luminosity is hence dominated by the initial breakout, which is likely to be
the detectable portion of the lightcurve in the luminous environment of an
AGN. In dense regions where the CSM mass until the breakout is much larger
than the ejecta mass $M_{\rm CSM}\gg M_{\rm ej}$, adiabatic cooling is the
main source of energy loss, while for models with lighter masses $M_{\rm
CSM}\lesssim M_{\rm ej}$ the optical depth is reduced, and radiative losses
will also contribute, and the decay will be steeper.

After some time, the photosphere will expand further out, the volume will
expand spherically, the optical depth will decrease, and deeper layers will be
exposed, which makes the lightcurve decay with a much shallower slope. Nakar &
Sari (2010) found a slope of $\approx-0.35$ for a radiative atmosphere
profile.

The uncertainty of the peak luminosity in terms of the coefficient $C_{L}$ and
the complex structure make it difficult to predict the peak luminosity and
subsequent lightcurve. Recent detailed modelling of explosions close to a
stellar surface has provided analytical expressions which were corroborated
numerically (YM19). We extend the modelling of YM19 to apply it to our disc
geometry under certain approximations in sec. 3.4 and derive an approximate
analytic expression for the lightcurves in sec. 3.5.

In sec. 4, we present the different initial conditions for our models, which
include different internal properties and spatial locations for the
explosions, different vertical profiles and SMBH masses, as well as discs with
reduced densities. We also describe the numerical techniques used in this
paper. We devote sec. 5 to describing the obtained lightcurves for the
aforementioned models, both analytically and numerically. We discuss which
models are in good agreement between the analytic and numerical models, and
which differ and why.

Finally, we discuss the observational aspects of our results in sec. 6. We
discuss which of the events may be observable and could be distinguished from
the underlying AGN luminosity. We then estimate the total rates of these
events, their observed fraction and the origins of their large uncertainty. We
also briefly discuss the response of the disc and comment on AGN variability
and potential past and future observations of SNe explosions in AGNs.

We make a summary of our key points in sec. 7. Readers interested in the key
takeaway messages without the technical details can directly skip to this
section.

##  2 AGN disc structure and properties

AGN discs are powered by accretion and are expected to be geometrically thin
but optically thick (although there may be optically thin windows, see sec.
2.3). In order to describe AGN disc structures, we use several parameters. The
radial coordinate $r$ is expressed in units of the Schwarzschild radius
$r_{s}\equiv 2GM_{\bullet}/c^{2}\approx 1.974M_{8}\ \rm AU$, where
$M_{8}\equiv M_{\bullet}/10^{8}\textnormal{M}_{\odot}$ is the SMBH mass
normalized to $10^{8}\,\textnormal{M}_{\odot}$. The vertical structure is
governed by the scale height $H$ which is generally small compared to the size
of the disc ($H\ll r$).

###  2.1 Radial structure models

AGN discs are generally modelled as viscous accretion discs. On the largest,
galactic scales, discs are expected to be cold, gravitationally unstable and
fragment into stars, with the effective Toomre $Q\equiv c_{s}\Omega/\pi
G\Sigma\ll 1$, where $c_{s}$ is the sound speed, $\Omega$ is the orbital
frequency and $\Sigma$ is the surface density. However, only a fraction of the
gas is fragmented and the remaining gaseous component is still retained. The
transition from the unstable galactic discs to accretion discs around a
central object is still unclear. Angular momentum transport via global
torques, rather than local viscosity, may keep the disc marginally stable and
avoid fragmentation, which is the mechanism that is assumed in the Thompson et
al. (2005) model throughout the disc, even at smaller radii.

The structure of AGN discs can be inferred from the optical/UV spectral energy
distribution (SED) in the inner regions, and from mass inflow and star
formation in the outer regions. Sirko & Goodman (2003, hereafter SG03)
obtained an accretion disc model that fits the SED observations, while
Thompson et al. (2005) constructed a model regulated by mass inflow and star
formation. The SG03 model best describes the disc structure in the inner
regions, up to $\lesssim 10^{5}r_{s}$, which is the radial region of interest
in our study

The radial gas density structure of the SG03 model may be inferred from the
marginal gravitational instability condition (Eq. 15 of SG03),

$\rho(r)=\frac{\Omega^{2}}{2\pi GQ_{\rm min}}\approx 1.22\cdot
10^{-9}M_{8}^{-2}\left(\frac{r}{10^{3}r_{s}}\right)^{-3}{\rm g\ cm^{-3}}.$ (2)

The latter equation is valid for $r\geq 10^{3}r_{s}$, where $Q_{\rm
min}\approx 1$ in the SG03 model. The disc is expected to be dominated by
radiation pressure $p_{\rm rad}=aT^{4}/3$, where $a=4\sigma_{\rm SB}/c$ is the
radiation constant and $\sigma_{\rm SB}$ is the Stefan-Boltzmann constant. In
the SG03 model, radiation pressure indeed dominates in most parts of the disc.
However, the gas pressure may have a significant contribution to the total
pressure near $r\approx 10^{3}r_{s}$. Nevertheless, we use the radiation
pressure throughout the paper, and the adjustments needed to include other
sources of pressure are straightforward.

The aspect ratio of a thin disc should be small, i.e. $H/r\ll 1$. In the SG03
model, the aspect ratio is approximately given by

$\frac{H}{r}=8\cdot 10^{-3}\left(\frac{r}{10^{3}r_{s}}\right)^{1/2},$ (3)

in our range of interest.

###  2.2 Vertical structure models

The vertical structure of an AGN disc and its stability are also under debate
(see the review of Davis & Tchekhovskoy, 2020 for discussion and further
references). While gas pressure dominated discs are expected to be vertically
stable, the stability of radiation dominated discs is undecided. Theory
predicts that radiation dominated discs will be unstable (Shakura & Sunyaev,
1976). Contrary to this prediction, early radiative magneto-hydrodynamical
(MHD) shearing box simulations, such as those with the ZEUS code, found that
the disc is stable over many thermal timescales (Hirose et al., 2009).
However, these results were later challenged by simulations with the ATHENA
code. Although small box sizes reproduced the stability of Hirose et al.
(2009), larger box size simulations resulted in an eventual runaway (where the
disc either expands or collapses), although on much longer timescales than the
thermal timescale (Jiang et al., 2013). Iron bound-bound opacity may revert
the situation again and make the vertical structure stable around temperatures
of $1.8\cdot 10^{5}{\rm K}$ (Jiang et al., 2016).

Given the uncertainties involved in the stability and the detailed vertical
structure, we will remain agnostic about the true structure. Nevertheless, a
different vertical structure can lead to different observational signatures
for the same initial conditions of the explosion. Here we consider several
possible models for the vertical structure as described below, where their
main common feature is that they all describe a thin disc. Thus the vertical
density profile should fall off sharply where the vertical height $h$ is
larger than the scale-height, $h\gtrsim H$. Note that it is rather different
from a $\rho\propto r^{-2}$ distribution that is often assumed in studies of
interaction-powered SNe. For our canonical model, we assume a Gaussian density
profile

$\rho_{\rm
gas}(h)=\rho_{0}\exp\left(-\frac{h^{2}}{2H^{2}}\right);\quad\frac{z}{h}\in\mathbb{R},$
(4)

which is applicable for gas-dominated discs and is stable. We use hydrostatic
equilibrium to derive the density profile for radiation dominated discs (if
they are stable) in Appendix A,

$\rho_{\rm
rad}(h)=\rho_{0}\left(1-\frac{h^{2}}{6H^{2}}\right)^{3};\quad\left|\frac{h}{H}\right|<\sqrt{6},$
(5)

and present some of the results with this structure, as well as a simplified
step profile of constant $\rho_{\rm step}=\rho_{0}$ up to a scale height $H$
away from the midplane. We assume that the disc is stable and static with the
respective vertical structure prior to the explosion. We compare the results
for different vertical structures in sec. 5.3.

###  2.3 Opacity

At high temperatures, $T>10^{4}\rm K$, the gas is mostly ionized and the
opacity is dominated by electron scattering and bremsstrahlung (free-free). A
typical value of the opacity in this range can be the electron scattering
opacity $\kappa_{\rm es}=0.2(1+X)\ \rm cm^{2}\ g^{-1}$ where $X$ is the
Hydrogen fraction. The opacities may be enhanced by $\approx 2$ orders of
magnitude when detailed opacity tables are considered (Iglesias & Rogers,
1996, see also Fig. 3 of Davis & Tchekhovskoy, 2020). At temperatures below
$\sim 100\ \rm K$, dust absorption is the dominant effect that increases
opacity, which scales as $\kappa_{\rm dust}\approx 2.4(T/100{\rm K})^{2}\ \rm
cm^{2}\ g^{-1}$, while the opacity is roughly constant for $T\approx
10^{2}-10^{3}\ \rm K$ in the range of $\kappa\approx 1-10\ \rm cm^{2}\ g^{-1}$
(see Fig 1. of Thompson et al., 2005, which is based on more recent opacity
models of Semenov et al., 2003). For intermediate temperatures around
$10^{3}-10^{4}\rm K$, the gas is neutral, and most of the dust is already
sublimated, thus the opacity is extremely low. Depending on the radial
temperature profile, an optically thin window is expected to occur in AGN
discs. The discs are optically thin at $r\gtrsim 10^{4}r_{s}$ in the models we
adopted. Explosions in the optically thin region will look like standard SNe
as the emitted radiation will be unaffected by the surrounding disc, while the
associated kinetic outflows will be choked by the enormous mass. Explosions in
the optically thick and massive regions ($\approx 10^{4}r_{s}$ for
$M_{\bullet}=10^{7}\textnormal{M}_{\odot}$) will be completely choked (see
sections 4 and 6). For our analytical modeling, we use a constant electron
scattering opacity of solar composition ($X=0.7$), $\kappa_{\rm es}=0.34\ \rm
cm^{2}\ g^{-1}$. We use realistic OPAL tables (Iglesias & Rogers, 1996) for
our radiative transfer simulations (see sec. 4 for details).

###  2.4 Optical depth

For a constant opacity $\kappa$, the optical depth at height $h$, in the
direction normal to the disc, is given by

$\tau(h)=\kappa\intop_{h}^{\infty}\rho(h^{\prime})dh^{\prime}\equiv\tau_{0}\mathcal{F}(\zeta),$
(6)

where we defined for convenience the dimensionless quantities
$\tau_{0}\equiv\kappa\rho_{0}H$, and $\zeta\equiv h/H$. The functional form of
$\mathcal{F}(\zeta)$ depends on the vertical structure.

For a uniform step density, the integral is straightforward, and the solution
is $\tau(h)=\kappa\rho_{0}(H-h)=\tau_{0}(1-\zeta)$ for $0\leq h\leq H$ and
$\tau=0$ for $h>H$ since there is no material there. This leads to
$\mathcal{F}_{\rm step}(\zeta)=1-\zeta$.

For the Gaussian profile, applicable for the gas dominated disc, we have

$\mathcal{F}_{\rm
gas}(\zeta)=\frac{1}{H}\intop_{h}^{\infty}\exp\left(\frac{-h^{\prime
2}}{2H^{2}}\right)dh^{\prime}=\sqrt{\frac{\pi}{2}}{\rm
erfc}\left(\frac{\zeta}{\sqrt{2}}\right)$ (7)

where ${\rm erfc}(x)$ is the complimentary error function.

For the radiation-dominated profile, we have

$\mathcal{F}_{\rm
rad}(\zeta)=\frac{1}{H}\intop_{h}^{\sqrt{6}H}\left(1-\frac{h^{\prime
2}}{6H^{2}}\right)^{3}dh^{\prime}=\intop_{\zeta}^{\sqrt{6}}\left(1-\frac{\zeta^{\prime
2}}{6}\right)^{3}d\zeta^{\prime}=p(\zeta),$ (8)

where $p(\zeta)=16\sqrt{6}/35-\zeta+\zeta^{3}/6-\zeta^{5}/60+\zeta^{7}/1512$.

##  3 Analytic lightcurve modeling

Here we model the lightcurve of typical explosions. We first discuss the
assumptions of the model. We then proceed to discuss the velocity of the shock
front, and finally calculate the the lightcurve by drawing analogies between
our disc explosions and the recently developed theory of surface explosions.
Although our analysis is generic, when applicable, we refer to our canonical
model (1) of an explosion of energy $E_{0}=10^{51}\ \rm erg$ and ejecta mass
$M_{\rm ej}=1.3\textnormal{M}_{\odot}$, which serves as a proxy for a standard
type Ia SN explosion. The radial location of the explosion is at
$r=10^{3}r_{s}$ around an SMBH of mass $10^{7}\textnormal{M}_{\odot}$, in the
disc midplane. The disc has density $\rho_{0}=1.22\cdot 10^{-7}\ \rm g\
cm^{-3}$ and scale height $H=8\cdot 10^{-3}r=2.36\cdot 10^{13}\ \rm cm$ with a
gas-dominated Gaussian vertical profile. For the opacity, we make use of a
constant opacity $\kappa=0.34\ \rm cm^{2}\ g^{-1}$, which corresponds to
electron scattering at Solar compositions. The optical depth is then
$\tau=\sqrt{\pi/2}\kappa\rho_{0}H=1.2\cdot 10^{6}$. Other models are described
in sec. 4 and their initial conditions are summarized in table 1.

###  3.1 Assumptions

Here we describe the assumptions made in constructing our solution for shock
propagation, breakout and the resulting lightcurve. Consider an explosion
going off in dense material. If the material is sufficiently optically thick
and the prompt energy deposition is sufficiently energetic, the primary mode
of energy transport, at least initially, will be kinetic outflows rather than
radiation transport (photon diffusion). The initial disc temperature is
$T_{0}$, and the explosion is spherical. For collimated outflows, not
considered here, one would need to introduce an additional solid angle
$\Omega_{s}$ spanned by the outflow.

i) Is the outflow relativistic? The characteristic velocity of a non-
relativistic outflow is $v_{0}=\sqrt{E_{0}/M_{\rm ej}}$, so we require that
$(v_{0}/c)^{2}=E_{0}/(M_{\rm ej}c^{2})\ll 1$ for an outflow to be non-
relativistic. And indeed, for the canonical model, representative of type Ia
SNe, $v_{0}/c=0.02$. Similarly, core-collapse SNe also produce non-
relativistic outflows.

ii) Is the kinetic outflow the dominant form of energy transport? For this,
the velocity must be faster than the photon diffusion speed, at least
initially, or, equivalently, the depth must be sufficiently large. The optical
depth at height $z$ is formally given by Eq. 6. If at the explosion site the
local density is $\rho_{0}$ and the shock propagates over a typical distance,
$R_{0}$, the optical depth is approximately estimated as
$\tau\sim\kappa\rho_{0}R_{0}$, where $\kappa$ is the opacity, and the
diffusion speed is $c/\tau$. For kinetic outflow to dominate, we require
$v_{0}\gg c/\tau$, or $\tau^{2}E_{0}/M_{\rm ej}c^{2}\gg 1$. Combining
conditions i) and ii), the optical depth must also be large, $\tau\gg 1$. Note
that this defines a hierarchy of velocities $c/\tau\ll v_{0}\ll c$. For our
canonical model, the explosion is at the midplane, and the midplane optical
depth ($z=0$) is estimated as $\tau\gtrsim 10^{6}$, so these conditions are
initially satisfied.

iii) Does the shock leave radiation-dominated material behind? To answer this
question, we require that the radiation energy density $aT^{4}$, which may be
approximated by $\sim\rho v_{0}^{2}$ near shock, will be much larger than the
gas energy density $\rho k_{\rm B}T/m_{\rm p}$, where $k_{\rm B}$ is
Boltzmann’s constant and $m_{\rm p}$ is the proton mass. Eliminating the
temperature, this condition becomes equivalent to $(\rho
v_{0}^{2}/a)\gg(k_{\rm B}\rho/m_{\rm p}a)^{4/3}$, or

$v_{0}^{2}=\frac{E_{0}}{M_{\rm ej}}\gg\left(\frac{k_{\rm B}^{4}\rho}{m_{\rm
p}^{4}a}\right)^{1/3}.$ (9)

The latter in only an estimate within an order of magnitude.

Note that the latter three conditions are equivalent to the assumptions made
by Yalinewich & Matzner (2019, hereafter: YM19) in the context of surface
explosions. We return to the detailed modeling of YM19 in our lightcurve
modelling is sections 3.4-3.5.

###  3.2 Shock velocity profile

Here we discuss several models for the shock velocity profile. The description
is mostly generic, although sometimes, for a concrete example, we use the
canonical model (model (1) in table 1). We later investigate how the variation
of some of the disc and explosion properties affect our results in sec. 5.

For an explosion with energy $E_{0}$ and ejecta mass $M_{\rm ej}$, the typical
velocity is $v_{0}\equiv\sqrt{E_{0}/M_{\rm ej}}$. The outflow in the plane of
the disc will be quenched by the disc material, and the shock will break out
in the vertical direction. Therefore, we focus on the vertical height $z$
above the explosion. Generally speaking, a shock wave is always accompanied by
three waves: forward shock, reverse shock (or rarefaction), and the contact
discontinuity. In our context, the forward shock compresses the CSM, while the
reverse shock compresses the ejecta. The contact discontinuity separates the
ejecta and the CSM matter. For most of our models, the CSM matter is more
massive than the ejecta, so the reverse shock quickly traverses the ejecta
while the forward shock is still propagating outwards. We therefore focus only
on the forward shock, unless stated otherwise explicitly.

The shock velocity will depend on the shock front radius $z(t)$ and on the
density at this radius. Generally, three qualitatively different regimes will
be evident:

i) Free expansion regime: In this case, the swept mass $M(z)$ is much smaller
than the ejected mass, i.e., $M(z)\sim\rho(z_{\rm min})(z(t)-z_{\rm
min})^{3}\ll M_{\rm ej}$, where $z_{\rm min}$ is the vertical location of the
explosion site, $z(t)$ is the height reached by the shock at time $t$, and the
shock velocity is close to $v_{0}$.  ii) Sedov-Taylor deceleration regime:
Once the swept mass becomes comparable to the ejecta mass, $M(z)\gtrsim M_{\rm
ej}$ the shock front will decelerate. Dimensional analysis leads us to a
length scale $z(t)=\beta(E_{0}/\rho(z_{\rm min}))^{1/5}t^{2/5}$ and a velocity
scale

$v_{\rm ST}(t)=\frac{dz}{dt}=\frac{2\beta}{5}\left(\frac{E_{0}}{\rho(z_{\rm
min})}\right)^{1/5}t^{-3/5}.$ (10)

Here, $\beta$ is an order unity parameter which can in principle be estimated
from the energy equation. We take for simplicity $\beta=1$. The latter scaling
is the self-similar Sedov-Taylor (ST) solution (Sedov, 1946; Taylor, 1950),
where the shock front expands as $R\propto t^{2/5}$ and the shock velocity
decelerates $v_{\rm ST}\propto t^{-3/5}$.

iii) Sakurai accelerating regime: The ST picture is correct if the density is
uniform. The vertical structure of the disc is generally far from uniform.
Close to the disc edge (either the physical edge or the photosphere, if the
disc is formally infinite), the density gradient is quite steep, and the shock
propagates more easily, which causes it to accelerate with decreasing density
as $v\propto\rho^{-\mu}$, which is known as the Sakurai law (Sakurai, 1960).
For radiation dominated material, $\mu\approx 0.19$ fits well with previous
estimates (see also Matzner & McKee, 1999, hereafter MM99, their sec. 4.2 and
references therein for more details and discussion). Note that the latter
picture is valid if the swept mass $M(z)\gtrsim M_{\rm ej}$. Some of our
models have the swept mass comparable to or smaller than $M_{\rm ej}$, and in
these cases the ST phase is skipped and the shock may only accelerate or
directly break out.

Close to the centre, the density is almost constant, and the evolution is as
in the ST regime. Close to the edge, a small change in the radial scale leads
to a large change in the density, and therefore the shock will follow the
Sakurai law. It is possible to write down a general formula that includes both
cases, as done in MM99:

$v(t)=v_{0}\left(\frac{t_{1}}{t}\right)^{3/5}\left(\frac{\rho(z(t))}{\rho(z_{\rm
min})}\right)^{-\mu},$ (11)

where $t_{1}$ is the time when a transition from the free expansion to the ST
deceleration occurs and is given by $v_{\rm ST}(t_{1})=v_{0}$ from Eq. 10.
Plugging $t_{1}$ into the ST length scale leads to $z(t_{1})=z_{1}=(4M_{\rm
ej}/25\rho(z_{\rm min}))^{1/3}$. The transition occurs when the swept mass is
$4M_{\rm ej}/25=0.16M_{\rm ej}$. For the ejecta mass of
$1.3\textnormal{M}_{\odot}$, as in our model 1, the transition occurs at
$0.21\textnormal{M}_{\odot}$, or at radius $z\approx 1.5\cdot 10^{13}\
\mathrm{cm}=0.63H$.

The transition from the ST decelerating velocity to Sakurai accelerating
velocity occurs when $v(t)$ is at a minimum. From the ST solution, we can
invert $z(t)$ to $t(z)\propto z^{5/2}$ and find the location of the minimum
$z_{2}$ given by $d\ln\rho/d\ln z|_{z_{2}}=-3/(2\mu)$. For a Gaussian profile,
the velocity behaves as $v\propto z^{-3/2}\exp(\mu z^{2}/(2H^{2}))$, the
logarithmic derivative of the density is $d\ln\rho/d\ln z=-(z/H)^{2}$, and the
minimum is achieved at $z_{2}/H=(3/(2\mu))^{1/2}\approx 2.8$.

The three regimes and their transitions described above allow us to construct
two possible solutions. The first is the piecewise velocity solution. Here we
divide the behaviour into three regimes, where we first begin with free
expansion and transform to ST deceleration. The ST approach is not accurate
since it was developed for a uniform density. The density change is negligible
for $z\ll H$, but it is already significant at $z_{1}$. We use an effective
scaling where we replace $\rho(z_{\rm min})$ in the ST solution by $\rho(z)$.
Although not entirely self-consistent, we can write $v$ as in Eq. (11) but
with an effective density power law $\mu^{\prime}=\mu+1/5$, where the $1/5$ is
coming from the ST scaling of $z\propto\rho^{-1/5}$. The piecewise velocity is

$\displaystyle{v_{\rm p}(z)=v_{0}\begin{cases}1&z<z_{1}\\\
\left(\frac{z_{1}}{z}\right)^{3/2}\left[\frac{\rho(z)}{\rho_{1}}\right]^{-(\mu+1/5)}&z_{1}<z<z_{2}\\\
\left(\frac{z_{1}}{z}\right)^{3/2}\left(\frac{\rho_{2}}{\rho_{1}}\right)^{-(\mu+1/5)}\left[\frac{\rho(z)}{\rho_{2}}\right]^{-\mu}&z>z_{2}\end{cases}}$
(12)

where $\rho_{i}=\rho(z_{i})$ for $i=1,2$.

Another solution is the direct MM99 shock velocity

$v_{\rm MM99}(z)=\left(\frac{E_{0}}{M_{{\rm
ej}}+M(z)}\right)^{1/2}\left(\frac{\rho(z)}{\rho(z_{\rm min})}\right)^{-\mu},$
(13)

Note that Eq. (13) captures well all of the regimes. For $M_{\rm ej}\gg M(z)$
the density and hence the velocity are essentially constant. When $M(z)\gg
M_{\rm ej}$ but $\rho(z)\sim\rho(z_{\rm min})$, the velocity is $\propto
z^{-3/2}\propto t^{-3/5}$ as expected from a ST solution. Finally, once
$\rho(z)\ll\rho(z_{\rm min})$, the Sakurai term dominates the acceleration.

Figure 1: Shock velocity prescriptions for model 1 in the AGN disc, in units
of the speed of light $c$. The solid lines are the theoretical shock
velocities; the piecewise model (blue), Eq. (12), and the MM99 model (green),
Eq. (13). The dashed lines are the numerical shock velocities. The
hydrodynamical code HORMONE is used to obtain the matter velocity $v_{m}$ just
behind the shocked CSM, and the shock velocity is $v_{m}(\gamma+1)/2$ in the
disc geometry (dashed black) and spherical geometry (dashed red). The result
from the spherically-symmetric radiative transfer code SNEC is also presented
in magenta. For completeness, the photon diffusion velocity, which is also the
inverse of the optical depth, is shown in dot-dashed grey lines. The shock
breakout occurs when the gray line crosses the shock velocity, which occurs
close to $z\sim 4H$ for all shock velocity estimates for model 1.

Figure 1 shows the different velocity models for the canonical (model 1)
explosion together with the shock velocity obtained numerically. The detailed
description of the hydrodynamical code HORMONE and the radiative transfer code
SNEC appears in sec. 4. The three velocity regimes are clearly evident. The
flat free expansion phase lasts until $z_{1}/H\approx 0.63$ for model 1, where
the velocity starts to decelerate. The transition to Sakurai acceleration
occurs somewhere around $z/H\sim 2-2.8$. We see that the matter-derived shock
velocity $v_{m}(\gamma+1)/2$ better agrees with the piecewise model, while the
SNEC shock velocity is somewhere in between the piecewise and MM99 model. In
either case, the overall acceleration in the low-density outskirts of the disc
is underestimated.

In summary, there are two main theoretical shock velocity estimates for an
explosion inside the AGN disc: a piecewise fit explicitly taking account the
three distinct propagation regimes (Eq. 12), and a direct continuous measure
of the swept mass (MM99 \- Eq. 13). The numerical shock velocity falls between
the two theoretical shock velocity models.

###  3.3 Breakout time

Even though the shock is accelerating in the Sakurai regime, the photon
diffusion velocity rises more steeply. Eventually, the photon diffusion speed
$c/\tau(z)$ will intersect the shock velocity $v(z)$. The two velocities are
equal when the grey dot-dashed line intersects the velocity model of choice in
Fig. 1. The location of this intersection defines the breakout shell $z_{\rm
bo}$, which is implicitly given by $\tau(z_{\rm bo})=c/v(z_{\rm bo})$ and is
around $\sim 4H$ for model 1, regardless of the velocity model.

The breakout time $t_{\rm bo}$ is the time elapsed for the first light to be
beyond the breakout shell from the time of the explosion. In order to finally
escape, the photons need to diffuse up to the photosphere. This diffusion time
is much shorter than the breakout time.

The breakout time can be calculated by

$t_{\rm bo}=\intop_{z_{\rm min}}^{z_{\rm bo}}\frac{dz}{v(z)}$ (14)

The exact expressions are given either by special functions or integrated
numerically. We provide them in Appendix B, together with the expressions for
off-plane explosions.

###  3.4 Modelling AGN explosions as surface explosions

Surface explosions are explosions that occur in a relatively sparse material,
close to a surface of a medium and far away from the dense core of the medium
(e.g. a massive star). YM19 studied in detail such surface explosion in
massive stellar envelopes. Our explosion is analogous to a surface explosion
in the following way. Due to the planar geometry of the disc and its large
spatial extent $r\gg H$, an explosion in the midplane of the disc cannot
propagate in the in-plane direction, and the breakout occurs effectively on
the vertical surface of the disc. This is analogous to an explosion at a
distance $l\ll R_{\star}$ below the surface of a star, the model for which was
recently put forward by YM19. The shock propagates only upwards. Here we
parametrise our explosion as a surface explosion, with some limitations.

####  3.4.1 Choosing length scales

The YM19 model assumes that the density vanishes at some boundary $z_{b}$ and
that the density scales as a power law profile of the form
$\rho=\rho_{0}(x/l)^{\omega}$ where $x$ is the distance from the edge and $l$
is a typical length scale of order the distance of the explosion hotspot from
the edge. When applied to our extended model, $l$ is generally on the order of
H, but its exact definition has some ambiguities as we discuss below.

The surface explosion modelling in YM19 is encapsulated by two dimensionless
parameters, $\Gamma\equiv E_{0}/(c^{2}\rho_{0}l^{3)}$ and
$\tau_{0}=\kappa\rho_{0}l$. Initially, the only length scale we have is $H$.
In YM19, $l$ has several different roles:

i) The distance from the explosion to the edge.

ii) The transition between the end of the Sedov-Taylor phase and the beginning
of the shock acceleration phase.

iii) The value in the expression for the velocity at the end of the Sedov-
Taylor phase, $v_{{\rm ST}}=\sqrt{E_{0}/\rho_{0}l^{3}}$.

These three definitions appear to be contradictory. First, the free expansion
phase is ignored, or a point-source explosion is assumed. Second, after
passing a length $l$ (as in i.), formally, we should be at the edge, which is
definitely after the breakout. However, according to ii. and iii., it is only
the phase where the shock starts to accelerate and is assumed to be optically
thick.

These differences are small in the stellar atmosphere with its steep density
profile towards the edge. However, we saw that for model 1, the differences
between the transition to accelerating shock at $z_{2}=2.81H$, the breakout
shell at $z_{{\rm bo}}=3.99H$ for the piecewise model and $z_{{\rm bo}}=3.96H$
for the MM99 model, and the photosphere at $z_{{\rm ph}}=4.93H$, are around a
scale height away from each other. Moreover, free expansion phase ends at
$z_{1}=0.64H$, which is also comparable to the scale height.

The velocity at the end of the ST phase is $v_{\rm ST}\sim c\Gamma^{1/2}$.
Hence the effective length scale $l_{\rm eff}$ in this case should replicate
the total range of the swept mass from the end of the free expansion phase
until the breakout shell. We therefore set $l_{\rm eff}=z_{\rm bo}-z_{1}$,
unless otherwise specified.

$\Gamma=\frac{E_{0}}{c^{2}(M_{{\rm ej}}+\rho_{0}l_{{\rm eff}}^{3})}.$ (15)

Here, the inclusion of $M_{\rm ej}$ in the denominator extends the analysis of
YM19 to include the free expansion phase and also avoids unphysically large
values of $\Gamma$ for low densities or small $l_{\rm eff}$.

The optical depth $\tau_{0}$ should be comparable with the midplane optical
depth, hence we keep $\tau_{0}$ as $\kappa\rho_{0}H$.

####  3.4.2 Effective power law index $\omega_{\rm eff}$

Once we set $l_{\rm eff}$ and $\Gamma$ for the explosion and also have
$\rho_{0},H$ and $\tau_{0}$ from the AGN disc model, we can assign an
effective power law index $\omega_{\rm eff}$ in the following way:
$\omega_{\rm eff}$ is a local value and formally varies along the vertical
shock propagation.

For a local density profile

$\rho(x)=\rho_{a}(x/l_{\rm eff})^{\omega_{\rm eff}},$ (16)

where $\rho_{a}$ is the density at a distance $l_{\rm eff}$ from the edge, and
$0\leq x\leq l_{\rm eff}$ is measured from the edge. We impose a matching
condition on $d\ln\rho/dz$ between the assumed disc density profile and the
power law model in order to fix $\omega$:

$\rho_{a}=\rho(z_{{\rm ph}}-l_{\rm eff});\quad\omega=-l_{\rm
eff}\left.\frac{d\ln\rho}{dz}\right|_{z_{\rm ph}-l_{\rm eff}}\ .$ (17)

For the slab profile, the derivative is zero and thus $\omega_{\rm eff}=0$
everywhere. For the Gaussian profile, $d\ln\rho/dz=-z/H^{2}$, and

$\omega_{{\rm gas}}(l_{{\rm eff}})=\frac{l_{{\rm eff}}(z_{\rm ph}-l_{{\rm
eff}})}{H^{2}}.$ (18)

For the radiation dominated profile the power law is

$\omega_{{\rm rad}}(l_{{\rm eff}})=6\frac{\sqrt{6}H-l_{{\rm
eff}}}{2\sqrt{6}H-l_{{\rm eff}}},$ (19)

where we used the actual edge of the disc $z_{\rm edge}=\sqrt{6}H$ instead of
the photosphere.

We note that $\rho_{a}$, $l_{\rm eff}$ and $\Gamma$ are determined by the
extent of the mass swept by the outflow prior to the breakout. $\tau_{0}$ is
also well defined. Although Eq. 17 suggests that determining $l_{\rm eff}$
uniquely determines $\omega_{\rm eff}$, this is misleading, since in the YM19
model of Eq. (16), the density profile is a global power with index
$\omega_{\rm eff}$, while in our case the power law is only a local fit to the
density profile. Therefore, we have the freedom to decide where to match the
density profile to the power law. We opted to do so at the depth of $l_{\rm
eff}$ beneath the breakout shell in equations 18 and 19. We could
alternatively choose to do so at the breakout radius $z_{\mathrm{bo}}$ in an
attempt to model the immediate post-breakout behavior more accurately; we
consider this alternative choice in section 5.4. Regardless of the local
choice of $\omega_{\rm eff}$, the dependence of the luminosity on $\omega$ is
expected to be much weaker than on the other parameters.

###  3.5 Lightcurves

As mentioned in 1.2, the slab geometry leads to adiabatic expansion in the
optically thick regime and rapid decay. The setting is similar to an explosion
close to a stellar surface studied by YM19. They applied the latter arguments
for a surface explosion using the dimensionless parameters $\Gamma$ and
$\tau_{0}$ as discussed in sec. 3.4, where an additional phase, where the
hotspot of the explosion becomes optically thin enough, which allows depletion
of material deeper than the hotspot and forming a crater. The time of the
transition where the hotspot is exposed is found to be $t_{\rm
sph}=(H/c)\tau_{0}^{1/2}\Gamma^{-1/4}$ (Eq. 24 of YM19).

Once we find $l_{\rm eff}$, $\omega_{\rm eff}$, followed by $\Gamma$ and
$\tau_{0}$ as describe in sec. 3.4, we can use YM19’s modelling for the
lightcurve. From now until the end of the section, we drop the ’eff’ subscript
and simply write $\omega$ to avoid cumbersome notation. The resulting
luminosity is

$\frac{L(t)}{E_{0}c/H}=\begin{cases}\Gamma^{(\omega\mu-2\omega/3-5/6)\delta_{-}}\tau_{0}^{(5\omega\mu/3-\omega-4/3)\delta_{-}}\left(\frac{ct}{H}\right)^{-4/3}&\rm{pl}\\\
\Gamma^{(-\omega\mu+1/6)\delta_{+}}\tau_{0}^{(\omega\mu/-\omega-4/3)\delta_{+}}\left(\frac{ct}{H}\right)^{(-4\omega\mu+2/3)\delta_{+}}&\rm{sph}\end{cases}$
(20)

in the planar and spherical phases, respectively, where
$\delta_{\pm}\equiv(1\pm\omega\mu+\omega)^{-1}$. For concreteness, we choose
$\omega(l_{\rm eff})$ as in equations 18 and 19.

In certain cases, we will use an expansion of $\omega_{\rm eff}$ close to the
edge, where we use $\omega^{\rm edge}_{\rm eff}=\omega_{\rm eff}(l_{\rm
edge})$ with $l_{\rm edge}=z_{\rm edge}-z_{\rm bo}$, where the edge is either
the physical edge of the disc (for radiative and slab models), or the
photosphere (for Gaussian models). We discuss this choice and its implications
in the results section.

Figure 2: Analytical lightcurve for model 1. We show the analytic bolometric
luminosity (i.e. the lightcurve) as a function of time with the solid blue
curve. The numerical SNEC result is shown as a dashed red curve. The analytic
curve is the maximum of the planar and spherical phases in Eq. 20. The
vertical dot-dashed lines are the breakout times for the two velocity models
given in Eq. 33 and 34, where the rightmost one is for the piecewise model
since it has a lower velocity than the MM99 model. The gray area between them
is the expected range of the breakout times. The dotted blue vertical line is
the end of the spherical expansion phase, $t_{\rm sph}$.

Fig. 2 shows the analytic lightcurve for model 1 in solid blue and the
numerical SNEC lightcurve in dashed red. We shifted the analytical curve in
time so that the breakout falls on top of the maximal luminosity of the
numerical curve. After the peak, SNEC model luminosity is initially more
shallow than the $L\propto t^{-4/3}$ decay expected in the planar phase from
the simple arguments of section 1.2, but subsequently steepens before becoming
more shallow again at the transition to the spherical phase. After the initial
breakout with $L\sim 10^{44}\ \rm erg\ s^{-1}$, the luminosity rapidly decays
by one order of magnitude over $\sim 0.4\ \rm d$ from breakout and then decays
more shallowly during the spherical phase. The end of the spherical phase is
indicated by the dotted vertical line around $\sim 150\ \rm d$. The gray area
in the region between $t_{\rm bo}|_{\rm MM99}<t<t_{\rm bo}|_{\rm pw}$ is given
by Eqs. 34 and 33, respectively, and indicates the range of times when the
breakout occurs. In addition, the SNEC lightcurve is also plotted in dashed
red, where the peak luminosity of the analytic model is set to match the peak
luminosity of the SNEC lightcurve, which falls in the margins of the breakout
time. Overall, there is a good correspondence, and the bumps in the SNEC
lightcurve can be attributed to the reverse shock initially and varying
opacity effects at later times, as discussed in detail in section 5.

We note that in order to set $z_{\rm bo}$ (and then $l_{\rm eff}$, $\omega$,
$\Gamma$ and the lightcurve), we need to solve for the breakout shell
$\tau(z_{\rm bo})=c/v_{\rm bo}$. This is done numerically using the standard
fsolve function of Python’s scipy.optimize package. We also need to choose a
model for the velocity. We use the earlier derived piecewise model, Eq. 12,
since the numerical breakout time is close to the breakout time evaluated
using this velocity prescription. We also tried the MM99 prescription, and it
gave almost indistinguishable results in most cases.

##  4 Numerical Simulations

Here we describe the numerical approach and initial conditions for our
simulations and analysis.

###  4.1 Initial conditions

We already introduced the initial conditions for our baseline model at the
beginning of section 3. In order to explore the sensitivity to other
parameters of the disc and the explosion, we also studied several other
models, varying one or two parameters at a time. For changes in the global
disc structure with the SMBH mass, we adopt the following scalings from the
SG03 model. At radial locations scaled by the gravitational radius of the
SMBH, i.e. at $r\propto r_{s}\propto M_{\bullet}$, the midplane density is
$\rho\propto\Omega^{2}\propto M_{\bullet}/r^{3}\propto M_{\bullet}^{-2}$. The
rescaled interacting mass is $\rho H^{3}\propto M_{\bullet}$. The optical
depth is also rescaled as $\tau\propto\rho H\propto M_{\bullet}^{-1}$. The
aspect ratio is taken from Eq. 3.

# | $M_{\bullet}$ | $r/r_{s}$ | $E_{51}$ | $\rho(z)$ | $z/H$ | $M_{\rm CSM}/\textnormal{M}_{\odot}$ | Remarks
---|---|---|---|---|---|---|---
1 | $10^{7}$ | $10^{3}$ | $1$ | Gas | $0$ | $12.8$ | a
2 | $10^{7}$ | $10^{4}$ | $1$ | Gas | $0$ | $390$ | b
3 | $10^{8}$ | $10^{3}$ | $1$ | Gas | $0$ | $128$ | c
4 | $10^{6}$ | $10^{4}$ | 1 | Gas | 0 | $39$ | b,c
5 | $10^{7}$ | $10^{3}$ | $10$ | Gas | $0$ | $12.8$ | d
6 | $10^{7}$ | $10^{3}$ | $0.1$ | Gas | $0$ | $12.8$ | d
7 | $10^{7}$ | $10^{3}$ | $1$ | Rad | $0$ | $7.6$ | e
8 | $10^{7}$ | $10^{3}$ | $1$ | Step | $0$ | $3.4$ | e
9 | $10^{7}$ | $10^{3}$ | $1$ | Gas | $1$ | $1.92$ | f
10 | $10^{7}$ | $10^{3}$ | $1$ | Gas | $2$ | $0.15$ | f
11 | $10^{7}$ | $10^{3}$ | $1$ | Rad | $1$ | $0.5$ | e,f
12 | $10^{7}$ | $10^{3}$ | $1$ | Rad | $2$ | $6.4\cdot 10^{-4}$ | e,f
13 | $10^{6}$ | $10^{3}$ | 1 | Gas | 0 | $1.28$ | c
14 | $10^{7}$ | $10^{3}$ | $1$ | Gas | $0$ | $1.28$ | g
15 | $10^{7}$ | $10^{3}$ | $1$ | Gas | $0$ | $12.8$ | h
Table 1: Initial conditions for the simulations. The columns are: simulation
number, central BH mass (in $M_{\odot}$), radial location, explosion energy
(in $10^{51}\rm\ erg$), the vertical height above the midplane, CSM mass (Eq.
21), and remarks describing the parameters being varied: a) The canonical
model. b) Model with varied radial distance $r$. c) Model with varied mass of
the SMBH. d) Model with varied explosion energy $E_{0}$. e) Model with varied
vertical density profile $\rho(z)$. f) Model with varied vertical distance
$z$. g) Model with a reduced density $\rho_{0}=1.23\cdot 10^{-8}\ \rm g\
cm^{-3}$ to mimic a starved AGN. h) Model with a larger ejecta mass, $M_{\rm
ej}=10\textnormal{M}_{\odot}$, to mimic a core-collapse SN.

Table 1 lists the grid of initial conditions for the numerical simulations.
Model 1 is our baseline canonical model, representing a typical type Ia SN
explosion. Models 2-4 and 13 explore the impact of the SMBH masses and the
radial location. Models 5-6 explore the impact of the explosion energy. Models
7-8 and 11-12 explore the impact of the vertical profile. Models 9-12 explore
the impact of the vertical location. Model 14 mimics a starved AGN disc with
reduced density, while Model 15 mimics a core collapse SN by considering more
massive ejecta than in other models, $M_{\rm ej}=10\ M_{\odot}$. The explosion
energy of SN 1987A is estimated at $(1.5\pm 0.12)\times 10^{51}$ erg (Utrobin,
2007). The observed explosion energies of core-collapse SNe range from a few
times $10^{50}$ erg to $\sim 10^{53}$ erg for hypernovae with massive
progenitors (Janka, 2012). We therefore keep the explosion energy at $10^{51}$
erg for model 15 as most representative for core-collapse SNe.

The CSM mass in the penultimate column of table 1 is defined as

$M_{\rm CSM}=4\pi\intop_{z_{\rm min}}^{z_{\rm max}}\rho(z)(z-z_{\rm
min})^{2}dz,$ (21)

which depends on the density profile, and the maximal and minimal vertical
locations. In Eq. 21, we assumed a spherical approximation, namely that the
most significant interaction of the ejecta with the disc material will be in
the vertical direction, allowing us to approximate the density profile at
radius $r$ as $\rho(r)=\rho(z=r)$. This is true in our geometry since the in-
plane directions are choked. We perform the explicit calculation of the CSM
mass for each model in Appendix C.

###  4.2 Spherically symmetric lightcurve modelling

To model the lightcurves arising from the CSM interaction numerically and
verify our analytic models, we perform detailed simulations with the supernova
explosion code SNEC (Morozova et al., 2015a; Morozova et al., 2015b). SNEC is
a spherically-symmetric Lagrangian radiative-hydrodynamics code that accounts
for shock propagation through the use of artificial viscosity. Compared to our
analytic models, SNEC directly solves the equations of radiative hydrodynamics
and, therefore, allows us to assess the validity of our approximations.
Additionally, it allows us to gain some intuition into the importance of
ionisation, detailed opacities and heating from ${}^{56}{\rm Ni}$ for forming
the final lightcurves.

SNEC uses OPAL opacities (Iglesias & Rogers, 1996) in the high-temperature
regime (at $3.75<\log_{10}T/{\rm K}<8.7$) and the solar-scaled Ferguson et al.
(2005) opacities, which accounts for molecular contributions, in the low-
temperature regime (at $2.7<\log_{10}T/{\rm K}<4.5$). In the overlapping
regions, the low-temperature opacity is preferred. The code accounts for
ionisation states by solving the equilibrium Saha equation and models
radiative transfer through flux-limited equilibrium photon diffusion. SNEC
also models the contribution to the lightcurve from the radioactive decay of
${}^{56}{\rm Ni}$ and 56Co. The decay leads to the deposition of gamma-rays
and positrons, which heat the ejecta, with gamma-rays dominating the heat
deposition. The code solves for gamma-ray propagation and absorption through
solving for radiative transfer equation with grey opacities.

We set up explosions by initialising balls of material with an exponential
density profile with a radial scale height set so that the density at the
boundary of the ball is $1\,{\rm g}/{\rm cm}^{3}$. Such a setup allows for a
mild density contrast near the CSM boundary, which improves the resolution of
the shock. For models 1-14, having a type Ia-like engine, we set the ball
radius to $0.1\textnormal{R}_{\odot}$ and assign it
$1.3\textnormal{M}_{\odot}$ of C-O material. For model 15, we set a
$10\textnormal{M}_{\odot}$ ball of $2\textnormal{R}_{\odot}$ radius, thus
mimicking a core-collapse SN potentially stripped by a companion. The engine
acquires the explosion energy split equally into kinetic energy and thermal
energy. Furthermore, the kinetic energy is distributed to provide comparable
escape velocities for all the velocity bins of ejecta. We have verified that
details of the setup, such as the fraction of the explosion energy injected as
thermal energy or the initial radius of the explosion ball, have a negligible
effect on the lightcurves. We simulate all the models twice, with and without
56Ni. To compare our results to the analytic models of CSM interaction, we use
the simulations without 56Ni, and we examine the qualitative effects of 56Ni
in the results section 5. For models with 56Ni, we assume
$0.6\,\textnormal{M}_{\odot}$ of the radioactive material, typical for type Ia
explosions (apart from model 15, wherein we assumed
$0.2\,\textnormal{M}_{\odot}$ of 56Ni). The CSM is set up in concentric
shells, following the density and composition profiles for the different AGN
models, with the full models containing between $1000$ and $2000$ grid cells.
We sample the lightcurves produced in SNEC with a step of $5$ minutes to
resolve the sharp hour-long peaks in some of the models. Overall, the
spherically-symmetric models are initialised similarly to how it is done in
the analytic modelling.

The code allows us to trace the detailed properties of the CSM and the shock
over time and calculate the lightcurves for the models.

###  4.3 Effects of the geometry on hydrodynamics and lightcurves

While the SNEC code can deal with most of the relevant physics such as
radiative transfer, detailed opacity tables, nuclear energy deposition from
56Ni and ionisation, it cannot model non-spherical geometries. To understand
how the supernova ejecta interact with the non-spherical AGN disc material, we
carry out additional 2D axisymmetric hydrodynamic simulations of the
interaction. For this, we use the hydrodynamic code HORMONE, which is a grid-
based code that solves the hydrodynamic equations through a Godunov-type
scheme (Hirai et al., 2016). We employ an equation of state with contributions
from ideal gas and radiation. A spherical coordinate system is used, where we
assume axisymmetry with the symmetry axis taken perpendicular to the disc and
the coordinate origin placed at the centre of the supernova explosion. In our
axisymmetric treatment we ignore any motion of gas in the disc prior to the
explosion, such as shearing motion expected in Keplerian discs, and instead
treat the disc as an initially static slab of material. We also could not
apply the gravity from the central SMBH, so we ignore all gravitational
forces, including self-gravity. Therefore, we neglect the initial thermal
energy in the disc material too, since it will lead to artificial expansion
without the gravitational forces that keep it bound. These assumptions are
still valid as long as the time-scale for the shock to reach the surface of
the disc is much shorter than the local Keplerian period, and the shock energy
is larger than the pre-explosion thermal energy of the disc.

We set up the simulation by placing a slab of material to represent the AGN
disc. At the coordinate origin, we place a homologously expanding supernova
ejecta model with an exponential density profile as

$\rho_{\textnormal{ej}}(r)=\frac{3\sqrt{6}M_{\textnormal{ej}}}{4\pi
v_{0}^{3}t_{0}^{3}}\exp{\left(-\frac{\sqrt{6}r}{v_{0}t_{0}}\right)},$ (22)

where $t_{0}$ is the time since explosion. Such exponential profiles are
commonly seen in type Ia SN ejecta models (e.g. Nomoto et al., 1984). The
velocity profile is set to be $v(r)=r/t_{0}$. When integrated from 0 to
infinity, the total mass is $\int_{0}^{\infty}4\pi
r^{2}\rho_{\mathrm{ej}}dr=M_{\mathrm{ej}}$ and the total kinetic energy
becomes $\int_{0}^{\infty}\frac{1}{2}\rho_{\mathrm{ej}}(r)v(r)^{2}\cdot 4\pi
r^{2}dr=E_{0}$. We cut off the explosion profile at a radius
$r=R_{\textnormal{exp}}$ and the time since explosion is set through
$t_{0}=R_{\textnormal{exp}}/(4v_{0})$. The factor 4 is chosen so that the
total integrated ejecta mass and energy are still very close to
$M_{\mathrm{ej}}$ and $E_{0}$. We set a small enough value of
$R_{\textnormal{exp}}$ such that the mass of the disc material within that
radius is negligible compared to the ejecta mass
($\int_{0}^{R_{\textnormal{exp}}}4\pi r^{2}\rho_{\textnormal{disk}}dr\ll
M_{\textnormal{exp}}$), where we can safely assume that the ejecta have freely
expanded up to that radius. For all the models presented here, we use an
ejecta mass of $M_{\textnormal{ej}}=1.3~{}\textnormal{M}_{\odot}$, which
represents the ejecta mass for typical type Ia SNe, except for model 15, which
has a representative mass of $M_{\rm ej}=10\textnormal{M}_{\odot}$, typical
for core-collapse SNe.

The ejecta are covered with at least 30 grid points, and the radial grid size
is increased as a geometrical series as it goes out. The outer boundary is
taken at $\sim 100~{}H$, and we use $\sim 1000$ grid points in the radial
direction. We divide the polar direction into 400 cells spaced uniformly in
$\cos\theta$ so that the solid angle of each cell is equal and the disc is
properly resolved.

The lightcurves are estimated through a post-process ray-tracing procedure
similar to that of Suzuki et al. (2016). See Appendix D for more details of
the lightcurve calculation method.

##  5 Results

Here we present the analytic and numerical lightcurves for the different
models.

###  5.1 Peak luminosities

| $L_{44}^{{\rm YM19}}$ | $L_{44}^{{\rm SNEC}}$ | $C_{L}^{{\rm YM19}}$ | $C_{L}^{{\rm SNEC}}$
---|---|---|---|---
1 | 1.48 | 1.7 | 7.97 | 9.12
2 | 0.12 | 0.14 | 5.54 | <6.48
3 | 0.26 | 0.28 | 4.34 | 4.69
4 | 0.62 | 0.93 | 5.8 | 8.6
5 | 29.38 | 27.68 | 10 | 9.44
6 | 0.08 | 0.13 | 6.63 | 10.87
7 | $4.24^{a}$ | 4.42 | 6.91 | 7.21
8 | $24.70^{b}$ | 33.5 | 1.15 | 1.56
9 | $3.96^{a}$ | 7.86 | 10.65 | 30.38
10 | $12.77^{a}$ | 14.90 | 40.97 | 26.84
11 | $7.82^{a}$ | 21.93 | 10.68 | 26.82
12 | $22.55^{a}$ | 5.43 | 47.2 | 11.53
13 | 7.85 | 2.65 | 21.05 | 7.1
14 | 8.89 | 14.16 | 14.69 | 23.4
15 | 1.45 | 1.51 | 13.45 | 14.04
Table 2: Peak luminosities and their coefficients. The first column is the
model number. The luminosities are normalized to $L_{44}=L/(10^{44}$ erg/s).
The second column is $L^{\rm YM19}$, the peak luminosity derived from the
maximal value of Eq. 20. The third column is the maximum luminosity obtained
in SNEC simulation. The forth column is the ratio $L^{\rm YM19}/(\rho
z^{2}v^{3})_{\rm bo}$, which is an estimate for $C_{L}$ from Eq. 1. The
denominator is evaluated at the breakout shell. The last column is the same,
only for the numerical maximal SNEC luminosity. a) Evaluated $\omega$ at the
breakout radius for the analytical fit. b) Used $\omega=1$.

In table 2, we show the peak bolometric luminosities from the analytical YM19
modelling and from SNEC simulations. The last two columns also show the
dimensionless prefactor $C_{L}$ for which the estimate in Eq. 1 matches the
relevant luminosity. We see that the peak luminosities from the YM19 modelling
correspond well with the SNEC results for most models. The most discrepant
models are 11-13, for which the luminosities differ by factors of 3 to 4.

The estimate of $C_{L}$ expressed from Eq. 1 shows a discrepancy of up to two
orders of magnitude compared to $C_{L}$ calculated with the previous method.
There is no clear trend, and larger values are obtained for models of compact
structure and low density, namely models which explode off-plane (model 9-12)
or with a reduced SMBH mass and scale height (model 13) or starved AGN (model
14). Excluding them and the slab model (model 8), which has discontinuous
density and is challenging to simulate, leads to a factor of two spread in
$C_{L}$ values, $C_{L}\approx 7-14$ when compared against both YM19 models and
SNEC simulations (models 1,5-7,15). We note that for a spherical model such as
the SNEC one, the spherical geometry suggests the value $C_{L}=4\pi\approx
12.6$ for thermalized matter, which is within this range.

We devote the remainder of this section to exploring predicted model
lightcurves and their unique features, including comparing numerical and
analytical results.

###  5.2 The lightcurve shape and the effect of changing the energy and ejecta
mass

Fig. 3 shows the lightcurves of models 1, 5, 6, which correspond to different
explosion energies ($10^{51}$, $10^{52}$ and $10^{50}{\rm erg}$, respectively)
and also the lightcurve of model 15, which is identical to model 1, but with a
larger ejecta mass of $M_{\rm ej}=10\ M_{\odot}$, which represents a fiducial
model of core-collapse SN. We also ran a model with $10M_{\odot}$ ejecta and a
lower energy of $10^{50}$ erg. The resulting lightcurve peak luminosity is
very similar to the low energy of model 6, and is not shown here. Overall we
see a good correspondence between the analytical lightcurve (solid) and the
numerical results based on the SNEC code (dashed). We see that the larger the
energy is, the faster is the breakout time, and the brighter is the peak. The
dependence of the peak luminosity on explosion energy is slightly steeper than
linear, $L_{\rm peak}\propto E_{0}\Gamma^{(1/2+\mu\omega)\delta_{-}}\propto
E_{0}^{1+(1/2+\mu\omega)\delta_{-}}$, which gives the value of $1.28-1.29$ for
our range of density power laws $\omega_{i}$, where $i={1,5,6}$. $\Gamma$ also
depends on $l_{\rm eff}$, but $l_{\rm eff}$ varies only within $10\%$ between
the different models. The additional corrections from $l_{\rm eff}$ slightly
change the scaling to $1.23-1.24$.

Thus, the peak luminosity reaches $\sim 3\times 10^{45}\ \rm erg\ s^{-1}$ for
model 5, while it is slightly less than $10^{43}\rm\ erg\ s^{-1}$ for model 6.
The different $\omega$’s arise due to the different breakout shells in each
model: the less energetic models have a lower velocity by a factor of
$E_{0}^{1/2}$ and hence breakout at later times. The more massive ejecta model
15 has a smaller velocity but has a longer phase of free expansion. Thus its
effective power law is $\omega_{15}=6.04$, which is the largest one in our
models. This is due to the fact that in this model, $l_{\rm eff}=2.67H$ is
very close to the maximal value of $\omega_{\rm 15,max}=6.076$ obtained at
half the photosphere $z_{\rm ph}/2H=4.93/2=2.465$.

In the spherical phase, the dependence of the luminosity on $\Gamma$ is a very
shallow power law, while the dependence on time since the breakout is
$L\propto t^{-x}$ with $x=0.44-0.47$. Note that this is steeper than the
maximal slope of $-0.35$ derived by Nakar & Sari (2010) for radiation
dominated material, which is obtained for $\omega=3$. This is because our
extended disc model allows larger values of the effective power law index.

Figure 3: Probing different energies and ejecta masses. The solid lines are
the analytical lightcurves where the semi-transparent dashed lines are the
SNEC simulations. Model 1 is our canonical model (blue), model 5 has a $10$
times higher explosion energy of $E_{0}=10^{52}\ \rm erg$ (red), model 6 has a
$10$ times lower explosion energy of $E_{0}=10^{50}\ \rm erg$ (green). Model
15 (black) has the same explosion energy as model 1, but with larger ejecta
mass of $M_{\rm ej}=10\textnormal{M}_{\odot}$. The coloured transparent
regions indicate the range of possible breakout times as in Fig. 2, with each
transparent colour matching the correspondingly coloured model. The vertical
dotted lines correspond to the end of the spherical phase for each model shown
in their respective colours, also similar to Fig. 2. The effective power laws
for the density for each model are $\omega_{1}=5.3$ (the same as in Fig. 2),
$\omega_{5}=4.7$, and $\omega_{6}=5.75$, and $\omega_{15}=6.04$.

###  5.3 Changing the SMBH mass and midplane density

Figure 4: Probing different AGN disc locations and properties. Models 4 and 13
correspond to a $10^{6}\textnormal{M}_{\odot}$ SMBH, while model 14
corresponds to a $10^{7}\textnormal{M}_{\odot}$ SMBH. Furthermore, the
explosion in model 4 is placed at $10^{3}r_{s}$, while for models 13 and 14,
the explosion is placed at $10^{4}r_{s}$. The line style and shaded areas have
the same meaning as in Fig. 3. The effective power laws are $\omega_{4}=4.97$,
$\omega_{13}=6.7$, and $\omega_{14}=4.98$.

Fig. 4 shows the lightcurves of models 4, 13, 14 where the AGN disc structure
is altered. In models 4 and 13, the reduced SMBH mass of
$10^{6}\textnormal{M}_{\odot}$, compared to $10^{7}\textnormal{M}_{\odot}$ for
model 13, induces an increase in the midplane density, but the scale height is
smaller. Model 4 is exploding at a more distant radial location of
$10^{4}r_{s}$, compared to the other two models that explode at $10^{3}r_{s}$,
thus making the swept CSM mass larger. While the analytical and SNEC
lightcurves broadly agree for model 4, they differ more significantly for
models 13 and 14. The general trend is that the SNEC lightcurves decrease more
rapidly compared to the analytic ones. The reason may lie in the fact that
models 13 and 14 have a reduced CSM interaction mass, which is slightly less
than the ejecta mass. Radiative cooling may contribute to the extra drop in
the luminosity in the simulations, not accounted for in the analytical model,
as discussed in sec. 3.5. Moreover, the peak luminosities in model 13 also
differ by a factor of a few. The difference in the peak luminosities and the
wiggles seen in the SNEC lightcurves for model 13 may be related to the
reverse shock travelling back through the ejecta. Once the reverse shock
reaches the origin of the spherical model, it reflects and starts travelling
forward. The reverse shock arrives at the breakout region with a delay
compared to the initial forward shock and is typically weaker. The reason the
reverse shock is prominent in model 13 is likely because the CSM radius is
small compared to other models, and the CSM mass is small enough to allow fast
breakout.

###  5.4 Changing vertical profile

Fig. 5 shows the dependence of the lightcurve on the vertical structure. We
see that the radiation-dominated profile (model 7) and the uniform step
profile (model 8) lead to progressively earlier and brighter breakouts. This
is due to the fact that the respective edges of the different vertical
profiles (and also their breakout shells) are more compact than the standard
Gaussian profile in model 1.

In Fig. 5, for the calculation of $\omega$ we used $l_{\rm edge}=z_{\rm
edge}-z_{\rm bo}$ in Eq. 19, where $z_{\rm edge}$ is the physical edge of the
disc for model 7, and the photosphere $z_{\rm phot}$ for model 1. The spatial
extent of the explosion which enters into the calculation of $\Gamma$ is
unchanged, $l_{\rm eff}=z_{\rm bo}-z_{1}$. The resulting values for these
models are $\omega_{1}^{\rm edge}=3.76$ for model 1, and a higher value of
$\omega_{7}^{\rm edge}=2.9$ for model 7 (compared to $\omega_{7}=1.46$ in the
default choice). For model 7, since $z_{\rm bo}$ is closer to the disc edge,
$\omega_{7}^{\rm edge}$ is closer to its limiting value of $3$, which is also
the power-law suitable for spherical radiative atmospheres. For model 8, the
profile is uniform and discontinuous, hence $\omega_{8}=0$ throughout the disc
and infinite at the edge. Since $\omega_{8}$ is ill-defined for model 8, we
express the range of possible lightcurves by shading a grey area, where the
upper limit is constructed by choosing $\omega_{8}=100$, and the lower limit
is given by $\omega_{8}=0$. For comparison, we also plot the analytic fit with
$\omega_{8}=1$, which gives a roughly flat lightcurve, similar to the SNEC
lightcurve (though we caution that a numerical treatment will also struggle
with a discontinuous density profile, and the SNEC models are likely not
robust in this case).

For model 1, this fit is comparable to the one with the default choice (as in
Figs. 2 and 3). This is also a representative case of all the gas dominated
Gaussian profiles, namely that they are not very sensitive to $\omega$. For
radiation dominated profiles, however, the latter choice of $\omega_{7}^{\rm
edge}$ gives a much better match to the SNEC results. We speculate that the
reason may be due to the fact that radiative dominated and slab profiles are
more compact and have an actual edge.

Figure 5: Probing different density profiles. The line style and shaded areas
have the same meaning as in Fig. 3. For models 1 and 7, we chose a different
location for evaluating $\omega$, at the breakout shell, namely using $l_{\rm
edge}=z_{\rm edge}-z_{\rm bo}$, in Eq. 19, where $z_{\rm edge}$ is either the
physical edge of the disc (e.g. radiation-dominated profile as in model 7) or
the photosphere, if the disc is formally infinite (Gaussian profile as in
model 1). The spatial extent $l_{\rm eff}$ which is used to calculate $\Gamma$
is unchanged. Model 1 is our canonical model (blue), model 8 has a step
profile (green), model 7 has a radiation dominated profile (red). The power-
law index for model 1 is now $\omega_{1}^{\rm edge}=3.76$, and for model 7
$\omega_{7}^{\rm edge}=2.9$. For model 8, the density profile is infinite at
the edge and zero elsewhere. We show the analytic fit for $\omega_{8}=1$
(green). The grey area indicates the possible range of lightcurve for model 8
between the two extremes: $\omega_{8}=0$ is the bottom limit, and
$\omega_{8}=100$ as the top limit. The grey area extends up to the peak of the
lightcurve.

###  5.5 Changing vertical location

Fig. 6 shows the lightcurves of the off-plane explosion models 9-12. Among the
Gaussian density profiles, the further out explosion (model 10) at $z_{0}=2H$
explodes earlier and is brighter, as expected. On the other hand, for
radiative density profiles, the further explosion (model 12) explodes earlier,
but is dimmer than the $z=1H$ case (model 11). The reason lies in the very low
CSM mass for model 12 ($M_{\rm CSM}<10^{-3}M_{\odot}$).

The predicted breakout times are compatible with the SNEC models, although for
the $z_{0}=2~{}H$ cases, the breakout occurs slightly earlier than expected.
The reason may lie again in the low CSM mass, which limits the thermalization
of the CSM material. This could account for the large discrepancy between the
analytic and numeric lightcurves in model 12.

Figure 6: Changing the vertical location of the explosion. Gaussian profiles
are for $z_{0}=1H$ (blue, model 9) and $z_{0}=2H$ (red, model 10). Radiative
profiles are for $z_{0}=1H$ (green, model 11) and $z_{0}=2H$ (black, model
12). The line style and shaded areas have the same meaning as in Fig. 3. We
use the $\omega^{\rm edge}$ choice as in fig. 5. The effective power laws are
$\omega^{\rm edge}_{9}=3.6$, $\omega^{\rm edge}_{10}=3.17$, $\omega^{\rm
edge}_{11}=2.93$, $\omega^{\rm edge}_{12}=2.87$.

###  5.6 Comparison between the codes

Figure 7:  Comparison between the SNEC code (dashed blue) and the HORMONE code
(solid black) for all sampled models. The red curve is the analytic fit, but
starting with the peak at the average breakout time of the two velocity models
calculated in Appendix B, $t_{\rm bo}=(t_{\rm bo}|_{\rm pw}+t_{\rm bo}|_{\rm
MM99})/2$. The grey area is where the AGN luminosity can obscure the
transient. The maximal AGN luminosity is set according to Eq. 23. The model
number appears in bold black for each panel. Note that the ordinate is the
same, but the abscissa is different for each panel.

Fig. 7 shows the lightcurves of different models and compares predictions from
the Lagrangian spherically-symmetric radiative hydrodynamics SNEC code (dashed
blue) to those from the Eulerian 2D hydrodynamics HORMONE code (solid black).
The goal is to understand how important are the effects of the more realistic
2D disc geometry, including its structure evolution and morphology, since our
analytic and SNEC models assume spherical symmetry. The HORMONE simulations do
not include radiative transfer and thus radiative cooling and the transfer of
energy from the inner parts of the CSM outwards are not properly modelled.
However, this only becomes important in the later phases and the dynamics up
to and close to shock breakout should be mostly adiabatic. The vertical shock
propagation in the 2D HORMONE simulations closely resembles that of the 1D
SNEC simulations as shown in Fig. 1.

The lightcurves for the HORMONE models in Fig. 7 are all computed for a face-
on viewing angle. This is expected to be the angle where the shock breakout is
observed with the brightest luminosity. Unlike the spherically symmetric
models, the shock breakout in the 2D models starts from a small patch and
spreads out radially due to the different shock propagation times along each
polar angle. Furthermore, the effective CSM mass the shock has to interact
with before breaking out also increases with the polar angle, decreasing the
velocity at the breakout shell and making the luminosity at corresponding
annuli lower. These differences may likely explain the lower luminosities
reached in HORMONE models in Fig. 7. The peak is also smeared out over a
slightly longer duration because of the delay in shock breakout at lower
latitude angles.

For the off-plane explosions, the difference in interacting mass along each
polar angle is smaller compared to the mid-plane explosion models. Therefore,
the opening angle of the shock breakout region is wider, leading to more
similar peak luminosities between 2D and spherical models (Fig. 7 panels 9,
10).

###  5.7 Summary of the results

In summary, the typical lightcurve of a SN explosion will peak at luminosities
around $10^{44}-10^{45}\rm erg\ s^{-1}$. The earliest peak will usually be the
brightest one. Radiation-dominated or slab profiles will result in more
compact and less massive vertical layers and could produce stronger observable
events when compared to gas pressure-dominated (Gaussian) profiles. The type
of the SN (or the ejecta mass) will play a relatively minor role in the
observational signature, with the main parameter being the explosion energy.
This is the case so long as most of the kinetic outflow encounters enough CSM
material to reprocess the kinetic energy into outgoing radiation.

The strongest peak will be reached if the CSM mass is comparable to or smaller
than the ejecta mass (models 8, 12, 13, 14), although models with CSM mass
exceeding the ejecta mass by a factor of ten are also observable. Models with
too much CSM mass ($\sim 100M_{\rm ej}$) will be choked and will not be
observable. For these reasons, off-plane explosions, explosions in the inner
($\sim 10^{3}r_{s}$) regions of AGN discs and explosions in the discs of
starved AGNs will lead to the brightest transients.

##  6 Discussion

Figure 8: SNEC-based ugriz lightcurves for models 1, 7 and 14. The models
represent the CSM interaction signature for a type Ia supernova in a typical
AGN around a $10^{7}\,{\rm M}_{\odot}$ black hole for a pressure-dominated
(model 1) and radiation-dominated (model 7) atmospheres, and the same
supernova in a starved AGN (model 14). The initial blue/UV half-day-long spike
due to shock breakout (which produces the brightest signal in extreme UV /
soft X-rays) is followed by the week-long blue/white peak, followed by an
extended red/IR emission from the deposited shock energy reprocessed by the
CSM. The shape and fall-off times of the lightcurves after the main peak may
be inferred from the bolometric lightcurves.

###  6.1 Observational picture

Although the SN explosions in AGN disks are typically more luminous than
identical explosion in the absence of circumstellar material, they still need
to outshine the luminous AGN environment. Hubeny et al. (2001) studied the
spectra of AGN accretion discs and found that for a
$10^{7}\textnormal{M}_{\odot}$ disc with a viscosity parameter $\alpha=0.01$
and luminosity

$L=0.3L_{\rm Edd}=4.4\cdot
10^{44}\left(\frac{M_{\bullet}}{10^{7}\textnormal{M}_{\odot}}\right)\rm erg\
s^{-1},$ (23)

the spectrum peaks at $\nu F_{\nu}\approx 10^{44}\ \rm erg\ s^{-1}$ in the
extreme UV / soft X-ray band around $0.1\ \rm keV$. In Eq. 23, the Eddington
luminosity $L_{\rm Edd}$ assumes the electron scattering opacity.

Bolometric luminosity: For our suite of models, models 1, 6, 9 and 15 peak
around $\sim 10^{44}\rm erg\ s^{-1}$ or lower and will likely be obscured by
the AGN luminosity in the absence of clear features that make it possible to
distinguish AGN spectra from the explosion spectra. Model 4 also peaks around
$\sim 10^{44}\ \rm erg\ s^{-1}$, but the host AGN is less massive. Thus it is
a factor of $\sim 10$ more luminous than the AGN background. Models 5, 8, 10,
11 and 14 peak with luminosities $L\gtrsim 10^{45}\ \rm erg\ s^{-1}$ and are
thus most likely to be observed. Models 13 and 14 are especially promising
since they are located in a low-density environment, either in a less massive
galaxy or in a starved AGN environment, and could be more luminous than their
respective background by a factor of $\gtrsim 100$. We note that this is a
conservative estimate, since the observed luminosity can be much lower ($\sim
0.01L_{\rm Edd}$, Fabian et al., 2009), which will make all models besides
3,4, and 6 observable.

In summary, the events with the best chances to be observed will be either in
a radiation-dominated disc, and/or in a low-density environment. There are
three possibilities to achieve this among the set of models we considered: i)
in a low mass galaxy ($M\approx 10^{6}\textnormal{M}_{\odot}$, ii) in a
starved AGN of reduced density, or iii) explosions away from the midplane.

Multiband lightcurves: We show the typical multi-band lightcurves in Fig. 8.
For all the models, the first breakout will produce a short UV/blue transient,
followed by a red/IR tail. As mentioned above, the initial breakout peaks in
extreme UV / soft X-rays, so these lightcurves understate the peak bolometric
luminosity. The emission is mostly optical after the early breakout peak,
which may help to separate these events from the AGN background, which peaks
in UV/X-rays. Furthermore, the models producing sharp peaks, e.g. models 8,
12, 14, reach temperatures high enough to give rise to X-ray flares. In this
case, we expect approximately hour-long flares reaching up to $10^{43}\rm
erg/s$ in the $0.3$ – $10\,\rm keV$ band. The exceptions to this picture are
models 2 and 3, in which the shock stalls due to the large CSM mass. In this
case, the thermalised energy is emitted over hundreds of days in red/IR bands.

The presence of 56Ni in real transients will lead to additional energy
deposition on a month timescale. However, since the energy is deposited mostly
in the ejecta material located behind the CSM, the luminosity contribution
from 56Ni will be delayed by the photon diffusion time and eventually emitted
over several months or longer in red/IR bands. Since the total energy yield in
56Ni is comparable to that of the regular nuclear transients in the field and
since this energy is emitted over longer times, we find that the contribution
from radiative decay typically does not exceed $10^{43}\rm erg/s$. The AGN
will typically outshine such red/IR contributions. Conversely, in models 12-14
with low CSM mass, the CSM lightcurve after the peak becomes quickly dominated
by the contribution from the nuclear decay. In this case, the late lightcurve
resembles that of isolated supernovae.

The explosion and its ejecta are expected to be non-relativistic. It is
possible that some fraction of the electrons, especially where the temperature
is large enough, will not be in thermal equilibrium and will be accelerated to
relativistic velocities producing $\gamma$-rays, either in GRBs (Zhu et al.,
2020; Perna et al., 2021b) or in hyper-Eddington accretion-induced Bondi
explosions (Wang et al., 2021), where due to inefficient energy transport the
material is heated to very large temperatures, and relativistic shocks make a
large cavity in the AGN disc. While the GRB radiation is a prompt emission,
the radiation from Bondi explosions will be visible only in the broad-line
region, $\sim 1\ \rm yr$ after the explosion. For GRBs, following the cocoon
breakout, non relativistic ejecta also break out at later times. We find
similar peak luminosities to the estimates of Zhu et al. (2020).

###  6.2 Rates of SNe in AGN discs

Figure 9: Density snapshots from the 2D simulations for model 1. Each panel
shows a different time snapshot. Note that the left panel has a smaller box
size. Figure 10: Same as Fig. 9 but for model 9.

The radiative efficiency of an accreting AGN is $\eta=L/(\dot{M}c^{2})<1$. If
the SMBH is accreting at a fraction $L=fL_{\rm Edd}$ of the Eddington
luminosity, and the same fraction $f$ of the mass accretion rate, the SMBH
mass doubling time is

$t_{\rm doubling}=\frac{M_{\bullet}}{\dot{M}}=\frac{\eta}{f}\frac{\kappa
c}{4\pi
G}=128\left(\frac{\eta}{0.1}\right)\left(\frac{f}{0.3}\right)^{-1}\rm{Myr}.$
(24)

For a disc lifetime of $t_{\rm AGN}=10^{8}\ \rm yr$, the mass that flows into
the disc is $M_{\bullet}t_{\rm AGN}/t_{\rm doubling}\approx M_{\bullet}$.

The type Ia SN rate is 1 per few hundred solar masses of star formation in the
field (Maoz & Graur, 2017). Assuming a similar but slightly higher efficiency
for AGN discs due to an enhanced binary merger rate (e.g., one SN event per
$100\textnormal{M}_{\odot}$), and further assuming a fraction $f_{1}$ is
converted into stars, we would have $f_{1}10^{5}$ SNe for a SMBH mass of
$10^{7}M_{\odot}$ over an AGN lifetime $t_{\rm AGN}=10^{8}\ \rm yr$. Thus, the
rate per AGN is $f_{1}10^{-3}\ \rm yr^{-1}$ per AGN disc. Roughly 1% galaxies
in the local Universe have an AGN disc, so the AGN density is $10^{-4}\ \rm
Mpc^{-3}$. Thus, we expect a total of $\sim f_{1}100$ AGN SNe per $\rm
Gpc^{3}$ per year. A similar estimate for the rate of core-collapse SNe leads
to a slightly higher but otherwise comparable rate of events. This rate is
also slightly larger than the the upper limit of the binary black hole merger
rate in AGN discs, $\mathcal{R}=0.002-18\ \rm Gpc^{-3}\ yr^{-1}$, which had
been recently estimated by Gröbner et al. (2020).

A proper rate estimate should include an integral over the mass function of
the AGN in the local universe. Moreover, indirect observations of the AGN
phase duration ($10^{6}-10^{9}\ \rm yr$), and hence the duty cycle
($10^{-1}-10^{-3}$), could vary by around one order of magnitude (Martini &
Weinberg, 2001; Marconi et al., 2004). Furthermore, AGNs may not be accreting
steadily but rather experience events of rapid accretion and quiescence every
$10^{5}$ years (Schawinski et al., 2015) which could affect the star formation
rate.

###  6.3 Observable event rates

What fraction of the SNe will be observable? In order to address this
question, we need to discuss the spatial distribution of stars in AGN discs,
the mass function of the formed stars and its dependence on accretion.

Star formation in AGN discs predominantly occurs in the outer regions of the
disc at low temperatures, where the disc is gravitationally unstable. The star
formation efficiency is especially large when the midplane temperature is
between $10^{3}-10^{4}\ \rm K$, and the opacity is very low since dust is
still vaporized, but the material is mostly neutral, so electron scattering is
low.

Migration: The interaction of the stars with the massive gaseous AGN discs can
induce torques on them and cause the stars to migrate (e.g. Artymowicz et al.,
1993; McKernan et al., 2012; Bellovary et al., 2016, and references therein).
Migration can be stopped if the torque changes sign. Extensive work had been
done on migration in protoplanetary discs (Goldreich & Tremaine, 1980; Tanaka
et al., 2002; Paardekooper & Mellema, 2006; Paardekooper et al., 2010), and
migration traps were numerically reproduced (Lyra et al., 2010). Bellovary et
al. (2016) found that the outermost migration traps for AGN discs lie at
around $\sim 300r_{g}$. If most stars end up in this region, the binary
mergers and SN explosions there will be both more frequent and more
observable.

However, AGN discs are not necessarily as well behaved as protoplanetary
discs. The immense radiation and mass make them prone to flaring and
instabilities. It is uncertain whether the disc structure of SG03 and Thompson
et al. (2005) is in a steady state, and the fitting formulae for the torque as
in Paardekooper et al. (2010) may not be applicable. Moreover, Bellovary et
al. (2016) used an adiabatic index of $\gamma=5/3$. Radiation-dominated discs
with $\gamma=4/3$ will result in a different location for the migration traps,
which may be shifted outwards.

Another issue is the timescale for the migration to take place. For type I
migration, in which viscous torques are unable to form a gap in the disc, the
migration time is (Tanaka et al., 2002)

$\tau_{I}\sim\frac{M_{\bullet}^{2}}{\textnormal{M}_{\odot}\Sigma(r)r^{2}}\left(\frac{H}{r}\right)^{2}\Omega^{-1}=7\left(\frac{M_{\bullet}}{10^{7}\textnormal{M}_{\odot}}\right)^{2}\left(\frac{r}{10^{4}r_{s}}\right)^{2}\
\rm Myr$ (25)

Stars that formed beyond $10^{5}r_{s}$ or in more massive discs would require
longer than the lifetime of the disc to migrate, while less massive discs
enable faster, more efficient migration. In some very low mass discs, a gap
may be carved, and the migration timescale will be of type II, which will not
depend on the mass of the disc (or the SMBH mass). Stars may still be captured
from the nuclear star cluster onto the AGN disc.

Finally, the dynamical friction experienced by the stars could be
significantly different from these simple estimates if the accretion onto the
stars produces strong outflows (Gruzinov et al., 2020), which is almost
inevitable for the AGN environment where the accretion rates can be well above
the Eddington limit (also see below).

The spatial distribution of the progenitors of SN explosions will depend on
the migration history. Since the shocks from explosions at radial distances
beyond $r\gtrsim 10^{3}r_{s}$ will be choked, the number of observed
progenitors in the disc will be $N_{3}=N(r\lesssim 10^{3}\ r_{s})$. If
migration is effective, most of the progenitors may be close to the migration
traps, yielding a fraction of observed events $f_{2}=N_{3}/N_{\rm tot}\sim 1$,
where $N_{\rm tot}$ is the total number of the progenitors. If no migration
occurs, then assuming an underlying log-uniform distribution between
$10^{2}-10^{5}\ r_{s}$, only $f_{2}\sim 0.3$ of the events will be observable.
If most stars are either stuck in the outer regions or migrate outward, then
$f_{2}\ll 1$ could be rather low.

Accretion: Accretion could be important for changing the initial mass function
of stars in AGN discs and affecting the type Ia / core-collapse SN branching
ratios. However, neither the details of the expected super-Eddington accretion
nor the initial mass function are certain. Artymowicz et al. (1993) first
studied the growth and accretion of stars in AGN discs. More recently, Davies
& Lin (2020); Cantiello et al. (2020) suggested that significant growth can
occur. Leigh et al. (2013) also had a similar discussion, with possible
implications for depleting the gas through accretion. However, a recent
follow-up study showed that rotation and radiative feedback could limit the
final mass to be around $10\ \rm M_{\odot}$ (Jermyn et al., 2021), especially
for the innermost regions, which we mainly consider. Moreover, the extremely
optically thick environment may block the energy release via radiation and
inhibit the accretion rate. A similar situation could lead to "Bondi
explosions" (Wang et al., 2021).

Accretion onto white dwarfs might, under appropriate conditions, allow for
their growth to close to the Chandrasekhar limit and eventual explosion as
type Ia SNe Ostriker (1983), thus increasing the rate for such events.
However, here too, the exact accretion efficiency is difficult to assess.

The frequency of off-plane explosions is also quite uncertain, but probably
low. Gas dynamical friction tends to damp the orbital inclination much faster
than the migration timescale (Rein, 2012; Grishin & Perets, 2015) if the
orbital inclination with respect to the disc midplane is not too large. Highly
inclined orbits $(i\gg H/r)$ may take much longer to dissipate into the plane
of the disc, but they spend most of their time outside of the disc, and are
therefore likely to explode in relative isolation without significant CSM
interaction. Accretion feedback may change this picture (Gruzinov et al.,
2020), while close encounters with neighbouring stars may excite the
progenitors’ orbits.

To summarize, the upper limit for the rates of detectable events is
$\mathcal{R}\lesssim f_{1}f_{2}100\ \rm yr\ Gpc^{-3}$, where $f_{1}$ is the
efficiency of forming AGN gas into stars, and $f_{2}$ is the fraction of
observable events. The fraction of core-collapse SNe is uncertain and could
constrain the efficiency of accretion or the underlying mass function.

Remaining agnostic to the efficiency of accretion and its final fate, we
comment that for a Salpeter mass function $dN/dm\propto m^{-2.35}$, around
$1\%$ of stars are born massive enough to produce core-collapse SN explosions,
while for a top-heavy mass function $dN/dm\propto m^{-0.45}$ (Bartko et al.,
2010), the fraction of massive stars increases to $\sim\ 10\%$, thus
increasing the rates by an additional factor of ten.

###  6.4 Disc feedback

The 2D hydrodynamical simulations allow us to explore the morphology of the
ejecta and the impact on the local disc morphology around the explosion site.
Fig. 9 shows density snapshots of model 1. We see that at the early time of
$60\ \rm hr$, shortly before shock breakout (left panel), the blast wave is
still rather spherical, expanding within the disc. The shock front is
gradually deformed as the in-plane directions are pinched, and the shock front
propagates faster in the vertical directions, making the blast wave take a
prolate form. Later, the shock breaks out from the disc surface, and a local
hole is carved in the disc. The ejecta flow out almost spherically from the
location of shock breakout.

The overall features are qualitatively similar to recent 3D numerical
simulations of an SN explosion in a more distant radial location
($M_{\bullet}=10^{8}~{}\textnormal{M}_{\odot},r=2\cdot 10^{4}~{}r_{s}$;
Moranchel-Basurto et al., 2021). Due to the larger spatial extent and higher
interaction mass, the more distant event will not be observable. However, SN
explosions may contribute to the global viscous evolution of the AGN disc. The
2D simulations we present do not take into account the shearing motions of the
disc. In reality, the hole excavated by the SN will be elongated in the disc
rotation direction (Moranchel-Basurto et al., 2021).

Fig. 10 shows the density snapshots of the off-plane explosion, model 9. We
see that the shock breaks out earlier from the top edge, while the downward
shock propagates longer as it climbs up the density gradient. This downwards
CSM interaction creates a strong upwards reverse shock, which quickly flows
out through the upper hole. The downwards shock also eventually breaks out
from the other side of the disc. The breakout from this side of the disc is
much weaker than the upper side, and the luminosity should be much dimmer.

###  6.5 AGN flaring and variability

The significant variability observed for AGNs over a wide range of timescales
gives rise to difficulties in identifying nuclear transients and separating
such transients from various forms of inherent AGN variability.

As discussed above, the observational signatures of SNe in AGN discs resulting
from our models suggest that the SN transients will often be challenging to
observe above the background luminosity of the AGN. However, in some cases,
these AGN SNe can surpass the AGN background during their rapid rise to a
peak. The timescale for such initial high-luminosity flaring is of the order
of hours up to 1-2 days at most (see sec. 5). Identifying such short-term
transients requires high-cadence transient surveys. Moreover, many transient
surveys do not focus on AGN nuclear regions where the AGN variability
challenges identification. To the best of our knowledge, the only high-cadence
sky surveys of AGNs were done using the Kepler mission. Interestingly, these
surveys did identify one fast, days-timescale flare event KIC 1606852 (Smith
et al., 2018), which could potentially be related to the AGN SNe we discuss.
However, it might be too long or be the result of other explosive processes
(Smith et al., 2018). The timescale of other nuclear transients, such as
candidate tidal disruption events and changing-look AGNs, are far longer than
the peak timescale of AGN SNe and are not likely to be related. Future high-
cadence transient surveys may be able to identify AGN SNe.

##  7 Summary

The AGN disc environment around SMBHs has dense stellar populations, many of
which are expected to be embedded in the disc. As stars end their lives in
supernova explosions, the dense AGN environment makes a fertile ground for
unique transient events. Understanding such explosions and their properties
opens a window to the physics of AGN discs and their interactions with the
supernova progenitors.

Motivated by superluminous explosions in circumstellar material (CSM), we have
developed an analytical model for the expected time of breakout and subsequent
lightcurve, depending on the disc and progenitor properties. We validated this
model with extensive numerical simulations. We found that the typical peak
luminosity for such events may reach a few times $10^{45}\ \rm erg\ s^{-1}$.
The most energetic events are also the quickest to break out, ranging from
hours to several days, until they become too faint to be detected.

The most luminous explosions are generally found either in radiation-dominated
discs, at larger explosion energies, or in regions with reduced density, such
as off-plane explosions, low mass SMBH, or starved AGNs with reduced density.
The latter two (i.e. low SMBH mass and starved AGN discs) also have better
chances of being observed due to reduced background AGN luminosity. The
initial breakout events should be dominated by the blue bands, while at later
times, the lightcurves will be dominated by the red and infrared bands. The
upper limit for the event rate is $\mathcal{R}\lesssim 100\ \rm yr\ Gpc^{-3}$,
where optimal star formation and observational conditions are assumed. The
actual rate could be lower by orders of magnitude.

In concise form, our study may be summarised as follows: We exploded stars in
an AGN. The photons can’t get out, the outflow can. The gas converts shocks
into brighter SN. We may see a few hundred, if we frequently scan.

## Data availability

The simulations underlying this article will be shared on reasonable request
to the corresponding author.

## Acknowledgements

We thank the referee, Pablo Fabián Velázquez, for valuable comments on the
manuscript. We thank Almog Yalinewich, Ari Laor and Yossef Zenati for
stimulating discussions. EG and HBP acknowledge support for this project from
the European Union’s Horizon 2020 research and innovation program under grant
agreement No 865932-ERC-SNeX. IM is a recipient of the Australian Research
Council Future Fellowship FT190100574.

## References

  * Artymowicz et al. (1993) Artymowicz P., Lin D. N. C., Wampler E. J., 1993, ApJ, 409, 592
  * Bartko et al. (2010) Bartko H., et al., 2010, ApJ, 708, 834
  * Baruteau et al. (2011) Baruteau C., Cuadra J., Lin D. N. C., 2011, ApJ, 726, 28
  * Bellovary et al. (2016) Bellovary J. M., Mac Low M.-M., McKernan B., Ford K. E. S., 2016, ApJ, 819, L17
  * Cantiello et al. (2020) Cantiello M., Jermyn A. S., Lin D. N. C., 2020, arXiv e-prints, p. arXiv:2009.03936
  * Cheng & Wang (1999) Cheng K. S., Wang J.-M., 1999, ApJ, 521, 502
  * Chevalier & Irwin (2011) Chevalier R. A., Irwin C. M., 2011, ApJ, 729, L6
  * Davies & Lin (2020) Davies M. B., Lin D. N. C., 2020, MNRAS, 498, 3452
  * Davis & Tchekhovskoy (2020) Davis S. W., Tchekhovskoy A., 2020, ARA&A, 58, 407
  * Dessart et al. (2015) Dessart L., Audit E., Hillier D. J., 2015, MNRAS, 449, 4304
  * Dittmann & Miller (2020) Dittmann A. J., Miller M. C., 2020, MNRAS, 493, 3732
  * Fabian et al. (2009) Fabian A. C., Vasudevan R. V., Mushotzky R. F., Winter L. M., Reynolds C. S., 2009, MNRAS, 394, L89
  * Ferguson et al. (2005) Ferguson J. W., Alexander D. R., Allard F., Barman T., Bodnarik J. G., Hauschildt P. H., Heffner-Wong A., Tamanai A., 2005, ApJ, 623, 585
  * Ferrara & Tolstoy (2000) Ferrara A., Tolstoy E., 2000, MNRAS, 313, 291
  * Gal-Yam (2019) Gal-Yam A., 2019, ARA&A, 57, 305
  * Ginzburg & Balberg (2012) Ginzburg S., Balberg S., 2012, ApJ, 757, 178
  * Goldreich & Tremaine (1980) Goldreich P., Tremaine S., 1980, ApJ, 241, 425
  * Goldreich et al. (2002) Goldreich P., Lithwick Y., Sari R., 2002, Nature, 420, 643
  * Grishin & Perets (2015) Grishin E., Perets H. B., 2015, ApJ, 811, 54
  * Gröbner et al. (2020) Gröbner M., Ishibashi W., Tiwari S., Haney M., Jetzer P., 2020, A&A, 638, A119
  * Gruzinov et al. (2020) Gruzinov A., Levin Y., Matzner C. D., 2020, MNRAS, 492, 2755
  * Hirai et al. (2016) Hirai R., Nagakura H., Okawa H., Fujisawa K., 2016, Phys. Rev. D, 93, 083006
  * Hirose et al. (2009) Hirose S., Krolik J. H., Blaes O., 2009, ApJ, 691, 16
  * Hubeny et al. (2001) Hubeny I., Blaes O., Krolik J. H., Agol E., 2001, ApJ, 559, 680
  * Iglesias & Rogers (1996) Iglesias C. A., Rogers F. J., 1996, ApJ, 464, 943
  * Janka (2012) Janka H.-T., 2012, Annual Review of Nuclear and Particle Science, 62, 407
  * Jermyn et al. (2021) Jermyn A. S., Dittmann A. J., Cantiello M., Perna R., 2021, arXiv e-prints, p. arXiv:2102.13114
  * Jiang et al. (2013) Jiang Y.-F., Stone J. M., Davis S. W., 2013, ApJ, 778, 65
  * Jiang et al. (2016) Jiang Y.-F., Davis S. W., Stone J. M., 2016, ApJ, 827, 10
  * Leigh et al. (2013) Leigh N. W. C., Böker T., Maccarone T. J., Perets H. B., 2013, MNRAS, 429, 2997
  * Lynden-Bell (1969) Lynden-Bell D., 1969, Nature, 223, 690
  * Lyra et al. (2010) Lyra W., Paardekooper S.-J., Mac Low M.-M., 2010, ApJ, 715, L68
  * Maoz & Graur (2017) Maoz D., Graur O., 2017, ApJ, 848, 25
  * Marconi et al. (2004) Marconi A., Risaliti G., Gilli R., Hunt L. K., Maiolino R., Salvati M., 2004, MNRAS, 351, 169
  * Martini & Weinberg (2001) Martini P., Weinberg D. H., 2001, ApJ, 547, 12
  * Matzner & McKee (1999) Matzner C. D., McKee C. F., 1999, ApJ, 510, 379
  * McKernan et al. (2012) McKernan B., Ford K. E. S., Lyra W., Perets H. B., 2012, MNRAS, 425, 460
  * McKernan et al. (2020) McKernan B., Ford K. E. S., O’Shaughnessy R., 2020, MNRAS, 498, 4088
  * Miralda-Escudé & Gould (2000) Miralda-Escudé J., Gould A., 2000, ApJ, 545, 847
  * Moranchel-Basurto et al. (2021) Moranchel-Basurto A., Sánchez-Salcedo F. J., Chametla R. O., Velázquez P. F., 2021, ApJ, 906, 15
  * Moriya et al. (2013) Moriya T. J., Maeda K., Taddia F., Sollerman J., Blinnikov S. I., Sorokina E. I., 2013, MNRAS, 435, 1520
  * Morozova et al. (2015a) Morozova V., Ott C. D., Piro A. L., 2015a, SNEC: SuperNova Explosion Code (ascl:1505.033)
  * Morozova et al. (2015b) Morozova V., Piro A. L., Renzo M., Ott C. D., Clausen D., Couch S. M., Ellis J., Roberts L. F., 2015b, ApJ, 814, 63
  * Morozova et al. (2017) Morozova V., Piro A. L., Valenti S., 2017, ApJ, 838, 28
  * Nakar & Piro (2014) Nakar E., Piro A. L., 2014, ApJ, 788, 193
  * Nakar & Sari (2010) Nakar E., Sari R., 2010, ApJ, 725, 904
  * Nomoto et al. (1984) Nomoto K., Thielemann F. K., Yokoi K., 1984, ApJ, 286, 644
  * Ostriker (1983) Ostriker J. P., 1983, ApJ, 273, 99
  * Paardekooper & Mellema (2006) Paardekooper S. J., Mellema G., 2006, A&A, 459, L17
  * Paardekooper et al. (2010) Paardekooper S. J., Baruteau C., Crida A., Kley W., 2010, MNRAS, 401, 1950
  * Paczynski (1978) Paczynski B., 1978, Acta Astron., 28, 91
  * Perets et al. (2016) Perets H. B., Li Z., Lombardi James C. J., Milcarek Stephen R. J., 2016, ApJ, 823, 113
  * Perna et al. (2021a) Perna R., Tagawa H., Haiman Z., Bartos I., 2021a, arXiv e-prints, p. arXiv:2103.10963
  * Perna et al. (2021b) Perna R., Lazzati D., Cantiello M., 2021b, ApJ, 906, L7
  * Piro et al. (2020) Piro A. L., Haynie A., Yao Y., 2020, arXiv e-prints, p. arXiv:2007.08543
  * Rein (2012) Rein H., 2012, MNRAS, 422, 3611
  * Rozyczka et al. (1995) Rozyczka M., Bodenheimer P., Lin D. N. C., 1995, MNRAS, 276, 597
  * Sakurai (1960) Sakurai A., 1960, Communications on Pure and Applied Mathematics, 13, 353
  * Schawinski et al. (2015) Schawinski K., Koss M., Berney S., Sartori L. F., 2015, MNRAS, 451, 2517
  * Sedov (1946) Sedov L. I., 1946, Prikl. Mat. Mekh, 10, 241
  * Semenov et al. (2003) Semenov D., Henning T., Helling C., Ilgner M., Sedlmayr E., 2003, A&A, 410, 611
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 500, 33
  * Shakura & Sunyaev (1976) Shakura N., Sunyaev R., 1976, Monthly Notices of the Royal Astronomical Society, 175, 613
  * Sirko & Goodman (2003) Sirko E., Goodman J., 2003, MNRAS, 341, 501
  * Smith et al. (2018) Smith K. L., Mushotzky R. F., Boyd P. T., Malkan M., Howell S. B., Gelino D. M., 2018, ApJ, 857, 141
  * Suzuki et al. (2016) Suzuki A., Maeda K., Shigeyama T., 2016, ApJ, 825, 92
  * Suzuki et al. (2019) Suzuki A., Moriya T. J., Takiwaki T., 2019, ApJ, 887, 249
  * Suzuki et al. (2021) Suzuki A., Nicholl M., Moriya T. J., Takiwaki T., 2021, ApJ, 908, 99
  * Syer et al. (1991) Syer D., Clarke C. J., Rees M. J., 1991, MNRAS, 250, 505
  * Tagawa et al. (2020) Tagawa H., Haiman Z., Kocsis B., 2020, ApJ, 898, 25
  * Takei & Shigeyama (2020) Takei Y., Shigeyama T., 2020, PASJ, 72, 67
  * Tanaka et al. (2002) Tanaka H., Takeuchi T., Ward W. R., 2002, ApJ, 565, 1257
  * Taylor (1950) Taylor G., 1950, Proceedings of the Royal Society of London Series A, 201, 159
  * Thompson et al. (2005) Thompson T. A., Quataert E., Murray N., 2005, ApJ, 630, 167
  * Tsuna et al. (2019) Tsuna D., Kashiyama K., Shigeyama T., 2019, ApJ, 884, 87
  * Utrobin (2007) Utrobin V. P., 2007, in Immler S., Weiler K., McCray R., eds, American Institute of Physics Conference Series Vol. 937, Supernova 1987A: 20 Years After: Supernovae and Gamma-Ray Bursters. pp 25–32, doi:10.1063/1.3682879
  * Wang et al. (2021) Wang J.-M., Liu J.-R., Ho L. C., Du P., 2021, arXiv e-prints, p. arXiv:2103.07708
  * Woosley et al. (2007) Woosley S. E., Blinnikov S., Heger A., 2007, Nature, 450, 390
  * Yalinewich & Matzner (2019) Yalinewich A., Matzner C. D., 2019, MNRAS, 490, 312
  * Zhu et al. (2020) Zhu J.-P., Zhang B., Yu Y.-W., Gao H., 2020, arXiv e-prints, p. arXiv:2011.08428
  * Zhu et al. (2021a) Zhu J.-P., Wang K., Zhang B., Yang Y.-P., Yu Y.-W., Gao H., 2021a, arXiv e-prints, p. arXiv:2103.00789
  * Zhu et al. (2021b) Zhu J.-P., Yang Y.-P., Zhang B., Liu L.-D., Yu Y.-W., Gao H., 2021b, arXiv e-prints, p. arXiv:2104.09389

##  Appendix A Vertical profile for radiation-dominated pressure

Consider the radiation pressure $P=aT^{4}/3=K\rho^{4/3}$, where
$a=4\sigma_{\rm SB}/c$, $\sigma_{\rm SB}$ is the Stefan-Boltzmann constant and
$c$ is the speed of light. We assume a polytropic equation of state where
$P=K\rho^{\gamma}$. The sound speed is given by
$c_{s}=\sqrt{dP/d\rho}=\sqrt{\gamma K\rho^{\gamma-1}}=\sqrt{\gamma P/\rho}$.
The vertical hydrostatic equilibrium equation is

$\frac{dP}{dh}=-\rho g_{z}=K\gamma\rho^{\gamma-1}\frac{d\rho}{dh}$ (26)

where $g_{z}=GM_{\rm SMBH}h/(r^{2}+h^{2})^{3/2}$ is the $z$ component of the
tidal gravitational acceleration. Since the disc is thin ($h\ll R$), we can
approximate it as $g_{z}\approx\Omega^{2}h$ where $\Omega^{2}=GM_{\rm
SMBH}/r^{3}$ is the Keplerian frequency. We have the ordinary differential
equation

$\rho^{\gamma-2}d\rho=-\frac{\Omega^{2}}{\gamma K}.$ (27)

Integrating, we have

$\frac{\rho^{\gamma-1}}{\gamma-1}=-\frac{\Omega^{2}}{2\gamma K}z^{2}+C.$ (28)

To determine $C$ we use the initial condition at $h=0$

$\rho^{\gamma-1}=(\gamma-1)C=\rho_{0}^{\gamma-1}\implies
C=\rho_{0}^{\gamma-1}/(\gamma-1).$ (29)

Finally, the vertical profile is

$\rho(h)=\rho_{0}\left[1-(\gamma-1)\frac{\Omega^{2}}{2\gamma
K\rho_{0}^{\gamma-1}}h^{2}\right]^{1/(\gamma-1)}.$ (30)

Note that the scale height $H$ is defined as

$H^{-2}=\frac{\Omega^{2}}{c_{s}^{2}}=\frac{\Omega^{2}}{\gamma
K\rho_{0}^{\gamma-1}},$ (31)

which leads to the final expression

$\rho_{\rm
rad}(h)=\rho_{0}\left[1-(\gamma-1)\frac{h^{2}}{2H^{2}}\right]^{1/(\gamma-1)}.$
(32)

Note that the radiation pressure-dominated profile is steeper than the
Gaussian profile $\rho_{\rm gas}\propto\exp(-h^{2}/2H^{2})$ for gas dominated
pressure $P\propto\rho$, and goes to zero at $h/H=\sqrt{2/(\gamma-1)}$, or
$h/H=\sqrt{6}\approx 2.45$.

##  Appendix B Breakout times in AGN discs

Gaussian density profile: For the piecewise velocity (Eq. 12) in the Gaussian
density profile (Eq. 4) the integral in Eq. 14 can be split into three
integrals. The first one over the free expansion phase is trivial, while the
other two involve integrals of the form $\int x^{3/2}e^{-\alpha
x^{2}/2}dx=-2^{1/4}\Gamma(5/4,\alpha x^{2}/2)/\alpha^{5/4}$, where
$\Gamma(s,x)=\int_{x}^{\infty}y^{s-1}e^{-y}dy$ is the upper incomplete Gamma
function. We get

$\displaystyle\left.\frac{t_{\rm bo}}{H/v_{0}}\right|_{\rm pw}=$
$\displaystyle\zeta_{1}+\frac{2^{1/4}\mathcal{A}_{1}}{\zeta_{1}^{3/2}\mu^{\prime
5/4}}\left[\Gamma\left(\frac{5}{4},\frac{\mu^{\prime}\zeta_{1}^{2}}{2}\right)-\Gamma\left(\frac{5}{4},\frac{\mu^{\prime}\zeta_{2}^{2}}{2}\right)\right]$
$\displaystyle+$
$\displaystyle\frac{2^{1/4}\mathcal{A}_{2}}{\zeta_{1}^{3/2}\mu{}^{5/4}}\left[\Gamma\left(\frac{5}{4},\frac{\mu\zeta_{2}^{2}}{2}\right)-\Gamma\left(\frac{5}{4},\frac{\mu\zeta_{{\rm
bo}}^{2}}{2}\right)\right]$ (33)

where $\zeta_{i},\zeta_{\rm bo}=z_{i}/H,z_{\rm bo}/H$ are dimensionless length
scales, $\mathcal{A}_{1}=(\rho_{0}/\rho_{1})^{\mu^{\prime}}$ and
$\mathcal{A}_{2}=(\rho_{2}/\rho_{1})^{\mu^{\prime}}(\rho_{0}/\rho_{1})^{\mu}$
are coefficients that depend on the density at each location, which comes from
the continuity requiremtns of the piecewise velocity in Eq. (12). For the
parameters of model 1, we get $t_{\rm bo}=3.17\ \rm d$.

For the MM99 velocity model (Eq. 13) in the Gaussian profile, the breakout
time is given by

$\left.t_{\rm bo}\right|_{\rm MM99}=\frac{H}{v_{0}}\intop_{0}^{\zeta_{{\rm
bo}}}\left(1+\frac{\rho_{0}H^{3}}{M_{{\rm
ej}}}\zeta^{3}\right)^{1/2}e^{-\mu\zeta^{2}/2}d\zeta.$ (34)

This is not an analytic expression, but can be evaluated numerically. For
model 1, the breakout time is around $2.55\ \rm d$, which is shorter by a
factor of $\sim 17\%$. Generally, as seen in Fig. 1, the shock velocity
obtained by simulations is somewhere in between the piecewise model and the
MM99 model. We compare to the numerical results in sec. 4.

Radiation-dominated density profile: For the piecewise velocity in the
radiation-dominated profile (Eq. 5), the integrals are similar. They have the
form

$J(y,\alpha)=\intop
y^{3/2}\left(1-\frac{y^{2}}{6}\right)^{3\alpha}dy=\frac{2}{5}y^{2/5}{}_{2}F_{1}\left(\frac{5}{4},-3\alpha,\frac{9}{4},\frac{y^{2}}{6}\right),$
(35)

where ${}_{2}F_{1}(a,b,c,x)$ is the hypergeometric function. The breakout time
is

$\displaystyle\left.\frac{t_{\rm bo}}{H/v_{0}}\right|_{\rm pw}=$
$\displaystyle\zeta_{1}+\frac{\mathcal{A}_{1}}{\zeta_{1}^{3/2}}\left[J(\zeta_{2},3\mu^{\prime})-J(\zeta_{1},3\mu^{\prime})\right]$
$\displaystyle+$
$\displaystyle\frac{\mathcal{A}_{2}}{\zeta_{1}^{3/2}}\left[J(\zeta_{\rm
bo},3\mu)-J(\zeta_{2},3\mu)\right].$ (36)

For the MM99 velocity model in the radiation-dominated profile, the breakout
time is given by

$\left.t_{\rm bo}\right|_{\rm MM99}=\frac{H}{v_{0}}\intop_{0}^{\zeta_{{\rm
bo}}}\left(1+\frac{\rho_{0}H^{3}}{M_{{\rm
ej}}}\zeta^{3}\right)^{1/2}\left(1-\frac{\zeta^{2}}{6}\right)^{3\mu}d\zeta,$
(37)

which is, again, not analytic but can be estimated numerically.

Off-plane explosions: For off-plane explosions the integrals are the same but
the integration domain is $(\zeta_{0},\zeta_{\rm bo})$, where $\zeta_{0}$ is
the explosion site above the midplane. For the MM99 model this is the only
change.

For the piecewise model, we divide the domain into three regions,
corresponding to the piecewise velocity solution in Eq. 12, namely

$t_{\rm bo}=\frac{H}{v_{0}}(I_{1}+I_{2}+I_{3}).$ (38)

Here, $I_{1}=\zeta_{1}-\zeta_{0}$, is the integral over the (constant velocity
$v_{0}$) free expansion phase; $I_{2}$ is the integral between $\zeta_{1}$ and
$\zeta_{2}$, where the velocity scales with the density with a power law
$v\propto\rho^{-\mu^{\prime}}$, where $\mu^{\prime}=\mu+1/5$; and $I_{3}$ is
the integral between $\zeta_{2}$ and $\zeta_{\rm bo}$, where the velocity
scales with the density with a power law $v\propto\rho^{-\mu}$.

For $\zeta_{\rm bo}>\zeta_{2}>\zeta_{1}>\zeta_{0}$, the result is essentially
the same as in Eq. 33, but with $I_{1}$ replacing $\zeta_{1}$ in the first
term. This is the case for model 9.

For $\zeta_{1}<\zeta_{\rm bo}<\zeta_{2}$, the shock breaks out before the
Sakurai acceleration begins. This leads to $I_{3}=0$ and the domain of
integration of $I_{2}$ is $(\zeta_{1},\zeta_{\rm bo})$. This is the case for
models 10 and 11.

For $\zeta_{\rm bo}<\zeta_{1}$, the shock breaks out before the free expansion
phase ends. In this case, we only have $I_{1}=\zeta_{\rm bo}-\zeta_{1}$. This
is the case for model 12.

##  Appendix C CSM mass

Here we explicitly calculate the integral in Eq. (21) for various models. We
will start with mid-plane explosions and hence fix $z_{\rm min}=0$.

For a flat profile, $z_{\rm max}=H$ and the CSM mass is simply $M_{\rm
CSM}=4\pi\rho_{0}H^{3}/3$.

For the Gaussian profile we have

$\displaystyle M_{\rm CSM}(z_{\rm max})$
$\displaystyle=4\pi\rho_{0}\intop_{0}^{z_{\rm max}}e^{-z^{2}/2H^{2}}z^{2}dz$
$\displaystyle=4\pi\rho_{0}H^{3}\left[\sqrt{\frac{\pi}{2}}{\rm
erf}\left(\zeta_{\rm max}/\sqrt{2}\right)-\zeta_{\rm max}e^{-\zeta_{\rm
max}^{2}/2}\right]$ (39)

where $\zeta_{\rm max}=z_{\rm max}/H.$ Note that for $\zeta_{\rm
max}\to\infty$ we get $M_{\rm{CSM},\infty}=(2\pi)^{3/2}\rho_{0}H^{3}$ while
for small $\zeta_{\rm max}\ll 1$ it is essentially as the uniform density
case, $M_{\rm CSM}(z_{\rm max}\ll H)\approx 4\pi\rho_{0}z_{\rm max}^{3}/3$.
For all practical purposes, the value at the photosphere is close enough to
the limiting value.

For the radiation-dominated profile, we have

$\displaystyle M_{\rm CSM}(z_{\rm max})$
$\displaystyle=4\pi\rho_{0}\intop_{0}^{z_{\rm
max}}\left(1-\frac{z^{2}}{6H^{2}}\right)^{3}z^{2}dz$
$\displaystyle=4\pi\rho_{0}H^{3}\left[\frac{\zeta_{\rm
max}^{3}}{3}-\frac{\zeta_{\rm max}^{5}}{10}+\frac{\zeta_{\rm
max}^{7}}{84}-\frac{\zeta^{9}_{\rm max}}{1944}\right]$ (40)

For the maximal value $\zeta_{\rm max}=\sqrt{6}$ we have $M(z_{\rm
max})=4\pi\rho_{0}H^{3}\times\frac{32}{35}\sqrt{\frac{2}{3}}=9.38\rho_{0}H^{3}$.

Off plane explosions: For an off plane explosion which occurs at $z_{\rm\rm
min}$, the density ranges from $\rho(z_{\rm min})$ to
$\rho(z_{\rm{max}}-z_{\rm min})$, which the vertical scale increases from $0$
to $z_{\rm max}-z_{\rm min}$. The resulting integral is

$M_{\rm CSM}(z_{\rm min},z_{\rm max})=4\pi\intop_{0}^{z_{\rm max}-z_{\rm
min}}\rho(z+z_{\rm min})z^{2}dz.$ (41)

For the Gaussian profile we have

$\displaystyle\frac{M_{\rm CSM}}{4\pi\rho_{0}H^{3}}=$
$\displaystyle\sqrt{\frac{\pi}{2}}(1+\zeta_{{\rm min}}^{2})\left[{\rm
erf}\left(\frac{\zeta_{{\rm max}}}{\sqrt{2}}\right)-{\rm
erf}\left(\frac{\zeta_{{\rm min}}}{\sqrt{2}}\right)\right]$ $\displaystyle-$
$\displaystyle e^{-\zeta_{{\rm min}}^{2}/2}\zeta_{{\rm min}}-e^{-\zeta_{{\rm
max}}^{2}/2}(\zeta_{{\rm max}}-2\zeta_{{\rm min}})$ (42)

where $\zeta_{\rm max/min}=z_{\rm max/min}/H$. For model 9 we use $\zeta_{\rm
max}=\infty$ and $\zeta_{\rm min}=1$, and the resulting mass is $M_{\rm
CSM}=1.92M_{\odot}$. For model 10 we use $\zeta_{\rm max}=\infty$ and
$\zeta_{\rm min}=2$, and the resulting mass is $M_{\rm CSM}=0.15M_{\odot}$.

For the radiation dominated profile we have

$\displaystyle\frac{M_{\rm CSM}}{4\pi\rho_{0}H^{3}}=$
$\displaystyle\intop_{\zeta_{\rm min}}^{\sqrt{6}}(\zeta-\zeta_{{\rm
min}})^{2}\left(1-\frac{\zeta^{2}}{6}\right)^{3}d\zeta=\frac{\zeta_{\rm
min}^{9}}{54432}-\frac{\zeta_{\rm min}^{7}}{1260}$ $\displaystyle+$
$\displaystyle\frac{\zeta_{\rm min}^{5}}{60}-\frac{\zeta_{\rm
min}^{3}}{3}+\frac{16\sqrt{6}\zeta_{\rm min}^{2}}{35}-\frac{3\zeta_{\rm
min}}{2}+\frac{32\sqrt{2}}{35\sqrt{3}}.$ (43)

For model 11 we use $\zeta_{\rm min}=1$, and the resulting mass is $M_{\rm
CSM}=0.49M_{\odot}$. For model 12 we use $\zeta_{\rm min}=2$, and the
resulting mass is $M_{\rm CSM}=6.4\cdot 10^{-4}M_{\odot}$.

##  Appendix D lightcurve computation in the 2D models

We first define an observer plane placed above the disc, parallel to the disc
plane. From behind the disc, we solve the simple radiative transfer equations
along parallel lines of sight to obtain the intensity at each position on the
observer plane by computing the sequence

$I(z+\Delta z)=I(z)e^{-\Delta\tau}+S(1-e^{-\Delta\tau}),$ (44)

where $I$ is the intensity, $S$ is the source function and
$\Delta\tau=\rho\kappa\Delta z$ is the optical depth along a small distance
$\Delta z$. We assume black body radiation for the source function ($S=\sigma
T^{4}/\pi$) and only compute the bolometric luminosity. For opacity, we only
assume electron scattering ($\kappa=0.2(1+X)$ cm2 g-1, where $X$ is the
hydrogen fraction). We do not take into account the effect of recombination in
order to avoid making the disc unreasonably optically thin due to the way we
set the initial disc temperature. This choice should not affect the early
evolution of the lightcurve, at least in the planar phase. The intensity is
then integrated over the whole observer plane and multiplied by $4\pi$ to
obtain the luminosity. We refer the reader to Suzuki et al. (2016) for more
details of the idea behind the calculation procedure.

Because radiation transport is not coupled to the hydrodynamics in our 2D
simulations, this methodology does not account for the diffusion and radiative
cooling in the ejecta. Therefore, the lightcurves obtained this way may not be
energetically self-consistent. Additionally, the temperature in the outer
parts of the ejecta may be overestimated. If the outer parts of the ejecta
have sufficiently cooled, the atoms may recombine and lower the opacity,
affecting the later lightcurve. Nevertheless, this estimate should give us a
rough idea of how the shock breakout may be observed.