# Electromagnetic Cascade Emission from Neutrino-Coincident Tidal Disruption
Events

Chengchao Yuan (袁成超) Deutsches Elektronen-Synchrotron DESY, Platanenallee 6,
15738 Zeuthen, Germany  Walter Winter Deutsches Elektronen-Synchrotron DESY,
Platanenallee 6, 15738 Zeuthen, Germany

###### Abstract

The potential association between Tidal Disruption Events (TDEs) and high-
energy astrophysical neutrinos implies the acceleration of cosmic rays. These
accelerated particles will initiate electromagnetic (EM) cascades spanning
from keV to GeV energies by the processes related to neutrino production. We
model the EM cascade and neutrino emissions by numerically solving the time-
dependent transport equations and discuss the implications for AT2019dsg and
AT2019fdr in the X-ray and $\gamma$-ray bands. We show that the $\gamma$-ray
constraints from _Fermi_ can constrain the size of the radiation zone and the
maximum energy of injected protons, and that the corresponding expected
neutrino event numbers in follow-up searches are limited to be less than about
0.1. Depending on the efficiency of $p\gamma$ interactions, the X-ray and
$\gamma$-ray signals can be expected closer to the peak of the optical-
ultraviolet (OUV) luminosity, or to the time of the neutrino production.

Tidal disruption; Radiative processes; Neutrino astronomy

##  1 Introduction

Recent observations have revealed that TDEs can produce intense flares of
radiation lasting from months to years, powered by infalling material from a
tidally disrupted star. Observationally, TDEs commonly exhibit thermal
emissions predominantly in the OUV range, with sub-populations also observed
in the X-ray and radio bands (see, e.g., Saxton et al., 2020; Alexander et
al., 2020; van Velzen et al., 2021a). The TDE catalog keeps rapidly expanding
as more and more TDEs have been continuously identified by the Zwicky
Transient Facility (ZTF, e.g., Hammerstein et al., 2023).

Three TDE candidates, AT2019dsg (Stein et al., 2021), AT2019fdr (Reusch et
al., 2022) and AT2019aalc (van Velzen et al., 2021b), are potentially
correlated to the three IceCube astrophysical neutrino events, IC191001A,
IC200530A and IC19119A, respectively. Two of these neutrino associations
(AT2019dsg and AT2019fdr) have been found by follow-up searches. It has been
pointed out in van Velzen et al. (2021b) that these associations have a strong
dust echo in common – which is delayed re-processed radiation from the OUV and
X-ray ranges into the infrared (IR) range by surrounding dust; this has led to
the identification of the third candidate (AT2019aalc). Remarkably, these
three TDEs are located in the 90% containment boxes of the corresponding
neutrino events, confirmed in the recent catalog of IceCube neutrino tracks
(Abbasi et al., 2023), which consolidates the multi-messenger correlations.
Apart from a strong dust echoes and OUV luminosities, AT2019dsg and AT2019fdr
have been observed in X-rays, and the neutrinos have come delayed by half a
year to a year with respect to the OUV peak. The neutrino counterparts of
these three TDEs have a profound impact on TDE models and bring significant
momentum to related researches since astrophysical neutrinos are primarily
generated in hadronic processes requiring the acceleration of cosmic-ray
primaries.

Extensive efforts have been made on theoretical modeling, simulations and
observations to establish a unified physical picture for TDEs. In this
picture, the accretion disks, sub-relativistic outflows or winds, and possibly
a jet and a dust torus are typically included (see, e.g., Dai et al., 2018).
Motivated by the observation of the luminous jetted TDE Swift J1644+57 in 2011
(Burrows et al., 2011), relativistic jets have been proposed as one promising
origin of the non-thermal emissions of TDEs (Wang et al., 2011; Wang & Liu,
2016; Dai & Fang, 2017; Lunardini & Winter, 2017; Senno et al., 2017). In the
jetted models, the energy released is beamed in a narrow region where the
power is sufficient high to explain the intense emissions. Another jetted TDE
candidate, AT2022cmc (Andreoni et al., 2022), has been recently reported to
exhibit a very high non-thermal X-ray luminosity $\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}10^{47}-10^{48}~{}\rm
erg\ s^{-1}$, which is likely produced by a collimated jet with an opening
angle $\theta_{j}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}1^{\circ}$ (Pasham et al., 2023). It was
also pointed out that only around one percent of TDEs are expected to have
relativistic jets; however, this small fraction could be consistent with the
fraction of neutrino-emitting TDEs from diffuse flux considerations (Winter &
Lunardini, 2023).

Apart from jets, quasi-isotropic emissions from accretion disks (Hayasaki &
Yamazaki, 2019), wide-angle outflows/winds (Fang et al., 2020) or tidal stream
interactions (Dai et al., 2015; Hayasaki & Yamazaki, 2019) have been
considered, particularly for TDEs that do not exhibit direct jet signatures,
such as radio signals. The TDE jets, disks, and winds can be efficient cosmic-
ray (CR) accelerators, potentially producing ultra-high-energy cosmic rays
(UHECRs, Farrar & Gruzinov, 2009; Farrar & Piran, 2014; Zhang et al., 2017;
Guépin et al., 2018; Biehl et al., 2018b). Consequently, these sites may also
serve as promising neutrino emitters. For AT2019dsg, various models such as
jets, disks, corona and hidden winds and outflow-cloud interactions have been
proposed to explain the neutrino counterparts (Liu et al., 2020; Winter &
Lunardini, 2021; Wu et al., 2022; Murase et al., 2020; Hayasaki, 2021).
Similarly, for AT2019fdr, Reusch et al. (2022) proposed the corona, hidden
wind and jet models, whereas a disk model for all three neutrino-coincident
TDEs has been proposed by van Velzen et al. (2021b). Motivated by the absence
of definitive evidence for jet signatures (e.g., Mohan et al., 2022), Winter &
Lunardini (2023) have recently presented unified time-dependent
interpretations of the neutrino emissions by considering three models (named
M-IR, M-OUV, and M-X) where thermal IR, OUV and X-ray photons respectively
dominate the neutrino productions depending on the maximally available proton
acceleration energy. An alternative model, proposed by Zheng et al. (2022),
has considered that the jets are choked inside the ejecta.

High-energy astrophysical neutrinos typically originate from charged pions
produced in cosmic ray interactions with matter or radiation. While the
charged pion decay chain also leads to electrons and positrons as decay
products, neutral pions co-produced with the charged pions directly decay into
gamma-rays. Therefore, a similar amount of energy is expected to be injected
into neutrinos and the electromagnetic (EM) cascade in source, where it is a
priori not clear where in energy the EM signatures dominate, and what the
time-dependent behavior of the EM cascade is – these questions require
dedicated theoretical modeling, which we perform in this study. A joint
analysis of the non-thermal X-ray and $\gamma$-ray observations with the
neutrino detection will then provide further insights: The X-ray light curve
for AT2019dsg was measured using Swift X-ray Telescope (XRT, Burrows et al.,
2005), XMM-Newton (XMM, Jansen et al., 2001), and NICER (Stein et al., 2021;
Cannizzaro et al., 2021). As for AT2019fdr, upper limits were obtained from
XRT and eROSITA (Saxton et al., 2020; Reusch et al., 2022; Predehl et al.,
2021). The _Fermi_ Large Area Telescope (Atwood et al., 2009) did not detect
significant $\gamma$-ray signals from both of them, resulting in a time-
dependent up limit of $10^{-9}-10^{-8}~{}\rm GeV~{}cm^{-2}~{}s^{-1}$ (Stein et
al., 2021). Given the relatively incomplete X-ray and $\gamma$-ray datasets
for AT2019aalc (van Velzen et al., 2021b), we choose AT2019dsg and AT2019fdr
as two prototypes to study the implications of the EM cascades.

We adopt the quasi-isotropic model proposed by Winter & Lunardini (2023) and
specifically study the time-dependent EM cascade emissions in radiation zones
characterized generically by their sizes and the maximum energies of injected
protons. Compared to Winter & Lunardini (2023), where the simulations were
based on the NeuCosmA (Neutrinos from Cosmic Accelerators) code (Hümmer et
al., 2010; Biehl et al., 2018a), we change technology to self-consistently the
electromagnetic cascade in-source. We use the AM3 (Astrophysical Modeling with
Multiple Messengers) software (Gao et al., 2017), which has been successfully
used for lepto-hadronic models of Active Galactic Nuclei blazars (see, e.g.,
Gao et al., 2019; Rodrigues et al., 2019, 2021) and Gamma-Ray Bursts (Rudolph
et al., 2022, 2023), and apply it for the first time to TDEs. We address
several key questions regarding EM cascade emissions in TDEs, e.g., what is
the spectral shape of EM cascade emissions? What kind of observational
features can be attributed to EM cascades? At what timescales and in which
energy ranges can we expect to observe these EM cascades? What insights can be
derived from X-ray and $\gamma$-ray observations, combined with neutrino
counterparts, to constrain the parameter space of the model? For most of the
results, we focus on the cases M-IR and M-OUV, where the thermal IR and OUV
photons, respectively, dominate the neutrino production. The M-X case, for
which X-ray targets dominate, will only be included in some of the results as
the qualitative conclusions are similar to M-IR.

The paper is organized as follows: in Sec. 2 we provide a brief description of
the TDE models, particle interactions and numerical methods. EM cascade
spectra and light curves for each TDE for the cases M-IR and M-OUV will be
presented in Sec. 3. In Sec. 4, we show the $\gamma$-ray constraints on model
parameters and predicted neutrino numbers for IceCube. We discuss our results
and summarize our conclusions in Sec. 5 & 6.

##  2 Model description

In this section, we describe the TDE model in this work, and the numerical
approaches employed for computing neutrino and EM cascade emissions.

###  2.1 Overview

We focus on the description of a spherical radiation zone of radius $R$
leading to quasi-isotropic emission of neutrinos and photons, following Winter
& Lunardini (2023). The acceleration is assumed to occur inside this region at
a radius $R_{\rm acc}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}R$.111This relationship applies to models
M-OUV and M-IR (see below), whereas $R_{\rm acc}\sim R$ for model M-X. We do
not simulate the accelerator explicitly, we instead parameterize the
acceleration zone by the maximal proton energy $E_{p,\mathrm{max}}$ and the
dissipation efficiency $\varepsilon_{\rm diss}$, which describes the
conversion efficiency from the mass accretion rate into non-thermal protons.
Because of the quasi-isotropic emission, large values of $\varepsilon_{\rm
diss}\simeq 0.05-0.2$ are required to obtain reasonable neutrino event rates;
we will use a value on the upper end of this range to saturate the EM cascade
bounds. Note that its precise value is not so relevant for the conclusions of
this work, as it scales the neutrino and EM emissions in the same way. The
accelerated protons are assumed to be magnetically confined in magnetic fields
of strength $B$ as long as the Larmor radius is smaller than the size of the
region, but they can diffuse out in the Bohm limit.

Protons in the radiation zone may interact with IR, OUV, and X-ray radiation
parameterized as thermal spectra with temperatures $T_{\mathrm{IR}}$,
$T_{\mathrm{OUV}}$, and $T_{\mathrm{X}}$, respectively; the higher the target
photon energy is, the lower is the required $E_{p,\mathrm{max}}$.
Consequently, the models are labeled M-X (lowest $E_{p,\mathrm{max}}$,
interactions with X-ray target only), M-OUV (intermediate
$E_{p,\mathrm{max}}$, interactions dominated by OUV target), and M-IR (highest
$E_{p,\mathrm{max}}$, interactions with IR target possible). Our interaction
targets are based solely on observations, which means that the do not add any
other theoretically motivated ingredients which have not been directly
observed. Our models are fully time-dependent, taking into account the time-
dependence of the proton injections assumed to follow the OUV light curves,
and the time-dependence of the target photons. This is illustrated in Fig. 1,
where the time-dependent proton injection lumimosities and the time evolutions
of the photon targets are shown for AT2019dsg and AT2019fdr; further model
details are listed in Tab. 1. It is especially noteworthy that the IR photons
are assumed to be be produced by the dust as so-called “dust echoes”, where
dust heated by X-ray and OUV radiation may re-emit IR photons (van Velzen et
al., 2016, 2021b; Reusch et al., 2022). We further on use the dust echo IR
light curves obtained by Winter & Lunardini (2023).

