# Hydrodynamic Entropy and Emergence of Order in Two-dimensional Euler
Turbulence

Mahendra K. Verma<EMAIL_ADDRESS>Department of Physics, Indian Institute of
Technology Kanpur, Kanpur 208016, India  Soumyadeep Chatterjee
<EMAIL_ADDRESS>Department of Physics, Indian Institute of Technology
Kanpur, Kanpur 208016, India

###### Abstract

Using numerical simulations, we show that the asymptotic states of two-
dimensional (2D) Euler turbulence exhibit large-scale flow structures due to
nonzero energy transfers among small wavenumber modes. These asymptotic
states, which depend on the initial conditions, are out of equilibrium, and
they are different from the predictions of Onsager and Kraichnan. We propose
“hydrodynamic entropy” to quantify order in 2D Euler turbulence; we show that
this entropy decreases with time, even though the system is isolated with no
dissipation and no contact with a heat bath.

##  I Introduction

Euler turbulence remains unsolved till date. In this paper, we address the
energy flux and entropy of two-dimensional (2D) Euler turbulence. The
equations of incompressible Euler flow are

$\displaystyle\partial_{t}{\bf u}+{\bf u}\cdot\nabla{\bf u}=-\nabla
p;~{}~{}~{}\nabla\cdot{\bf u}=0,$ (1)

where ${\bf u},p$ are the velocity and pressure fields respectively [1, 2].
The above system is isolated, as it lacks external force and viscous
dissipation. Consequently, the thermodynamic entropy of an Euler flow remains
constant [1]. However, when we solve Euler equations with an ordered initial
condition, structurally, three-dimensional (3D) Euler turbulence becomes more
random during its evolution [3], whereas 2D Euler turbulence tends to become
more orderly [4, 5]. Hence, the entropy of Euler turbulence needs a
reexamination.

Onsager [6] modelled 2D Euler flow using a collection of point vortices
interacting via logarithmic potential. Onsager showed that for large energy,
2D Euler turbulence exhibits “negative temperature” and a large cluster of
same-circulation vortices. Recently, Gauthier  _et al._ [7] observed such
giant vortices in an experiment involving 2D quantum fluid, thus providing an
experimental verification of Onsager’s theory. Billam  _et al._ [8] developed
a first-principles realization of Onsager’s vortex model in a 2D superfluid.
Miller [9] and Robert [10] extended Onsager’s theory to continuum version of
2D Euler turbulence and computed entropy for the flow. Using tools of
equilibrium statistical mechanics, Bouchet and Simonnet [11], and Bouchet and
Venaille [12] derived multiple stationary states, namely a dipole and
unidirectional flow (shear layer), for 2D Euler turbulence. Pakter and Levin
[13] provide a contrary viewpoint and showed that 2D Euler turbulence is out
of equilibrium, and that a system of interacting vortices becomes trapped in a
nonequilibrium stationary state; these results deviate from the predictions of
Onsager [6]. Refer to review articles by Eyink and Sreenivasan [14], and
Bouchet and Venaille [12] for an extensive discussion.

Lee [15] and Kraichnan [16] provide an alternative framework for Euler
turbulence. They showed that the evolution of Fourier modes of Euler equation
follows Liouville’s theorem, and that the equilibrium solutions of Euler
turbulence are

$\displaystyle E(k)=\frac{k^{2}}{\beta-\gamma k^{2}}~{}\mathrm{for~{}3D};$ (2)
$\displaystyle E(k)=\frac{k}{\beta+\gamma k^{2}}~{}\mathrm{for~{}2D},$ (3)

where $\beta$ and $\gamma$ are constants. Here, the Fourier modes form a
microcanonical ensemble. The derivation of Eqs. (2,3) involves two competing
conservation laws: kinetic energy ($\int d{\bf r}u^{2}/2$) and kinetic
helicity ($\int d{\bf r}({\bf u}\cdot\bm{\omega})$) in 3D, and kinetic energy
and enstrophy ($\int d{\bf r}\omega^{2}/2$) in 2D, where $\bm{\omega}$ is the
vorticity field. For some combinations of energy and enstrophy, 2D Euler
turbulence yields $\beta<0$ or “negative temperature” [4, 17]. Kraichnan and
Montgomery [17] provide a detailed review of 2D Euler turbulence.

For a $\delta$-correlated random velocity field as an initial condition, both
2D and 3D Euler turbulence follow the energy spectra of Eqs. (3, 2) with
$\gamma\approx 0$ [18, 19]. In addition, for an initial condition with large-
scale structures, 3D Euler turbulence asymptotes to $E(k)$ of Eq. (2) [3].
However, $E(k)$ of 2D Euler turbulence differs from Eq. (3) for coherent
velocity field as an initial condition. For example, for enstrophy-dominated
2D Euler turbulence, Fox and Orszag [4] reported deviations from Eq. (3) at
small wavenumbers. For parameters where $\beta+\gamma k^{2}\approx 0$, Seyler
_et al._ [5] observed large vortex structures, similar to those in a discrete
vortex system [20]. Dritschel  _et al._ [21] studied the unsteady nature of 2D
flow structures on a sphere. Robert and Sommeria [22], and Bouchet and
Venaille [12] analyzed such structures in the framework of equilibrium
statistical mechanics.

The works of Fox and Orszag [4], Seyler  _et al._ [5], Pakter and Levin [13],
Bouchet and Simonnet [11], Dritschel  _et al._ [21], and Modin and Viviani
[23] indicate that 2D Euler turbulence is out of equilibrium, contrary to the
assumptions of Onsager [6] and Kraichnan [16]. Bouchet and Venaille [12]
argued that even though nonequilibrium steady states of 2D Euler turbulence
often break detailed balance, under weak force and zero viscosity, they may be
described by microcanonical measures and entropy functional. Bouchet and
coworkers [11, 12, 24], and Modin and Viviani [23] explained the structures of
2D Euler turbulence in this framework. In this paper, we advance this theme by
carefully examining the energy transfers and energy flux of 2D Euler
turbulence. We show that the small wavenumber modes exhibit nonzero energy
transfers, hence break the detailed balance, which is a stringent criterion
for equilibrium. Thus, we demonstrate the nonequilibrium nature of 2D Euler
turbulence. We also quantify the order of the structures using hydrodynamic
entropy.

The outline of the paper is as follows. In Sec. II, we describe nonequilibrium
nature of 2D Euler turbulence. We propose hydrodynamic entropy in Sec. III to
quantify this nature. We conclude the paper in Sec. IV.

##  II Nonequilibrium Nature of 2D Euler Turbulence

Prior to a detailed discussion on 2D Euler turbulence, we summarize the energy
spectrum and flux of 3D Euler turbulence. Cichowlas  _et al._ [3] simulated 3D
Euler turbulence with Taylor-Green vortex as an initial condition. For such
simulations, in the early phase, the energy flows from large scales to small
scales, and the energy flux is positive. After several eddy turnover times,
the system approaches equilibrium with vanishing energy flux. Refer to
Appendix A for details.

For 2D Euler turbulence, we performed pseudo-spectral simulations on a
$(2\pi)^{2}$ box with a $M^{2}$ grid. Here, $M=512$. We dealiase the code by
setting all the modes outside the sphere of radius $M/3$ to zero. To conserve
the total energy, we time evolve Eq. (1) using position-extended Forest-Ruth-
like (PEFRL) scheme [25, 26] with time step = $10^{-4}$. We carried out three
runs with the following initial conditions:

  1. 1.

Run A: The initial velocity profile is taken as $(\sin 11x\cos
11y+\eta_{x},-\cos 11x\sin 11y+\eta_{y})$, where $(\eta_{x},\eta_{y})$ is
random noise. We take $|\eta_{x}|\ll 1$ and $|\eta_{y}|\ll 1$.

  2. 2.

Run B: The initial nonzero velocity Fourier modes are ${\bf u}(1,0)=(0,1)$,
${\bf u}(0,1)=(1,0)$, and ${\bf u}(1,1)=(-i,i)$.

  3. 3.

Run C: The initial nonzero velocity Fourier modes are ${\bf u}(10,0)=(0,1)$,
${\bf u}(0,10)=(1,0)$, and ${\bf u}(10,10)=(-i,i)$.

Figure 1: For Runs A, B, and C of 2D Euler turbulence: (a,c,e) the initial
states, (b,d,f) the final states respectively. Here we plot the velocity field
over the density plots of the vorticity field.

We time advance Runs A, B, C up to 170, 30, and 120 turnover times
($2\pi/U_{\mathrm{rms}}$) respectively. We observe Runs B and C reach the
steady flow profiles after 10 to 20 eddy turnover times, whereas Run A reaches
the steady state after 100 eddy turnover times. We label these states as
asymptotic states. In this paper, we show that these states are out of
equilibrium. The dependence of the asymptotic states on initial condition is
consistent with earlier works, e.g., Celani  _et al._ [27]. The top and bottom
panels of Fig. 1 illustrate respectively the initial and asymptotic states of
the three runs. Here, the velocity field is superposed over the density plots
of the vorticity field. Runs A, B, C asymptote respectively to a vortex-
antivortex pair [6], a unidirectional flow (shear layer), and four vortex-
antivortex pairs. The above large-scale flow structures are embedded in small-
scale noisy flow. Similar structures have been observed by Bouchet and
Simonnet [11], and Dritschel  _et al._ [21].

For Runs A, B, and C, total energy $E=\int d\textbf{r}(u^{2}/2)/\int
d\textbf{r}=0.2500954,4,4$, while total enstrophy $\Omega=\int
d\textbf{r}(\omega^{2}/2)/\int d\textbf{r}=62.17,6,600$, respectively; and
these quantities are conserved. For validation of our numerical codes, we plot
the time series of the total energy and total enstrophy, along with their
relative errors in Figs. 2 and 3. The relative errors for energy and enstrophy
are defined as

$\epsilon(t)=\frac{|E(t)-E(t=0)|}{E(t=0)};~{}~{}~{}~{}\epsilon_{\Omega}(t)=\frac{|\Omega(t)-\Omega(t=0)|}{\Omega(t=0)}.$
(4)

Since the three runs take different amounts of time to reach their respective
asymptotic states, in the plots, we employ normalized time,
$t^{\prime}=3t/17,t,t/4$, for the Runs A, B, and C respectively. For the
asymptotic states of Runs A, B, and C, the relative errors in energy
($\epsilon$) are $7.3\times 10^{-11}$, $1.9\times 10^{-11}$ and $4.5\times
10^{-8}$ respectively, while those in entropy ($\epsilon_{\Omega}$) are
$6.3\times 10^{-9}$, $2.7\times 10^{-7}$, and $6.5\times 10^{-6}$
respectively.

