# Big-Bounce in projectively invariant Nieh-Yan models: the Bianchi I case

Flavio Bombacigno  Departament de Física Teòrica and IFIC, Centro Mixto
Universitat de València - CSIC, Universitat de València, Burjassot 46100,
València, Spain
E-mail<EMAIL_ADDRESS>Simon Boudet  Dipartimento di Fisica,
Università di Trento,
Via Sommarive 14, I-38123 Povo (TN), Italy
and
Trento Institute for Fundamental Physics and Applications (TIFPA)-INFN,
Via Sommarive 14, I-38123 Povo (TN), Italy
E-mail<EMAIL_ADDRESS>Gonzalo J. Olmo  Departament de Física Teòrica
and IFIC, Centro Mixto Universitat de València - CSIC, Universitat de
València, Burjassot 46100, València, Spain
<EMAIL_ADDRESS>Giovanni Montani  Physics Department, “Sapienza”
University of Rome,
P.le Aldo Moro 5, 00185 (Roma), Italy
and
ENEA, Fusion and Nuclear Safety Department,
C. R. Frascati, Via E. Fermi 45, 00044 Frascati (Roma), Italy
E-mail<EMAIL_ADDRESS>
###### Abstract

We show that the Nieh-Yan topological invariant breaks projective symmetry and
loses its topological character in presence of non vanishing nonmetricity. The
notion of the Nieh-Yan topological invariant is then extended to the generic
metric-affine case, defining a generalized Nieh-Yan term, which allows to
recover topologicity and projective invariance, independently. As a concrete
example a class of modified theories of gravity is considered and its
dynamical properties are investigated in a cosmological setting. In
particular, bouncing cosmological solutions in Bianchi I models are derived.
Finite time singularities affecting these solutions are analysed, showing that
the geodesic completeness and the regular behavior of scalar perturbations in
these space-times are not spoiled.

###### keywords:

Nieh-Yan, metric-affine, Bianchi I, big-bounce

##  1 Introduction

The theory of General Relativity (GR) [1, 2] relies on the geometric
interpretation of the gravitational field, described in terms of a metric
tensor and a connection on a pseudo-Riemannian manifold. Both GR and many
alternative theories of gravity are based on a metric formulation, in which
the connection is given by the symmetric and metric compatible Levi-Civita
connection, which is completely determined by the metric and its derivatives.
An alternative formulation for geometric theories of gravity consists in
adopting the metric-affine paradigm, in which the metric tensor and the
connection are considered as independent variables. In this approach, symmetry
and metric compatibility of the connection are not imposed a priori, resulting
in the presence of torsion and nonmetricity, respectively. Well known examples
of metric-affine theories are Ricci based gravity [3, 4], Palatini $f(R)$
theory [5], quadratic gravity [6], Born-Infeld-type models [7], general
teleparallel models [8], generalized hybrid metric-Palatini gravity [9, 10,
11, 12, 13] and metric-affine extension of higher order theories [14, 15, 16,
17].
The metric-affine approach plays a crucial role also in one of the current
attempts to quantize gravity, i.e. loop quantum gravity (LQG) [18, 19], where
GR is reformulated in terms of a gauge $SU(2)$ connection (Ashtekar-Barbero-
Immirzi connection) and its conjugate momentum, the densitized triad [20, 21,
22, 23]. This formulation, indeed, can be derived [24] by including an
additional contribution to the first order (Palatini) action of GR, namely the
Nieh-Yan (NY) topological invariant [25, 26] (the Holst term [27] can be used
as well). The NY term was discovered in the context of Riemann-Cartan theory
(where nonmetricity is set to zero) and its main property is topologicity: it
reduces to a boundary term without affecting the field equations at all. This
additional term is driven by the so called Immirzi parameter $\beta$ [28, 29],
which concurs in the definition of the Ashtekar variables and is related to a
quantization ambiguity [30]. Attempts to address this issue led to the
proposal of considering the Immirzi parameter as a new fundamental field [31,
32, 33], an idea that has been later developed within several different
contexts [34, 35, 36, 37, 38, 39, 40, 41, 31, 33, 42, 43, 44, 45, 46, 33, 40].
The promotion of such constant parameter to a dynamical field is usually
pursued “by hand”, substituting $\beta\rightarrow\beta(x)$ in the Lagrangian
and possibly adding a potential term $V(\beta)$.
More recently, beside LQG the NY term has been studied in the context of
teleparallel gravity [47] and in condensed matter physics [48, 49, 50, 51].
Another important property we will focus on, is projective invariance [52,
53], which has recently been shown to be of crucial importance in metric-
affine theories since the breaking of this symmetry can give rise to dynamical
instabilities [54]. In this regard, we want to stress that the NY term breaks
this symmetry. This feature has always been neglected in literature and a
revision of previous formulations seems necessary. Moreover, as will be shown
in the following, the topological character of the NY term is also lost when
nonmetricity is included.
The approach followed in this note is grounded on the choice of recovering
these features from the very beginning in the action, without imposing any
restriction on the affine connection. After a formal discussion, we will
implement the gravitational model in a cosmological setting. In particular, we
investigate Bianchi I models [55, 56, 57], focusing our attention on the
emergence of a classical bouncing cosmology [58, 59, 60, 61, 62, 42, 63, 64].

##  2 The role of nonmetricity in the Nieh-Yan term

In Einstein-Cartan theory the NY term [25, 26] is explicitly defined as

$\displaystyle NY$

where

$R\indices{{}^{\rho}_{\mu\sigma\nu}}=\partial_{\sigma}\Gamma\indices{{}^{\rho}_{\mu\nu}}-\partial_{\nu}\Gamma\indices{{}^{\rho}_{\mu\sigma}}+\Gamma\indices{{}^{\rho}_{\tau\sigma}}\Gamma\indices{{}^{\tau}_{\mu\nu}}-\Gamma\indices{{}^{\rho}_{\tau\nu}}\Gamma\indices{{}^{\tau}_{\mu\sigma}},$
(1)

is the Riemann tensor built with the independent connection and

$T\indices{{}^{\mu}_{\nu\rho}}=\Gamma\indices{{}^{\mu}_{\nu\rho}}-\Gamma\indices{{}^{\mu}_{\rho\nu}}$
(2)

is the torsion tensor. The starting point of our discussion is the observation
that for a non-vanishing nonmetricity tensor
$Q_{\mu\nu\rho}=-\nabla_{\mu}g_{\nu\rho}$, the NY term (LABEL:NY0) is spoilt
of its topological character. Indeed, extracting the nonriemmanian part of the
Riemann tensor leads to[65]

$\displaystyle NY$ $\displaystyle=-\frac{1}{2}\bar{\nabla}\cdot
S-\frac{1}{2}\varepsilon^{\mu\nu\rho\sigma}T\indices{{}^{\lambda}_{\mu\nu}}Q\indices{{}_{\rho\sigma\lambda}},$
(3)

where $\bar{\nabla}_{\mu}$ is built with the Levi-Civita connection
$\bar{\Gamma}\indices{{}^{\mu}_{\nu\rho}}$ and

$S_{\mu}\equiv\varepsilon_{\mu\nu\rho\sigma}T^{\nu\rho\sigma}.$ (4)

Therefore, when $Q_{\rho\mu\nu}\neq 0$, the Nieh-Yan term does not simply
reduce to the divergence of a vector, and the appearance of nonmetricity
spoils the topologicity. Let us now consider the behavior of (LABEL:NY0) under
projective transformations of the connection, namely

$\tilde{\Gamma}\indices{{}^{\rho}_{\mu\nu}}=\Gamma\indices{{}^{\rho}_{\mu\nu}}+\delta\indices{{}^{\rho}_{\mu}}\xi_{\nu},$
(5)

It can easily be shown that (LABEL:NY0) is also not invariant under projective
transformations, since

$\frac{1}{4}\varepsilon^{\mu\nu\rho\sigma}\tilde{T}\indices{{}^{\lambda}_{\mu\nu}}\tilde{T}_{\lambda\rho\sigma}-\frac{1}{4}\varepsilon^{\mu\nu\rho\sigma}T\indices{{}^{\lambda}_{\mu\nu}}T_{\lambda\rho\sigma}=-S_{\mu}\xi^{\mu}.$
(6)

Now, by looking at (3), we point out that a newly topological Nieh-Yan term
can be recovered by simply setting

$NY^{*}\equiv
NY+\frac{1}{2}\varepsilon^{\mu\nu\rho\sigma}T\indices{{}^{\lambda}_{\mu\nu}}Q\indices{{}_{\rho\sigma\lambda}}.$
(7)

We note that projective invariance is now enclosed as well, since

$\frac{1}{2}\varepsilon^{\mu\nu\rho\sigma}\tilde{T}\indices{{}^{\lambda}_{\mu\nu}}\tilde{Q}\indices{{}_{\rho\sigma\lambda}}-\frac{1}{2}\varepsilon^{\mu\nu\rho\sigma}T\indices{{}^{\lambda}_{\mu\nu}}Q\indices{{}_{\rho\sigma\lambda}}=+S_{\mu}\xi^{\mu},$
(8)

which exactly cancels out (6). We stress, however, that projective invariance
is not strictly related to topologicity, and suitable generalizations of (7)
breaking up only with the latter can be actually formulated. Let us consider,
for instance, the following modified Nieh-Yan term

$NY_{gen}\equiv\frac{1}{2}\varepsilon^{\mu\nu\rho\sigma}\left(\frac{\lambda_{1}}{2}T\indices{{}^{\lambda}_{\mu\nu}}T_{\lambda\rho\sigma}+\lambda_{2}\,T\indices{{}^{\lambda}_{\mu\nu}}Q\indices{{}_{\rho\sigma\lambda}}-R_{\mu\nu\rho\sigma}\right),$
(9)

where we introduced the real parameters $\lambda_{1},\,\lambda_{2}$. In this
case the term (9) transforms under a projective transformation as

$NY_{gen}\rightarrow NY_{gen}-(\lambda_{1}-\lambda_{2})\xi^{\mu}S_{\mu},$ (10)

so that by setting $\lambda_{1}=\lambda_{2}$ we can recover again projective
invariance, despite topologicity being in general violated if
$\lambda_{1}=\lambda_{2}\neq 1$, since

$NY_{gen}=-\frac{1}{2}\bar{\nabla}\cdot
S+\frac{(\lambda_{1}-1)}{4}\varepsilon^{\mu\nu\rho\sigma}T\indices{{}^{\lambda}_{\mu\nu}}T_{\lambda\rho\sigma}+\frac{(\lambda_{2}-1)}{2}\varepsilon^{\mu\nu\rho\sigma}T\indices{{}^{\lambda}_{\mu\nu}}Q\indices{{}_{\rho\sigma\lambda}}.$
(11)

In the following, we will consider the general form (9), which by a suitable
choice of the parameters $\lambda_{1,2}$ can reproduce all known actions
usually studied in literature, as the Holst ($\lambda_{1}=\lambda_{2}=0$) or
the standard Nieh-Yan (3) ($\lambda_{1}=1,\,\lambda_{2}=0$) terms.

##  3 Generalized Nieh-Yan models

