# Multi-messenger Implications of Sub-PeV Diffuse Galactic Gamma-Ray Emission

Ke Fang  Department of Physics, Wisconsin IceCube Particle Astrophysics
Center, University of Wisconsin, Madison, WI, 53706  Kohta Murase  Department
of Physics, The Pennsylvania State University, University Park, Pennsylvania
16802, USA  Department of Astronomy & Astrophysics, The Pennsylvania State
University, University Park, Pennsylvania 16802, USA  Center for
Multimessenger Astrophysics, Institute for Gravitation and the Cosmos, The
Pennsylvania State University, University Park, Pennsylvania 16802, USA
Center for Gravitational Physics, Yukawa Institute for Theoretical Physics,
Kyoto, Kyoto 606-8502 Japan

###### Abstract

The diffuse Galactic gamma-ray flux between 0.1 and 1 PeV has recently been
measured by the Tibet AS$\gamma$ Collaboration. The flux and spectrum are
consistent with the decay of neutral pions from hadronuclear interactions
between Galactic cosmic rays and the interstellar medium (ISM). We derive the
flux of the Galactic diffuse neutrino emission from the same interaction
process that produces the gamma rays. Our calculation accounts for the effect
of gamma-ray attenuation inside the Milky Way and uncertainties due to the
spectrum and distribution of cosmic rays, gas density, and infrared emission
of the ISM. We find that the contribution from the Galactic plane to the all-
sky neutrino flux is $\lesssim 5-10$% around 100 TeV. The Galactic and
extragalactic neutrino intensities are comparable in the Galactic plane
region. Our results are consistent with the upper limit reported by the
IceCube and ANTARES Collaborations, and predict that next-generation neutrino
experiments may observe the Galactic component. We also show that the Tibet
AS$\gamma$ data imply either an additional component in the cosmic-ray nucleon
spectrum or contribution from discrete sources, including Pevatrons such as
superbubbles and hypernova remnants, and PeV electron accelerators. Future
multi-messenger observations between 1 TeV and 1 PeV are crucial to
decomposing the origin of sub-PeV gamma rays.

Galactic cosmic rays (567), Neutrino astronomy (1100), High-energy cosmic
radiation (731)

##  1 Introduction

Diffuse gamma rays with energies between 100 TeV and 1 PeV have recently been
detected by the Tibet air-shower gamma (Tibet AS$\gamma$) experiment. The
arrival directions of the sub-PeV gamma rays extend over the Galactic plane
(GP) and are consistent with being from diffuse cosmic rays in the Galaxy
(Amenomori et al., 2021). The origin of cosmic rays has been an enigma
especially around and beyond the knee energy at $E_{\rm knee}\sim 3-4$ PeV.
Various sources, including superbubbles, hypernovae, supernovae in dense
circumstellar material, merger remnants, pulsar wind nebulae, and the Galactic
center, have been proposed as potential contributors. Dissecting the diffuse
gamma-ray emission at MeV-GeV energies (Hunter et al., 1997; Ackermann et al.,
2012a) to PeV energies is important to understanding the cosmic-ray sources
and propagation physics.

Sub-PeV gamma rays play crucial roles in multi-messenger astrophysics. The
detection of $\sim 0.1-1$ PeV neutrinos enables us to probe $\sim 3-30$ PeV
protons beyond the knee energy. The IceCube experiment has measured diffuse
TeV-PeV neutrinos that are consistent with an extragalactic origin (Aartsen et
al., 2013a, b, 2015, 2020; Stettner, 2019). The neutrino and gamma-ray
connection has been shown to be powerful for revealing both Galactic and
extragalactic sources of neutrinos (Murase et al., 2013; Ahlers & Murase,
2014; Murase et al., 2016).

To date diffuse Galactic neutrino emission has not been discovered (Aartsen et
al., 2017; Albert et al., 2018), despite that it has been predicted to exist
for decades (Stecker, 1979; Kheirandish, 2020, for a recent review) and
studied in light of the IceCube measurements (Ahlers & Murase, 2014;
Anchordoqui et al., 2014; Neronov et al., 2014; Joshi et al., 2013; Kachelrieß
& Ostapchenko, 2014; Spurio, 2014; Gaggero et al., 2015; Palladino & Vissani,
2016; Denton et al., 2017). In particular, primarily based on previous sub-PeV
gamma-ray limits posed by the CASA-MIA (Borione et al., 1998) and KASCADE
experiments, Ahlers & Murase (2014) showed that the Galactic contribution to
IceCube neutrinos is subdominant (see also Murase et al., 2016), and the GP
may give $\sim 3-10$% of the $10-100$ TeV all-sky neutrino flux with
$E_{\nu}^{2}\Phi_{\nu}^{\rm IC}\sim(5-10)\times{10}^{-8}~{}{\rm GeV}~{}{\rm
cm}^{-2}~{}{\rm s}^{-1}~{}{\rm sr}^{-1}$ (Aartsen et al., 2015, 2020). This is
also consistent with neutrino constraints (e.g., Spurio, 2014; Ahlers et al.,
2016) as well as the latest IceCube and ANTARES results, $\lesssim 8.5$ %
(Albert et al., 2018).

The Galactic disk is magnetized, and can be regarded as a cosmic-ray
“reservoir” with a typical escape time scale of $t_{\rm esc}\sim 30~{}{\rm
Myr}~{}{(R/1~{}{\rm GV})}^{-\delta}$ (e.g., Murase & Fukugita, 2019), where
$R$ is the rigidity and $\delta\sim 0.3-0.5$. Pionic gamma rays and neutrinos
are co-produced when cosmic-ray ions interact with gas and dust particles in
the Galaxy (Hayakawa, 1952; Stecker, 1977). The detection of diffuse Galactic
gamma rays therefore provides a solid reference to the flux level of the
diffuse Galactic neutrinos as well as the cosmic-ray confinement and injection
around the knee energy. In this work, by taking into account details of the
gamma-ray attenuation along the line of sight, we evaluate the diffuse
Galactic neutrino flux using the sub-PeV gamma-ray flux observed by the Tibet
AS$\gamma$ experiment. We show that the Galactic contribution can be
comparable to the extragalactic flux in the direction to the GP especially
below $E_{\nu}\lesssim 100$ TeV, and may be observed by the next-generation
neutrino telescopes. Our conclusion also applies to neutrino emission from
unresolved sources in the Galaxy.

This letter is organized as follows. We study the diffuse neutrino emission in
$\S\ref{sec:diffuse}$. We discuss the potential contribution to sub-PeV gamma
rays from discrete sources, including Cygnus Cocoon
($\S\ref{sec:CygnusCocoon}$) and unresolved hadronic ($\S\ref{sec:HNR}$) and
leptonic ($\S\ref{sec:leptonic}$) sources. We conclude in
$\S\ref{sec:discussion}$.

##  2 Multimessenger Connection in Galactic Diffuse Emission

In hadronuclear ($pp$) scenarios, the differential spectrum of neutrinos and
gamma rays per interaction are related by (e.g., Murase et al., 2013; Ahlers &
Murase, 2014),

$E_{\nu}^{2}\frac{dN_{\nu}}{dE_{\nu}}\approx\frac{3}{2}\left.\left(E_{\gamma}^{2}\frac{dN_{\gamma}}{dE_{\gamma}}\right)\right|_{E_{\nu}\approx
E_{\gamma}/2},$ (1)

considering that the ratio of charged and neutral pions is approximately $2:1$
at high energies and each neutrino carries $\sim 1/4$ of the pion energy.
Neutrinos with characteristic energy carry $\sim 3-5$% of the parent nucleon
energy, i.e., $E_{\nu}\sim(0.03-0.05)E_{p}$. Equation (1) is further subject
to the gamma-ray attenuation. This is because PeV gamma rays travel only $\sim
10$ kpc. In this work, we take into account the effect of the gamma-ray
absorption in detail. We do not consider electromagnetic cascades since the
effect is small for steep spectra with $\alpha_{\nu}>2$ (see, e.g., Murase &
Beacom, 2012).

###  2.1 Methods

Below $\vec{x}$ indicates a cylindrical coordinate system with the Galactic
center (GC) at the origin and the GP on the $xy$ plane, referred as the GC
frame. The projected distance of $\vec{x}$ to the GC on the $xy$ plane is
noted as $r$ and the distance to the plane is $z$. This coordinate system is
suitable for the description of the diffuse infrared emission and cosmic-ray
source distribution in the Galaxy, which are approximately cylindrically
symmetric. For observation of neutrinos and gamma rays at the solar
neighborhood we will use the Galactic coordinate $\vec{x}_{g}$. Taking the
direction toward the Galactic Center as the $x$-axis, we have
$\vec{x}_{g}=(s\cos b\,\cos l,s\cos b\,\sin l,s\,\sin b)$, where $s$ is the
distance to the observer (the Sun), $l$ and $b$ are the Galactic longitude and
latitude, respectively. $\vec{x}_{g}$ and $\vec{x}$ are converted by
$\vec{x}=\vec{x}_{g}+\vec{x}_{\rm obs}$, where $\vec{x}_{\rm obs}$ is the
coordinate of the observer in the GC frame.

The inverse of the mean free path for a gamma ray of energy $E_{\gamma}$ and
direction $\hat{u}$ at a space point $\vec{x}$ is,

$\lambda_{\gamma\gamma}^{-1}(E_{\gamma},\hat{u},\vec{x})=\int
d\Omega(1-\hat{u}\cdot\hat{k})\int d\varepsilon\frac{dn}{d\varepsilon
d\Omega}(\vec{x})\sigma_{\gamma\gamma}(E_{\gamma}\varepsilon(1-\hat{u}\cdot\hat{k})).$
(2)

In the equation $\sigma_{\gamma\gamma}$ is the cross section for pair
production $\gamma\gamma\rightarrow e^{+}e^{-}$, Here $dn/d\varepsilon
d\Omega(\vec{x})$ is the number density of target photons per unit energy per
unit solid angle at the position $\vec{x}$ and direction $\hat{k}$. For
photons above 10 TeV, the main target photons are the cosmic microwave
background (CMB) and the infrared emission by dust. The former is uniform
within the Galaxy while the latter depends on $r$ and $z$. We describe the
calculation of $dn/d\varepsilon d\Omega(\vec{x})$ and discuss the attenuation
by different radiation fields in Appendix A.

The ${\gamma\gamma}$ optical depth $\tau_{\gamma\gamma}$ for a photon
traveling from an initial position $\vec{x}_{0}$ to an observer at
$\vec{x}_{\rm ob}$ is

$\tau_{\gamma\gamma}(E_{\gamma},\vec{x}_{0},\vec{x}_{\rm
ob})=\int_{0}^{|\vec{x}_{\rm
ob}-\vec{x}_{0}|}ds~{}\lambda_{\gamma\gamma}^{-1}(E_{\gamma},\hat{u},\vec{x}_{0}+s\hat{u})$
(3)

with $\hat{u}=(\vec{x}_{\rm ob}-\vec{x}_{0})/|\vec{x}_{\rm ob}-\vec{x}_{0}|$.
The probability for a photon to survive from the pair production is (Vernetto
& Lipari, 2016),

$P_{\gamma,\rm surv}(E_{\gamma},\vec{x}_{0},\vec{x}_{\rm
ob})=\exp\left(-\tau_{\gamma\gamma}(E_{\gamma},\vec{x}_{0},\vec{x}_{\rm
ob})\right)$ (4)

The averaged gamma-ray intensity from a region of solid angle $\Delta\Omega$
is the sum of the photons that have survived from all the sources in that
area,

$\displaystyle E_{\gamma}^{2}\Phi_{\gamma}^{\Omega}$ $\displaystyle\approx$
$\displaystyle\frac{1}{\Delta\Omega}\int d\Omega~{}\int ds\,n_{\rm
CR}n_{N}(\vec{x}_{0})\sigma_{pp}c$ $\displaystyle\times$
$\displaystyle\left(E_{\gamma}^{2}\frac{dN_{\gamma}}{dE_{\gamma}}\right)\frac{1}{4\pi|\vec{x}_{0}-\vec{x}_{\rm
ob}|^{2}}P_{\gamma,\rm surv}(E_{\gamma},\vec{x}_{0},\vec{x}_{\rm ob})$