Figure 2: (a) For Runs A, B, C of 2D Euler turbulence, (a) time series of the
total energy, $E$, and (b) its relative error, $\epsilon$. We employ
normalized time $t^{\prime}=3t/17,t,t/4$ for the three runs.  Figure 3:  For
Runs A, B, C of 2D Euler turbulence, (a) time series of total enstrophy,
$\Omega$, and (b) its relative error, $\epsilon_{\Omega}$, as a function of
$t^{\prime}$, where $t^{\prime}=3t/17,t,t/4$ for the three runs.

In Fig. 4(a,b), we plot the averaged energy spectra and fluxes of the three
runs. Here, we average 6500, 1000, and 10000 frames of Runs A, B, and C in the
respective time intervals $(105,170),(20,30)$ and $(20,120)$. For intermediate
and large wavenumbers, $E(k)$ of the three runs follow Kraichnan’s predictions
[Eq. (3)]:

$\displaystyle\mathrm{Run~{}A:}$ $\displaystyle
E(k)=\frac{k}{-3230+237k^{2}}~{}\mathrm{for}~{}k>10,$ (5)
$\displaystyle\mathrm{Run~{}B:}$ $\displaystyle
E(k)=\frac{k}{-6357840+9361k^{2}}~{}\mathrm{for}~{}k>40,$ (6)
$\displaystyle\mathrm{Run~{}C:}$ $\displaystyle
E(k)=\frac{k}{888+26k^{2}}~{}\mathrm{for}~{}k>10,$ (7)

with constants having significant errors. We construct these functions using
non-linear least squares fit of the numerical data of $k/E(k)$ to
$\beta+\gamma k^{2}$. We employ the standard scipy.optimize.curve_fit function
of Python for the computation.

Figure 4: For the Runs A, B, C of 2D Euler turbulence: plots of the averaged
energy spectra, $E(k)$, and fluxes, $\Pi(k)$, of the asymptotic states. In
(a), the best-fit curves of Eqs. (5-7) are shown as dashed curves. Figure (b)
exhibits $\Pi(k)$ for small $k$’s, while inset shows $\Pi(k)$ for the whole
range. Figure 5:  For the Runs A, B, C of 2D Euler turbulence, plots of the
averaged energy spectra for small wavenumbers $k\in[1,10]$. The black-dashed
line represents $k^{-5/3}$ spectrum. Clearly, Runs B and C do not exhibit
$k^{-5/3}$ spectrum, while Run A matches with $k^{-5/3}$ for very small
region.

Based on Eqs. (5-7) and Figure 4(a), we claim that the large wavenumber modes
are in equilibrium. Fjørtoft [28] and Nazarenko [29] showed that for 2D
hydrodynamic turbulence, $k_{E}\leq\sqrt{\Omega/E}\leq k_{\Omega}$, where
$k_{E}$ and $k_{\Omega}$ are the centroids of energy and enstrophy
respectively. For Runs A, B, and C, $\sqrt{\Omega/E}=15.8,1.2,12.2$. Hence,
the wavenumbers far beyond $k_{\Omega}$ are dominated by enstrophy. In this
regime, our simulations reveal that $E(k)\propto k^{-1}$, which corresponds to
an equipartition of enstrophy. That is, modal enstrophy, $E_{\omega}({\bf
k})=E({\bf k})k^{2}=\mathrm{constant}$, which leads to $E(k)\sim 2\pi
k/k^{2}\sim k^{-1}$. Moreover, Runs A and B exhibit $\beta<0$ or “negative
temperature” [6, 16, 17], which is related to the emergence of large scale
structures. Note that Onsager [6] and Kraichnan [16] assumed the 2D Euler flow
to be in equilibrium.

Equations (5-7) and Figure 4(a) reveal that the small wavenumber modes deviate
strongly from Kraichnan’s predictions for equilibrium Euler turbulence,
consistent with the works of Fox and Orszag [4]. In the following discussion,
we will show that the small-wavenumber modes of 2D Euler are out of
Equilibrium from various perspectives. As illustrated in Fig. 5, for small
wavenumbers, $E(k)$’s do not follow $k^{-5/3}$ energy spectrum, except for a
small range of wavenumbers for Run A. This is expected because Euler
turbulence is very different from 2D hydrodynamic turbulence, which exhibits
$E(k)\propto k^{-5/3}$ when the intermediate scales are forced [30].

We complement the spectral studies with a quantification of energy transfers
and fluxes. Euler turbulence lacks external forcing and dissipation, hence the
temporal evolution of $E(k)$ is given by [2, 31, 32, 33, 34]

$\frac{dE(k,t)}{dt}=T(k,t),$ (8)

where $T(k,t)$ is the nonlinear energy transfer to wavenumber shell $k$. If
Euler turbulence were to be in equilibrium, then for all $k$’s,
$dE(k,t)/dt=0$, implying that $T(k,t)=0$. Consequently, the energy flux,
$\Pi(k,t)=-\int_{0}^{k}T(k^{\prime},t)dk^{\prime}$ would vanish for all
wavenumber spheres. This test is an alternative one to the entropy
maximization principle [33]. Energy conserving systems, as well as Hamiltonian
systems, are expected to approach equilibrium asymptotically [35, 31]. It has
been shown that 3D Euler turbulence asymptotically reaches the equilibrium
state with zero energy flux [3]. But, we show below that $T(k,t)$ and
$\Pi(k,t)$ for 2D Euler turbulence are nonzero for small $k$’s, thus 2D Euler
turbulence is out of equilibrium.

For 2D Euler turbulence, we illustrate $T(k,t)$ for the dominant wavenumber
shells ($k$ = 1 to 4) in Fig. 6. Clearly, these $T(k,t)$’s fluctuate
significantly, and $dE(k,t)/dt\approx T(k,t)$, thus validating Eq. (8). For
small $k$’s, the nonzero $T(k,t)$’s yield negative energy flux. In particular,
$\min[\left\langle\Pi(k)\right\rangle]\approx-3\times
10^{-4},-10^{-3},-8.4\times 10^{-4}$ for Runs A, B, C respectively. Note that
the finite $T(k,t)$ and $\Pi(k,t)$ for 2D turbulence are much larger than the
corresponding quantities of 3D Euler turbulence (see Appendix A). The nonzero
$T(k,t)$ and $\Pi(k,t)$ break the detailed balance of energy transfers and
indicate nonequilibrium nature of 2D Euler turbulence.

Figure 6:  For Runs A, B, C, plots of the time evolutions of $dE(k)/dt$ (solid
lines) and $T(k)$ (dashed lines) for three different wavenumbers.

Based on the above observations, we conclude that the intermediate and large
wavenumber shells are in equilibrium, but the small wavenumber shells are out
of equilibrium. However, there is an exception to the above rule. For
$\delta$-correlated velocity field as an initial condition, 2D Euler
turbulence exhibits equilibrium solution with $E(k)\propto k$ [19]. For this
particular case, $\gamma\approx 0$, and $k_{\Omega}$ exceeds the grid size.

In classical literature of hydrodynamics, the thermodynamic entropy of Euler
turbulence is taken to be constant [1]. However, as discussed above, the
disorder in Euler turbulence varies with time. Hence, the thermodynamic
entropy cannot capture the disorder in Euler turbulence. In Section III, we
define “hydrodynamic entropy” that can describe the disorder in Euler
turbulence.

##  III Hydrodynamic Entropy

Consider two gaseous systems shown in Fig. 7. In Fig. 7(a), the molecules are
in thermal equilibrium, and they move randomly with thermal speed $C_{s}$.
Thermodynamic entropy provides a good measure for the disorder in such a
system.

Figure 7: (a) In a 2D gas in thermal equilibrium, the molecules move randomly
with the sound speed $C_{s}$. (b) A schematic diagram illustrating the
velocities of molecules in a hydrodynamic vortex.

In Fig. 7(b), the molecules co-move in a hydrodynamic vortex. This system is
out of equilibrium because 2D hydrodynamic vortices coalesce to form large
vortices via inverse cascade of energy. To quantify disorder in an
hydrodynamic flow, we scale separate the microscopic thermal processes and the
macroscopic fluid processes, and define two different entropies for them [36,
18]. Thermodynamic entropy is employed for the microscopic processes [37, 38],
whereas “hydrodynamic entropy”, to be defined below, for describing
macroscopic order. For an Euler flow, only hydrodynamic entropy would be
meaningful because thermodynamic entropy is zero for this case.

We employ Shannon entropy [39] to quantify the disorder of the fluid
structures at macroscopic level. For the same, we postulate that the
probability of occurrence of a Fourier mode with wavenumber ${\bf k}$ is
$p_{\bf k}=E({\bf k})/E$, where $E({\bf k})$ is the modal energy, and $E$ is
the total energy. Now, the hydrodynamic entropy of the flow is defined as

$S=-\sum_{\bf k}p_{\bf k}\log_{2}(p_{\bf k}).$ (9)

The above entropy, which is defined for a snapshot, can be applied to any
fluid flow. Note that hydrodynamic entropy is very different from
thermodynamic entropy, which depends on temperature and volume of the system.
At present, it is not apparent if hydrodynamic entropy and thermodynamic
entropy could be put on a similar footing. It would be interesting to explore
whether we can relate the two in a way Landauer [40] connected computation to
entropy.

It is important to keep in mind that Shannon entropy [Eq. (9)] has been used
earlier to quantify order in images, music, messages, DNA, as as well as in
variety of flow phenomena, e.g., boundary layers, transitions, etc. [41].
There are many other measures of entropy. For example, Clark  _et al._ [42]
employed Kolmogorov-Sinai entropy for quantifying 2D hydrodynamic flow,
whereas Drivas and Elgindi [43] quantified entropy of Euler flow in terms of
accessible phase space.

Euler flow is an ideal case where thermal processes are ignored (due to zero
viscosity). Therefore, Euler flow has only hydrodynamic entropy. The entropies
of the vortex solutions of Onsager [6] and Miller [9] are hydrodynamic,
however, these measures may differ from those computed using Eq. (9).
Interestingly, Kraichnan’s equilibrium solution with $\gamma=0$ corresponds to
$\delta$–correlated velocity field [15, 16, 18]. For such a flow, $p_{\bf
k}=1/M$ and the hydrodynamic entropy $S=\log_{2}(M)$, where $M$ is the number
of modes of the system; this is the maximum possible hydrodynamic entropy for
a flow with $M$ degrees of freedom.