We focus on two scenarios, M-IR and M-OUV first, as they are representative
for our results. However, we will later consider extended radiation zones
(e.g., located close to the dust radius) and compact radiation zones (e.g.,
located close to the OUV-BB), respectively, as we will use $R$ as a free
parameter in Sec. 4. Further model details are given in the following
subsections.

###  2.2 Proton injection

Table 1: Observational and TDE modeling parameters for AT2019dsg and
AT2019fdr. In all scenarios, the universal values of energy dissipation
efficiency $\varepsilon_{\rm diss}=0.2$ and magnetic field strength
$B=0.1~{}\rm G$ are used.

| AT2019dsga | AT2019fdrb
---|---|---
|  $z=0.051,~{}M=5\times 10^{6}M_{\odot}$, $t_{\rm dyn}=670~{}\rm d$ | $z=0.267,~{}M=1.3\times 10^{7}M_{\odot},~{}t_{\rm dyn}=1730~{}\rm d$
$k_{B}T_{\rm X,~{}OUV,~{}IR}$ | 72 eV, 3.4 eV, 0.16 eV | 56 eV, 1.2 eV, 0.14 eV
$E_{\nu}$ | 217 TeV (IC191001A) | 82 TeV (IC200530A)
$t_{\nu}-t_{\rm pk}$ | 154 d | 324 d
$N_{\nu}(\rm GFU)^{c}$ | $0.008-0.76$ | $0.007-0.13$
Scenario | M-IR | M-OUV | M-IR | M-OUV
$R$ [cm] | $5.0\times 10^{16}$ | $5.0\times 10^{14}$ | $2.5\times 10^{17}$ | 5.0$\times 10^{15}$
$E_{p,\rm max}$ [GeV] | $5.0\times 10^{9}$ | $1.0\times 10^{8}$ | $5.0\times 10^{9}$ | 1.0$\times 10^{8}$

  * (1)

aAT2019dsg data references: redshift $z$, expected neutrino number via IceCube
GFU searches $N_{\nu}$(GFU), $T_{\rm OUV}$ and $T_{\rm X}$ (Stein et al.,
2021); SMBH mass $M$ (van Velzen et al., 2021b); peak time of OUV light curve
$t_{\rm pk}$ (Stein et al., 2021); Neutrino energy $E_{\nu}$ (IceCube
Collaboration, 2019a); $T_{\rm IR}$ (Winter & Lunardini, 2023).

  * (2)

bAT2019fdr data references: $z$, $t_{\rm pk}$, $N_{\nu}$(GFU), $T_{\rm OUV}$,
$T_{\rm X}$ and $T_{\rm IR}$ (Reusch et al., 2022); $M$ (van Velzen et al.,
2021b); $E_{\nu}$ (IceCube Collaboration, 2019b).

  * (3)

cExpected neutrino number from IceCube gamma-ray follow up (GFU) searches.

Figure 1: Time-dependent evolution of IR, OUV, X-ray, and injected proton
luminosities, shown as the orange, red, blue, and black solid curves
respectively. The dissipation efficiency $\varepsilon_{\rm
diss}=L_{p}/(\dot{M}c^{2})\simeq 0.2$ is used to obtain the proton injection
luminosity $L_{p}$ for AT2019dsg (left panel) and AT2019fdr (right panel). The
horizontal green dashed curves represents the Eddington luminosity $L_{\rm
Edd}$. The OUV luminosity peaks ($t_{\rm pk}$) and neutrino detection times
($t_{\nu}-t_{\rm pk}$) are illustrated as the vertical gray solid and magenta
dotted curves, respectively.

Possible locations for the acceleration zone include the accretion disk
corona, relativistic jets, or the disk winds; correspondingly, $R_{\rm acc}$
may vary significantly. Near the X-ray photosphere or the disk corona, the
size can be $\sim(10-30)R_{s}$ (e.g., Murase et al., 2020; van Velzen et al.,
2021b), where $R_{s}=2GM/c^{2}$ is the Schwarzschild radius of a SMBH with
mass $M$. Inside the jet, particle acceleration is likely to occur at a radius
of a few hundred to 1000 $R_{s}$ (e.g., Dai et al., 2018), depending on the
location of the shocks. Alternatively, if particle acceleration in the wind is
considered, $R_{\rm acc}$ can extend to the boundary of the wind envelope,
e.g., $R_{\rm acc}\sim 10^{16}-10^{17}~{}\rm cm$ (Stein et al., 2021), in
proximity to the dust radius.

In this work, we focus on the description of the radiation zone characterized
by the radius $R\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}R_{\rm acc}$ without explicitly specifying
the particle acceleration site, indicating that the physical conditions of the
radiation and acceleration regions, such as magnetic field strength, may not
necessarily be identical. Note that even though the proton emission from the
accelerator may be anisotropic, the magnetic field in the radiation zone may
isotropize the protons in the calorimetric limit; an example would be an off-
axis jet as accelerator in a large radiation zone (model M-IR). Using an
acceleration rate $t_{\rm acc}^{-1}=\eta_{\rm acc}(c/R_{L,\rm
acc})\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}c/R_{L,\rm acc}$ (Hillas, 1984), the
protons could be accelerated to $10^{8}-10^{9}~{}\rm GeV$ near the accretion
disk corona ($R_{\rm acc}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}10^{14}~{}\rm cm$) with a strong magnetic
field strength or in the extended disk wind ($R_{\rm acc}\sim 10^{17}~{}\rm
cm$), where $R_{L,\rm acc}=E_{p}/(eB_{\rm acc})$ is the Larmor radius of
protons in the acceleration zone characterized by the magnetic field strength
$B_{\rm acc}$, $\eta_{\rm acc}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}1$ is the acceleration efficiency and
$t_{\rm acc}$ is the acceleration time required to energize protons to $E_{p}$
(see discussion in Winter & Lunardini, 2023).

We adopt an energy-differential power-law proton injection spectrum (per unit
volume per unit time) $Q_{p}(E_{p},t)\propto E_{p}^{-2}\exp(-E_{p}/E_{p,\rm
max})$ where $E_{p,\rm max}$ is the maximum energy of injected proton spectrum
and is treated as a free parameter. The spectrum can be normalized via $\int
E_{p}Q_{p}dE_{p}=L_{p}/V$, where $L_{p}$ denotes the time-dependent proton
luminosity and $V=4\pi R^{3}/3$ is the volume of the radiation region.
Typically, the minimum proton energy can be written as $E_{p,\rm
min}=\Gamma_{\rm rel}m_{p}c^{2}$, where $\Gamma_{\rm rel}$ is the relative
Lorentz factor between the acceleration site and the radiation zone. Since the
absence of the specification of the acceleration site, the conservative value
of $E_{p,\rm min}=1~{}\rm GeV$ is used. As depicted by the red curves in Fig.
1, OUV observations reveal that TDEs generically exhibit characteristic light
curves, with an initial rapid increase and a subsequent slower decay ($\propto
t^{-5/3}$) that is consistent with the predicted mass fallback rate. This
behavior suggests the OUV luminosity ($L_{\rm OUV}$) reflects the accretion
history. Henceforth, following the methodology in Winter & Lunardini (2023),
we assume that a fraction of the accreted energy onto the SMBH is dissipated
into protons and the accretion rate aligns with the observed temporal
evolution of OUV light curves. We write down $L_{p}$ in terms of the Eddington
accretion rate $\dot{M}_{\rm Edd}$, which represents the accretion rate
required for a black hole to radiate at the Eddington luminosity $L_{\rm
Edd}\simeq 1.3\times 10^{45}{~{}\rm erg~{}s^{-1}}M/(10^{7}M_{\odot})$, as

$L_{p}=\varepsilon_{\rm diss}\dot{M}c^{2}=\varepsilon_{\rm
diss}\zeta\dot{M}_{\rm Edd}c^{2}\frac{L_{\rm OUV}}{L_{\rm OUV}(t_{\rm pk})}.$
(1)

In this expression, $t_{\rm pk}$ is the time when OUV light curve reaches its
peak $L_{\rm OUV}(t_{\rm pk})$, and $\zeta\equiv\dot{M}(t_{\rm
pk})/\dot{M}_{\rm Edd}$ is the ratio of the peak accretion rate to the
Eddington accretion rate $\dot{M}_{\rm Edd}$. Generally, $\dot{M}_{\rm Edd}$
can be estimated from the radiation efficiency $\eta_{\rm rad}\sim 0.01-0.1$
(e.g., McKinney et al., 2015), that represents the fraction of accreted energy
reprocessed into radiation energy at the level of $L_{\rm Edd}$, via
$\dot{M}_{\rm Edd}=L_{\rm Edd}/(\eta_{\rm rad}c^{2})$. Considering the super-
Eddington accretion nature of TDEs during the peak accretion phase, e.g.,
$\zeta\sim 10-100$ (e.g., Dai et al., 2018), we establish the relationship,
$\dot{M}(t_{\rm pk})=\zeta\dot{M}_{\rm Edd}=(\zeta/\eta_{\rm rad})L_{\rm
Edd}/c^{2}$. We adopt $\zeta/\eta_{\rm rad}=100$ in the following analysis;
this value can be also justified by comparing the time-integrated accretion
mass with the the mass of disrupted stars (Winter & Lunardini, 2023). Under
these parameterizations, the thick black curves in the left and right panels
of Fig. 1 depict the injected proton luminosity for AT2019dsg and AT2019fdr.
We show also the Eddington luminosity as green dashed lines for reference.

###  2.3 Interactions and numerical methods

In the radiation zone, neutrinos can be generated through photomeson
($p\gamma$) and hadronuclear ($pp$) interactions. In the $p\gamma$ process,
the quasi-thermal X-ray, OUV, and IR emissions serve as the target photon
fields. With the time-dependent luminosity $L$ obtained from direct
observations of the IR/OUV/X-ray components, we calculate the target photon
density through

$\int\varepsilon_{\gamma}n(\varepsilon_{\gamma})d\varepsilon_{\gamma}=\frac{L}{4\pi
R^{2}c},$ (2)

where

$n\propto\frac{\varepsilon_{\gamma}^{2}}{\exp\left(\frac{\varepsilon_{\gamma}}{k_{B}T}\right)-1}$
(3)

is the black body spectrum (in the units of $\rm eV^{-1}~{}cm^{-3}$) and $T$
is the characteristic temperature determined by fitting the observed spectra.
For OUV components, the photon density can be straightforwardly obtained from
Eqs. (2) and (3) using the observed $L_{\rm OUV}$ and $T_{\rm OUV}$. In the
context of X-ray photons originated from accretion disks, we assume that the
in-source X-ray luminosity remains relatively stable (e.g., Wen et al., 2020).
The time-dependent fluctuations observed in X-ray observations could be
attributed to absorption and geometrical obscuration. Therefore, we consider
the high-level observed luminosity as the representative X-ray luminosity
within the radiation zone, as described by Winter & Lunardini (2023). Note
that this assumption is only relevant in Sec. 4 for small enough proton
energies. Regarding the IR emissions measured subsequent to the OUV peak, the
dust echo interpretation suggests that echo IR photons are generated by dust
located at a distance of approximately $R_{\rm IR}\sim 10^{16}-10^{17}~{}\rm
cm$. The dust heated by the X-ray and OUV radiation will re-emit IR photons
with a temperature of approximately $k_{B}T_{\rm IR}\simeq 0.16~{}\rm eV$ (van
Velzen et al., 2016, 2021b; Reusch et al., 2022). Additionally, the deflection
of these photons with respect to the line of sight explains the observed time
delay. Reusch et al. (2022) have demonstrated that the IR light curve can be
explained by convolving the OUV luminosity with a normalized box function.
Using this method, Winter & Lunardini (2023) performed a least-square fitting
of the IR observations and obtained the IR light curves for AT2019dsg and
AT2019fdr. Building upon the work of Winter & Lunardini (2023), we incorporate
the IR light curves presented in their study directly into our modeling of the
target photon density. Figure 1 illustrates the time dependence of IR ($L_{\rm
IR}$; orange lines), OUV ($L_{\rm OUV}$; red lines), and X-ray ($L_{X}$; blue
lines) luminosities. The vertical purple dotted lines represent the time of
neutrino detection, $t_{\nu}-t_{\rm pk}$. The time axes in this figure are
defined with the OUV luminosity peak time $t_{\rm pk}$ as the origin. The IR,
OUV, and X-ray target photon densities are determined by Eq. (2) and Eq. (3),
using the respective luminosities. The measured temperatures of the IR, OUV,
and X-ray black body emissions are provided in Tab. 1.