As a specific gravitational model featuring the generalized NY term we
consider an action defined by a general function of two arguments, the Ricci
scalar and the generalized NY term (9):

$S_{g}=\frac{1}{2\kappa}\int d^{4}x\sqrt{-g}\,F(R,NY_{gen}),$ (12)

Now, performing the transformation to the Jordan frame leads to the scalar
tensor representation

$S_{g}=\frac{1}{2\kappa}\int d^{4}x\sqrt{-g}\left(\phi R+\beta
NY_{gen}-W(\phi,\beta)\right),$ (13)

with $\phi\equiv\frac{\partial F}{\partial R},\,\beta\equiv\frac{\partial
F}{\partial NY_{gen}}$ and $W\equiv\phi R(\phi,\beta)+\beta
NY_{gen}(\phi,\beta)-F(\phi,\beta)$. The scalar field $\beta$ can be
identified with the Immirzi field, which acquires in this way a dynamical
character without the need of introducing this feature by hand in the action.
Moreover, this formulation offers a viable mechanism to produce an interaction
term $W(\phi,\beta)$ as well. Now, the field equation for the connection are
obtained varying (13) with respect to $\Gamma\indices{{}^{\mu}_{\nu\rho}}$.
For the full set of equations the reader may cosnult [65], while here we are
interest in the following contraction

$\delta^{\lambda}_{\nu}\frac{\delta
S_{g}}{\delta\Gamma\indices{{}^{\lambda}_{\nu\mu}}}=0$ (14)

which leads to

$(\lambda_{1}-\lambda_{2})S^{\mu}=0.$ (15)

This implies that the features of the solutions depend on the parameters
$\lambda_{1}$ and $\lambda_{2}$, and when projective invariance is broken
($\lambda_{1}\neq\lambda_{2}$) one is compelled to set $S_{\mu}=0$. In this
case, (9) can be re-expressed as

$NY_{gen}=-\frac{1}{2}\bar{\nabla}\cdot
S-\frac{(1-\lambda_{1})}{3}S^{\mu}T\indices{{}^{\rho}_{\mu\rho}}-\frac{(1-\lambda_{2})}{2}S^{\mu}Q\indices{{}^{\rho}_{\mu\rho}},$
(16)

implying that the generalized Nieh-Yan term (9) is identically vanishing on
half-shell. In other words, terms violating projective invariance are harmless
along the dynamics. This can be further appreciated deriving the effective
scalar tensor action stemming from (13), when the solutions of the full set of
connection field equations are plugged in it. Explicit calculations (see [65]
for details) lead to

$S=\frac{1}{2\kappa}\int
d^{4}x\sqrt{-g}\left(\phi\bar{R}+\frac{3}{2\phi}\bar{\nabla}_{\mu}\phi\bar{\nabla}^{\mu}\phi-W(\phi,\beta)\right),$
(17)

where $\bar{R}$ is the Ricci scalar of the Levi-CIvita connection. This
resembles the form of a Palatini $f(R)$ theory, with a potential depending on
the Immirzi field as well. The equation for the latter, i.e.

$\frac{\partial W(\phi,\beta)}{\partial\beta}=0,$ (18)

fixes its form in terms of the remaining scalar field: $\beta=\beta(\phi)$.
Then, using the trace of the metric field equations, variation of (17) with
respect to $\phi$ results in the usual structural equation featuring Palatini
$f(R)$ theories [5], i.e.

$\left[2W(\phi,\beta)-\phi\frac{\partial
W(\phi,\beta)}{\partial\phi}\right]_{\beta=\beta(\phi)}=\kappa T,$ (19)

where $T$ is the trace of the stress energy tensor of matter. This implies
that the dynamics of the scalaron $\phi$ is frozen as well, and completely
determined by $T$. Conditions (18) and (19) then guarantee that the scalar
fields $\phi,\;\beta$ are not propagating degrees of freedom, and reduce to
constants in vacuum, where the theory is stable and the breaking of projective
invariance does not lead to ghost instabilities, in contrast to [54]. If
$\lambda_{1}=\lambda_{2}\equiv\lambda$, instead, the projective invariance of
the model can be used to get rid of one vector degree of freedom, which can be
set to zero properly choosing the vector $\xi_{\mu}$. A convenient choice
consists in setting $\xi_{\mu}=-\frac{1}{2}Q\indices{{}^{\rho}_{\mu\rho}}$,
which allows to deal only with torsion in the connection field equations. The
effective action stemming from (13) then reads (see [65])

$S=\frac{1}{2\kappa}\int
d^{4}x\sqrt{-g}\left(\phi\bar{R}+\frac{3}{2\phi}\bar{\nabla}_{\mu}\phi\bar{\nabla}^{\mu}\phi-\frac{3\phi}{2}\frac{1}{\phi^{2\lambda}+(1-\lambda)^{2}\psi^{2}}\bar{\nabla}_{\mu}\psi\bar{\nabla}^{\mu}\psi-V(\phi,\psi)\right),$
(20)

where we used the transformation $\psi\equiv\beta\phi^{\lambda-1}$ and
redefined the potential as $V(\phi,\psi)=W(\phi,\psi\phi^{1-\lambda})$. In
general, the Immirzi field is expected to be a well-behaved dynamical degree
of freedom, since in the Einstein frame action, defined by the conformal
rescaling $\tilde{g}_{\mu\nu}=\phi\,g_{\mu\nu}$, the kinetic term for the
Immirzi field takes the form

$-\frac{3}{2}\frac{\tilde{g}^{\mu\nu}\nabla_{\mu}\psi\nabla_{\nu}\psi}{\phi^{2\lambda}+(1-\lambda)^{2}\psi^{2}}.$
(21)

Since $\phi^{2\lambda}+(1-\lambda)^{2}\psi^{2}>0$ for every value of
$\phi,\,\psi$ and $\lambda$, (21) has always the correct sign and no ghost
instability arise [66]. Let us end this section with a remark on how previous
results with vanishing nonmetricity can properly be recovered. In particular,
the Riemann-Cartan structure of [32, 67, 68, 69, 46, 33, 37, 36] can be
replicated by inserting in (12) the condition of vanishing nonmetricity with a
Lagrange multiplier, i.e. adding to the Lagrangian a term
$l^{\rho\mu\nu}Q_{\rho\mu\nu}$, with $l^{\rho\mu\nu}=l^{\rho\nu\mu}$. Then,
results of [32, 67, 68, 69, 46, 33, 37, 36] are simply obtained by setting
$\lambda_{1}=1$. The fact that the usual Einstein-Cartan NY invariant and
related models are recovered in this way, supports the correctness of our
generalization of the NY term, with respect to other possible generalizations.

##  4 Big bounce in Bianchi I cosmology

In this section we consider dynamical models, i.e. those described by (20),
which are characterized by a dynamical Immirzi field and look for cosmological
solutions in Bianchi I spacetimes. In particular, we will be interested in
obtaining solutions characterized by a bouncing behavior for the universe
volume, thanks to which the big bang singularity is regularized in favour of a
big bounce scenario.

Let us start from the equations of motion for the metric and scalar fields
which are given by

$\displaystyle\bar{G}_{\mu\nu}$
$\displaystyle=\frac{\kappa}{\phi}T_{\mu\nu}+\frac{1}{\phi}\left(\bar{\nabla}_{\mu}\bar{\nabla}_{\nu}-g_{\mu\nu}\bar{\Box}\right)\phi-\frac{3}{2\phi^{2}}\bar{\nabla}_{\mu}\phi\bar{\nabla}_{\nu}\phi+\frac{3}{2}\frac{\bar{\nabla}_{\mu}\psi\bar{\nabla}_{\nu}\psi}{\phi^{2\lambda}+(1-\lambda)^{2}\psi^{2}}$
$\displaystyle+\frac{1}{2}g_{\mu\nu}\left(\frac{3(\bar{\nabla}\phi)^{2}}{2\phi^{2}}-\frac{3}{2}\frac{(\bar{\nabla}\psi)^{2}}{\phi^{2\lambda}+(1-\lambda)^{2}\psi^{2}}-\frac{V(\phi,\psi)}{\phi}\right),$
(22) $\displaystyle 2V(\phi,\psi)-\phi\frac{\partial
V(\phi,\psi)}{\partial\phi}+\frac{3\lambda\phi^{2\lambda+1}}{\left(\phi^{2\lambda}+(1-\lambda)^{2}\psi^{2}\right)^{2}}(\bar{\nabla}\psi)^{2}=\kappa
T,$ (23)
$\displaystyle\bar{\Box}\psi-\frac{(1-\lambda)^{2}\psi}{\phi^{2\lambda}+(1-\lambda)^{2}\psi^{2}}(\bar{\nabla}\psi)^{2}+\left(1-\frac{2\lambda\phi^{2\lambda}}{\phi^{2\lambda}+(1-\lambda)^{2}\psi^{2}}\right)\bar{\nabla}_{\mu}\ln\phi\bar{\nabla}^{\mu}\psi=\frac{\partial
V(\phi,\psi)}{3\partial\psi}.$ (24)

As will be shown, cosmological solutions can be found for projective invariant
models ($\lambda=1$), and restricting to potentials of the form
$V(\phi,\psi)=V(\phi)$. To this end, we consider the metric for a Bianchi I
flat spacetime, i.e.

$ds^{2}=-dt^{2}+a(t)^{2}dx^{2}+b(t)^{2}dy^{2}+c(t)^{2}dz^{2},$ (25)

which represents a homogeneous but anisotropic spacetime, with three different
scale factors $a(t),b(t),c(t)$. We include matter in the form of a perfect
fluid with stress energy tensor

$T_{\mu\nu}=\text{diag}(\rho,a^{2}p,b^{2}p,c^{2}p),$ (26)

where $\rho$ is the energy density and $p$ the pressure. Assuming the equation
of state $p=w\rho$, conservation of the stress energy tensor implies

$\rho(t)=\frac{\mu^{2}}{(abc)^{w+1}},$ (27)

where $\mu^{2}$ is a constant. Now, considering a Starobinsky quadratic
potential [70]

$V(\phi)=\frac{1}{\alpha}\left(\frac{\phi-1}{2}\right)^{2},$ (28)

it can be shown that the equations for the scalar fields yield111A dot denotes
derivatives with respect to time.

$\dot{\psi}=\frac{k_{0}\phi}{v},\qquad\qquad\qquad\phi=\frac{v^{2}f(v)}{6\alpha
k_{0}^{2}+v^{2}},$ (29)

in terms of the volume-like variable $v\equiv abc$ and the function
$f(v)\equiv 1-2\alpha\kappa(3w-1)\rho(v)$, while $k_{0}$ is an integration
constant. Regarding the metric field equations, lengthy computations[65] allow
to rewrite the tt component in the form of a modified Friedmann equation:

$H^{2}\equiv\left(\frac{\dot{v}}{3v}\right)^{2}=\frac{\frac{\kappa}{3}\left(\frac{\mu_{I}^{2}}{v^{2}}+\frac{\rho}{\phi}+\frac{\mu_{AN}^{2}}{\phi^{2}v^{2}}\right)+\frac{V(\phi)}{6\phi}}{\left(1+\frac{3v}{2}\frac{d}{dv}\ln\phi\right)^{2}},$
(30)