Before reporting the hydrodynamic entropy of 2D Euler turbulence, we describe
the hydrodynamic entropy of 3D Euler turbulence. For the 3D Euler flow
described in Appendix A, we compute the entropy using Eq. (9) and plot its
time series in Fig. 8. We derive the functional dependences of $S(t)$ using
the non-linear least squares method. In the initial phase, the hydrodynamic
entropy increases exponentially, after which it approaches the maximum
possible value $S=\log_{2}(M)=18.3$, where $M=(4\pi/3)(128/3)^{3}$ is the
number of modes of the system. We observe that for 3D Euler turbulence, the
approach to equilibrium is slow. Since $S$ increases monotonically in time, we
claim that 3D Euler turbulence evolves from order to disorder.

Figure 8: In 3D Euler turbulence, hydrodynamic entropy exhibits monotonic
growth, with exponential increase in the beginning and saturation in the end.

Now, we compute the hydrodynamic entropies of the 2D Euler flows of Runs A, B,
and C, and plot the entropy time series in Fig. 9. The duration and time
scales of the three runs are quite different. Hence, for a proper comparison,
we normalize the time appropriately, i.e., $t^{\prime}=3t/17,t,t/4$ for Runs
A, B, and C respectively. We observe that for each case, the entropy
fluctuates in the early stages, after which it decreases exponentially to an
asymptotic value. Here too, we employ the non-linear least squares method to
derive the functions that describe the decrease in $S$. The asymptotic
entropies for Runs A, B, and C are 4.9, 1.2, and 3.1 respectively. These
values are smaller than the maximum possible value, which is
$\log_{2}(M)\approx 16.5$, where $M=\pi(512/3)^{2}$ is the number of active
dealiased modes. Note, however, that $S$ for all the runs exhibit small
fluctuations in the asymptotic regime due to the dynamic nature of the flow.
Recall the temporal fluctuations of $T(k,t)$ and $E(k,t)$ described in Section
II. It is also important to note that the evolution of hydrodynamic entropy
for 2D Euler turbulence is very different that of 3D Euler turbulence.

Figure 9: For Runs A, B, C of 2D Euler turbulence, the temporal evolution of
hydrodynamic entropies with $t^{\prime}=3t/17,t,t/4$ respectively. In each
case, after initial transients, the entropy decreases with time and asymptotes
to an approximate constant value.

The asymptotic state of Run B is the most ordered one (having least entropy)
among the three runs. In the asymptotic regime of Run B, the Fourier modes
${\bf u}(0,\pm 1)$ contains 99% of the total energy, hence, their
probabilities are approximately 1/2 each. Therefore, these two modes yield
$S\approx\log_{2}(2)=1$, while the other modes contribute the rest (0.2). To
quantify the evolution of these modes, in Fig. 10, we plot the time series of
$\Re[{u_{x}(0,1)}]$ and $\Im[{u_{x}(0,1)}]$ of Run B. As shown in the figure,
$\Im[{u_{x}(0,1)}]=0$, while $\Re[{u_{x}(0,1)}]$ approaches nearly a constant
value for $t>10$. The fluctuations in $u_{x}(0,1)$ is consistent with those in
$T(k,t)$ discussed in Section II.

Figure 10: For Run B of 2D Euler turbulence, time series of
$\Re[{u_{x}(0,1)}]$ and $\Im[{u_{x}(0,1)}]$. In this run, ${\bf u}(0,1)$ is
the most dominant Fourier mode.

In contrast, for Run C, the seven most energetic modes yield an approximate
entropy of 2.3, while the rest of them contribute the remaining 0.8. The
Fourier modes of Run A are more wide-spread than those of Runs B and C, hence
Run A has the maximum entropy among the three runs. Note that the shear layer
of Run B has the least entropy, but the vortex-antivortex pair of Run A has
the maximum entropy, which is still smaller than the maximum possible value of
16.5.

Thus, we show that the hydrodynamic entropy of 2D Euler turbulence decreases
with time for a significant duration, even though the system is isolated.
Thus, 2D Euler turbulence is a rare isolated system that exhibits evolution
from disorder to order.

##  IV Discussions and Conclusions

Kraichnan [16] and Onsager [6] assumed that 2D Euler turbulence reaches
equilibrium. Kraichnan argued that 2D Euler turbulence exhibits the energy
spectrum of Eq. (3), whereas Onsager [6] advocated a large cluster of same-
circulation vortices. We performed numerical simulations to test whether 2D
Euler turbulence is in equilibrium or out of equilibrium.

Our numerical simulations, as well as several past ones, report that at small
wavenumbers, the energy spectrum and flux of 2D Euler turbulence differ from
those predicted by Kraichnan [16]. These differences arise due to the
nonequilibrium nature of 2D Euler turbulence. It is important to note that 2D
Euler turbulence is one of the few energy-conserving systems that do not
thermalize (or approach equilibrium). Note, however, that several experiments
on quantum fluids report consistency with Onsager’s theory [7]. Hence, the
relationship between Euler turbulence and quantum fluids needs to be examined
carefully. We also remark that the flow structures of 2D Euler turbulence do
not exhibit long-range order, as in a ferromagnet. Instead, the flow
structures resemble vortex-antivortex pairs of Berezinskii–Kosterlitz–Thouless
transition [44].

More importantly, we propose hydrodynamic entropy to quantify disorder in
Euler turbulence. We observe that the hydrodynamic entropy of 3D Euler
turbulence increases monotonically with time, whereas it decreases for 2D
Euler turbulence. Thus, 2D Euler turbulence is a unique isolated system that
exhibits evolution form disorder to order. This feature arises due to the
inverse energy cascade, which is a property of 2D hydrodynamics [30, 45, 2,
32]. Hence, the emergence of hydrodynamic order in 2D Euler turbulence has a
dynamic origin. Note that the thermodynamic entropy of Euler turbulence
remains constant throughout its evolution. Therefore, the decrease in
hydrodynamic entropy with time does not violate second law of thermodynamics.

Euler equation is time reversible due to an absence of viscous dissipation
[2]. This is the reason why thermodynamic entropy of Euler turbulence is
constant [1]. Note, however, that the solutions of Euler equation exhibits
irreversibility due to their chaotic and nonequilibrium nature. As we show in
this paper, 2D Euler turbulence is a special energy-conserving system that
exhibit nonequilibrium behaviour. However, 3D Euler turbulence reaches
equilibrium asymptotically where detailed balance is preserved statistically
[3]. Interestingly, the hydrodynamic entropy captures the irreversibility and
disorder of 2D and 3D Euler turbulence quite well.

We conclude this paper by emphasizing that thermal and hydrodynamic processes,
with or without viscosity, are multiscale phenomena where hierarchical energy
transfers play a critical role [2, 34]. Recently, Verma [36, 18] attributed
irreversibility in a turbulent flow to the asymmetric energy transfers (e.g.,
forward cascade in 3D Navier-Stokes equation), which is a hydrodynamic
property. Following a similar approach, in this paper, we propose hydrodynamic
entropy that successfully captures the evolution of 2D Euler turbulence from
disorder to order.

## Acknowledgements

The authors thank Arul Lakshminarayan, Siva Chandran, Shashwat Bhattacharya,
Pankaj Mishra, Sagar Chakraborty, and Anurag Gupta for useful discussions.
This work is supported by the project 6104-1 from the Indo-French Centre for
the Promotion of Advanced Research (IFCPAR/CEFIPRA). Soumyadeep Chatterjee is
supported by INSPIRE fellowship (IF180094) from Department of Science &
Technology, India.

##  Appendix A Evolution of 3D Euler Turbulence

In this Appendix, we summarize the energy spectrum and flux of 3D Euler
turbulence. Cichowlas  _et al._ [3] simulated 3D Euler turbulence with Taylor-
Green vortex as an initial condition. For such simulations, in the early
phase, the energy flows from large scales to small scales, and the energy flux
is positive. After several eddy turnover times, the system approaches
equilibrium with vanishing energy flux and $E(k)$ given by Eq. (2).

To compute the hydrodynamic entropy for 3D Euler turbulence, we simulated 3D
Euler flow. We performed our simulation on a $(2\pi)^{3}$ box with a $M^{3}$
grid. Here, $M=128$. We dealiase the code by setting all the modes outside the
sphere of radius $M/3$ to zero. We time evolve Eq. (1) of the main text using
position-extended Forest-Ruth-like (PEFRL) scheme [25, 26] to ensure energy
conservation. As in Cichowlas  _et al._ [3], we take Taylor-Green vortex
($k_{0}=1$) as an initial condition, and time evolve the system till 120 eddy
turnover time ($2\pi/U_{\mathrm{rms}}$) with a constant $dt=10^{-4}$. The
total energy per unit volume, $E=\int d\textbf{r}(u^{2}/2)/\int
d\textbf{r}=0.125$, and it is conserved up to 12 significant digits till the
final time $t=120$. The system reaches equilibrium in approximately 100 eddy
turnover time, consistent with the estimates of Cichowlas  _et al._ [3] and
Verma  _et al._ [19]. In the asymptotic state, the velocity field is
$\delta$-correlated, and it is as random as that in a thermodynamic gas [19].

Figure 11: For 3D Euler turbulence, plots of (a) the energy spectra $E(k)$,
and (b) the energy fluxes $\Pi(k)$ at $t=0,10,30,120$. At $t=120$,
$E(k)\propto k^{2}$.

Using the numerical data, we compute the energy spectrum and flux for the flow
at times $t=0,10,30,120$. For the computation of energy flux, we employ the
algorithm outlined in Dar  _et al._ [46] and Verma [47]. In Fig. 11(a,b), we
exhibit the energy spectra and energy fluxes at $t=0,10,30,120$. Clearly, with
time, the energy spectrum $E(k)$ spreads from low wavenumbers to high
wavenumbers, and asymptotes to $k^{2}$, consistent with the predictions of
Kraichnan [16] [Eq. (2) with $\gamma=0$].