where $n_{\rm CR}n_{N}(\vec{x}_{0})$ is the product of the CR and
gas/molecular densities. $\sigma_{pp}$ is the inelastic $pp$ cross section,
which moderately increases from $\sim 40$ mb at $E_{p}=100$ TeV to $\sim 70$
mb at $E_{p}=10$ PeV (e.g., Tanabashi et al., 2018). This energy dependence
impacts the shapes of the intrinsic gamma-ray spectrum and the neutrino
spectrum in the same way.

The all-flavor neutrino flux can be computed in the same way,

$\displaystyle E_{\nu}^{2}\Phi_{\nu}^{\Omega}\approx\frac{1}{\Delta\Omega}\int
d\Omega\int ds\frac{n_{\rm
CR}n_{N}(\vec{x}_{0})\sigma_{pp}c}{4\pi|\vec{x}_{0}-\vec{x}_{\rm
ob}|^{2}}\left(E_{\nu}^{2}\frac{dN_{\nu}}{dE_{\nu}}\right).$ (6)

Notice that the chance of a neutrino interaction inside the Galaxy is
negligible. In other words, the survival probability of high-energy neutrinos
is always 1,

Using equation (1), we can write the neutrino intensity using the gamma-ray
intensity in the Galactic coordinate,

$\displaystyle
E_{\nu}^{2}\Phi_{\nu}^{\Omega}\approx\frac{3}{2}\left.\left(E_{\gamma}^{2}\Phi_{\gamma}^{\Omega}\right)\right|_{E_{\gamma}=2E_{\nu}}$
$\displaystyle\times$ $\displaystyle\frac{\int ds\int\cos{b}\,db\int
dl\,n_{\rm CR}n_{N}(s,b,l)}{\int ds\int\cos{b}\,db\int dl\,n_{\rm
CR}n_{N}(s,b,l)P_{\gamma,\rm surv}(E_{\gamma}=2E_{\nu},s,b,l)},$

which is improved compared to equation 2 of Ahlers & Murase (2014).

To account for the uncertainty caused by the $n_{\rm CR}n_{N}$ distribution,
which depends on details of the CR propagation and gas density distribution,
we consider the two limits. In the first model, as in the leaky box model,
cosmic rays are assumed to be uniformly distributed within the disk (that is
smaller than the cosmic-ray halo),

$n_{\rm CR}n_{N}(r,z)\propto\begin{cases}1&r<R_{\rm
disk}\,\mbox{and}\,|z|<z_{\rm disk}\\\ 0&\mbox{otherwise }\end{cases},$ (8)

with $R_{\rm disk}=15$ kpc and $z_{\rm disk}=0.2$ kpc (that is compatible with
the scale height of HI gas). This gives a conservative estimate on the
Galactic diffuse emission (see, e.g., Ahlers & Murase, 2014; Ahlers et al.,
2016). In the second model, we assume that it follows the spatial distribution
of supernova remnants (SNRs),

$n_{\rm
CR}n_{N}(r,z)\propto\left(\frac{r}{R_{\odot}}\right)^{\zeta}\exp{\left[-\eta\left(\frac{r-R_{\odot}}{R_{\odot}}\right)-\frac{|z|}{z_{g}}\right]}.$
(9)

where $R_{\odot}=8.5$ kpc is the solar distance from the GC and the following
parameter values are adopted, $\zeta=1.09$, $\eta=3.87$ (Green, 2015) and
$z_{g}=0.2$ kpc (Blasi & Amato, 2012). This model is appropriate when we
consider GP emission from a number of discrete sources. But this would give a
optimistic estimate on the Galactic diffuse emission since cosmic rays diffuse
out from the sources.

Finally, the total neutrino (or gamma-ray) flux from the region of interest
with $\Omega$ can be derived by the flux from the observed sky region,
$\Omega_{\rm obs}$, through

$\displaystyle E_{i}^{2}F_{i,\rm GP}$ $\displaystyle=$ $\displaystyle
E_{i}^{2}F_{i,\rm obs}$ (10) $\displaystyle\times$ $\displaystyle\frac{\int
ds\int^{\Omega}d\Omega\,n_{\rm CR}n_{N}(s,b,l)P_{i,\rm surv}(s,b,l)}{\int
ds\int^{\Omega_{\rm obs}}d\Omega\,n_{\rm CR}n_{N}(s,b,l)P_{i,\rm
surv}(s,b,l)}\,\,\,\,\,$ $\displaystyle=$ $\displaystyle\frac{\Delta\Omega
n_{\rm CR}^{\Omega}\tau_{pp}^{\Omega}\tilde{P}_{i,\rm
surv}^{\Omega}}{\Delta\Omega_{\rm obs}n_{\rm CR}^{\Omega_{\rm
obs}}\tau_{pp}^{\Omega_{\rm obs}}\tilde{P}_{i,\rm surv}^{\Omega_{\rm obs}}}$
(11)

where $i$ indicates either $\nu$ or $\gamma$, $n_{\rm CR}^{\Omega}$ is the
averaged CR density and $\tau_{pp}^{\Omega}$ is the averaged $pp$ optical
depth (Ahlers & Murase, 2014).

Figure 1: The diffuse Galactic gamma-ray intensity from two sky regions,
region A: $25^{\circ}<l<100^{\circ}$, $|b|<5^{\circ}$, and region B:
$50^{\circ}<l<200^{\circ}$, $|b|<5^{\circ}$. The red and orange data points
are the Tibet AS$\gamma$ measurement of the diffuse $\gamma$-ray emission from
the two regions (Amenomori et al., 2021). In the last energy bin, the fainter
data points indicate the residual intensity after removing events relevant to
Cygnus Cocoon. The red and orange bands are the best-fit $\gamma$-ray models
derived in this work, accounting for uncertainties in the gamma-ray
attenuation and cosmic-ray models. The brown long and short dashed curves
indicate the diffuse gamma-ray spectra for the GP and unresolved hypernova
remnants, respectively, which are taken from Ahlers & Murase (2014) for
$|b|<5^{\circ}$. The thin solid curve shows the sum of the two components,
which demonstrates that the Tibet AS$\gamma$ data are also consistent with a
prediction with source emission.  Figure 2:  All-sky-averaged intensity of all
flavor diffuse neutrinos from the GP, compared to neutrino observations. The
GP neutrino intensity, $E_{\nu}^{2}\Phi_{\nu}$, (blue shaded band) is derived
with the best-fit gamma-ray intensities in Figure 1. The model is consistent
with the combined upper limits at 90% confidence level posed by ANTARES and
IceCube (sky blue dashed and dash-dotted curves; Albert et al., 2018), the 90%
limits with 7-year IceCube data (blue dotted curve; Aartsen et al., 2017), and
the upper limits on neutrinos from the GP (blue downward arrows), which are
derived from the CASA-MIA gamma-ray limits in region B, assuming that sources
follow the SNR distribution (cyan downward arrows; Borione et al., 1998). The
hatched band shows the intensity $E_{\nu}^{2}\Phi^{\Omega}_{\nu}$ of the
$|b|<5^{\circ}$ region, which is comparable to the isotropic neutrino
background from the IceCube Cascade (light blue data points; Aartsen et al.,
2020) and muon neutrino (light blue shaded area; Stettner, 2019) data below
$\sim 100$ TeV.

###  2.2 Results

Two neutrino spectral models are used to account for the uncertainty in the
cosmic-ray nucleon spectrum (see Appendix B). The gamma-ray intensity
$\Phi_{\gamma}^{\Omega}$ from a sky region of solid angle $\Omega$ is obtained
from equations (2.1, B1, B2) and normalized by fitting to the Tibet AS$\gamma$
data using a $\chi^{2}$ statistic,

$\displaystyle\chi^{2}=\sum_{i,j}\frac{(\Phi^{\Omega_{j}}_{\gamma,\rm
obs}(E_{i})-\Phi_{\gamma}^{\Omega_{j}}(E_{i}))^{2}}{\left(\sigma^{\Omega_{j}}_{\gamma,\rm
obs}(E_{i})\right)^{2}}.$ (12)

In this equation, $\Phi^{\Omega_{j}}_{\gamma,\rm obs}(E_{i})$ and
$\sigma^{\Omega_{j}}_{\gamma,\rm obs}(E_{i})$ are the observed intensity and
uncertainty of diffuse gamma rays from sky region $j$ in energy bin $i$,
respectively. Two sky regions are considered, namely, region A with
$25^{\circ}<l<100^{\circ}$ and $|b|<5^{\circ}$, and region B with
$50^{\circ}<l<200^{\circ}$ and $|b|<5^{\circ}$. In each energy bin, the upper
(lower) error is used if the model is above (below) the mean of the
measurement. Amenomori et al. (2021) notes that $40\%$ events above 398 TeV in
each of the two sky regions are close to the Cygnus Cocoon (see our discussion
in Section 3.1.1). Therefore, our fit in the highest energy bin uses $60\%$ of
the measured values. The fit uses six flux points from Tibet AS$\gamma$ and
one free parameter, the flux norm of equation B1 or 8, and thus has a total of
five degrees of freedom. The resulted chi-square per degree of freedom is
$\chi^{2}\sim 1.6-2.0$ and $\chi^{2}\sim 0.7-1.0$ for the uniform and SNR
source distribution models, respectively.

The best-fit gamma-ray intensities for the two sky regions measured by Tibet
AS$\gamma$ are shown as colored bands in Figure 1. The boundaries of the bands
are decided by the minimum and maximum values from the four cases in use,
which include the two gamma-ray emissivity distribution models and the two
cosmic-ray nucleon spectra. In particular, in these shaded bands, upper
boundaries below $\sim 150$ TeV are determined by the nucleon spectrum Model A
(equation B1), whereas those above are governed by the nucleon spectrum Model
B (equation B2). See also Appendix B. We find that the sub-PeV gamma-ray
spectral shape is barely affected by uncertainties in the gamma-ray
attenuation. Rather, it depends on the shape of the cosmic-ray nucleon
spectrum. For example, if the cosmic-ray nucleon spectrum has a break at
$E_{*}=0.9$ PeV (Gaisser et al., 2013), the gamma-ray spectrum should be
steepened at $\sim 0.08E_{*}\sim 80$ TeV, which causes a tension with the
Tibet AS$\gamma$ data at the highest-energy bin. The tension can readily
alleviated if discrete sources make a significant contribution at sub-PeV
energies as we discuss in Section 3.

Now, let us consider the GP contribution to IceCube neutrinos. Figure 2 shows
the $4\pi$-averaged intensity of diffuse Galactic neutrinos (shaded blue band)
calculated using the best-fit gamma-ray models in Figure 1, equations (10) and
(2.1). The GP contribution to the all-sky neutrino flux depends on neutrino
energy, and is found to be $\lesssim 5-10$% in the 100 TeV range. We caution
that the exact value depends on the all-sky IceCube flux and varies with
energy depending on the spectra. Analytically, noting that the all-sky-
averaged intensity is denoted as $\Phi_{\nu}\equiv(\int
d\Omega~{}\Phi_{\nu}^{\Omega}/4\pi)=(\Delta\Omega_{\rm
GP}/4\pi)\Phi_{\nu}^{\rm GP}$, we may write

$\displaystyle\frac{\Delta\Omega_{\rm GP}E_{\nu}^{2}\Phi_{\nu}^{\rm GP}}{4\pi
E_{\nu}^{2}\Phi_{\nu}^{\rm IC}}\sim 5\%{\left(\frac{E_{\nu}^{2}\Phi_{\nu}^{\rm
IC}}{5\times{10}^{-8}~{}{\rm GeV}{\rm cm}^{-2}~{}{\rm s}^{-1}~{}{\rm
sr}^{-1}}\right)}^{-1}$ $\displaystyle\times\left(\frac{\Delta\Omega_{\rm
GP}}{1.1~{}\rm sr}\right)\left(\frac{E_{\gamma}^{2}\Phi_{\gamma}^{\rm
GP}|_{E_{\gamma}=2E_{\nu}}}{2\times{10}^{-8}~{}{\rm GeV}{\rm cm}^{-2}~{}{\rm
s}^{-1}~{}{\rm sr}^{-1}}\right).\,\,\,\,\,\,\,\,\,\,$ (13)