where $\mu_{A}^{2}$ and $\mu^{2}_{I}$ are constants, representing the
anisotropy density parameter and the energy density parameter of the Immirzi
field, respectively. We note that the r.h.s is a rational function of the
volume alone. In the following we will take into account dust and radiation as
matter contributions, specified by the choices $w=0$, $\mu=\mu_{D}$ and
$w=1/3$, $\mu=\mu_{R}$ in (27), respectively. Then, bouncing solutions can be
derived for $\alpha<0$ integrating Eq.(30) for $v(t)$, which then yields the
scalar fields behavior via (29).

###  4.1 Vacuum case

It is convenient to first focus on the vacuum case, where $f(v)=1$ and
$\mu_{D}=\mu_{R}=0$. In this case we found the big bounce in the volume
variable depicted in Fig. 1, where also the behavior of each scale factor is
shown.

(a)

(b)

Figure 1: Numerical solutions for $\alpha=-5/3$, $\mu_{I}=\sqrt{3}$,
$\mu_{A}=0.2\mu_{I}$ as a function of $t/t_{Pl}$. Dotted and dashed lines
represent where bounce and future time singularity happen, respectively. The
bounce is centered at the origin of time for convenience, and the values of
the parameters are chosen in order to yield graphs that display features in a
clear fashion. (a) Universe volume normalised to $v_{B}$. (b) Scale factors.

We note that the volume is affected by a future finite-time singularity [71,
72] where the Hubble function diverges, while the scale factors are always
finite and nonvanishing. Such singular points will be carefully investigated
in the next section studying the geodesic completeness and scalar
perturbations across them. Here we just note that, in general, quantum effects
of particle creation [73, 74, 75, 76] can produce additional effective terms
in the Friedman equation, able to regularize singularities of the Hubble
function. Finally, the scalar fields behaviour is shown in Fig. 2 where the
field $\phi$ asymptotically reaches unity as $t\to\infty$, while the Immirzi
field relaxes to a constant Immirzi parameter.

Figure 2: Scalaron $\phi$ and Immirzi field derivative $\dot{\psi}$ for
$\alpha=-5/3$, $\mu_{I}=\sqrt{3}$, $\mu_{A}=0.2\mu_{I}$ as a function of
$t/t_{Pl}$.

###  4.2 Radiation and dust

The above analysis can be extended to the non vacuum case, including the
energy density of radiation and dust. A first difference concerns the late
time region, since both radiation and dust are able to provide an
isotropization mechanism for the universe, as is apparent studying the
anisotropy degree

$A(t)=\frac{\left(H_{A}^{2}+H_{B}^{2}+H_{C}^{2}\right)}{3H^{2}}-1,$ (31)

whose behavior is shown in Fig. 3.

Figure 3: Asymptotic behavior of the anisotropy degree $A$ as a function of
$t/t_{Pl}$ after the finite time singularity for various values of $\mu_{R}$,
$\mu_{D}$ and $\alpha=-5/3$, $\mu_{I}=\sqrt{3}$, $\mu_{A}=0.2\mu_{I}$.

Regarding the universe in its early phase, instead, two different scenarios
may occur, depending on the value of the parameter $\alpha$. If
$\bar{\alpha}<\alpha<0$, where $\bar{\alpha}=-2\mu_{I}^{2}/\mu_{D}^{4}$, the
solutions are qualitatively equivalent to the vacuum case. If
$\bar{\alpha}>\alpha$, instead, we still obtain a bouncing behavior for the
volume (See Fig. 4), which however is now devoid of singularities. However,
now the scale factors either diverge or vanish at some critical time $t_{c}$.

(a)

(b)

Figure 4: Numerical solutions as a function of $t/t_{Pl}$ for $\mu_{I}=0.057$,
$\mu_{A}=2.4$, $\mu_{D}=0.365$, $\mu_{R}=1.56$ and
$\alpha=-8.42<\bar{\alpha}$. The dashed lines represent where the scalaron
vanishes. (a) Volume normalised to $v_{B}$, $\phi$ and derivative of the
Immirzi field. (b) Scale factors.

##  5 Physical implications of curvature divergences

As described in the previous section there are two classes of solutions,
characterized by a singularity in the Hubble function and non vanishing and
finite scale factors (Fig. 1) or with regular Hubble function but vanishing or
divergent scale factors (Fig. 4). In this section we will analyse such
singularities studying the behavior of null geodesics and of scalar
perturbations near the critical time $t_{c}$, at which the singularity is
located.

Regarding null geodesics with tangent vector $u^{\alpha}=dx^{\alpha}/ds$, it
can be shown that they admit a first integral of the form [77, 78]

$x^{\prime}=\frac{k_{a}}{a^{2}},\qquad y^{\prime}=\frac{k_{b}}{b^{2}},\qquad
z^{\prime}=\frac{k_{c}}{c^{2}},\qquad
t^{\prime}=\left(\frac{k^{2}_{a}}{a^{2}}+\frac{k^{2}_{b}}{b^{2}}+\frac{k^{2}_{c}}{c^{2}}\right)^{1/2}+C_{0},$
(32)

where prime denotes derivative with respect to the affine parameter $s$ and
${k_{a},k_{b},k_{c}}$ and $C_{0}$ are integration constants. It follows that
if $a(t)$, $b(t)$, and $c(t)$ are continuous and non-vanishing, as in Fig. 1,
the tangent vector to the geodesics will be unique and well defined.
Therefore, such cases are geodesically complete, a result that holds both in
the anisotropic and in the isotropic case [79].

In the other class of solutions (Fig. 4), instead, we see that the volume
remains finite despite the vanishing/divergence of some scale factors. The
divergence of individual scale factors does not affect the geodesics, but the
vanishing of some of them may spoil the continuity and lead to the
impossibility of a unique extension across $t_{c}$. Indeed, let us consider
the case in which one scale factor vanishes at some affine parameter $s_{c}$.
In particular, suppose that

$a(s)=a_{0}(s-s_{c})^{\gamma},$ (33)

with $\gamma>0$. Then, integrating the relevant equations yields

$\displaystyle x(s)$
$\displaystyle=x_{c}+\frac{k_{a}(s-s_{c})^{1-2\gamma}}{a_{0}^{2}(1-2\gamma)},$
(34) $\displaystyle t(s)$
$\displaystyle=t_{c}+C_{0}(s-s_{c})+\frac{k_{a}(s-s_{c})^{1-\gamma}}{a_{0}(1-\gamma)},$
(35)

which are smooth if $0<\gamma<1/2$ and $0<\gamma<1$, respectively. Note that
if $1/2<\gamma<1$, then $x(s)\xrightarrow[]{s\to s_{c}}\pm\infty$, which would
imply reaching infinity in finite coordinate time. Conversely, if
$0<\gamma<1/2$, then the geodesic path will span the range
$\\{t,x\\}\in\;(-\infty,\infty)$, and the geodesics would be complete, despite
the vanishing of some scale factors. However, this result is dependent on how
rapidly the zero is reached, i.e. we have to determine the value of $\gamma$
relative to the solutions in question. The results are shown in Fig.5 and
prove that the solution approaches zero too rapidly, corresponding to a value
of $\gamma$ larger than $1/2$. We are thus forced to conclude that the example
shown in Fig. 4 does represent a geodesically incomplete space-time.

Figure 5: Outcomes of null geodesics test for $\alpha<\bar{\alpha}$. Scale
factor $a(s(t))$ for different values of $\gamma$. The dashed-black line
represent the numerical solution $a(t)$ reported in Fig. 4.

We now turn the attention to the behavior of scalar field perturbations. For a
scalar mode given by
$\sigma_{\vec{k}}(t,\vec{x})=\Theta(t)e^{i\vec{k}\cdot\vec{x}}$, one finds[65]
an equation of the form

$\ddot{\Theta}+h(v)\frac{\dot{v}}{v}\dot{\Theta}+\left(\frac{k_{x}^{2}}{a^{2}}+\frac{k_{y}^{2}}{b^{2}}+\frac{k_{z}^{2}}{c^{2}}\right)\Theta=0\
,$ (36)

where $h(v)$ represents some regular function of the volume $v$ and
$\vec{k}=(k_{x},k_{y},k_{z})$ represents a set of constants. From this
expression, we see that scalar modes are sensible to the presence of the
individual scale factors $a,b,$ and $c$, and of the Hubble function
$H=\dot{v}/3v$. If the scalar factors do not vanish anywhere, then any
potential problems should come only from the term involving the Hubble
function, which diverges in $t_{c}$. In particular, in vacuum one finds[65]
the following approximate solution near the singularity

$\Theta(t)\approx\Theta_{c}+\frac{\dot{\Theta}_{c}}{2\tilde{h}_{c}^{2}}e^{\mp
2\tilde{h}_{c}|t-t_{c}|^{1/2}}\left(1\pm 2\tilde{h}_{c}|t-t_{c}|^{1/2}\right)\
,$ (37)

with $\Theta_{c},\dot{\Theta}_{c},\tilde{h}_{c}$ constants, from which we see
that scalar field perturbations remain bounded around $t_{c}$ despite the
divergence in the Hubble function. A similar result holds in presence of dust
and radiation, since

$\displaystyle\Theta\approx$
$\displaystyle\Theta_{c}-\frac{\dot{\Theta}_{c}}{32\tilde{h}_{c}^{4}}e^{\mp
4\tilde{h}_{c}|t-t_{c}|^{1/4}}\left(3\pm
12\tilde{h}_{c}|t-t_{c}|^{1/4}\right.$
$\displaystyle\left.+24\tilde{h}_{c}^{2}|t-t_{c}|^{1/2}\pm
32\tilde{h}_{c}^{3}|t-t_{c}|^{3/4}\right),$ (38)

which is easy to see to be again bounded. On the other hand, for a regular
Hubble function but vanishing scale factors, equation (36) describes a
harmonic oscillator with a time dependent frequency,

$\ddot{\Theta}(t)+\frac{k_{x}^{2}}{a^{2}(t)}\Theta(t)\approx 0,$ (39)

which diverges as $a(t)\to 0$. Therefore, for the values of $a(t)$ obtained
numerically in the previous section neither geodesics nor scalar perturbations
are well behaved.

##  6 Conclusion