The energy flux is nonzero for the intermediate configurations, but for $t\geq
30$, $\Pi(k)\approx 0$ for small $k$’s. Note, however, that $\Pi(k)$ is of the
order of $10^{-4}$ for large wavenumbers. The fluctuations in the energy flux
is suppressed significantly on averaging over 2000 frames in time interval
(100,120) (see Fig. 12). We expect these fluctuations to subside at a later
time when the system has thermalized fully.

Figure 12: For 3D Euler turbulence, plots of (a) averaged $T(k,t)$, and (b)
averaged $\Pi(k,t)$. We average over 2000 frames in time interval (100, 120).

We also remark that for $\delta$-correlated initial velocity profile,
$E(k)\propto k^{2}$ [Eq. (2) with $\gamma=0$], and energy flux vanishes from
the beginning itself [19]. Thus, Euler turbulence remains thermalized
throughout for $\delta$-correlated initial condition.

## References

  * Landau and Lifshitz [1987] L. D. Landau and E. M. Lifshitz,  _Fluid Mechanics_ , 2nd ed., Course of Theoretical Physics (Elsevier, Oxford, 1987).
  * Frisch [1995] U. Frisch,  _Turbulence: The Legacy of A. N. Kolmogorov_ (Cambridge University Press, Cambridge, 1995).
  * Cichowlas  _et al._ [2005] C. Cichowlas, P. Bonaïti, F. Debbasch, and M. E. Brachet, Effective Dissipation and Turbulence in Spectrally Truncated Euler Flows, Phys. Rev. Lett. 95, 264502 (2005).
  * Fox and Orszag [1973] D. G. Fox and S. A. Orszag, Inviscid dynamics of two-dimensional turbulence, Phys. Fluids 16, 169 (1973).
  * Seyler  _et al._ [1975] C. E. Seyler, Y. Salu, D. Montgomery, and G. Knorr, Two-dimensional turbulence in inviscid fluids or guiding center plasmas, Physics of Fluids 18, 803 (1975).
  * Onsager [1949] L. Onsager, Statistical hydrodynamics, Il Nuovo Cimento 6, 279 (1949).
  * Gauthier  _et al._ [2019] G. Gauthier, M. T. Reeves, X. Yu, A. S. Bradley, M. A. Baker, T. A. Bell, H. Rubinsztein-Dunlop, M. J. Davis, and T. W. Neely, Giant vortex clusters in a two-dimensional quantum fluid, Science 364, 1264 (2019).
  * Billam  _et al._ [2014] T. P. Billam, M. T. Reeves, B. P. Anderson, and A. S. Bradley, Onsager-Kraichnan Condensation in Decaying Two-Dimensional Quantum Turbulence, Phys. Rev. Lett. 112, 145301 (2014).
  * Miller [1990] J. Miller, Statistical mechanics of Euler equations in two dimensions, Phy. Rev. Lett. 65, 2137 (1990).
  * Robert [1991] R. Robert, A maximum-entropy principle for two-dimensional perfect fluid dynamics, J. Stat. Phys. 65, 531 (1991).
  * Bouchet and Simonnet [2009] F. Bouchet and E. Simonnet, Random Changes of Flow Topology in Two-Dimensional and Geophysical Turbulence, Phys. Rev. Lett. 102, 094504 (2009), 0804.2231 .
  * Bouchet and Venaille [2012] F. Bouchet and A. Venaille, Statistical mechanics of two-dimensional and geophysical flows, Phys. Rep. 515, 227 (2012).
  * Pakter and Levin [2018] R. Pakter and Y. Levin, Nonequilibrium Statistical Mechanics of Two-Dimensional Vortices, Phys. Rev. Lett. 121, 020602 (2018).
  * Eyink and Sreenivasan [2006] G. L. Eyink and K. R. Sreenivasan, Onsager and the theory of hydrodynamic turbulence, Rev. Mod. Phys. 78, 87 (2006).
  * Lee [1952] T. D. Lee, On some statistical properties of hydrodynamical and magneto-hydrodynamical fields, Quart. Appl. Math. 10, 69 (1952).
  * Kraichnan [1973] R. H. Kraichnan, Helical turbulence and absolute equilibrium, J. Fluid Mech. 59, 745 (1973).
  * Kraichnan and Montgomery [1980] R. H. Kraichnan and D. C. Montgomery, Two-dimensional turbulence, Rep. Prog. Phys. 43, 547 (1980).
  * Verma [2020] M. K. Verma, Boltzmann equation and hydrodynamic equations: their equilibrium and non-equilibrium behaviour, Phil. Trans. R. Soc. A. 378, 20190470 (2020).
  * Verma  _et al._ [2020] M. K. Verma, S. Bhattacharya, and S. Chatterjee, Euler Turbulence and thermodynamic equilibrium , arXiv , arXiv:2004.09053 (2020).
  * Joyce and Montgomery [1973] G. Joyce and D. Montgomery, Negative temperature states for the two-dimensional guiding-centre plasma, J. Plasma Phys. 10, 107 (1973).
  * Dritschel  _et al._ [2015] D. G. Dritschel, W. Qi, and J. B. Marston, On the late-time behaviour of a bounded, inviscid two-dimensional flow, J. Fluid Mech. 783, 1 (2015).
  * Robert and Sommeria [1991] R. Robert and J. Sommeria, Statistical equilibrium states for two-dimensional flows, J. Fluid Mech. 229, 291 (1991).
  * Modin and Viviani [2021] K. Modin and M. Viviani, Canonical scale separation in 2d incompressible hydrodynamics, arXiv , 2102.01451 (2021).
  * van Kan  _et al._ [2022] A. van Kan, A. Alexakis, and M. Brachet, Geometric microcanonical theory of two-dimensional truncated Euler flows, Phil. Trans. R. Soc. A 380, 20210049 (2022).
  * Omelyan  _et al._ [2002] I. P. Omelyan, I. M. Mryglod, and R. Folk, Optimized Forest-Ruth- and Suzuki-like algorithms for integration of motion in many-body systems, Computer Physics Communications 146, 188 (2002).
  * Forest and Ruth [1990] E. Forest and R. D. Ruth, Fourth-order symplectic integration, Physica D 43, 105 (1990).
  * Celani  _et al._ [2010] A. Celani, S. Musacchio, and D. Vincenzi, Turbulence in More than Two and Less than Three Dimensions, Phys. Rev. Lett. 104, 184506 (2010).
  * Fjørtoft [1953] R. Fjørtoft, On the Changes in the Spectral Distribution of Kinetic Energy for Twodimensional, Nondivergent Flow, Tellus 5, 225 (1953).
  * Nazarenko [2011] S. V. Nazarenko,  _Wave Turbulence_ (Springer-Varlog, Berlin, 2011).
  * Kraichnan [1971] R. H. Kraichnan, Inertial-range transfer in two-and three-dimensional turbulence, J. Fluid Mech. 47, 525 (1971).
  * Majda [2015] A. J. Majda, Statistical energy conservation principle for inhomogeneous turbulent dynamical systems., PNAS 112, 8937 (2015).
  * Verma [2019a] M. K. Verma,  _Energy transfers in Fluid Flows: Multiscale and Spectral Perspectives_ (Cambridge University Press, Cambridge, 2019).
  * Verma  _et al._ [2022] M. K. Verma, S. Chatterjee, A. Sharma, and A. Mohapatra, Equilibrium states of Burgers and Korteweg–de Vries equations, Phys. Rev. E 105, 034121 (2022).
  * Verma [2022] M. K. Verma, Variable energy flux in turbulence, Journal of Physics A: Mathematical and Theoretical 55, 013002 (2022).
  * Landau and Lifshitz [1980] L. D. Landau and E. M. Lifshitz,  _Statistical Physics_ , 3rd ed., Course of Theoretical Physics (Elsevier, Oxford, 1980).
  * Verma [2019b] M. K. Verma, Asymmetric energy transfers in driven nonequilibrium systems and arrow of time, Eur. Phys. J. B 92, 190 (2019b).
  * Kundu  _et al._ [2015] P. K. Kundu, I. M. Cohen, and D. R. Dowling,  _Fluid Mechanics_ , 6th ed. (Academic Press, San Diego, 2015).
  * Choudhuri [2010] A. R. Choudhuri,  _Astrophysics for Physicists_ (Cambridge University Press, Cambridge, 2010).
  * Shannon [1948] C. E. Shannon, A mathematical theory of communication, Bell Labs Tech. J. 27, 379 (1948).
  * Landauer [2000] R. Landauer, Irreversibility and heat generation in the computing process, IBM Journal of Research and Development 44, 261 (2000).
  * Aubry  _et al._ [1991] N. Aubry, R. Guyonnet, and R. Lima, Spatiotemporal analysis of complex signals: Theory and applications, Journal of Statistical Physics 64, 683 (1991).
  * Clark  _et al._ [2020] D. Clark, L. Tarra, and A. Berera, Chaos and information in two-dimensional turbulence, Phys. Rev. Fluids 5, 064608 (2020), 2003.08159 .
  * Drivas and Elgindi [2022] T. D. Drivas and T. M. Elgindi, Singularity formation in the incompressible Euler equation in finite and infinite time, arXiv (2022), 2203.17221 .
  * Chaikin and Lubensky [2000] P. M. Chaikin and T. C. Lubensky,  _Principles of Condensed Matter Physics_ (Cambridge University Press, New York, 2000).
  * Lesieur [2008] M. Lesieur,  _Turbulence in Fluids_ (Springer-Verlag, Dordrecht, 2008).
  * Dar  _et al._ [2001] G. Dar, M. K. Verma, and V. Eswaran, Energy transfer in two-dimensional magnetohydrodynamic turbulence: formalism and numerical results, Physica D 157, 207 (2001).
  * Verma [2004] M. K. Verma, Statistical theory of magnetohydrodynamic turbulence: recent results, Phys. Rep. 401, 229 (2004).


# Light scalar explanation for 18 TeV GRB 221009A

Shyam Balaji<EMAIL_ADDRESS>Laboratoire de Physique Théorique et
Hautes Energies (LPTHE),
UMR 7589 CNRS & Sorbonne Université, 4 Place Jussieu, F-75252, Paris, France
Institut d’Astrophysique de Paris, UMR 7095 CNRS & Sorbonne Université, 98 bis
boulevard Arago, F-75014 Paris, France  Maura E. Ramirez-Quezada
<EMAIL_ADDRESS>Department of Physics, University of
Tokyo, Bunkyo-ku, Tokyo 113–0033, Japan  Dual CP Institute of High Energy
Physics, C.P. 28045, Colima, México  Joseph Silk<EMAIL_ADDRESS>Institut
d’Astrophysique de Paris, UMR 7095 CNRS & Sorbonne Université, 98 bis
boulevard Arago, F-75014 Paris, France  Yongchao Zhang
<EMAIL_ADDRESS>School of Physics, Southeast University, Nanjing
211189, China