Here the effective solid angle of the GP is assumed to be $\Delta\Omega_{\rm
GP}=2\sin b\Delta l\simeq 1.1~{}{\rm sr}~{}{(|b|/5^{\circ}}){(\Delta
l/360^{\circ})}$ and $E_{\gamma}^{2}\Phi_{\gamma}^{\rm GP}\sim
E_{\gamma}^{2}\Phi_{\gamma}^{\Omega_{\rm GP}}$ 111The all-sky estimate may
depend on the distribution in the vertical direction of the disk. For a a
scale height of $0.1$ kpc, we have $E_{\gamma}^{2}\Phi_{\gamma}^{\rm
GP}/E_{\gamma}^{2}\Phi_{\gamma}^{\Omega_{\rm GP}}\simeq 1.6$ (Ahlers & Murase,
2014).. This value is also consistent with previous results (Ahlers & Murase,
2014; Ahlers et al., 2016). Following Ahlers & Murase (2014), we also show
limits on GP neutrino emission by converting the CASA-MIA upper limits through
equations 2.1 and 10, assuming that sources follow SNR distribution. As
clearly seen in Figure 2, our results obtained with the new Tibet AS$\gamma$
data are also consistent with the upper limits posed by Albert et al. (2018)
and Aartsen et al. (2017) through independent neutrino observations. This is
not surprising because IceCube neutrinos are consistent with an isotropic
distribution, and they mostly come from the region outside the GP
(Kheirandish, 2020, for a review).

The hatched gray band in Figure 2 shows the neutrino intensity of the GP with
$|b|<5^{\circ}$. The Tibet AS$\gamma$ data imply that the intensity of diffuse
Galactic neutrinos can be comparable to the isotropic diffuse neutrino
background in the GP region in the $\sim 10-100$ TeV range, i.e.,

$E_{\nu}^{2}\Phi_{\nu}^{{\Omega}_{\rm GP}}\sim E_{\nu}^{2}\Phi_{\nu}^{\rm
IC}.$ (14)

Our result suggests that the diffuse Galactic emission below $\sim 100$ TeV
can be promisingly observed by next-generation neutrino telescopes such as
IceCube-Gen2 (The IceCube-Gen2 Collaboration et al., 2020), KM3Net (Adrián-
Martínez et al., 2016), and Baikal-GVD (Baikal-GVD Collaboration et al.,
2018). At higher energies, the source contribution can be relevant (see the
HNR curve in Figure 2 obtained by averaging the flux within $b=10^{\circ}$),
which is also encouraging for next-generation high-energy neutrino telescopes.
Even non-detections would also be useful, because as discussed in the next
section, the diffuse gamma rays may potentially come from discrete sources.

##  3 Contribution of Discrete sources to Sub-PeV Gamma Rays and Neutrinos

The sub-PeV diffuse emission discussed in Section 2 comes from ions that were
injected into the Galaxy by sources $\sim 10^{4}-10^{6}$ yr ago, propagate in
the Galactic magnetic field, and interact with gas and dust in the Milky Way.
In this section we focus on the contribution of discrete sources to the sub-
PeV gamma rays observed by the Tibet AS$\gamma$ experiment.

Sub-PeV gamma rays from discrete sources are naturally expected for two main
reasons. First, as noted in Section 2, a nucleon spectrum with a break energy
lower than $\sim 1$ PeV, which is suggested by the composition modeling around
the knee, is in tension with the last Tibet AS$\gamma$ data point. Second,
theoretically, it is natural that gamma rays are produced inside or in the
vicinity of Pevatrons, including the recently detected TeV source Cygnus
Cocoon and other unresolved hadronic and leptonic source populations. Harder
spectra of source neutrinos and gamma rays may dominate over the steep spectra
of diffuse secondaries at the highest energies.

###  3.1 Hadronic Sources

####  3.1.1 Cygnus Cocoon

Figure 3:  Spectral energy distribution of the Cygnus Cocoon measured by
Fermi-LAT (Abdollahi et al., 2020), ARGO-YBJ (Bartoli et al., 2014), and HAWC
(Abeysekara et al., 2021). The light pink and orange flux points indicate 40%
of the Tibet AS$\gamma$ flux of regions A and B (Amenomori et al., 2021). The
thick pink and orange markers additionally scale the fluxes to the HAWC size
of the Cygnus Cocoon. The blue dotted curve shows the limit on the
$\gamma$-ray flux based on the non-detection of neutrinos from the region by
IceCube (Kheirandish & Wood, 2019). The two $\gamma$-ray emission models from
Abeysekara et al. (2021) are shown for comparison. A significant detection of
the Cygnus Cocoon at the estimated flux level may favor the burst model and
the presence of a Pevatron.

Amenomori et al. (2021) indicates that $\sim 40\%$ of the events in their
highest-energy bin are detected within $4^{\circ}$ around the Cygnus Cocoon.
Extended $1-200$ TeV gamma-ray emission from the Cygnus Cocoon has recently
been reported by the HAWC Observatory, with emission above 100 TeV detected at
$\sim 2.4\sigma$ significance level (Abeysekara et al., 2021). The gamma-ray
spectrum can be explained by protons that either have been continuously
injected over the lifetime of the stellar clusters (a few Myr), or were
produced by a recent (sub-Myr) burst-like event. The latter scenario invokes
the presence of a Pevatron, specifically, a PeV proton accelerator with a hard
spectrum $dN_{p}/dE_{p}\propto E_{p}^{-2.1}$.

The Cygnus OB2 association (e.g., Wright, 2020) has been among the most
promising sites for cosmic-ray acceleration for many years (see Bykov, 2014,
for a review). Superbubbles may accelerate cosmic rays up to PeV energies or
beyond via multiple shocks and turbulence (Bykov & Fleishman, 1992; Klepach et
al., 2000). It has been shown that cosmic rays escaping from star clusters or
superbubbles may explain cosmic rays above the knee energy (Murase & Fukugita,
2019; Zhang et al., 2020). The Cygnus region has also been of much interest as
the promising source of high-energy neutrinos (e.g., Anchordoqui et al., 2007;
Beacom & Kistler, 2007; Halzen & O Murchadha, 2007; Halzen et al., 2017).

Figure 3 suggests that observations of the Cygnus Cocoon region above 200 TeV
may disentangle the continuous and burst scenarios of Abeysekara et al.
(2021). The light pink and orange data points correspond to 40% of the Tibet
AS$\gamma$ flux above 398 TeV from regions A and B. Since the radius of the
Cocoon is measured to be $\sim 2^{\circ}$ at $1-100$ TeV by HAWC, the thick
pink and orange markers show a more conservative estimation of the Tibet
AS$\gamma$ flux using the average event number within the HAWC radius. We
caveat that the pink and orange flux points in Figure 3 are approximate. The
actual flux depends on the $\gamma$-ray morphology and the detector exposure.
No high-energy neutrino emission has been detected from the Cygnus Cocoon. The
blue dotted curve shows the gamma-ray upper limit converted from the IceCube
limit on this source (Kheirandish & Wood, 2019). For comparison, we overlay
the continuous model and the burst model from Abeysekara et al. (2021). In
particular, we here update the maximum proton energy in the burst model from 2
PeV to 10 PeV to accommodate the estimated Tibet AS$\gamma$ flux. The other
model parameters remain the same. We find that the burst scenario, hence a
Pevatron, would be favored if the flux above 400 TeV reaches $\gtrsim 3\times
10^{-13}\,\rm TeV\,cm^{-2}\,s^{-1}$ assuming $\sim$30% measurement
uncertainty.

####  3.1.2 Hypernova Remnants

Recent optical observations have revealed that energetic supernovae with a
kinetic energy of ${\mathcal{E}}_{\rm ej}\gtrsim 10^{52}$ erg are not
negligible as the cosmic-ray energy budget (e.g., Murase & Fukugita, 2019),
and their rate is about $\sim 3$% of the core-collapse supernova rate that is
$\sim 3$ per century. Energetic supernovae so-called hypernovae (that are
mostly broad-line Type Ibc supernovae) and trans-relativistic supernovae,
which are often associated with low-luminosity gamma-ray bursts, have been
widely discussed as cosmic-ray accelerators, and they may accelerate cosmic
rays up to $\sim 10-100$ PeV energies (e.g., Sveshnikova, 2003; Murase et al.,
2013; Senno et al., 2015). It has been argued that the X-ray emission from the
Cygnus region can be attributed to a hypernova remnant (Kimura et al., 2013;
Bluem et al., 2020). The burst model discussed in the previous subsection is
consistent with such a model. The required cosmic-ray input, $\sim 10^{51}$
erg, is comparable to the energy amount of cosmic rays accelerated by a
hypernova. Dozens of hypernova remnants are expected to exist in the Milky
Way, and a fraction of IceCube neutrinos may come from them (Fox et al.,
2013). There may be a few hypernova remnants in region A and B, and one of
them could be in the Cygnus region. As shown in Figure 1, the model (Ahlers &
Murase, 2014) is consistent with the Tibet AS$\gamma$ data. This demonstrates
the potential relevance of contributions from discrete sources, and we stress
that other candidate sources are also possible.

###  3.2 Leptonic Sources

Figure 4:  Demonstration of a hybrid $\gamma$-ray emission model, in which the
inverse Compton of relativistic electrons (grey dashed curve) explains the
Tibet AS$\gamma$ measurement in the region $25^{\circ}<l<100^{\circ}$ (red
round data points), and $\pi^{0}$ decay by Galactic diffuse protons (grey
dash-dotted curve) explains the lower-energy observations of the same region
by EGRET (black plus markers; Hunter et al., 1997), Fermi-LAT (brown shaded
region, scaled from Ackermann et al., 2012b to the EGRET flux), and ARGO-YBJ
(pink triangle data points; Bartoli et al., 2015). The electrons are assumed
to have an intrinsic spectrum $dN/dE_{e}\propto E_{e}^{-2}$ and maximum energy
$E_{e,\rm max}=3$ PeV.

While the sub-PeV $\gamma$-ray emission can be plausibly explained by the
decay of neutral pions from hadronuclear interactions (Amenomori et al.,
2021), a leptonic origin may not yet be excluded. Figure 4 demonstrates such a
scenario, where the Tibet AS$\gamma$ data can in principle explained by
electrons that upscatter the CMB. We here assume that relativistic electrons
are injected by discrete Galactic sources such as pulsar wind nebulae,
confined close to the vicinity of the emission region, while being cooled via
synchrotron radiation in the Galactic magnetic field and inverse-Compton
scattering with the CMB. The steady-state electron distribution is calculated
by solving the transport equation, with $B=3~{}\mu$G and $u_{\rm
CMB}=0.26\,\rm eV\,cm^{-3}$ for the energy density of magnetic field and the
CMB. We note that $B$ near the sources could be higher than the average ISM
field strength that we take. Besides, our example model does not account for
the IR radiation field at the emission sites which could further contribute to
gamma-ray production below $\sim 10$ TeV. At tens to hundreds of TeV electron
energies the cooling timescale is much shorter than the diffusion timescale,
therefore the diffusion process is negligible for the calculation. Assuming
that $Q_{e}\propto E_{e}^{-2}$ up to a maximum energy $E_{e,\rm max}=3~{}\,\rm
PeV$ and source emission time $\sim 0.1$ Myr (just for demonstrative
purposes), we find that the total electron power
$L_{e}=\int_{m_{e}c^{2}}^{E_{e,\rm max}}\sim 10^{37}\,(d/5\,{\rm
kpc})^{2}\,\rm erg\,s^{-1}$ is sufficient to explain the Tibet AS$\gamma$
flux, where $d$ is the average source distance. Although the mechanism of PeV
electron acceleration remains an open question, leptonic sources with such
hard spectra and high $E_{e,\rm max}$ have been previously observed (see,
e.g., HAWC Collaboration et al., 2018, 2019; Fang et al., 2020).