We proposed a generalization of the Nieh-Yan term to metric-affine gravity, by
including an additional term featuring nonmetricity and inserting two
parameters ($\lambda_{1}$, $\lambda_{2}$), which allow to recover the
projective invariance and the topological character, otherwise lost in
presence of nonmetricity. In particular, projective invariance can be
independently obtained by setting $\lambda_{1}=\lambda_{2}=\lambda$, whereas
topologicity is only guaranteed for $\lambda=1$. As an explicit example, we
considered a model with Lagrangian $F(R,NY_{gen})$, which re-expressed in the
Jordan frame features two scalar fields. We identified these scalar degrees as
the $f(R)$-like scalaron $\phi$ and the Immirzi field $\beta$ and showed that
the latter acquires dynamical character and a potential term in a more natural
way than in previous treatments, where these features are introduced by hand
in the action. Depending on the values of $\lambda_{1}$ and $\lambda_{2}$, we
found two different scalar tensor theories. Models with
$\lambda_{1}\neq\lambda_{2}$, are non-dynamical, and the scalar fields are
frozen to constant values in vacuum. In the projective invariant case
($\lambda_{1}=\lambda_{2}$), instead, the theory is endowed with one
additional dynamical degree of freedom, the Immirzi field, while the field
$\phi$ is algebraically related to the latter via a modified structural
equation. Then, we considered in more detail the dynamical models, looking for
cosmological solutions in Bianchi I spacetimes. We found numerical solutions
in which the big bang singularity is replaced by a big bounce scenario, in
which the universe volume undergoes a contraction up to a minimum value and
then bounces back, re-expanding in another branch. The scalaron and the
Immirzi field reach a maximum during the bounce and relax to constant values
at later times, where the standard LQG picture, with $\phi=1$ and a constant
Immirzi parameter $\beta=\beta_{0}$, is recovered. Moreover, the inclusion of
dust and radiation turns out to provide an isotropization mechanism at late
times, a feature that is absent in the vacuum case.

Such solutions are characterized by future finite time singularities after the
bounce, either in the Hubble function or in the individual scale factors. In
the former case, we showed that null geodesics are still well behaved and
scalar perturbations bounded, which allows us to conclude that the solution is
physically acceptable. In the latter case, instead, the study of null
geodesics shows that they cannot be extended across the singular point, where
also scalar perturbations grow in time, leading us to regard such solutions as
unphysical.

## References

  * [1] Wald R. M., General Relativity (University of Chicago Press, 1984).
  * [2] T. K. S. W. J. A. Misner C. W. and K. D.I., Gravitation (Princeton University Press, 2017).
  * [3] V. I. Afonso, G. J. Olmo and D. Rubiera-Garcia, Mapping Ricci-based theories of gravity into general relativity, Phys. Rev. D 97, p. 021503 (2018).
  * [4] V. I. Afonso, G. J. Olmo, E. Orazi and D. Rubiera-Garcia, Correspondence between modified gravity and general relativity with scalar fields, Phys. Rev. D 99, p. 044040 (2019).
  * [5] G. J. Olmo, PALATINI APPROACH TO MODIFIED GRAVITY: f(R) THEORIES AND BEYOND, International Journal of Modern Physics D 20, 413 (2011).
  * [6] F. S. N. Lobo, G. J. Olmo and D. Rubiera-Garcia, Semiclassical geons as solitonic black hole remnants, JCAP 07, p. 011 (2013).
  * [7] J. Beltran Jimenez, L. Heisenberg, G. J. Olmo and D. Rubiera-Garcia, Born–Infeld inspired modifications of gravity, Phys. Rept. 727, 1 (2018).
  * [8] J. Beltrán Jiménez, L. Heisenberg, D. Iosifidis, A. Jiménez-Cano and T. S. Koivisto, General teleparallel quadratic gravity, Physics Letters B 805, p. 135422 (2020).
  * [9] J. a. L. Rosa, F. S. N. Lobo and G. J. Olmo, Weak-field regime of the generalized hybrid metric-Palatini gravity (4 2021).
  * [10] T. Harko, T. S. Koivisto, F. S. N. Lobo and G. J. Olmo, Metric-Palatini gravity unifying local constraints and late-time cosmic acceleration, Phys. Rev. D 85, p. 084016 (2012).
  * [11] S. Capozziello, T. Harko, T. S. Koivisto, F. S. N. Lobo and G. J. Olmo, Hybrid metric-Palatini gravity, Universe 1, 199 (2015).
  * [12] F. Bombacigno, F. Moretti and G. Montani, Scalar modes in extended hybrid metric-Palatini gravity: weak field phenomenology, Phys. Rev. D 100, p. 124036 (2019).
  * [13] K. A. Bronnikov, Spherically symmetric black holes and wormholes in hybrid metric-Palatini gravity, Grav. Cosmol. 25, 331 (2019).
  * [14] M. Borunda, B. Janssen and M. Bastero-Gil, Palatini versus metric formulation in higher curvature gravity, JCAP 11, p. 008 (2008).
  * [15] B. Janssen, A. Jiménez-Cano and J. A. Orejuela, A non-trivial connection for the metric-affine Gauss-Bonnet theory in $D=4$, Phys. Lett. B 795, 42 (2019).
  * [16] B. Janssen and A. Jiménez-Cano, On the topological character of metric-affine Lovelock Lagrangians in critical dimensions, Phys. Lett. B 798, p. 134996 (2019).
  * [17] R. Percacci and E. Sezgin, New class of ghost- and tachyon-free metric affine gravities, Phys. Rev. D 101, p. 084040 (2020).
  * [18] C. Rovelli, Quantum Gravity (Cambridge University Press, Cambridge, England, 2004).
  * [19] T. Thiemann, Modern Canonical Quantum General Relativity (Cambridge University Press, Cambridge, England, 2007).
  * [20] A. Ashtekar, New Variables for Classical and Quantum Gravity, Phys. Rev. Lett. 57, 2244 (nov 1986).
  * [21] A. Ashtekar, New Hamiltonian formulation of general relativity, Physical Review D 36, 1587 (sep 1987).
  * [22] A. Ashtekar, J. D. Romano and R. S. Tate, New variables for gravity: Inclusion of matter, Physical Review D 40, 2572 (oct 1989).
  * [23] A. Ashtekar and C. J. Isham, Representations of the holonomy algebras of gravity and nonAbelian gauge theories, Classical and Quantum Gravity 9, 1433 (jun 1992).
  * [24] G. Date, R. K. Kaul and S. Sengupta, Topological interpretation of Barbero-Immirzi parameter, Phys. Rev. D 79, p. 44008 (feb 2009).
  * [25] H. T. Nieh and M. L. Yan, An identity in Riemann–Cartan geometry, Journal of Mathematical Physics 23, 373 (1982).
  * [26] H. T. Nieh, A TORSIONAL TOPOLOGICAL INVARIANT, International Journal of Modern Physics A 22, 5237 (2007).
  * [27] S. Holst, Barbero’s Hamiltonian derived from a generalized Hilbert-Palatini action, Phys. Rev. D 53, 5966 (may 1996).
  * [28] G. Immirzi, Real and complex connections for canonical gravity, Class. Quant. Grav. 14, L177 (1997).
  * [29] G. Immirzi, Real and complex connections for canonical gravity, Classical and Quantum Gravity 14, L177 (oct 1997).
  * [30] C. Rovelli and T. Thiemann, The Immirzi parameter in quantum general relativity, Phys. Rev. D 57, 1009 (1998).
  * [31] V. Taveras and N. Yunes, Barbero-Immirzi parameter as a scalar field: K-inflation from loop quantum gravity?, Phys. Rev. D 78, p. 64070 (sep 2008).
  * [32] G. Calcagni and S. Mercuri, Barbero-Immirzi field in canonical formalism of pure gravity, Phys. Rev. D 79, p. 84004 (apr 2009).
  * [33] F. Bombacigno, F. Cianfrani and G. Montani, Big-bounce cosmology in the presence of Immirzi field, Phys. Rev. D 94, p. 64021 (sep 2016).
  * [34] O. J. Veraguth and C. H. Wang, Immirzi parameter without Immirzi ambiguity: Conformal loop quantization of scalar-tensor gravity, Physical Review D 96, p. 084011 (oct 2017).
  * [35] C. H. Wang and D. P. Rodrigues, Closing the gaps in quantum space and time: Conformally augmented gauge structure of gravitation, Physical Review D 98, p. 124041 (dec 2018).
  * [36] F. Bombacigno and G. Montani, Implications of the Holst term in a f(R) theory with torsion, Phys. Rev. D 99, p. 64016 (mar 2019).
  * [37] F. Bombacigno and G. Montani, f(R) gravity with torsion and the Immirzi field: Signature for gravitational wave detection, Phys. Rev. D 97, p. 124066 (jun 2018).
  * [38] C. H.-T. Wang and M. Stankiewicz, Quantization of time and the big bang via scale-invariant loop gravity, Physics Letters B 800, p. 135106 (jan 2020).
  * [39] D. Iosifidis and L. Ravera, Parity Violating Metric-Affine Gravity Theories (9 2020).
  * [40] F. Bombacigno, S. Boudet and G. Montani, Generalized ashtekar variables for palatini f(r) models, Nuclear Physics B 963, p. 115281 (2021).
  * [41] M. Långvik, J.-M. Ojanperä, S. Raatikainen and S. Räsänen, Higgs inflation with the Holst and the Nieh–Yan term, Phys. Rev. D 103, p. 083514 (2021).
  * [42] F. Bombacigno and G. Montani, Big bounce cosmology for Palatini $R^{2}$ gravity with a Nieh-Yan term, The European Physical Journal C 79, p. 405 (may 2019).
  * [43] S. Boudet, F. Bombacigno, G. Montani and M. Rinaldi, Superentropic black hole with immirzi hair, Phys. Rev. D 103, p. 084034 (Apr 2021).
  * [44] M. Lattanzi and S. Mercuri, A solution of the strong $cp$ problem via the peccei-quinn mechanism through the nieh-yan modified gravity and cosmological implications, Phys. Rev. D 81, p. 125015 (Jun 2010).
  * [45] O. Castillo-Felisola, C. Corral, S. Kovalenko, I. Schmidt and V. E. Lyubovitskij, Axions in gravity with torsion, Phys. Rev. D 91, p. 085017 (Apr 2015).
  * [46] S. Mercuri, Peccei-Quinn Mechanism in Gravity and the Nature of the Barbero-Immirzi Parameter, Phys. Rev. Lett. 103, p. 81302 (aug 2009).
  * [47] M. Li, H. Rao and D. Zhao, A simple parity violating gravity model without ghost instability, Journal of Cosmology and Astroparticle Physics 2020, 023 (nov 2020).
  * [48] J. Nissinen and G. E. Volovik, Thermal nieh-yan anomaly in weyl superfluids, Phys. Rev. Research 2, p. 033269 (Aug 2020).
  * [49] J. Nissinen and G. E. Volovik, On Thermal Nieh-Yan Anomaly in Topological Weyl Materials, JETP Letters 110, 789 (2019).
  * [50] Z.-M. Huang, B. Han and M. Stone, Nieh-yan anomaly: Torsional landau levels, central charge, and anomalous thermal hall effect, Phys. Rev. B 101, p. 125201 (Mar 2020).
  * [51] C.-X. Liu, Phonon helicity and nieh-yan anomaly in the kramers-weyl semimetals of chiral crystals (2021).
  * [52] V. I. Afonso, C. Bejarano, J. Beltran Jimenez, G. J. Olmo and E. Orazi, The trivial role of torsion in projective invariant theories of gravity with non-minimally coupled matter fields, Class. Quant. Grav. 34, p. 235003 (2017).
  * [53] D. Iosifidis, Linear Transformations on Affine-Connections, Class. Quant. Grav. 37, p. 085010 (2020).
  * [54] J. Beltrán Jiménez and A. Delhom, Ghosts in metric-affine higher order curvature gravity, The European Physical Journal C 79, p. 656 (2019).
  * [55] G. Montani, M. V. Battisti, R. Benini and G. Imponente, Primordial cosmology (World Scientific, Singapore, 2009).
  * [56] A. Corichi and P. Singh, A Geometric perspective on singularity resolution and uniqueness in loop quantum cosmology, Phys. Rev. D 80, p. 044024 (2009).
  * [57] A. Y. Kamenshchik, E. O. Pozdeeva, A. A. Starobinsky, A. Tronconi, G. Venturi and S. Y. Vernov, Induced gravity, and minimally and conformally coupled scalar fields in Bianchi-I cosmological models, Phys. Rev. D 97, p. 023536 (2018).
  * [58] G. Montani and R. Chiovoloni, A Scenario for a Singularity-free Generic Cosmological Solution (3 2021).
  * [59] G. Montani, A. Marchi and R. Moriconi, Bianchi I model as a prototype for a cyclical Universe, Phys. Lett. B 777, 191 (2018).
  * [60] F. Cianfrani, A. Marchini and G. Montani, The picture of the Bianchi I model via gauge fixing in Loop Quantum Gravity, EPL 99, p. 10003 (2012).
  * [61] R. Moriconi and G. Montani, Behavior of the Universe anisotropy in a big-bounce cosmology, Phys. Rev. D 95, p. 123533 (2017).
  * [62] E. Giovannetti, G. Montani and S. Schiattarella, Semiclassical and quantum features of the Bianchi I cosmology in the polymer representation (5 2021).
  * [63] C. Barragan and G. J. Olmo, Isotropic and Anisotropic Bouncing Cosmologies in Palatini Gravity, Phys. Rev. D 82, p. 084015 (2010).
  * [64] C. Barragan, G. J. Olmo and H. Sanchis-Alepuz, Bouncing Cosmologies in Palatini f(R) Gravity, Phys. Rev. D 80, p. 024016 (2009).
  * [65] F. Bombacigno, S. Boudet, G. J. Olmo and G. Montani, Big bounce and future time singularity resolution in bianchi i cosmologies: The projective invariant nieh-yan case, Phys. Rev. D 103, p. 124031 (Jun 2021).
  * [66] G. J. Olmo, Post-Newtonian constraints on f(R) cosmologies in metric and Palatini formalism, Phys. Rev. D 72, p. 083505 (2005).
  * [67] S. Mercuri, Fermions in the Ashtekar-Barbero connection formalism for arbitrary values of the Immirzi parameter, Phys. Rev. D 73, p. 84016 (apr 2006).
  * [68] S. Mercuri, From the Einstein-Cartan to the Ashtekar-Barbero canonical constraints, passing through the Nieh-Yan functional, Physical Review D \- Particles, Fields, Gravitation and Cosmology 77, p. 024036 (jan 2008).
  * [69] S. Mercuri and V. Taveras, Interaction of the Barbero-Immirzi field with matter and pseudoscalar perturbations, Phys. Rev. D 80, p. 104007 (nov 2009).
  * [70] A. Starobinsky, A new type of isotropic cosmological models without singularity, Physics Letters B 91, 99 (1980).
  * [71] S. Nojiri, S. D. Odintsov and S. Tsujikawa, Properties of singularities in (phantom) dark energy universe, Phys. Rev. D 71, p. 063004 (2005).
  * [72] S. D. Odintsov and V. K. Oikonomou, Dynamical Systems Perspective of Cosmological Finite-time Singularities in $f(R)$ Gravity and Interacting Multifluid Cosmology, Phys. Rev. D 98, p. 024013 (2018).
  * [73] G. Montani, Influence of the particles creation on the flat and negative curved FLRW universes, Class. Quant. Grav. 18, 193 (2001).
  * [74] K. Dimopoulos, Is the Big Rip unreachable?, Phys. Lett. B 785, 132 (2018).
  * [75] L. H. Ford, Gravitational Particle Creation and Inflation, Phys. Rev. D 35, p. 2955 (1987).
  * [76] F. Contreras, N. Cruz, E. Elizalde, E. González and S. Odintsov, Linking little rip cosmologies with regular early universes, Phys. Rev. D 98, p. 123520 (2018).
  * [77] P. Singh, Curvature invariants, geodesics and the strength of singularities in Bianchi-I loop quantum cosmology, Phys. Rev. D 85, p. 104011 (2012).
  * [78] K. Nomura and D. Yoshida, Past extendibility and initial singularity in Friedmann-Lemaître-Robertson-Walker and Bianchi I spacetimes (5 2021).
  * [79] J. Beltran Jimenez, R. Lazkoz, D. Saez-Gomez and V. Salzano, Observational constraints on cosmological future singularities, Eur. Phys. J. C 76, p. 631 (2016).