###### Abstract

The reecent astrophysical transient Swift J1913.1+1946 may be associated with
the $\gamma$-ray burst GRB 221009A. The redshift of this event is $z\simeq
0.151$. Very high-energy $\gamma$-rays (up to 18 TeV) followed the transient
and were observed by LHAASO; additionally, Carpet-2 detected a photon-like air
shower of 251 TeV. Photons of such high energy are expected to readily
annihilate with the diffuse extragalactic background light (EBL) before
reaching Earth. If the $\gamma$-ray identification and redshift measurements
are correct, new physics could be necessary to explain these measurements.
This letter provides the first CP-even scalar explanation of the most
energetic 18 TeV event reported by LHAASO. In this minimal scenario, the light
scalar singlet $S$ mixes with the Standard Model (SM) Higgs boson $h$. The
highly boosted $S$ particles are produced in the GRB and then undergo the
radiative decay di-photon $S\rightarrow\gamma\gamma$ while propagating to
Earth. The resulting photons may thus be produced at a remote region without
being nullified by the EBL. Hence, the usual exponential reduction of
$\gamma$-rays is lifted due to an attenuation that is inverse in the optical
depth, which becomes much larger due to the scalar carriers.

##  I Introduction

Gamma-ray burst GRB 221009 A events over energies from $500\,\rm GeV$
extending up to $18\,\rm TeV$ have been recently detected. The events are
reported to have occurred at a redshift of $z\approx 0.15$ corresponding to a
distance of $d=645\,\rm Mpc$ de Ugarte Postigo  _et al._ (2022). The initial
detection was first reported by the Swift-Burst Alert Telescope Kennea and
Williams (2022), the Fermi satellite Veres  _et al._ (2022) and CARPET-2
Dzhappuev  _et al._ (2022). The number of photons detected by LHAASO Huang
_et al._ (2022) in an interval of $2000$ s was $\mathcal{O}(5000)$ events.
However, observing such energetic photons is extremely unlikely since the flux
is expected to be rapidly attenuated when propagating and interacting with
extragalactic particlesGould and Schréder (1966); Nikishov (1961); Fazio and
Stecker (1970); Balaji (2022). The unlikeliness for photons in the upper band
of the $500\,{\rm GeV}$-$18\,\rm TeV$ range to survive in this manner while
propagating to the Earth motivated several studies beyond the Standard Model
(SM) as possible solutions. These beyond-the-SM models include explanations
from sterile neutrinos Smirnov and Trautner (2022); Brdar and Li (2022); Huang
_et al._ (2022); Guo  _et al._ (2023), axion-like particlesBaktash  _et al._
(2022); Galanti  _et al._ (2022); Lin and Yanagida (2022); Troitsky (2022);
Nakagawa  _et al._ (2022) and axion-photon conversion scenarios Zhang and Ma
(2022). Even SM explanations have recently been discussed, where external
inverse-Compton interactions with the cosmological radiation fields and proton
synchrotron emissions could generate the ultra-energetic $\gamma$-rays Zhang
_et al._ (2022); Alves Batista (2022); Gonzalez  _et al._ (2022); Sahu  _et
al._ (2023).

In this work, we consider another possible scenario that can explain the high
energy $\gamma$-ray emission. Here, we consider a generic scalar $S$,
corresponding to the inclusion of a gauge singlet to the SM scalar sector.
This new field couples to the SM particles through mixing with the SM Higgs
field $h$, and the interaction is parameterised by the mixing parameter
$\sin\theta$ with angle $\theta$. The singlet extension of the SM scalar
sector is one of the simplest underlying models, and has a host of motivations
including aiding in stabilizing the SM vacuum Gonderinger  _et al._ (2010,
2012); Lebedev (2012); Elias-Miro  _et al._ (2012); Khan and Rakshit (2014);
Falkowski  _et al._ (2015); Ghorbani (2021), addressing the hierarchy problem
in relaxion models Graham  _et al._ (2015); Flacke  _et al._ (2017); Frugiuele
_et al._ (2018); Banerjee  _et al._ (2020); Brax and Burrage (2021),
generating the observed baryon asymmetry in the Universe via baryogenesis
Espinosa and Quiros (1993); Choi and Volkas (1993); Ham  _et al._ (2005);
Profumo  _et al._ (2007); Espinosa  _et al._ (2012); Barger  _et al._ (2012);
Profumo  _et al._ (2015); Curtin  _et al._ (2014); Kotwal  _et al._ (2016);
Chen  _et al._ (2017); Balaji  _et al._ (2021), addressing the cosmological
constant problem through radiative breaking of classical scale invariance Foot
and Kobakhidze (2015); Heikinheimo  _et al._ (2014); Wang  _et al._ (2016),
and mediating the interactions between dark matter (DM) and the SM sector
Pospelov  _et al._ (2008); Baek  _et al._ (2012a, b, 2013); Schmidt-Hoberg
_et al._ (2013); Krnjaic (2016); Beniwal  _et al._ (2016); Balaji and
Kobakhidze (2018).

To accommodate the aforementioned high energy photon observations, the $S$
particles corresponding to this new field are produced in large numbers during
the GRB through hadronic scattering and then undergo the radiative decay
$S\rightarrow\gamma\gamma$ while propagating to Earth. The resulting decay
photons may therefore be produced at a region with less interstellar medium
and avoid being nullified. Hence, the $S$ particles provide an effective means
for the survival of the photons to Earth. More explicitly, this corresponds to
the usual exponential reduction of $\gamma$-rays being lifted due to an
attenuation that is inverse in the optical depth.

In this work, we will use the 18 TeV GRB $\gamma$-ray flux observation and the
number of the detected events in the $2000$ second observation window to set
new astrophysical bounds. These limits complement the existing bounds derived
from stellar cooling and luminosity measurements of the Sun, white dwarfs, red
giants and supernovae Balaji  _et al._ (2022); Dev  _et al._ (2021, 2022).

The outline of this paper is as follows. In Sec. II we introduce the
simplified extension to the SM scalar sector and show the decay rates of the
scalar $S$ into various states. In Sec. III we compute the flux of
$\gamma$-rays produced by the S-decay. Finally, we summarise and discuss the
GRB bounds on the scalar mass and mixing in Sec. IV and conclude in Sec. V.

##  II Scalar decay

This section reviews the decay rates of the CP-even scalar $S$ into
$\gamma$-rays. Therefore, we first introduce the scalar interactions with the
SM particles. At tree-level, the $S$ interacts with leptons and pions via
mixing with the SM Higgs. The interaction Lagrangian for these specific
interactions is given by

$\mathcal{L}=-\sin\theta
S\left[A_{\pi}(\pi^{0}\pi^{0}+\pi^{+}\pi^{-})+\frac{m_{\ell}}{v}\bar{\ell}\ell\right]\,,$
(1)

where $A_{\pi}$ is the effective coupling of $S$ to pions Voloshin (1986);
Donoghue  _et al._ (1990). $\sin\theta$ is the $S$-mixing with the SM Higgs
$h$ and $\ell$ denotes leptons.111Technically $\ell$ represents electrons,
muons and $\tau$ leptons, however in this work, we only consider electrons but
do not consider scalars heavy enough to decay into muons and $\tau$ leptons.
Hence, the scalar can decay into photons at one-loop level as well as
electrons, muons and pions at tree-level from mixing with the Higgs boson.

For $S$ scalars, with masses below twice the electron mass ($m_{S}<2m_{e}$),
the dominant channel is the decay to two photons ($S\to\gamma\gamma$). The
decay rate for this process is given by Dev  _et al._ (2017),

$\Gamma_{S\gamma\gamma}=\frac{121}{9}\frac{\alpha^{2}m_{S}^{3}\sin^{2}\theta}{512\pi^{3}v^{2}}\,,$
(2)

where $\alpha=e^{2}/4\pi$ is the fine structure-constant and $v=246\,\rm GeV$
is the Higgs vacuum expectation value.

When the $S$ scalar masses are above twice the lepton mass threshold
($m_{S}\gtrsim 2m_{\ell}$), the scalar can decay into two oppositely charged
leptons respectively. In this case, the decay rate ($S\to\ell^{+}\ell^{-}$) is
given by Dev  _et al._ (2020a),

$\displaystyle\Gamma_{S\ell^{+}\ell^{-}}$
$\displaystyle=\frac{m_{S}m_{\ell}^{2}\sin^{2}\theta}{8\pi
v^{2}}\left(1-\frac{4m_{\ell}^{2}}{m_{S}^{2}}\right)^{3/2}$ (3)

where $m_{\ell}$ is either the electron or muon mass depending on the decay
channel of interest. Here, we will only consider the scalar massive enough to
decay into a pair of electrons $e$.

##  III Photon flux

We now compute the $\gamma$-ray flux $\Phi_{\gamma}$ from the $S$ scalar
decay. First, we introduce the $\gamma$-ray flux $\Phi_{\gamma}^{0}$ of GRB
221009A, which corresponds to the flux where all photons survive their
interactions with the extragalactic medium on their way to the Earth,

$\displaystyle\Phi_{\gamma}^{0}(E_{\gamma})=$ $\displaystyle 2.1\times
10^{-6}\,{\rm TeV}^{-1}\,{\rm cm}^{-2}\,{\rm s}^{-1}$
$\displaystyle\times\left(\frac{E_{\gamma}}{\rm TeV}\right)^{-1.87\pm 0.04}.$
(4)

This is the so-called unattenuated flux. As discussed in Ref. Smirnov and
Trautner (2022), this is ascertained by extrapolating the flux measured by
Fermi-LAT (GCN 32658) in [0.1,1] GeV to TeV scale Baktash  _et al._ (2022).

We define the decay-production probabilities that a single particle $S$ decays
into two photons, and they reach the Earth. If the $S$ scalar decay occurs at
a distance interval of $[x,x+dx]$, the decay-production probability is then
given by Smirnov and Trautner (2022),

$P_{\rm{decay}}=B_{\gamma}e^{-x/\lambda_{S}}\times\frac{dx}{\lambda_{S}}e^{-(d-x)/\lambda_{\gamma}}.$
(5)