While propagating inside the radiation zone, the injected protons will also
interact with the thermal protons in the sub-relativistic outflow or wind.
Given the wind mass rate $\dot{M}_{w}=\eta_{w}\dot{M}$ (Dai et al., 2018; Yuan
et al., 2021), where $\eta_{w}\equiv\dot{M}_{w}/\dot{M}$ is the fraction of
accreted mass converted to wind, we estimate the average target proton number
density $n_{w}\sim\eta_{w}\dot{M}/(4\pi m_{p}v_{w}R^{2})$ and the $pp$
interaction rate

$t_{pp}^{-1}\simeq
cn_{w}\sigma_{pp}\sim\frac{1}{4\pi}\sigma_{pp}\beta_{w}^{-1}\eta_{w}\frac{\dot{M}}{R^{2}m_{p}}\,.$
(4)

Here $v_{w}=\beta_{w}c\sim\mathcal{O}(0.1c)$ is the wind velocity and
$\sigma_{pp}$ is the cross section for inelastic $pp$ collisions. For super-
Eddington accretors, $\eta_{w}$ typical ranges from $10^{-2}$ to $10^{-1}$ in
the case of the efficient magnetically arrested disks (see e.g., Ohsuga et
al., 2009; Jiang et al., 2019; Yuan et al., 2020). We will demonstrate that
the $p\gamma$ process is significantly more efficient than $pp$, even if an
optimistic value of $\eta_{w}\sim 0.1-0.2$ is used.

We also consider the electromagnetic cascade emission resulting from the
secondary particles generated by the $p\gamma$ and $pp$ processes. Charged
pions produced from $p\gamma$ and $pp$ interactions decay, yielding neutrinos,
electrons, and positrons, while neutral pions directly decay into photons
(denoted as ‘$\pi^{0}$ decay’). The energy deposited in the secondary
electrons and positrons can then be redistributed to the radiation field
through synchrotron (SY) and inverse Compton (IC) radiation in the presence of
magnetic fields. We denote this process as ‘$pp/p\gamma$-SY/IC’. As a
competing channel to the $p\gamma$ process, protons can also lose energy
through Bethe-Heitler (BH) interactions, leading to the production of
electron-positron ($e^{\pm}$) pairs. Additionally, the $\gamma\gamma$
attenuatios between the EM cascade emissions and the abundant thermal IR, OUV
and X-ray photons give rise to $e^{\pm}$ pairs. These $e^{\pm}$ pairs from BH
and $\gamma\gamma$ processes are expected to contribute to the EM cascade
radiation field, represented as ‘BH-SY/IC’ and ‘$\gamma\gamma$-SY/IC’
respectively.

Figure 2: Proton interaction rates at neutrino detection time $t_{\nu}$ for
AT2019dsg, as a function of proton energy (in the observer’s frame). The left
and right panels correspond to the M-IR and M-OUV scenarios. In each figure,
the $p\gamma$, $pp$ with $\eta_{w}=0.2$, Bethe-Heitler (BH), escape, and
dynamical rates are shown as the blue, green, orange, red, and magenta curves.
The horizontal curves depict the free streaming rate $t_{\rm
fs}^{-1}=(R/c)^{-1}$, whereas the vertical gray curves show the maximum
injected proton energies $E_{p,\rm max}$.

To obtain the neutrino and EM cascade spectra, we numerically solve the
coupled time-dependent transport equations for all relevant particle species

$\frac{\partial n_{i}}{\partial t}=Q_{i,\rm ext}+\sum_{k}Q_{k\rightarrow
i}-\frac{\partial}{\partial E}(\dot{E}\cdot n_{i})-\frac{n_{i}}{t_{i,\rm
esc}}$ (5)

in a self-consistent way using AM3 software (e.g., Gao et al., 2017). In Eq.
(5), $i$ is the label of the particle, $n_{i}\equiv dN_{i}/dVdE$ is the in-
source particle density differential in energy and volume, $N_{i}$ is the
total number, $Q_{i,\rm ext}$ is the external particle injection rate, e.g.,
$Q_{p,\rm ext}=Q_{p}(E_{p},t)$, $\sum_{k}Q_{k\rightarrow i}$ represents the
in-source particle injection from other radiation processes, $\dot{E}=dE/dt$
is the total energy loss rate, and $t_{i,\rm esc}$ is the particle escape
time. Since we focus on the EM cascade emission induced by accelerated protons
and do not consider the primary leptonic loading in the acceleration zone, the
external electron/positron injection rate is set to zero, e.g., $Q_{e,\rm
ext}=0$. The influence of primary electron injections will be discussed in
Sec. 5.

Regarding the particle escape times, the photons, neutrinos, and neutrons
produced in $p\gamma$ processes escape the radiation zone with the free-
streaming time $t_{\rm esc}=t_{\rm fs}=R/c$, whereas the escape rate of
charged particles is determined by the diffusion rate. We use as proton escape
rate $t_{p,\rm esc}^{-1}=\min[t_{\rm fs}^{-1},t_{p,\rm diff}^{-1}]$ (Baerwald
et al., 2013), where $t_{p,\rm diff}^{-1}\equiv D/R^{2}$ is the diffusion
rate, $D\sim R_{L}c$ is the diffusion coefficient in Bohm limit, and
$R_{L}=E_{p}/(eB)$ is the Larmor radius of protons with energy $E_{p}$ in the
magnetic field $B$. For simplicity, we choose the magnetic field strength
$B=0.1$ G as the fiducial value for both AT2019dsg and AT2019fdr as in Winter
& Lunardini (2023). This choice of magnetic field is consistent with the case
of outflows from AGNs and is sufficiently strong to confine protons with
energies $E_{p}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}1$ PeV in the radiation zone. We will
investigate the impacts on neutrino and EM cascade emissions by varying $B$
from 0.1 G to 1.0 G in Sec. 4.

We have checked the consistency of our results with the earlier computations
in Winter & Lunardini (2023) using NeuCosmA. While our results overall agree
very well, we find small differences e.g. in the shapes of the neutrino
spectra at the level expected from different softwares using different solver
and particle interaction implementations (see e.g. Cerruti et al., 2022); for
example, these come from the different implementations of photohadronic
interactions: while AM3 uses Hümmer et al. (2010), NeuCosmA is based on Biehl
et al. (2018a).

###  2.4 M-IR and M-OUV scenarios

Based on the formulation described before, the physical properties of the
radiation zone can be characterized by two parameters: $R$ and $E_{p,\rm
max}$, assuming a fixed magnetic field strength of $B=0.1$ G. We investigate
the EM cascade emissions of two scenarios, M-IR and M-OUV, for AT2019dsg and
AT2019fdr. The corresponding parameter sets are listed in Tab. 1.

_M-IR scenario._ As for $p\gamma$ process, the proton threshold energy is

$E_{p,\rm th}=\frac{m_{\pi}(2m_{p}+m_{\pi})c^{4}}{4E_{\gamma}}\simeq 7.1\times
10^{8}~{}{\rm GeV}\left(\frac{E_{\gamma}}{0.1~{}\rm eV}\right)^{-1},$ (6)

where $E_{\gamma}$ is target photon energy and $m_{\pi}\simeq 140~{}\rm
MeV/c^{2}$ is the pion mass. We choose $R=5\times 10^{16}~{}\rm
cm\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}R_{\rm IR}$ and $E_{p,\rm max}=5.0\times
10^{9}~{}{\rm GeV}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}E_{p,\rm th}$ for the M-IR scenario of
AT2019dsg to ensure that $p\gamma$ interactions with thermal IR photons with
energy $E_{\gamma,\rm IR}=k_{B}T_{\rm IR}\simeq 0.16~{}\rm eV$ dominate the
neutrino production. The left panel of Fig. 2 shows the proton interaction
rates at neutrino detection time $t_{\nu}$ due to $p\gamma$
($t_{p\gamma}^{-1}$, blue curve), $pp$ ($t_{pp}^{-1}$ with $\eta_{w}=0.2$,
green line) and BH ($t_{\rm BH}^{-1}$, orange curve) interactions. The
horizontal black line is the particle free-streaming rate $t_{\rm fs}^{-1}$,
whereas the magenta horizontal line depicts the reciprocal of the dynamic
time, e.g., $t_{\rm dyn}^{-1}$. We define the dynamic time ($t_{\rm dyn}$) as
the duration of the super-Eddington phase, e.g., $\dot{M}c^{2}>L_{\rm Edd}$.
The energy-dependent escape rate is illustrated as the red line. The vertical
gray line and the shaded region describe the maximum proton energy $E_{p,\rm
max}$.

From this figure, we can deduce that the system is $p\gamma$ optically thin,
indicated by $\tau_{p\gamma}^{\rm fs}=t_{p\gamma}^{-1}/t_{\rm fs}^{-1}<1$, but
nearly calorimetric at ultra-high energies, as $\tau_{p\gamma}^{\rm
cal}=t_{p\gamma}^{-1}/t_{\rm esc}^{-1}\sim 1$ around $10^{9}$ GeV. This means
that in the M-IR scenario, the TDE can be a promising very-high-energy
neutrino emitter, but a time delay in the scale of $t_{p\gamma}$ of the
$p\gamma/pp-$SY/IC cascade channel is typically expected. Furthermore, the
$p\gamma$ process dominates the neutrino production as $t_{p\gamma}^{-1}\gg
t_{pp}^{-1}$.

_M-OUV scenario._ According to Eq. (6), the proton threshold energy decreases
with increasing energy of the target photon. We therefore consider a compact
radiation zone with $R=5\times 10^{14}~{}\rm cm\ll R_{\rm IR}$ and a
relatively lower maximum proton energy $E_{p,\rm max}=1.0\times 10^{8}~{}\rm
GeV$ for the M-OUV scenario of AT2019dsg. In this case, $p\gamma$ interactions
with OUV target photons dominate the neutrino flux. As shown in the right
panel of Fig. 2, the IR bump in the $p\gamma$ interaction rate (blue curve) is
suppressed due to the higher density of OUV photons ($n_{\rm OUV}\propto
R^{-2}$) compared to IR photons ($n_{\rm IR}\propto R_{\rm IR}^{-2}$).
Notably, in this scenario, the source can be $p\gamma$ optically thick from
$t_{\rm pk}$ to $t_{\nu}$, as indicated by the right panel of Fig. 2, where
$\tau_{p\gamma}^{\rm fs}=t_{p\gamma}^{-1}/t_{\rm fs}^{-1}>1$ is satisfied at
$t_{\nu}$. The $p\gamma$ interactions can be much faster at $t_{\rm pk}$ where
the OUV luminosity is approximately one order of magnitude higher (see the red
curve in the left panel of Fig. 1).

We use AT2019dsg as an example to illustrate the specifications of the M-IR
and M-OUV scenarios. Similar conclusions can be drawn for AT2019fdr using the
parameters listed in Tab. 1.