# Could nearby star-forming galaxies light up the point-like neutrino sky?

Antonio Ambrosone111ciao<EMAIL_ADDRESS>Dipartimento di Fisica “Ettore
Pancini”, Università degli studi di Napoli “Federico II”, Complesso Univ.
Monte S. Angelo, I-80126 Napoli, Italy  INFN - Sezione di Napoli, Complesso
Univ. Monte S. Angelo, I-80126 Napoli, Italy  Marco Chianese  Dipartimento di
Fisica “Ettore Pancini”, Università degli studi di Napoli “Federico II”,
Complesso Univ. Monte S. Angelo, I-80126 Napoli, Italy  INFN - Sezione di
Napoli, Complesso Univ. Monte S. Angelo, I-80126 Napoli, Italy  Damiano F.G.
Fiorillo  Dipartimento di Fisica “Ettore Pancini”, Università degli studi di
Napoli “Federico II”, Complesso Univ. Monte S. Angelo, I-80126 Napoli, Italy
INFN - Sezione di Napoli, Complesso Univ. Monte S. Angelo, I-80126 Napoli,
Italy  Antonio Marinelli  INFN - Sezione di Napoli, Complesso Univ. Monte S.
Angelo, I-80126 Napoli, Italy  Gennaro Miele  Dipartimento di Fisica “Ettore
Pancini”, Università degli studi di Napoli “Federico II”, Complesso Univ.
Monte S. Angelo, I-80126 Napoli, Italy  INFN - Sezione di Napoli, Complesso
Univ. Monte S. Angelo, I-80126 Napoli, Italy  Scuola Superiore Meridionale,
Università degli studi di Napoli “Federico II”, Largo San Marcellino 10, 80138
Napoli, Italy

###### Abstract

Star-forming and starburst galaxies, which are well-known cosmic-rays
reservoirs, are expected to emit gamma-rays and neutrinos predominantly via
hadronic collisions. In this Letter, we analyze the 10-year Fermi-LAT spectral
energy distributions of 13 nearby galaxies by means of a physical model which
accounts for high-energy proton transport in starburst nuclei and includes the
contribution of primary and secondary electrons. In particular, we test the
hypothesis that the observed gamma-ray fluxes are mostly due to star-forming
activity, in agreement with the available star formation rates coming from IR
and UV observations. Through this observation-based approach, we determine the
most-likely neutrino counterpart from star-forming and starburst galaxies and
quantitatively assess the ability of current and upcoming neutrino telescopes
to detect them as point-like sources. Remarkably, we find that the cores of
the Small Magellanic Cloud and the Circinus galaxy are potentially observable
by KM3NeT/ARCA with 6 years of observation. Moreover, most of the nearby
galaxies are likely to be just a factor of a few below the KM3NeT and IceCube-
Gen2 point-like sensitivities. After investigating the prospects for detection
of gamma-rays above TeV energies from these sources, we conclude that the
joint observations of high-energy neutrinos and gamma-rays with upcoming
telescopes will be an objective test for our emission model and may provide
compelling evidence of star-forming activity as a tracer of neutrino
production.

galaxies: starburst — stars: formation — gamma rays — neutrinos

##  1 Introduction

IceCube has been measuring high-energy neutrinos from tens of TeV up to PeVs
(Aartsen et al., 2020; Albert et al., 2020) during the last decade and yet
their origin is still unclear. These measurements certainly reveal the
presence of neutrino emitters in our Universe, even though the only clear
multi-messenger indication of a neutrino source has been related to the blazar
TXS 0506+056 Aartsen et al. (2018). In this context, star-forming and
starburst galaxies (SFGs and SBGs) are abundant sources which could be
responsible for at least a part of such observations (Ambrosone et al., 2021).
Indeed, SFGs and SBGs are galaxies endowed with an intense star-forming
activity, and therefore they have a particularly high gas density $(n_{\rm
ISM}\propto\ 10$–$100\ {\rm cm^{-3}})$, strong magnetic fields $(B\sim
10^{2}{\rm\mu G})$, and supernova explosion rates of the order of $0.01$–$1\
{\rm yr^{-1}}$ (Thompson et al., 2006). For these reasons, they are supposed
to exhibit strong turbulence which should be able to trap high-energy cosmic-
rays inside their environment, giving them an enhanced probability to produce
high-energy neutrinos and gamma-rays through hadronic collisions (Peretti et
al., 2019). In this regard, many authors have studied the properties of the
diffuse gamma-ray and neutrino fluxes of these sources (Loeb & Waxman, 2006;
Murase et al., 2013; Tamborra et al., 2014; Senno et al., 2015; Bechtol et
al., 2017; Murase & Waxman, 2016; Sudoh et al., 2018; Peretti et al., 2020;
Ajello et al., 2020; Ambrosone et al., 2021; Blanco & Linden, 2021; Peretti et
al., 2021). Different cosmic-ray transport scenarios have been used to
describe the observed gamma-ray emission from SFGs and SBGs, disentangling a
calorimetric behaviour from a diffusion-dominated regime (Peretti et al.,
2019; Krumholz et al., 2020; Ha et al., 2020; Müller et al., 2020; Shimono et
al., 2021; Xiang et al., 2021b; Owen et al., 2021; Werhahn et al., 2021a, b).
However, in order to explain a sizeable portion of measured high-energy
neutrino flux with SFGs and SBGs, the deep Universe must be considered up to
redshift $\sim$ 4–5, because of their dimness (Peretti et al., 2020; Ajello et
al., 2020; Ambrosone et al., 2021). The low gamma-ray luminosity of SFGs and
SBGs typically represents a bound for their contribution to the observed
astrophysical neutrinos as a point-like component. Currently, only a dozen of
these sources have been catalogued as gamma-ray point-like sources using the
Fermi-LAT data and only few of them have been observed through Imaging
Cherenkov telescopes.

In this Letter, for the first time we employ a multi-messenger and multi-
wavelength approach to assess the ability of current and upcoming neutrino
telescopes to observe such galaxies as neutrino point-like sources. This is a
question of utmost importance after the IceCube collaboration reported a
$2.9\sigma$ excess of neutrino events coming from the direction of NGC 1068
analyzing $10$-year data (Aartsen et al., 2020). Even though additional Active
Galactic Nuclei (AGN) activity could be present, observations like this one
suggest that star-forming activity could trace the hadronic emission and
produce point-like excesses in the TeV sky observed by IceCube as well as in
the future skymap of the upcoming KM3NeT/ARCA and IceCube-Gen2 telescopes.