$d$ is the distance at which the flux was produced. $B_{\gamma}$ is the
branching ratio for the scalar decay into photons. The second exponential
factor is the probability that the produced $\gamma$-rays survive their
trajectory and is given in terms of the photon absorption length (or optical
depth) $\tau=d/\lambda_{\gamma}$ where $\lambda_{\gamma}$ is the mean-free-
path of the photon. In this expression, the $\lambda_{S}$ is the mean free
path (MFP) of the $S$ scalar in terms of the $S$ decay rate $\Gamma_{S}$ and
energy of the scalar $E_{S}$ in the Earth’s rest frame

$\lambda_{S}=\frac{E_{S}}{m_{S}\Gamma_{S}}.$ (6)

To compute the $\gamma$-ray flux produced from the $S$ scalar decay, we
integrate the expression given in Eq. (5) over $x$ and multiply by the
incoming $S$ scalar flux at Earth $\Phi_{S}$ from the GRB source, that would
be expected if the scalar did not decay

$\Phi_{\gamma}=\frac{2\Phi_{S}B_{\gamma}}{\lambda_{S}/\lambda_{\gamma}-1}\left(e^{-d/\lambda_{S}}-e^{-d/\lambda_{\gamma}}\right),$
(7)

the factor of 2 is included to account for the two photons produced in the
$S$-decay. Rewriting this expression in terms of the photon absorption length
$\tau(E_{\gamma})$ leads to

$\Phi_{\gamma}=\frac{2\Phi_{S}B_{\gamma}}{\tau\lambda_{S}/d-1}\left(e^{-d/\lambda_{S}}-e^{-\tau}\right).$
(8)

Figure 1: Photon energy flux from GRB 221009A at Earth as a function of photon
energy $E_{\gamma}$ for a scalar flux of $\Phi_{S}=5\times
10^{-23}\,\rm(eVcm^{2}s)^{-1}$ for the cases of 50 events (cyan), 500 events
(magenta) and 5000 events (blue). Flux extrapolated from Fermi-LAT (GCN 32658)
is $\Phi_{\gamma}^{0}$ from Eq. (4) (red-dashed). The energy of the emitted
scalars is taken to be $E_{S}=2\,E_{\gamma}$ TeV

We show the $\gamma$-ray flux from the $S$-decay as a function of the photon
energy $E_{\gamma}$ in Fig. 1, the $S$ scalar energy is $E_{S}=2\,E_{\gamma}$
with a mass of $m_{S}=0.5\,\rm MeV$. The mixing has been fixed by requiring
three different scenarios, one where the scalar decay leads to 5000 events
(blue), 500 events (magenta), and 50 events (cyan). We use the full energy
dependence of $\tau(E_{\gamma})$ from Ref. Domínguez  _et al._ (2011) for the
flux calculation.

Since there are large uncertainties in the hadronic environment and energy
profile of particles at the source, we do not have exact knowledge about the
number of high-energy nucleons that can be produced or preexist in that
region. However, we know that the dominant mode for light scalar production is
nucleon-nucleon bremsstrahlung via pion exchange. Krnjaic (2016); Ishizuka and
Yoshimura (1990); Arndt and Fox (2003); Diener and Burgess (2013); Tu and Ng
(2017); Lee (2018); Dev  _et al._ (2020a, 2022). Possible sources for very
high energy scalars in high nucleon density environments could be via the
interaction of AGN jets with obstacles, such as supernova remnants Torres-Albà
(2020). The collision of the supernova ejecta with the high-energy jets could
lead to the formation of an interaction region, in which particles can be
accelerated and produce high-energy scalar emission Vieyro  _et al._ (2019).
Scalars can also be produced by heavy meson decays. With present information
about the GRB, the exact fraction of energy that can be extracted from the GRB
in exotic states is unknown. However we can make some conservative estimates
and adjust them in them in the event of more refined measurements in the
future. The amount of energy emitted in a supernova explosion is typically
about $10^{44}\,\textrm{J}\approx 6.2\times 10^{62}$ eV. In general, almost
all of this energy is carried away by neutrinos. In the case of the GRB of
interest here, the energy output is estimated at $2\times
10^{54}\,\textrm{erg}=1.2\times 10^{66}$ eV Cheung (2022). Assuming that
almost all the energy is lost and emitted as neutrinos, and using 0.3% of this
energy as an upper bound for scalar emission with a rectangular energy
distribution in the range $[1,36]$ TeV corresponding to an average energy of
$\overline{E}_{S}=18.5$ TeV 222Note that the rectangular spectrum in the
energy of the scalar flux at the GRB source between $[1,36]$ TeV is chosen
because this results in photon energies of $E_{\gamma}\simeq E_{S}/2$ which
lies in the $[0.5,18]$ TeV range LHAASO observations window. In the absence of
more information about the GRB source energy profile, this is the simplest
possible choice. , we get a scalar flux (assuming no decays) of
$\Phi_{S}\simeq 5\times
10^{-23}\,\textrm{eV}^{-1}\,\textrm{cm}^{-2}\,\textrm{s}^{-1}$ at Earth over
the 2000 s duration. It should be noted that this is very conservative since
when considering supernova bounds, typically, the flux of the new particle is
bounded by 10% of the neutrino luminosity, see, e.g. Refs. Balaji  _et al._
(2022); Arndt and Fox (2003); Lee (2018). It should be noted that
$E_{S}=2\,E_{\gamma}$ is fixed simply due to kinematics of two-body decays.
Since the max energy for the GRB photons is observed at 18 TeV, the energy of
the parent scalar has to be at least twice this energy. Other, more exotic
beyond the SM scenarios can also be invoked to generate the scalar flux from
the source, but we do not explore them in this work.

Finally, we can observe in Fig. 1 that the flux of $\gamma$-rays from the $S$
scalar decay increases at high energies, approaching the unattenuated
$\gamma$-ray flux. This indicates that the photons can survive interactions
with the extragalactic background.

##  IV Bounds on the scalar mass-mixing parameter space

The flux in the previous section provides us with a means to set bounds on the
scalar mass-mixing parameter space. We require the number of photons produced
by the scalar decay to be $N_{\gamma}^{S}=50,\,500$ and 5000 respectively.
Here $N_{\gamma}^{S}$ is obtained by integrating Eq. (8) over the photon
energy $E_{\gamma}$ over the range of [$0.5,18$] TeV as well as integrating
over the time $t$ from 0 s to 2000 s and the detector area of $1\,\rm{km}^{2}$
Cui  _et al._ (2014); Ma  _et al._ (2022).

We consider scalars with masses up to 2 MeV. Therefore, these scalars can only
decay directly into photons and electrons. It is easy to find and extrapolate
the relationship between $\sin\theta$ and $m_{S}$ as a function of the number
of events expected in 2000 s. This mixing-mass relationship providing upper
bounds is found to be given by the relation

$\sin\theta\simeq\left(\frac{m_{S}^{2}}{1.57\times
10^{-8}\,\rm{eV^{2}}}\right)^{-1}\sqrt{\frac{N_{\gamma}^{S}}{\Phi_{S}/(\rm{eV\,cm}^{2}\,\rm{s})^{-1}}}.$
(9)

Figure 2: Limits for light scalar decays into photons to explain GRB 221009 A
in the scalar mass-mixing ($m_{S},\sin\theta$) plane requiring
$N_{\gamma}^{S}\geq 50,\,500$ and 5000 number of events. The shaded regions
correspond to supernova exclusions from Ref. Balaji  _et al._ (2022) from
various models such as Fischer 11$M_{\odot}$ (blue), Fischer 18$M_{\odot}$
(cyan) and Nakazato 13$M_{\odot}$ (magenta). Additionally, we add
complementary limits from SN1987A and collider searches Dev  _et al._ (2017);
Egana-Ugrinovic  _et al._ (2020); Dev  _et al._ (2020b).

This relationship also depends on the scalar flux $\Phi_{S}$ and assumes an
even distribution of $\Phi_{S}$ over the $E_{S}$ range $[1,36]$ TeV.

In Fig. 2, we show the upper bounds for light scalar decaying into photons to
explain GRB 221009 A in the scalar mixing-mass plane. This was computed
numerically using the expression in Eq. (8) considering three benchmarks for
the number of events $N_{\gamma}^{S}=50,\,500$ and 5000. We also show the
stellar exclusions in Ref. Balaji  _et al._ (2022) given by the supernova
SN1987A using different numerical profiles Fischer 11.8M⊙ (blue), Fischer 18M⊙
(cyan) Fischer  _et al._ (2016) and Nakazato 13M⊙ (magenta) Nakazato  _et al._
(2013). We also show the combined collider searches (gray) Dev  _et al._
(2017); Egana-Ugrinovic  _et al._ (2020); Dev  _et al._ (2020b). Note that
when there is an order of magnitude difference between the expected number of
events, the difference between the bounds is a factor of a few. For the case
of two orders of magnitude difference in $N_{\gamma}^{S}$, we find an order of
magnitude difference between the limits. This is easily understood from the
proportionality between the mixing and the number of events expected scales
like $\sin\theta\propto\sqrt{N_{\gamma}^{S}}$. Additionally, we check the MFP
of the scalar at an average energy of $\overline{E}_{S}=18.5$ TeV for the
three different expected event boundaries shown in Fig. 2, these are on the
order of $\simeq 10^{30}\,\rm cm$ (50 events), $\simeq 10^{29}\,\rm cm$ (500
events) and $\simeq 4\times 10^{26}$-$5\times 10^{27}\,\rm cm$ (5000 events)
which is comparable to the estimated GRB source distance of $\simeq 10^{27}$
cm.

The bounds obtained from GRB 221009 A by requiring $N_{\gamma}^{S}=50,\,500$
and $5000$ events complement those given by supernovae Balaji  _et al._ (2022)
and collider exclusions Dev  _et al._ (2017); Egana-Ugrinovic  _et al._
(2020); Dev  _et al._ (2020b). Note that considering regions with $>500$ or
$>50$ events, there is no upper bound on mixing. However, there is a small
band in Fig. 2 for $>5000$ events because, in this region, $\lambda_{S}\simeq
d$ and hence the photon flux is locally maximised. This can be understood as
follows: At higher mixing, the scalars decay closer to the source, and more of
the resulting photons get attenuated by the EBL. However, even if the scalars
decay very quickly near the source, some minimum number of lower energy
photons in the $[0.5,18]$ TeV range will still survive their journey to Earth,
albeit with a smaller number of events, as shown. The lower bounds on mixing
are because the scalars are too long-lived. On average, they may decay well
beyond the Earth, reducing the observed photon flux and hence the event yield.