Figure 3: Left column: Electromagnetic cascade and neutrino spectra at the
time of the neutrino observation $t_{\nu}$. The black and magenta curves show
the spectra of EM cascade and single-flavor neutrino emissions. The black body
IR, OUV and X-ray spectra are depicted as the thin cyan dotted curves. The
thin blue dashed, orange dash-dotted, green dotted, and red curves represent
the components contributing to the EM cascade. The blue, red, and green shaded
areas depict the energy ranges of XRT, _Fermi_ -LAT, and HAWC observations,
whereas the gray-shaded areas indicate where $\gamma\gamma$ attenuation with
CMB and EBL prohibits the observation of the emitted gamma rays. Right column:
Neutrino, $\gamma$-ray and X-ray light curves. The the magenta, green, orange,
and blue solid curves illustrate the single-flavor neutrino ($L_{\nu}/4\pi
d_{L}^{2}$), very-high-energy (VHE) $\gamma$-ray (0.3 - 100 TeV), $\gamma$-ray
(0.1 - 800 GeV), and X-ray light curves, respectively. The flux up limits from
HAWC, _Fermi_ , and XRT/XMM/NICER or eROSITA (Stein et al., 2021; Cannizzaro
et al., 2021; Predehl et al., 2021; Reusch et al., 2022) are shown as the
green, orange, and blue points, respectively. Furthermore, the neutrino
detection times are marked as the vertical magenta dotted curves. The OUV
fluxes are also plotted as the gray dashed curves for reference. Light blue
shaded areas and curves correspond to XRT (0.3 - 10 keV), whereas dark blue
shaded areas and curves to eROSITA (0.3 - 2 keV). The upper and lower rows
correspond to AT2019dsg and AT2019fdr, respectively. All panels are shown for
model M-IR in the observer’s frame.

##  3 Neutrino and EM cascade emissions

In this section, we present the spectral energy distributions (SEDs), i.e.,
$E_{i}^{2}\mathcal{F}_{i}=t_{\rm fs}^{-1}E_{i}^{2}dN_{i}/dE_{i}/(4\pi
d_{L}^{2})$, of neutrino and electromagnetic cascade emissions for M-IR (Sec.
3.1) and M-OUV (Sec. 3.2) scenarios of AT2019dsg and AT2019fdr, where $d_{L}$
is the luminosity distance to the TDE and the red shift corrections to fluxes
and energies are not explicitly shown in this equation for simplicity. To
establish a connection between our models and observations, we provide also
light curves in the relevant X-ray and gamma-ray energy bands, and compare our
predictions with measured data points or upper limits.

###  3.1 M-IR: extended radiation zone ($R\simeq R_{\rm IR}$)

Using the methodology described in Section 2 and the parameter values for M-IR
scenario provided in Tab. 1, we present in Fig. 3 the SEDs and light curves
for AT2019dsg (upper row) and AT2019fdr (lower row) in the observer’s frame.
The left column shows the SEDs at neutrino detection time $t_{\nu}$, depicting
the thermal IR, OUV and X-ray spectra (in blue dotted curves denoted as ‘BB’),
as well as the single-flavor neutrino (e.g., $\nu_{\mu}$, in magenta curves)
and EM cascade spectra spanning from MeV to 100 PeV. The black curves
represent the cascade spectra without taking into account the $\gamma\gamma$
attenuation with extragalactic background light (EBL) and cosmic microwave
background (CMB). The different components of the electromagnetic cascades,
including ‘$pp/p\gamma-$SY/IC’, ‘$\gamma\gamma-$SY/IC’, ‘BH-SY/IC’, and
‘$\pi^{0}$ decay’, are depicted as the blue dashed, orange dash-dotted, green
dotted, and red solid curves, respectively. The gray shaded areas in the
figure highlight the energy ranges of cosmic $\gamma\gamma$ attenuation with
EBL and CMB. One notable feature of the EM cascade spectra is the pronounced
dips in the energy range $\sim 1~{}\rm GeV-1~{}PeV$. These dips can be
attributed to the in-source $\gamma\gamma$ absorption with thermal photons.
Considering the characteristic energy of thermal photons $E_{\rm BB}$, the
$\gamma\gamma$ attenuation occurs at
$E_{\gamma\gamma}\sim(m_{e}c^{2})^{2}/E_{\rm BB}\simeq 2.6{~{}\rm
GeV~{}}(E_{\rm BB}/100~{}{\rm eV})^{-1}$ in the laboratory frame. For
instance, X-ray photons with energy of $k_{B}T_{X}=72$ eV will lead to the
first absorption dip at $\sim 4$ GeV in the black curve of the upper left
panel.

For AT2019dsg, the blue, red, and green shaded areas in the upper left panel
illustrate the energy ranges for the observations or limits using
XRT/XMM/NICER (0.3-10 keV), _Fermi_ LAT (0.1-800 GeV) and High-Altitude Water
Cherenkov Observatory (HAWC, 0.3-100 TeV) (see Ayala & HAWC Collaboration,
2019; Stein et al., 2021; Cannizzaro et al., 2021; Reusch et al., 2022;
Predehl et al., 2021, and references therein). The corresponding light curves
are shown in the upper right panel. The orange and green data points are the
energy-integrated upper limits obtained by _Fermi_ and HAWC, whereas the
XRT/XMM/NICER observations are plotted as blue points. The _Fermi_
$\gamma$-ray up limits are obtained in three time slots, $t-t_{\rm
pk}\sim-20-100{~{}\rm d},~{}{100-200~{}\rm d}$ and $200-500~{}\rm d$. The
upper limit in the first time slot is more stringent according to the upper
right panel of Fig. 3 and will be used in Sec. 4 for the constraints on the
$R-E_{p,\rm max}$ plane. The multiwavelength light curves indicate that the
$\gamma$-ray upper limits are satisfied in the M-IR scenario, although the
limits are not far way from the predictions. This also implies that much
higher predicted neutrino fluxes would contradict the $\gamma$-ray limits.

For reference purposes, we include the thermal OUV light curve, e.g., $L_{\rm
OUV}/4\pi d_{L}^{2}$, as the black dashed curve. The peaks of EM cascade light
curves appear roughly $30-50$ days after the OUV peak, which is commonly
referred to as the time delay. In the M-IR scenario, where the radiation zone
is calorimetric but optically thin to $p\gamma$ interactions
($\tau_{p\gamma}^{\rm fs}<1$), the time delay observed in the EM cascade
emissions, primarily initiated via $p\gamma$ processes, can be interpreted as
the $p\gamma$ time scale, e.g., $t_{p\gamma}\sim 30-50$ days. This reflects
the time-dependent nature of EM cascade as a redistribution of energy via a
series of successive radiation processes. Hence, if future observed TDEs
exhibit time delays in their X-ray and $\gamma$-ray light curves, one can
employ $p\gamma$ optically-thin models to interpret and understand such
temporal signatures.

The rapid (exponential) decay of observed X-ray light curve in the time
interval $0\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}t-t_{\rm pk}\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}50~{}\rm d$
(illustrated as the blue dashed line in the upper right panel of Fig. 3) might
be caused by the cooling of the accretion disk or by dust obscuration.
Interestingly, there is an additional X-ray (Swift-XRT) data point around
$t-t_{\rm pk}=100$ d which was identified in Cannizzaro et al. (2021) and does
not fit this rapid decay picture. Notably, the X-ray cascade emission can
describe that data point, which means that Swift may have actually seen an
electromagnetic cascade signature there. One may speculate that if the
populations of neutrino-emitting TDEs and X-ray bright TDEs have overlap, the
often complicated and puzzling behavior of X-ray observations may in some
cases be explained by the contribution of the EM cascade, which can be
substantially delayed with respect to the BB peak (see e.g. example for
AT2019fdr below). The identification of a spectral hardening in the X-ray
range (see upper left panel) may provide future evidence.

Our calculation also provide the neutrino spectra and light curves, as shown
in the magenta curves in Fig. 3. Since $p\gamma$ interactions dominate the
neutrino production compared to the $pp$ process, we show only the single-
flavor neutrino spectra obtained from the $p\gamma$ channel in the left
column. The black arrows indicate the most-likely energies of the coincident
neutrinos. From the upper right panel, we observe that the energy-integrated
neutrino flux, $L_{\nu}/4\pi d_{L}^{2}$, aligns well with the OUV light curve.
The explanation is that, in the M-IR scenario of AT2019dsg, the radiation zone
is nearly calorimetric to UHECRs, e.g., $\tau_{p\gamma}^{\rm
cal}=t_{p\gamma}^{-1}/t_{p,\rm esc}^{-1}\sim 1$ and the $p\gamma$ interaction
rate is close to the free streaming rate (see the left panel of Fig. 2). This
implies that the $p\gamma$ interactions are relatively efficient and fast (see
Winter & Lunardini, 2023, for a more detailed discussion).

Figure 4: Electromagnetic cascade spectra and neutrino spectra at $t_{\nu}$
(left column), and neutrino, $\gamma$-ray and X-ray light curves (right
column), for the M-OUV model. The curves, areas, and data points have the same
meaning as in Fig. 3.

Correspondingly, the results for AT2019fdr are displayed in the lower row of
Fig. 3. The areas, data points, and curves have the same interpretation as in
the AT2019dsg cases. However, in this case, the deep blue areas, data points,
and curves specifically represent the observations and predictions in the
energy range of 0.3$-$2 keV, associated with the eROSITA upper limits. The
upper limit on $\gamma$-ray emissions is only measured within the time
interval $-20~{}\rm{d}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}t-t_{\rm pk}\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}300~{}\rm{d}$.
Additionally, X-ray follow-up observations using Swift XRT provide an upper
bound around the time of the OUV peak, while the later time constraints
(represented by the green points) are measured by eROSITA. The M-IR scenario
of AT2019fdr is consistent with the X-ray and $\gamma$-ray observational
limits. However, it is challenging for the EM cascade models to explain the
late-time X-ray peak observed at $t-t_{\rm pk}\simeq 600$ days by eROSITA.
This particular signature may be produced by other mechanisms or regions that
are distinct from the neutrino-emitting zones discussed here, especially
considering that it appears roughly 300 days after the detection of neutrinos.

The time delays in X-ray and $\gamma$-ray light curves are also observed and
can be attributed to the $p\gamma$ interaction times. Since a larger radius of
the radiation zone, e.g., $R=2.5\times 10^{17}$ cm, is used for AT2019fdr,
compared to that of AT2019dsg ($R=5\times 10^{16}$ cm), a prolonged time delay
of approximately 200 days is expected due to the decreasing $p\gamma$
interaction rate. Unlike the case of AT2019dsg, the neutrino light curve of
AT2019fdr no longer simply follows the time dependence of OUV luminosity. The
reason is that the neutrino production is determined not only by the $p\gamma$
interaction rate but also the in-source proton distributions. For AT2019fdr,
the protons are more strongly confined within the radiation zone in comparison
to AT2019dsg. As a consequence, there is a buildup or “piling-up” of protons,
which enhances the late-time neutrino productions and lead to the time
delay222Here, we only provide a qualitative explanation for the time
dependence of neutrino emission, as this paper primarily focuses on the EM
cascades. A more quantitative discussion for the neutrino light curves of
AT2019dsg and AT2019fdr can be found in Winter & Lunardini (2023)..

In summary, in the M-IR scenario where the radiation zone can extend to the
dust radius, the $p\gamma$ optical depth is typically optically thin. This
leads to time delays in the EM cascade light curves. One advantage of this
scenario is that the predicted X-ray and $\gamma$-ray emissions obey well with
the observational upper limits, as shown in the right column of Fig. 3.
Remarkably, the M-IR scenario of AT2019dsg is capable of explaining the X-ray
observations during the time interval $50~{}{\rm d}\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}t-t_{\rm
pk}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}100~{}\rm d$.

###  3.2 M-OUV: compact radiation zone ($R\ll R_{\rm IR}$)