The Galactic diffuse emission at GeV-TeV energies is expected to be
contributed by the $\pi^{0}$ decay of hadronic cosmic rays and the inverse-
Compton emission by diffuse electrons and electrons from pulsars (Ackermann et
al., 2012b; Linden & Buckman, 2018). The dash-dotted curve in Figure 4 shows
such a component. It shows one of our benchmark case (cosmic-ray model I and
uniform source distribution) but with 50% of the best-fit normalization.

##  4 Conclusions and Discussion

The flux level of Galactic neutrinos has been a mystery and remains undetected
by the current-generation neutrino telescopes. Using the recently measured
Galactic diffuse gamma-ray intensity, we derived the neutrino flux from the
GP. We took into account uncertainties in the gamma-ray attenuation and
cosmic-ray spectrum, and the all-sky-averaged neutrino intensity from the GP
is estimated to be $\lesssim(3-6)\times 10^{-9}\,\rm
GeV\,cm^{-2}\,s^{-1}\,sr^{-1}$ at 100 TeV. Our results are consistent with
both previous gamma-ray and neutrino constraints. Our calculation relies on
the simple connection between neutrino and gamma rays produced by $pp$
interactions, so the results on the sub-PeV neutrino flux are not much
sensitive to the chemical composition of cosmic rays as well as details of the
propagation and source distribution.

The diffuse neutrino intensity toward the GP may be comparable to the
extragalactic neutrino intensity in the same sky region. Given that the GP is
$\sim 10$% of the sky, the detection is promising for next-generation
telescopes such as KM3Net and IceCube-Gen2, which may also find a large-scale
anisotropy due to the Galactic component (Becker Tjus & Merten, 2020). The
seven-year all-flavor IceCube GP sensitivity around 100 TeV is $\sim
10^{-8}\,\rm GeV\,cm^{-2}\,s^{-1}\,sr^{-1}$. With a factor of $\sim 5$
increase in the neutrino effective area (The IceCube-Gen2 Collaboration et
al., 2020), the GP sensitivity of IceCube-Gen2 may reach $\sim 3\times
10^{-9}\,\rm GeV\,cm^{-2}\,s^{-1}\,sr^{-1}$ in ten years, while KM3Net (Aiello
et al., 2019) would reach $\sim(3-6)\times 10^{-9}\,\rm
GeV\,cm^{-2}\,s^{-1}\,sr^{-1}$ at 100 TeV for the diffuse Galactic neutrino
spectrum, although the detailed value depends on the spectral shape. These
would be sufficient to detect the diffuse neutrino flux indicated by the lower
bound of the blue shaded region in Figure 2.

The origin of neutrinos below 100 TeV has emerged as a new mystery (Aartsen et
al., 2015, 2020). The measured spectral index ($\sim 2.53$) is compatible with
that of diffuse GP emission. With these similar indices, the Tibet AS$\gamma$
data imply that the GP contribution remains $\lesssim 10$% even at $\sim
10-100$ energies, further supporting the manifestation of extragalactic
origins (Murase et al., 2016). We do not exclude possibilities that a fraction
of neutrinos come from other regions such as the Galactic halo (see Section 2
of Ahlers & Murase, 2014), although the Tibet AS$\gamma$ off-source data imply
a Galactic halo contribution lower than $E_{\nu}^{2}\Phi_{\nu}\lesssim
2\times{10}^{-9}~{}{\rm GeV}~{}{\rm cm}^{-2}~{}{\rm s}^{-1}~{}{\rm sr}^{-1}$.
Other Galactic sources, such as the Fermi bubbles, have also been constrained
(Ahlers & Murase, 2014; Lunardini et al., 2014; Fang et al., 2017). These are
consistent with the fact there is no significant northern-southern asymmetry
in the neutrino sky (Aartsen et al., 2020).

While the diffuse Galactic interpretation of the Tibet AS$\gamma$ data seems
the most natural, discrete sources may still significantly contribute
especially at the highest energies. This is especially the case if the cosmic-
ray nucleon spectrum is as steep as $E^{-2.7}$ with a break energy of
$\lesssim 1$ PeV. If a crucial fraction of the highest-energy events detected
by the Tibet AS$\gamma$ experiment is associated with the Cygnus Cocoon, the
presence of an efficient Pevaron would be supported. The Tibet AS$\gamma$ data
can also be explained by unresolved Pevatrons such as hypernova remnants in
the Cygnus region and (or) other part of the Galaxy. Finally, the leptonic
scenario is not excluded. Future multi-messenger observations by not only
neutrino telescopes but also near-future gamma-ray experiments such as LHAASO,
ALPACA, and SWGO are necessary to discriminate between these scenarios. The
spatial distribution would give us crucial information, and observations in
the southern sky are relevant (Ahlers & Murase, 2014; Huentemeyer et al.,
2019). A few or dozens of sources are sufficient to explain the sub-PeV gamma-
ray intensity, which is promising for source identification.

While preparing the manuscript, we became aware that Dzhatdoev (2021); Qiao et
al. (2021); Liu & Wang (2021) appeared on arXiv. Our work was carried out
independently.

We thank Markus Ahlers, Julia Becker Tjus, Kazumasa Kawata, and Walter Winter
for useful comments and communications. The work of K.F. is supported by the
Office of the Vice Chancellor for Research and Graduate Education at the
University of Wisconsin–Madison with funding from the Wisconsin Alumni
Research Foundation. The work of K.M. is supported by NSF Grant No.
AST-1908689, and KAKENHI No. 20H01901 and No. 20H05852.

## References

  * Aartsen et al. (2013a) Aartsen, M., et al. 2013a, Phys.Rev.Lett., 111, 021103
  * Aartsen et al. (2013b) —. 2013b, Science, 342, 1242856
  * Aartsen et al. (2015) Aartsen, M. G., et al. 2015, Astrophys. J., 809, 98
  * Aartsen et al. (2017) Aartsen, M. G., Ackermann, M., Adams, J., et al. 2017, The Astrophysical Journal, 849, 67.  https://doi.org/10.3847/1538-4357/aa8dfb
  * Aartsen et al. (2019) Aartsen, M. G., Ackermann, M., Adams, J., et al. 2019, Phys. Rev. D, 100, 082002
  * Aartsen et al. (2020) Aartsen, M. G., et al. 2020, Phys. Rev. Lett., 125, 121104
  * Abdollahi et al. (2020) Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33
  * Abeysekara et al. (2021) Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2021, Nature Astronomy, arXiv:2103.06820
  * Acero et al. (2016) Acero, F., et al. 2016, Astrophys. J. Suppl., 223, 26
  * Ackermann et al. (2012a) Ackermann, M., et al. 2012a, Astrophys.J., 750, 3
  * Ackermann et al. (2012b) Ackermann, M., Ajello, M., Atwood, W. B., et al. 2012b, The Astrophysical Journal, 750, 3.  https://doi.org/10.1088/0004-637x/750/1/3
  * Adrián-Martínez et al. (2016) Adrián-Martínez, S., Ageron, M., Aharonian, F., et al. 2016, Journal of Physics G Nuclear Physics, 43, 084001
  * Ahlers et al. (2016) Ahlers, M., Bai, Y., Barger, V., & Lu, R. 2016, Phys. Rev. D, 93, 013009
  * Ahlers & Murase (2014) Ahlers, M., & Murase, K. 2014, Phys.Rev., D90, 023010
  * Aiello et al. (2019) Aiello, S., Akrame, S., Ameli, F., et al. 2019, Astroparticle Physics, 111, 100\.  https://www.sciencedirect.com/science/article/pii/S0927650518302809
  * Albert et al. (2018) Albert, A., et al. 2018, Astrophys. J. Lett., 868, L20
  * Albert et al. (2018) Albert, A., André, M., Anghinolfi, M., et al. 2018, ApJ, 868, L20
  * Amenomori et al. (2021) Amenomori, M., Bao, Y. W., Bi, X. J., et al. 2021, Phys. Rev. Lett., 126, 141101\.  https://link.aps.org/doi/10.1103/PhysRevLett.126.141101
  * Anchordoqui et al. (2007) Anchordoqui, L. A., Beacom, J. F., Goldberg, H., Palomares-Ruiz, S., & Weiler, T. J. 2007, Phys. Rev. D, 75, 063001
  * Anchordoqui et al. (2014) Anchordoqui, L. A., Goldberg, H., Lynch, M. H., et al. 2014, Phys.Rev., D89, 083003
  * Baikal-GVD Collaboration et al. (2018) Baikal-GVD Collaboration, :, Avrorin, A. D., et al. 2018, arXiv e-prints, arXiv:1808.10353
  * Bartoli et al. (2014) Bartoli, B., Bernardini, P., Bi, X. J., et al. 2014, ApJ, 790, 152
  * Bartoli et al. (2015) Bartoli, B., Bernardini, P., Bi, X. J., et al. 2015, The Astrophysical Journal, 806, 20.  https://doi.org/10.1088/0004-637x/806/1/20
  * Beacom & Kistler (2007) Beacom, J. F., & Kistler, M. D. 2007, Phys. Rev. D, 75, 083001
  * Becker Tjus & Merten (2020) Becker Tjus, J., & Merten, L. 2020, Phys. Rep., 872, 1
  * Blasi & Amato (2012) Blasi, P., & Amato, E. 2012, J. Cosmology Astropart. Phys, 2012, 011
  * Bluem et al. (2020) Bluem, J., Kaaret, P., Fuelberth, W., et al. 2020, Astrophys. J., 905, 91
  * Borione et al. (1998) Borione, A., Catanese, M. A., Chantell, M. C., et al. 1998, The Astrophysical Journal, 493, 175.  https://doi.org/10.1086/305096
  * Bykov (2014) Bykov, A. M. 2014, Astron. Astrophys. Rev., 22, 77
  * Bykov & Fleishman (1992) Bykov, A. M., & Fleishman, G. D. 1992, MNRAS, 255, 269
  * Denton et al. (2017) Denton, P. B., Marfatia, D., & Weiler, T. J. 2017, JCAP, 08, 033
  * Draine (2003) Draine, B. T. 2003, ARA&A, 41, 241
  * Dzhatdoev (2021) Dzhatdoev, T. 2021, arXiv e-prints, arXiv:2104.02838
  * Fang et al. (2020) Fang, K., Charles, E., & Blandford, R. D. 2020, ApJ, 889, L5
  * Fang et al. (2017) Fang, K., Su, M., Linden, T., & Murase, K. 2017, Phys. Rev. D, 96, 123007
  * Fox et al. (2013) Fox, D., Kashiyama, K., & Mészáros, P. 2013, Astrophys.J., 774, 74
  * Gaggero et al. (2015) Gaggero, D., Grasso, D., Marinelli, A., Urbano, A., & Valli, M. 2015, ApJ, 815, L25
  * Gaisser (2012) Gaisser, T. K. 2012, Astropart. Phys., 35, 801
  * Gaisser (2013) Gaisser, T. K. 2013, in European Physical Journal Web of Conferences, Vol. 52, European Physical Journal Web of Conferences, 09004
  * Gaisser et al. (2013) Gaisser, T. K., Stanev, T., & Tilav, S. 2013, Front. Phys. China, 8, 748
  * Green (2015) Green, D. A. 2015, MNRAS, 454, 1517
  * Halzen et al. (2017) Halzen, F., Kheirandish, A., & Niro, V. 2017, Astropart. Phys., 86, 46
  * Halzen & O Murchadha (2007) Halzen, F., & O Murchadha, A. 2007, Phys. Rev. D, 76, 123003
  * HAWC Collaboration et al. (2018) HAWC Collaboration, Abeysekara, A. U., Albert, A., et al. 2018, arXiv e-prints, arXiv:1810.01892
  * HAWC Collaboration et al. (2019) —. 2019, arXiv e-prints, arXiv:1909.08609
  * Hayakawa (1952) Hayakawa, S. 1952, Progress of Theoretical Physics, 8, 571
  * Hörandel (2003) Hörandel, J. R. 2003, Astroparticle Physics, 19, 193
  * Huentemeyer et al. (2019) Huentemeyer, P., BenZvi, S., Dingus, B., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 109
  * Hunter et al. (1997) Hunter, S. D., et al. 1997, Astrophys. J., 481, 205
  * Hunter et al. (1997) Hunter, S. D., Bertsch, D. L., Catelli, J. R., et al. 1997, ApJ, 481, 205
  * Joshi et al. (2013) Joshi, J. C., Winter, W., & Gupta, N. 2013, Mon.Not.Roy.Astron.Soc., 3414, 3419
  * Kachelrieß & Ostapchenko (2014) Kachelrieß, M., & Ostapchenko, S. 2014, Phys.Rev., D90, 083002
  * Kamae et al. (2006) Kamae, T., Karlsson, N., Mizuno, T., Abe, T., & Koi, T. 2006, ApJ, 647, 692
  * Kelner et al. (2006) Kelner, S. R., Aharonian, F. A., & Bugayov, V. V. 2006, Phys. Rev., D74, 034018, [Erratum: Phys. Rev. D79, 039901(E) (2009)]
  * Kheirandish (2020) Kheirandish, A. 2020, Astrophys. Space Sci., 365, 108
  * Kheirandish & Wood (2019) Kheirandish, A., & Wood, J. 2019, in International Cosmic Ray Conference, Vol. 36, 36th International Cosmic Ray Conference (ICRC2019), 932
  * Kimura et al. (2013) Kimura, M., Tsunemi, H., Tomida, H., et al. 2013, PASJ, 65, 14
  * Klepach et al. (2000) Klepach, E. G., Ptuskin, V. S., & Zirakashvili, V. N. 2000, Astroparticle Physics, 13, 161
  * Linden & Buckman (2018) Linden, T., & Buckman, B. J. 2018, Phys. Rev. Lett., 120, 121101.  https://link.aps.org/doi/10.1103/PhysRevLett.120.121101
  * Lipari & Vernetto (2018) Lipari, P., & Vernetto, S. 2018, Phys. Rev. D, 98, 043003
  * Lipari & Vernetto (2020) Lipari, P., & Vernetto, S. 2020, Astroparticle Physics, 120, 102441
  * Liu & Wang (2021) Liu, R.-Y., & Wang, X.-Y. 2021, arXiv e-prints, arXiv:2104.05609
  * Lunardini et al. (2014) Lunardini, C., Razzaque, S., Theodoseau, K. T., & Yang, L. 2014, Phys.Rev., D90, 023016
  * Murase et al. (2013) Murase, K., Ahlers, M., & Lacki, B. C. 2013, Phys.Rev., D88, 121301(R)
  * Murase & Beacom (2012) Murase, K., & Beacom, J. F. 2012, JCAP, 1210, 043
  * Murase & Fukugita (2019) Murase, K., & Fukugita, M. 2019, Phys. Rev., D99, 063012
  * Murase et al. (2016) Murase, K., Guetta, D., & Ahlers, M. 2016, Phys. Rev. Lett., 116, 071101
  * Neronov et al. (2014) Neronov, A., Semikoz, D., & Tchernin, C. 2014, Phys.Rev., D89, 103002
  * Palladino & Vissani (2016) Palladino, A., & Vissani, F. 2016, Astrophys. J., 826, 185
  * Qiao et al. (2021) Qiao, B.-Q., Liu, W., Zhao, M.-J., Bi, X.-J., & Guo, Y.-Q. 2021, arXiv e-prints, arXiv:2104.03729
  * Senno et al. (2015) Senno, N., Mészáros, P., Murase, K., Baerwald, P., & Rees, M. J. 2015, Astrophys.J., 806, 24
  * Spurio (2014) Spurio, M. 2014, Phys. Rev. D, 90, 103004
  * Stecker (1977) Stecker, F. W. 1977, Astrophys. J., 212, 60
  * Stecker (1979) —. 1979, Astrophys.J., 228, 919
  * Stettner (2019) Stettner, J. 2019, in International Cosmic Ray Conference, Vol. 36, 36th International Cosmic Ray Conference (ICRC2019), 1017
  * Sveshnikova (2003) Sveshnikova, L. G. 2003, Astron. Astrophys., 409, 799
  * Tanabashi et al. (2018) Tanabashi, M., Hagiwara, K., Hikasa, K., et al. 2018, Phys. Rev. D, 98, 030001\.  https://link.aps.org/doi/10.1103/PhysRevD.98.030001
  * The IceCube-Gen2 Collaboration et al. (2020) The IceCube-Gen2 Collaboration, :, Aartsen, M. G., et al. 2020, arXiv e-prints, arXiv:2008.04323
  * The KASCADE-Grande Collaboration et al. (2013) The KASCADE-Grande Collaboration, :, Apel, W. D., et al. 2013, arXiv e-prints, arXiv:1306.6283
  * Thoudam et al. (2016) Thoudam, S., Rachen, J. P., van Vliet, A., et al. 2016, A&A, 595, A33
  * Vernetto & Lipari (2016) Vernetto, S., & Lipari, P. 2016, Phys. Rev. D, 94, 063009
  * Wright (2020) Wright, N. J. 2020, New A Rev., 90, 101549
  * Zhang et al. (2020) Zhang, Z., Murase, K., & Mészáros, P. 2020, Mon. Not. Roy. Astron. Soc., 492, 2250