While supernova bounds excluded parameter space across the whole scalar mass
range between mixing of the order of $10^{-7}$-$10^{-5}$, the excluded area
from GRB 221009 A intersects the excluded area from the supernovae limits, at
masses between $\simeq{\rm a\,few}\,10$-$10^{3}\,\rm keV$, in the same mixing
interval. However, GRB 221009 A excludes a large parameter of space above a
mixing of $10^{-5}$. It is important to notice that the GRB 221009 A gives a
better constraint than the supernova for a scalar with a mass of
$m_{S}=10^{3}$ keV at a fixed mixing of $\sin\theta\simeq 10^{-8}$.

The time delay $\Delta t$ in arrival (or time of flight difference) to travel
a distance $d$ between a massive particle of mass $m_{S}$ and a massless
particle is given by the simple formula $\Delta
t=\frac{m_{S}^{2}}{2E_{S}^{2}}d$. We require that the time delay is smaller
than the GRB duration of 2000s. We can see from this that the largest time
delay is controlled by the higher mass scalars with lower total energy. To
approximate the time delay, we consider the limiting mass for an average
energy in the scalar spectrum of $\overline{E}_{S}$, which yields $m_{S}\simeq
4.5\times 10^{3}$ keV, well outside our allowed region in Fig. 2 which is cut-
off at $10^{3}$ keV. Hence the allowed region is not affected by the time
delay and complies easily with the 2000s requirement. If, however, additional
information arises enabling the identification of the lowest photon energy in
the range as low (or lower) than 0.5 TeV from $S$ decays, this would
strengthen the bound to lower masses, although such an identification is
unlikely given the large background from conventional photons in this region.
Conversely, if a very high energy photon is identified with a long time delay
event, this would indicate higher $m_{S}$. For massless particles in the final
state produced by the decay of a massive boosted scalar, like in our scenario,
the angular dispersion $\Theta$ in extremely relativistic limit in the lab
frame is given in general by $\Theta=\frac{2m_{S}}{E_{S}}$. Taking the maximum
mass of interest as a limiting scenario and the average scalar energy, we get
$\Theta\simeq 10^{-7}$. Hence, as long as the GRB jet opening angle exceeds
this tiny number, there is no additional suppression of the $\gamma$-rays from
$S$ at the Earth.

##  V Conclusion

In this work, we have considered a generic scalar CP-even $S$ scalar
corresponding to a CP-even singlet extension of the SM scalar sector as an
explanation for GRB 221009 A. We consider a scenario wherein the $S$ particles
are produced in the GRB through hadronic scattering and then undergo the
radiative loop-level decay $S\rightarrow\gamma\gamma$ while propagating to
Earth. If the photons are produced without being nullified by the
extragalactic background, it could explain the anomalously high observed
$\gamma-$ray emission at 18 TeV.

We calculated the flux of $\gamma$-rays from this process for different
benchmarks of the expected number of events on Earth. We found that the
$\gamma$-ray flux produced increases at higher photon energies and approaches
the unattenuated $\gamma$-ray flux.

We have also computed bounds on the mixing-mass parameter space of the $S$
scalar for different fractions of the observed LHAASO events. We found that at
masses between $\simeq{\rm a\,few}10$-$10^{3}\,\rm keV$, it overlaps with the
excluded parameter space set by the SN1987A luminosity bounds. However, an
upper bound exists in the parameter space of large mixing for all masses
between $\simeq 10^{1}$-$10^{3}\,\rm keV$ that is not excluded by SN1987A or
collider measurements. We also found that the strongest $S$ constraint for GRB
221009 A is for a scalar with a mass of $10^{3}$ keV and a mixing of $\simeq
10^{-8}$ for 50 events.

##  VI Acknowledgements

S.B. would like to thank Andrea Caputo for helpful discussion. S.B. is
supported by funding from the European Union’s Horizon 2020 research and
innovation programme under grant agreement No. 101002846 (ERC CoG
“CosmoChart”) as well as support from the Initiative Physique des Infinis
(IPI), a research training program of the Idex SUPER at Sorbonne Université.
M.R.Q. is supported by the JSPS KAKENHI Grant Number 20H01897. Y.Z. is
supported by the National Natural Science Foundation of China under grant No.
12175039, the 2021 Jiangsu Shuangchuang (Mass Innovation and Entrepreneurship)
Talent Program No. JSSCBS20210144, and the “Fundamental Research Funds for the
Central Universities.”