Here we consider more compact radiation zones in the M-OUV scenarios and
discuss the EM cascade signatures in the $p\gamma$ optically thick cases, for
AT2019dsg and AT2019fdr. Similar to the Fig. 3, the SEDs and light curves for
the M-OUV scenarios are demonstrated in Fig. 4. The parameters used here are
summarized in Tab. 1. We find that the neutrino peak fluxes in the M-OUV
scenarios for both AT2019dsg and AT2019fdr are significantly higher compared
to the M-IR scenarios. This enhancement is attributed to the efficient
$p\gamma$ interactions in the M-OUV scenarios, as indicated in the right panel
of Fig. 2, facilitated by the much denser target photons. As a result, the
neutrino fluxes in the M-OUV scenarios are enhanced by a factor of $2-3$
compared to those in the M-IR scenarios. Additionally, the peaks of the
neutrino spectra closely align with the energies of the coincident neutrinos.
This, in combination with the high rate of efficient $p\gamma$ interactions,
makes TDEs in M-OUV scenarios as the promising neutrino emitters. From the
right column, we find that the neutrino light curves exhibit the same shape as
the OUV luminosities ($L_{\rm OUV}/4\pi d_{L}^{2}$ is not shown for
simplicity). This can be straightforwardly interpreted by the optically thick
nature of the radiation zones, i.e., the protons interact very quickly before
they can escape even compared to the size of the region, where the condition
$\tau_{p\gamma}^{\rm fs}=t_{p\gamma}^{-1}/t_{p,\rm fs}^{-1}>1$ is satisfied,
as depicted in the right panel of Fig. 2. Consequently, no significant time
delays are observed.

The EM cascade fluxes, similar to the neutrinos, are enhanced as well,
particularly at the peak time $t_{\rm pk}$, and there are no significant time
delays with respect to the OUV curve. For AT2019dsg, the upper right panel of
Fig. 4 demonstrates that the _Fermi_ $\gamma$-ray upper limits are violated in
the M-OUV scenarios with the radius of $R=5\times 10^{14}~{}\rm cm$. It is
worth noting that when the radius is increased to $10^{15}-10^{16}~{}\rm cm$,
the peak flux shows a significant drop, and the upper limits from _Fermi_
observations are no longer violated. This is observed in the M-OUV scenario of
AT2019fdr, as shown in the lower right panel. In this discussion, we provide a
intuitive conclusion regarding the _Fermi_ constraints. A more detailed
analysis of the implications of these limits on the $R-E_{p,\rm max}$ plane
will be presented in Sec. 4.

One common feature of the EM cascade SEDs in the M-IR and M-OUV scenarios is
that the ‘$pp/p\gamma-$SY/IC’ (blue dashed curves) and ‘$\gamma\gamma-$SY/IC’
(orange dashed-dotted curves) components in the left columns of Figs. 3 and 4
give rise to the low-energy ($\sim 0.1-10$ MeV) and high-energy ($\sim 0.1-10$
GeV) humps, respectively. Noting that in $p\gamma$ channel around $5\%$ of the
proton energy is inherited by secondary leptons, we estimate the
characteristic Lorentz factor of electrons and positrons to be
$\gamma_{e,p\gamma}\simeq(E_{p,\rm cal}/20)/(m_{e}c^{2})\sim 10^{8}(E_{p,\rm
cal}/(10^{6}\,\mathrm{GeV}))$, where $E_{p,\rm cal}\simeq 10^{5}-10^{6}$ GeV
is the proton energy at which the radiation zone becomes calorimetric, e.g.,
$t_{p\gamma}^{-1}\sim t_{\rm esc}^{-1}$. We then infer the peak energy of
synchrotron radiation from $p\gamma$ electrons and positrons

$E_{p\gamma,\rm
sy}=\frac{3}{4\pi}h\gamma_{e,p\gamma}^{2}\frac{eB}{m_{e}c}\simeq 17{~{}\rm
MeV}\left(\frac{E_{p,\rm cal}}{10^{6}~{}\rm GeV}\right)^{2},$ (7)

where the magnetic field strength $B=0.1$ G is used. The corresponding inverse
Compton components are typically invisible due to $\gamma\gamma$ absorption.
Regarding the ‘$\gamma\gamma-$SY/IC’ channel, the interpretation is more
complicated since the $\gamma\gamma$ cross section exhibits slow decay beyond
the threshold, suggesting a wide energy range of absorption. From the EM
cascade spectra, we observe that the strongest absorption occurs close to the
peak energies of $\gamma$-rays from $\pi_{0}$ decays, as predicted by
$E_{\gamma\gamma,\rm max}\sim E_{p,\rm th}/10$, where $E_{p,\rm th}$ is
defined in Eq. (6) and we consider that approximately 10% of the proton energy
is transferred to $\gamma$-rays via the $\pi_{0}\rightarrow\gamma\gamma$
channel. The peak energy of the resulting synchrotron emission, produced by
electrons of $E_{\gamma\gamma,\rm max}/2$, is given by $E_{\gamma\gamma,\rm
sy}\sim\mathcal{O}(1{~{}\rm TeV})(E_{\gamma}/0.1\rm eV)^{-2}$, where
$E_{\gamma}$ is the energy of target photons. In this case, the synchrotron
emission in the ‘$\gamma\gamma-$SY/IC’ channel suffers from attenuation before
reaching its maximum $E_{\gamma\gamma,\rm sy}$, which explains the sharp
spikes appearing in the EM cascade spectra. This semi-theoretical estimation
justifies the numerical results shown in Figs. 3 and 4.

Figure 5: Number of expected single-flavor neutrino events as a function of
$R$ and $E_{p,\rm max}$ as given by the color scales; we use the IceCube
Gamma-Ray Follow-up (GFU) effective areas to compute the neutrino numbers. The
panels correspond to AT2019dsg (upper row) and AT2019fdr (lower row), and
magnetic field strengths $B=0.1$ G (left column) and $1$ G are (right column).
The _Fermi_ $\gamma$-ray upper limits on the $R-E_{p,\rm max}$ plane are shown
as red curves (the parameters in the upper left corners are excluded), and the
regions where the optical thickness $\tau_{p\gamma}^{\rm fs}(E_{p\gamma,\rm
max})=1$ are shown as blue dashed curves (the parameters in the upper left
corners are optically thick); here the $p\gamma$ optical depth is evaluated at
$t_{\nu}$ and $E_{p\gamma,\rm max}$ and the definition of $E_{p\gamma,\rm
max}$ can be found in the main text. Moreover, the models M-X, M-IR and M-OUV
models are marked by stars.

##  4 $\gamma$-ray constraints on $E_{p,\rm max}$ and $R$

In the previous sections, we examined EM cascade and neutrino emissions in two
special scenarios, M-IR and M-OUV, for AT2019dsg and AT2019fdr. It is useful
to investigate the model-predicted neutrino numbers as a function of $R$ and
$E_{p,\rm max}$ and explore the potential constraints on source parameters
using _Fermi_ $\gamma$-ray up limits.

We choose $R\in[5\times 10^{14}~{}{\rm cm},10^{18}~{}{\rm cm}]$
($[10^{15}~{}{\rm cm},10^{18}~{}{\rm cm}]$ for AT2019fdr) and $E_{p,\rm
max}\in[10^{6}~{}{\rm GeV},10^{10}~{}{\rm GeV}]$ as the fiducial ranges to
allow for a continuous interpolation among the different models. From the left
panels of Figs. 3 and 4, we notice that the _Fermi_ $\gamma$-ray constraints
are more stringent compared to the X-ray and HAWC $\gamma$-ray observations.
Hence, we use the _Fermi_ constraints in the energy range $0.1-800$ GeV to
study the permissible parameter sets that are compatible with observations.

Our main results for parameter constraints are illustrated in Fig. 5. The
upper and lower rows correspond to the TDEs AT2019dsg and AT2019fdr, whereas
the magnetic fields of $B=0.1$ G and $1.0$ G are used in the left and right
columns respectively to demonstrate the impact of different magnetic field
strengths. The M-OUV and M-IR scenarios are indicated by blue and red stars,
respectively. We also include the parameters (e.g., $R=5\times 10^{15}{~{}\rm
cm},~{}E_{p,\rm max}=5\times 10^{6}{~{}\rm GeV}$ and $B=1.0$ G) for the M-X
scenarios discussed in Winter & Lunardini (2023) as the magenta stars. The red
contours represent the _Fermi_ upper limits (denoted as ‘ _Fermi_ UL’). For a
TDE radiation zone described by the parameter sets above the red contours, the
EM cascade emission will overshoot the _Fermi_ upper limits. This conclusion
is consistent with the results in Sec. 3, which show that the compact
radiation zones violate the upper limits.

To understand the compactness of the radiation zone, we present the critical
conditions for the source to be $p\gamma$ optically thick in the dashed blue
contours as well. These conditions are obtained through $\tau_{p\gamma}^{\rm
fs}(E_{p\gamma,\rm max})=1$, where $E_{p\gamma,\rm max}<E_{p,\rm max}$ is the
proton energy at which $t_{p\gamma}^{-1}$ reaches its maximum. We find that
there exists a critical radius beyond which the system can no longer be
$p\gamma$ optically thick even for UHE protons. The reason is that
$t_{p\gamma}^{-1}$ is limited by the maximum target photon density and the
optical depth decreases with respect to $R$, e.g., $\tau_{p\gamma}^{\rm
fs}=t_{p\gamma}^{-1}/t_{\rm fs}^{-1}\propto R^{-1}$, as
$t_{p\gamma}^{-1}\propto R^{-2}$ and $t_{\rm fs}^{-1}\propto R^{-1}$. Above
the blue dashed lines, the EM cascade emissions do not exhibit significant
time delays due to the faster secondary radiation processes compared to photon
free streaming. However, below the blue dashed line, a time delay on the scale
of $t_{p\gamma}$ can be expected.

To discuss the implications of _Fermi_ $\gamma$-ray constraints on neutrino
detection, we show the model-predicted neutrino numbers as a function of $R$
and $E_{p,\rm max}$ in the form of meshed color maps with the red contours in
Fig. 5. The IceCube gamma-ray follow-up (GFU)333The GFU neutrino numbers,
estimated using IceCube GFU effective areas (IceCube Collaboration et al.,
2016), are more suitable when comparing the model predictions to actual
follow-up observations, whereas point source (PS) neutrino numbers are
typically used for independent point-source analyses. neutrino numbers can be
estimated via

$N_{\nu,\mu}({\rm GFU})=\int dE_{\nu}A_{\rm
eff}(E_{\nu})F_{\nu}(E_{\nu},t_{\nu})$ (8)

where $F_{\nu}(E_{\nu},t_{\nu})=\int^{t_{\nu}}dt\mathcal{F}_{\nu}(E_{\nu},t)$
is the cumulative neutrino fluence until $t_{\nu}$ and $A_{\rm eff}$ is the
IceCube GFU effective area. The black contours in Fig. 5 indicate the single-
flavor neutrino numbers $N_{\nu,\mu}$(GFU)$=1,~{}0.1$, and 0.01. From this
figure, we find that to avoid overshooting the _Fermi_ $\gamma$-ray up limits,
the model-predicted neutrino numbers detected by IceCube are constrained to be
$N_{\nu,\mu}({\rm GFU})\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}0.1$. This conclusion applies to both
cases with $B=0.1$ G (left column) and 1.0 G (right column). Specifically, our
findings indicate that the M-X scenarios share similarities with the M-IR
cases. Both scenarios are situated in the $p\gamma$ optically thin regime, and
they both predict a similar neutrino event rate and satisfy the upper limits
set by the _Fermi_ observations.

When comparing the left column with the right column, it is worth noting that
a stronger magnetic field would result in a higher predicted neutrino number,
assuming fixed source parameters. Additionally, a stronger magnetic field
would also increase the likelihood of exceeding the $\gamma$-ray limits due to
enhanced confinement of charged particles within the radiation zones. This
enhanced confinement would lead to more efficient $p\gamma$ interactions and
subsequent EM cascades.

From a multi-messenger perspective, we observe a correlation among three
regions: high optical thickness for $p\gamma$ interactions, high rates of
neutrino events, and high fluxes of $\gamma$-rays. Specifically, high neutrino
event rates indicate the presence of intense EM cascades, which are
facilitated by efficient and rapid $p\gamma$ interactions. The fact that time
delays in neutrino signals with respect to OUV peaks are observed while no
$\gamma$-ray induced by EM cascades has been seen implies that the neutrino
production cannot be too efficient, e.g., the predicted upper limit for
neutrino event rates from _Fermi_ constraints ranges between 0.01 and 0.1
events per TDE, as one can read off from the Fig. 5, and the system is
probably optically thin for $p\gamma$ interactions. Future multi-messenger
observations, in particularly incorporating the $\gamma$-rays measurements,
could shed more lights onto the radiation processes and physical conditions in
the radiation zones.