##  Appendix A Optical Depth for Sub-PeV Gamma Rays

Figure 5:  Survival probability of gamma rays travelling to an observer in the
solar neighborhood. The attenuation due to the CMB is indicated by dotted
curves and that due to the infrared dust emission is indicated by dashed
curves. Left: the survival probability as a function of gamma-ray energy for
photons starting from the Galactic Center (GC, grey) and a certain direction
($l=90^{\circ},\,b=0^{\circ}$, $s=8.5\,\rm kpc$) inside the sky region where
diffuse gamma rays are observed (light blue). Right: $P_{\rm surv}$ as a
function of distance to the Sun. Gamma-ray energy is indicated in color (red
corresponds to $E_{\gamma}=200$ TeV and blue corresponds to $E_{\gamma}=1$
PeV) and the direction is indicated by the thickness of the curves (thin
dashed indicates the GC direction and thick dashed indicates
$l=90^{\circ},b=0^{\circ}$).

Gamma rays above a few 100 TeV mostly pair-produce with the cosmic microwave
background (CMB). The differential number density of the isotropic black body
emission is

$\left.\frac{dn}{d\varepsilon d\Omega}\right|_{\rm
CMB}=\frac{2\varepsilon^{2}}{h^{3}c^{3}}\frac{1}{e^{\varepsilon/(k_{B}T_{\rm
CMB})}-1},$ (A1)

where $h$ is the Planck constant, $k_{B}$ is the Boltzmann constant, and
$T_{\rm CMB}=2.73$ K is the CMB temperature.

Gamma rays between $\sim$10 TeV and $\sim 1$ PeV also interact with the dust
emission with wavelength $\lambda\gtrsim 50\,\mu$m. We follow Vernetto &
Lipari (2016) for the calculation of the infrared density field of the Galaxy.
The intensity of the infrared emission at the location $\vec{x}$ and direction
$\hat{k}$ is obtained by integrating the dust emissivity along the line-of-
sight $s$,

$I_{\lambda}(\vec{x})=\int_{0}^{\infty}ds\,\eta_{\lambda}(\vec{x}+\hat{k}s).$
(A2)

Here $\eta_{\lambda}$ is the power emitted per unit volume, unit solid angle
and unit wavelength by the dust,

$\eta_{\lambda}=\rho_{d}\kappa_{\lambda}B_{\lambda}(T)$ (A3)

$\kappa_{\lambda}$ is the absorption cross section per mass of dust, for which
we have used the $\kappa_{\lambda}$ values from Draine (2003) (the $R_{V}=3.1$
model). $B_{\lambda}(T)$ is the spectral radiance of a black body,

$B_{\lambda}(T)=\frac{2hc^{2}}{\lambda^{5}}\frac{1}{e^{hc/(\lambda
k_{B}T)}-1}.$ (A4)

The infrared emission comes from cold and warm dust components. Their density
and temperature profiles are assumed to follow the disk structure,

$\rho_{c,w}\propto\exp{\left[-\frac{r}{r^{0}_{c,w}}-\frac{|z|}{z^{0}_{c,w}}\right]}.$
(A5)

We adopt the $\rho^{0}_{c,w}$ and $T^{0}_{c,w}$ expressions described by
equations 16 and 17 of Vernetto & Lipari (2016). The number density of the
dust emission is computed by

$\left.\frac{dn}{d\varepsilon d\Omega}\right|_{\rm
IR}(\vec{x})=\frac{1}{c\,\varepsilon}I_{\varepsilon}(\vec{x})$ (A6)

where $I_{\varepsilon}=I_{\lambda}\lambda/\varepsilon$. The extragalactic
background light (EBL) may also interact with gamma rays. However, it is sub-
dominant comparing to the Galactic dust emission and therefore ignored in our
calculation.

Figure 5 compares the attenuation effect of CMB and infrared photons for gamma
rays traveling from the different directions at various energies. Gamma rays
above $\sim$500 TeV are mostly absorbed by the CMB. Gamma rays between $\sim
50$ TeV and $\sim 500$ TeV may pair produce with the dust emission. The
survival probability depends on the direction and distance of the gamma-ray
source. Due to the spatial distribution of the dust, gamma rays from the inner
Galaxy is more absorbed than those from the other parts of the Galaxy.

##  Appendix B Diffuse Galactic Neutrino Spectrum

Figure 6:  All-sky-averaged intensities of all flavor diffuse neutrinos from
the GP, for the two neutrino models described by equations B1 (Model A) and B2
(Model B) and the two source distribution models described by equations 8
(uniform) and 9 (SNR). The lowest and highest fluxes in each energy bin are
used as the boundary of the blue filled area in Figure 2. The intensity
normalization is set by the Tibet AS$\gamma$ data.

The cosmic-ray nucleon spectrum between 10 TeV and 10 PeV has been studied by
different groups based on the measured cosmic-ray spectrum and modeling of the
chemical composition (e.g., Hörandel, 2003; Gaisser, 2012, 2013; Thoudam et
al., 2016; Lipari & Vernetto, 2020). The fitting to observation relies on the
assumption of the spectral model, including the dependence of the break energy
on the charge and the mass number of each chemical group. It also depends on
the choice of hadronic interaction models. For example, based on direct
measurements of cosmic rays, Hörandel (2003) presented the polygonato model
with the rigidity-dependence hypothesis, where the rigidity cutoff is set to 4
PeV (see also Gaisser, 2012). However, the proton spectral break may exist
well below the knee energy, and nuclei such as helium may be dominant around
the knee (see also Gaisser et al., 2013; Lipari & Vernetto, 2020). More
recently, Lipari & Vernetto (2020) showed that different proton spectral
indices are needed to explain the KASCADE (The KASCADE-Grande Collaboration et
al., 2013) and IceTop/IceCube data (Aartsen et al., 2019) that are interpreted
with QGSJet and Sibyll models. To account for uncertainties in the cosmic-ray
spectrum and composition, we consider two models for the diffuse Galactic
neutrino spectrum. In Model A, following Ahlers & Murase (2014), we adopt

$\frac{dN}{dE_{\nu}}\propto
E_{\nu}^{-\alpha_{\nu}}\left[1+\left(\frac{E_{\nu}}{E^{*}}\right)^{w}\right]^{-\delta/w}$
(B1)

with $\alpha_{\nu}=2.54$, $\delta=0.67$, $E_{\nu,*}=0.04\times 0.9$ PeV and
$w=3.0$. The spectral index is based on the nucleon spectrum obtained by
Gaisser (2013) with a $0.1$ spectral hardening to account for the energy
dependence of the inelastic $pp$ cross section (e.g., Kamae et al., 2006;
Kelner et al., 2006). The nucleon spectral cutoff, 0.9 PeV, is motivated by
the limit that cosmic rays around the knee energy are dominated by helium
nuclei with $A=4$ (Gaisser, 2012). The factor of $\sim 0.04$ converts nucleon
energy to neutrino energy, since leading neutrinos from $pp$ interactions
carry $\sim(3-4)\%$ of the nucleon energy in this energy range.

In Model B, the neutrino spectrum is assumed to be a broken power law,

