1 Introduction↩︎

Stars radiate energy by emitting neutrinos, in addition to electromagnetic emission. By means of nuclear fusion in the core of stars, light atomic nuclei combine to form heavier ones. Thermonuclear neutrinos are produced as a by-product with \(\mathcal{O}(0.1\)–\(10\)) MeV energy (see, e.g., Refs. [1]–[3] for an overview of thermonuclear processes). Stars also emit thermal neutrinos from thermal processes that strongly depend on the temperature, density, and composition of the medium; thermal neutrinos have typical energies of \(\mathcal{O}(10\)–\(100)\) keV (see e.g. Refs. [3]–[5] for an introduction to thermal processes in stars). Hence, all stars in our Universe are neutrino sources, including those in our own Galaxy–the Milky Way. The flux of Galactic neutrinos from stars was first estimated in Ref. [6], relying on the FRANEC evolutionary code [7] and a Galactic model from Ref. [8]. Similar results for thermonuclear neutrinos were presented in Ref. [9] for other models accounting for the structure and evolution of our Galaxy, relying on the stellar evolution results and neutrino emission properties presented in Ref. [6].

Our understanding of stellar evolution and neutrino production in stellar environments has greatly improved in the past decade; we refer the reader to, e.g., Refs. [10], [11] for recent reviews on the topic. Solar neutrinos from different thermonuclear channels, such as the pp-chain [12]–[15] and the CNO cycle [16], have been experimentally identified. In this regard, solar neutrino spectroscopy has been fundamental to improve our understanding of the Sun, see e.g. Refs. [17]–[21], as well as to learn about neutrino properties [22], [23]. In addition, we have gained a deeper understanding of the dependence of neutrino emission on stellar mass and metallicity [24]–[26].

In this work, we forecast of the Galactic stellar neutrino flux (GS\(\nu\)F): the cumulative neutrino flux from all stars in the Milky Way expected at Earth. Our findings are summarized in Fig. 1. We rely on Ref. [27], that employs the Gaia Data Release 2, for modeling the spatial distribution of stars in our Galaxy and use new insights from galactic archaeology [28] to model the Galactic star formation history (SFH). Moreover, we compute the neutrino emission from stars with zero-age-main-sequence (ZAMS) mass between \(0.08\;M_\odot\) and \(100\;M_\odot\), following the stellar evolution from the pre-main sequence stage to the final evolutionary stages.

This paper is organized as follows. In Sec. 2, we describe our modelling of the Milky Way stellar population, including the spatial distribution of stars and the SFH. In Sec. 3, we outline our modelling of the evolution of each star in our suite. Section 4 describes the approach adopted to compute the GS\(\nu\)F. Section 5 presents our results for low-mass stars (\(M < 0.9\;M_\odot\)), intermediate-mass stars (\(0.9\;M_\odot \lesssim M < 8\;M_\odot\)), and high-mass stars (\(M \geq 8\;M_\odot\)); the main GS\(\nu\)F features are also discussed. Finally, we summarize our findings and discuss the GS\(\nu\)F detectability prospects in Sec. 6.

2 Galactic star formation rate↩︎

In this section, we present the model of the Galactic stellar population and our parametric description of the birthrate of stars.

2.1 Spatial distribution of stars↩︎

As sketched in the top left panel of Fig. 2, we divide the baryonic mass distribution in our Galaxy in three independent components: the Galactic buldge, a thin disk, and a thick disk. We describe the mass distribution of stars in the Milky Way in Galactocentric coordinates \((R, z)\).