##  5 Discussion

###  5.1 EM cascades from TDEs

As noted before in Sec. 3, the EM cascade emissions exhibit some commonalities
across different scenarios and TDEs. One prominent feature of the SEDs is the
deep dip in the energy range approximately from 10-100 GeV to 1-10 PeV, which
is a result of the in-source $\gamma\gamma$ attenuation with thermal IR, OUV,
and X-ray photons. Therefore, only a part of the whole EM cascade spectra,
ranging from approximately 1-10 keV to 1-10 GeV, can be measured by eROSITA,
Swift XRT, XMM-Newton, NICER, and _Fermi_ LAT in the observer’s frame.
Additionally, the multi-wavelength light curves of the EM cascades shown in
the right columns of Figs. 3 and 4 also display unique signatures with respect
to the IR, OUV and X-ray BB emissions. For instance, in M-IR cases where the
$p\gamma$ interactions are slow, the peaks of EM cascade light curves usually
emerge after the OUV peak, leading to the time delays. More generally, similar
time delays are also expected if the radiation zone has a large radius or the
maximum energy of injected protons is low, as indicated by the regions below
the blue dashed lines in Fig. 5. In this case, X-ray and $\gamma$-ray follow-
ups after the identification of OUV peaks would reveal these delayed
emissions. The fully time-dependent treatment in this paper provides a more
accurate modeling for the neutrino and EM cascade emissions compared to
steady-state approximations, particularly in cases where the interactions are
efficient (calorimetric) but not sufficiently fast ($p\gamma$ optically thin).
On the other hand, if the radiation zone is compact and $E_{p,\rm max}$ is
high (see the regions above the blue dashed lines in Fig. 5), the system can
be $p\gamma$ optically thick and the light curves of EM cascade emissions
reach their peak roughly at the same time with the neutrino and thermal OUV
light curves. In both cases, the X-ray and $\gamma$-ray fluxes can reach the
detectable levels, making the follow-up observations possible.

In addition to the neutrino-coincident TDEs AT2019dsg and AT2019fdr, whose
neutrino counterparts are identified via follow-up searches, AT2019aalc stands
out with the highest dust echo flux among all ZTF transients and is also
likely to be associated with the neutrino event IC 191119A. The SMBH mass and
AT2019aalc is comparable to that of AT2019fdr, while the _Fermi_ limits and
the time difference of neutrino arrival measured relative to the OUV peak are
similar to those of AT2019dsg. However, no significant X-ray limits or
observations for AT2019aalc has been reported. Winter & Lunardini (2023)
systematically studied the neutrino emissions from all three neutrino-
coincident TDEs and confirmed the similarities between them. Despite not shown
explicitly in this paper, we investigated the EM cascade emissions from
AT2019aalc as well and got similar conclusions with AT2019dsg and AT2019fdr.
For example, high neutrino event rates imply strong EM cascades and violation
of _Fermi_ upper limits. For instance, the M-OUV scenario of AT2019aalc would
also overshoot the _Fermi_ upper limits obtained by 207 days of observations
from the discovery of the optical emission. This on the other hand
demonstrates that our model can be widely applied to TDEs.

Figure 6: Comparison of EM cascade emission (black solid curve) with the
leptonic contribution (dashed curves). The shaded areas have the same meaning
as those in the upper left panel of Fig. 3. To obtain the EM cascade spectrum,
the parameters used in M-IR scenario of AT2019dsg are applied (e.g., $B=0.1$
G). From the upper to lower dashed curves, the ratio of injected electron
luminosity to the proton luminosity, $K_{e/p}\equiv L_{e}/L_{p}$, varies from
1 to $10^{-4}$. The blue and orange dotted curves show the synchrotron and
inverse Compton components of the $L_{e}/L_{p}=1$ case.

###  5.2 Leptonic loading in TDEs

The EM cascade emissions discussed in this work are initiated by hadronic
processes without taking into account the injection of primary electrons or
positrons. In practice, we also expect the acceleration of leptons at the
acceleration sites where protons are energized. However, the accretion power
converted to leptons is typically less than that of protons. The leptonic
loading factor ($K_{e/p}$), defined as the ratio of the luminosities of
injected electrons/positrons to protons, is typically less than unity, e.g.,
$K_{e/p}\equiv L_{e}/L_{p}<1$. Simulations of non-relativistic diffusive shock
acceleration in supernova remnants have demonstrated that the efficiency of
electron energy deposition, compared to proton acceleration, is at least two
orders of magnitude lower, e.g., $K_{e/p}\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}10^{-2}$ (see, e.g.,
Jones, 2011; Morlino & Caprioli, 2012). Furthermore, the particle-in-cell
simulation gives a similar value, e.g., $K_{e/p}\simeq 10^{-3}$ (Katz &
Waxman, 2008; Caprioli & Spitkovsky, 2014; Park et al., 2015). For
relativistic shocks, a leptonic loading factor on the order of
$K_{e/p}\sim\mathcal{O}(10^{-4}-1)$ is widely used in active galactic nucleus
(AGN) models (roughly corresponding to the baryon loading factor $\xi_{\rm
cr}\sim 1-10^{4}$, see e.g., Murase et al., 2014; Kang et al., 2014; Gao et
al., 2019; Petropoulou et al., 2020). Therefore, in order to obtain a more
comprehensive understanding of the radiation zone, we compare the EM cascade
spectra with the synchrotron and inverse Compton emissions (in the synchrotron
self-Compton picture) from leptonic loading by varying $K_{e/p}$ from
$10^{-4}$ to the extreme value $K_{e/p}=1.0$. We consider a power law
electron/positron injection
$Q_{e}\propto\gamma_{e}^{-2}\exp(-\gamma_{e}/\gamma_{e,\rm max})$ with the
minimum and maximum Lorentz factors, $\gamma_{e,\rm min}=300$ and
$\gamma_{e,\rm max}=10^{5}$, comparable to AGNs (e.g., Böttcher et al., 2013).
Similarly, the injected electron distributions can be normalized via
$\int(\gamma_{e}m_{e}c^{2})Q_{e}d\gamma_{e}=L_{e}/V$. Given the magnetic field
strength $B=0.1$, we estimate the electron cooling Lorentz factor to be
$\gamma_{e,c}\simeq 6\pi m_{e}c/(t_{\rm fs}\sigma_{T}B^{2})\sim 4\times
10^{4}$, where $\sigma_{T}$ is the Thomson cross section. The cooling Lorentz
factor is significantly larger than $\gamma_{e,\rm min}$, which indicates that
the electrons are in the slow-cooling regime. The presence of slow-cooling
electrons together with a relatively weak magnetic field would result in an
enhanced inverse Compton component compared to the synchrotron emission.

In the M-IR case of AT2019dsg, the EM emissions from leptonic injections and
EM cascades are illustrated as the dashed and black solid curves in Fig. 6
(see captions for the meaning of each curve). We find that the leptonic
injection predicts a more broad spectrum from radio to $\gamma$-ray bands.
Moreover, the leptonic contribution (especially the inverse Compton component)
is comparable to the EM cascade for a moderate leptonic loading factor
$K_{e/p}\simeq 0.01$. The less efficient leptonic loading
($K_{e/p}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}0.01$) implies that the EM cascade
initiated by hadronic processes could dominate the EM emissions444This
conclusion is consistent with the secondary emissions in merging galaxies
(Yuan et al., 2019).. Meanwhile, the contributions from primary
electron/positron injections could become increasingly important for efficient
leptonic loading with $K_{e/p}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}0.01$. If we take the _Fermi_ upper limits
into account, the leptonic loading factor is constrained to be
$K_{e/p}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}0.01$, which falls within the expected
ranges from particle acceleration simulations or the AGN models. Near this
upper limit, relying solely on leptonic models or EM cascade models with
purely hadronic origins is insufficient to comprehend the multi-messenger
emissions from TDEs, emphasizing the necessity of a lepto-hadronic model. One
caveat in this discussion is that the limits on $K_{e/p}$ depend on the
assumptions of the electron/positron distributions. If a higher $\gamma_{e,\rm
min}$ or a stronger magnetic field is used, the electrons may be in the fast-
cooling regime, and the synchrotron component would become increasingly
important, leading to more stringent X-ray constraints. A more quantitative
discussion of leptonic injections is beyond the scope of this work and will be
provided in future lepto-hadronic models of TDEs.

##  6 Summary and conclusions

In this work, we have presented a comprehensive examination of the quasi-
isotropic EM cascade emissions in TDE radiation zones in a fully time-
dependent manner and investigated the interconnections between the neutrino
and EM cascade counterparts. We have demonstrated that in the observer’s
frame, EM cascades can be measured in the X-ray to $\gamma$-ray ranges at the
OUV peak time or exhibit a time delay depending on the $p\gamma$ optical
thickness. For instance, if the system is $p\gamma$ optically thin, the peak
of EM cascade emissions appears later than the OUV peak, and the time delay is
approximately in the order of the $p\gamma$ interaction time scale, ranging
from tenths to hundreds of days.

For AT2019dsg, an X-ray data point measured around 100 days after the OUV
peak, which is potentially incompatible with the early exponential decay of
the X-ray flux, can be described by the X-ray cascade emission in the M-IR
scenario, suggesting that the EM cascade emissions could be detectable for
Swift-XRT, XMM-Newton, and NICER. Additionally, we find that if no
$\gamma$-ray from the TDEs is seen, the source parameters, such as the radii
of the radiation zones and the maximum energies of injected protons, can be
stringently constrained. Specifically, for AT2019dsg and AT2019fdr, the
_Fermi_ $\gamma$-ray upper limits imply that the source may not be a very
efficient neutrino emitter, and the radiation zone is likely $p\gamma$
optically thin, where the predicted neutrino event rate (i.e., 0.01-0.1
neutrino events per TDE) is consistent with the expected neutrino numbers from
GFU searches (Stein et al., 2021; Reusch et al., 2022).

The joint analysis of the EM cascades and neutrino counterparts presented in
this paper provides an intriguing template for the multi-messenger studies of
TDEs. Future observations of the SEDs and light curves of X-ray, $\gamma$-ray,
and neutrino emissions will allow us to infer the physical conditions of the
radiation zones and test our EM cascade models. In addition to the TDEs
considered in this paper, the quasi-isotropic EM cascade models can be widely
applied to more TDEs that do not exhibit a direct connection to relativistic
jets. Moreover, EM cascades could play a significant role in jetted TDEs as
well, thereby highlighting the utility of coherent lepto-hadronic modeling.
Detailed discussions on this topic will be presented in our forthcoming works.

We would like to thank Maria Petropoulou for useful comments during the Lepto-
Hadronic Workshop in Bochum, Germany, and Claire Guepin for questions
motivating this work. We also extend our appreciation to Steffen Hallmann and
Xin-Yue Shi for their thorough internal review.