$\frac{dN}{dE_{\nu}}\propto\begin{cases}E_{\nu}^{-\alpha_{\nu
1}}&E_{\nu}<E_{\nu,\rm bk}\\\ E_{\nu}^{-\alpha_{\nu 2}}&E_{\nu}>E_{\nu,\rm
bk}\end{cases}$ (B2)

where $\alpha_{\nu 1}=2.45$, $\alpha_{\nu 2}=2.85$, and $E_{\nu,\rm
bk}=0.04\times 3$ PeV. This model is based on modeling of the all-particle
cosmic-ray spectrum with a nucleon spectral cutoff energy of $\sim 3-4$ PeV
comparable to the knee energy (e.g., Hörandel, 2003; Gaisser, 2012; Lipari &
Vernetto, 2020), but we further take possible hardening of the cosmic-ray
spectral index due to the spatial inhomogeneity (by 0.15; Acero et al., 2016;
Lipari & Vernetto, 2018). The spectrum around the knee is believed to be
largely contributed by helium nuclei, so the break energy is likely to be
lower. This choice of the break energy and spectral index should be regarded
as the most optimistic case.


00footnotetext: $\bf{2010\ Mathematics\ Subject\ Classification\>}.$ 53D10,
53D15, 53D25, 53D35.
Key words and phrases: Lorentzian manifolds, perfect fluid spacetime, $GRW$
spacetimes, pseudosymmetric spacetimes

# Sufficient conditions for a pseudosymmetric spacetime to be a perfect fluid
spacetime

Peibiao Zhao  Uday Chand De  Bülent Ünal and Krishnendu De  Department of
Applied Mathematics, Nanjing University of Science and Technology, Nanjing
210094, P. R. China<EMAIL_ADDRESS>Department of Pure Mathematics,
University of Calcutta 35, Ballygaunge Circular Road Kolkata 700019, West
Bengal, India<EMAIL_ADDRESS>Department of Mathematics, Bilkent University,
Bilkent, 06800 Ankara, Turkey<EMAIL_ADDRESS>Department of Mathematics,
Kabi Sukanta Mahavidyalaya, The University of Burdwan. Bhadreswar, P.O.-Angus,
Hooghly, Pin 712221, West Bengal, India<EMAIL_ADDRESS>
###### Abstract

The aim of the present paper is to obtain the condition under which a
pseudosymmetric spacetime to be a perfect fluid spacetime. It is proven that a
pseudosymmetric generalized Robertson-Walker spacetime is a perfect fluid
spacetime. Moreover, we establish that a conformally flat pseudosymmetric
spacetime is a generalized Robertson-Walker spacetime. Next, it is shown that
a pseudosymmetric dust fluid with constant scalar curvature satisfying
Einstein’s field equations without cosmological constant is vacuum. Finally,
we construct a non-trivial example of pseudosymmetric spacetime.

##  1 Introduction

Let $M$ be a semi-Riemannian manifold of dimension $n\geq 2$ equipped with a
semi-Riemannian metric $g$ of signature $(m,p)$, where $m+p=n$. If $g$ is a
Lorentzian metric of signature $(1,n-1)$ or $(n-1,1)$, then $M$ endowed with
$g$ is said to be an $n$-dimensional Lorentzian manifold [23]. If
$M=-I\times_{\mathfrak{f}}\,\mathcal{M}$, where $I$ is an open interval of
$\mathbb{R}$ (set of real numbers), $\mathcal{M}$ is a Riemannian manifold of
dimension $n-1$ and $\mathfrak{f}>0$ denotes a smooth function, named as
warping function or scale factor, then $M$ is said to be a generalized
Robertson-Walker $(\text{briefly},GRW)$ spacetime ([1], [2]). Particularly, if
we suppose that $\mathcal{M}$ is a Riemannian manifold of dimension $3$ with
constant scalar curvature, then the $GRW$ spacetime becomes a Robertson-Walker
($RW$) spacetime. This infers that $GRW$ spacetimes are the natural extension
of $RW$ spacetimes. It is well-known that the Lorentzian Minkowski spacetime,
the Friedmann cosmological models, the static Einstein spacetime, the
Einstein-de Sitter spacetime, the de Sitter spacetime are included in the
$GRW$ spacetimes [27]. They have also many applications in inhomogeneous
spacetimes admitting an isotropic radiation. Several remarkable results on
$GRW$ spacetimes are investigated in ([3], [8], [9], [18], [20], [21]) and by
others.

A vector field $\zeta$ on an $n$-dimensional semi-Riemannian manifold $M$ is
said to be concircular if it satisfies $\nabla_{X}\zeta=\aleph\,X$ for some
smooth function $\aleph$, where $\nabla_{X}$ denotes the covariant derivative
along the smooth vector field $X$ of $M$. The notion of concircular vector
field on Riemannian manifold has been introduced by Failkow [15], and further
studied by [10] and others.

Chen [9] has given a local characterization of a Lorentzian manifold admitting
a timelike concircular vector field:
Theorem A. [9] A Lorentzian manifold $M$ of dimension $n\geq 3$ is a $GRW$
spacetime if and only if it admits a timelike concircular vector field such
that $X_{i,j}=\rho g_{ij}$ for some smooth function $\rho$.
The concept of torse-forming vector field is weaker than the concircular one,
but with an additional condition, the following local characterization of
$GRW$ spacetimes is possible.
Theorem B. [20] A Lorentzian manifold of dimension $n\geq 3$ is a $GRW$
spacetime if and only if it admits a unit timelike vector,
$u_{i,j}=\varphi(g_{ij}+u_{i}u_{j})$, that is also an eigenvector of the Ricci
tensor.
An $n$-dimensional spacetime $M$ with a non-vanishing Ricci tensor $R_{ij}$ is
said to be a perfect fluid spacetime if $R_{ij}$ satisfies the relation

$R_{ij}=\alpha g_{ij}+\beta A_{i}A_{j}$ (1.1)

for some smooth functions $\alpha$ and $\beta$ on $M$. Here $g$ denotes the
Lorentzian metric of $M$ and $A_{i}$ is the velocity vector such that
$g_{ij}A^{i}A^{j}=-1$ and $A_{i}=g_{ij}A^{j}$. It is noticed that a $RW$
spacetime is a perfect fluid spacetime ([23], Theorem 12.11). Every $GRW$
spacetime with $n=4$ is a perfect fluid spacetime if and only if it is a $RW$
spacetime. We cite ([4], [5], [13], [16], [19], [20], [22], [31]) and its
references for some deep results of the perfect fluid spacetimes.

During the investigation of conformally flat Riemannian manifolds of class
one, Sen and Chaki [28] found that the covariant derivative of the curvature
tensor $R$ of type $(0,4)$ satisfies

$R_{hijk,l}=2v_{l}R_{hijk}+v_{h}R_{lijk}+v_{i}R_{hljk}+v_{j}R_{hilk}+v_{k}R_{hijl},$
(1.2)

where ”, (comma)” denotes the covariant derivative with respect to the metric
tensor, $R_{hijk}$ are the components of the curvature tensor and $v_{l}$ is a
non-zero covector. Later, Chaki [6] named a Riemannian manifold to be
pseudosymmetric if the curvature tensor satisfies the condition (1.2). Then
Chaki and De [7] examined the Riemannian manifolds with the above condition.
If the covector $v_{l}$ vanishes, then the manifold becomes a locally
symmetric manifold in the sense of Cartan. An $n$-dimensional pseudosymmetric
manifold is denoted by $(PS)_{n}$. In $1992$, Deszcz [14] introduced the
notion of a pseudosymmetric manifold which is different from Chaki’s notion.
According to Deszcz, a pseudo-symmetric manifold is said to be pseudosymmetric
if $Riem\cdot Riem$ and $Q(g,Riem)$ are linearly dependent at each point of
the manifold, where $Riem$ and $Q$ denote the Riemannian curvature and Ricci
operator, respectively. Thus we have $Riem\cdot Riem=\mathfrak{h}Q(g,Riem)$
for some smooth function $\mathfrak{h}$. This paper deals with Chaki’s notion
of pseudosymmetric. Pseudosymmetric manifolds have been studied by several
authors such as ([11], [12], [17], [25]) and others.

In [26], Prvanović called such a manifold as a generalized recurrent manifold
and proved that the associated covector $v_{l}$ is gradient, that is,
irrotational.

On the other hand, a spacetime of general relativity is regarded as a
$4$-dimensional time-oriented Lorentzian manifold $(M,g)$. In general
relativity, the matter content of the spacetime is described by the energy
momentum tensor. The matter content is assumed to be fluid having density,
pressure and possessing dynamical and kinematical quantities like velocity,
acceleration, vorticity, shear and expansion. The fluid is called perfect
because of the absence of heat conduction terms and stress tensor
corresponding to viscosity. For a perfect fluid, the energy momentum tensor
$T_{ij}$ is given by

$T_{ij}=(\sigma+p)A_{i}A_{j}+pg_{ij},$ (1.3)

where $p$ and $\sigma$ are the isotropic pressure and energy density of the
perfect fluid, respectively [23].

The Einstein’s field equations without cosmological constant is given by

$R_{ij}-\frac{R}{2}g_{ij}=\kappa T_{ij},$ (1.4)

where $R_{ij}=R^{h}_{ijh}$ is the Ricci tensor, $R=g^{ij}R_{ij}$ denotes the
scalar curvature and $\kappa$ is the gravitational constant. The Einstein’s
field equations, in fact, connect the distribution of mass (represented by the
energy momentum tensor) with the curvature of the spacetime (represented by
the Einstein tensor).

In [21], Mantica, Molinari and De proved that a $GRW$ spacetime with
divergence free conformal curvature tensor ($C^{m}_{jkl,m}=0$) is a perfect
fluid spacetime.
Our aim is to improve the above result of [21]. We will prove the following:

######  Theorem 1.1.

Every pseudosymmetric $GRW$ spacetime is a perfect fluid spacetime.

The converse of Theorem 1.1 is not true in general. In this series, Mantica,
Molinari and De [21] have established the following theorem:
Theorem C. [21] A perfect fluid spacetime together with $div\,C=0$ is a
generalized Robertson-Walker spacetime with Einstein fiber, provided that the
velocity vector field of the perfect fluid is irrotational.

Motivated from Theorem C, we are going to prove the following theorem:

######  Theorem 1.2.

Every conformally flat $(PS)_{n}$ spacetime is a $GRW$ spacetime.

##  2 Pseudosymmetric $GRW$ spacetimes

This section is dedicated to give the condition for which a pseudosymmetric
spacetime to be a perfect fluid spacetime.
Proof of Theorem 1.1. Let us suppose that the pseudosymmetric spacetime is a
$GRW$ spacetime, then from Theorem B we have

$v_{i,j}=\varphi(g_{ij}+v_{i}v_{j}),$ (2.1)

where $\varphi$ is a scalar. Again the covariant derivative of equation (2.1)
gives

$v_{i,jk}=\varphi_{k}(g_{ij}+v_{i}v_{j})+\varphi\\{\varphi(g_{ik}+v_{i}v_{k})v_{j}+\varphi(g_{jk}+v_{j}v_{k})v_{i}\\},$
(2.2)

where $\varphi_{k}=\varphi_{,k}$ (the covariant derivative of the smooth
function $\varphi$). Interchanging $j$ and $k$ in the above equation and then
subtracting the foregoing equation from (2.1) and using Ricci identity, we
infer

$\displaystyle
v_{h}R^{h}_{ijk}=(\varphi_{k}g_{ij}-\varphi_{j}g_{ik})+v_{i}(\varphi_{k}v_{j}-\varphi_{j}v_{k})+\varphi^{2}(g_{ik}v_{j}-g_{ij}v_{k})$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,=(\varphi_{k}-\varphi^{2}v_{k})g_{ij}-(\varphi_{j}-\varphi^{2}v_{j})g_{ik}+v_{i}\varphi_{k}v_{j}-\varphi_{j}v_{i}v_{k}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,=w_{k}(g_{ij}+v_{i}v_{j})-w_{j}(g_{ik}+v_{i}v_{k}),$
(2.3)