For the bulge, we adopt a parameterization of the baryonic density that does not account for the complex stellar distribution in the bulge and neglects the bar structure [29], [30]: \[\begin{align} \rho_{\rm bulge} = \frac{\rho_{0, \textrm{bulge}}}{(1 + r'/r_0)^\alpha} \exp \left[-(r'/r_{\rm cut})^2\right]\, , \label{eq:bulge} \end{align}\tag{1}\] with \(r' = \sqrt{R^2 + (z/q)^2}\) ¹. Here, \(q\), \(r_0\), \(r_{\rm cut}\), and \(\alpha\) are phenomenological parameters determining the three-dimensional shape and size of the bulge, and \(\rho_{0,{\rm bulge}}\) is the baryon density of the bulge.

The baryon density of the thin and thick disks is modeled in three dimensions as follows [31], [32]: \[\begin{align} \rho_{d} = \frac{\Sigma_0}{2z_d} \exp\left( -\frac{|z|}{z_d} - \frac{R}{R_d}\right) \, ; \label{eq:disk} \end{align}\tag{2}\] the parameters \(z_d\) and \(R_d\) determine the size of each of the disks, and \(\Sigma_0\) is the corresponding density.

Table 1 summarizes the values adopted for the parameters introduced in Eqs. 1 and 2 , following Ref. [27]. The parameterization proposed in Ref. [27] relies on Gaia Data Release 2 [33], assumes a Navarro-Frenk-White dark matter distribution [34], an axisymmetric model of the Milky Way with no spiral arms, and accounts for two additional gas disks as well as the circumgalactic medium. Under these assumptions, the baryonic mass of the bulge is \(M_{\rm bulge} = 9.2 \times 10^9~M_\odot\), the baryonic mass of the thin disk is \(M_{\rm d, thin} = 3.98 \times 10^{10}~M_\odot\), and the one of the thick disk is \(M_{\rm d, thick} = 1.07 \times 10^{10}~M_\odot\) [27].

We implement the arm structure of the Milky Way of the thick and thin disks, specifically the arms Scutum–Centaurus, Sagittarius–Carina, Perseus, and Norma-Outer following Ref. [35]. However, we do not account for the gas disks and the intergalactic medium. Those additional components have a small density of stars and, therefore, are expected to provide a negligible contribution to the GS\(\nu\)F.

2.2 Birthrate function↩︎

The stellar birthrate function, \(B\), provides the number of stars born in the time interval \((t, t+{\rm d}t)\) and mass interval \((M, M+{\rm d}M)\) in the Milky Way. It can be expressed as the product of the SFH (\(\zeta(t)\)) and the initial mass function (IMF, \(\xi(M)\)): \[\begin{align} B(t,M)\,{\rm d}t \,{\rm d}M= \zeta(t)\,\xi(M)\,{\rm d}t \,{\rm d}M\, ; \end{align}\] note that we assume that the IMF is constant in time and the star formation rate does not depend on \(M\).

2.2.1 Star formation history↩︎

Different methods have been proposed to investigate the star formation rate of the Galaxy based on its chemical evolution, such as the serial and parallel approaches [28]; in the former, the formation of the Milky Way is modeled through continuous accretion of gas during which the halo, thick and thin disks are formed, or assuming several episodes of gas infall. The latter approach considers that Galaxy formation results from different episodes of gas accretion occurring in parallel, but at a different rate.

In this work, we adopt an effective two-infall model that parametrizes the star formation rate of the Milky Way following Ref. [36] ². The two-infall model assumes that an initial collapse led to the formation of the halo of the thick disk. Star formation continued at a very efficient rate (much larger than the current one) until the gas density dropped below a certain threshold. The gas lost by the halo accumulated in the center, forming the bulge (cf. left top panel of Fig. 2). Then, a second infall event is responsible for the thin disk. This was either the result of a merger event with a small galaxy or due to the longer time required for material with a high angular momentum to accrete. A recent analysis supports the occurrence of a delayed episode of gas infall [39].

We consider the first infall epoch, which resulted in the formation of the bulge and the thick disk, and a later infall forming the thin disk (see the top left panel of Fig. 2). We parametrize the SFH during the two-infall episodes as \[\begin{align} \zeta(t) \propto t^\kappa \exp(- t/\tau)\, , \label{eq:sfh} \end{align}\tag{3}\] and determine the effective parameters \(\kappa\) and \(\tau\) for each infall episode by fitting the two-infall best-fit model from Ref. [40]. For the SFH of the bulge and thick disk, we take \(\kappa =0.4\) and \(\tau = 0.5\) Gyr, while for the thin disk SFH, we assume \(\kappa=0.5\) and \(\tau = 6.2\) Gyr.

The top right panel of Fig. 2 shows the SFH for each of the two infalls, using the parameters above. Stars in the bulge and thick disk formed rapidly (\(12\)–\(13\) Gyr ago), followed by stars in the thin disk forming at a roughly constant rate for the past \(12\) Gyr. Hence, the bulge and thick disk host older stars, whereas the stellar population of the thin disk should exhibit a wider spread in stellar age.

2.2.2 Initial mass function↩︎

The IMF describes the distribution of stars at birth as a function of stellar mass (see, e.g., Ref. [41] for a recent review). We adopt a broken power-law for the IMF [42]: \[\begin{align} \label{eq:imf} \xi(M) &\propto M^{-\alpha} \\&\text{with}\begin{cases} \alpha = 1.3 \qquad 0.08\;M_\odot\leq M < 0.5\;M_\odot \\ \alpha = 2.3 \qquad 0.5\;M_\odot \,\,\,\leq M \leq 100\;M_\odot \\ \end{cases}\, . \end{align}\tag{4}\] The lower limit in mass of the IMF is chosen according to the fact that, if a proto-star forms with a mass smaller than \(0.08\;M_\odot\), its internal temperature never becomes high enough for thermonuclear fusion to begin, leading to a brown dwarf. The choice of the IMF upper limit in this work is arbitrary, since the IMF strongly suppresses the population of high-mass stars. The bottom panel of Fig. 2 illustrates the IMF in Eq. 4 .

3 Stellar evolution↩︎

In this section, we present our modeling of the evolution of stars with ZAMS mass of \([0.08, 100]\;M_\odot\). We also illustrate the approach adopted to reconstruct the spectral energy distribution of neutrinos.

3.1 Stellar models↩︎

3.1.1 Stellar population modeling↩︎

We construct a suite of stellar models, considering each star with fixed ZAMS mass between \([0.1, 100]\;M_\odot\) as a single, non-rotating object with solar metallicity. We divide the population of stars in three groups: low-mass stars with \(M < 0.9\;M_\odot\), intermediate-mass stars with \(0.9\,M_\odot\leq M < 8\,M_\odot\), and high-mass stars with \(M \geq 8\,M_\odot\).

As illustrated in the bottom panel of Fig. 2, we select a discrete subset of \(13\) progenitors with ZAMS \(M/M_\odot = 0.1\), \(0.2\), \(0.4\), \(0.6\), \(0.8\), \(1\), \(1.2\), \(2\), \(3\), \(12\), \(15\), \(20\), and \(40\). We evolve each of these models using Modules for Experiments in Stellar Astrophysics (MESA) 24.08.1 [43]–[47], adapting existing models from the MESA Test Suite. We assume that each of these progenitors is representative of a range of ZAMS masses, as shown in the bottom panel of Fig. 2. The models with mass \(\lesssim 1\;M_\odot\) are evolved from their pre-main sequence evolutionary stage until their terminal age main sequence (i.e., until hydrogen is approximately exhausted, since their main-sequence lifetime would otherwise be larger than the age of the Universe). Intermediate-mass progenitors (\(1\;M_\odot \leq M \leq 3\;M_\odot\)) are evolved until the white-dwarf phase, when the photon luminosity significantly exceeds the neutrino luminosity. Finally, the high-mass progenitors are evolved until pre-core-collapse conditions are reached.

3.1.2 Convection↩︎

For each stellar model, we adopt a mixing-length treatment of convection, following Refs. [43], [48]. We fix the mixing-length parameter (\(\alpha_{\rm ML}\)), relying on a solar calibration. Other parameters related to the mixing-length velocity and temperature gradient in a rising bubble are set to default as in Ref. [43]. We also account for convective overshooting with an exponential parametrization of the diffusion coefficient, as illustrated in Ref. [49], with \(f_{\rm{ov}} = 0.016\) [24] ³.

3.1.3 Solar calibration↩︎

We calibrate our suite of models to solar data to fix the initial hydrogen, helium, and heavy-element mass fraction (X, Y, and Z, respectively), as well as the mixing length parameter \(\alpha_{\rm ML}\). To this end, we use the MESA simplex built-in routine. We find the set of parameters that best reproduces the solar luminosity (\(L_\odot = 3.828 \times 10^{33}\) erg s\(^{-1}\) [50]), effective temperature (\(T_{\rm eff} = 5772\) K [50]), helium surface fraction (\(Y_{\rm surface} = 0.2485 \pm 0.0035\) [51]), the radius of the convective zone (\(R_{\rm cz} = (0.713 \pm 0.001) R_\odot\) [52]), and the measured sound-speed profile [53].

We perform calibration for two solar models with different predicted surface abundances of heavy elements: \((Z/X)_{\rm surface} = (1.8 \pm 0.1)\times 10^{-2}\) for the AGSS09 solar model [54] and \((Z/X)_{\rm surface} = (2.3 \pm 0.1)\times 10^{-2}\) for the GS98 solar model [55], for the estimated age of the Sun (\(4.61 \times 10^9\) Gyr) [56], [57]. We have tested that the predicted neutrino flux from the pp-chain, which dominates the solar neutrino emission, is in agreement with observations for both solar models. Hereafter, we adopt the calibration values for the GS98 model (\(X_0 = 0.7020\), \(Y_0 = 0.2778\), \(Z_0 = 0.202\), and \(\alpha_{\rm ML} = 1.96\)), which yield a solar neutrino luminosity \(L_{\nu, \odot} = 0.02469 L_{\odot}\), where \(L_{\odot}\) is the photon luminosity.

3.1.4 Nuclear network, thermonuclear and thermal processes↩︎

We rely on the built-in MESA nuclear reaction networks, mesa_49 for low- and intermediate-mass stars and mesa_206 for the high-mass stars. The nuclear reaction rates are taken from the Joint Institute for Nuclear Astrophysics REACLIB library [58] and account for plasma screening effects [59] (we refer the interested reader to Refs. [46], [60] for details on the calculation of the nuclear reaction rates). In addition, tabulated weak reaction rates are from Refs. [61]–[63].

Table [tab:reactions] summarizes some of the thermonuclear and thermal reactions, taking place during stellar evolution and that lead to neutrino production. We take these reaction channels into account for computing the GS\(\nu\)F, as illustrated below.

3.2 Neutrino energy spectrum↩︎

For each time snapshot of stellar evolution, MESA provides the reaction rates for each thermonuclear process, see Table [tab:reactions], as well as the total energy released. The neutrino energy spectrum is not provided as MESA output, but can be computed considering the kinematics of each reaction.

Neutrinos from \(\beta\) processes occurring as part of the pp and CNO cycles produce a continuous spectrum that can be computed knowing the Q-value of the reaction; whereas reactions with two bodies in the final state produce spectral lines. As for the thermonuclear component of the neutrino flux, we compute the energy spectrum for all pp-chain reactions (i.e. the continuous neutrino spectra from pp and hep neutrinos, \(^8\)B decay, and the spectral lines of pep neutrinos and \(^7\)Be neutrinos), the CNO-cycle (i.e. the energy spectra from the decay of \(^{13}\)N, \(^{15}\)O, and \(^{17}\)F), and the spectrum from \(^{18}\)F decay. For a detailed discussion on the modelling of the spectral shape of pp-chain and CNO-cycle reactions, we refer the reader to Ref. [10] and references therein. An analogous procedure approximately holds for \(\beta\) processes, such as \(^{18}\)F.

As for thermal processes, MESA outputs the reaction rates as functions of radius and for several time steps of the evolution, for plasmon decay, pair annihilation, bremsstrahlung, photoneutrinos, and neutrinos from recombination. For such processes, the neutrino emission strongly scales as a function of e.g. the electron temperature and electron number density (see e.g. Ref. [4]). Hence, we extract the radial profiles of the relevant thermodynamic quantities from MESA at the relevant time snapshots during stellar evolution; relying on such profiles, one can compute the energy spectrum of neutrinos from plasmon decay [64], [65], pair annihilations [66], [67], bremsstrahlung  [68], recombination [69], [70], and photoneutrinos [71]. As discussed in Sec. 5, we only account for plasmon decay as contributing to white dwarfs cooling, and we refer the reader to Sec. 5 for details on the relevance of the remaining thermal processes.

In thermonuclear reactions, only electron neutrinos are produced. However, thermal processes produce neutrino-antineutrino pairs. Photoneutrinos and recombination neutrinos are always of electron flavor, while all the other thermal processes produce a mixture of electron-type and non-electron type neutrino pairs, where the ratio is fixed by the structure of weak interactions. When traveling towards the outer layers of stars, and depending on their energy and the radial profile of the electron number density at each stellar evolution stage, neutrinos undergo flavor conversion [72], [73]. Flavor conversion in the source, together with loss of coherence en route to Earth, determines the flavor composition of neutrinos reaching Earth. We do not account for flavor conversion due to the complexity of integrating the energy-dependent and time-dependent flavor evolution for each progenitor. Therefore, hereafter we report the total neutrino and antineutrino flux for all six flavors.

4 Galactic stellar neutrino flux↩︎

The GS\(\nu\)F at Earth and at present time is: \[\begin{align} \Phi_\nu = \int \text{d}^3 \vec{r}\,&\frac{N_* (\vec{r})}{ 4\pi D_\oplus^2(\vec{r})}\int \text{d}t\,\zeta(t, \vec{r})\,\nonumber \\&\times \int \text{d}M\,\xi(M) \frac{{\rm d} \mathcal{\phi}_\nu}{{\rm d}E_\nu} (M, \tau)\, . \label{eq:nu-flux} \end{align}\tag{5}\] Here, \(N_*\) denotes the number of stars at the location \(\vec{r}\) expressed in Galactic coordinates and corresponding to the distance from Earth \(D^2_\oplus(\vec{r})\); \(\zeta\) denotes the SFH normalized to unity (introduced in Eq. 3 ), \(\xi\) is the IMF normalized to unity (defined as in Eq. 4 ) under the assumption of uniform IMF across the Milky Way, \(\tau = t_U -t\) with \(t_U = 13.5\) Gyr being the age of the Universe. The neutrino energy flux emitted from a star of ZAMS mass \(M\) and age \(\tau\) is \({\rm d}\mathcal{\phi}_\nu/{\rm d}E_\nu\).

For simplicity, we treat the bulge, thin disk, and thick disk separately. This allows us to break Eq. 5 into three integrals with independent SFH (cf. Sec. 2.2.1). As for the spatial distribution of stars, we construct a finely binned spatial grid for each of the Milky Way components. Under the assumption of isotropic IMF, the number density of stars at each discrete volume element of the 3D grid (\(V_k\)) is the ratio between the baryon density of each Galactic component (see Table 1) and the mean mass, \(\overline{M}\), according to the IMF, i.e. \[\begin{align} \overline{M} = \frac{\int_{0.08\,M_\odot}^{100\,M_\odot}\,{\rm d}M \, \xi(M) M }{\int_{0.08\,M_\odot}^{100\,M_\odot}\,{\rm d}M \,\xi(M)}\,. \end{align}\] For each 3D volume element, we draw a random position (\(\vec{r}_k\)) and assume that all stars in that volume element are located in \(\vec{r}_k\). Then, the spatial integral in each of the elements of the spatial grid is \[\begin{align} \int_{V_k} {\rm d}^3\vec{r} \frac{N_* (\vec{r})}{4\pi D^2_\oplus (\vec{r})} \simeq \frac{1}{4\pi D^2_{\oplus}(\vec{r}_k)} \frac{\int_{V_k} {\rm d}^3 \vec{r} \rho_{n}(\vec{r})}{\overline{M}} \equiv \varepsilon_{k}^{(n)}\, , \end{align}\] where \(\rho_n\) is the baryon density of the bulge (Eq. 1 ) and each of the disks (Eq. 2 ). The index \(n\) indicates each of the three components, i.e. \(n = \lbrace\)bulge, thin disk, thick disk\(\rbrace\).

We also consider a grid of stellar ages for which we compute the average neutrino emission. We use the normalized SFH distribution and IMF as weights: \[\begin{align} w_{ij} = \frac{\int_{M_{i,{\rm min}}}^{M_{i,{\rm max}}}\, {\rm d}M \, \xi(M) }{\int_{0.08\,M_\odot}^{100\,M_\odot}\, {\rm d}M\,\xi(M) } \frac{\int_{t_{j,{\rm min}}}^{t_{j,{\rm max}}}\,{\rm d}t\, \zeta(t_U-t) }{\int_{0}^{t_U}\,{\rm d}t\, {\zeta}(t_U- t) }\, . \end{align}\] Here \(M_{i,{\rm min}}\) and \(M_{i,{\rm max}}\) indicate the minimum and maximum mass of the \(i\)-th mass bin. Note that since the bulge, thin and thick disks have different SFHs (see Sec. 2.2.1), we compute such weights for each of the three Galactic components separately.

Using the discretization method outlined above, Eq. 5 becomes \[\begin{align} \label{eq:SnuB} \Phi^{(n)}_{\nu} = \sum_k \varepsilon^{(n)}_{k} \sum_{i,j} w^{(n)}_{ij} \frac{d\mathcal{\phi}_\nu}{dE_\nu}(M_i, \tau_j)\, . \end{align}\tag{6}\] The cumulative neutrino flux is the sum of the neutrino emission from the three Galaxy components.

In our numerical modeling, we consider \(4331\) volume elements for the bulge, and \(50601\) volume elements for each of the disks. For low- and intermediate-mass stars, our stellar-age grid has \(2700\) elements linearly spaced between \(5\) and \(13.5\) Gyr. For high-mass stars, we adopt \(4000\) bins in stellar age between \(5000\) yr and \(20\) Myr. We have performed resolution tests for the stellar age grid and the spatial grid to ensure convergence of our results (plots not shown here).

5 Results: stellar neutrino flux↩︎

In this section, we present our findings on the GS\(\nu\)F from low-, intermediate- and high-mass stars in our Galaxy. We discuss the reactions dominating the neutrino emission for each class and discriminate among neutrino and antineutrino emission from the bulge, the thin disk, and thick disk.

5.1 Low-mass stars↩︎

For low-mass stars (we refer the reader to e.g. Refs. [74], [75] for details on the evolution of low-mass stars including post main-sequence stages), we follow the evolution of stars with ZAMS mass of \(0.1\,M_\odot\), \(0.2\,M_\odot\), \(0.4\,M_\odot\), \(0.6\,M_\odot\) and \(0.8\,M_\odot\), as outlined in Sec. 3.1. We evolve such models until the terminal-age main sequence, i.e. until hydrogen is close to being exhausted and the central mass fraction of \(^1\)H reaches \(0.01\%\). We halt the evolution of these models before the end of the main sequence because the lifetime of these stars as main sequence stars would be larger than the age of the Universe. We then compute the neutrino emission by means of Eq. 6 . Our results concerning the evolution of low-mass progenitors is in agreement with the existing literature [25].

The top panel of Fig. 3 displays the GS\(\nu\)F at Earth from low-mass stars in our Galaxy. The neutrino luminosity of low-mass stars increases as a function of the ZAMS mass and as the star ages in its main sequence phase. Therefore, the main contribution to the GS\(\nu\)F from Galactic low-mass stars comes from the oldest and most massive stars of this subpopulation. The leading neutrino emission channel is the pp chain, although neutrinos from the CNO cycle tend to dominate towards the end of their main sequence lifetime. Lines emerge from the continuum spectrum in Fig. 3; such spectral lines are characteristic of the \(^7\)Be and pep reactions.

The GS\(\nu\)F from low-mass stars shown in the top panel of Fig. 3 is divided according to the contribution from the bulge, thin and thick disks, from left to right, respectively. The emission from the thin disk dominates since the latter has the largest stellar population. Notice, however, that a bump is expected in the spectrum at \(1\)–\(1.6\) MeV from the CNO cycle and due to the contributions from the bulge and thick disk, as visible from Fig. 1 (orange curve). The bulge and the thick disk are the oldest in the Milky Way, therefore a fraction of their stellar population has already reached the temperature threshold for the CNO cycle to become the dominant neutrino production channel.

5.2 Intermediate-mass stars↩︎

Main neutrino production channels↩︎

Our sample of intermediate-mass stars consists of progenitors with ZAMS mass \(M = 1\,M_\odot\), \(1.2\,M_\odot\), \(2\,M_\odot\), \(3\,M_\odot\). We explore the evolution of these stellar models employing MESA and then compute the neutrino energy spectrum for this population, as outlined in Sec. 3.1 and Eq. 6 . For these stars, the CNO cycle dominates the hydrogen burning. As for low-mass stars, the neutrino luminosity increases as a function of the ZAMS mass and stellar age during the main sequence.

During the post-main-sequence evolution, since the core temperature of these stars is not large enough (\(\lesssim 10^9\) K) to fuse carbon, an inert core of carbon and oxygen builds up. The star then loses its outer layers to stellar winds, forming a planetary nebula and a white dwarf. White dwarfs lack a source of heat; electron degeneracy supports them against gravitational collapse. Hence, first they cool down by means of neutrino emission (mainly via plasmon decay, which we compute following Ref. [64]); then, at later stages, radiative cooling takes over. We have tested that our stellar models are in agreement with existing literature, including recent work focusing on neutrino emission [24].

The middle panels of Fig. 3 show the GS\(\nu\)F from intermediate-mass stars in our Galaxy. For the oldest Galactic components (the bulge and the thick disk), the neutrino emission from the smaller ZAMS masses is the largest, since these are the only stars that are still in the main sequence to date. The thermal neutrino emission from the cooling of white dwarfs is also led by intermediate-mass stars with smaller mass due to their abundance, as implied by the IMF and shown in Fig. 2. For the thin disk, the contribution from thermonuclear reactions is similar for all masses in this range, since there is a significant part of the population in the main-sequence stage. These differences arise from the relative importance between the pp chain and the CNO cycle during hydrogen burning according to the core temperature. The thermal part of the neutrino energy spectrum is coming from plasmon decay and it is dominated by the highest ZAMS masses in the low-energy range.

Additional neutrino production channels↩︎

At the end of the main sequence, stars with mass \(\lesssim 2\;M_\odot\) have a degenerate core. When the mass of the degenerate core reaches \(\sim 0.5\;M_\odot\), the core temperature is about \(10^8\) K and helium ignition is activated via triple-\(\alpha\) process. Helium ignition is responsible for an electromagnetic flash. The helium flash is sometimes followed by subflashes; such flashes become less likely as the progenitor mass increases and burns a smaller amount of helium. As a consequence, the nuclear energy that is generated mainly goes into lifting electron degeneracy in the core. This process leads to the formation of a stable convective core of helium burning under non-degenerate conditions.

The slowest step in the CNO cycle is proton capture on \(^{14}\)N. As a result, when hydrogen burning is completed, one expects a significant abundance of \(^{14}\)N. During helium ignition, \(^{14}\)N is converted into \(^{22}\)Ne by the following set of reactions, \[\begin{align} &^{14}{\rm N} \,+\, \alpha \to \gamma \,+\, ^{18}{\rm F}\, ,\\ &^{18}{\rm F} \to e^+ + \mathbf{\nu_e} \,+ \,^{18}{\rm O}\, ,\\ & ^{18}{\rm O} \,+\, \alpha \to \gamma \,+\,^{22}{\rm Ne}\, . \end{align}\] Hence, the helium flash (and the subflashes) are correlated with the emission of a burst (or multiple bursts) of neutrinos from \(^{18}\)F decay (with a Q-value of \(1.655\) MeV) [76]. Despite their large neutrino luminosity, these flashes barely affect the time-integrated energy emitted in neutrinos from each intermediate-mass star during its life. In addition, the end point of the energy spectrum of neutrinos from \(^{18}\)F is very close to the one of neutrinos from \(^{15}\)O decay (\(E_\nu \simeq 1.73\) MeV) from the CNO cycle. Hence, the neutrino energy spectra from these two channels tend to overlap in energy.

For masses \(M > 2\;M_\odot\), helium ignition occurs in the core under non-degerate conditions. Hence, no helium flash (nor neutrino burst) is expected.

At later stages in the post-main-sequence evolution, during the asymptotic giant branch phase, hydrogen and helium burning occur in geometrically thin shells on top of the carbon-oxygen core. Thermal pulses occur when the helium shell ignites. These episodes also produce a burst of neutrinos from \(^{18}\)F decay. Such neutrino bursts are less luminous than the ones following helium flashes. We stress, however, that this phase of stellar evolution is very uncertain, therefore existing models are calibrated to observational data [77], [78]. Importantly, the number of thermal pulses depends on the treatment of mass loss [79]. The latter may halt pulses once the envelope has been ejected.

In our models, a parametric treatment of convection has been adopted, based on mixing-length theory, and calibrated to match observations in the Sun. However, for masses larger than \(1.5\)–\(2\;M_\odot\), stars develop a convective core and a radiative envelope, contrary to the case of solar-mass stars. Hence, the number of thermal pulses in our models, and even their existence, might be a result of the approximations intrinsic to our model. Despite that, given the short duration of such pulses, their contribution to the GS\(\nu\)F is negligible.

During the hydrogen- and helium-burning stages of evolution, the stellar-core temperature is such that neutrinos from plasmon decay and photoneutrinos are also emitted. The neutrino flux from these reaction channels was estimated to be comparable or smaller than the number of neutrinos emitted during white-dwarf cooling [6]. Our forecast of the GS\(\nu\)F does not account for these thermal processes in these stages of stellar evolution; in this sense, our prediction is conservative. These thermal processes would partially modify the thermal component of the spectrum (below \(100\) keV) in shape, enhancing it. Given the strong dependence of these processes on temperature, a detailed assessment of their impact is beyond the scope of this work; in fact, a more accurate treatment of convection in the core would be needed to this purpose.

5.3 High-mass stars↩︎

Main neutrino production channels↩︎

Our discrete sample of high-mass stars, evolved using MESA as outlined in Sec. 3, consists of four progenitors with ZAMS mass of \(12\;M_\odot\), \(15\;M_\odot\), \(20\;M_\odot\), and \(40\;M_\odot\). We follow the evolution of these progenitors until they reach conditions close to core collapse–namely when the material inside the iron core exceeds \(10^7\) cm/s.

During their life as main sequence stars, neutrino emission is dominated by the CNO cycle. At later stages, thermal neutrino emission, mainly from pair processes becomes relevant. This efficient energy-loss mechanism sets a rapid evolution (years to hours) of the pre-supernova stars during the advanced stages of nuclear fusion. Hence, neutrino emission during the main-sequence evolution dominates the contribution to the GS\(\nu\)F.

The bottom panels of Fig. 3 show the GS\(\nu\)F from high-mass stars, computed based on Eq. 6 . Despite the smaller abundance of high-mass stars dictated by the IMF (cf. Fig. 2), the neutrino emission grows with the ZAMS mass in the range below \(1.6\) MeV and is mainly due to the CNO cycle. For higher neutrino energies, the spectrum is shaped by pp-neutrinos, \(^8\)B and hep neutrinos, from progenitors with mass \(\lesssim 15\,M_\odot\) within this population, i.e. from the less massive stars of this subpopulation.

High mass stars have a very short main-sequence lifetime. Therefore, the contribution to the GS\(\nu\)F is smaller for the older components of the Milky Way (i.e., the bulge and thick disk). However, since star formation in the thin disk is a continuous process, occurring even in recent times, we find that there is a significant population of high-mass stars in the main-sequence stages of evolution at present; this population of main-sequence stars emits a copious amount of neutrinos.

We note that our treatment of convection calibrated to the Sun might not be accurate for high-mass stars with a convective core. This could affect the late stages of stellar evolution. However, such stages have a negligible impact on the GS\(\nu\)F due to their short duration and overall neutrino emission smaller than the one integrated during the remaining lifetime of the massive star.

Additional neutrino production channels↩︎

As the pre-supernova phase approaches, \(\beta\)-processes contribute to increase the electron-to-baryon ratio, making the medium more neutron rich and favoring the emission of a significant amount of MeV neutrinos [26], [80]–[82]. Thermal neutrino emission from the pre-supernova phase should be expected [67], [83]. Nevertheless, we do not account for the contribution of pre-supernova neutrinos to the GS\(\nu\)F. This is justified because about \(30\) supernova candidates have been identified at a distance smaller than \(1\) kpc from Earth [84]; their pre-supernova emission is expected to increase substantially only less than a day before the core collapse of the massive star. As a consequence, we estimate that the cumulative neutrino emission from these local supernova candidates is negligible with respect to the total neutrino emission from the remaining massive stars in our Galaxy.

Notice that our approach yields a local rate of \(2\)–\(3\) stars with \(M \gtrsim 8\, M_\odot\) undergoing core collapse per century; this is compatible with more sophisticated estimates [85]. Since the number of core-collapse supernovae in the Galaxy at present is \(\ll 1/\)s, we conclude that supernovae do not contribute to the GS\(\nu\)F. Instead, a local core-collapse supernova would be observed at neutrino telescopes as a transient signal [11], [86], [87].

For high-mass stars, it is known that photoneutrino production is the dominant thermal process during main-sequence evolution, especially during the helium burning phase. Previous work evaluated this contribution to be smaller than or comparable to the one resulting from the cooling of a white dwarf [6]. For reasons analogous to the ones discussed in Sec. 5.2, we do not include this contribution in our prediction.

5.4 Cumulative neutrino flux↩︎

The GS\(\nu\)F is displayed in Fig. 1. For readers wishing to reproduce their own version of the GS\(\nu\)F plot, we provide the tabulated fluxes of the different GS\(\nu\)F components displayed in Fig. 1 as Supplemental Material of this paper (and ancillary data files in the arXiv repository). In addition, Fig. 4 highlights the contributions to the GS\(\nu\)F according to the distance of stars from Earth and the overall contributions from the bulge, thin and thick disks (in the left and right panels, respectively). For energies smaller than \(100\) keV, the GS\(\nu\)F is dominated by thermal processes, in particular plasmon decay from white dwarfs marking the end of the evolution of intermediate-mass stars; see also Fig. 3. For larger energies, thermonuclear processes dominate the GS\(\nu\)F. Overall, the flux peaks at energies \(\lesssim 0.5\) MeV; this emission is mostly due to pp neutrinos produced in low- and intermediate-mass stars. For energies up to \(\sim 1.6\) MeV, the GS\(\nu\)F is shaped by CNO reactions occurring in the core of intermediate- and high-mass stars. Therefore, the spectral lines from \(^7\)Be and pep neutrinos displayed in Fig. 1 for low-mass stars (orange curve) do not leave any trace in the overall emission plotted as black dashed lines in Fig. 4. The GS\(\nu\)F is significantly reduced at high energies (\(\gtrsim 1.5\) MeV), with a spectral shape determined by \(^8\)B and hep neutrinos.

The left panel of Fig. 4 suggests that stars within \(6\) and \(10\) kpc from Earth are the main contributors to the GS\(\nu\)F. Notably, this distance range includes the Galactic Center (located at roughly \(8\) kpc from Earth), where the number density of stars is the largest.

The right panel of Fig. 4 highlights that the largest contribution to the GS\(\nu\)F comes from the thin disk of the Milky Way. The reason is twofold. First, the thin disk is the most massive component (and hence the one with the largest number of stars, under the assumption of isotropic IMF). Second, the thin disk is the youngest component of the Galaxy and therefore the one hosting the largest number of intermediate and high-mass stars in the main sequence (with the main sequence being the evolutionary stage dominating the neutrino emission).

Figures 1 and 4 suggest that the neutrino flux from stars in our Galaxy is about \(5\) orders of magnitude smaller than the solar neutrino flux and the diffuse supernova neutrino background (DSNB). Although the development of dedicated experimental techniques would be necessary to make the detection of stellar neutrinos possible, it is worth noting that the magnitude of the GS\(\nu\)F has never been larger than at present. The Earth formed when the Universe was about \(4.5\) Gyr old (after the first infall leading to star formation in the bulge and thick disk); star formation has taken place at a slightly decreasing rate since then. Hence, one can expect that the population of intermediate- and high-mass stars in the thin disk (which happens to be the most massive component of the Galaxy) has increased since then. Most of these are main-sequence stars, with neutrinos being abundantly produced.

Our forecast of the GS\(\nu\)F is in overall agreement with earlier work presented in Refs. [6], [9], both as for the spectral shape and normalization of the flux. However, our work is based on up-to-date 1D modelling of stellar evolution (with a complete nuclear reaction network) and accounts for a modeling of the Milky Way stellar population based on recent observational data. In principle, one could extend these results and estimate the diffuse cumulative flux from stars in other galaxies. Estimates of this sort can be found in Refs. [6], [9] based on simple assumptions.

5.5 Components of the GS\(\nu\)F↩︎

In order to compare the various components of the GS\(\nu\)F to the well-known solar neutrino flux, Table [tab:fluxes] summarizes the observed flux of neutrinos in all flavors from the thermonuclear reactions taking place in the Sun and the GS\(\nu\)F. We note that instead of relying on the solar neutrino flux from Ref. [10] (as otherwise done in this work), we quote the solar neutrino components from Ref. [88], obtained fitting all available solar neutrino data. For completeness, we also provide the \(Q\) value for each thermonuclear reaction channel, which corresponds to the maximum electron neutrino energy for the reaction channels leading to spectral energy distributions and the neutrino energy for the \(^7\)Be and pep channels responsible for spectral lines. We can see that the solar neutrino flux dominates over the GS\(\nu\)F across all neutrino-production channels. However, the ratio of the flux from pp neutrinos to CNO-cycle neutrinos is larger for the Sun than for the rest of the GS\(\nu\)F. This is because the neutrino emission from intermediate- and high-mass stars is dominated by the burning of hydrogen through the CNO cycle during their life as main sequence stars.

Figure 5 shows the rate per unit volume of pp, CNO, and \(^{18}\)F neutrinos as functions of the distance from Earth for the bulge, the thin disk, and the thick disk. For all Galactic components CNO neutrinos constitute the largest contribution to the overall neutrino flux. However, the radial evolution of each contribution depends on the matter distribution of each of the spatial components of our Milky-Way model. The largest neutrino flux at Earth should be expected from the Galactic Center, as highlighted by the vertical dashed line at \(8\) kpc. Of course, if the observer were to point in other directions, the GS\(\nu\)F should be expected to be significantly smaller.

Figure 6 represents the GS\(\nu\)F contributions of neutrinos and antineutrinos. The high-energy tail of the GS\(\nu\)F, for \(E_\nu \gtrsim 0.1\) MeV, comes from thermonuclear reactions and, therefore, is composed solely of neutrinos. However, the thermal component of the flux has equal contributions from neutrinos and antineutrinos since these are produced in pairs in thermal processes. For comparison, we also show other relevant neutrino and antineutrino backgrounds in the energy range of interest, such as the DSNB, solar neutrinos, geoneutrinos, and reactor antineutrinos [10]. We stress, however, that the reactor antineutrinos and geoneutrinos, although could be backgrounds to the GS\(\nu\)F detection, are strongly dependent on the neutrino telescope location; hence the fluxes show in Fig. 6 should only serve for orientation.

6 Concluding remarks↩︎

Stars produce neutrinos throughout their life because of thermal and thermonuclear processes. In addition to the routinely observed solar neutrinos, we should expect the GS\(\nu\)F between a few keV and tens of MeV from all stars in our Galaxy. In this work, we compute this neutrino signal. In order to do so, we consider three different components of our Galaxy: the bulge, a thin disk, and a thick disk. The SFH for each of these components is modelled according to the two-infall model. The baryon density distribution is inferred from the analysis presented in Ref. [27] which uses Gaia Data Release 2. In order to compute the neutrino energy spectrum expected from each star, we construct a suite of \(13\) progenitors with solar metallicity and ZAMS mass ranging between \(0.1\;M_\odot\) and \(40\;M_\odot\). We rely on the stellar evolution code MESA to evolve each of these stars from their pre-main sequence phase to their final fates. We assume that such models are representative of the Galactic population of stars with ZAMS in the range between \(0.08\;M_\odot\) and \(100\;M_\odot\). For each time snapshot during stellar evolution, MESA outputs the reaction rates and total energy released from the thermonuclear (pp chain, the CNO cycle, and \(^{18}\)F decay) and thermal processes (plasmon decay, photoneutrinos, pair annihilation, recombination, and bremsstrahlung) taking place in the core of each star. We then compute the neutrino energy spectrum for each of these processes during the evolution of the models in our suite.

We find that the bulk of the GS\(\nu\)F comes from thermal processes for neutrino energies \(\lesssim \mathcal{O}(0.1)\) MeV and thermonuclear ones otherwise. For energies in the range of \(0.1\, M_\odot\)–\(0.5\) MeV and above \(\sim1.6\) MeV, the GS\(\nu\)F is dominated by pp neutrinos from stars with ZAMS smaller than \(8\,M_\odot\); the remaining (intermediate) part of the spectrum for \(0.5 ~{\rm MeV} \lesssim E_\nu \lesssim 1.6\) MeV results from the CNO cycle occurring in stars with mass \(\gtrsim 1\; M_\odot\). The neutrino emission mainly comes from the thin disk of the Milky Way, since this is the youngest and most massive component of our Galaxy, with the largest density of stars with ZAMS mass \(\gtrsim 1\; M_\odot\) in their main-sequence evolutionary stage. The GS\(\nu\)F dominated by thermonuclear reactions is solely composed of neutrinos, whereas the thermal component of the GS\(\nu\)F has an equal number of neutrinos and antineutrinos (see Table [tab:reactions]). The bulk of the GS\(\nu\)F comes from stars within \(5\) and \(10\) kpc from Earth; this implies that we are favorably located to receive a copious amount of neutrinos from the Galactic Center where the density of stars is the highest.

The GS\(\nu\)F lies in the same energy region of solar neutrinos, reactor antineutrinos, geoneutrinos, the DSNB, and possibly the diffuse background of neutrinos from Thorne-Żytkow Objects [10], [89], although the spectral features of these fluxes are very different and the reactor and geoneutrino contributions are strongly dependent on the neutrino telescope location. In particular, the detection of the GS\(\nu\)F is challenged by the fact that this flux is about five orders of magnitude smaller than the one of solar neutrinos. It would be crucial to devote experimental efforts to overcome solar neutrinos and the DSNB, favoring a unique discovery for stellar archaeology. In fact, the detection of the GS\(\nu\)F would provide unique insight on the Galactic stellar population, especially for stars with ZAMS mass greater than \(1\;M_\odot\).

Leveraging on the directionality of the neutrino events detected at Earth could be crucial to discriminating between the GS\(\nu\)F and the solar and supernova backgrounds. In fact, we could combine our precise understanding of the solar neutrino flux with directional information on solar neutrino events to tag and remove the solar background. Moreover, the annual modulations of the solar neutrino flux could be employed as an additional discriminant factor to disentangle the solar neutrino background from the GS\(\nu\)F. We expect that a large fraction of GS\(\nu\)F comes from the Galactic Center. However, even with perfect directional reconstruction, the DSNB might only be reduced by approximately one order of magnitude (at most by a factor \(1/4 \pi\) due to the DSNB isotropy). Hence, work is strongly encouraged to improve the detection prospects of neutrinos with energy \(\lesssim 1.6\) MeV, where the DSNB is negligible.

As a final remark, we stress that a detailed knowledge of the GS\(\nu\)F can be of relevance to the discovery of physic scenarios beyond the Standard Model. For instance, in the presence of non-standard neutrino self-interactions, neutrinos from distant astrophysical sources would scatter off the GS\(\nu\)F, leading to potential observable features, such as spectral distortions or changes in the flavor composition. We refer the reader to Ref. [90]–and references therein–for analogous examples involving astrophysical neutrinos scattering off the cosmic neutrino background. Moreover, a high-energy neutrino flux could originate from the scattering of cosmic rays on the GS\(\nu\)F or the diffuse flux from stars in other galaxies, as proposed for cosmological relic neutrinos [91], [92]. If neutrinos were unstable, their emission from stars and decay en route to Earth could give rise to a detectable flux of photons or dark matter; this scenario was explored in Ref. [9] based on their prediction of the diffuse flux from all galaxies.

References↩︎