## References

  * Abbasi et al. (2023) Abbasi, R., Ackermann, M., Adams, J., et al. 2023, arXiv e-prints, arXiv:2304.01174, doi: 10.48550/arXiv.2304.01174
  * Alexander et al. (2020) Alexander, K. D., van Velzen, S., Horesh, A., & Zauderer, B. A. 2020, Space Sci. Rev., 216, 81, doi: 10.1007/s11214-020-00702-w
  * Andreoni et al. (2022) Andreoni, I., Coughlin, M. W., Perley, D. A., et al. 2022, Nature, 612, 430, doi: 10.1038/s41586-022-05465-8
  * Atwood et al. (2009) Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071, doi: 10.1088/0004-637X/697/2/1071
  * Ayala & HAWC Collaboration (2019) Ayala, H., & HAWC Collaboration. 2019, GRB Coordinates Network, 25936, 1
  * Baerwald et al. (2013) Baerwald, P., Bustamante, M., & Winter, W. 2013, ApJ, 768, 186, doi: 10.1088/0004-637X/768/2/186
  * Biehl et al. (2018a) Biehl, D., Boncioli, D., Fedynitch, A., & Winter, W. 2018a, A&A, 611, A101, doi: 10.1051/0004-6361/201731337
  * Biehl et al. (2018b) Biehl, D., Boncioli, D., Lunardini, C., & Winter, W. 2018b, Scientific Reports, 8, 10828, doi: 10.1038/s41598-018-29022-4
  * Böttcher et al. (2013) Böttcher, M., Reimer, A., Sweeney, K., & Prakash, A. 2013, ApJ, 768, 54, doi: 10.1088/0004-637X/768/1/54
  * Burrows et al. (2005) Burrows, D. N., Hill, J. E., Nousek, J. A., et al. 2005, Space Sci. Rev., 120, 165, doi: 10.1007/s11214-005-5097-2
  * Burrows et al. (2011) Burrows, D. N., Kennea, J. A., Ghisellini, G., et al. 2011, Nature, 476, 421, doi: 10.1038/nature10374
  * Cannizzaro et al. (2021) Cannizzaro, G., Wevers, T., Jonker, P. G., et al. 2021, MNRAS, 504, 792, doi: 10.1093/mnras/stab851
  * Caprioli & Spitkovsky (2014) Caprioli, D., & Spitkovsky, A. 2014, ApJ, 783, 91, doi: 10.1088/0004-637X/783/2/91
  * Cerruti et al. (2022) Cerruti, M., Kreter, M., Petropoulou, M., et al. 2022, in 37th International Cosmic Ray Conference, 979, doi: 10.22323/1.395.0979
  * Dai & Fang (2017) Dai, L., & Fang, K. 2017, MNRAS, 469, 1354, doi: 10.1093/mnras/stx863
  * Dai et al. (2015) Dai, L., McKinney, J. C., & Miller, M. C. 2015, ApJ, 812, L39, doi: 10.1088/2041-8205/812/2/L39
  * Dai et al. (2018) Dai, L., McKinney, J. C., Roth, N., Ramirez-Ruiz, E., & Miller, M. C. 2018, ApJ, 859, L20, doi: 10.3847/2041-8213/aab429
  * Fang et al. (2020) Fang, K., Metzger, B. D., Vurm, I., Aydi, E., & Chomiuk, L. 2020, ApJ, 904, 4, doi: 10.3847/1538-4357/abbc6e
  * Farrar & Gruzinov (2009) Farrar, G. R., & Gruzinov, A. 2009, ApJ, 693, 329, doi: 10.1088/0004-637X/693/1/329
  * Farrar & Piran (2014) Farrar, G. R., & Piran, T. 2014, arXiv e-prints, arXiv:1411.0704, doi: 10.48550/arXiv.1411.0704
  * Gao et al. (2019) Gao, S., Fedynitch, A., Winter, W., & Pohl, M. 2019, Nature Astronomy, 3, 88, doi: 10.1038/s41550-018-0610-1
  * Gao et al. (2017) Gao, S., Pohl, M., & Winter, W. 2017, ApJ, 843, 109, doi: 10.3847/1538-4357/aa7754
  * Guépin et al. (2018) Guépin, C., Kotera, K., Barausse, E., Fang, K., & Murase, K. 2018, A&A, 616, A179, doi: 10.1051/0004-6361/201732392
  * Hammerstein et al. (2023) Hammerstein, E., van Velzen, S., Gezari, S., et al. 2023, ApJ, 942, 9, doi: 10.3847/1538-4357/aca283
  * Hayasaki (2021) Hayasaki, K. 2021, Nature Astronomy, 5, 436, doi: 10.1038/s41550-021-01309-z
  * Hayasaki & Yamazaki (2019) Hayasaki, K., & Yamazaki, R. 2019, ApJ, 886, 114, doi: 10.3847/1538-4357/ab44ca
  * Hillas (1984) Hillas, A. M. 1984, ARA&A, 22, 425, doi: 10.1146/annurev.aa.22.090184.002233
  * Hümmer et al. (2010) Hümmer, S., Rüger, M., Spanier, F., & Winter, W. 2010, ApJ, 721, 630, doi: 10.1088/0004-637X/721/1/630
  * IceCube Collaboration (2019a) IceCube Collaboration. 2019a, GRB Coordinates Network, 25913, 1
  * IceCube Collaboration (2019b) —. 2019b, GRB Coordinates Network, 26258, 1
  * IceCube Collaboration et al. (2016) IceCube Collaboration, Aartsen, M. G., Abraham, K., et al. 2016, Journal of Instrumentation, 11, P11009, doi: 10.1088/1748-0221/11/11/P11009
  * Jansen et al. (2001) Jansen, F., Lumb, D., Altieri, B., et al. 2001, A&A, 365, L1, doi: 10.1051/0004-6361:20000036
  * Jiang et al. (2019) Jiang, Y.-F., Stone, J. M., & Davis, S. W. 2019, ApJ, 880, 67, doi: 10.3847/1538-4357/ab29ff
  * Jones (2011) Jones, T. W. 2011, Journal of Astrophysics and Astronomy, 32, 427, doi: 10.1007/s12036-011-9128-y
  * Kang et al. (2014) Kang, S.-J., Chen, L., & Wu, Q. 2014, ApJS, 215, 5, doi: 10.1088/0067-0049/215/1/5
  * Katz & Waxman (2008) Katz, B., & Waxman, E. 2008, J. Cosmology Astropart. Phys, 2008, 018, doi: 10.1088/1475-7516/2008/01/018
  * Liu et al. (2020) Liu, R.-Y., Xi, S.-Q., & Wang, X.-Y. 2020, Phys. Rev. D, 102, 083028, doi: 10.1103/PhysRevD.102.083028
  * Lunardini & Winter (2017) Lunardini, C., & Winter, W. 2017, Phys. Rev. D, 95, 123001, doi: 10.1103/PhysRevD.95.123001
  * McKinney et al. (2015) McKinney, J. C., Dai, L., & Avara, M. J. 2015, MNRAS, 454, L6, doi: 10.1093/mnrasl/slv115
  * Mohan et al. (2022) Mohan, P., An, T., Zhang, Y., et al. 2022, ApJ, 927, 74, doi: 10.3847/1538-4357/ac4cb2
  * Morlino & Caprioli (2012) Morlino, G., & Caprioli, D. 2012, A&A, 538, A81, doi: 10.1051/0004-6361/201117855
  * Murase et al. (2014) Murase, K., Inoue, Y., & Dermer, C. D. 2014, Phys. Rev. D, 90, 023007, doi: 10.1103/PhysRevD.90.023007
  * Murase et al. (2020) Murase, K., Kimura, S. S., Zhang, B. T., Oikonomou, F., & Petropoulou, M. 2020, ApJ, 902, 108, doi: 10.3847/1538-4357/abb3c0
  * Ohsuga et al. (2009) Ohsuga, K., Mineshige, S., Mori, M., & Kato, Y. 2009, PASJ, 61, L7, doi: 10.1093/pasj/61.3.L7
  * Park et al. (2015) Park, J., Caprioli, D., & Spitkovsky, A. 2015, Phys. Rev. Lett., 114, 085003, doi: 10.1103/PhysRevLett.114.085003
  * Pasham et al. (2023) Pasham, D. R., Lucchini, M., Laskar, T., et al. 2023, Nature Astronomy, 7, 88, doi: 10.1038/s41550-022-01820-x
  * Petropoulou et al. (2020) Petropoulou, M., Murase, K., Santander, M., et al. 2020, ApJ, 891, 115, doi: 10.3847/1538-4357/ab76d0
  * Predehl et al. (2021) Predehl, P., Andritschke, R., Arefiev, V., et al. 2021, A&A, 647, A1, doi: 10.1051/0004-6361/202039313
  * Reusch et al. (2022) Reusch, S., Stein, R., Kowalski, M., et al. 2022, Phys. Rev. Lett., 128, 221101, doi: 10.1103/PhysRevLett.128.221101
  * Rodrigues et al. (2019) Rodrigues, X., Gao, S., Fedynitch, A., Palladino, A., & Winter, W. 2019, ApJ, 874, L29, doi: 10.3847/2041-8213/ab1267
  * Rodrigues et al. (2021) Rodrigues, X., Garrappa, S., Gao, S., et al. 2021, ApJ, 912, 54, doi: 10.3847/1538-4357/abe87b
  * Rudolph et al. (2022) Rudolph, A., Petropoulou, M., Bošnjak, Ž., & Winter, W. 2022, arXiv e-prints, arXiv:2212.00765, doi: 10.48550/arXiv.2212.00765
  * Rudolph et al. (2023) Rudolph, A., Petropoulou, M., Winter, W., & Bošnjak, Ž. 2023, ApJ, 944, L34, doi: 10.3847/2041-8213/acb6d7
  * Saxton et al. (2020) Saxton, R., Komossa, S., Auchettl, K., & Jonker, P. G. 2020, Space Sci. Rev., 216, 85, doi: 10.1007/s11214-020-00708-4
  * Senno et al. (2017) Senno, N., Murase, K., & Mészáros, P. 2017, ApJ, 838, 3, doi: 10.3847/1538-4357/aa6344
  * Stein et al. (2021) Stein, R., van Velzen, S., Kowalski, M., et al. 2021, Nature Astronomy, 5, 510, doi: 10.1038/s41550-020-01295-8
  * van Velzen et al. (2016) van Velzen, S., Mendez, A. J., Krolik, J. H., & Gorjian, V. 2016, ApJ, 829, 19, doi: 10.3847/0004-637X/829/1/19
  * van Velzen et al. (2021a) van Velzen, S., Gezari, S., Hammerstein, E., et al. 2021a, ApJ, 908, 4, doi: 10.3847/1538-4357/abc258
  * van Velzen et al. (2021b) van Velzen, S., Stein, R., Gilfanov, M., et al. 2021b, arXiv e-prints, arXiv:2111.09391, doi: 10.48550/arXiv.2111.09391
  * Wang & Liu (2016) Wang, X.-Y., & Liu, R.-Y. 2016, Phys. Rev. D, 93, 083005, doi: 10.1103/PhysRevD.93.083005
  * Wang et al. (2011) Wang, X.-Y., Liu, R.-Y., Dai, Z.-G., & Cheng, K. S. 2011, Phys. Rev. D, 84, 081301, doi: 10.1103/PhysRevD.84.081301
  * Wen et al. (2020) Wen, S., Jonker, P. G., Stone, N. C., Zabludoff, A. I., & Psaltis, D. 2020, ApJ, 897, 80, doi: 10.3847/1538-4357/ab9817
  * Winter & Lunardini (2021) Winter, W., & Lunardini, C. 2021, Nature Astronomy, 5, 472, doi: 10.1038/s41550-021-01305-3
  * Winter & Lunardini (2023) —. 2023, ApJ, 948, 42, doi: 10.3847/1538-4357/acbe9e
  * Wu et al. (2022) Wu, H.-J., Mou, G., Wang, K., Wang, W., & Li, Z. 2022, MNRAS, 514, 4406, doi: 10.1093/mnras/stac1621
  * Yuan et al. (2020) Yuan, C., Murase, K., Kimura, S. S., & Mészáros, P. 2020, Phys. Rev. D, 102, 083013, doi: 10.1103/PhysRevD.102.083013
  * Yuan et al. (2019) Yuan, C., Murase, K., & Mészáros, P. 2019, ApJ, 878, 76, doi: 10.3847/1538-4357/ab1f06
  * Yuan et al. (2021) Yuan, C., Murase, K., Zhang, B. T., Kimura, S. S., & Mészáros, P. 2021, ApJ, 911, L15, doi: 10.3847/2041-8213/abee24
  * Zhang et al. (2017) Zhang, B. T., Murase, K., Oikonomou, F., & Li, Z. 2017, Phys. Rev. D, 96, 063007, doi: 10.1103/PhysRevD.96.063007
  * Zheng et al. (2022) Zheng, J.-H., Liu, R.-Y., & Wang, X.-Y. 2022, arXiv e-prints, arXiv:2209.11005, doi: 10.48550/arXiv.2209.11005


# Investigating the parton shower model in pythia8 with pp collision data at
$\mathrm{\surd{s}=13\,TeV}$ ††thanks: Presented at XXIV DAE-BRNS High Energy
Physics Symposium