where $w_{k}=\varphi_{k}-\varphi^{2}v_{k}$. Multiplying equation (2) with
$g^{ij}$, we get

$v_{h}R^{h}_{k}=(n-2)w_{k}-w_{j}v^{j}v_{k}.$ (2.4)

If $v_{k}$ is an eigenvector of the Ricci tensor, then $w_{k}=fv_{k}$ for some
scalar $f$. Hence equation (2.4) turns into

$v_{h}R^{h}_{k}=(n-1)fv_{k}.$

It is noted that the associated vector field $v_{i}$ is gradient [26]. Using
this fact in equation (1.2), we obtain

$\displaystyle
R_{hijk,lm}-R_{hijk,ml}=A_{hm}R_{lijk}+A_{im}R_{hljk}+A_{jm}R_{hilk}+A_{km}R_{hijl}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-A_{hl}R_{mijk}-A_{il}R_{hmjk}-A_{jl}R_{himk}-A_{kl}R_{hijm},$

where $A_{hm}=v_{h,m}-v_{h}v_{m}$. The foregoing equation together with
equation (2.1) and Ricci identity entail that

$\displaystyle
R_{pijk}R^{p}_{hlm}+R_{hpjk}R^{p}_{ilm}+R_{hipk}R^{p}_{jlm}+R_{hijp}R^{p}_{klm}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,=(\varphi
g_{hm}+(\varphi-1)v_{h}v_{m})R_{lijk}+(\varphi
g_{im}+(\varphi-1)v_{i}v_{m})R_{hljk}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+(\varphi
g_{jm}+(\varphi-1)v_{j}v_{m})R_{hilk}+(\varphi
g_{km}+(\varphi-1)v_{k}v_{m})R_{hijl}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-(\varphi
g_{hl}+(\varphi-1)v_{h}v_{l})R_{mijk}-(\varphi
g_{il}+(\varphi-1)v_{i}v_{l})R_{hmjk}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-(\varphi
g_{jl}+(\varphi-1)v_{j}v_{l})R_{himk}-(\varphi
g_{kl}+(\varphi-1)v_{k}v_{l})R_{hijm}.$ (2.5)

Multiplying equation (2) with $v^{h}$ and using (2), we conclude that

$\displaystyle L.H.S.$ $\displaystyle=$ $\displaystyle
R^{p}_{ijk}R_{phlm}v^{h}+v^{h}R_{hpjk}R^{p}_{ilm}+v^{h}R_{hipk}R^{p}_{jlm}+v^{h}R_{hijp}R^{p}_{klm}$
$\displaystyle=$ $\displaystyle-
R_{lijk}w_{m}-v_{p}R^{p}_{ijk}v_{l}w_{m}+w_{l}R_{mijk}+w_{l}v_{p}R^{p}_{ijk}v_{m}+w_{k}v_{j}v_{p}R^{p}_{ilm}$
$\displaystyle-
w_{j}v_{p}v_{k}R^{p}_{ilm}+w_{k}v_{i}v_{p}R^{p}_{jlm}-w_{p}R^{p}_{jlm}g_{ik}-w_{p}v_{i}v_{k}R^{p}_{jlm}$
$\displaystyle+w_{p}R^{p}_{klm}g_{ij}+w_{p}v_{i}v_{j}R^{p}_{klm}-w_{j}v_{i}v_{p}R^{p}_{klm}.$

Again, using equations (2) and $w_{k}=fv_{k}$ for some scalar $f$ in the
foregoing equation we have

$L.H.S.=-fv_{m}R_{lijk}+fv_{l}R_{mijk}-f^{2}v_{i}v_{k}(v_{m}g_{jl}-v_{l}g_{jm})+f^{2}(v_{m}g_{kl}-v_{l}g_{km})g_{ij}.$

Following the above process, the right-hand side of equation (2) assumes the
form

$\displaystyle
R.H.S.=v_{m}R_{lijk}-v_{l}R_{mijk}+f(\varphi-1)v_{i}v_{m}\\{v_{k}g_{lj}-v_{j}g_{lk}\\}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-f\varphi
g_{jm}v_{l}g_{ik}+f(\varphi-1)v_{j}v_{m}\\{v_{k}g_{il}-v_{l}g_{ik}\\}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+f(\varphi
g_{km}+(\varphi-1)v_{k}v_{m})(v_{l}g_{ij}-v_{j}g_{il})+f\varphi
v_{j}g_{il}g_{mk}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-f(\varphi-1)v_{i}v_{l}(v_{k}g_{mj}-v_{j}g_{mk})+f\varphi
g_{jl}v_{m}g_{ik}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-f(\varphi-1)v_{j}v_{l}\\{v_{k}g_{im}-v_{m}g_{ik}\\}-f\varphi
g_{kl}v_{m}g_{ij}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-f(\varphi-1)v_{k}v_{l}\\{v_{m}g_{ij}-v_{j}g_{im}\\}.$

The above equations infer that

$\displaystyle-
fv_{m}R_{lijk}+fv_{l}R_{mijk}-f^{2}v_{i}v_{k}(v_{m}g_{il}-v_{l}g_{jm})+f^{2}(v_{m}g_{kl}-v_{l}g_{km})g_{ij}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,=v_{m}R_{lijk}-v_{l}R_{mijk}+f(\varphi-1)v_{i}v_{m}\\{v_{k}g_{lj}-v_{j}g_{lk}\\}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-f\varphi
g_{jm}v_{l}g_{ik}+f(\varphi-1)v_{j}v_{m}\\{v_{k}g_{il}-v_{l}g_{ik}\\}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+f(\varphi
g_{km}+(\varphi-1)v_{k}v_{m})(v_{l}g_{ij}-v_{j}g_{il})+f\varphi
v_{j}g_{il}g_{mk}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-f(\varphi-1)v_{i}v_{l}(v_{k}g_{mj}-v_{j}g_{mk})+f\varphi
g_{jl}v_{m}g_{ik}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-f(\varphi-1)v_{j}v_{l}\\{v_{k}g_{im}-v_{m}g_{ik}\\}-f\varphi
g_{kl}v_{m}g_{ij}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-f(\varphi-1)v_{k}v_{l}\\{v_{m}g_{ij}-v_{j}g_{im}\\}.$
(2.6)

Multiplying equation (2) with $v^{l}$, and making use of equations (2) and
$w_{k}=fv_{k}$ we obtain

$\displaystyle-f^{2}v_{m}(v_{k}g_{ij}-v_{j}g_{ik})-R_{mijk}-f^{2}v_{i}v_{k}(v_{m}v_{j}+g_{jm})+f^{2}(v_{m}v_{k}+g_{mk})g_{ij}$
$\displaystyle=f^{2}v_{m}\\{v_{k}g_{ij}-v_{j}g_{ik}\\}+R_{mijk}+f(\varphi-1)v_{j}v_{m}\\{v_{k}v_{i}+g_{ik}\\}+f\varphi
g_{jm}g_{ik}$ $\displaystyle\,\,\,\,+f\\{\varphi
g_{jm}+(\varphi-1)v_{j}v_{m}\\}(v_{k}v_{i}+g_{ki})+f(\varphi-1)v_{i}(v_{k}g_{mj}-v_{j}g_{mk})$
$\displaystyle\,\,\,\,+f\varphi v_{j}v_{i}g_{mk}+f\varphi
v_{j}v_{m}g_{ik}+(\varphi-1)v_{j}\\{v_{k}g_{im}-v_{m}g_{ik}\\}$
$\displaystyle\,\,\,\,-f\varphi
v_{k}v_{m}g_{ij}+f(\varphi-1)v_{k}\\{v_{m}g_{ij}-v_{j}g_{im}\\}.$

Again multiplying the above equation with $g^{ij}$, we find

$R_{km}=\alpha g_{km}+\beta v_{k}v_{m},$ (2.7)

where $\beta=\frac{1}{2}[-f^{2}+f\\{(2\varphi-1)+n(\varphi-1)\\}]$ and
$\alpha=\frac{1}{2}[-nf^{2}+f(\varphi-1)]$. The last equation together with
(1.1) show that the pseudosymmetric $GRW$ spacetime is a perfect fluid
spacetime. This finishes the proof of the Theorem.

######  Remark.

In [21], Mantica, Molinari and first author proved that a perfect fluid
spacetime in dimension $n\geq 4$, with differentiable equation of state
$p=p(\sigma)$, $p+\sigma\neq 0$, and with null divergence of the Weyl
conformal curvature tensor $(C^{m}_{jkl,m}=0)$, is a $GRW$ spacetime. In
Theorem 1.1, we gave an affirmative answer of the question ”Under which
condition a $GRW$ spacetime to be a perfect fluid spacetime?”

In consequence of equations (1.3), (1.4) and (2.7), we notice that

$\left(\alpha-\frac{R}{2}-\kappa
p\right)g_{ij}+(\beta-\kappa(p+\sigma))v_{i}v_{j}=0,$

which gives

$\kappa\sigma=(n-1)\alpha\,\,\,\,\,\text{and}\,\,\,\,(n-1)p\kappa=(n-1)\alpha-\frac{(n-2)(n\alpha-\beta)}{2}.$

These reflect that the equation of state assumes the form

$\frac{p}{\sigma}=\frac{1}{n-1}-\frac{(n-2)(n\alpha-\beta)}{2\alpha(n-1)^{2}}.$
(2.8)

For $n=4$, the above equation with (2.7) take the form

$\frac{p}{\sigma}=\frac{1}{3}+\frac{15f+2\varphi-1}{\varphi-4f-1}.$

If $2\varphi=1-15f$, then the above equation becomes $3p=\sigma$, that is, the
equation of state represents the radiation era in the evolution of the
universe [30]. Now, we state our finding as:

######  Corollary 2.1.

Let a pseudosymmetric $GRW$ spacetime satisfy the Einstein’s field equations
without cosmological constant, then the equation of state is given by (2.8).
Also, if $n=4$ and $2\varphi=1-15f$, the matter of the spacetime represents
the radiation era.

##  3 Conformally flat $(PS)_{n}$ spacetimes

This section deals with the study of pseudosymmetric conformally flat
spacetimes.
Proof of Theorem 1.2. It is well-known that in a conformally flat $(PS)_{n}$
spacetime, the Ricci tensor $R_{ij}$ is of the form

$R_{ij}=\frac{R-t}{n-1}g_{ij}+\frac{nt-R}{(n-1)v_{p}v^{p}}v_{i}v_{j},$

where $R$ denotes the scalar curvature and $t$ is a scalar [29]. The above
equation can be rewritten as

$R_{ij}=\alpha g_{ij}+\beta\lambda_{i}\lambda_{j},$ (3.1)

where $\alpha=\frac{R-t}{n-1}$ and $\beta=\frac{nt-R}{(n-1)v_{p}v^{p}}$ are
two scalars and $\lambda_{i}=\frac{v_{i}}{\sqrt{-v_{j}v^{j}}}$ is a unit
timelike vector. Equation (3.1) reflects that the conformally flat $(PS)_{n}$
spacetimes are the perfect fluid spacetimes.

According to our hypothesis, the conformal curvature tensor $C$ of our ambient
space vanishes and therefore $div\,C=0$, where $div$ denotes the divergence.
In [21], Mantica, Molinari and De proved that a perfect fluid spacetime with
$div\,C=0$ is a $GRW$ spacetime. Thus, by considering equation (3.1) and the
above discussions, we obtain the required result.

##  4 Viscous fluid

We consider the pseudosymmetric spacetimes with constant scalar curvature.
From equation (1.2) we get

$R_{ij,l}=2v_{l}R_{ij}+v_{i}R_{lj}+v_{j}R_{il}+v^{k}R_{lijk}+v^{k}R_{kijl}.$

Considering an orthonormal frame field and contracting the above equation
entails that $R_{,l}=2v_{l}R+4R_{lh}v^{h}$. Since the scalar curvature $R$ is
constant, therefore we get $R_{lh}v^{h}=-\frac{R}{2}v_{l}$, which implies that
$-\frac{R}{2}$ is an eigenvalue of the Ricci tensor $R_{lh}$ corresponding to
the eigenvector $v_{l}$.