In this work, following the scenario proposed by Peretti et al. (2019) which
considers the transport of high-energy protons and electrons, we analyze the
10-year Fermi-LAT data provided by Ajello et al. (2020), taking into account
the spectral energy distributions (SEDs) of 13 SFGs and SBGs. We describe
these gamma-ray observations through the hadronic and leptonic processes at
work in the star-forming regions and determine the corresponding neutrino
fluxes. For each SBG, we require the star formation rate (SFR) to be
consistent with the one derived from infra-red (IR) and ultra-violet (UV)
observations (Kornecki et al., 2020) within a maximal level of discrepancy,
thus making our predictions more robust. Hence, we compare the most-likely
neutrino flux normalizations at $1\ \text{TeV}$ with IceCube, IceCube-Gen2,
and KM3NeT sensitivities (see Figure 1). Besides, for the brightest sources we
also show in Figure 2 the prospects of the Cherenkov Telescope Array (CTA).

##  2 Modelling the neutrino and gamma-ray emission

In this section, we summarize the emission model adopted for the neutrino and
gamma-ray spectra of SFGs and SBGs. We adapt the model provided by Peretti et
al. (2019) (see also Peretti et al. (2020); Ambrosone et al. (2021) for
further details). We assume that the star-forming activity is limited to a
small spherical region of the galaxy called starburst nucleus (SBN). Hence, we
solve the leaky-box equations to determine the energy distribution of high-
energy cosmic-rays (protons and electrons) injected by supernovae explosions.
The solution of the leaky-box equations can be written as

$F_{p,e}=Q_{p,e}\left(\frac{1}{T_{\text{adv}}}+\frac{1}{T_{\text{loss}}}+\frac{1}{T_{\text{diff}}}\right)^{-1}\,,$
(1)

where $F_{p,e}$ and $Q_{p,e}$ are respectively the distribution function and
the injection rate for protons and electrons. The cosmic-ray distribution
stems from a balance of different processes with characteristic timescales:
the advection timescale $T_{\text{adv}}=R_{\mathrm{SBN}}/v_{\text{wind}}$,
with $R_{\mathrm{SBN}}$ being the SBN radius and $v_{\text{wind}}$ the dust-
wind velocity, the energy losses timescale $T_{\text{loss}}$, and the
diffusion timescale $T_{\text{diff}}$. For protons, energy losses are
dominated by ionization, proton-proton, and Coulomb interactions, whereas for
electrons they are dominated by ionization, Inverse Compton scatterings, and
synchrotron interactions (Peretti et al., 2019). For the cosmic-ray diffusion,
we consider a Kolmogorov-like model with diffusion coefficient $D(E)\propto
E^{1/3}$, so that $T_{\text{diff}}(E)\propto E^{-1/3}$. The diffusion
timescale is typically larger than the other ones due to the large magnetic
fields and the high levels of turbulence of SBGs and SFGs. This implies that
the gamma-ray and neutrino emission is only marginally affected by the
specific choice of the diffusion model.

We assume that protons are injected with a power-law spectrum in momentum
space with spectral index $\Gamma+2$ and an exponential cutoff at
$10~{}\mathrm{PeV}$, in agreement with the combined fit of IceCube and Fermi-
LAT diffuse data (Ambrosone et al., 2021). The precise value of the cutoff
energy is irrelevant for our purposes, since the gamma-ray flux is suppressed
already at much smaller energies by internal and external absorption. The
proton spectrum is directly proportional to the star formation rate
$\dot{M}_{*}$ and normalized by requiring that each supernova releases into
protons $10\%$ of its total explosion kinetic energy ($\sim
10^{51}~{}\mathrm{erg}$). On the other hand, electrons have the same spectral
index, a Gaussian cutoff at $10~{}\mathrm{TeV}$, and a normalization fixed to
one fiftieth of the one of the protons, similar to what is inferred for our
Galaxy (Torres, 2004; Peretti et al., 2019).

Neutrinos are emitted through the decay of charged pions
($\pi\to\mu\,\nu_{\mu}$, $\mu\to e\,\nu_{e}\,\nu_{\mu}$) that are produced in
hadronic interactions of the injected protons with the interstellar gas. We
determine the neutrino production rate by assuming that pions always carry 17%
of their parent proton energy (Kelner et al., 2006). Gamma-rays are emitted in
hadronic processes through neutral pion decays ($\pi^{0}\to\gamma\gamma$) and
in leptonic processes through bremsstrahlung and Inverse Compton scatterings
of electrons. We also take into account the contribution of secondary
electrons produced by the decay of charged pions. We emphasize that both the
hadronic and leptonic components examined in this study are dictated by the
star formation rate. Investigating the possible presence of additional
components with a different origin (e.g. related to AGN activities) is beyond
the scope of this work. The gamma-ray emission from Inverse Compton
scatterings depends on the density of background photons acting as a target.
In Peretti et al. (2019), such a radiation is modeled as a superposition of
black-body spectra. However, as emphasized also there, we can approximate it
as a monochromatic spectrum with an energy equal to the peak energy
$\epsilon_{\rm peak}=0.1\ \text{eV}$. Therefore, the Inverse Compton spectrum
can be simply characterized by the energy density $U_{\mathrm{rad}}$ of the
background photons. Analogously to Ambrosone et al. (2021), we consider all
sources to have a similar spectral shape for the background photons equal to
the M82 best-fit background spectrum reported by Peretti et al. (2019). The
normalization for the different sources is self-consistently determined by the
radiation energy density of each source $U_{\text{rad}}$ (see below). Finally,
we account for internal and external gamma-ray absorption due to pair-
production processes with background photons. For the latter, we consider as a
target CMB photons as well as the extragalactic background light model
reported in Franceschini & Rodighiero (2017). We point out that gamma-ray
absorption is relevant at energies larger than a few TeV, and therefore it
affects only the forecast predictions for the CTA telescope.

The neutrino and gamma-ray emission mainly depends on the spectral index
$\Gamma$ and on the star formation rate. For the density $n_{\mathrm{ISM}}$ of
the interstellar gas (target of the proton interactions), we rely on the
Kennicutt relation (Kennicutt, 1998; Kennicutt & Evans, 2012; Kennicutt & De
Los Reyes, 2021). This relation connects the surface density of SFR,
$\Sigma_{\rm SFR}=\dot{M}_{*}/(\pi R_{\mathrm{SBN}}^{2})$, and the gas surface
density, $\Sigma_{\text{gas}}=m_{p}\,n_{\mathrm{ISM}}\,R_{\mathrm{SBN}}$. For
simplicity, we consider for all the galaxies
$R_{\mathrm{SBN}}=200~{}\mathrm{pc}$, which is an average size of the
circumnuclear regions (Peretti et al., 2019). By using the benchmark value
$n_{\mathrm{ISM}}=175~{}{\rm cm^{-3}}$ for $\dot{M}_{*}=5~{}{\rm
M_{\odot}\,yr^{-1}}$ (Peretti et al., 2019), we obtain the scaling relation

$n_{\rm ISM}=175\left(\frac{\dot{M}_{*}}{5~{}{\rm
M_{\odot}\,yr^{-1}}}\right)^{2/3}\,{\rm cm^{-3}}\,.$ (2)

For the background photon density $U_{\rm rad}$, we assume a direct
proportionality to the SFR, which is expected to be tightly related to the
infrared (IR) luminosity (Kennicutt, 1998; Inoue et al., 2000; Hirashita et
al., 2003; Yuan et al., 2011; Kennicutt & Evans, 2012; Kennicutt & De Los
Reyes, 2021). In particular, we consider the scaling relation

$U_{\mathrm{rad}}=2500\left(\frac{\dot{M}_{*}}{5\,{\rm
M_{\odot}\,yr^{-1}}}\right)\,{\rm eV\,cm^{-3}}\,,$ (3)

where the reference values, $U_{\mathrm{rad}}=2500\,{\rm eV\,cm^{-3}}$ and
$\dot{M}_{*}=5\,{\rm M_{\odot}\,yr^{-1}}$, are obtained from what has been
inferred for M82 by Peretti et al. (2019, 2020). Such a relation is generally
satisfied by high-SFR galaxies which behave as good calorimeters. On the other
hand, low-SFR galaxies have been shown to deviate from the assumed linear SFR-
IR relation, as also recently pointed out by Kornecki et al. (2020); Werhahn
et al. (2021b). Nevertheless, for these galaxies we have verified that
equation 3 is a good and conservative approximation: indeed the gamma-ray
emission depends only slightly on $U_{\mathrm{rad}}$ since the leptonic
component is in general subdominant in our model. Moreover, a smaller value
for $U_{\rm rad}$, as observed for low-SFR galaxies (Kornecki et al., 2020),
would further reduce the leptonic production in favor of the hadronic one.
Finally, for all the galaxies, we adopt the average values $B=200~{}\mu{\rm
G}$ (Thompson et al., 2006) and $v_{\rm wind}=500~{}{\rm km/s}$ (Peretti et
al., 2020) for the magnetic field and the dust-wind velocity, respectively.

Table 1: Results of the likelihood analysis of current gamma-ray data. The
columns report the source name, the SFR prior, the most-likely values of the
two parameters, the 68% maximum posterior density credible intervals of the
marginal distributions, and the reduced chi-squared values considered as an
estimate of the goodness of the fit.

Source | Uniform prior | Most-likely values | 68% credible intervals | $\chi^{2}/\mathrm{dof}$
---|---|---|---|---
| $\dot{M}_{*}$ | $(\dot{M}_{*},\,\Gamma)$ | $\dot{M}_{*}$ | $\Gamma$ |
M82 | 3.0 – 30 | $\left(4.5,\,2.30\right)$ | $\left[4.3,\,4.6\right]$ | $\left[2.27,\,2.33\right]$ | 1.24
NGC 253 | 1.4 – 17 | $\left(3.3,\,2.30\right)$ | $\left[3.14,\,3.40\right]$ | $\left[2.28,\,2.32\right]$ | 1.32
ARP 220 | 60 – 740 | $\left(740,\,2.66\right)$ | $\left[492,\,740\right]$ | $\left[2.51,\,2.68\right]$ | 1.52
NGC 4945 | 0.35 – 4.15 | $\left(4.15,\,2.30\right)$ | $\left[4.05,\,4.15\right]$ | $\left[2.23,\,2.32\right]$ | 1.52
NGC 1068 | 5 – 93 | $\left(16,\,2.52\right)$ | $\left[13,\,20\right]$ | $\left[2.45,\,2.65\right]$ | 0.65
NGC 2146 | 3 – 57 | $\left(15,\,2.50\right)$ | $\left[9,\,27\right]$ | $\left[2.44,\,2.88\right]$ | 0.50
ARP 299 | 28 – 333 | $\left(28,\,2.15\right)$ | $\left[28,\,200\right]$ | $\left[1.40,\,1.90\right]\cup\left[2.77,\,3.00\right]$ | 0.18
M31 | 0.09 – 0.90 | $\left(0.34,\,2.40\right)$ | $\left[0.31,\,0.40\right]$ | $\left[2.29,\,2.61\right]$ | 0.52
M33 | 0.09 – 0.90 | $\left(0.44,\,2.76\right)$ | $\left[0.19,\,0.56\right]$ | $\left[2.57,\,2.96\right]$ | 0.44
NGC 3424 | 0.4 – 5.4 | $\left(5.4,\,2.22\right)$ | $\left[2.5,\,5.4\right]$ | $\left[1.92,\,2.67\right]$ | 1.63
NGC 2403 | 0.1 – 1.2 | $\left(0.75,\,2.12\right)$ | $\left[0.58,\,0.96\right]$ | $\left[1.92,\,2.36\right]$ | 0.38
SMC | 0.008 – 0.090 | $\left(0.038,\,2.14\right)$ | $\left[0.037,\,0.039\right]$ | $\left[2.13,\,2.16\right]$ | 1.90
Circinus Galaxy | 0.1 – 8.1 | $\left(6.6,\,2.32\right)$ | $\left[6.2,\,7.8\right]$ | $\left[2.15,\,2.45\right]$ | 0.92