## References

  * de Ugarte Postigo  _et al._ (2022) A. de Ugarte Postigo, L. Izzo, G. Pugliese, D. Xu, B. Schneider, J. P. U. Fynbo, N. R. Tanvir, D. B. Malesani, A. Saccardi, D. A. Kann, K. Wiersema, B. P. Gompertz, C. C. Thoene, A. J. Levan, and Stargate Collaboration, GRB Coordinates Network 32648, 1 (2022).
  * Kennea and Williams (2022) J. A. Kennea and M. Williams, GCN CIRCULAR 32635, 1 (2022).
  * Veres  _et al._ (2022) P. Veres, E. Burns, E. Bissaldi, S. Lesage, O. Roberts, and Fermi GBM Team, GRB Coordinates Network 32636, 1 (2022).
  * Dzhappuev  _et al._ (2022) D. D. Dzhappuev, Y. Z. Afashokov, I. M. Dzaparova, T. A. Dzhatdoev, E. A. Gorbacheva, I. S. Karpikov, M. M. Khadzhiev, N. F. Klimenko, A. U. Kudzhaev, A. N. Kurenya, A. S. Lidvansky, O. I. Mikhailova, V. B. Petkov, E. I. Podlesnyi, N. A. Pozdnukhov, V. S. Romanenko, G. I. Rubtsov, S. V. Troitsky, I. B. Unatlokov, I. A. Vaiman, A. F. Yanin, and K. V. Zhuravleva, The Astronomer’s Telegram 15669, 1 (2022).
  * Huang  _et al._ (2022) Y. Huang, S. Hu, S. Chen, M. Zha, C. Liu, Z. Yao, Z. Cao, and T. L. Experiment, GRB Coordinates Network 32677, 1 (2022).
  * Gould and Schréder (1966) R. J. Gould and G. Schréder, Phys. Rev. Lett. 16, 252 (1966).
  * Nikishov (1961) A. I. Nikishov, Zhur. Eksptl’. i Teoret. Fiz. Vol: 41 (1961).
  * Fazio and Stecker (1970) G. G. Fazio and F. W. Stecker, Nature 226, 135 (1970).
  * Balaji (2022) S. Balaji, (2022), arXiv:2211.03807 [astro-ph.HE] .
  * Smirnov and Trautner (2022) A. Y. Smirnov and A. Trautner, (2022), arXiv:2211.00634 [hep-ph] .
  * Brdar and Li (2022) V. Brdar and Y.-Y. Li, (2022), arXiv:2211.02028 [hep-ph] .
  * Huang  _et al._ (2022) J. Huang, Y. Wang, B. Yu, and S. Zhou, (2022), arXiv:2212.03477 [hep-ph] .
  * Guo  _et al._ (2023) S.-Y. Guo, M. Khlopov, L. Wu, and B. Zhu, (2023), arXiv:2301.03523 [hep-ph] .
  * Baktash  _et al._ (2022) A. Baktash, D. Horns, and M. Meyer, (2022), arXiv:2210.07172 [astro-ph.HE] .
  * Galanti  _et al._ (2022) G. Galanti, M. Roncadelli, and F. Tavecchio, (2022), arXiv:2210.05659 [astro-ph.HE] .
  * Lin and Yanagida (2022) W. Lin and T. T. Yanagida, (2022), arXiv:2210.08841 [hep-ph] .
  * Troitsky (2022) S. V. Troitsky, (2022), arXiv:2210.09250 [astro-ph.HE] .
  * Nakagawa  _et al._ (2022) S. Nakagawa, F. Takahashi, M. Yamada, and W. Yin, (2022), arXiv:2210.10022 [hep-ph] .
  * Zhang and Ma (2022) G. Zhang and B.-Q. Ma, (2022), arXiv:2210.13120 [hep-ph] .
  * Zhang  _et al._ (2022) B. T. Zhang, K. Murase, K. Ioka, D. Song, C. Yuan, and P. Mészáros, (2022), arXiv:2211.05754 [astro-ph.HE] .
  * Alves Batista (2022) R. Alves Batista, (2022), arXiv:2210.12855 [astro-ph.HE] .
  * Gonzalez  _et al._ (2022) M. M. Gonzalez, D. A. Rojas, A. Pratts, S. Hernandez-Cadena, N. Fraija, R. Alfaro, Y. P. Araujo, and J. A. Montes, (2022), arXiv:2210.15857 [astro-ph.HE] .
  * Sahu  _et al._ (2023) S. Sahu, B. Medina-Carrillo, G. Sánchez-Colón, and S. Rajpoot, Astrophys. J. Lett. 942, L30 (2023), arXiv:2211.04057 [astro-ph.HE] .
  * Gonderinger  _et al._ (2010) M. Gonderinger, Y. Li, H. Patel, and M. J. Ramsey-Musolf, JHEP 01, 053 (2010), arXiv:0910.3167 [hep-ph] .
  * Gonderinger  _et al._ (2012) M. Gonderinger, H. Lim, and M. J. Ramsey-Musolf, Phys. Rev. D 86, 043511 (2012), arXiv:1202.1316 [hep-ph] .
  * Lebedev (2012) O. Lebedev, Eur. Phys. J. C 72, 2058 (2012), arXiv:1203.0156 [hep-ph] .
  * Elias-Miro  _et al._ (2012) J. Elias-Miro, J. R. Espinosa, G. F. Giudice, H. M. Lee, and A. Strumia, JHEP 06, 031 (2012), arXiv:1203.0237 [hep-ph] .
  * Khan and Rakshit (2014) N. Khan and S. Rakshit, Phys. Rev. D 90, 113008 (2014), arXiv:1407.6015 [hep-ph] .
  * Falkowski  _et al._ (2015) A. Falkowski, C. Gross, and O. Lebedev, JHEP 05, 057 (2015), arXiv:1502.01361 [hep-ph] .
  * Ghorbani (2021) P. Ghorbani, Nucl. Phys. B 971, 115533 (2021), arXiv:2104.09542 [hep-ph] .
  * Graham  _et al._ (2015) P. W. Graham, D. E. Kaplan, and S. Rajendran, Phys. Rev. Lett. 115, 221801 (2015), arXiv:1504.07551 [hep-ph] .
  * Flacke  _et al._ (2017) T. Flacke, C. Frugiuele, E. Fuchs, R. S. Gupta, and G. Perez, JHEP 06, 050 (2017), arXiv:1610.02025 [hep-ph] .
  * Frugiuele  _et al._ (2018) C. Frugiuele, E. Fuchs, G. Perez, and M. Schlaffer, JHEP 10, 151 (2018), arXiv:1807.10842 [hep-ph] .
  * Banerjee  _et al._ (2020) A. Banerjee, H. Kim, O. Matsedonskyi, G. Perez, and M. S. Safronova, JHEP 07, 153 (2020), arXiv:2004.02899 [hep-ph] .
  * Brax and Burrage (2021) P. Brax and C. Burrage, Phys. Rev. D 104, 015011 (2021), arXiv:2101.10693 [hep-ph] .
  * Espinosa and Quiros (1993) J. R. Espinosa and M. Quiros, Phys. Lett. B 305, 98 (1993), arXiv:hep-ph/9301285 .
  * Choi and Volkas (1993) J. Choi and R. R. Volkas, Phys. Lett. B 317, 385 (1993), arXiv:hep-ph/9308234 .
  * Ham  _et al._ (2005) S. W. Ham, Y. S. Jeong, and S. K. Oh, J. Phys. G 31, 857 (2005), arXiv:hep-ph/0411352 .
  * Profumo  _et al._ (2007) S. Profumo, M. J. Ramsey-Musolf, and G. Shaughnessy, JHEP 08, 010 (2007), arXiv:0705.2425 [hep-ph] .
  * Espinosa  _et al._ (2012) J. R. Espinosa, T. Konstandin, and F. Riva, Nucl. Phys. B 854, 592 (2012), arXiv:1107.5441 [hep-ph] .
  * Barger  _et al._ (2012) V. Barger, D. J. H. Chung, A. J. Long, and L.-T. Wang, Phys. Lett. B 710, 1 (2012), arXiv:1112.5460 [hep-ph] .
  * Profumo  _et al._ (2015) S. Profumo, M. J. Ramsey-Musolf, C. L. Wainwright, and P. Winslow, Phys. Rev. D 91, 035018 (2015), arXiv:1407.5342 [hep-ph] .
  * Curtin  _et al._ (2014) D. Curtin, P. Meade, and C.-T. Yu, JHEP 11, 127 (2014), arXiv:1409.0005 [hep-ph] .
  * Kotwal  _et al._ (2016) A. V. Kotwal, M. J. Ramsey-Musolf, J. M. No, and P. Winslow, Phys. Rev. D 94, 035022 (2016), arXiv:1605.06123 [hep-ph] .
  * Chen  _et al._ (2017) C.-Y. Chen, J. Kozaczuk, and I. M. Lewis, JHEP 08, 096 (2017), arXiv:1704.05844 [hep-ph] .
  * Balaji  _et al._ (2021) S. Balaji, M. Spannowsky, and C. Tamarit, JCAP 03, 051 (2021), arXiv:2010.08013 [hep-ph] .
  * Foot and Kobakhidze (2015) R. Foot and A. Kobakhidze, Int. J. Mod. Phys. A 30, 1550126 (2015), arXiv:1112.0607 [hep-ph] .
  * Heikinheimo  _et al._ (2014) M. Heikinheimo, A. Racioppi, M. Raidal, C. Spethmann, and K. Tuominen, Mod. Phys. Lett. A 29, 1450077 (2014), arXiv:1304.7006 [hep-ph] .
  * Wang  _et al._ (2016) Z.-W. Wang, T. G. Steele, T. Hanif, and R. B. Mann, JHEP 08, 065 (2016), arXiv:1510.04321 [hep-ph] .
  * Pospelov  _et al._ (2008) M. Pospelov, A. Ritz, and M. B. Voloshin, Phys. Lett. B 662, 53 (2008), arXiv:0711.4866 [hep-ph] .
  * Baek  _et al._ (2012a) S. Baek, P. Ko, and W.-I. Park, JHEP 02, 047 (2012a), arXiv:1112.1847 [hep-ph] .
  * Baek  _et al._ (2012b) S. Baek, P. Ko, W.-I. Park, and E. Senaha, JHEP 11, 116 (2012b), arXiv:1209.4163 [hep-ph] .
  * Baek  _et al._ (2013) S. Baek, P. Ko, W.-I. Park, and E. Senaha, JHEP 05, 036 (2013), arXiv:1212.2131 [hep-ph] .
  * Schmidt-Hoberg  _et al._ (2013) K. Schmidt-Hoberg, F. Staub, and M. W. Winkler, Phys. Lett. B 727, 506 (2013), arXiv:1310.6752 [hep-ph] .
  * Krnjaic (2016) G. Krnjaic, Phys. Rev. D 94, 073009 (2016), arXiv:1512.04119 [hep-ph] .
  * Beniwal  _et al._ (2016) A. Beniwal, F. Rajec, C. Savage, P. Scott, C. Weniger, M. White, and A. G. Williams, Phys. Rev. D 93, 115016 (2016), arXiv:1512.06458 [hep-ph] .
  * Balaji and Kobakhidze (2018) S. Balaji and A. Kobakhidze, (2018), arXiv:1812.10914 [hep-ph] .
  * Balaji  _et al._ (2022) S. Balaji, P. S. B. Dev, J. Silk, and Y. Zhang, (2022), arXiv:2205.01669 [hep-ph] .
  * Dev  _et al._ (2021) P. S. B. Dev, R. N. Mohapatra, and Y. Zhang, JCAP 05, 014 (2021), arXiv:2010.01124 [hep-ph] .
  * Dev  _et al._ (2022) P. S. B. Dev, J.-F. Fortin, S. P. Harris, K. Sinha, and Y. Zhang, JCAP 01, 006 (2022), arXiv:2111.05852 [hep-ph] .
  * Voloshin (1986) M. B. Voloshin, Sov. J. Nucl. Phys. 44, 478 (1986).
  * Donoghue  _et al._ (1990) J. F. Donoghue, J. Gasser, and H. Leutwyler, Nucl. Phys. B 343, 341 (1990).
  * Dev  _et al._ (2017) P. S. B. Dev, R. N. Mohapatra, and Y. Zhang, Nucl. Phys. B 923, 179 (2017), arXiv:1703.02471 [hep-ph] .
  * Dev  _et al._ (2020a) P. S. B. Dev, R. N. Mohapatra, and Y. Zhang, JCAP 08, 003 (2020a), [Erratum: JCAP 11, E01 (2020)], arXiv:2005.00490 [hep-ph] .
  * Domínguez  _et al._ (2011) A. Domínguez, J. R. Primack, D. J. Rosario, F. Prada, R. C. Gilmore, S. M. Faber, D. C. Koo, R. S. Somerville, M. A. Pérez-Torres, P. Pérez-González, J. S. Huang, M. Davis, P. Guhathakurta, P. Barmby, C. J. Conselice, M. Lozano, J. A. Newman, and M. C. Cooper, Mon. Not. Roy. Astron. Soc. 410, 2556 (2011), arXiv:1007.1459 [astro-ph.CO] .
  * Ishizuka and Yoshimura (1990) N. Ishizuka and M. Yoshimura, Prog. Theor. Phys. 84, 233 (1990).
  * Arndt and Fox (2003) D. Arndt and P. J. Fox, JHEP 02, 036 (2003), arXiv:hep-ph/0207098 .
  * Diener and Burgess (2013) R. Diener and C. P. Burgess, JHEP 05, 078 (2013), arXiv:1302.6486 [hep-ph] .
  * Tu and Ng (2017) H. Tu and K.-W. Ng, JHEP 07, 108 (2017), arXiv:1706.08340 [hep-ph] .
  * Lee (2018) J. S. Lee, (2018), arXiv:1808.10136 [hep-ph] .
  * Torres-Albà (2020) N. Torres-Albà, PoS HEPROVII, 005 (2020), arXiv:1911.00299 [astro-ph.HE] .
  * Vieyro  _et al._ (2019) F. L. Vieyro, V. Bosch-Ramon, and N. Torres-Albà, Astron. Astrophys. 622, A175 (2019), arXiv:1901.09003 [astro-ph.HE] .
  * Cheung (2022) K. Cheung, (2022), arXiv:2210.14178 [hep-ph] .
  * Cui  _et al._ (2014) S. Cui, Y. Liu, Y. Liu, and X. Ma, Astroparticle Physics 54, 86 (2014).
  * Ma  _et al._ (2022) X.-H. Ma, Y.-J. Bi, Z. Cao, M.-J. Chen, S.-Z. Chen, Y.-D. Cheng, G.-H. Gong, M.-H. Gu, H.-H. He, C. Hou, W.-H. Huang, X.-T. Huang, C. Liu, O. Shchegolev, X.-D. Sheng, Y. Stenkin, C.-Y. Wu, H.-R. Wu, S. Wu, G. Xiao, Z.-G. Yao, S.-S. Zhang, Y. Zhang, and X. Zuo, Chinese Physics C 46, 030001 (2022).
  * Egana-Ugrinovic  _et al._ (2020) D. Egana-Ugrinovic, S. Homiller, and P. Meade, Phys. Rev. Lett. 124, 191801 (2020), arXiv:1911.10203 [hep-ph] .
  * Dev  _et al._ (2020b) P. S. B. Dev, R. N. Mohapatra, and Y. Zhang, Phys. Rev. D 101, 075014 (2020b), arXiv:1911.12334 [hep-ph] .
  * Fischer  _et al._ (2016) T. Fischer, S. Chakraborty, M. Giannotti, A. Mirizzi, A. Payez, and A. Ringwald, Phys. Rev. D 94, 085012 (2016), arXiv:1605.08780 [astro-ph.HE] .
  * Nakazato  _et al._ (2013) K. Nakazato, K. Sumiyoshi, H. Suzuki, T. Totani, H. Umeda, and S. Yamada, Astrophys. J. Suppl. 205, 2 (2013), arXiv:1210.6841 [astro-ph.HE] .