S. K. Kundu1, T. Sarkar2, M. Maity3 sumankundu.rs@visva-
bharati.ac.intanmay.sarkar<EMAIL_ADDRESS>
(1,3 Visva-Bharati University, Santiniketan, India
2 National Central University ( NCU), Taiwan.
)

###### Abstract

Understanding the production of quarks and gluons in high energy collisions
and their evolution is a very active area of investigation. Monte carlo event
generator pythia8 uses the parton shower model to simulate such collisions and
is optimized using experimental observations. Recent measurements of event
shape variables and differential jet cross-sections in pp collisions at
$\mathrm{\surd{s}=13\,TeV}$ at the Large Hadron Collider have been used to
investigate further the parton shower model as used in pythia8.

##  1 Introduction

Matrix element calculations with fixed order treatment is not sufficient to
understand the production of quarks and gluons, collectively called partons,
in high energy collisions or their evolution into jets of hadrons. Comparison
with experimental results demand fully exclusive description of the final
states based on the shower evolution and hadronization. Such methods are
described through phenomenological models embedded in the shower Monte Carlo
(MC) codes.

pythia8 uses leading order(LO) calculations followed by ‘transverse momentum’
($p_{\perp}$) ordered parton shower[1] with $p^{2}_{\perp}$ as evolution
variable for the generation of $2\rightarrow n$ ($n\geq 2$) final states by
taking account initial (ISR) and final (FSR) state shower. Shower evolution
for a parton like $a\rightarrow bc$, is based on the standard (LO) DGLAP
splitting kernels and the branching probability expressed as:

$d\mathcal{P}_{a}=\frac{dp^{2}_{\perp}}{p^{2}_{\perp}}\sum_{b,c}\frac{\alpha_{s}(p^{2}_{\perp})}{2\pi}P_{a\rightarrow
bc}(z)dz$ (1.1)

where $\mathrm{P_{a\rightarrow bc}}$ is the DGLAP splitting function and
$p^{2}_{\perp}$ represents the scale of the branching; $z$ represents the
sharing of $p_{\perp}$ of $a$ between the two daughters, with $b$ taking a
fraction $z$ and $c$ the rest, $1-z$. Here the summation goes over all allowed
branchings, e.g. $q\rightarrow qg$ and $q\rightarrow q\gamma$ and etc. Now,
the divergence at $p^{2}_{\perp}\rightarrow 0$ is taken care of by introducing
a term
$\mathcal{P}^{no}_{a}(p^{2}_{{\perp}_{\textrm{max}}},p^{2}_{{\perp}_{\textrm{evol}}})$
known as _Sudakov form factor_ [2]. This Sudakov factor ensures that there
will be no emission between scale $p^{2}_{{\perp}_{\textrm{max}}}$ to a given
$p^{2}_{{\perp}_{\textrm{evol}}}$.

Considering lightcone kinematics, evolution variables
$\mathrm{p^{2}_{{\perp}_{\textrm{evol}}}}$ for $a\rightarrow bc$ at virtuality
scale $Q^{2}$ for space-like branching (ISR) and time-like branching (FSR) are
given by $(1-z)Q^{2}$ and $z(1-z)Q^{2}$ respectively. Finally, equations 1.2
and 1.3 describe the evolutions for ISR and FSR respectively [1].

$\displaystyle\textrm{d}\mathcal{P}_{b}$ $\displaystyle=$
$\displaystyle\frac{\textrm{d}p^{2}_{{\perp}_{\textrm{evol}}}}{p^{2}_{{\perp}_{\textrm{evol}}}}\frac{\alpha_{s}(p^{2}_{{\perp}_{\textrm{evol}}})}{2\pi}\frac{x^{\prime}f_{a}(x^{\prime},p^{2}_{{\perp}_{\textrm{evol}}})}{xf_{a}(x,p^{2}_{{\perp}_{\textrm{evol}}})}P_{a\rightarrow
bc}(z)dz\mathcal{P}^{\textrm{no}}_{b}(x,p^{2}_{{\perp}_{\textrm{max}}},p^{2}_{{\perp}_{\textrm{evol}}})$
(1.2) $\displaystyle\textrm{d}\mathcal{P}_{a}$ $\displaystyle=$
$\displaystyle\frac{\textrm{d}p^{2}_{{\perp}_{\textrm{evol}}}}{p^{2}_{{\perp}_{\textrm{evol}}}}\frac{\alpha_{s}(p^{2}_{{\perp}_{\textrm{evol}}})}{2\pi}P_{a\rightarrow
bc}(z)\textrm{d}z\mathcal{P}^{\textrm{no}}_{a}(p^{2}_{{\perp}_{\textrm{max}}},p^{2}_{{\perp}_{\textrm{evol}}})$
(1.3)

Currently both the running re-normalisation and factorisation shower scales,
i.e. the scales at which $\alpha_{s}$ and the PDFs are evaluated, are chosen
to be $p^{2}_{\perp_{\textrm{evol}}}$ [3]. The general methodology of pythia8
for ISR, FSR and MPI is to start from some maximum scale
$p^{2}_{{\perp}_{\textrm{max}}}$ and evolve downward in energy towards next
branching untill the daughter partons reach some cut-off.

pythia8 | Monash | Sampling range | Optimized
---|---|---|---
Parameters set | values |  | values
SpaceShower:alphaSvalue | 0.1365 | $0.1092-0.1638$ | $0.11409^{+0.00078}_{-0.00073}$
TimeShower:alphaSvalue | 0.1365 | $0.1092-0.1638$ | $0.15052^{+0.00077}_{-0.00076}$
SpaceShower:PTmaxFudge | 1.0 | $0.6-1.4$ | $0.9323^{+0.0065}_{-0.0064}$
Table 1: Optimized result of three parameters of pythia8 is shown along with
their default values in the Monash tune and the sampling range.

##  2 Optimizing the Parton Shower Model of PYTHIA8

CMS and ATLAS have done several tunning of pythia8 around its Monash tune [4]
for underlying event (UE), the strong coupling, and MPI related parameters [5,
6] [7]. In this study[8] Monash tune also used as default for pythia v8.235
with NNPDF2.3 PDF (LO) set to optimize with four event shapes [9] measurement
from CMS. These are - the complement of transverse thrust
($\mathrm{\tau_{\perp}}$), total jet mass ($\mathrm{\rho_{\text{Tot}}}$),
total transverse jet mass ($\mathrm{\rho^{\textrm{T}}_{\text{Tot}}}$) and
total jet broadening ($\mathrm{B_{\textrm{T}}}$).

Monash tune overestimates the multijet regions of these event shapes [9],
Hence we examined ISR and FSR utilising the provision that pythia8 allows the
use of separate values of $\mathrm{\alpha_{s}(M_{Z})}$ for the showering
frameworks used for these. The maximum evolution scale involved in the
showering is set to match the scale of the hard process itself. In pythia8 it
is set equal to the factorization scale, but allows its modification by
multiplicative factors SpaceShower:PTmaxFudge for ISR and
TimeShower:PTmaxFudge for FSR. The latter is seen not to have much effect on
the ESVs, it is excluded from the optimization.

For each point in the parameter space, resulting distributions have been
compared with data in terms of $\mathrm{\chi^{2}/NDF}$. Then professor v2.3.0
[10] along with RIVET v2.6 [11] has been used to optimize the complete set of
ESV distributions from pythia8 [9]. Post optimization, the new parameter set
is checked[8] using other relevant results from the CMS [12] and ATLAS [13].

Figure 1: Predictions of the optimized parameter set is compared with CMS data
and Monash tune for $H_{\textrm{T,2}}$ range $165<H_{\textrm{T,2}}<225$.
normalized distributions of the $\mathrm{\tau_{\perp}}$(top left),
$\mathrm{\rho_{\text{Tot}}}$(top right), $\mathrm{B_{\textrm{T}}}$(bottom
left) and $\mathrm{\rho^{\textrm{T}}_{\text{Tot}}}$(bottom right)

Figure 2: Normalized distributions of differential inclusive cross-section for
$\textrm{anti-}k_{\textrm{T}}$ jets (R=0.4) for CMS(left) and ATLAS(right) are
compared with the predictions of pythia8 with the optimized parameter set and
Monash tune.

##  3 Validation of results

The optimized values of the three parameters (see, table 1) are used to
calculate the ESVs. Agreement with data deteriorates slightly for
$\mathrm{\tau_{\perp}}$ and $\mathrm{\rho^{\textrm{T}}_{\text{Tot}}}$ (figure
1) compared to the good agreement with the Monash tune. But, there is
significant improvement in agreement with data for
$\mathrm{\rho_{\text{Tot}}}$ and $\mathrm{B_{\textrm{T}}}$ (figure 1)compared
to the Monash tune. Since $\mathrm{\rho_{\text{Tot}}}$ and
$\mathrm{B_{\textrm{T}}}$ had a rather poor agreement between data and the
Monash tune, overall this new set of parameters is better.

Inclusive jet cross-section measurements being sensitive to PDF of protons and
$\mathrm{\alpha_{s}}$ are also compared with those optimized values. CMS [12]
and ATLAS [13] studies with the 13 TeV data considered for this validation.
The CMS measurements of inclusive cross-sections for
$\textrm{anti-}k_{\textrm{T}}$ jets with R = 0.4, 0.7. Figures 2 show that the
new parameter set improves the agreement between data and the Monash tune of
Pythia8. Similar improvement is seen for the ATLAS measurement of
$\textrm{anti-}k_{\textrm{T}}$ jets with R = 0.4 (figure 2).

Since pythia8 is widely used, its optimization is important. This study shows
that certain aspects of the experimental observations can be better described
with this optimized set of parameters.

## References

  * [1] T. Sjostrand, P.Z. Skands, Eur. Phys. J. C39, 129 (2005).  DOI 10.1140/epjc/s2004-02084-y
  * [2] V.V. Sudakov, Sov. Phys. JETP 3, 65 (1956).  [Zh. Eksp. Teor. Fiz.30,87(1956)]
  * [3] R. Corke, T. Sjostrand, JHEP 03, 032 (2011).  DOI 10.1007/JHEP03(2011)032
  * [4] P. Skands, S. Carrazza, J. Rojo, Eur. Phys. J. C 74, 3024 (2014).  DOI 10.1140/epjc/s10052-014-3024-y
  * [5] V. Khachatryan, et al., Eur. Phys. J. C 76, 155 (2016).  DOI 10.1140/epjc/s10052-016-3988-x
  * [6] A.M. Sirunyan, et al., Eur. Phys. J. C 80(1), 4 (2020).  DOI 10.1140/epjc/s10052-019-7499-4
  * [7] A. Buckley, in _Proceedings of the Sixth International Workshop on Multiple Partonic Interactions at the Large Hadron Collider_. CERN (CERN, Geneva, 2014), p. 29
  * [8] S.K. Kundu, T. Sarkar, M. Maity, Int. J. Mod. Phys. A 34(33), 1950219 (2019).  DOI 10.1142/S0217751X19502191
  * [9] A.M. Sirunyan, et al., JHEP 12, 117 (2018).  DOI 10.1007/JHEP12(2018)117
  * [10] A. Buckley, H. Hoeth, H. Lacker, H. Schulz, J.E. von Seggern, Eur. Phys. J. C65, 331 (2010).  DOI 10.1140/epjc/s10052-009-1196-7
  * [11] A. Buckley, J. Butterworth, L. Lonnblad, D. Grellscheid, H. Hoeth, J. Monk, H. Schulz, F. Siegert, Comput. Phys. Commun. 184, 2803 (2013).  DOI 10.1016/j.cpc.2013.05.021
  * [12] V. Khachatryan, et al., Eur. Phys. J. C76(8), 451 (2016).  DOI 10.1140/epjc/s10052-016-4286-3
  * [13] Measurement of inclusive-jet cross-sections in proton-proton collisions at $\sqrt{s}=13\leavevmode\nobreak\ \mathrm{TeV}$ centre-of-mass energy with the ATLAS detector.  Tech. Rep. ATLAS-CONF-2016-092, CERN, Geneva (2016)