In viscous fluid, the energy momentum tensor [24] is given by

$T_{ij}=(p+\sigma)v_{i}v_{j}+pg_{ij}+D_{ij},$ (4.1)

where $D_{ij}$ is a $(0,2)$-type symmetric tensor, named as the anisotropic
pressure of the fluid. Also $D_{ij}$ satisfies $D_{ij}v^{j}=0$ and trace of
$D=0$. Suppose the spacetime under consideration obeys Einstein’s field
equation without cosmological constant, that is equation (1.4) is satisfied.
Then we have

$R_{ij}-\frac{R}{2}g_{ij}=\kappa[(p+\sigma)v_{i}v_{j}+pg_{ij}+D_{ij}].$ (4.2)

Multiplying equation (4.2) with $g^{ij}$ we find

$\frac{2-n}{2\kappa}R=(n-1)p-\sigma,$ (4.3)

since trace of $D=0.$ Now, multiplying equation (4.2) with $v^{j}$ we infer

$-Rv_{i}=-\kappa\sigma v_{i}\implies\kappa\sigma=R,$ (4.4)

since $D_{ij}v^{j}=0$ and $R_{ij}v^{j}=-\frac{R}{2}v_{i}$. Equations (4.3) and
(4.4) entail that

$p=\frac{4-n}{2(n-1)}\sigma,$ (4.5)

which means that the fluid is isentropic [18]. Thus, we can write our result
as:

######  Theorem 4.1.

Let the pseudosymmetric spacetime with the constant scalar curvature satisfy
the Einstein’s field equations. If the spacetime is filled with viscous fluid,
then the fluid is isentropic and the equation of state is given by (4.5).

######  Remark.

For $n=4$, equation (4.5) reflects that $p=0$, that is the isotropic pressure
of the viscous fluid is zero and the spacetime is filled with dust matter.

In case of dust or pressure-less fluid, the energy momentum tensor assumes the
form

$T_{ij}=\sigma v_{i}v_{j}.$ (4.6)

In consequence of equation (1.4), equation (4.6) takes the form

$R_{ij}-\frac{R}{2}g_{ij}=\kappa\sigma v_{i}v_{j}.$ (4.7)

Let us consider an orthonormal frame field and contracting the above equation
for $i$ and $j$, we get

$\frac{2-n}{2}R=-\kappa\sigma.$ (4.8)

Again, multiplying equation (4.7) with $v^{j}$ we lead to

$-Rv_{i}=-\kappa\sigma v_{i}\implies\kappa\sigma=R,$ (4.9)

since $R_{ij}v^{j}=-\frac{R}{2}v_{i}$. From equations (4.8) and (4.9) we infer
that $\sigma=0$ and hence equation (4.6) reduces to $T_{ij}=0$, which implies
that the spacetime is vacuum.

######  Theorem 4.2.

Let a pseudosymmetric spacetime with the constant scalar curvature satisfy the
Einstein’s field equation without a cosmological constant. If the spacetime is
filled with dust fluid, then it is vacuum.

##  5 An Example of Pseudosymmetric spacetime

In this section, we construct a non-trivial example to prove the existence of
a pseudosymmetric spacetime whose associated vector is irrotational.

Let us consider a Lorentzian metric $g$ on $\mathbb{R}^{4}$ by

$ds^{2}=g_{ij}dx^{i}dx^{j}=(dx^{1})^{2}+(x^{1})^{2}(dx^{2})^{2}+(x^{2})^{2}(dx^{3})^{2}-(dx^{4})^{2},$

where $(x^{1},x^{2},x^{3},x^{4})$ are the standard coordinate of
$\mathbb{R}^{4}$ and $i,j=1,2,3,4$. The only non-vanishing components of the
Christoffel symbols and curvature tensor are

$\Gamma^{1}_{22}=-x^{1},\,\,\,\Gamma^{2}_{33}=-\frac{x^{2}}{(x^{1})^{2}},\,\,\,\Gamma^{2}_{12}=\frac{1}{x^{1}},\,\,\,\Gamma^{3}_{23}=\frac{1}{x^{2}}\,\,\text{and}\,\,R_{1332}=-\frac{x^{2}}{x^{1}}.$

From equation
$R_{hijk,l}=2v_{l}R_{hijk}+v_{h}R_{lijk}+v_{i}R_{hljk}+v_{j}R_{hilk}+v_{k}R_{hijl}$
we have

$R_{1332,1}=2v_{1}R_{1332}+v_{1}R_{1332}+v_{3}R_{1132}+v_{3}R_{1312}+v_{2}R_{1331},$
(5.1)
$R_{1332,2}=2v_{2}R_{1332}+v_{1}R_{2332}+v_{3}R_{1232}+v_{3}R_{1322}+v_{2}R_{1332},$
(5.2)
$R_{1332,3}=2v_{3}R_{1332}+v_{1}R_{3332}+v_{3}R_{1332}+v_{3}R_{1332}+v_{2}R_{1333},$
(5.3)
$R_{1332,4}=2v_{4}R_{1332}+v_{1}R_{4332}+v_{3}R_{1432}+v_{3}R_{1342}+v_{2}R_{1334}.$
(5.4)

We choose the covariant vector $v_{i}$ as follows:

$v_{i}=\left\\{\begin{array}[]{rcl}-\frac{2}{3x^{1}},&\mbox{for}&i=1\\\ \\\
-\frac{1}{3x^{2}},&\mbox{for}&i=2\\\ \\\
0,&\mbox{for}&i=3,4.\end{array}\right.$

It can be easily shown that

$\displaystyle
v_{1,4}=v_{4,1}=0,\,\,v_{1,2}=v_{2,1}=\frac{1}{3x^{1}x^{2}},\,\,v_{1,3}=v_{3,1}=0,$
$\displaystyle\,\,\,\,\,v_{2,4}=v_{4,2}=0,\,\,v_{2,3}=v_{3,2}=0,\,\,v_{3,4}=v_{4,3}=0.$

From the above discussions, we can conclude that the equations (5.1)-(5.4) are
true. Thus, the manifold under consideration is a pseudosymmetric spacetime.
Moreover, equation (2.2) shows that the covector ($1$-form) $v$ is
irrotational (closed).

## References

  * [1] L. Alias, A. Romero and M. Sánchez, Uniqueness of complete spacelike hypersurfaces of constant mean curvature in generalized Robertson-Walker spacetimes, Gen. Relativ. Gravit. 27 (1995), 71-84.
  * [2] L. Alias, A. Romero and M. Sánchez, Compact spacelike hypersurfaces of constant mean curvature in generalized Robertson-Walker spacetimes, in Geometry and Topology of Submanifolds VII, River Edge NJ, USA, World Scientific Publ., 1995, 67-70.
  * [3] K. Arslan, R. Deszcz, R. Ezentaş, M. Hotloś and C. Murathan, On generalized Robertson-Walker spacetimes satisfying some curvature condition, Turk. J. Math. 38 (2) (2014), 353-373.
  * [4] A. M. Blaga, Solitons and geometrical structures in a perfect fluid spacetime, Rocky Mountain J. Math. 50(1) (2020), 41-53.
  * [5] A. M. Blaga, On harmonicity and Miao-Tam critical metrics in a perfect fluid spacetime, Bol. Soc. Mat. Mex. (3) 26(3) (2020), 1289-1299.
  * [6] M. C. Chaki, On pseudo symmetric manifolds, An. Ştiinţ. Univ. Al. I. Cuza Iaşi. Mat. (N.S.) 33 (1987), 53-58.
  * [7] M. C. Chaki and U. C. De, On pseudo symmetric spaces, Acta Math. Hungarica 54 (3-4) (1989), 185-190.
  * [8] S. K. Chaubey, Y. J. Suh and U. C. De, Characterizations of the Lorentzian manifolds admitting a type of semi-symmetric metric connection, Anal. Math. Phys. 10(4) (2020), Paper No. 61, 15 pp.
  * [9] B. Y. Chen, A simple characterization of generalized Robertson-Walker spacetimes, Gen. Relativ. Gravit. 46 (2014), 1833 (5 pages).
  * [10] B. Y. Chen, Concircular vector fields and pseudo-Kaehler manifolds, Kragujevac J. Math. 40 (1) (2016), 7-14.
  * [11] U. C. De, Y. J. Suh, S. K. Chaubey and S. Shenawy, On pseudo $H$-symmetric Lorentzian manifolds with applications to relativity, Filomat 34 (10) (2020), 3287-3297.
  * [12] U. C. De and A. K. Gazi, On almost pseudo conformally symmetric manifolds, Demonstr. Math. 42(4) (2009), 507-520.
  * [13] U. C. De, S. K. Chaubey and S. Shenawy, Perfect fluid spacetimes and Yamabe solitons, J. Math. Phys. 62, 032501 (2021); doi: 10.1063/5.0033967.
  * [14] R. Deszcz, On pseudosymmetric spaces, Bull. Belg. Math. Soc., Series A 44 (1992), 1-34.
  * [15] A. Fialkow, Conformal geodesics, Trans. Amer. Math. Soc. 45(3) (1939), 443-473.
  * [16] S. Mallick, U. C. De and Y. J. Suh, Spacetimes with different forms of energy-momentum tensor, J. Geom. Phys. 151 (2020), 103622 (8 pp).
  * [17] C. A. Mantica and Y. J. Suh, Pseudo $Z$-symmetric Riemannian manifolds with harmonic curvature tensors, Int. J. Geom. Meth. Mod. Phys. 9 (2012), 1250004 (21 pages).
  * [18] C. A. Mantica, Y. J. Suh and U. C. De, A note on generalized Robertson-Walker space-times, Int. J. Geom. Meth. Mod. Phys. 13 (6) (2016), 1650079 (9 pages).
  * [19] C. A. Mantica, U. C. De, Y. J. Suh and L. G. Molinari, Perfect fluid spacetimes with harmonic generalized curvature tensor, Osaka J. Math. 56 (2019), 173-182.
  * [20] C. A. Mantica and L. G. Molinari, Generalized Robertson-Walker spacetimes-A survey, Int. J. Geom. Meth. Mod. Phys. 14 (3) (2017), 1730001 (27 pages).
  * [21] C. A. Mantica, L. G. Molinari and U. C. De, A condition for a perfect-fluid space-time to be a generalized Robertson-Walker space-time, J. Math. Phys. 57 (2) (2016), 022508.
  * [22] C. A. Mantica, L. G. Molinari, Y. J. Suh and S. Shenawy, Perfect-fluid, generalized Robertson-Walker space-times, and Grays decomposition, J. Math. Phys. 60 (2019), 052506.
  * [23] B. ONeill, Semi-Riemannian Geometry with Applications to Relativity, Academic Press, New York, 1983.
  * [24] M. Novello and M. J. Reboucas, The stability of a rotating universe, The Astrophysics Journal, 225 (1978), 719-724.
  * [25] F. Ozen and S. Altay, Weakly and pseudo-symmetric Riemannian spaces, Indian J. Pure Appl. Math. 33 (2002), 1477-1488.
  * [26] M. Prvanović, Generalized recurrent Riemannian manifolds, Ann. St. Univ. ”AL.I.Cuza” Iasi 38 (1992), 423-434.
  * [27] M. Sánchez, On the geometry of generalized Robertson-Walker spacetimes: geodesics, Gen. Relativ. Gravit. 30 (1998), 915-932.
  * [28] R. N. Sen and M. C. Chaki, On curvature restrictions of a certain kind of conformally flat Riemannian space of class one, Proc. Nat. Inst. Sci. India Part A 33 (1967), 100-102.
  * [29] M. Tarafdar, On conformally flat pseudosymmetric manifolds, An. Ştiinţ. Univ. Al. I. Cuza Iaşi. Mat. (N.S.) 41 (1995), 237-241.
  * [30] S. Weinberg, Gravitation and cosmology: Principles and applications of the general theory of gravitation (Wiley, New York, 1972).
  * [31] F.O. Zengin and A.Y. Taşci, Spacetimes admitting the $Z$-symmetric tensor, Quaestiones Mathematicae 2020: 1-11. https://doi.org/10.2989/16073606.2020.1816587.