Note. — The star formation rate $\dot{M}_{*}$ is in units of
$\text{M}_{\odot}\,\text{yr}^{-1}$.

##  3 Analysis and Results

Our main purpose is to determine the high-energy neutrino emission from SFGs
and SBGs based on a likelihood analysis of gamma-ray data. A crucial
observation must be made: for some of these galaxies, a possible additional
source of gamma-rays may be related to AGN activity (Alonso-Herrero et al.,
2013; Yoast-Hull et al., 2017; Guo et al., 2019; Inoue et al., 2020; Murase et
al., 2020; Ajello et al., 2020; Xiang et al., 2021a; Kheirandish et al.,
2021), which is not included in the following analysis. Thus we test the
hypothesis that the observed gamma-ray SEDs are saturated by star-forming
activity only, and determine the most-likely neutrino emission under this
assumption. We analyze the gamma-ray SEDs of 13 galaxies observed by Fermi-LAT
with 10 years of observation (Ajello et al., 2020).222The data reported in
Ajello et al. (2020) are normalized to a constant value of
$10^{-9}\,\text{GeV}\,\text{cm}^{-2}\,\text{s}^{-1}$ at an energy of $1$ GeV.
We have suitably rescaled them using the reported best-fit parameters of the
power-law model. For M82 and NGC 253 we make also use of the data provided by
VERITAS (VERITAS Collaboration et al., 2009) and H.E.S.S. (H. E. S. S.
Collaboration et al., 2018), respectively. For each galaxy, we pursue a
Bayesian approach to assess the most-likely values for the two free parameters
of the model: the star formation rate $\dot{M}_{*}$ and the spectral index
$\Gamma$ of injected protons and electrons. We determine the posterior
distribution as

$p(\dot{M}_{*},\,\Gamma|{\rm SED})\propto p({\rm
SED}|\dot{M}_{*},\,\Gamma)\,p(\dot{M}_{*})\,p(\Gamma)\,,$ (4)

with a Gaussian likelihood function

$p({\rm
SED}|\dot{M}_{*},\,\Gamma)=e^{-\frac{1}{2}\sum_{i}\left(\frac{\text{SED}_{i}-E^{2}_{i}\Phi_{\gamma}(E_{i}|\dot{M}_{*},\,\Gamma)}{\sigma_{i}}\right)^{2}}\,.$
(5)

Here, $\text{SED}_{i}$ are the measured data, where $i$ runs over the energy
bins centered around the energy $E_{i}$, and $\sigma_{i}$ are the
observational uncertainties. The data are compared to the gamma-ray emission
$\Phi_{\gamma}(E_{i}|\dot{M}_{*},\,\Gamma)$ predicted by our model. We adopt
the source distances provided by Kornecki et al. (2020) and neglect the
uncertainty on such values. For all the galaxies we consider the same uniform
prior on the spectral index in the range 1.0–3.0. For the SFR, we instead
adopt a different uniform prior distribution for each source. We require
$\dot{M}_{*}$ to be consistent within a factor of 3 with the corresponding
values reported in Kornecki et al. (2020). This choice is made to cover the
wide variety of SFR estimates presented in the literature (Groves et al.,
2008; Bolatto et al., 2011; For et al., 2012; Rahmani et al., 2016; Yoast-Hull
et al., 2017; Peretti et al., 2020; Ajello et al., 2020).

The results of the analysis are summarized in Table 1. For M82, NGC 253, NGC
2146, NGC 2403, Circinus galaxy, M31, M33, and SMC, we find a reasonable
agreement both with the data and the measured SFR values. We emphasize that in
our model all the galaxies are approximated as point-like sources with a
dominant emission limited to the SBN. Hence, the SBN extension intrinsically
defines the extension of the gamma-ray source emission. Most of the sources
have not been resolved and therefore they are consistent with our starburst
nucleus approximation. On the other hand, SMC and M31 have been observed by
Ajello et al. (2020) as extended sources with an angular extension of
$1.5^{\circ}$ and $0.45^{\circ}$, respectively. Thus, for these sources we
have also performed our analysis by assuming a different core radius in
agreement with these observations (i.e. $800$ pc and $3000$ pc, respectively)
and have obtained results similar to the ones reported in Fig. 1. Indeed, the
most likely values of the neutrino flux normalizations change less than a
factor of two, thereby making our conclusions unchanged. We stress that our
main goal is to assess the possibility for future neutrino telescope of
singularly observing these sources, only using their gamma-ray data and the
star formation rate as constraints for their neutrino emissions. In order to
robustly constrain the parameters of these extended sources, it would be
necessary a numerical simulation of their astrophysical environment, which we
leave for future work.

For NGC 4945, NGC 3424, and ARP 220, we find a lower limit for the SFR,
implying the need for additional contributions to the gamma-ray production.
Indeed these sources host an AGN (Yoast-Hull et al., 2017; Ajello et al.,
2020; Xiang et al., 2021a). A special case is ARP 299, for which we find an
upper limit for the SFR. Most of the gamma-ray measurements of this source are
only flux upper limits, thus jeopardizing the likelihood procedure. Finally,
for NGC 1068, we find a good agreement both with the data and the SFR value
given by Kornecki et al. (2020). However, as will be discussed later, the
corresponding neutrino emission predicted by our model is much lower than the
one released by IceCube as a $2.9\sigma$ excess (Aartsen et al., 2020).
Therefore, while the gamma-ray flux could well be explained by the star-
forming activity of the source, the explanation of its potential neutrino
emission might require an intense AGN hot corona activity (Inoue et al., 2020;
Ajello et al., 2020; Kheirandish et al., 2021; Anchordoqui et al., 2021).

##  4 Prospects for multi-messenger observations

Figure 1: Muon neutrino flux normalizations at $1~{}\mathrm{TeV}$ predicted by
the source star-forming activity as a function of source declination. The
diamonds correspond to the most-likely values of the source parameters deduced
by current gamma-ray data (see Table 1), while the bands represent the 68%
credible intervals of the marginal flux distributions. The lines shown the
point-like sensitivity of different neutrino telescopes: 6-year KM3NeT/ARCA
(Aiello et al., 2019), 10-year IceCube (Aartsen et al., 2020), and 10-year
IceCube-Gen2 estimated according to Aartsen et al. (2021).

The results of the likelihood analysis of current gamma-ray data are hence
employed to estimate the high-energy muon neutrino flux from each source. The
normalization of such fluxes at $1~{}\text{TeV}$ are shown in Figure 1, where
the diamonds represent the predictions for the most-likely source parameters
and the bands cover the 68% interval of the marginal flux distributions. In
the plot, our predictions are compared with the IceCube point-like sensitivity
(Aartsen et al., 2020), as well as with the ones of the upcoming neutrino
telescopes KM3NeT/ARCA (Aiello et al., 2019) and IceCube-Gen2 with different
observation times. The latter is estimated to be at least five times better
than the current IceCube one (Aartsen et al., 2021). Remarkably, we find that
the cores of SMC and Circinus galaxy could be potentially detected by KM3NeT
thanks to its higher sensitivity to the southern sky. In the northern sky, a
promising source is ARP 299 which might be within reach of IceCube-Gen2. Other
galaxies such as NGC 4945, M31, NGC 2403 and M82, are instead just a factor of
$\sim 3$ below the KM3NeT and IceCube-Gen2 point-like sensitivities within
their 68% credible intervals. These considerations are likely to be
representative of the real detection potential for these sources. It is worth
noticing that the experimental sensitivities are determined for an $E^{-2}$
neutrino spectrum, and therefore the conclusions might change for softer
spectra. However, as reported in Table 1, the most-likely spectral index for
most galaxies is close to 2. Moreover, the largest values for the neutrino
flux normalization reported in Figure 1 are obtained in case of hard spectra,
thus supporting the comparison with the reported point-like sensitivities.

Figure 2: Most-likely gamma-ray spectral energy distributions of different
galaxies as a function of gamma-ray energy, in agreement with their star-
forming activities and current gamma-ray observations. The left (right) panel
collects the galaxies of the southern (northern) sky. The thick lines show the
CTA differential sensitivity for an observation time of 50 hours (Acharya et
al., 2018).

For the brightest neutrino galaxies, we also show in Figure 2 the most-likely
very-high-energy gamma-ray flux to be compared with the CTA differential
sensitivity, for which we consider an observation time of 50 hours towards the
direction of the galaxy (Acharya et al., 2018). The cutoff of the gamma-ray
flux is due both to internal and external absorption. A few sources are
expected to be above the CTA sensitivity in the range between $100$ GeV and
$10$ TeV; for SMC this range extends to $100$ TeV. This is particularly
relevant, since one of the main limitations in our study is the lack of gamma-
ray measurements in the range above 1 TeV (apart from the VERITAS and H.E.S.S.
data on M82 and NGC 253). Therefore, our results suggest that CTA will allow
us to draw more robust conclusions on the gamma-ray production in SFGs and
SBGs. It is worth emphasizing that, under the assumption of a dominant
hadronic production linked to the star-forming activity, at high energies the
gamma-ray and neutrino emission are directly related each other through the
SFR which is derived from the UV and IR observations.

##  5 Conclusions

In the near future, upcoming neutrino telescopes will potentially observe
nearby star-forming and starburst galaxies as point-like sources. In
particular, they could pose a solid link between the hadronic emission of
these sources, supposed to dominate the GeV-TeV gamma-rays, and the expected
intense star-forming activity as obtained from IR and UV observations.

As the brightest sources predicted by our emission model are in the southern
sky, a leading role will be played by neutrino telescopes located in the
northern hemisphere. They include KM3NeT/ARCA as well as the planned Baikal-
GVD (Avrorin et al., 2020) and P-ONE (Agostini et al., 2020) telescopes. On
the other hand, IceCube-Gen2 will have a better point-like sensitivity in the
northern sky, thus remaining a crucial observer for the sources positioned in
this portion of the sky. In any case, the advent of a Global Neutrino Network
would be necessary to observe most of the point-like emission from SFGs and
SBGs predicted in this work, by increasing the available worldwide effective
volume of neutrino telescopes.

While our results are mainly derived from the star-forming activity, we cannot
exclude a possible AGN emission counterpart for some of the galaxies selected.
However, an additional AGN activity would typically exhibit a flaring
behaviour, unless the related duty cycle is very high. Therefore, the neutrino
measurements over large observational periods would be dominated by the star-
forming steady component, thereby making our detection prospects more robust.
In this regard, a crucial role will be also played by CTA that will
potentially measure the gamma-ray emission above tens of TeV energies for some
of the galaxies investigated. This will allow us to place better constraints
on the sources production mechanism as well as on the parameters of the
emitting cores.

This work indicates that the next decade will be decisive to assess whether
the local star-forming activity can be a tracer of point-like neutrino
production, and it highlights the importance of pivoting the presented
scenarios through the low-energy thermal emissions.

## Acknowledgements

This work was partially supported by the research grant number 2017W4HA7S
“NAT-NET: Neutrino and Astroparticle Theory Network” under the program PRIN
2017 funded by the Italian Ministero dell’Università e della Ricerca (MUR).
The authors also acknowledge the support by the research project TAsP
(Theoretical Astroparticle Physics) funded by the Istituto Nazionale di Fisica
Nucleare (INFN).

## References

  * Aartsen et al. (2018) Aartsen, M., et al. 2018, Science, 361, eaat1378, doi: 10.1126/science.aat1378
  * Aartsen et al. (2020) —. 2020, Phys. Rev. Lett., 124, 051103, doi: 10.1103/PhysRevLett.124.051103
  * Aartsen et al. (2021) Aartsen, M. G., et al. 2021, J. Phys. G, 48, 060501, doi: 10.1088/1361-6471/abbd48
  * Acharya et al. (2018) Acharya, B. S., et al. 2018, Science with the Cherenkov Telescope Array (WSP), doi: 10.1142/10986
  * Agostini et al. (2020) Agostini, M., et al. 2020, Nature Astron., 4, 913, doi: 10.1038/s41550-020-1182-4
  * Aiello et al. (2019) Aiello, S., et al. 2019, Astropart. Phys., 111, 100, doi: 10.1016/j.astropartphys.2019.04.002
  * Ajello et al. (2020) Ajello, M., Di Mauro, M., Paliya, V., & Garrappa, S. 2020, Astrophys. J., 894, 88, doi: 10.3847/1538-4357/ab86a6
  * Albert et al. (2020) Albert, A., et al. 2020, Astrophys. J., 892, 92, doi: 10.3847/1538-4357/ab7afb
  * Alonso-Herrero et al. (2013) Alonso-Herrero, A., et al. 2013, Astrophys. J. Lett., 779, L14, doi: 10.1088/2041-8205/779/1/L14
  * Ambrosone et al. (2021) Ambrosone, A., Chianese, M., Fiorillo, D. F. G., et al. 2021, Mon. Not. Roy. Astron. Soc., 503, 4032, doi: 10.1093/mnras/stab659
  * Anchordoqui et al. (2021) Anchordoqui, L. A., Krizmanic, J. F., & Stecker, F. W. 2021.  https://arxiv.org/abs/2102.12409
  * Avrorin et al. (2020) Avrorin, A., et al. 2020, PoS, ICRC2019, 1011, doi: 10.22323/1.358.1011
  * Bechtol et al. (2017) Bechtol, K., Ahlers, M., Di Mauro, M., Ajello, M., & Vandenbroucke, J. 2017, Astrophys. J., 836, 47, doi: 10.3847/1538-4357/836/1/47
  * Blanco & Linden (2021) Blanco, C., & Linden, T. 2021.  https://arxiv.org/abs/2104.03315
  * Bolatto et al. (2011) Bolatto, A. D., Leroy, A. K., Jameson, K., et al. 2011, The Astrophysical Journal, 741, 12, doi: 10.1088/0004-637x/741/1/12
  * For et al. (2012) For, B.-Q., Koribalski, B. S., & Jarrett, T. H. 2012, Monthly Notices of the Royal Astronomical Society, 425, 1934–1950, doi: 10.1111/j.1365-2966.2012.21416.x
  * Franceschini & Rodighiero (2017) Franceschini, A., & Rodighiero, G. 2017, Astron. Astrophys., 603, A34, doi: 10.1051/0004-6361/201629684
  * Groves et al. (2008) Groves, B., Dopita, M. A., Sutherland, R. S., et al. 2008, The Astrophysical Journal Supplement Series, 176, 438, doi: 10.1086/528711
  * Guo et al. (2019) Guo, X.-L., Xin, Y.-L., Liao, N.-H., & Fan, Y.-Z. 2019, doi: 10.3847/1538-4357/ab44ba
  * H. E. S. S. Collaboration et al. (2018) H. E. S. S. Collaboration, Abdalla, H., Aharonian, F., et al. 2018, A&A, 617, A73, doi: 10.1051/0004-6361/201833202
  * Ha et al. (2020) Ha, J.-H., Ryu, D., & Kang, H. 2020.  https://arxiv.org/abs/2008.06650
  * Hirashita et al. (2003) Hirashita, H., Buat, V., & Inoue, A. K. 2003, Astron. Astrophys., 410, 83, doi: 10.1051/0004-6361:20031144
  * Inoue et al. (2000) Inoue, A. K., Hirashita, H., & Kamaya, H. 2000, Publ. Astron. Soc. Jap., 52, 539, doi: 10.1093/pasj/52.3.539
  * Inoue et al. (2020) Inoue, Y., Khangulyan, D., & Doi, A. 2020, Astrophys. J. Lett., 891, L33, doi: 10.3847/2041-8213/ab7661
  * Kelner et al. (2006) Kelner, S., Aharonian, F. A., & Bugayov, V. 2006, Phys. Rev. D, 74, 034018, doi: 10.1103/PhysRevD.74.034018
  * Kennicutt (1998) Kennicutt, Robert C., J. 1998, Ann. Rev. Astron. Astrophys., 36, 189, doi: 10.1146/annurev.astro.36.1.189
  * Kennicutt & De Los Reyes (2021) Kennicutt, R. C., & De Los Reyes, M. A. C. 2021, The Astrophysical Journal, 908, 61, doi: 10.3847/1538-4357/abd3a2
  * Kennicutt & Evans (2012) Kennicutt, Jr., R. C., & Evans, II, N. J. 2012, Ann. Rev. Astron. Astrophys., 50, 531, doi: 10.1146/annurev-astro-081811-125610
  * Kheirandish et al. (2021) Kheirandish, A., Murase, K., & Kimura, S. S. 2021.  https://arxiv.org/abs/2102.04475
  * Kornecki et al. (2020) Kornecki, P., Pellizza, L. J., del Palacio, S., et al. 2020, Astron. Astrophys., 641, A147, doi: 10.1051/0004-6361/202038428
  * Krumholz et al. (2020) Krumholz, M. R., Crocker, R. M., Xu, S., et al. 2020, Mon. Not. Roy. Astron. Soc., 493, 2817, doi: 10.1093/mnras/staa493
  * Loeb & Waxman (2006) Loeb, A., & Waxman, E. 2006, JCAP, 05, 003, doi: 10.1088/1475-7516/2006/05/003
  * Müller et al. (2020) Müller, A. L., Romero, G. E., & Roth, M. 2020, Mon. Not. Roy. Astron. Soc., 496, 2474, doi: 10.1093/mnras/staa1720
  * Murase et al. (2013) Murase, K., Ahlers, M., & Lacki, B. C. 2013, Phys. Rev. D, 88, 121301, doi: 10.1103/PhysRevD.88.121301
  * Murase et al. (2020) Murase, K., Kimura, S. S., & Meszaros, P. 2020, Phys. Rev. Lett., 125, 011101, doi: 10.1103/PhysRevLett.125.011101
  * Murase & Waxman (2016) Murase, K., & Waxman, E. 2016, Phys. Rev. D, 94, 103006, doi: 10.1103/PhysRevD.94.103006
  * Owen et al. (2021) Owen, E. R., Lee, K.-G., & Kong, A. K. H. 2021, doi: 10.1093/mnras/stab1707
  * Peretti et al. (2019) Peretti, E., Blasi, P., Aharonian, F., & Morlino, G. 2019, Mon. Not. Roy. Astron. Soc., 487, 168, doi: 10.1093/mnras/stz1161
  * Peretti et al. (2020) Peretti, E., Blasi, P., Aharonian, F., Morlino, G., & Cristofari, P. 2020, Mon. Not. Roy. Astron. Soc., 493, doi: 10.1093/mnras/staa698
  * Peretti et al. (2021) Peretti, E., Morlino, G., Blasi, P., & Cristofari, P. 2021.  https://arxiv.org/abs/2104.10978
  * Rahmani et al. (2016) Rahmani, S., Lianou, S., & Barmby, P. 2016, Monthly Notices of the Royal Astronomical Society, 456, 4128, doi: 10.1093/mnras/stv2951
  * Senno et al. (2015) Senno, N., Mészáros, P., Murase, K., Baerwald, P., & Rees, M. J. 2015, Astrophys. J., 806, 24, doi: 10.1088/0004-637X/806/1/24
  * Shimono et al. (2021) Shimono, N., Totani, T., & Sudoh, T. 2021.  https://arxiv.org/abs/2103.08287
  * Sudoh et al. (2018) Sudoh, T., Totani, T., & Kawanaka, N. 2018, Publ. Astron. Soc. Jap., 70, Publications of the Astronomical Society of Japan, Volume 70, Issue 3, 1 June 2018, 49, https://doi.org/10.1093/pasj/psy039, doi: 10.1093/pasj/psy039
  * Tamborra et al. (2014) Tamborra, I., Ando, S., & Murase, K. 2014, JCAP, 09, 043, doi: 10.1088/1475-7516/2014/09/043
  * Thompson et al. (2006) Thompson, T. A., Quataert, E., Waxman, E., Murray, N., & Martin, C. L. 2006, Astrophys. J., 645, 186, doi: 10.1086/504035
  * Torres (2004) Torres, D. F. 2004, Astrophys. J., 617, 966, doi: 10.1086/425415
  * VERITAS Collaboration et al. (2009) VERITAS Collaboration, Acciari, V. A., Aliu, E., et al. 2009, Nature, 462, 770, doi: 10.1038/nature08557
  * Werhahn et al. (2021a) Werhahn, M., Pfrommer, C., & Girichidis, P. 2021a.  https://arxiv.org/abs/2105.12134
  * Werhahn et al. (2021b) Werhahn, M., Pfrommer, C., Girichidis, P., & Winner, G. 2021b, doi: 10.1093/mnras/stab1325
  * Xiang et al. (2021a) Xiang, Y.-c., Deng, J.-h., & Jiang, Z.-j. 2021a.  https://arxiv.org/abs/2105.08240
  * Xiang et al. (2021b) Xiang, Y.-c., Jiang, Z.-j., & Tang, Y.-y. 2021b.  https://arxiv.org/abs/2105.02594
  * Yoast-Hull et al. (2017) Yoast-Hull, T. M., Gallagher, J. S., Aalto, S., & Varenius, E. 2017, Mon. Not. Roy. Astron. Soc., 469, L89, doi: 10.1093/mnrasl/slx054
  * Yuan et al. (2011) Yuan, F.-T., Takeuchi, T. T., Buat, V., et al. 2011, PASJ, 63, 1207, doi: 10.1093/pasj/63.6.1207