# Origin of $^{60}$Fe nuclei in cosmic rays: the contribution of local OB associations

Nicolas de Séréville, Vincent Tatischeff, Pierre Cristofari, Stefano Gabici, Roland Diehl

► To cite this version:

Nicolas de Séréville, Vincent Tatischeff, Pierre Cristofari, Stefano Gabici, Roland Diehl. Origin of $^{60}$Fe nuclei in cosmic rays: the contribution of local OB associations. Mon.Not.Roy.Astron.Soc., 2024, 530 (1), pp.684-698. 10.1093/mnras/stae336. hal-04446659

HAL Id: hal-04446659

https://hal.science/hal-04446659

Submitted on 16 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Origin of $^{60}$Fe nuclei in cosmic rays: the contribution of local OB associations

Nicolas de Séréville, Vincent Tatischeff, Pierre Cristofari, Stefano Gabici, Roland Diehl

► To cite this version:

Nicolas de Séréville, Vincent Tatischeff, Pierre Cristofari, Stefano Gabici, Roland Diehl. Origin of $^{60}$Fe nuclei in cosmic rays: the contribution of local OB associations. Monthly Notices of the Royal Astronomical Society, 2024, 530 (1), pp.684-698. 10.1093/mnras/stae336. hal-04446659

HAL Id: hal-04446659

https://hal.science/hal-04446659v1

Submitted on 16 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Origin of $^{60}$Fe nuclei in cosmic rays: the contribution of local OB associations

Nicolas de Séréville,<fnref n="1" />* Vincent Tatischeff,<fnref n="1" /> Pierre Cristofari,<fnref n="1" />,<fnref n="2" /> Stefano Gabici<fnref n="3" /> and Roland Diehl<fnref n="4" />,<fnref n="5" />

Accepted 2024 January 26. Received 2023 December 24; in original form 2023 July 27

## ABSTRACT

The presence of live $^{60}$Fe nuclei (lifetime of 3.8 Myr) in cosmic rays detected by the Advanced Composition Explorer/Cosmic Ray Isotope Spectrometer instrument suggests a nearby nucleosynthesis source. $^{60}$Fe is primarily produced in core-collapse supernovae, and we aim to clarify whether the detected $^{60}$Fe nuclei can be associated with a particular local supernova. We consider 25 OB associations and subgroups located within 1 kpc of the Solar system based on recent Gaia census. A model is developed that combines stellar population synthesis within these OB associations, cosmic ray acceleration within associated superbubbles, and cosmic ray transport to the Solar system. The most critical model parameter impacting $^{60}$Fe cosmic ray production is the explodability criterion, which determines if a massive star ends its life as a supernova. Our study points to the Scorpius–Centaurus (Sco–Cen) OB association as the most probable origin of the observed $^{60}$Fe nuclei, particularly suggesting they were accelerated in the Sco–Cen superbubble by a young supernova aged $\le 500$ kyr with a progenitor mass of approximately 13–20 $M_\odot$. A less likely source is the supernova at the origin of the Geminga pulsar 342 kyr ago, if the progenitor originated in the Orion OB1 association. The contribution of local OB associations to the cosmic ray density of stable $^{56}$Fe is estimated to be around 20 per cent, with some sensitivity to cosmic ray acceleration efficiency and diffusion coefficient. These findings shed light on the origins of cosmic ray nuclei, connecting them to nucleosynthesis events within our local cosmic neighbourhood.

**Key words:** nuclear reactions, nucleosynthesis, abundances – cosmic rays – acceleration of particles – open clusters and associations: general – pulsars: individual: Geminga software: simulations.

## 1 INTRODUCTION

Cosmic rays (CRs) are believed to be a common component of the interstellar medium (ISM) in galaxies, with an energy density that is comparable to the energy densities of other ISM components, such as the kinetic energy of bulk atomic or molecular gas motions, the thermal energy of hot plasma, and the magnetic energy of regular and turbulent fields (Blasi 2013; Gabici et al. 2019). The consensual picture is that strong shock waves (of Mach number $\gg 1$) accelerate CRs through diffusive shock acceleration (DSA; Blandford & Ostriker 1978; Berezhko & Ellison 1999; Lee, Ellison & Nagataki 2012). Such shock waves have typically been associated with massive star winds, supernovae (SNe), and/or their remnants (Drury 2012; Blasi 2013). However, several questions about CR acceleration remain poorly understood, including the spectra and maximum energy achieved by the CR particles at their sources and the efficiency of the acceleration process (Gabici et al. 2019).

The propagation of CRs in the Galaxy after escaping from their sources also remains an important topic of research. The mean CR

lifetime in the Milky Way ($\tau_{CR} \sim 15$ Myr) is much longer than the light-crossing time ($< 0.1$ Myr), which is explained by diffusive confinement of the non-thermal particles by scattering on small-scale electromagnetic turbulence. Both pre-existing magnetohydrodynamic (MHD) turbulence (Lazarian & Xu 2021; Lazarian, Xu & Hu 2023) and plasma waves self-generated by the CR streaming instability (Kulsrud & Pearce 1969; Farmer & Goldreich 2004) are considered as scattering centres, but their relative importance for CR transport strongly depends on local plasma conditions in the multiphase ISM (Kempski & Quataert 2022). Effective diffusion models are commonly used to describe CR propagation (e.g. Evoli, Aloisio & Blasi 2019), but the diffusion coefficient is hard to determine from first principles and may significantly vary within the Galaxy.

Recent gamma-ray observations of CR interactions with interstellar matter report significant variations of CR densities in specific regions, such as the Central Molecular Zone (HESS Collaboration 2016), the inner Galaxy region between 1.5 and 4.5 kpc from the Galactic Centre (Peron et al. 2021), and the Cygnus region at a distance of 2–3 kpc from Earth (Ackermann et al. 2011; Astiasarain et al. 2023, see also discussions in Aharonian, Yang & de Oña

<footnote n="1">¹Université Paris-Saclay, CNRS/IN2P3, IJCLab, F-91405 Orsay, France</footnote>

<footnote n="2">²Laboratoire Univers et Théories, Université de Paris, Observatoire de Paris, Université PSL, CNRS, F-92190 Meudon, France</footnote>

<footnote n="3">³Université de Paris, CNRS, Astroparticule et Cosmologie, F-75006 Paris, France</footnote>

<footnote n="4">⁴Max Planck Institut für extraterrestrische Physik, Giessenbachstraße 1, D-85748 Garching, Germany</footnote>

<footnote n="5">⁵Excellence Cluster Origins, Boltzmannstraße 2, D-85748 Garching, Germany</footnote>

<footnote>* E-mail: nicolas.de-sereville@ijclab.in2p3.fr</footnote>
# Origin of $^{60}$Fe nuclei in cosmic rays: the contribution of local OB associations

Nicolas de Séréville,<fnref n="1" />* Vincent Tatischeff,<fnref n="1" /> Pierre Cristofari,<fnref n="1" />,<fnref n="2" /> Stefano Gabici<fnref n="3" /> and Roland Diehl<fnref n="4" />,<fnref n="5" />

Accepted 2024 January 26. Received 2023 December 24; in original form 2023 July 27

## ABSTRACT

The presence of live $^{60}$Fe nuclei (lifetime of 3.8 Myr) in cosmic rays detected by the *Advanced Composition Explorer/Cosmic Ray Isotope Spectrometer* instrument suggests a nearby nucleosynthesis source. $^{60}$Fe is primarily produced in core-collapse supernovae, and we aim to clarify whether the detected $^{60}$Fe nuclei can be associated with a particular local supernova. We consider 25 OB associations and subgroups located within 1 kpc of the Solar system based on recent *Gaia* census. A model is developed that combines stellar population synthesis within these OB associations, cosmic ray acceleration within associated superbubbles, and cosmic ray transport to the Solar system. The most critical model parameter impacting $^{60}$Fe cosmic ray production is the explodability criterion, which determines if a massive star ends its life as a supernova. Our study points to the Scorpius–Centaurus (Sco–Cen) OB association as the most probable origin of the observed $^{60}$Fe nuclei, particularly suggesting they were accelerated in the Sco–Cen superbubble by a young supernova aged $\le 500$ kyr with a progenitor mass of approximately 13–20 $M_{\odot}$. A less likely source is the supernova at the origin of the Geminga pulsar 342 kyr ago, if the progenitor originated in the Orion OB1 association. The contribution of local OB associations to the cosmic ray density of stable $^{56}$Fe is estimated to be around 20 per cent, with some sensitivity to cosmic ray acceleration efficiency and diffusion coefficient. These findings shed light on the origins of cosmic ray nuclei, connecting them to nucleosynthesis events within our local cosmic neighbourhood.

**Key words:** nuclear reactions, nucleosynthesis, abundances – cosmic rays – acceleration of particles – open clusters and associations: general – pulsars: individual: Geminga software: simulations.

## 1 INTRODUCTION

Cosmic rays (CRs) are believed to be a common component of the interstellar medium (ISM) in galaxies, with an energy density that is comparable to the energy densities of other ISM components, such as the kinetic energy of bulk atomic or molecular gas motions, the thermal energy of hot plasma, and the magnetic energy of regular and turbulent fields (Blasi 2013; Gabici et al. 2019). The consensual picture is that strong shock waves (of Mach number $\gg 1$) accelerate CRs through diffusive shock acceleration (DSA; Blandford & Ostriker 1978; Berezhko & Ellison 1999; Lee, Ellison & Nagataki 2012). Such shock waves have typically been associated with massive star winds, supernovae (SNe), and/or their remnants (Drury 2012; Blasi 2013). However, several questions about CR acceleration remain poorly understood, including the spectra and maximum energy achieved by the CR particles at their sources and the efficiency of the acceleration process (Gabici et al. 2019).

The propagation of CRs in the Galaxy after escaping from their sources also remains an important topic of research. The mean CR

lifetime in the Milky Way ($\tau_{CR} \sim 15$ Myr) is much longer than the light-crossing time ($< 0.1$ Myr), which is explained by diffusive confinement of the non-thermal particles by scattering on small-scale electromagnetic turbulence. Both pre-existing magnetohydrodynamic (MHD) turbulence (Lazarian & Xu 2021; Lazarian, Xu & Hu 2023) and plasma waves self-generated by the CR streaming instability (Kulsrud & Pearce 1969; Farmer & Goldreich 2004) are considered as scattering centres, but their relative importance for CR transport strongly depends on local plasma conditions in the multiphase ISM (Kempski & Quataert 2022). Effective diffusion models are commonly used to describe CR propagation (e.g. Evoli, Aloisio & Blasi 2019), but the diffusion coefficient is hard to determine from first principles and may significantly vary within the Galaxy.

Recent gamma-ray observations of CR interactions with interstellar matter report significant variations of CR densities in specific regions, such as the Central Molecular Zone (HESS Collaboration 2016), the inner Galaxy region between 1.5 and 4.5 kpc from the Galactic Centre (Peron et al. 2021), and the Cygnus region at a distance of 2–3 kpc from Earth (Ackermann et al. 2011; Astiasarain et al. 2023, see also discussions in Aharonian, Yang & de Oña

<footnote n="1">¹Université Paris-Saclay, CNRS/IN2P3, IJCLab, F-91405 Orsay, France</footnote>

<footnote n="2">²Laboratoire Univers et Théories, Université de Paris, Observatoire de Paris, Université PSL, CNRS, F-92190 Meudon, France</footnote>

<footnote n="3">³Université de Paris, CNRS, Astroparticule et Cosmologie, F-75006 Paris, France</footnote>

<footnote n="4">⁴Max Planck Institut für extraterrestrische Physik, Giessenbachstraße 1, D-85748 Garching, Germany</footnote>

<footnote n="5">⁵Excellence Cluster Origins, Boltzmannstraße 2, D-85748 Garching, Germany</footnote>

<footnote>* E-mail: nicolas.de-sereville@ijclab.in2p3.fr</footnote>
Wilhelmi 2019; Bykov & Kalyashova 2022). Significant variations of the measured CR-induced ionization rate in molecular clouds also point to variations in the density of low-energy CRs throughout the Galaxy (Indriolo & McCall 2012; Gabici 2022; Phan et al. 2023). In particular, the local spectrum of MeV CRs measured by the Voyager probes may not be representative of the low-energy CR spectrum elsewhere in the Galaxy (Phan et al. 2021). In addition, according to Kachelrieß, Neronov & Semikoz (2018), the unexpected hardness of the CR positron and antiproton spectra above ~100 GeV can be explained by a significant contribution to the CR flux of particles accelerated in a local SN some 2–3 Myr ago.

The detection of $^{60}$Fe nuclei in CRs with the Cosmic Ray Isotope Spectrometer (CRIS) on the Advanced Composition Explorer (ACE) spacecraft (Binns et al. 2016) offers a unique opportunity to study the contribution of localized and nearby sources to the CR population seen here, hence addressing CR source and transport simultaneously. $^{60}$Fe is a primary CR, i.e. it is not produced to any significant extent by nuclear spallation of heavier CRs in the ISM. It is thought to be synthesized mainly in core-collapse supernovae (ccSNe) of massive stars. Its radioactive lifetime of 3.8 Myr is sufficiently long such that it can potentially survive the time interval between nucleosynthesis and detection on Earth. But the $^{60}$Fe lifetime is significantly shorter than $\tau_{CR}$, which suggests that nucleosynthesis sites far out in the Galaxy are plausibly beyond reach for $^{60}$Fe CRs surviving such a journey.

$^{60}$Fe has also been found in sediments from the Pacific ocean floor (Knie et al. 2004), complemented by findings in other sediments across Earth and even on the Moon (Wallner et al. 2016, 2021). Its live presence on Earth, combined with its radioactive decay time, and with typical velocities for the transport of interstellar matter (transport of $^{60}$Fe to Earth generally assumes adsorption on dust grains travelling at velocities of the order of ~10 km s⁻¹), suggested that it may be due to recent nucleosynthesis activity near the Solar system.

In parallel to the CR measurements, and to the recent data obtained on $^{60}$Fe in sediments and on the Moon, our knowledge of the distribution of stars, and especially massive stars and OB associations in our local environment within a few kiloparsec, is rapidly increasing, as recently illustrated with Gaia observations (Zucker et al. 2022, 2023). In the problem of the origin of CRs, OB associations are especially relevant, since they are expected to substantially enrich the ISM, injecting nuclear material through their winds and when exploding. The potential important contribution of OB associations in the CR content has been discussed in several works (Parizot et al. 2004; Binns et al. 2007; Murphy et al. 2016; Tatischeff et al. 2021). The recent $^{60}$Fe data and ever-increasing knowledge on the local OB associations provide an opportunity for probing the contribution of OB associations to CRs.

In this paper, we aim to set up a bottom-up model for the origin of $^{60}$Fe in CRs near Earth, based on modelling both the plausible nearby massive star groups as sources of the nucleosynthesis ejecta including $^{60}$Fe, together with modelling the acceleration near the sources, and the transport through the specifics of ISM trajectories from the sources to near-Earth space. We rely on Monte Carlo simulations, developing a model combining a description of the OB stellar population, accounting for CR acceleration and transport, and confront it to available $^{60}$Fe data. The model also allows us to discuss the origin of other CR nuclei such as $^{56}$Fe and $^{26}$Al.

This paper is organized as follows. First, we convert the measurement data of $^{60}$Fe in CRs into interstellar fluxes (Section 2). Then, we present our population synthesis model for the determination of

time-dependent production of $^{60}$Fe, followed by CR acceleration and transport (Section 3). We apply this to nearby massive star groups (Section 4), and evaluate these results towards constraints for locally found $^{60}$Fe CRs (Section 5). We conclude with a discussion of the sensitivity of our findings to various assumptions and ingredients of this bottom-up modelling.

## 2 DENSITY OF $^{60}$FE AND $^{26}$AL CRS FROM ACE/CRIS MEASUREMENTS

ACE/CRIS collected $^{56}$Fe and $^{60}$Fe CR nuclei between ~195 and ~500 MeV nucleon⁻¹, reporting 15 $^{60}$Fe CR nuclei (Binns et al. 2016). The reconstructed mean energy at the top of the CRIS instrument is 340 MeV nucleon⁻¹ for $^{56}$Fe and 327 MeV nucleon⁻¹ for $^{60}$Fe. According to Binns et al. (2016), the CR modulation inside the Solar system during the 17 yr period of the data taking can be accounted for with an average force-field potential $\phi = 453$ MV, corresponding to an energy loss of 210 MeV nucleon⁻¹ for $^{56}$Fe and 196 MeV nucleon⁻¹ for $^{60}$Fe. Thus, the mean energies in the local interstellar space are 550 MeV nucleon⁻¹ for $^{56}$Fe and 523 MeV nucleon⁻¹ for $^{60}$Fe, and the corresponding velocities are 0.778c and 0.768c (c is the speed of light).

The measured iron isotopic ratio near Earth is $(^{60}\text{Fe}/^{56}\text{Fe})_{\text{CRIS}} = (4.6 \pm 1.7) \times 10^{-5}$ (Binns et al. 2016). The flux ratio in the local ISM (LISM) can be estimated from the force-field approximation to the transport equation describing the CR modulation in the heliosphere (Gleeson & Axford 1968). In this simple model, the CR flux in the LISM is related to the one measured near Earth by a shift in particle momentum, which gives for the Fe isotopic ratio:

$$(^{60}\text{Fe}/^{56}\text{Fe})_{\text{LISM}} = (^{60}\text{Fe}/^{56}\text{Fe})_{\text{CRIS}} \times (p_{60,\text{LISM}}/p_{60,\text{CRIS}})^2 \times (p_{56,\text{CRIS}}/p_{56,\text{LISM}})^2, \quad (1)$$

where $p_{56,CRIS} = 48.5$ GeV/c, $p_{60,CRIS} = 50.8$ GeV/c, $p_{56,LISM} = 64.5$ GeV/c, and $p_{60,LISM} = 67.0$ GeV/c are the $^{56}$Fe and $^{60}$Fe mean momenta at the top of the CRIS instrument and in the LISM. Thus, $(^{60}\text{Fe}/^{56}\text{Fe})_{\text{LISM}} = (4.5 \pm 1.7) \times 10^{-5}$.

The spectrum of $^{56}$Fe CRs in the LISM can be estimated from the work of Boschini et al. (2021), who used recent the Alpha Magnetic Spectrometer 02 (AMS-02) results (Aguilar et al. 2021), together with Voyager 1 and ACE/CRIS data, to study the origin of Fe in the CR population. Their calculations are based on the GALPROP code to model the CR propagation in the ISM (Strong & Moskalenko 1998) and the HELMOD model to describe the particle transport within the heliosphere (Boschini et al. 2019). Integrating the iron spectrum given by these authors in the energy range from 400 to 700 MeV nucleon⁻¹, which approximately corresponds to the range of the CRIS measurements, we find $I_{\text{LISM}}(\text{Fe}) = 2.9 \times 10^{-5} \text{ cm}^{-2} \text{ s}^{-1} \text{ sr}^{-1}$ and the density $n_{\text{LISM}}(\text{Fe}) = I_{\text{LISM}} \times 4\pi/v = 1.6 \times 10^{-14} \text{ CR cm}^{-3}$. The intensity of $^{60}$Fe CRs between 400 and 700 MeV nucleon⁻¹ in the LISM is $I_{\text{LISM}}(^{60}\text{Fe}) = I_{\text{LISM}}(\text{Fe}) \times (^{60}\text{Fe}/^{56}\text{Fe})_{\text{LISM}} = (1.3 \pm 0.5) \times 10^{-9} \text{ cm}^{-2} \text{ s}^{-1} \text{ sr}^{-1}$ and the density in this energy range is $n_{\text{LISM}}(^{60}\text{Fe}) = (7.1 \pm 2.7) \times 10^{-19} \text{ CR cm}^{-3}$.

Recently, Boschini et al. (2022) found that the aluminium CR spectrum measured by AMS-02 presents a significant excess in the rigidity range from 2 to 7 GV compared to the spectrum predicted with the GALPROP-HELMOD framework from spallation of $^{28}$Si CRs and heavier nuclei. They suggested that this excess could be attributed to a source of primary CRs of radioactive $^{26}$Al (half-life $T_{1/2} = 7.17 \times 10^5$ yr) possibly related to the well-known $^{22}$Ne excess in the CR composition. The latter is interpreted as arising from the acceleration of massive star wind material in OB associations (see
Wilhelmi 2019; Bykov & Kalyashova 2022). Significant variations of the measured CR-induced ionization rate in molecular clouds also point to variations in the density of low-energy CRs throughout the Galaxy (Indriolo & McCall 2012; Gabici 2022; Phan et al. 2023). In particular, the local spectrum of MeV CRs measured by the Voyager probes may not be representative of the low-energy CR spectrum elsewhere in the Galaxy (Phan et al. 2021). In addition, according to Kachelrieß, Neronov & Semikoz (2018), the unexpected hardness of the CR positron and antiproton spectra above ~100 GeV can be explained by a significant contribution to the CR flux of particles accelerated in a local SN some 2–3 Myr ago.

The detection of $^{60}$Fe nuclei in CRs with the Cosmic Ray Isotope Spectrometer (CRIS) on the Advanced Composition Explorer (ACE) spacecraft (Binns et al. 2016) offers a unique opportunity to study the contribution of localized and nearby sources to the CR population seen here, hence addressing CR source and transport simultaneously. $^{60}$Fe is a primary CR, i.e. it is not produced to any significant extent by nuclear spallation of heavier CRs in the ISM. It is thought to be synthesized mainly in core-collapse supernovae (ccSNe) of massive stars. Its radioactive lifetime of 3.8 Myr is sufficiently long such that it can potentially survive the time interval between nucleosynthesis and detection on Earth. But the $^{60}$Fe lifetime is significantly shorter than $\tau_{CR}$, which suggests that nucleosynthesis sites far out in the Galaxy are plausibly beyond reach for $^{60}$Fe CRs surviving such a journey.

$^{60}$Fe has also been found in sediments from the Pacific ocean floor (Knie et al. 2004), complemented by findings in other sediments across Earth and even on the Moon (Wallner et al. 2016, 2021). Its live presence on Earth, combined with its radioactive decay time, and with typical velocities for the transport of interstellar matter (transport of $^{60}$Fe to Earth generally assumes adsorption on dust grains travelling at velocities of the order of ~10 km s⁻¹), suggested that it may be due to recent nucleosynthesis activity near the Solar system.

In parallel to the CR measurements, and to the recent data obtained on $^{60}$Fe in sediments and on the Moon, our knowledge of the distribution of stars, and especially massive stars and OB associations in our local environment within a few kiloparsec, is rapidly increasing, as recently illustrated with Gaia observations (Zucker et al. 2022, 2023). In the problem of the origin of CRs, OB associations are especially relevant, since they are expected to substantially enrich the ISM, injecting nuclear material through their winds and when exploding. The potential important contribution of OB associations in the CR content has been discussed in several works (Parizot et al. 2004; Binns et al. 2007; Murphy et al. 2016; Tatischeff et al. 2021). The recent $^{60}$Fe data and ever-increasing knowledge on the local OB associations provide an opportunity for probing the contribution of OB associations to CRs.

In this paper, we aim to set up a bottom-up model for the origin of $^{60}$Fe in CRs near Earth, based on modelling both the plausible nearby massive star groups as sources of the nucleosynthesis ejecta including $^{60}$Fe, together with modelling the acceleration near the sources, and the transport through the specifics of ISM trajectories from the sources to near-Earth space. We rely on Monte Carlo simulations, developing a model combining a description of the OB stellar population, accounting for CR acceleration and transport, and confront it to available $^{60}$Fe data. The model also allows us to discuss the origin of other CR nuclei such as $^{56}$Fe and $^{26}$Al.

This paper is organized as follows. First, we convert the measurement data of $^{60}$Fe in CRs into interstellar fluxes (Section 2). Then, we present our population synthesis model for the determination of

time-dependent production of $^{60}$Fe, followed by CR acceleration and transport (Section 3). We apply this to nearby massive star groups (Section 4), and evaluate these results towards constraints for locally found $^{60}$Fe CRs (Section 5). We conclude with a discussion of the sensitivity of our findings to various assumptions and ingredients of this bottom-up modelling.

## 2 DENSITY OF $^{60}$FE AND $^{26}$AL CRS FROM ACE/CRIS MEASUREMENTS

ACE/CRIS collected $^{56}$Fe and $^{60}$Fe CR nuclei between ~195 and ~500 MeV nucleon⁻¹, reporting 15 $^{60}$Fe CR nuclei (Binns et al. 2016). The reconstructed mean energy at the top of the CRIS instrument is 340 MeV nucleon⁻¹ for $^{56}$Fe and 327 MeV nucleon⁻¹ for $^{60}$Fe. According to Binns et al. (2016), the CR modulation inside the Solar system during the 17 yr period of the data taking can be accounted for with an average force-field potential $\phi = 453$ MV, corresponding to an energy loss of 210 MeV nucleon⁻¹ for $^{56}$Fe and 196 MeV nucleon⁻¹ for $^{60}$Fe. Thus, the mean energies in the local interstellar space are 550 MeV nucleon⁻¹ for $^{56}$Fe and 523 MeV nucleon⁻¹ for $^{60}$Fe, and the corresponding velocities are 0.778c and 0.768c (c is the speed of light).

The measured iron isotopic ratio near Earth is $(^{60}\text{Fe}/^{56}\text{Fe})_{\text{CRIS}} = (4.6 \pm 1.7) \times 10^{-5}$ (Binns et al. 2016). The flux ratio in the local ISM (LISM) can be estimated from the force-field approximation to the transport equation describing the CR modulation in the heliosphere (Gleeson & Axford 1968). In this simple model, the CR flux in the LISM is related to the one measured near Earth by a shift in particle momentum, which gives for the Fe isotopic ratio:

$$(^{60}\text{Fe}/^{56}\text{Fe})_{\text{LISM}} = (^{60}\text{Fe}/^{56}\text{Fe})_{\text{CRIS}} \times (p_{60,\text{LISM}}/p_{60,\text{CRIS}})^2 \times (p_{56,\text{CRIS}}/p_{56,\text{LISM}})^2, \quad (1)$$

where $p_{56,\text{CRIS}} = 48.5$ GeV/c, $p_{60,\text{CRIS}} = 50.8$ GeV/c, $p_{56,\text{LISM}} = 64.5$ GeV/c, and $p_{60,\text{LISM}} = 67.0$ GeV/c are the $^{56}$Fe and $^{60}$Fe mean momenta at the top of the CRIS instrument and in the LISM. Thus, $(^{60}\text{Fe}/^{56}\text{Fe})_{\text{LISM}} = (4.5 \pm 1.7) \times 10^{-5}$.

The spectrum of $^{56}$Fe CRs in the LISM can be estimated from the work of Boschini et al. (2021), who used recent the Alpha Magnetic Spectrometer 02 (AMS-02) results (Aguilar et al. 2021), together with Voyager 1 and ACE/CRIS data, to study the origin of Fe in the CR population. Their calculations are based on the GALPROP code to model the CR propagation in the ISM (Strong & Moskalenko 1998) and the HELMOD model to describe the particle transport within the heliosphere (Boschini et al. 2019). Integrating the iron spectrum given by these authors in the energy range from 400 to 700 MeV nucleon⁻¹, which approximately corresponds to the range of the CRIS measurements, we find $I_{\text{LISM}}(\text{Fe}) = 2.9 \times 10^{-5} \text{ cm}^{-2} \text{ s}^{-1} \text{ sr}^{-1}$ and the density $n_{\text{LISM}}(\text{Fe}) = I_{\text{LISM}} \times 4\pi/v = 1.6 \times 10^{-14} \text{ CR cm}^{-3}$. The intensity of $^{60}$Fe CRs between 400 and 700 MeV nucleon⁻¹ in the LISM is $I_{\text{LISM}}(^{60}\text{Fe}) = I_{\text{LISM}}(\text{Fe}) \times (^{60}\text{Fe}/^{56}\text{Fe})_{\text{LISM}} = (1.3 \pm 0.5) \times 10^{-9} \text{ cm}^{-2} \text{ s}^{-1} \text{ sr}^{-1}$ and the density in this energy range is $n_{\text{LISM}}(^{60}\text{Fe}) = (7.1 \pm 2.7) \times 10^{-19} \text{ CR cm}^{-3}$.

Recently, Boschini et al. (2022) found that the aluminium CR spectrum measured by AMS-02 presents a significant excess in the rigidity range from 2 to 7 GV compared to the spectrum predicted with the GALPROP-HELMOD framework from spallation of $^{28}$Si CRs and heavier nuclei. They suggested that this excess could be attributed to a source of primary CRs of radioactive $^{26}$Al (half-life $T_{1/2} = 7.17 \times 10^5$ yr) possibly related to the well-known $^{22}$Ne excess in the CR composition. The latter is interpreted as arising from the acceleration of massive star wind material in OB associations (see
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Illustration of our bottom-up CR source model. The population synthesis of a massive star cluster (left) evaluates stellar evolution and SN models to determine the abundance of particular isotopes within the cluster gas and its superbubble. This gas is accelerated to CRs, which are transported through galactic ISM towards the Sun (right).</figcaption></figure>

Tatischeff et al. 2021, and references therein). Here, we study the contribution of primary $^{26}$Al CRs originating together with $^{60}$Fe from local OB associations.

ACE/CRIS measured ($^{26}$Al/$^{27}$Al)CRIS = 0.042 ± 0.002 between 125 and 300 MeV nucleon$^{-1}$, corresponding to the LISM energy range 190–560 MeV nucleon$^{-1}$ (Yanasak et al. 2001). From the Al spectrum in the LISM computed by Boschini et al. (2022), we find the mean energy of Al CRs in the LISM to be 355 MeV nucleon$^{-1}$ ($v$ = 0.690c) and the LISM density of Al CRs between 190 and 560 MeV nucleon$^{-1}$ to be $n_{\text{LISM}}(\text{Al}) = 1.0 \times 10^{-14}$ CR cm$^{-3}$. The $^{26}$Al CRs density in this energy range is then $n_{\text{LISM}}(^{26}\text{Al}) = (4.2 \pm 0.2) \times 10^{-16}$ CR cm$^{-3}$.

## 3 CR POPULATION SYNTHESIS AND TRANSPORT

We developed a bottom-up model (Fig. 1) for the CR flux at the Solar system, integrating contributions from the presumed sources of radioactive nuclei within massive star clusters. The basic ingredients are the yields of ejecta from stars and SNe. For each cluster, its age and richness are used together with a generic initial mass distribution to determine proper weighting, thus building a time profile of interstellar nuclide abundances for and within each specific cluster. With plausible assumptions about CR acceleration efficiency within such a massive star group and the likely superbubble configuration resulting from the clustered stellar and SN activity, we derive a CR source density for each star cluster, as it varies with time. Propagation of these CRs towards the Solar system requires a CR transport model that accounts for the location of the source within the Galaxy and its distance from the Solar system, accounting for specifics of CR transport in the solar neighbourhood. Integrating contributions of all sources from which CRs could have reached instruments near Earth in our present epoch, we thus obtain a bottom-up determination of the local CR flux in terms of model parameters based on stars and SNe.

Our model is similar to and builds on those of Gounelle et al. (2009), Voss et al. (2009), and Young (2014) for example, and we focus on the activity of massive stars ($M \ge 8$ M$_{\odot}$) in OB associations. The novelty of our work is to couple the nucleosynthesis output of a massive star group to a CR transport model, which then allows the prediction of the flux of CRs near Earth. Adjusting the parameters of our model to best match CR data

taken near Earth, we can therefore constrain the origins of locally observed CRs, backtracing them to the contributing massive star groups.

### 3.1 Radioisotope production at CR sources

We aim to know the production of radioactive isotopes from the ensemble of stars and SNe in the nearby Galaxy. Population synthesis is the tool commonly used to predict the integrated outcome and properties of stellar populations (Cerviño & Luridiana 2006; Cerviño 2013). This approach has been used in particular to predict the radionuclide enrichment of the ISM near OB associations using such a bottom-up approach that implements our knowledge about stars, their evolution, and their nucleosynthesis yields (Voss et al. 2009). In the following, we describe key aspects of the stellar population synthesis part of our model.

#### 3.1.1 Population synthesis ingredients

Initial mass function (IMF). A population of stars that formed simultaneously and within the same environment, such as in a cluster, is characterized by the distribution in mass of the stars after having been formed, the IMF. Observationally, the stellar population seen within a cluster reflects the current mass distribution. From this, one may estimate an initial mass distribution by corrections for the stars of high mass that already may have disappeared, when the cluster age is known, or can reliably be estimated. There is considerable debate about how generic the initial mass distribution may be, or how it may depend on the feedback properties for different stellar density and interstellar gas density (e.g. Kroupa 2020). However, the widely observed similarity of the power-law shape of the mass distribution (Kroupa 2001) suggests that the mass distribution of newly formed stars is a result of the physical processes during star formation, as it may be inhibited or modified by energetic feedback from the newly formed stars. The IMF was initially described for intermediate to large stellar masses by a single power-law function by Salpeter (1955). Toward the low-mass end down to the brown dwarf limit the IMF flattens and can be described by a lognormal shape (Miller & Scalo 1979), or a broken power law (Kroupa 2001). Our model implements any IMF described by a multipart power law, and we use as default the parameters given by Kroupa (2001, equation 6); this gives an average mass of the association members of 0.21 M$_{\odot}$ and a fraction of stars having a mass greater than 8 M$_{\odot}$ of $1.6 \times 10^{-3}$. The stellar content of specific known OB associations is mostly derived from a census of bright stars such as O and B stars ($M \ge 2.8$ M$_{\odot}$; Habets & Heintze 1981), thus only the high-mass end of the IMF is relevant. The upper end of the mass distribution for massive stars is debated (e.g. Heger et al. 2003; Vanbeveren 2009; Schneider et al. 2018). Theoretical uncertainties derive from the star formation processes for very massive stars as nuclear burning sets in during the mass accretion phase, but also from the late evolution of massive stars towards core collapse that may be inhibited by pair instability. Observationally, stars with masses up to 300 M$_{\odot}$ have been claimed to exist in the Large Magellanic Cloud's 30 Dor region (Schneider et al. 2018). In our model, we consider an upper limit of 120 M$_{\odot}$, which seems reasonable compared to the observational upper limit for single stars in our Galaxy of ~150 M$_{\odot}$ (Maíz Apellániz et al. 2007). This allows us to use the full range of the mass grid for stellar yields from Limongi & Chieffi (2018).

Stellar yields. Massive stars contribute significantly to the enrichment of the ISM by releasing nuclear-processed material through
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Illustration of our bottom-up CR source model. The population synthesis of a massive star cluster (left) evaluates stellar evolution and SN models to determine the abundance of particular isotopes within the cluster gas and its superbubble. This gas is accelerated to CRs, which are transported through galactic ISM towards the Sun (right).</figcaption></figure>

Tatischeff et al. 2021, and references therein). Here, we study the contribution of primary $^{26}$Al CRs originating together with $^{60}$Fe from local OB associations.

ACE/CRIS measured ($^{26}$Al/$^{27}$Al)CRIS = 0.042 ± 0.002 between 125 and 300 MeV nucleon$^{-1}$, corresponding to the LISM energy range 190–560 MeV nucleon$^{-1}$ (Yanasak et al. 2001). From the Al spectrum in the LISM computed by Boschini et al. (2022), we find the mean energy of Al CRs in the LISM to be 355 MeV nucleon$^{-1}$ ($v$ = 0.690c) and the LISM density of Al CRs between 190 and 560 MeV nucleon$^{-1}$ to be $n_{\text{LISM}}(\text{Al}) = 1.0 \times 10^{-14}$ CR cm$^{-3}$. The $^{26}$Al CRs density in this energy range is then $n_{\text{LISM}}(^{26}\text{Al}) = (4.2 \pm 0.2) \times 10^{-16}$ CR cm$^{-3}$.

## 3 CR POPULATION SYNTHESIS AND TRANSPORT

We developed a bottom-up model (Fig. 1) for the CR flux at the Solar system, integrating contributions from the presumed sources of radioactive nuclei within massive star clusters. The basic ingredients are the yields of ejecta from stars and SNe. For each cluster, its age and richness are used together with a generic initial mass distribution to determine proper weighting, thus building a time profile of interstellar nuclide abundances for and within each specific cluster. With plausible assumptions about CR acceleration efficiency within such a massive star group and the likely superbubble configuration resulting from the clustered stellar and SN activity, we derive a CR source density for each star cluster, as it varies with time. Propagation of these CRs towards the Solar system requires a CR transport model that accounts for the location of the source within the Galaxy and its distance from the Solar system, accounting for specifics of CR transport in the solar neighbourhood. Integrating contributions of all sources from which CRs could have reached instruments near Earth in our present epoch, we thus obtain a bottom-up determination of the local CR flux in terms of model parameters based on stars and SNe.

Our model is similar to and builds on those of Gounelle et al. (2009), Voss et al. (2009), and Young (2014) for example, and we focus on the activity of massive stars ($M \ge 8$ M$_{\odot}$) in OB associations. The novelty of our work is to couple the nucleosynthesis output of a massive star group to a CR transport model, which then allows the prediction of the flux of CRs near Earth. Adjusting the parameters of our model to best match CR data

taken near Earth, we can therefore constrain the origins of locally observed CRs, backtracing them to the contributing massive star groups.

### 3.1 Radioisotope production at CR sources

We aim to know the production of radioactive isotopes from the ensemble of stars and SNe in the nearby Galaxy. Population synthesis is the tool commonly used to predict the integrated outcome and properties of stellar populations (Cerviño & Luridiana 2006; Cerviño 2013). This approach has been used in particular to predict the radionuclide enrichment of the ISM near OB associations using such a bottom-up approach that implements our knowledge about stars, their evolution, and their nucleosynthesis yields (Voss et al. 2009). In the following, we describe key aspects of the stellar population synthesis part of our model.

#### 3.1.1 Population synthesis ingredients

Initial mass function (IMF). A population of stars that formed simultaneously and within the same environment, such as in a cluster, is characterized by the distribution in mass of the stars after having been formed, the IMF. Observationally, the stellar population seen within a cluster reflects the current mass distribution. From this, one may estimate an initial mass distribution by corrections for the stars of high mass that already may have disappeared, when the cluster age is known, or can reliably be estimated. There is considerable debate about how generic the initial mass distribution may be, or how it may depend on the feedback properties for different stellar density and interstellar gas density (e.g. Kroupa 2020). However, the widely observed similarity of the power-law shape of the mass distribution (Kroupa 2001) suggests that the mass distribution of newly formed stars is a result of the physical processes during star formation, as it may be inhibited or modified by energetic feedback from the newly formed stars. The IMF was initially described for intermediate to large stellar masses by a single power-law function by Salpeter (1955). Toward the low-mass end down to the brown dwarf limit the IMF flattens and can be described by a lognormal shape (Miller & Scalo 1979), or a broken power law (Kroupa 2001). Our model implements any IMF described by a multipart power law, and we use as default the parameters given by Kroupa (2001, equation 6); this gives an average mass of the association members of 0.21 M$_{\odot}$ and a fraction of stars having a mass greater than 8 M$_{\odot}$ of $1.6 \times 10^{-3}$. The stellar content of specific known OB associations is mostly derived from a census of bright stars such as O and B stars ($M \ge 2.8$ M$_{\odot}$; Habets & Heintze 1981), thus only the high-mass end of the IMF is relevant. The upper end of the mass distribution for massive stars is debated (e.g. Heger et al. 2003; Vanbeveren 2009; Schneider et al. 2018). Theoretical uncertainties derive from the star formation processes for very massive stars as nuclear burning sets in during the mass accretion phase, but also from the late evolution of massive stars towards core collapse that may be inhibited by pair instability. Observationally, stars with masses up to 300 M$_{\odot}$ have been claimed to exist in the Large Magellanic Cloud's 30 Dor region (Schneider et al. 2018). In our model, we consider an upper limit of 120 M$_{\odot}$, which seems reasonable compared to the observational upper limit for single stars in our Galaxy of ~150 M$_{\odot}$ (Maíz Apellániz et al. 2007). This allows us to use the full range of the mass grid for stellar yields from Limongi & Chieffi (2018).

Stellar yields. Massive stars contribute significantly to the enrichment of the ISM by releasing nuclear-processed material through
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Wind (dashed lines) and explosion (solid lines) yields for $^{26}$Al (top) and $^{60}$Fe (bottom) nuclides produced by solar-metallicity stars. Yields for stellar models from Limongi & Chieffi (2018, LC18; set M) and Ebinger et al. (2019, PUSH; based on the pre-explosion models of Woosley & Heger 2007) are displayed for non-rotating stars and rotating ones with an initial rotational velocity of 150 and 300 km s$^{-1}$.</figcaption></figure>

stellar winds and during their explosive phase. Models of stellar evolution that include a detailed nucleosynthesis network trace stars through all evolutionary stages and thus predict the nuclide yields both from the stellar winds and the SN explosion phase. An example of $^{26}$Al and $^{60}$Fe yields as a function of the initial stellar mass is presented in Fig. 2 for models from Limongi & Chieffi (2018) and Ebinger et al. (2019). Comparing the yields for non-rotating stars gives an idea of the systematic uncertainties of these models. The contribution of explosive nucleosynthesis (solid lines) typically amounts to ~$10^{-5}$ to $10^{-3}$ M$_{\odot}$ of ejected mass for both nuclides, with a mild dependence on the initial stellar mass. The wind contribution is completely negligible for $^{60}$Fe, quite in contrast with the case of $^{26}$Al where this is very significant for the high-end massive stars, even comparable to the contribution from the explosion. In Fig. 2, yields are also given as a function of the star's initial rotational velocity. Massive stars are known to be rotating objects (Głębocki & Gnaciński 2005), and yields are affected through mixing processes stimulated by stellar rotation. Indeed, rotation induces a slow mixing of both fresh fuel from the envelope into the burning core and of freshly synthesized material from the burning H-core into the envelope of the star. Stellar rotation also enhances the ejection into the ISM with stronger winds (Meynet & Maeder 2000). This leads to a larger $^{26}$Al wind contribution for stellar models that include such rotation (see Fig. 2). The large difference observed between rotating and non-rotating models from Limongi & Chieffi (2018) in the low-mass range ($\le 30$ M$_{\odot}$) is due to the treatment of rotational mixing and

the impact of a dust-driven wind on the stellar mass loss (Chieffi & Limongi 2013). The effect of rotation on the explosive yields is more difficult to assess, with no clear enhancement for rotating models except for the case of $^{60}$Fe for stars up to 20 M$_{\odot}$, where the yields are about 10 times larger than for the non-rotating models. In our population synthesis model, we follow the prescription of Prantzos et al. (2018) where the initial distribution of rotational velocity of stars is constrained from the study of the production of s-elements with a galactic chemical evolution model. We therefore consider that the probability for solar metallicity OB stars to have an initial rotational velocity of 0, 150, and 300 km s$^{-1}$ is 67, 32, and 1 per cent, respectively.

**Stellar explodability.** The nucleosynthetic output from a massive star strongly depends on its fate during its gravitational collapse at the end of its evolution. Massive stars that collapse and form black holes, either directly or through fallback, are not expected to enrich the ISM, while their successful explosion will disseminate freshly synthesized nuclear material into the ISM. Which star of a specific mass may experience a successful explosion and for which stellar mass this fails is an actively debated question (e.g. Foglizzo et al. 2015; Sukhbold et al. 2016). Even though there are several observations of ccSNe with indication of the progenitor mass (Ebinger et al. 2019, appendix A), the initial mass uncertainty and the low rate of ccSNe make it difficult to constrain the explodability from observations only. For simplicity, some models assume that massive stars collapse directly to black holes when their initial mass is greater than 25 M$_{\odot}$ (Limongi & Chieffi 2018), while others claim a transition mass in the range 100 $\le M \le 140$ M$_{\odot}$ (Janka 2012). Detailed numerical treatment of the explosion of massive stars has suggested that their explodability depends on the compactness in their pre-SN phase (O'Connor & Ott 2011), which leads to irregular gaps within the range of the stellar initial masses where massive stars undergo a successful explosion (e.g. Sukhbold et al. 2016; Ebinger et al. 2019).

**Fiducial model.** In this work, we use as a nominal set of parameters the IMF from Kroupa (2001, see also Kroupa 2002) with the stellar evolution prescription and yields from Limongi & Chieffi (2018). For explodability we assume that only stars below $\le 25$ M$_{\odot}$ explode as ccSN and subsequently release ejecta in the ISM. This corresponds to the case of set R defined in Limongi & Chieffi (2018), which is equivalent to set M (displayed in Fig. 2) where the explosion yield is set to zero above 25 M$_{\odot}$. Both stellar yields and lifetime depend on the star metallicity. However, since most OB associations are relatively young with typical ages below 50 Myr (Wright 2020), we adopt stellar yields and lifetimes for solar-metallicity stars. Concerning the initial distribution of rotational velocities of stars, we follow the prescription of Prantzos et al. (2018). The flexibility of our model allows us to switch for different IMFs, stellar and explosion yields, and explodability criteria very easily. We investigate the impact of changing these input parameters and describe this in Section 6.1.

### 3.1.2 Nuclide enrichment of the gas in OB association

As the starting point of our population synthesis model, we sample the IMF to generate the masses of the OB association members. We use the IMF function described in Kroupa (2001, 2002), and only massive stars ($M \ge 8$ M$_{\odot}$) are considered. We use random sampling for simplicity, considering the difference to optimal sampling (see e.g. Yan, Jerabkova & Kroupa 2023) rather insignificant for our purposes. The sampling procedure is repeated until a given total stellar content of the OB association is reproduced. This content can be deduced from the observations (see Section 4) or it can be
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Wind (dashed lines) and explosion (solid lines) yields for $^{26}$Al (top) and $^{60}$Fe (bottom) nuclides produced by solar-metallicity stars. Yields for stellar models from Limongi & Chieffi (2018, LC18; set M) and Ebinger et al. (2019, PUSH; based on the pre-explosion models of Woosley & Heger 2007) are displayed for non-rotating stars and rotating ones with an initial rotational velocity of 150 and 300 km s$^{-1}$.</figcaption></figure>

stellar winds and during their explosive phase. Models of stellar evolution that include a detailed nucleosynthesis network trace stars through all evolutionary stages and thus predict the nuclide yields both from the stellar winds and the SN explosion phase. An example of $^{26}$Al and $^{60}$Fe yields as a function of the initial stellar mass is presented in Fig. 2 for models from Limongi & Chieffi (2018) and Ebinger et al. (2019). Comparing the yields for non-rotating stars gives an idea of the systematic uncertainties of these models. The contribution of explosive nucleosynthesis (solid lines) typically amounts to ~$10^{-5}$ to $10^{-3}$ M$_{\odot}$ of ejected mass for both nuclides, with a mild dependence on the initial stellar mass. The wind contribution is completely negligible for $^{60}$Fe, quite in contrast with the case of $^{26}$Al where this is very significant for the high-end massive stars, even comparable to the contribution from the explosion. In Fig. 2, yields are also given as a function of the star's initial rotational velocity. Massive stars are known to be rotating objects (Głębocki & Gnaciński 2005), and yields are affected through mixing processes stimulated by stellar rotation. Indeed, rotation induces a slow mixing of both fresh fuel from the envelope into the burning core and of freshly synthesized material from the burning H-core into the envelope of the star. Stellar rotation also enhances the ejection into the ISM with stronger winds (Meynet & Maeder 2000). This leads to a larger $^{26}$Al wind contribution for stellar models that include such rotation (see Fig. 2). The large difference observed between rotating and non-rotating models from Limongi & Chieffi (2018) in the low-mass range ($\le 30$ M$_{\odot}$) is due to the treatment of rotational mixing and

the impact of a dust-driven wind on the stellar mass loss (Chieffi & Limongi 2013). The effect of rotation on the explosive yields is more difficult to assess, with no clear enhancement for rotating models except for the case of $^{60}$Fe for stars up to 20 M$_{\odot}$, where the yields are about 10 times larger than for the non-rotating models. In our population synthesis model, we follow the prescription of Prantzos et al. (2018) where the initial distribution of rotational velocity of stars is constrained from the study of the production of s-elements with a galactic chemical evolution model. We therefore consider that the probability for solar metallicity OB stars to have an initial rotational velocity of 0, 150, and 300 km s$^{-1}$ is 67, 32, and 1 per cent, respectively.

**Stellar explodability.** The nucleosynthetic output from a massive star strongly depends on its fate during its gravitational collapse at the end of its evolution. Massive stars that collapse and form black holes, either directly or through fallback, are not expected to enrich the ISM, while their successful explosion will disseminate freshly synthesized nuclear material into the ISM. Which star of a specific mass may experience a successful explosion and for which stellar mass this fails is an actively debated question (e.g. Foglizzo et al. 2015; Sukhbold et al. 2016). Even though there are several observations of ccSNe with indication of the progenitor mass (Ebinger et al. 2019, appendix A), the initial mass uncertainty and the low rate of ccSNe make it difficult to constrain the explodability from observations only. For simplicity, some models assume that massive stars collapse directly to black holes when their initial mass is greater than 25 M$_{\odot}$ (Limongi & Chieffi 2018), while others claim a transition mass in the range 100 $\le M \le 140$ M$_{\odot}$ (Janka 2012). Detailed numerical treatment of the explosion of massive stars has suggested that their explodability depends on the compactness in their pre-SN phase (O'Connor & Ott 2011), which leads to irregular gaps within the range of the stellar initial masses where massive stars undergo a successful explosion (e.g. Sukhbold et al. 2016; Ebinger et al. 2019).

**Fiducial model.** In this work, we use as a nominal set of parameters the IMF from Kroupa (2001, see also Kroupa 2002) with the stellar evolution prescription and yields from Limongi & Chieffi (2018). For explodability we assume that only stars below $\le 25$ M$_{\odot}$ explode as ccSN and subsequently release ejecta in the ISM. This corresponds to the case of set R defined in Limongi & Chieffi (2018), which is equivalent to set M (displayed in Fig. 2) where the explosion yield is set to zero above 25 M$_{\odot}$. Both stellar yields and lifetime depend on the star metallicity. However, since most OB associations are relatively young with typical ages below 50 Myr (Wright 2020), we adopt stellar yields and lifetimes for solar-metallicity stars. Concerning the initial distribution of rotational velocities of stars, we follow the prescription of Prantzos et al. (2018). The flexibility of our model allows us to switch for different IMFs, stellar and explosion yields, and explodability criteria very easily. We investigate the impact of changing these input parameters and describe this in Section 6.1.

### 3.1.2 Nuclide enrichment of the gas in OB association

As the starting point of our population synthesis model, we sample the IMF to generate the masses of the OB association members. We use the IMF function described in Kroupa (2001, 2002), and only massive stars ($M \ge 8$ M$_{\odot}$) are considered. We use random sampling for simplicity, considering the difference to optimal sampling (see e.g. Yan, Jerabkova & Kroupa 2023) rather insignificant for our purposes. The sampling procedure is repeated until a given total stellar content of the OB association is reproduced. This content can be deduced from the observations (see Section 4) or it can be
specified a priori as a total number of stars integrated over the full IMF mass range. For each massive star, an initial rotational velocity is randomly generated. Then, the lifetime of the star, which depends on the initial stellar mass and rotational velocity, is determined from stellar evolution models, and their nucleosynthesis yields during evolution are assembled. The contribution of stellar winds of massive stars is also taken into account in our model: for simplicity, stellar winds are assumed to be released at the end of the star lifetime since we are mainly interested in $^{60}$Fe that is not significantly produced by stellar winds. For massive stars ending their lives as ccSNe as controlled by the explodability criterion, their ejecta are released at their time of explosion.

It is usually assumed that massive stars within a stellar cluster can be considered as a coeval population (Lada 2005). Thus, the temporal evolution of the mass of a radionuclide $M(j, t)$ in the gas of an OB association is calculated as the sum of the individual contributions associated with each massive star:

$$M(j, t) = \sum_{i=1}^{n} \left[ Y_{i}^{\text{wind}}(j) + \eta_{i} \times Y_{i}^{\text{expl}}(j) \right] e^{-(t-t_{i})/\tau_{j}}, \quad (2)$$

where $Y_{i}^{\text{wind}}(j)$ and $Y_{i}^{\text{expl}}(j)$ are the wind and explosive yields, respectively, for nuclide j associated with the ith massive star with stellar lifetime $t_i$. $\eta_i$ is a parameter taking value of 0 or 1 whether the considered star explodes as a ccSN or not depending on the adopted explodability criterion. The exponential term reflects the free radioactive decay of nuclide j according to its corresponding lifetime $\tau_j$. This term should be set to 1 in case of stable nuclides.

In order to account for the stochastic nature of forming an OB association, our population synthesis model of an OB association is typically repeated 4000 times. This ensures to obtain a meaningful average for the temporal evolution of the nuclides abundance. As an example, the temporal evolution of the abundance of a few nuclides relevant to this work ($^{60}$Fe, $^{56}$Fe, and $^{26}$Al) is presented in Fig. 3 for two OB associations having a total stellar mass of $10^4 M_{\odot}$ (red) and $10^5 M_{\odot}$ (blue). For clarity sake, only five Monte Carlo realizations are shown in black solid lines for each case. The temporal evolution of the average mass of each nuclide computed for all realizations is represented as a solid coloured line. As expected, the total mass of a given nuclide scales linearly with the stellar content of the OB association, and the variance of the nuclide mass distribution is larger for the OB association with the lowest stellar content (red case).

The temporal evolution of the mass of a nuclide in the gas of the OB association shows distinct behaviour depending on its lifetime. In the case of stable nuclides (e.g. $^{56}$Fe) the abundance increases monotonically with time as a result of the cumulative effect of successive nucleosynthetic events. For radioactive nuclides a typical saw-tooth pattern is observed where sudden rises, corresponding to the enrichment of the OB association gas by the release of the wind and SNe yields, are followed by the radioactive decay until another nucleosynthetic event builds up on top of the previous one. The obtained pattern depends on how the radionuclide lifetime compares with the mean time between two successive ccSN explosions $\Delta t_{\text{ccSN}}$ (Côté et al. 2019). In the case of the $10^4 M_{\odot}$ OB association $\Delta t_{\text{ccSN}} \approx 0.5$ Myr. This is similar to the $^{26}$Al lifetime ($\tau = 1.03$ Myr) and the temporal variation of its mass exhibits a larger scatter than for $^{60}$Fe that has a longer lifetime ($\tau = 3.78$ Myr). Since the stellar content of the $10^5 M_{\odot}$ OB association is higher, the mean time between two successive SNe is lower ($\Delta t_{\text{ccSN}} \approx 0.1$ Myr), and much smaller than both the $^{26}$Al and $^{60}$Fe lifetimes. In that case, the deviation between individual realizations (black curve) and the average (blue curve) is significantly reduced.

<figure><img src="image_4.png" /><figcaption>Figure 3. Time evolution of the mass of $^{60}$Fe, $^{56}$Fe, and $^{26}$Al produced by an OB association with a stellar content of $10^4 M_{\odot}$ (red) and $10^5 M_{\odot}$ (blue). The time origin corresponds to the birth of the OB association that is evolved during 55 Myr. Five Monte Carlo realizations are represented in each case (black lines), and the associated average is shown by the coloured solid line.</figcaption></figure>

The nucleosynthetic enrichment of the gas of the OB association for a given nuclide may start at different epochs, as shown in Fig. 3. When a nuclide is produced significantly by stellar winds (e.g. $^{26}$Al and $^{56}$Fe), it enriches the OB association gas at early times. Since the wind contribution is released at explosion time in our model, the earliest possible release time occurs at ~3.5 Myr that corresponds to the stellar lifetime of the most massive stars of our model (120 $M_{\odot}$). In the case of nuclides that are produced during the SN phase only (e.g. $^{60}$Fe), the first contributing stars are the exploding stars with the shortest lifetime. This depends on the explodability criterion, which, in the present calculation, is such that stars with $M > 25 M_{\odot}$ directly collapse to form black holes with no explosive contribution to the nucleosynthesis. The earliest release time in that case is ~7.7 Myr that corresponds to the lifetime of a $25 M_{\odot}$ star.
specified a priori as a total number of stars integrated over the full IMF mass range. For each massive star, an initial rotational velocity is randomly generated. Then, the lifetime of the star, which depends on the initial stellar mass and rotational velocity, is determined from stellar evolution models, and their nucleosynthesis yields during evolution are assembled. The contribution of stellar winds of massive stars is also taken into account in our model: for simplicity, stellar winds are assumed to be released at the end of the star lifetime since we are mainly interested in $^{60}$Fe that is not significantly produced by stellar winds. For massive stars ending their lives as ccSNe as controlled by the explodability criterion, their ejecta are released at their time of explosion.

It is usually assumed that massive stars within a stellar cluster can be considered as a coeval population (Lada 2005). Thus, the temporal evolution of the mass of a radionuclide $M(j, t)$ in the gas of an OB association is calculated as the sum of the individual contributions associated with each massive star:

$$M(j, t) = \sum_{i=1}^{n} \left[ Y_{i}^{\text{wind}}(j) + \eta_{i} \times Y_{i}^{\text{expl}}(j) \right] e^{-(t-t_{i})/\tau_{j}}, \quad (2)$$

where $Y_{i}^{\text{wind}}(j)$ and $Y_{i}^{\text{expl}}(j)$ are the wind and explosive yields, respectively, for nuclide j associated with the ith massive star with stellar lifetime $t_i$. $\eta_i$ is a parameter taking value of 0 or 1 whether the considered star explodes as a ccSN or not depending on the adopted explodability criterion. The exponential term reflects the free radioactive decay of nuclide j according to its corresponding lifetime $\tau_j$. This term should be set to 1 in case of stable nuclides.

In order to account for the stochastic nature of forming an OB association, our population synthesis model of an OB association is typically repeated 4000 times. This ensures to obtain a meaningful average for the temporal evolution of the nuclides abundance. As an example, the temporal evolution of the abundance of a few nuclides relevant to this work ($^{60}$Fe, $^{56}$Fe, and $^{26}$Al) is presented in Fig. 3 for two OB associations having a total stellar mass of $10^4 M_{\odot}$ (red) and $10^5 M_{\odot}$ (blue). For clarity sake, only five Monte Carlo realizations are shown in black solid lines for each case. The temporal evolution of the average mass of each nuclide computed for all realizations is represented as a solid coloured line. As expected, the total mass of a given nuclide scales linearly with the stellar content of the OB association, and the variance of the nuclide mass distribution is larger for the OB association with the lowest stellar content (red case).

The temporal evolution of the mass of a nuclide in the gas of the OB association shows distinct behaviour depending on its lifetime. In the case of stable nuclides (e.g. $^{56}$Fe) the abundance increases monotonically with time as a result of the cumulative effect of successive nucleosynthetic events. For radioactive nuclides a typical saw-tooth pattern is observed where sudden rises, corresponding to the enrichment of the OB association gas by the release of the wind and SNe yields, are followed by the radioactive decay until another nucleosynthetic event builds up on top of the previous one. The obtained pattern depends on how the radionuclide lifetime compares with the mean time between two successive ccSN explosions $\Delta t_{ccSN}$ (Côté et al. 2019). In the case of the $10^4 M_{\odot}$ OB association $\Delta t_{ccSN} \approx 0.5$ Myr. This is similar to the $^{26}$Al lifetime ($\tau = 1.03$ Myr) and the temporal variation of its mass exhibits a larger scatter than for $^{60}$Fe that has a longer lifetime ($\tau = 3.78$ Myr). Since the stellar content of the $10^5 M_{\odot}$ OB association is higher, the mean time between two successive SNe is lower ($\Delta t_{ccSN} \approx 0.1$ Myr), and much smaller than both the $^{26}$Al and $^{60}$Fe lifetimes. In that case, the deviation between individual realizations (black curve) and the average (blue curve) is significantly reduced.

<figure><img src="image_4.png" /><figcaption>Figure 3. Time evolution of the mass of $^{60}$Fe, $^{56}$Fe, and $^{26}$Al produced by an OB association with a stellar content of $10^4 M_{\odot}$ (red) and $10^5 M_{\odot}$ (blue). The time origin corresponds to the birth of the OB association that is evolved during 55 Myr. Five Monte Carlo realizations are represented in each case (black lines), and the associated average is shown by the coloured solid line.</figcaption></figure>

The nucleosynthetic enrichment of the gas of the OB association for a given nuclide may start at different epochs, as shown in Fig. 3. When a nuclide is produced significantly by stellar winds (e.g. $^{26}$Al and $^{56}$Fe), it enriches the OB association gas at early times. Since the wind contribution is released at explosion time in our model, the earliest possible release time occurs at ~3.5 Myr that corresponds to the stellar lifetime of the most massive stars of our model (120 $M_{\odot}$). In the case of nuclides that are produced during the SN phase only (e.g. $^{60}$Fe), the first contributing stars are the exploding stars with the shortest lifetime. This depends on the explodability criterion, which, in the present calculation, is such that stars with $M > 25 M_{\odot}$ directly collapse to form black holes with no explosive contribution to the nucleosynthesis. The earliest release time in that case is ~7.7 Myr that corresponds to the lifetime of a $25 M_{\odot}$ star.
## 3.2 CR production and transport

Having assembled the interstellar content of $^{60}$Fe nuclei within a group of stars, we proceed to determine the fraction ending up in locally accelerated CRs, and propagate these then from the source through ISM toward the Solar system.

### 3.2.1 CR acceleration efficiency

Galactic CRs are widely believed to be produced by the DSA process in SN remnants, but alternative sources such as massive star clusters, pulsar wind nebulae, and the supermassive black hole at the Galactic Centre may also contribute to the CR population (see Gabici et al. 2019, and references therein). The DSA theory predicts that a fraction of interstellar particles of about $10^{-5}$ to $10^{-3}$ swept-up by a SN shock during the free expansion and the Sedov–Taylor phases become non-thermal CR particles (e.g. Blasi 2013). The CR-related gamma-ray luminosity of the Milky Way (Strong et al. 2010) suggests that the acceleration efficiency of protons, alpha particles, and other volatile elements is relatively low, of the order of $10^{-5}$ (Tatischeff et al. 2021). But refractory elements such as Al and Fe are significantly more abundant than volatile ones in the CR composition compared to the Solar system composition (Meyer, Drury & Ellison 1997), which requires an acceleration efficiency of the former of the order of a few $10^{-4}$. Such higher efficiency could plausibly be explained by a more efficient injection of dust grains than ions into the DSA process, due to the higher rigidity of the former (Ellison, Drury & Meyer 1997).

Massive star winds and SN ejecta within an OB association leave their sources in the form of hot, fast gas. As they expand, dust may form in dense clumps of stellar ejecta and condense a significant fraction of the refractory material. This has been suggested from infrared observations of SN 1987A (e.g. Matsuura et al. 2019). However, some or all of this dust could be efficiently destroyed by thermal sputtering in the SN reverse shock. This is suggested from the paucity of pre-solar grains with characteristic signatures of ccSNe as analysed in meteoritic materials (Nittler et al. 1996; Hoppe et al. 2019). Subsequently, stellar ejecta are expected to be diluted in the hot superbubble plasma encompassing the stellar association. However, in a young and compact star cluster embedded in a molecular cloud, a fraction of the ejecta could be rapidly incorporated in cold molecular gas (Vasileiadis, Nordlund & Bizzarro 2013). Gamma-ray observations of $^{26}$Al decay in nearby sources, such as the Scorpius–Centaurus (Sco–Cen) and the Orion–Eridanus superbubbles, provide a unique way of studying the interstellar transport of massive star ejecta (see Diehl et al. 2021, and references therein).

The acceleration efficiency of massive star ejecta by SN shocks propagating into the superbubble plasma thus depends in theory on several parameters including the size and age of the parent OB association, as well as on the efficiencies of dust production in stellar ejecta and destruction by thermal sputtering. In our model, all these poorly known processes are included in a single efficiency factor $\epsilon_{acc}$, which we vary from $10^{-6}$ to $10^{-4}$.

### 3.2.2 CR propagation

The general formalism of CR transport in the Galaxy includes particle diffusion, advection, ionization losses, spallation, and the radioactive decay of unstable nuclei (Ginzburg & Syrovatskii 1964). The specific transport of $^{60}$Fe CRs has been recently studied by Morlino & Amato (2020) within the framework of a disc–halo diffusion model. They

used for the CR diffusion coefficient, assumed to be the same in the disc and the halo (see also Evoli et al. 2019):

$$D(R) = \beta D_0 \frac{(R/GV)^\delta}{[1 + (R/R_b)^{\Delta\delta/s}]^\delta}, \quad (3)$$

where R is the particle rigidity, $D_0 = 3.08 \times 10^{28} \text{ cm}^2 \text{ s}^{-1}$, $\delta = 0.54$, $\Delta\delta = 0.2$, $s = 0.1$, and $R_b = 312$ GV. For $^{60}$Fe CRs of $\approx 523$ MeV nucleon$^{-1}$ (the mean LISM energy of the $^{60}$Fe nuclei detected by ACE/CRIS; see Section 2), we have $D \approx 4.0 \times 10^{28} \text{ cm}^2 \text{ s}^{-1}$.

However, the diffusion coefficient in the LISM is very uncertain. It depends in particular on the structure of the interstellar magnetic field between the nearby sources and the Solar system. In addition, the spatial diffusion coefficient in an active superbubble environment is expected to be lower than that in the average ISM ($D_0$ in the range $10^{27}-10^{28} \text{ cm}^2 \text{ s}^{-1}$; see Vieu et al. 2022). Moreover, in order to escape from a superbubble, CRs must diffuse mainly perpendicularly to the compressed magnetic field in the supershell, which could enhance the particle confinement in the hot plasma. Detailed modelling of these effects is beyond the scope of this paper. Here, we assume as a nominal value the same diffusion coefficient as Evoli et al. (2019) and Morlino & Amato (2020), and study in Section 6.1 the impact on the results of reducing D by an order of magnitude.

We now compare the time-scales for the various processes involved in the transport of $^{60}$Fe ions in the Galactic disc, assuming the half-thickness of the disc to be $h = 150$ pc. With $D \approx 4.0 \times 10^{28} \text{ cm}^2 \text{ s}^{-1}$ (as obtained from equation 3), the diffusion time-scale of $^{60}$Fe CRs over this distance is

$$\tau_{\text{diff}} = h^2/D = 1.7 \times 10^5 \text{ yr}, \quad (4)$$

which is significantly shorter than the CR advection time-scale:

$$\tau_{\text{adv}} = h/u_0 = 2.9 \times 10^7 \text{ yr}, \quad (5)$$

where $u_0 = 5 \text{ km s}^{-1}$ is the typical CR advection velocity (Morlino & Amato 2020). Advection can thus be neglected.

The time-scale for the catastrophic losses of $^{60}$Fe nuclei by nuclear spallation reactions in the ISM can be estimated as

$$\begin{aligned} \tau_{\text{spal}} &= \frac{1}{n_H v [\sigma_H + (n_{\text{He}}/n_H)\sigma_{\text{He}}]} \\ &= 1.54 \times 10^7 (n_H/0.1 \text{ cm}^{-3})^{-1} \text{ yr}, \end{aligned} \quad (6)$$

where $n_H$ is the average ISM density into which $^{60}$Fe CRs propagate from their sources to the Solar system, and $\sigma_H$ and $\sigma_{\text{He}}$ are the total reaction cross-sections for fast ions propagating in interstellar H and He, respectively (we assume 90 per cent H and 10 per cent He by number). We used for these cross-sections the universal parametrization of Tripathi, Cucinotta & Wilson (1996, 1999). We then found for the interaction mean free path of $^{60}$Fe nuclei of $523$ MeV nucleon$^{-1}$ in the ISM $\lambda_{\text{spal}} = 2.69 \text{ g cm}^{-2}$, which is 5 per cent above the value reported by Binns et al. (2016), $\lambda_{\text{spal}} = 2.56 \text{ g cm}^{-2}$. The total loss time-scale of $^{60}$Fe CRs in the ISM is given by

$$\tau_{\text{loss}} = \left( \frac{1}{\tau_{\text{spal}}} + \frac{1}{\tau_{\text{decay}}} \right)^{-1}, \quad (7)$$

where

$$\tau_{\text{decay}} = \gamma \tau_{\text{decay},0} = (5.90 \pm 0.09) \times 10^6 \text{ yr}. \quad (8)$$

Here $\gamma = 1.56$ is the Lorentz factor and $\tau_{\text{decay},0} = 3.78 \pm 0.06$ Myr is the mean lifetime for radioactive decay of $^{60}$Fe at rest.

$^{60}$Fe ions originating from the younger subgroups of the nearby Sco–Cen OB association are expected to have propagated mainly in
## 3.2 CR production and transport

Having assembled the interstellar content of $^{60}$Fe nuclei within a group of stars, we proceed to determine the fraction ending up in locally accelerated CRs, and propagate these then from the source through ISM toward the Solar system.

### 3.2.1 CR acceleration efficiency

Galactic CRs are widely believed to be produced by the DSA process in SN remnants, but alternative sources such as massive star clusters, pulsar wind nebulae, and the supermassive black hole at the Galactic Centre may also contribute to the CR population (see Gabici et al. 2019, and references therein). The DSA theory predicts that a fraction of interstellar particles of about $10^{-5}$ to $10^{-3}$ swept-up by a SN shock during the free expansion and the Sedov–Taylor phases become non-thermal CR particles (e.g. Blasi 2013). The CR-related gamma-ray luminosity of the Milky Way (Strong et al. 2010) suggests that the acceleration efficiency of protons, alpha particles, and other volatile elements is relatively low, of the order of $10^{-5}$ (Tatischeff et al. 2021). But refractory elements such as Al and Fe are significantly more abundant than volatile ones in the CR composition compared to the Solar system composition (Meyer, Drury & Ellison 1997), which requires an acceleration efficiency of the former of the order of a few $10^{-4}$. Such higher efficiency could plausibly be explained by a more efficient injection of dust grains than ions into the DSA process, due to the higher rigidity of the former (Ellison, Drury & Meyer 1997).

Massive star winds and SN ejecta within an OB association leave their sources in the form of hot, fast gas. As they expand, dust may form in dense clumps of stellar ejecta and condense a significant fraction of the refractory material. This has been suggested from infrared observations of SN 1987A (e.g. Matsuura et al. 2019). However, some or all of this dust could be efficiently destroyed by thermal sputtering in the SN reverse shock. This is suggested from the paucity of pre-solar grains with characteristic signatures of ccSNe as analysed in meteoritic materials (Nittler et al. 1996; Hoppe et al. 2019). Subsequently, stellar ejecta are expected to be diluted in the hot superbubble plasma encompassing the stellar association. However, in a young and compact star cluster embedded in a molecular cloud, a fraction of the ejecta could be rapidly incorporated in cold molecular gas (Vasileiadis, Nordlund & Bizzarro 2013). Gamma-ray observations of $^{26}$Al decay in nearby sources, such as the Scorpius–Centaurus (Sco–Cen) and the Orion–Eridanus superbubbles, provide a unique way of studying the interstellar transport of massive star ejecta (see Diehl et al. 2021, and references therein).

The acceleration efficiency of massive star ejecta by SN shocks propagating into the superbubble plasma thus depends in theory on several parameters including the size and age of the parent OB association, as well as on the efficiencies of dust production in stellar ejecta and destruction by thermal sputtering. In our model, all these poorly known processes are included in a single efficiency factor $\epsilon_{acc}$, which we vary from $10^{-6}$ to $10^{-4}$.

### 3.2.2 CR propagation

The general formalism of CR transport in the Galaxy includes particle diffusion, advection, ionization losses, spallation, and the radioactive decay of unstable nuclei (Ginzburg & Syrovatskii 1964). The specific transport of $^{60}$Fe CRs has been recently studied by Morlino & Amato (2020) within the framework of a disc–halo diffusion model. They

used for the CR diffusion coefficient, assumed to be the same in the disc and the halo (see also Evoli et al. 2019):

$$D(R) = \beta D_0 \frac{(R/GV)^\delta}{[1 + (R/R_b)^{\Delta\delta/s}]^s} \quad (3)$$

where R is the particle rigidity, $D_0 = 3.08 \times 10^{28} \text{ cm}^2 \text{ s}^{-1}$, $\delta = 0.54$, $\Delta\delta = 0.2$, $s = 0.1$, and $R_b = 312$ GV. For $^{60}$Fe CRs of $\approx 523$ MeV nucleon$^{-1}$ (the mean LISM energy of the $^{60}$Fe nuclei detected by ACE/CRIS; see Section 2), we have $D \approx 4.0 \times 10^{28} \text{ cm}^2 \text{ s}^{-1}$.

However, the diffusion coefficient in the LISM is very uncertain. It depends in particular on the structure of the interstellar magnetic field between the nearby sources and the Solar system. In addition, the spatial diffusion coefficient in an active superbubble environment is expected to be lower than that in the average ISM ($D_0$ in the range $10^{27}-10^{28} \text{ cm}^2 \text{ s}^{-1}$; see Vieu et al. 2022). Moreover, in order to escape from a superbubble, CRs must diffuse mainly perpendicularly to the compressed magnetic field in the supershell, which could enhance the particle confinement in the hot plasma. Detailed modelling of these effects is beyond the scope of this paper. Here, we assume as a nominal value the same diffusion coefficient as Evoli et al. (2019) and Morlino & Amato (2020), and study in Section 6.1 the impact on the results of reducing D by an order of magnitude.

We now compare the time-scales for the various processes involved in the transport of $^{60}$Fe ions in the Galactic disc, assuming the half-thickness of the disc to be $h = 150$ pc. With $D \approx 4.0 \times 10^{28} \text{ cm}^2 \text{ s}^{-1}$ (as obtained from equation 3), the diffusion time-scale of $^{60}$Fe CRs over this distance is

$$\tau_{\text{diff}} = h^2/D = 1.7 \times 10^5 \text{ yr}, \quad (4)$$

which is significantly shorter than the CR advection time-scale:

$$\tau_{\text{adv}} = h/u_0 = 2.9 \times 10^7 \text{ yr}, \quad (5)$$

where $u_0 = 5 \text{ km s}^{-1}$ is the typical CR advection velocity (Morlino & Amato 2020). Advection can thus be neglected.

The time-scale for the catastrophic losses of $^{60}$Fe nuclei by nuclear spallation reactions in the ISM can be estimated as

$$\begin{aligned} \tau_{\text{spal}} &= \frac{1}{n_H v [\sigma_H + (n_{\text{He}}/n_H)\sigma_{\text{He}}]} \\ &= 1.54 \times 10^7 (n_H/0.1 \text{ cm}^{-3})^{-1} \text{ yr}, \end{aligned} \quad (6)$$

where $n_H$ is the average ISM density into which $^{60}$Fe CRs propagate from their sources to the Solar system, and $\sigma_H$ and $\sigma_{\text{He}}$ are the total reaction cross-sections for fast ions propagating in interstellar H and He, respectively (we assume 90 per cent H and 10 per cent He by number). We used for these cross-sections the universal parametrization of Tripathi, Cucinotta & Wilson (1996, 1999). We then found for the interaction mean free path of $^{60}$Fe nuclei of 523 MeV nucleon$^{-1}$ in the ISM $\lambda_{\text{spal}} = 2.69 \text{ g cm}^{-2}$, which is 5 per cent above the value reported by Binns et al. (2016), $\lambda_{\text{spal}} = 2.56 \text{ g cm}^{-2}$. The total loss time-scale of $^{60}$Fe CRs in the ISM is given by

$$\tau_{\text{loss}} = \left( \frac{1}{\tau_{\text{spal}}} + \frac{1}{\tau_{\text{decay}}} \right)^{-1}, \quad (7)$$

where

$$\tau_{\text{decay}} = \gamma \tau_{\text{decay},0} = (5.90 \pm 0.09) \times 10^6 \text{ yr}. \quad (8)$$

Here $\gamma = 1.56$ is the Lorentz factor and $\tau_{\text{decay},0} = 3.78 \pm 0.06$ Myr is the mean lifetime for radioactive decay of $^{60}$Fe at rest.

$^{60}$Fe ions originating from the younger subgroups of the nearby Sco–Cen OB association are expected to have propagated mainly in
the low-density gas ($n_H \sim 0.1 \text{ cm}^{-3}$) filling the Local Hot Bubble (Zucker et al. 2022), and thus have suffered negligible catastrophic losses. But $^{60}$Fe ions coming from more distant OB associations (e.g. Orion, Cygnus OB2, etc.) and diffusing in the Galactic disc could have passed through denser regions (superbubble shells in particular) and seen on average ISM densities of $n_H \sim 1 \text{ cm}^{-3}$. However, $^{60}$Fe ions produced in distant associations should have mainly propagated in the low-density halo of the Galaxy before reaching the Solar system and then seen $n_H \sim 0.1 \text{ cm}^{-3}$ (see Morlino & Amato 2020). We thus adopt $n_H = 0.1 \text{ cm}^{-3}$ as the nominal value in our model, and will discuss the effect of changing the density parameter in Section 6.1. For $n_H = 0.1 \text{ cm}^{-3}$, $\tau_{\text{loss}} = 4.27$ Myr.

The ionization energy loss time-scale for $^{60}$Fe ions of kinetic energy $E = 523$ MeV nucleon$^{-1}$ is

$$ \tau_{\text{ion}} = E/(dE/dt)_{\text{ion}} = 5.7 \times 10^7 (n_{\text{H}}/0.1 \text{ cm}^{-3})^{-1} \text{ yr}, \quad (9) $$

where $(dE/dt)_{\text{ion}}$ is the ionization energy-loss rate, which is calculated from Mannheim & Schlickeiser (1994, equation 4.24). As for the catastrophic energy losses, the significance of the ionization energy losses could depend on the OB association from which the $^{60}$Fe CRs originate. However, we see from equations (6) and 9 that $\tau_{\text{ion}} > \tau_{\text{spal}}$ whatever $n_H$, so that the ionization losses can always be neglected in front of the catastrophic losses.

So finally we consider a simple propagation model where accelerated ions, when escaping from their source, diffuse isotropically in the ISM and suffer both catastrophic and radioactive losses. We use as nominal set of input parameters $D_0 = 3.08 \times 10^{28} \text{ cm}^2 \text{ s}^{-1}$ (equation 3), $n_H = 0.1 \text{ cm}^{-3}$, and $\epsilon_{\text{acc}} = 10^{-5}$ (Section 3.2.1), and we will study the impact of changing these parameters in Section 6.1. Future work could take into account in more detail the specific locations of the local OB associations and consider non-isotropic diffusion from MHD modelling of the LISM, but this is beyond the scope of this paper.

### 3.2.3 CR density in the LISM

In order to compare the observed density of CRs by ACE/CRIS with our model, we need to compute the number density of CRs $n(j)$ for a given nuclide $j$. This is computed as the sum of the contribution of each ccSN explosion from our model, where each SN accelerates with an efficiency $\epsilon_{\text{acc}}$ (Section 3.2.1), the number of atoms $N_0(j, t_i)$ of nuclide $j$ present at the explosion time $t_i$ of the ith ccSN counted from the birth time of the OB association. This number of atoms is deduced from the temporal evolution of the mass of $j$ given in equation (2) assuming that ccSNe do accelerate their own winds, since they are released prior to the collapse, but not their own ejecta (Wiedenbeck et al. 1999).

CRIS measurements of $^{56}$Fe and $^{60}$Fe CRs were performed between ~195 and ~500 MeV nucleon$^{-1}$, corresponding to $\Delta E \sim 400-700$ MeV nucleon$^{-1}$ in the LISM (Section 2). From the Fe source spectrum obtained by Boschini et al. (2021), we find the fraction of Fe nuclei released with energies in the range $\Delta E$ to be $\epsilon_{\Delta E} \approx 5\%$. This quantity slightly depends on the assumed minimum CR energy used to calculate the total number of accelerated Fe. Thus, we have $\epsilon_{\Delta E} = 4.6\%$ and $7.0\%$ for $E_{\text{min}} = 1$ and $10$ MeV nucleon$^{-1}$, respectively.

The resulting CR population must then diffuse across the distance $d$ between the OB association and the Solar system during a time $\Delta t_i = t_{\text{OB}} - t_i$, where $t_{\text{OB}}$ is the age of the association. The contribution of the ith ccSN to the total number density is obtained from the solution

<figure><img src="image_5.png" /><figcaption>Figure 4. CR density for three different nuclides ($^{56}$Fe, $^{60}$Fe, and $^{26}$Al) as a function of propagation time for OB associations located at different distances $d$. The calculations are performed for a single ccSN event accelerating the same number of atoms for each nuclide corresponding to $1.4 \times 10^{-4} \text{ M}_\odot$ of $^{60}$Fe. All calculations are performed with a diffusion coefficient $D = 4 \times 10^{28} \text{ cm}^2 \text{ s}^{-1}$, a total efficiency $\epsilon_{\text{acc}} \times \epsilon_{\Delta E} = 5 \times 10^{-7}$, and an average density of the ISM $n_H = 0.1 \text{ cm}^{-3}$.</figcaption></figure>

of the diffusion equation and reads

$$ n_i(j) = \frac{N_0(j, t_i)\epsilon_{\text{acc}}\epsilon_{\Delta E}(j)}{(4\pi D \Delta t_i)^{3/2}} \exp\left(-\frac{d^2}{4D\Delta t_i}\right) \\ \times \exp\left(-\frac{\Delta t_i}{\tau_{\text{loss}}(j)}\right), \quad (10) $$

where the last exponential decay term accounts for the catastrophic and radioactive losses (when $j$ is a radioactive species).

The CR density obtained from equation (10) is displayed in Fig. 4 as a function of the propagation time $\Delta t_i$ for three different distances $d$ of the parent OB association. Calculations are performed for one stable nuclide ($^{56}$Fe) and two radionuclides ($^{60}$Fe and $^{26}$Al). For all cases we consider for illustration purpose the same number of atoms in the parent superbubble plasma, $N_0 = 2.8 \times 10^{51}$, which corresponds to $1.4 \times 10^{-4} \text{ M}_\odot$ of $^{60}$Fe. This value is obtained from an average of the Limongi & Chieffi (2018) yields over the IMF from Kroupa (2001) and the initial rotational velocity from Prantzos et al. (2018). Fig. 4 exhibits the expected time evolution of the CR number density at the Solar system location from sources at various distances, with a sharp rise and a longer decay. For the low average ISM density considered, $n_H = 0.1 \text{ cm}^{-3}$, the catastrophic losses are negligible with respect to the radioactive decay losses of $^{26}$Al and $^{60}$Fe, which explains why for $\Delta t_i \gtrsim 1$ Myr the density of $^{26}$Al ($\tau_{\text{decay},0} = 1.03$ Myr) decreases faster than that of $^{60}$Fe ($\tau_{\text{decay},0} = 3.78$ Myr) and $^{56}$Fe.

We see in Fig. 4 that the CR density at maximum varies a lot with the source distance, e.g. by more than three orders of magnitude from $d = 140$ pc (the approximate distance of the Sco–Cen association) to $d = 1.5$ kpc (the approximate distance of Cyg OB2). The time when the CR density reaches its maximum can be obtained after cancelling the derivative of $n_i(j)$ from equation (10):

$$ t_{\text{max}} = \frac{\sqrt{D\tau(4d^2 + 9D\tau)} - 3D\tau}{4D}. \quad (11) $$

For stable nuclei and when catastrophic losses are negligible, we retrieve the well-known formula

$$ t_{\text{max}} = \frac{d^2}{6D}. \quad (12) $$
the low-density gas ($n_H \sim 0.1 \text{ cm}^{-3}$) filling the Local Hot Bubble (Zucker et al. 2022), and thus have suffered negligible catastrophic losses. But $^{60}$Fe ions coming from more distant OB associations (e.g. Orion, Cygnus OB2, etc.) and diffusing in the Galactic disc could have passed through denser regions (superbubble shells in particular) and seen on average ISM densities of $n_H \sim 1 \text{ cm}^{-3}$. However, $^{60}$Fe ions produced in distant associations should have mainly propagated in the low-density halo of the Galaxy before reaching the Solar system and then seen $n_H \sim 0.1 \text{ cm}^{-3}$ (see Morlino & Amato 2020). We thus adopt $n_H = 0.1 \text{ cm}^{-3}$ as the nominal value in our model, and will discuss the effect of changing the density parameter in Section 6.1. For $n_H = 0.1 \text{ cm}^{-3}$, $\tau_{\text{loss}} = 4.27$ Myr.

The ionization energy loss time-scale for $^{60}$Fe ions of kinetic energy $E = 523$ MeV nucleon$^{-1}$ is

$$ \tau_{\text{ion}} = E/(dE/dt)_{\text{ion}} = 5.7 \times 10^7 (n_{\text{H}}/0.1 \text{ cm}^{-3})^{-1} \text{ yr}, \quad (9) $$

where $(dE/dt)_{\text{ion}}$ is the ionization energy-loss rate, which is calculated from Mannheim & Schlickeiser (1994, equation 4.24). As for the catastrophic energy losses, the significance of the ionization energy losses could depend on the OB association from which the $^{60}$Fe CRs originate. However, we see from equations (6) and 9 that $\tau_{\text{ion}} > \tau_{\text{spal}}$ whatever $n_H$, so that the ionization losses can always be neglected in front of the catastrophic losses.

So finally we consider a simple propagation model where accelerated ions, when escaping from their source, diffuse isotropically in the ISM and suffer both catastrophic and radioactive losses. We use as nominal set of input parameters $D_0 = 3.08 \times 10^{28} \text{ cm}^2 \text{ s}^{-1}$ (equation 3), $n_H = 0.1 \text{ cm}^{-3}$, and $\epsilon_{\text{acc}} = 10^{-5}$ (Section 3.2.1), and we will study the impact of changing these parameters in Section 6.1. Future work could take into account in more detail the specific locations of the local OB associations and consider non-isotropic diffusion from MHD modelling of the LISM, but this is beyond the scope of this paper.

### 3.2.3 CR density in the LISM

In order to compare the observed density of CRs by ACE/CRIS with our model, we need to compute the number density of CRs $n(j)$ for a given nuclide $j$. This is computed as the sum of the contribution of each ccSN explosion from our model, where each SN accelerates with an efficiency $\epsilon_{\text{acc}}$ (Section 3.2.1), the number of atoms $N_0(j, t_i)$ of nuclide $j$ present at the explosion time $t_i$ of the $i$th ccSN counted from the birth time of the OB association. This number of atoms is deduced from the temporal evolution of the mass of $j$ given in equation (2) assuming that ccSNe do accelerate their own winds, since they are released prior to the collapse, but not their own ejecta (Wiedenbeck et al. 1999).

CRIS measurements of $^{56}$Fe and $^{60}$Fe CRs were performed between ~195 and ~500 MeV nucleon$^{-1}$, corresponding to $\Delta E \sim 400-700$ MeV nucleon$^{-1}$ in the LISM (Section 2). From the Fe source spectrum obtained by Boschini et al. (2021), we find the fraction of Fe nuclei released with energies in the range $\Delta E$ to be $\epsilon_{\Delta E} \approx 5\%$. This quantity slightly depends on the assumed minimum CR energy used to calculate the total number of accelerated Fe. Thus, we have $\epsilon_{\Delta E} = 4.6\%$ and $7.0\%$ for $E_{\text{min}} = 1$ and $10$ MeV nucleon$^{-1}$, respectively.

The resulting CR population must then diffuse across the distance $d$ between the OB association and the Solar system during a time $\Delta t_i = t_{\text{OB}} - t_i$, where $t_{\text{OB}}$ is the age of the association. The contribution of the $i$th ccSN to the total number density is obtained from the solution

<figure><img src="image_5.png" /><figcaption>**Figure 4.** CR density for three different nuclides ($^{56}$Fe, $^{60}$Fe, and $^{26}$Al) as a function of propagation time for OB associations located at different distances $d$. The calculations are performed for a single ccSN event accelerating the same number of atoms for each nuclide corresponding to $1.4 \times 10^{-4} \text{ M}_\odot$ of $^{60}$Fe. All calculations are performed with a diffusion coefficient $D = 4 \times 10^{28} \text{ cm}^2 \text{ s}^{-1}$, a total efficiency $\epsilon_{\text{acc}} \times \epsilon_{\Delta E} = 5 \times 10^{-7}$, and an average density of the ISM $n_H = 0.1 \text{ cm}^{-3}$.</figcaption></figure>

of the diffusion equation and reads

$$ n_i(j) = \frac{N_0(j, t_i)\epsilon_{\text{acc}}\epsilon_{\Delta E}(j)}{(4\pi D \Delta t_i)^{3/2}} \exp\left(-\frac{d^2}{4D\Delta t_i}\right) \\ \times \exp\left(-\frac{\Delta t_i}{\tau_{\text{loss}}(j)}\right), \quad (10) $$

where the last exponential decay term accounts for the catastrophic and radioactive losses (when $j$ is a radioactive species).

The CR density obtained from equation (10) is displayed in Fig. 4 as a function of the propagation time $\Delta t_i$ for three different distances $d$ of the parent OB association. Calculations are performed for one stable nuclide ($^{56}$Fe) and two radionuclides ($^{60}$Fe and $^{26}$Al). For all cases we consider for illustration purpose the same number of atoms in the parent superbubble plasma, $N_0 = 2.8 \times 10^{51}$, which corresponds to $1.4 \times 10^{-4} \text{ M}_\odot$ of $^{60}$Fe. This value is obtained from an average of the Limongi & Chieffi (2018) yields over the IMF from Kroupa (2001) and the initial rotational velocity from Prantzos et al. (2018). Fig. 4 exhibits the expected time evolution of the CR number density at the Solar system location from sources at various distances, with a sharp rise and a longer decay. For the low average ISM density considered, $n_H = 0.1 \text{ cm}^{-3}$, the catastrophic losses are negligible with respect to the radioactive decay losses of $^{26}$Al and $^{60}$Fe, which explains why for $\Delta t_i \gtrsim 1$ Myr the density of $^{26}$Al ($\tau_{\text{decay},0} = 1.03$ Myr) decreases faster than that of $^{60}$Fe ($\tau_{\text{decay},0} = 3.78$ Myr) and $^{56}$Fe.

We see in Fig. 4 that the CR density at maximum varies a lot with the source distance, e.g. by more than three orders of magnitude from $d = 140$ pc (the approximate distance of the Sco–Cen association) to $d = 1.5$ kpc (the approximate distance of Cyg OB2). The time when the CR density reaches its maximum can be obtained after cancelling the derivative of $n_i(j)$ from equation (10):

$$ t_{\text{max}} = \frac{\sqrt{D\tau(4d^2 + 9D\tau)} - 3D\tau}{4D}. \quad (11) $$

For stable nuclei and when catastrophic losses are negligible, we retrieve the well-known formula

$$ t_{\text{max}} = \frac{d^2}{6D}. \quad (12) $$
<table><caption>**Table 1.** Properties of local OB associations considered in this work. Distances and ages come from Wright (2020). The richness ($N_{*} \ge 8 M_{\odot}$ at formation time) and the number of core-collapse supernovae having already exploded ($N_{ccSN}$) are calculated from the number of observed stars and the OB association age (see text). $N_{ccSN}$ is reported for two different explodability criteria.</caption><thead><tr><th>Association</th><th>Distance (pc)</th><th>Age (Myr)</th><th>Number of observed stars</th><th>Ref.</th><th>N<sub>*</sub> (≥ 8 M<sub>☉</sub>)</th><th>N<sub>ccSN</sub><sup>LC18</sup></th><th>N<sub>ccSN</sub><sup>Sukhbold</sup></th></tr></thead><tbody><tr><td>Sco-Cen: US</td><td>143 ± 6</td><td>10 ± 7</td><td>N(OB) = 49</td><td>DZ99</td><td>13.9 ± 3.5</td><td>1.8 ± 1.8</td><td>1.6 ± 1.4</td></tr><tr><td>Sco-Cen: UCL</td><td>136 ± 5</td><td>16 ± 7</td><td>N(OB) = 66</td><td>DZ99</td><td>19.4 ± 4.3</td><td>5.1 ± 3.2</td><td>4.2 ± 2.7</td></tr><tr><td>Sco-Cen: LCC</td><td>115 ± 4</td><td>15 ± 6</td><td>N(OB) = 42</td><td>DZ99</td><td>12.4 ± 3.3</td><td>3.1 ± 2.1</td><td>2.5 ± 1.7</td></tr><tr><td>Ori OB1a</td><td>~360</td><td>8–12</td><td>N(4–15 M<sub>☉</sub>) = 53</td><td>B94</td><td>26.5 ± 4.9</td><td>2.2 ± 1.6</td><td>2.2 ± 1.4</td></tr><tr><td>Ori OB1b</td><td>360–420</td><td>2–8</td><td>N(4–120 M<sub>☉</sub>) = 45</td><td>B94</td><td>19.4 ± 3.6</td><td>0.5 ± 0.1</td><td>0.8 ± 0.6</td></tr><tr><td>Ori OB1c</td><td>~385</td><td>2–6</td><td>N(7–36 M<sub>☉</sub>) = 23</td><td>B94</td><td>22.7 ± 2.5</td><td>0</td><td>0.8 ± 0.6</td></tr><tr><td>Ori OB1d</td><td>~380</td><td>1–2</td><td>N(&gt; 1 M<sub>☉</sub>) = 145</td><td>H97</td><td>10.0 ± 3.0</td><td>0</td><td>0</td></tr><tr><td>Per OB1</td><td>~1830</td><td>8–11</td><td>N(OB) = 133</td><td>MD17</td><td>36.1 ± 5.6</td><td>2.2 ± 1.6</td><td>2.6 ± 1.1</td></tr><tr><td>Per OB2</td><td>296 ± 17</td><td>1–10</td><td>N(1–17 M<sub>☉</sub>) = 800</td><td>B02</td><td>53.9 ± 7.2</td><td>0.9 ± 1.1</td><td>1.8 ± 1.5</td></tr><tr><td>Per OB3</td><td>175 ± 3</td><td>50</td><td>N(OB) = 30</td><td>DZ99</td><td>10.8 ± 3.7</td><td>7.9 ± 3.4</td><td>7.1 ± 3.2</td></tr><tr><td>Cyg OB2</td><td>1350–1750</td><td>1–7</td><td>N(O) = 78</td><td>B20</td><td>241.2 ± 47.8</td><td>0</td><td>3.8 ± 3.4</td></tr><tr><td>Cyg OB4</td><td>~800</td><td>~8.3</td><td>N(OB) = 2</td><td>MD17</td><td>1.3 ± 0.8</td><td>0.5 ± 0.2</td><td>0.5 ± 0.2</td></tr><tr><td>Cyg OB7</td><td>~630</td><td>1–13</td><td>N(OB) = 25</td><td>MD17</td><td>7.2 ± 2.4</td><td>0.7 ± 0.6</td><td>0.8 ± 0.6</td></tr><tr><td>Cyg OB9</td><td>~960</td><td>2–4</td><td>N(OB) = 31</td><td>MD17</td><td>8.5 ± 2.5</td><td>0</td><td>0.5 ± 0.2</td></tr><tr><td>Vel OB2</td><td>411 ± 12</td><td>10–30</td><td>N(&gt; 2.5 M<sub>☉</sub>) = 72</td><td>A18</td><td>18.4 ± 4.4</td><td>6.8 ± 3.9</td><td>5.6 ± 3.5</td></tr><tr><td>Trumpler 10</td><td>372 ± 23</td><td>45–50</td><td>N(OB) = 22</td><td>DZ99</td><td>8.3 ± 3.2</td><td>6.7 ± 2.8</td><td>6.1 ± 2.6</td></tr><tr><td>Cas-Tau</td><td>125–300</td><td>~50</td><td>N(OB) = 83</td><td>DZ99</td><td>29.4 ± 6.6</td><td>22.1 ± 6.5</td><td>20.0 ± 6.2</td></tr><tr><td>Lac OB1</td><td>368 ± 17</td><td>2–25</td><td>N(OB) = 36</td><td>DZ99</td><td>10.7 ± 3.1</td><td>2.5 ± 2.3</td><td>2.2 ± 2.0</td></tr><tr><td>Cep OB2</td><td>~730</td><td>5</td><td>N(OB) = 56</td><td>DZ99</td><td>15.5 ± 3.4</td><td>0</td><td>0.8 ± 0.5</td></tr><tr><td>Cep OB3</td><td>~700</td><td>5–8</td><td>N(OB) = 25</td><td>MD17</td><td>7.2 ± 2.4</td><td>0.5 ± 0.1</td><td>0.7 ± 0.4</td></tr><tr><td>Cep OB4</td><td>~660</td><td>1–6</td><td>N(OB) = 7</td><td>MD17</td><td>2.5 ± 1.2</td><td>0</td><td>0.5 ± 0.2</td></tr><tr><td>Cep OB6</td><td>270 ± 12</td><td>~50</td><td>N(OB) = 6</td><td>DZ99</td><td>2.9 ± 1.8</td><td>2.3 ± 1.5</td><td>2.1 ± 1.4</td></tr><tr><td>Collinder 121</td><td>543 ± 23</td><td>5</td><td>N(OB) = 87</td><td>DZ99</td><td>22.9 ± 4.1</td><td>0</td><td>0.9 ± 8.8</td></tr><tr><td>Cam OB1</td><td>~800</td><td>7–14</td><td>N(OB) = 45</td><td>MD17</td><td>12.9 ± 3.3</td><td>1.6 ± 1.3</td><td>1.5 ± 1.1</td></tr><tr><td>Mon OB1</td><td>~580</td><td>1–10</td><td>N(OB) = 6</td><td>MD17</td><td>2.3 ± 1.2</td><td>0.5 ± 0.1</td><td>0.5 ± 0.2</td></tr></tbody></table>

Note. References: A18 – Armstrong, Wright & Jeffries (2018); B02 – Belikov et al. (2002); B20 – Berlanas et al. (2020); B94 – Brown, de Geus & de Zeeuw (1994); DZ99 – de Zeeuw et al. (1999); H97 – Hillenbrand (1997); and MD17 – Mel'nik & Dambis (2017).

For $d = 140$ pc and $D = 4 \times 10^{28}$ cm$^2$ s$^{-1}$, we get $t_{max} = 25$ kyr, which is much shorter than $\tau_{loss}$ for both $^{26}$Al, $^{60}$Fe and $^{56}$Fe. But for $d = 1.5$ kpc, $t_{max}$ is comparable to the radioactive lifetime of $^{26}$Al and $^{60}$Fe, which thus have time to decay before reaching the Solar system.

## 4 NEARBY OB ASSOCIATIONS

The nearest OB associations have been identified and studied for a long time (e.g. Blaauw 1964). The catalogue from de Zeeuw et al. (1999) based on Hipparcos positions, proper motions, and parallaxes provides a census of the stellar content of the OB associations within 1 kpc from the Sun. With improved astrometry, Gaia allows a better determination of the membership of stars belonging to OB associations, and the identification of new subgroups (Zucker et al. 2023). Recent compilations of O and B stars (Pantaleoni González et al. 2021) and OB associations (Wright 2020) make use of Gaia's results.

In this work, we consider all the well-studied OB associations listed in Wright (2020) and all the high-confidence OB associations at less than 1 kpc. Properties of these OB associations are summarized in Table 1 where their distance and age come from the review of Wright (2020). The numbers of observed stars mainly come from the catalogues of de Zeeuw et al. (1999) and Mel'nik & Dambis (2017) except for a few OB associations for which the star census

has been extensively studied such as Orion OB1 (Brown, de Geus & de Zeeuw 1994; Hillenbrand 1997), Perseus OB2 (Belikov et al. 2002), and Vela OB2 (Armstrong, Wright & Jeffries 2018).

We compute in the sixth column of Table 1 the richness of the OB association (or subgroup) that we define as the number of massive stars (≥ 8 M$_{☉}$) present when the association is formed. It is estimated considering the number of observed stars and the age of the OB association. For each realization of our population synthesis model an OB association is first given an age obtained by uniformly sampling the range of adopted ages (Table 1, third column). We assume 50 per cent uncertainty on the age when this is not specified. In a second step, the IMF (from Kroupa 2001) is sampled until the number of observed stars is reproduced taking into account the star lifetime from Limongi & Chieffi (2018). The number of massive stars is then recorded for each realization; the richness and associated standard deviation are obtained for typically 4000 realizations.

The determination of the richness depends on the mass range associated with the number of observed stars. However, this is only reported for very few OB associations: Orion OB1, Perseus OB1 and OB3, and Vela OB2. When the number of O and B stars is given instead, we use for the latest B-type stars a mass of 2.8 M$_{☉}$ obtained from a study of binary systems (Habets & Heintze 1981). A similar value is obtained using the evolutionary tracks from Palla & Stahler (1999) for pre-main-sequence models as shown in Preibisch
<table><caption>**Table 1.** Properties of local OB associations considered in this work. Distances and ages come from Wright (2020). The richness ($N_{*} \ge 8 M_{\odot}$ at formation time) and the number of core-collapse supernovae having already exploded ($N_{ccSN}$) are calculated from the number of observed stars and the OB association age (see text). $N_{ccSN}$ is reported for two different explodability criteria.</caption><thead><tr><th>Association</th><th>Distance (pc)</th><th>Age (Myr)</th><th>Number of observed stars</th><th>Ref.</th><th>N<sub>*</sub> (≥ 8 M<sub>☉</sub>)</th><th>N<sub>ccSN</sub><sup>LC18</sup></th><th>N<sub>ccSN</sub><sup>Sukhbold</sup></th></tr></thead><tbody><tr><td>Sco-Cen: US</td><td>143 ± 6</td><td>10 ± 7</td><td>N(OB) = 49</td><td>DZ99</td><td>13.9 ± 3.5</td><td>1.8 ± 1.8</td><td>1.6 ± 1.4</td></tr><tr><td>Sco-Cen: UCL</td><td>136 ± 5</td><td>16 ± 7</td><td>N(OB) = 66</td><td>DZ99</td><td>19.4 ± 4.3</td><td>5.1 ± 3.2</td><td>4.2 ± 2.7</td></tr><tr><td>Sco-Cen: LCC</td><td>115 ± 4</td><td>15 ± 6</td><td>N(OB) = 42</td><td>DZ99</td><td>12.4 ± 3.3</td><td>3.1 ± 2.1</td><td>2.5 ± 1.7</td></tr><tr><td>Ori OB1a</td><td>~360</td><td>8–12</td><td>N(4–15 M<sub>☉</sub>) = 53</td><td>B94</td><td>26.5 ± 4.9</td><td>2.2 ± 1.6</td><td>2.2 ± 1.4</td></tr><tr><td>Ori OB1b</td><td>360–420</td><td>2–8</td><td>N(4–120 M<sub>☉</sub>) = 45</td><td>B94</td><td>19.4 ± 3.6</td><td>0.5 ± 0.1</td><td>0.8 ± 0.6</td></tr><tr><td>Ori OB1c</td><td>~385</td><td>2–6</td><td>N(7–36 M<sub>☉</sub>) = 23</td><td>B94</td><td>22.7 ± 2.5</td><td>0</td><td>0.8 ± 0.6</td></tr><tr><td>Ori OB1d</td><td>~380</td><td>1–2</td><td>N(&gt; 1 M<sub>☉</sub>) = 145</td><td>H97</td><td>10.0 ± 3.0</td><td>0</td><td>0</td></tr><tr><td>Per OB1</td><td>~1830</td><td>8–11</td><td>N(OB) = 133</td><td>MD17</td><td>36.1 ± 5.6</td><td>2.2 ± 1.6</td><td>2.6 ± 1.1</td></tr><tr><td>Per OB2</td><td>296 ± 17</td><td>1–10</td><td>N(1–17 M<sub>☉</sub>) = 800</td><td>B02</td><td>53.9 ± 7.2</td><td>0.9 ± 1.1</td><td>1.8 ± 1.5</td></tr><tr><td>Per OB3</td><td>175 ± 3</td><td>50</td><td>N(OB) = 30</td><td>DZ99</td><td>10.8 ± 3.7</td><td>7.9 ± 3.4</td><td>7.1 ± 3.2</td></tr><tr><td>Cyg OB2</td><td>1350–1750</td><td>1–7</td><td>N(O) = 78</td><td>B20</td><td>241.2 ± 47.8</td><td>0</td><td>3.8 ± 3.4</td></tr><tr><td>Cyg OB4</td><td>~800</td><td>~8.3</td><td>N(OB) = 2</td><td>MD17</td><td>1.3 ± 0.8</td><td>0.5 ± 0.2</td><td>0.5 ± 0.2</td></tr><tr><td>Cyg OB7</td><td>~630</td><td>1–13</td><td>N(OB) = 25</td><td>MD17</td><td>7.2 ± 2.4</td><td>0.7 ± 0.6</td><td>0.8 ± 0.6</td></tr><tr><td>Cyg OB9</td><td>~960</td><td>2–4</td><td>N(OB) = 31</td><td>MD17</td><td>8.5 ± 2.5</td><td>0</td><td>0.5 ± 0.2</td></tr><tr><td>Vel OB2</td><td>411 ± 12</td><td>10–30</td><td>N(&gt; 2.5 M<sub>☉</sub>) = 72</td><td>A18</td><td>18.4 ± 4.4</td><td>6.8 ± 3.9</td><td>5.6 ± 3.5</td></tr><tr><td>Trumpler 10</td><td>372 ± 23</td><td>45–50</td><td>N(OB) = 22</td><td>DZ99</td><td>8.3 ± 3.2</td><td>6.7 ± 2.8</td><td>6.1 ± 2.6</td></tr><tr><td>Cas-Tau</td><td>125–300</td><td>~50</td><td>N(OB) = 83</td><td>DZ99</td><td>29.4 ± 6.6</td><td>22.1 ± 6.5</td><td>20.0 ± 6.2</td></tr><tr><td>Lac OB1</td><td>368 ± 17</td><td>2–25</td><td>N(OB) = 36</td><td>DZ99</td><td>10.7 ± 3.1</td><td>2.5 ± 2.3</td><td>2.2 ± 2.0</td></tr><tr><td>Cep OB2</td><td>~730</td><td>5</td><td>N(OB) = 56</td><td>DZ99</td><td>15.5 ± 3.4</td><td>0</td><td>0.8 ± 0.5</td></tr><tr><td>Cep OB3</td><td>~700</td><td>5–8</td><td>N(OB) = 25</td><td>MD17</td><td>7.2 ± 2.4</td><td>0.5 ± 0.1</td><td>0.7 ± 0.4</td></tr><tr><td>Cep OB4</td><td>~660</td><td>1–6</td><td>N(OB) = 7</td><td>MD17</td><td>2.5 ± 1.2</td><td>0</td><td>0.5 ± 0.2</td></tr><tr><td>Cep OB6</td><td>270 ± 12</td><td>~50</td><td>N(OB) = 6</td><td>DZ99</td><td>2.9 ± 1.8</td><td>2.3 ± 1.5</td><td>2.1 ± 1.4</td></tr><tr><td>Collinder 121</td><td>543 ± 23</td><td>5</td><td>N(OB) = 87</td><td>DZ99</td><td>22.9 ± 4.1</td><td>0</td><td>0.9 ± 8.8</td></tr><tr><td>Cam OB1</td><td>~800</td><td>7–14</td><td>N(OB) = 45</td><td>MD17</td><td>12.9 ± 3.3</td><td>1.6 ± 1.3</td><td>1.5 ± 1.1</td></tr><tr><td>Mon OB1</td><td>~580</td><td>1–10</td><td>N(OB) = 6</td><td>MD17</td><td>2.3 ± 1.2</td><td>0.5 ± 0.1</td><td>0.5 ± 0.2</td></tr></tbody></table>

Note. References: A18 – Armstrong, Wright & Jeffries (2018); B02 – Belikov et al. (2002); B20 – Berlanas et al. (2020); B94 – Brown, de Geus & de Zeeuw (1994); DZ99 – de Zeeuw et al. (1999); H97 – Hillenbrand (1997); and MD17 – Mel'nik & Dambis (2017).

For $d = 140$ pc and $D = 4 \times 10^{28}$ cm$^2$ s$^{-1}$, we get $t_{max} = 25$ kyr, which is much shorter than $\tau_{loss}$ for both $^{26}$Al, $^{60}$Fe and $^{56}$Fe. But for $d = 1.5$ kpc, $t_{max}$ is comparable to the radioactive lifetime of $^{26}$Al and $^{60}$Fe, which thus have time to decay before reaching the Solar system.

## 4 NEARBY OB ASSOCIATIONS

The nearest OB associations have been identified and studied for a long time (e.g. Blaauw 1964). The catalogue from de Zeeuw et al. (1999) based on Hipparcos positions, proper motions, and parallaxes provides a census of the stellar content of the OB associations within 1 kpc from the Sun. With improved astrometry, Gaia allows a better determination of the membership of stars belonging to OB associations, and the identification of new subgroups (Zucker et al. 2023). Recent compilations of O and B stars (Pantaleoni González et al. 2021) and OB associations (Wright 2020) make use of Gaia's results.

In this work, we consider all the well-studied OB associations listed in Wright (2020) and all the high-confidence OB associations at less than 1 kpc. Properties of these OB associations are summarized in Table 1 where their distance and age come from the review of Wright (2020). The numbers of observed stars mainly come from the catalogues of de Zeeuw et al. (1999) and Mel'nik & Dambis (2017) except for a few OB associations for which the star census

has been extensively studied such as Orion OB1 (Brown, de Geus & de Zeeuw 1994; Hillenbrand 1997), Perseus OB2 (Belikov et al. 2002), and Vela OB2 (Armstrong, Wright & Jeffries 2018).

We compute in the sixth column of Table 1 the richness of the OB association (or subgroup) that we define as the number of massive stars ($\ge 8 M_{\odot}$) present when the association is formed. It is estimated considering the number of observed stars and the age of the OB association. For each realization of our population synthesis model an OB association is first given an age obtained by uniformly sampling the range of adopted ages (Table 1, third column). We assume 50 per cent uncertainty on the age when this is not specified. In a second step, the IMF (from Kroupa 2001) is sampled until the number of observed stars is reproduced taking into account the star lifetime from Limongi & Chieffi (2018). The number of massive stars is then recorded for each realization; the richness and associated standard deviation are obtained for typically 4000 realizations.

The determination of the richness depends on the mass range associated with the number of observed stars. However, this is only reported for very few OB associations: Orion OB1, Perseus OB1 and OB3, and Vela OB2. When the number of O and B stars is given instead, we use for the latest B-type stars a mass of $2.8 M_{\odot}$ obtained from a study of binary systems (Habets & Heintze 1981). A similar value is obtained using the evolutionary tracks from Palla & Stahler (1999) for pre-main-sequence models as shown in Preibisch
<figure><img src="image_6.png" /><figcaption>**Figure 5.** 3D representation of the local OB associations in the solar neighbourhood. The volume of each sphere is proportional to the richness of the association. the colour accounts for the typical age of the association (see Table 1). The star symbol indicates the position of the Sun at (X,Y,Z) = (0,0,0) and the arrow points towards the Galactic Centre.</figcaption></figure>

et al. (2002). In these conditions, the total number of stars we obtain for Upper Scorpius when normalizing the IMF from Preibisch et al. (2002) to the 49 B stars reported in de Zeeuw et al. (1999) is 2590, in very good agreement with the 2525 stars reported by Preibisch et al. (2002). For the latest O-type stars we consider that they have masses of 16 $M_\odot$ or larger (Habets & Heintze 1981; Weidner & Vink 2010).

The past nucleosynthetic activity of an OB association is related to the number of ccSN that have exploded so far ($N_{ccSN}$). In the same calculation as for the richness, the number of exploding massive stars with a stellar lifetime shorter than the age of the OB association is recorded. We use by default the explodability criterion from Limongi & Chieffi (2018), i.e. $M < 25 M_\odot$, and the corresponding $N_{ccSN}$ is reported in the seventh column of Table 1. For young OB associations with ages smaller than the first ccSN explosion time (occurring at about 7.8 Myr, and corresponding to the lifetime of a non-rotating 25 $M_\odot$ star) no ccSN has exploded yet. In these cases the enrichment of the gas of the OB association mainly comes from stellar winds that cannot be accelerated as CRs because no SN exploded yet. Hence, the OB association is not expected to contribute to the CR density budget even though it may be a high richness association (e.g. Cyg OB2, Orion OB1c, and Collinder 121). On the contrary, for rather old OB associations with ages greater than the last massive star explosion time (occurring at about 40 Myr, and corresponding to the lifetime of a non-rotating 8 $M_\odot$ star) all massive stars may have exploded. However, even in the case of a high richness OB association (e.g. Cas-Tau), the present-day enrichment in short-lived radionuclides (e.g. $^{26}$Al and $^{60}$Fe) of the associated superbubble gas will most likely be negligible because of the smaller yields for the low-end massive stars and the free radionuclide decay after the last massive star explosion (see Fig. 3).

The number of past ccSN for a given OB association significantly depends on the explodability criterion that is considered. In the last column of Table 1, we compute $N_{ccSN}$ using Sukhbold et al. (2016) explodability criterion. In this case, and at variance with the case using Limongi & Chieffi (2018) explodability criterion, some stars having initial masses greater than 25 $M_\odot$ explode as ccSN. These stars have lifetime smaller than 7.8 Myr, so younger OB associations will have a nucleosynthetic activity while this is not the case with Limongi & Chieffi (2018) explodability criterion (e.g. Orion OB1c). On the other hand, for older OB associations the nucleosynthetic activity may be reduced when considering the Sukhbold et al. (2016) explosion criterion (e.g. Sco–Cen) since some stars in the low-end of the massive range ($< 25 M_\odot$) may not explode as ccSN.

The number of past ccSN we obtain for Sco–Cen is 10 ± 7.1 or 8.3 ± 5.8 depending on the explodability criterion, which is in reasonable agreement with the number of past SNe, between 14 and 20, needed to excavate the Local Bubble (Fuchs et al. 2006; Breitschwerdt et al. 2016).

A 3D representation of the OB associations in our solar neighbourhood is presented in Fig. 5 with the volume of each OB association proportional to its richness while the age of the association is colour coded.

# 5 ORIGIN OF $^{56,60}$FE AND $^{26}$AL IN CRS

## 5.1 CR density distribution and observations

Our CR population synthesis model was used to compute the LISM CR density of $^{60}$Fe, $^{56}$Fe, and $^{26}$Al resulting from the contribution of all OB associations listed in Table 1. We take as nominal parameters for the CRs acceleration and propagation a mean
<figure><img src="image_6.png" /><figcaption>**Figure 5.** 3D representation of the local OB associations in the solar neighbourhood. The volume of each sphere is proportional to the richness of the association. the colour accounts for the typical age of the association (see Table 1). The star symbol indicates the position of the Sun at (X,Y,Z) = (0,0,0) and the arrow points towards the Galactic Centre.</figcaption></figure>

et al. (2002). In these conditions, the total number of stars we obtain for Upper Scorpius when normalizing the IMF from Preibisch et al. (2002) to the 49 B stars reported in de Zeeuw et al. (1999) is 2590, in very good agreement with the 2525 stars reported by Preibisch et al. (2002). For the latest O-type stars we consider that they have masses of 16 $M_\odot$ or larger (Habets & Heintze 1981; Weidner & Vink 2010).

The past nucleosynthetic activity of an OB association is related to the number of ccSN that have exploded so far ($N_{\text{ccSN}}$). In the same calculation as for the richness, the number of exploding massive stars with a stellar lifetime shorter than the age of the OB association is recorded. We use by default the explodability criterion from Limongi & Chieffi (2018), i.e. $M < 25 M_\odot$, and the corresponding $N_{\text{ccSN}}$ is reported in the seventh column of Table 1. For young OB associations with ages smaller than the first ccSN explosion time (occurring at about 7.8 Myr, and corresponding to the lifetime of a non-rotating 25 $M_\odot$ star) no ccSN has exploded yet. In these cases the enrichment of the gas of the OB association mainly comes from stellar winds that cannot be accelerated as CRs because no SN exploded yet. Hence, the OB association is not expected to contribute to the CR density budget even though it may be a high richness association (e.g. Cyg OB2, Orion OB1c, and Collinder 121). On the contrary, for rather old OB associations with ages greater than the last massive star explosion time (occurring at about 40 Myr, and corresponding to the lifetime of a non-rotating 8 $M_\odot$ star) all massive stars may have exploded. However, even in the case of a high richness OB association (e.g. Cas-Tau), the present-day enrichment in short-lived radionuclides (e.g. $^{26}$Al and $^{60}$Fe) of the associated superbubble gas will most likely be negligible because of the smaller yields for the low-end massive stars and the free radionuclide decay after the last massive star explosion (see Fig. 3).

The number of past ccSN for a given OB association significantly depends on the explodability criterion that is considered. In the last column of Table 1, we compute $N_{\text{ccSN}}$ using Sukhbold et al. (2016) explodability criterion. In this case, and at variance with the case using Limongi & Chieffi (2018) explodability criterion, some stars having initial masses greater than 25 $M_\odot$ explode as ccSN. These stars have lifetime smaller than 7.8 Myr, so younger OB associations will have a nucleosynthetic activity while this is not the case with Limongi & Chieffi (2018) explodability criterion (e.g. Orion OB1c). On the other hand, for older OB associations the nucleosynthetic activity may be reduced when considering the Sukhbold et al. (2016) explosion criterion (e.g. Sco–Cen) since some stars in the low-end of the massive range ($< 25 M_\odot$) may not explode as ccSN.

The number of past ccSN we obtain for Sco–Cen is 10 $\pm$ 7.1 or 8.3 $\pm$ 5.8 depending on the explodability criterion, which is in reasonable agreement with the number of past SNe, between 14 and 20, needed to excavate the Local Bubble (Fuchs et al. 2006; Breitschwerdt et al. 2016).

A 3D representation of the OB associations in our solar neighbourhood is presented in Fig. 5 with the volume of each OB association proportional to its richness while the age of the association is colour coded.

# 5 ORIGIN OF $^{56,60}$FE AND $^{26}$AL IN CRS

## 5.1 CR density distribution and observations

Our CR population synthesis model was used to compute the LISM CR density of $^{60}$Fe, $^{56}$Fe, and $^{26}$Al resulting from the contribution of all OB associations listed in Table 1. We take as nominal parameters for the CRs acceleration and propagation a mean
acceleration efficiency $\epsilon_{acc} = 10^{-5}$, a mean diffusion coefficient $D_0 = 3.08 \times 10^{28} \text{ cm}^2 \text{ s}^{-1}$ (equation 3), and a mean ISM average density $n_H = 0.1 \text{ cm}^{-3}$. The impact of these parameters will be discussed in Section 6.1. A realization of our CR population synthesis model is defined as the sampling of the IMF until the richness of each considered OB association is reproduced. For each realization, the distance of an OB association needed to calculate the CR densities is obtained by uniformly sampling the adopted distances (Table 1, second column). We assume a 15 per cent error on the distance when the uncertainty is not specified. The total CR density distribution for a given nuclide, obtained as the sum of the contribution of each OB association, is shown in Fig. 6 (blue histogram) for 4000 realizations, and compared to the ACE/CRIS measurements (hatched and solid vertical red lines). The median of the total $^{60}$Fe CR density distribution is indicated with the solid brown vertical lines and the 16th and 84th percentiles, defining a 68 per cent probability coverage, are the dashed brown vertical lines.

The calculated $^{60}$Fe CR density has a rather broad distribution with a mean of $\approx 2.5 \times 10^{-19}$ atoms cm$^{-3}$. This is about 2.8 times smaller than what is deduced from the ACE/CRIS observations (Binns et al. 2016). However, the observations are well within the calculated distribution at slightly more than $1\sigma$ from the median. This indicates that the observed density of $^{60}$Fe CRs in the LISM is not exceptional, indeed it represents $\approx 8$ per cent of the simulated cases. The spread of the distribution arises from the stochastic nature of the IMF, the different $^{60}$Fe yields as a function of the stellar initial mass, the competition between the $^{60}$Fe lifetime and the mean time between two successive ccSN, and the contribution from the different OB associations.

The calculated $^{56}$Fe CR density distribution (blue) is not as broad as in the case of $^{60}$Fe that is due to the similar $^{56}$Fe yield for each ccSN and the stable nature of $^{56}$Fe. The CR distribution for the realizations matching the $^{60}$Fe observations is also displayed as a brown histogram. On average, the calculated $^{56}$Fe CR density represents about 20 per cent of the observed value. This suggests that a non-negligible fraction of the $^{56}$Fe CR density in the LISM comes from local sources (see further discussion in Section 6.3).

Concerning the calculated CR density for $^{26}$Al it is lower by more than one order of magnitude than the ACE/CRIS observations. This is expected since $^{26}$Al is mostly produced by CRs spallation (Yanasak et al. 2001) and that our CR population synthesis model only computes the primary component of CRs. This result suggests that the excess in the Al CR spectrum found by Boschini et al. (2022) is not produced by a contribution of primary $^{26}$Al.

## 5.2 The role of Sco–Cen

For each realization of our CR population synthesis model it is interesting to know which OB association is contributing the most to the total CR density of $^{60}$Fe (shown in Fig. 6). This is what is represented by the blue histogram in Fig. 7 that indicates that in $\approx 70$ per cent of the cases Sco–Cen is the main contributor to the predicted total CR density of $^{60}$Fe, followed by Vela OB2, Orion OB1, and Cas-Tau at the 10 per cent level. However, this does not tell anything about whether, realization by realization, the most contributing OB association dominates largely the other ones or whether its contribution is more equally shared. The inset in Fig. 7 shows the distribution of the fraction of total CR density for the most contributing OB associations. The distributions are very different between Sco–Cen and the other associations. For Sco–Cen the distribution is peaked for large fractions meaning that when Sco–

Cen is the main contributing association this is by far the dominant one. Specifically, Sco–Cen contributes by more than 80 per cent to the total CR density in 64 per cent of the cases. For Vela OB2, Orion OB1, and Cas-Tau the fraction distributions are rather flat with a maximum at about 50–60 per cent indicating that the contribution to the total CR density is much more equally shared between the participating associations.

If we now only consider the realizations compatible with the ACE/CRIS observations, it appears that the 15 detected $^{60}$Fe nuclei are nearly always coming from the Sco–Cen association as shown by the red hatched histogram in Fig. 7. The configuration of these realizations is quite specific since they all involve at least one massive star having exploded recently. Fig. 8 shows the explosion time distribution for the realizations where a single SN in Sco–Cen accelerates more than 50 per cent of the total $^{60}$Fe CR density. The mean explosion time is 146 kyr with an rms of 90 kyr, and in 93 per cent of the cases the age of the SN is smaller than 300 kyr. The initial mass distribution of these SNe at the origin of the acceleration of the $^{60}$Fe present in the enriched gas of the OB association is represented in the inset of Fig. 8. It is characterized by a mean stellar initial mass of 16.1 $\text{M}_\odot$ with an rms of 3.2 $\text{M}_\odot$.

# 6 DISCUSSION

## 6.1 Sensitivity to input parameters

The results presented so far are obtained with a nominal set of input parameters. However, the underlying physics of both the population synthesis and the acceleration and propagation of CRs is far from being under control. It is thus important to explore the robustness of our findings by varying the main ingredients of our CR population synthesis model within uncertainties. Concerning the population synthesis part of our model, we first investigate a case where the IMF is determined from the observed Upper Scorpius OB association members (Preibisch et al. 2002) rather than from the model developed by Kroupa (2001). We also investigate a case with the explodability criterion taken from Sukhbold et al. (2016) and not from Limongi & Chieffi (2018). Finally, we perform a simulation with the yields obtained from the PUSH model (Curtis et al. 2019) based on the pre-explosion models of Woosley & Heger (2007) for non-rotating solar-metallicity stars.

Concerning the acceleration and propagation of CRs, we explore four different cases by varying the relevant input parameters. First, we consider a case where the diffusion coefficient is reduced by an order of magnitude, i.e. taking $D_0 = 3.08 \times 10^{27} \text{ cm}^2 \text{ s}^{-1}$ in equation (3), which may be expected if CRs spent most of their time in an active superbubble environment (Section 3.2.2). We also consider a case where the average density of the ambient medium is increased by an order of magnitude, i.e. $n_H = 1 \text{ cm}^{-3}$, which may be typical if CRs diffuse for a significant time in superbubble shells and/or in the ISM of the Galactic disc outside superbubbles. We also study the impact of increasing the CR acceleration efficiency by an order of magnitude, i.e. $\epsilon_{acc} = 10^{-4}$, to take into account that refractory elements such as Fe may be more efficiently accelerated by the DSA process than volatile elements (Section 3.2.1). Finally, we compute a case where $D_0$, $n_H$, and $\epsilon_{acc}$ are independently determined for each model realization from a lognormal distribution with a factor uncertainty of 2. This case is intended to take into account that the acceleration of CRs and their propagation to the Solar system could depend on the individual properties of the nearby OB associations and their specific location in the ISM. In particular, $\epsilon_{acc}$ could depend on the size and age of the parent OB association, whereas $D_0$ and $n_H$
acceleration efficiency $\epsilon_{acc} = 10^{-5}$, a mean diffusion coefficient $D_0 = 3.08 \times 10^{28} \text{ cm}^2 \text{ s}^{-1}$ (equation 3), and a mean ISM average density $n_H = 0.1 \text{ cm}^{-3}$. The impact of these parameters will be discussed in Section 6.1. A realization of our CR population synthesis model is defined as the sampling of the IMF until the richness of each considered OB association is reproduced. For each realization, the distance of an OB association needed to calculate the CR densities is obtained by uniformly sampling the adopted distances (Table 1, second column). We assume a 15 per cent error on the distance when the uncertainty is not specified. The total CR density distribution for a given nuclide, obtained as the sum of the contribution of each OB association, is shown in Fig. 6 (blue histogram) for 4000 realizations, and compared to the ACE/CRIS measurements (hatched and solid vertical red lines). The median of the total $^{60}$Fe CR density distribution is indicated with the solid brown vertical lines and the 16th and 84th percentiles, defining a 68 per cent probability coverage, are the dashed brown vertical lines.

The calculated $^{60}$Fe CR density has a rather broad distribution with a mean of $\approx 2.5 \times 10^{-19}$ atoms cm$^{-3}$. This is about 2.8 times smaller than what is deduced from the ACE/CRIS observations (Binns et al. 2016). However, the observations are well within the calculated distribution at slightly more than $1\sigma$ from the median. This indicates that the observed density of $^{60}$Fe CRs in the LISM is not exceptional, indeed it represents $\approx 8$ per cent of the simulated cases. The spread of the distribution arises from the stochastic nature of the IMF, the different $^{60}$Fe yields as a function of the stellar initial mass, the competition between the $^{60}$Fe lifetime and the mean time between two successive ccSN, and the contribution from the different OB associations.

The calculated $^{56}$Fe CR density distribution (blue) is not as broad as in the case of $^{60}$Fe that is due to the similar $^{56}$Fe yield for each ccSN and the stable nature of $^{56}$Fe. The CR distribution for the realizations matching the $^{60}$Fe observations is also displayed as a brown histogram. On average, the calculated $^{56}$Fe CR density represents about 20 per cent of the observed value. This suggests that a non-negligible fraction of the $^{56}$Fe CR density in the LISM comes from local sources (see further discussion in Section 6.3).

Concerning the calculated CR density for $^{26}$Al it is lower by more than one order of magnitude than the ACE/CRIS observations. This is expected since $^{26}$Al is mostly produced by CRs spallation (Yanasak et al. 2001) and that our CR population synthesis model only computes the primary component of CRs. This result suggests that the excess in the Al CR spectrum found by Boschini et al. (2022) is not produced by a contribution of primary $^{26}$Al.

## 5.2 The role of Sco–Cen

For each realization of our CR population synthesis model it is interesting to know which OB association is contributing the most to the total CR density of $^{60}$Fe (shown in Fig. 6). This is what is represented by the blue histogram in Fig. 7 that indicates that in $\approx 70$ per cent of the cases Sco–Cen is the main contributor to the predicted total CR density of $^{60}$Fe, followed by Vela OB2, Orion OB1, and Cas-Tau at the 10 per cent level. However, this does not tell anything about whether, realization by realization, the most contributing OB association dominates largely the other ones or whether its contribution is more equally shared. The inset in Fig. 7 shows the distribution of the fraction of total CR density for the most contributing OB associations. The distributions are very different between Sco–Cen and the other associations. For Sco–Cen the distribution is peaked for large fractions meaning that when Sco–

Cen is the main contributing association this is by far the dominant one. Specifically, Sco–Cen contributes by more than 80 per cent to the total CR density in 64 per cent of the cases. For Vela OB2, Orion OB1, and Cas-Tau the fraction distributions are rather flat with a maximum at about 50–60 per cent indicating that the contribution to the total CR density is much more equally shared between the participating associations.

If we now only consider the realizations compatible with the ACE/CRIS observations, it appears that the 15 detected $^{60}$Fe nuclei are nearly always coming from the Sco–Cen association as shown by the red hatched histogram in Fig. 7. The configuration of these realizations is quite specific since they all involve at least one massive star having exploded recently. Fig. 8 shows the explosion time distribution for the realizations where a single SN in Sco–Cen accelerates more than 50 per cent of the total $^{60}$Fe CR density. The mean explosion time is 146 kyr with an rms of 90 kyr, and in 93 per cent of the cases the age of the SN is smaller than 300 kyr. The initial mass distribution of these SNe at the origin of the acceleration of the $^{60}$Fe present in the enriched gas of the OB association is represented in the inset of Fig. 8. It is characterized by a mean stellar initial mass of 16.1 $\text{M}_\odot$ with an rms of 3.2 $\text{M}_\odot$.

# 6 DISCUSSION

## 6.1 Sensitivity to input parameters

The results presented so far are obtained with a nominal set of input parameters. However, the underlying physics of both the population synthesis and the acceleration and propagation of CRs is far from being under control. It is thus important to explore the robustness of our findings by varying the main ingredients of our CR population synthesis model within uncertainties. Concerning the population synthesis part of our model, we first investigate a case where the IMF is determined from the observed Upper Scorpius OB association members (Preibisch et al. 2002) rather than from the model developed by Kroupa (2001). We also investigate a case with the explodability criterion taken from Sukhbold et al. (2016) and not from Limongi & Chieffi (2018). Finally, we perform a simulation with the yields obtained from the PUSH model (Curtis et al. 2019) based on the pre-explosion models of Woosley & Heger (2007) for non-rotating solar-metallicity stars.

Concerning the acceleration and propagation of CRs, we explore four different cases by varying the relevant input parameters. First, we consider a case where the diffusion coefficient is reduced by an order of magnitude, i.e. taking $D_0 = 3.08 \times 10^{27} \text{ cm}^2 \text{ s}^{-1}$ in equation (3), which may be expected if CRs spent most of their time in an active superbubble environment (Section 3.2.2). We also consider a case where the average density of the ambient medium is increased by an order of magnitude, i.e. $n_H = 1 \text{ cm}^{-3}$, which may be typical if CRs diffuse for a significant time in superbubble shells and/or in the ISM of the Galactic disc outside superbubbles. We also study the impact of increasing the CR acceleration efficiency by an order of magnitude, i.e. $\epsilon_{acc} = 10^{-4}$, to take into account that refractory elements such as Fe may be more efficiently accelerated by the DSA process than volatile elements (Section 3.2.1). Finally, we compute a case where $D_0$, $n_H$, and $\epsilon_{acc}$ are independently determined for each model realization from a lognormal distribution with a factor uncertainty of 2. This case is intended to take into account that the acceleration of CRs and their propagation to the Solar system could depend on the individual properties of the nearby OB associations and their specific location in the ISM. In particular, $\epsilon_{acc}$ could depend on the size and age of the parent OB association, whereas $D_0$ and $n_H$
<figure><img src="image_7.png" /></figure>

**Figure 6.** Calculated CR density distribution in the LISM (blue histograms) for *N* = 4000 realizations of our model including all OB associations listed in Table 1 for <sup>60</sup>Fe (left panel), <sup>56</sup>Fe (middle panel), and <sup>26</sup>Al (right panel). Vertical brown thick lines correspond to the 16th, 50th, and 84th percentile of the <sup>60</sup>Fe CR density distribution. Hatched and solid vertical red lines correspond to the CR densities derived from the ACE/CRIS measurements. Brown histograms represent cases for which the predicted <sup>60</sup>Fe CR density matches the observations. Nuclides in the gas of the OB associations are injected in CRs with an efficiency ε<sub>acc</sub> = 10<sup>-5</sup>. CRs are assumed to propagate in the ISM of average density n<sub>H</sub> = 0.1 cm<sup>-3</sup> with a mean diffusion coefficient D<sub>0</sub> = 3.08 × 10<sup>28</sup> cm<sup>2</sup> s<sup>-1</sup> (equation 3).

<figure><img src="image_8.png" /></figure>

**Figure 7.** Occurrence of the most contributing OB association to the total CR density of <sup>60</sup>Fe in the LISM (solid histogram) for *N* = 4000 realizations of our CR population synthesis model. The hatched histogram is the same but for realizations matching the <sup>60</sup>Fe ACE/CRIS observations. The distribution of the fraction of total CR density for the main contributing OB associations is given in the inset; the distribution for OB associations other than Sco–Cen is displayed with a scaling factor of 5.

<figure><img src="image_9.png" /></figure>

**Figure 8.** Age and stellar initial mass distributions for the Sco–Cen SN accelerating more than 50 per cent of the observed <sup>60</sup>Fe CR density measured by ACE/CRIS.

<figure><img src="image_10.png" /></figure>

**Figure 9.** Effect of the input parameters of our CR population synthesis model when we consider the realizations matching the <sup>60</sup>Fe observations from ACE/CRIS. Top: Occurrence of Sco–Cen as the most contributing OB association (filled squares), and the ratio of the calculated mean <sup>56</sup>Fe CR density to the observed one (open star symbols). Bottom: Age (down-pointing triangles) and mass (up-pointing triangles) of the Sco–Cen SN accelerating more than 50 per cent of the observed <sup>60</sup>Fe CR density.

could be related to the distance of the source and its position with respect to the magnetic field lines passing near the Solar system. Results are gathered in Fig. 9 for the realizations matching the LISM <sup>60</sup>Fe CR density determined from the ACE/CRIS observations.

The red full squares represent the probability that the observed LISM <sup>60</sup>Fe CR density can be explained by the contribution of
<figure><img src="image_7.png" /></figure>

**Figure 6.** Calculated CR density distribution in the LISM (blue histograms) for *N* = 4000 realizations of our model including all OB associations listed in Table 1 for <sup>60</sup>Fe (left panel), <sup>56</sup>Fe (middle panel), and <sup>26</sup>Al (right panel). Vertical brown thick lines correspond to the 16th, 50th, and 84th percentile of the <sup>60</sup>Fe CR density distribution. Hatched and solid vertical red lines correspond to the CR densities derived from the ACE/CRIS measurements. Brown histograms represent cases for which the predicted <sup>60</sup>Fe CR density matches the observations. Nuclides in the gas of the OB associations are injected in CRs with an efficiency ε<sub>acc</sub> = 10<sup>-5</sup>. CRs are assumed to propagate in the ISM of average density n<sub>H</sub> = 0.1 cm<sup>-3</sup> with a mean diffusion coefficient D<sub>0</sub> = 3.08 × 10<sup>28</sup> cm<sup>2</sup> s<sup>-1</sup> (equation 3).

<figure><img src="image_8.png" /></figure>

**Figure 7.** Occurrence of the most contributing OB association to the total CR density of <sup>60</sup>Fe in the LISM (solid histogram) for *N* = 4000 realizations of our CR population synthesis model. The hatched histogram is the same but for realizations matching the <sup>60</sup>Fe ACE/CRIS observations. The distribution of the fraction of total CR density for the main contributing OB associations is given in the inset; the distribution for OB associations other than Sco–Cen is displayed with a scaling factor of 5.

<figure><img src="image_9.png" /></figure>

**Figure 8.** Age and stellar initial mass distributions for the Sco–Cen SN accelerating more than 50 per cent of the observed <sup>60</sup>Fe CR density measured by ACE/CRIS.

<figure><img src="image_10.png" /></figure>

**Figure 9.** Effect of the input parameters of our CR population synthesis model when we consider the realizations matching the <sup>60</sup>Fe observations from ACE/CRIS. Top: Occurrence of Sco–Cen as the most contributing OB association (filled squares), and the ratio of the calculated mean <sup>56</sup>Fe CR density to the observed one (open star symbols). Bottom: Age (down-pointing triangles) and mass (up-pointing triangles) of the Sco–Cen SN accelerating more than 50 per cent of the observed <sup>60</sup>Fe CR density.

could be related to the distance of the source and its position with respect to the magnetic field lines passing near the Solar system. Results are gathered in Fig. 9 for the realizations matching the LISM <sup>60</sup>Fe CR density determined from the ACE/CRIS observations.

The red full squares represent the probability that the observed LISM <sup>60</sup>Fe CR density can be explained by the contribution of
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Same as Fig. 7 except that the explodability criterion is from Sukhbold et al. (2016).</figcaption></figure>

the Sco–Cen OB association only. We see that the results are largely independent of the assumed input parameters except for the acceleration efficiency and the explodability criterion. In the latter case the Orion OB1 and Perseus OB2 associations are also able, on their own, to produce the observed $^{60}$Fe CR density (see Fig. 10). This can be semiquantitatively explained by comparing the number of past SNe ($N_{\text{ccSN}}$) reported in the last two columns of Table 1. Indeed, the ratio $N_{\text{ccSN}}(\text{Sco}-\text{Cen})/N_{\text{ccSN}}(\text{Orion OB1})$ decreases from 10.0/2.7 = 3.7 to 8.3/3.8 = 2.2 when considering the Limongi & Chieffi (2018) or Sukhbold et al. (2016) explodability criterion, respectively. This indicates that the relative contribution of Orion OB1 should increase with respect to Sco–Cen when the Sukhbold et al. (2016) explosion criterion is considered, as observed in Fig. 10. In a similar way, the number of past SNe is higher (lower) for Perseus OB2 (Vela OB2) when considering the Sukhbold et al. (2016) explodability criterion, which consequently increases (reduces) the importance of these OB associations with respect to the reference case using the Limongi & Chieffi (2018) explodability criterion.

In the case where the acceleration efficiency is multiplied by a factor of 10, the Sco–Cen association still dominates over the other associations, but with a reduced probability of about 60 per cent. Indeed, with a higher acceleration efficiency more material is accelerated resulting in an increased predicted $^{60}$Fe CR density. It follows that OB associations other than Sco–Cen are then also able to account for the observed $^{60}$Fe CR density. There are then a larger number of realizations where a smaller contribution from Sco–Cen due to the initial draw of the massive star population is compensated by a larger contribution from associations such as Orion OB1 or Perseus OB2. This consequently reduces the occurrence of Sco–Cen as the most contributing OB association.

The green open stars in Fig. 9 represent the ratio between the average value of the calculated $^{56}$Fe CR density for the realizations, which account for the $^{60}$Fe observations, and the observed $^{56}$Fe CR density value. This can be taken as an indicator of the $^{56}$Fe CR density fraction that may come from a local source. This fraction is typically 20 per cent except for the cases where the CR diffusion coefficient and acceleration efficiency are varied by a factor of 10; for such cases the fraction is about 50 per cent. Since the nucleosynthesis activity for these two cases is the same, the higher predicted $^{56}$Fe CR density is correlated to an increase of the $^{60}$Fe CR density. The reason for such an increase when the CR acceleration efficiency is higher has previously been discussed. When the CR diffusion coefficient is

decreased, the propagation time needed to reach the maximum CR density for a given OB association distance is higher (see equation 12). Thus, there are possibly more ccSN that can explode during this time lapse, leading then to a higher predicted $^{56,60}$Fe CR density.

Since in most cases Sco–Cen is the OB association contributing the most to the LISM $^{60}$Fe CR density, it is interesting to explore the impact of the input parameters on the properties of the SNe accelerating the material present in the gas of the superbubble. Fig. 9 (bottom) shows the age (black) and mass (blue) of the Sco–Cen ccSN accelerating more than 50 per cent of the observed $^{60}$Fe CR density. The triangle markers represent the median (50th percentile) of the distribution while the lower (upper) bound of the error bar corresponds to the 16th (84th) percentile, defining a 68 per cent coverage probability. We see that most of the Sco–Cen realizations host a young SN with an age smaller than 300 kyr, except for two cases for which the explosion time distribution of the SN extends to higher values. This arises from different effects depending on the case. For smaller values of the CR diffusion coefficient, as discussed previously, one expects the age of the SNe to be greater and to span a larger range since the propagation time for CRs is longer. In the case where the CR acceleration efficiency is increased, older SNe, which would have contributed negligibly otherwise, may now contribute significantly to the observed $^{60}$Fe CR density budget. Interestingly, the mass of the Sco–Cen ccSN accelerating more than 50 per cent of the observed $^{60}$Fe CR density is nearly independent of the different test cases that we explored: the median value is $\approx 15–16 M_{\odot}$ and the $1\sigma$ range is $\approx 13–20 M_{\odot}$.

## 6.2 Geminga

The Geminga pulsar is currently located in the constellation Gemini at a distance of about $157_{-34}^{+59}$ pc with a proper motion of $170 \pm 4$ mas yr$^{-1}$ (Caraveo et al. 1996). Its spin-down age deduced from the pulsar period and period derivative (Bignami & Caraveo 1992) is 342 kyr, which can be considered as representative of its true age (Pellizza et al. 2005). Even if these properties are well established, the place of birth of Geminga is not clearly identified yet. By tracing back the space motion of Geminga, Pellizza et al. (2005) find that Geminga was born at 90–240 pc from the Sun, most probably inside the Orion OB1a association or the Cas-Tau OB association. Moreover, these authors conclude that the Geminga progenitor mass should not be greater than $15 M_{\odot}$.

One of the main conclusions of this study is that the Sco–Cen OB association plays a specific role in explaining the observed LISM $^{60}$Fe CR density. However, our results also show (i) that the Cas-Tau OB association is able to reproduce the observations, even though this is much more unlikely, and (ii) that the occurrence of Orion OB1 as the most contributing association can be significant when the Sukhbold et al. (2016) explodability criterion is used (see Fig. 10).

Here, we investigate whether the Geminga progenitor could be the SN that accelerated the $^{60}$Fe nuclei observed by ACE/CRIS, and whether this SN could be associated with the Orion OB1 or Cas-Tau associations.

In the following, we consider Orion OB1 and Cas-Tau as two independent OB associations. We first compute the nucleosynthetic activity as a function of time of each OB association, and we estimate, for each realization of our model ($N = 4000$), the amount of $^{60}$Fe present in the associated superbubble gas when the Geminga progenitor exploded 342 kyr ago. In a second step, and for each realization, the distance of Geminga to the Solar system $d$ is uniformly sampled up to the distance of the considered OB associations. CRs are then accelerated and propagated across the distance $d$ during a time
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Same as Fig. 7 except that the explodability criterion is from Sukhbold et al. (2016).</figcaption></figure>

the Sco–Cen OB association only. We see that the results are largely independent of the assumed input parameters except for the acceleration efficiency and the explodability criterion. In the latter case the Orion OB1 and Perseus OB2 associations are also able, on their own, to produce the observed $^{60}$Fe CR density (see Fig. 10). This can be semiquantitatively explained by comparing the number of past SNe ($N_{\text{ccSN}}$) reported in the last two columns of Table 1. Indeed, the ratio $N_{\text{ccSN}}(\text{Sco}-\text{Cen})/N_{\text{ccSN}}(\text{Orion OB1})$ decreases from 10.0/2.7 = 3.7 to 8.3/3.8 = 2.2 when considering the Limongi & Chieffi (2018) or Sukhbold et al. (2016) explodability criterion, respectively. This indicates that the relative contribution of Orion OB1 should increase with respect to Sco–Cen when the Sukhbold et al. (2016) explosion criterion is considered, as observed in Fig. 10. In a similar way, the number of past SNe is higher (lower) for Perseus OB2 (Vela OB2) when considering the Sukhbold et al. (2016) explodability criterion, which consequently increases (reduces) the importance of these OB associations with respect to the reference case using the Limongi & Chieffi (2018) explodability criterion.

In the case where the acceleration efficiency is multiplied by a factor of 10, the Sco–Cen association still dominates over the other associations, but with a reduced probability of about 60 per cent. Indeed, with a higher acceleration efficiency more material is accelerated resulting in an increased predicted $^{60}$Fe CR density. It follows that OB associations other than Sco–Cen are then also able to account for the observed $^{60}$Fe CR density. There are then a larger number of realizations where a smaller contribution from Sco–Cen due to the initial draw of the massive star population is compensated by a larger contribution from associations such as Orion OB1 or Perseus OB2. This consequently reduces the occurrence of Sco–Cen as the most contributing OB association.

The green open stars in Fig. 9 represent the ratio between the average value of the calculated $^{56}$Fe CR density for the realizations, which account for the $^{60}$Fe observations, and the observed $^{56}$Fe CR density value. This can be taken as an indicator of the $^{56}$Fe CR density fraction that may come from a local source. This fraction is typically 20 per cent except for the cases where the CR diffusion coefficient and acceleration efficiency are varied by a factor of 10; for such cases the fraction is about 50 per cent. Since the nucleosynthesis activity for these two cases is the same, the higher predicted $^{56}$Fe CR density is correlated to an increase of the $^{60}$Fe CR density. The reason for such an increase when the CR acceleration efficiency is higher has previously been discussed. When the CR diffusion coefficient is

decreased, the propagation time needed to reach the maximum CR density for a given OB association distance is higher (see equation 12). Thus, there are possibly more ccSN that can explode during this time lapse, leading then to a higher predicted $^{56,60}$Fe CR density.

Since in most cases Sco–Cen is the OB association contributing the most to the LISM $^{60}$Fe CR density, it is interesting to explore the impact of the input parameters on the properties of the SNe accelerating the material present in the gas of the superbubble. Fig. 9 (bottom) shows the age (black) and mass (blue) of the Sco–Cen ccSN accelerating more than 50 per cent of the observed $^{60}$Fe CR density. The triangle markers represent the median (50th percentile) of the distribution while the lower (upper) bound of the error bar corresponds to the 16th (84th) percentile, defining a 68 per cent coverage probability. We see that most of the Sco–Cen realizations host a young SN with an age smaller than 300 kyr, except for two cases for which the explosion time distribution of the SN extends to higher values. This arises from different effects depending on the case. For smaller values of the CR diffusion coefficient, as discussed previously, one expects the age of the SNe to be greater and to span a larger range since the propagation time for CRs is longer. In the case where the CR acceleration efficiency is increased, older SNe, which would have contributed negligibly otherwise, may now contribute significantly to the observed $^{60}$Fe CR density budget. Interestingly, the mass of the Sco–Cen ccSN accelerating more than 50 per cent of the observed $^{60}$Fe CR density is nearly independent of the different test cases that we explored: the median value is $\approx 15–16 M_{\odot}$ and the $1\sigma$ range is $\approx 13–20 M_{\odot}$.

## 6.2 Geminga

The Geminga pulsar is currently located in the constellation Gemini at a distance of about $157_{-34}^{+59}$ pc with a proper motion of $170 \pm 4$ mas yr$^{-1}$ (Caraveo et al. 1996). Its spin-down age deduced from the pulsar period and period derivative (Bignami & Caraveo 1992) is 342 kyr, which can be considered as representative of its true age (Pellizza et al. 2005). Even if these properties are well established, the place of birth of Geminga is not clearly identified yet. By tracing back the space motion of Geminga, Pellizza et al. (2005) find that Geminga was born at 90–240 pc from the Sun, most probably inside the Orion OB1a association or the Cas-Tau OB association. Moreover, these authors conclude that the Geminga progenitor mass should not be greater than $15 M_{\odot}$.

One of the main conclusions of this study is that the Sco–Cen OB association plays a specific role in explaining the observed LISM $^{60}$Fe CR density. However, our results also show (i) that the Cas-Tau OB association is able to reproduce the observations, even though this is much more unlikely, and (ii) that the occurrence of Orion OB1 as the most contributing association can be significant when the Sukhbold et al. (2016) explodability criterion is used (see Fig. 10).

Here, we investigate whether the Geminga progenitor could be the SN that accelerated the $^{60}$Fe nuclei observed by ACE/CRIS, and whether this SN could be associated with the Orion OB1 or Cas-Tau associations.

In the following, we consider Orion OB1 and Cas-Tau as two independent OB associations. We first compute the nucleosynthetic activity as a function of time of each OB association, and we estimate, for each realization of our model ($N = 4000$), the amount of $^{60}$Fe present in the associated superbubble gas when the Geminga progenitor exploded 342 kyr ago. In a second step, and for each realization, the distance of Geminga to the Solar system *d* is uniformly sampled up to the distance of the considered OB associations. CRs are then accelerated and propagated across the distance *d* during a time
corresponding to the age of the Geminga pulsar, and their density is calculated using equation (10).

In the case of the Orion OB1 association, we find that in ≈6.5 per cent of the realizations (Sukhbold et al. 2016 explodability criterion) a ccSN exploding 342 kyr ago at a distance between 90 and 240 pc from the Sun is able to accelerate the 15 $^{60}$Fe nuclei detected by ACE/CRIS. On the contrary, we find that for the Cas-Tau OB association on its own no SN is able to reproduce the ACE/CRIS observations (Limongi & Chieffi 2018 explodability criterion). From these results the Geminga progenitor could be the SN that accelerated the $^{60}$Fe observed in the LISM if and only if it is associated with the Orion OB1 association. For these cases the $^{56}$Fe CR density is at least 10 times lower than the density measured in the LISM.

## 6.3 Locality of $^{56}$Fe CRs

An intriguing result of our work is that a substantial fraction of the $^{56}$Fe in the CR composition is found to be of local origin: of the order of ~20 per cent in most cases and up to ~50 per cent for special input parameters (see Fig. 9).

This fact can be understood by computing the maximum distance that a CR nucleus of $^{56}$Fe can diffuse in the ISM before undergoing catastrophic losses (spallation). The characteristic spallation time for $^{56}$Fe isotopes of energy equal to 550 MeV nucleon$^{-1}$ is $\tau_{spal} \sim 16 (n_H/0.1 \text{ cm}^{-3})^{-1}$ Myr, which translates into a maximum diffusion distance of

$$d_{\max} = (6 D \tau_{\text{spal}})^{1/2} \\ \sim 3.6 \left( \frac{D}{4 \times 10^{28} \text{ cm}^2 \text{ s}^{-1}} \right)^{1/2} \left( \frac{n_H}{0.1 \text{ cm}^{-3}} \right)^{-1/2} \text{kpc}, \quad (13)$$

where the diffusion coefficient has been normalized to the appropriate value for $^{56}$Fe CRs, and the gas density $n_H$ corresponds to the effective density experienced by CRs during their trip to the Solar system, while propagating both through the halo and disc (for sources located at a distance larger than the thickness of the gaseous disc, CRs spend a sizeable fraction of the propagation time in the halo). This means that the $^{56}$Fe CRs that we observe on the Earth have been produced at sources located at a distance smaller than $d_{\max}$.

More quantitatively, let us consider a situation where CRs are produced at a constant (in both space and time) rate at any location on the Galactic disc. Let $q_{CR}$ be the rate at which CRs of a given species and of a given energy are produced within an infinitesimal surface of the disc $d\sigma$. Then, an observer on the Earth would measure a density of CRs coming from a region $d\sigma$ located at a distance $R$ equal to $dn_{CR} = q_{CR}d\sigma/4\pi D R$ if $R < d_{\max}$ and $dn_{CR} \sim 0$ otherwise. Integrating over the entire surface of the disc one can see that the local density of CRs produced within a distance $R$ scales as $n_{CR}(\langle R \rangle = q_{CR}R/2D)$. Therefore, the fraction $\eta$ of observed CRs coming from within a distance $R$ is simply

$$\eta(\langle R \rangle) \sim \left(\frac{R}{d_{\max}}\right) \\ \sim 0.28 \left(\frac{R}{\text{kpc}}\right) \left(\frac{D}{4 \times 10^{28} \text{ cm}^2 \text{ s}^{-1}}\right)^{-1/2} \left(\frac{n_H}{0.1 \text{ cm}^{-3}}\right)^{1/2}. \quad (14)$$

This is true provided that $d_{\max}$ is smaller than the size of the CR halo $H$, otherwise $\eta \sim R/H$. Equation (14) shows that a fraction of about 30 per cent of the observed $^{56}$Fe is indeed expected to be produced in the LISM, where the star clusters listed in Table 1 are located. In fact, CRs are not injected at a constant rate at any location within the

disc, but are rather associated with SN explosions. The discreteness and stochasticity of stellar explosions plays a crucial role, and results in inhomogeneities in the spatial and temporal distribution of low-energy CRs (Phan et al. 2021). Therefore, the result obtained by means of equation (14) should be considered only as an indicative estimate.

## 7 SUMMARY AND CONCLUSIONS

Live $^{60}$Fe CRs have been detected in near-Earth space by the ACE/CRIS instrument over 17 yr of operation (Binns et al. 2016). The $^{60}$Fe radioactive lifetime of 3.8 Myr is sufficiently long such that an origin from a nearby nucleosynthesis site is plausible, and short enough so that the nucleosynthesis sites far out in the Galaxy are plausibly beyond reach for $^{60}$Fe surviving such a journey. In this paper, we thus investigated the possible local sources that may have accelerated the observed $^{60}$Fe nuclei.

We developed a bottom-up model computing the CR flux at the Solar system where the nucleosynthetic output from a massive star group is coupled to a CR transport model. The population synthesis part of our model relies on the yields from stars and SNe, which are properly weighted by an IMF using a Monte Carlo approach, addressing statistical fluctuations of stellar and star group parameters. The time profile of any nuclide abundance has thus been obtained in the gas of the superbubble that is excavated by the massive star cluster activity. We find that among the different ingredients of the population synthesis model the explodability criterion, which determines whether a massive star ends its life as an SN or avoids explosion, has the largest impact on the nuclide abundance in the superbubble.

Once the superbubble content in $^{60}$Fe is evaluated, we determine the fraction ending up in locally accelerated CRs, and propagate these from their source through the ISM toward the Solar system. We consider a simple acceleration and propagation model where the advection and ionization energy losses can be neglected, and where accelerated ions, when escaping from their source, diffuse isotropically in the ISM and suffer both catastrophic and radioactive losses. Both the CR acceleration efficiency and diffusion coefficient are very uncertain, in part because of the structure of the magnetic field and the superbubble environment (diffusion coefficient), and the efficiency of dust production and its destruction by thermal sputtering (acceleration efficiency).

When applying our CR population synthesis and transport model to all the OB associations within 1 kpc of our Solar system (Wright 2020) we find that the 15 nuclei of $^{60}$Fe detected by the ACE/CRIS instrument most probably originate from the Sco-Cen OB association. Moreover, we find that a young SN (age $\le 500$ kyr) with a progenitor mass of $\approx 13-20$ M$_{\odot}$ might be the source of acceleration of the observed $^{60}$Fe nuclei. These results are largely independent of the assumed input parameters of our model except for the explodability criterion. When the Sukhbold et al. (2016) criterion is used, the Orion OB1 association may also contribute significantly to the observed $^{60}$Fe CR density in the LISM.

The Orion OB1 association and the Cas-Tau OB association are both possible birthplaces of the Geminga pulsar (Pellizza et al. 2005). We investigate the possibility that the observed $^{60}$Fe nuclei were accelerated by the SN explosion that gave birth to the Geminga pulsar, and we show that a ccSN exploding 342 kyr ago (age of Geminga) at a distance between 90 and 240 pc from the Sun (presumed distance of Geminga at its birth) can account for the observed $^{60}$Fe CR density in the LISM if and only if the progenitor of Geminga is located in the
corresponding to the age of the Geminga pulsar, and their density is calculated using equation (10).

In the case of the Orion OB1 association, we find that in ≈6.5 per cent of the realizations (Sukhbold et al. 2016 explodability criterion) a ccSN exploding 342 kyr ago at a distance between 90 and 240 pc from the Sun is able to accelerate the 15 $^{60}$Fe nuclei detected by ACE/CRIS. On the contrary, we find that for the Cas-Tau OB association on its own no SN is able to reproduce the ACE/CRIS observations (Limongi & Chieffi 2018 explodability criterion). From these results the Geminga progenitor could be the SN that accelerated the $^{60}$Fe observed in the LISM if and only if it is associated with the Orion OB1 association. For these cases the $^{56}$Fe CR density is at least 10 times lower than the density measured in the LISM.

## 6.3 Locality of $^{56}$Fe CRs

An intriguing result of our work is that a substantial fraction of the $^{56}$Fe in the CR composition is found to be of local origin: of the order of ~20 per cent in most cases and up to ~50 per cent for special input parameters (see Fig. 9).

This fact can be understood by computing the maximum distance that a CR nucleus of $^{56}$Fe can diffuse in the ISM before undergoing catastrophic losses (spallation). The characteristic spallation time for $^{56}$Fe isotopes of energy equal to 550 MeV nucleon$^{-1}$ is $\tau_{spal} \sim 16 (n_H/0.1 \text{ cm}^{-3})^{-1}$ Myr, which translates into a maximum diffusion distance of

$$d_{\max} = (6 D \tau_{\text{spal}})^{1/2} \\ \sim 3.6 \left( \frac{D}{4 \times 10^{28} \text{ cm}^2 \text{ s}^{-1}} \right)^{1/2} \left( \frac{n_H}{0.1 \text{ cm}^{-3}} \right)^{-1/2} \text{kpc}, \quad (13)$$

where the diffusion coefficient has been normalized to the appropriate value for $^{56}$Fe CRs, and the gas density $n_H$ corresponds to the effective density experienced by CRs during their trip to the Solar system, while propagating both through the halo and disc (for sources located at a distance larger than the thickness of the gaseous disc, CRs spend a sizeable fraction of the propagation time in the halo). This means that the $^{56}$Fe CRs that we observe on the Earth have been produced at sources located at a distance smaller than $d_{\max}$.

More quantitatively, let us consider a situation where CRs are produced at a constant (in both space and time) rate at any location on the Galactic disc. Let $q_{CR}$ be the rate at which CRs of a given species and of a given energy are produced within an infinitesimal surface of the disc $d\sigma$. Then, an observer on the Earth would measure a density of CRs coming from a region $d\sigma$ located at a distance $R$ equal to $dn_{CR} = q_{CR}d\sigma/4\pi D R$ if $R < d_{\max}$ and $dn_{CR} \sim 0$ otherwise. Integrating over the entire surface of the disc one can see that the local density of CRs produced within a distance $R$ scales as $n_{CR}(\langle R \rangle = q_{CR}R/2D)$. Therefore, the fraction $\eta$ of observed CRs coming from within a distance $R$ is simply

$$\eta(\langle R \rangle) \sim \left( \frac{R}{d_{\max}} \right) \\ \sim 0.28 \left( \frac{R}{\text{kpc}} \right) \left( \frac{D}{4 \times 10^{28} \text{ cm}^2 \text{ s}^{-1}} \right)^{-1/2} \left( \frac{n_H}{0.1 \text{ cm}^{-3}} \right)^{1/2} . \quad (14)$$

This is true provided that $d_{\max}$ is smaller than the size of the CR halo $H$, otherwise $\eta \sim R/H$. Equation (14) shows that a fraction of about 30 per cent of the observed $^{56}$Fe is indeed expected to be produced in the LISM, where the star clusters listed in Table 1 are located. In fact, CRs are not injected at a constant rate at any location within the

disc, but are rather associated with SN explosions. The discreteness and stochasticity of stellar explosions plays a crucial role, and results in inhomogeneities in the spatial and temporal distribution of low-energy CRs (Phan et al. 2021). Therefore, the result obtained by means of equation (14) should be considered only as an indicative estimate.

## 7 SUMMARY AND CONCLUSIONS

Live $^{60}$Fe CRs have been detected in near-Earth space by the ACE/CRIS instrument over 17 yr of operation (Binns et al. 2016). The $^{60}$Fe radioactive lifetime of 3.8 Myr is sufficiently long such that an origin from a nearby nucleosynthesis site is plausible, and short enough so that the nucleosynthesis sites far out in the Galaxy are plausibly beyond reach for $^{60}$Fe surviving such a journey. In this paper, we thus investigated the possible local sources that may have accelerated the observed $^{60}$Fe nuclei.

We developed a bottom-up model computing the CR flux at the Solar system where the nucleosynthetic output from a massive star group is coupled to a CR transport model. The population synthesis part of our model relies on the yields from stars and SNe, which are properly weighted by an IMF using a Monte Carlo approach, addressing statistical fluctuations of stellar and star group parameters. The time profile of any nuclide abundance has thus been obtained in the gas of the superbubble that is excavated by the massive star cluster activity. We find that among the different ingredients of the population synthesis model the explodability criterion, which determines whether a massive star ends its life as an SN or avoids explosion, has the largest impact on the nuclide abundance in the superbubble.

Once the superbubble content in $^{60}$Fe is evaluated, we determine the fraction ending up in locally accelerated CRs, and propagate these from their source through the ISM toward the Solar system. We consider a simple acceleration and propagation model where the advection and ionization energy losses can be neglected, and where accelerated ions, when escaping from their source, diffuse isotropically in the ISM and suffer both catastrophic and radioactive losses. Both the CR acceleration efficiency and diffusion coefficient are very uncertain, in part because of the structure of the magnetic field and the superbubble environment (diffusion coefficient), and the efficiency of dust production and its destruction by thermal sputtering (acceleration efficiency).

When applying our CR population synthesis and transport model to all the OB associations within 1 kpc of our Solar system (Wright 2020) we find that the 15 nuclei of $^{60}$Fe detected by the ACE/CRIS instrument most probably originate from the Sco-Cen OB association. Moreover, we find that a young SN (age $\le 500$ kyr) with a progenitor mass of $\approx 13-20$ M$_{\odot}$ might be the source of acceleration of the observed $^{60}$Fe nuclei. These results are largely independent of the assumed input parameters of our model except for the explodability criterion. When the Sukhbold et al. (2016) criterion is used, the Orion OB1 association may also contribute significantly to the observed $^{60}$Fe CR density in the LISM.

The Orion OB1 association and the Cas-Tau OB association are both possible birthplaces of the Geminga pulsar (Pellizza et al. 2005). We investigate the possibility that the observed $^{60}$Fe nuclei were accelerated by the SN explosion that gave birth to the Geminga pulsar, and we show that a ccSN exploding 342 kyr ago (age of Geminga) at a distance between 90 and 240 pc from the Sun (presumed distance of Geminga at its birth) can account for the observed $^{60}$Fe CR density in the LISM if and only if the progenitor of Geminga is located in the
Orion OB1 association. The associated probability for such a case is about 6–7 per cent.

The origin of the live $^{60}$Fe nuclei detected by the ACE/CRIS instrum-
strument could be traced back to the closest nearby OB associations.
With the same formalism we computed the CR density of radioactive
$^{26}$Al and stable $^{56}$Fe nuclei in the LISM. We find that the $^{26}$Al
density calculated from local OB associations is more than an order
of magnitude lower than that deduced from ACR/CRIS observations,
which confirms that $^{26}$Al in CRs is mainly a secondary species
produced by spallation of heavier nuclei (mainly $^{28}$Si). However,
we also find that about 20 per cent of the observed $^{56}$Fe density
can be accounted for by local OB associations located at less than
~1 kpc from the Solar system. These results are independent of the
population synthesis parameters (IMF, yields, and explodability), but
do show a sensitivity to the CR acceleration efficiency and diffusion
coefficient. Varying by a factor of 10 down and up the CR acceleration
efficiency and the diffusion coefficient, respectively, the $^{56}$Fe density
calculated from local OB associations can represent up to 50 per cent
of the observed value. Overall, the calculated contribution of local
sources to the $^{56}$Fe CR population appears to be consistent with a
simple estimate assuming homogeneous CR production at a constant
rate across the Galactic disc.

ACKNOWLEDGEMENTS

SG acknowledges support from Agence Nationale de la Recherche
(grant ANR-21-CE31-0028).

DATA AVAILABILITY

The stellar yields used in this work are publicly available from the
works of Limongi & Chieffi (2018) and Ebinger et al. (2019). Data
arising from this work are available on reasonable request to the
corresponding author.

REFERENCES

Ackermann M. et al., 2011, *Science*, **334**, 1103

Aguilar M. et al., 2021, *Phys. Rev. Lett.*, 126, 041104

Aharonian F., Yang R., de Oña Wilhelmi E., 2019, *Nat. Astron.*, 3, 561

Armstrong J. J., Wright N. J., Jeffries R. D., 2018, MNRAS, 480, L121

Astiasarain X., Tibaldo L., Martin P., Knödlseder J., Remy Q., 2023, A&A,
671, A47

Belikov A. N., Kharchenko N. V., Piskunov A. E., Schilbach E., Scholz R. D., 2002, A&A, 387, 117

Berezhko E. G., Ellison D. C., 1999, ApJ, 526, 385

Berlanas S. R. et al., 2020, A&A, 642, A168

Bignami G. F., Caraveo P. A., 1992, *Nature*, **357**, 287

Binns W. R. et al., 2007, in von Steiger R., Gloeckler G., Mason G. M., eds, The Composition of Matter: Symposium Honouring Johannes Geiss on the Occasion of His 80th Birthday. Springer, New York, p. 439

Binns W. R. et al., 2016, *Science*, **352**, 677

Blaauw A., 1964, ARA&A, **2**, 213

Blandford R. D., Ostriker J. P., 1978, ApJ, 221, L29

Blasi P., 2013, A&AR, 21, 70

Boschini M. J., Della Torre S., Gervasi M., La Vacca G., Rancoita P. G., 2019, *Adv. Space Res.*, **64**, 2459

Boschini M. J. et al., 2021, ApJ, 913, 5

Boschini M. J. et al., 2022, ApJ, 933, 147

Breitschwerdt D., Feige J., Schulreich M. M., Avillez M. A. D., Dettbarn C., Fuchs B., 2016, *Nature*, **532**, 73

Brown A. G. A., de Geus E. J., de Zeeuw P. T., 1994, A&A, 289, 101

Bykov A. M., Kalyashova M. E., 2022, *Adv. Space Res.*, **70**, 2685

Caraveo P. A., Bignami G. F., Mignani R., Taff L. G., 1996, ApJ, 461, L91

Cerviño M., 2013, *New Astron. Rev.*, 57, 123

Cerviño M., Luridiana V., 2006, A&A, 451, 475

Chieffi A., Limongi M., 2013, ApJ, 764, 21

Côté B., Yagüe A., Világos B., Lugaro M., 2019, ApJ, 887, 213

Curtis S., Ebinger K., Fröhlich C., Hempel M., Perego A., Liebendörfer M., Thielemann F.-K., 2019, ApJ, 870, 2

de Zeeuw P. T., Hoogerwerf R., de Bruijne J. H. J., Brown A. G. A., Blaauw A., 1999, AJ, 117, 354

Diehl R. et al., 2021, *Publ. Astron. Soc. Aust.*, 38, e062

Drury L. O. C., 2012, *Astropart. Phys.*, 39, 52

Ebinger K., Curtis S., Fröhlich C., Hempel M., Perego A., Liebendörfer M., Thielemann F.-K., 2019, ApJ, 870, 1

Ellison D. C., Drury L. O., Meyer J.-P., 1997, ApJ, 487, 197

Evoli C., Aloisio R., Blasi P., 2019, *Phys. Rev. D*, **99**, 103023

Farmer A. J., Goldreich P., 2004, ApJ, 604, 671

Foglizzo T. et al., 2015, *Publ. Astron. Soc. Aust.*, 32, e009

Fuchs B., Breitschwerdt D., de Avillez M. A., Dettbarn C., Flynn C., 2006, MNRAS, 373, 993

Gabici S., 2022, A&AR, 30, 4

Gabici S., Evoli C., Gaggero D., Lipari P., Mertsch P., Orlando E., Strong A., Vittino A., 2019, *Int. J. Mod. Phys. D*, **28**, 1930022

Ginzburg V. L., Syrovatskii S. I., 1964, The Origin of Cosmic Rays. Pergamon Press, Oxford

Głěbocki R., Graciński P., 2005, in Favata F., Hussain G. A. J., Bat-trick B., eds, Proceedings of the 13th Cambridge Workshop on Cool Stars, Stellar Systems and the Sun (ESA SP-560). ESA, Noordwijk, p. 571

Gleeson L. J., Axford W. I., 1968, ApJ, 154, 1011

Gounelle M., Meibom A., Hennebelle P., Inutsuka S.-i., 2009, ApJ, 694, L1

Habets G. M. H. J., Heintze J. R. W., 1981, A&AS, 46, 193

Heger A., Fryer C. L., Woosley S. E., Langer N., Hartmann D. H., 2003, ApJ,
591, 288

HESS Collaboration, 2016, *Nature*, **531**, 476

Hillenbrand L. A., 1997, AJ, 113, 1733

Hoppe P., Stancliffe R. J., Pignatari M., Amari S., 2019, ApJ, 887, 8

Indriolo N., McCall B. J., 2012, ApJ, 745, 91

Janka H.-T., 2012, *Annu. Rev. Nucl. Part. Sci.*, **62**, 407

Kachelrieß M., Neronov A., Semikoz D. V., 2018, Phys. Rev. D, 97, 063011

Kempski P., Quataert E., 2022, MNRAS, 514, 657

Knie K., Korschinek G., Faestermann T., Dorfi E. A., Rugel G., Wallner A., 2004, Phys. Rev. Lett., 93, 171103

Kroupa P., 2001, MNRAS, 322, 231

Kroupa P., 2002, in Grebel E. K., Brandner W., eds, ASP Conf. Ser. Vol. 285, Modes of Star Formation and the Origin of Field Populations. Astron. Soc. Pac., San Francisco, p. 86

Kroupa P., 2020, in Bragglia A., Davies M., Sills A., Vesperini E., eds, Proc. IAU Symp. Vol. 351, Star Clusters: From the Milky Way to the Early Universe. Cambridge Univ. Press, Cambridge, p. 117

Kulsrud R., Pearce W. P., 1969, ApJ, 156, 445

Lada C. J., 2005, Prog. Theor. Phys. Suppl., 158, 1

Lazarian A., Xu S., 2021, ApJ, 923, 53

Lazarian A., Xu S., Hu Y., 2023, *Frontiers Astron. Space Sci.*, 10, 1154760

Lee S.-H., Ellison D. C., Nagataki S., 2012, ApJ, 750, 156

Limongi M., Chieffi A., 2018, ApJS, 237, 13

Mańcz Appellániz J., Walborn N. R., Morrell N. I., Niemela V. S., Nelan E. P., 2007, ApJ, 660, 1480

Mannheim K., Schlickeiser R., 1994, A&A, 286, 983

Matsuura M. et al., 2019, MNRAS, 482, 1715

Mel'nik A. M., Dambis A. K., 2017, MNRAS, 472, 3887

Meyer J.-P., Drury L. O., Ellison D. C., 1997, ApJ, 487, 182

Meynet G., Maeder A., 2000, A&A, 361, 101

Miller G. E., Scalo J. M., 1979, ApJS, 41, 513

Morlino G., Amato E., 2020, Phys. Rev. D, 101, 083017

Murphy R. P. et al., 2016, ApJ, 831, 148

Nittler L. R., Amari S., Zinner E., Woosley S. E., Lewis R. S., 1996, ApJ,
462, L31

O'Connor E., Ott C. D., 2011, ApJ, 730, 70

Palla F., Stahler S. W., 1999, ApJ, 525, 772
Orion OB1 association. The associated probability for such a case is about 6–7 per cent.

The origin of the live $^{60}$Fe nuclei detected by the ACE/CRIS instrum-
strument could be traced back to the closest nearby OB associations.
With the same formalism we computed the CR density of radioactive
$^{26}$Al and stable $^{56}$Fe nuclei in the LISM. We find that the $^{26}$Al
density calculated from local OB associations is more than an order
of magnitude lower than that deduced from ACR/CRIS observations,
which confirms that $^{26}$Al in CRs is mainly a secondary species
produced by spallation of heavier nuclei (mainly $^{28}$Si). However,
we also find that about 20 per cent of the observed $^{56}$Fe density
can be accounted for by local OB associations located at less than
~1 kpc from the Solar system. These results are independent of the
population synthesis parameters (IMF, yields, and explodability), but
do show a sensitivity to the CR acceleration efficiency and diffusion
coefficient. Varying by a factor of 10 down and up the CR acceleration
efficiency and the diffusion coefficient, respectively, the $^{56}$Fe density
calculated from local OB associations can represent up to 50 per cent
of the observed value. Overall, the calculated contribution of local
sources to the $^{56}$Fe CR population appears to be consistent with a
simple estimate assuming homogeneous CR production at a constant
rate across the Galactic disc.

ACKNOWLEDGEMENTS

SG acknowledges support from Agence Nationale de la Recherche
(grant ANR-21-CE31-0028).

DATA AVAILABILITY

The stellar yields used in this work are publicly available from the
works of Limongi & Chieffi (2018) and Ebinger et al. (2019). Data
arising from this work are available on reasonable request to the
corresponding author.

REFERENCES

Ackermann M. et al., 2011, *Science*, **334**, 1103

Aguilar M. et al., 2021, *Phys. Rev. Lett.*, 126, 041104

Aharonian F., Yang R., de Oña Wilhelmi E., 2019, *Nat. Astron.*, 3, 561

Armstrong J. J., Wright N. J., Jeffries R. D., 2018, MNRAS, 480, L121

Astiasarain X., Tibaldo L., Martin P., Knödlseder J., Remy Q., 2023, A&A,
671, A47

Belikov A. N., Kharchenko N. V., Piskunov A. E., Schilbach E., Scholz R. D., 2002, A&A, 387, 117

Berezhko E. G., Ellison D. C., 1999, ApJ, 526, 385

Berlanas S. R. et al., 2020, A&A, 642, A168

Bignami G. F., Caraveo P. A., 1992, *Nature*, **357**, 287

Binns W. R. et al., 2007, in von Steiger R., Gloeckler G., Mason G. M., eds, The Composition of Matter: Symposium Honouring Johannes Geiss on the Occasion of His 80th Birthday. Springer, New York, p. 439

Binns W. R. et al., 2016, *Science*, **352**, 677

Blaauw A., 1964, ARA&A, **2**, 213

Blandford R. D., Ostriker J. P., 1978, ApJ, 221, L29

Blasi P., 2013, A&AR, 21, 70

Boschini M. J., Della Torre S., Gervasi M., La Vacca G., Rancoita P. G., 2019, *Adv. Space Res.*, **64**, 2459

Boschini M. J. et al., 2021, ApJ, 913, 5

Boschini M. J. et al., 2022, ApJ, 933, 147

Breitschwerdt D., Feige J., Schulreich M. M., Avillez M. A. D., Dettbarn C., Fuchs B., 2016, *Nature*, **532**, 73

Brown A. G. A., de Geus E. J., de Zeeuw P. T., 1994, A&A, 289, 101

Bykov A. M., Kalyashova M. E., 2022, *Adv. Space Res.*, **70**, 2685

Caraveo P. A., Bignami G. F., Mignani R., Taff L. G., 1996, ApJ, 461, L91

Cerviño M., 2013, *New Astron. Rev.*, 57, 123

Cerviño M., Luridiana V., 2006, A&A, 451, 475

Chieffi A., Limongi M., 2013, ApJ, 764, 21

Côté B., Yagüe A., Világos B., Lugaro M., 2019, ApJ, 887, 213

Curtis S., Ebinger K., Fröhlich C., Hempel M., Perego A., Liebendörfer M., Thielemann F.-K., 2019, ApJ, 870, 2

de Zeeuw P. T., Hoogerwerf R., de Bruijne J. H. J., Brown A. G. A., Blaauw
A., 1999, AJ, 117, 354

Diehl R. et al., 2021, *Publ. Astron. Soc. Aust.*, 38, e062

Drury L. O. C., 2012, *Astropart. Phys.*, 39, 52

Ebinger K., Curtis S., Fröhlich C., Hempel M., Perego A., Liebendörfer M., Thielemann F.-K., 2019, ApJ, 870, 1

Ellison D. C., Drury L. O., Meyer J.-P., 1997, ApJ, 487, 197

Evoli C., Aloisio R., Blasi P., 2019, *Phys. Rev. D*, **99**, 103023

Farmer A. J., Goldreich P., 2004, ApJ, 604, 671

Foglizzo T. et al., 2015, *Publ. Astron. Soc. Aust.*, 32, e009

Fuchs B., Breitschwerdt D., de Avillez M. A., Dettbarn C., Flynn C., 2006, MNRAS, 373, 993

Gabici S., 2022, A&AR, 30, 4

Gabici S., Evoli C., Gaggero D., Lipari P., Mertsch P., Orlando E., Strong A., Vittino A., 2019, *Int. J. Mod. Phys. D*, **28**, 1930022

Ginzburg V. L., Syrovatskii S. I., 1964, The Origin of Cosmic Rays. Pergamon Press, Oxford

Głěbocki R., Graciński P., 2005, in Favata F., Hussain G. A. J., Bat-trick B., eds, Proceedings of the 13th Cambridge Workshop on Cool Stars, Stellar Systems and the Sun (ESA SP-560). ESA, Noordwijk, p. 571

Gleeson L. J., Axford W. I., 1968, ApJ, 154, 1011

Gounelle M., Meibom A., Hennebelle P., Inutsuka S.-i., 2009, ApJ, 694, L1

Habets G. M. H. J., Heintze J. R. W., 1981, A&AS, 46, 193

Heger A., Fryer C. L., Woosley S. E., Langer N., Hartmann D. H., 2003, ApJ,
591, 288

HESS Collaboration, 2016, *Nature*, 531, 476

Hillenbrand L. A., 1997, AJ, 113, 1733

Hoppe P., Stancliffe R. J., Pignatari M., Amari S., 2019, ApJ, 887, 8

Indriolo N., McCall B. J., 2012, ApJ, 745, 91

Janka H.-T., 2012, *Annu. Rev. Nucl. Part. Sci.*, **62**, 407

Kachelrieß M., Neronov A., Semikoz D. V., 2018, Phys. Rev. D, 97, 063011

Kempski P., Quataert E., 2022, MNRAS, 514, 657

Knie K., Korschinek G., Faestermann T., Dorfi E. A., Rugel G., Wallner A., 2004, Phys. Rev. Lett., 93, 171103

Kroupa P., 2001, MNRAS, 322, 231

Kroupa P., 2002, in Grebel E. K., Brandner W., eds, ASP Conf. Ser. Vol. 285, Modes of Star Formation and the Origin of Field Populations. Astron. Soc. Pac., San Francisco, p. 86

Kroupa P., 2020, in Bragglia A., Davies M., Sills A., Vesperini E., eds, Proc. IAU Symp. Vol. 351, Star Clusters: From the Milky Way to the Early Universe. Cambridge Univ. Press, Cambridge, p. 117

Kulsrud R., Pearce W. P., 1969, ApJ, 156, 445

Lada C. J., 2005, Prog. Theor. Phys. Suppl., 158, 1

Lazarian A., Xu S., 2021, ApJ, 923, 53

Lazarian A., Xu S., Hu Y., 2023, *Frontiers Astron. Space Sci.*, 10, 1154760

Lee S.-H., Ellison D. C., Nagataki S., 2012, ApJ, 750, 156

Limongi M., Chieffi A., 2018, ApJS, 237, 13

Mańcz Appellániz J., Walborn N. R., Morrell N. I., Niemela V. S., Nelan E. P., 2007, ApJ, 660, 1480

Mannheim K., Schlickeiser R., 1994, A&A, 286, 983

Matsuura M. et al., 2019, MNRAS, 482, 1715

Mel'nik A. M., Dambis A. K., 2017, MNRAS, 472, 3887

Meyer J.-P., Drury L. O., Ellison D. C., 1997, ApJ, 487, 182

Meynet G., Maeder A., 2000, A&A, 361, 101

Miller G. E., Scalo J. M., 1979, ApJS, 41, 513

Morlino G., Amato E., 2020, Phys. Rev. D, 101, 083017

Murphy R. P. et al., 2016, ApJ, 831, 148

Nittler L. R., Amari S., Zinner E., Woosley S. E., Lewis R. S., 1996, ApJ,
462, L31

O'Connor E., Ott C. D., 2011, ApJ, 730, 70

Palla F., Stahler S. W., 1999, ApJ, 525, 772
Pantaleoni González M., Maíz Apellániz J., Barbá R. H., Reed B. C., 2021, MNRAS, 504, 2968

Parizot E., Marcowith A., van der Swaluw E., Bykov A. M., Tatischeff V., 2004, A&A, 424, 747

Pellizza L. J., Mignani R. P., Grenier I. A., Mirabel I. F., 2005, A&A, 435, 625

Peron G., Aharonian F., Casanova S., Yang R., Zanin R., 2021, ApJ, 907, L11

Phan V. H. M., Schulze F., Mertsch P., Recchia S., Gabici S., 2021, Phys. Rev. Lett., 127, 141101

Phan V. H. M., Recchia S., Mertsch P., Gabici S., 2023, Phys. Rev. D, 107, 123006

Prantzos N., Abia C., Limongi M., Chieffi A., Cristallo S., 2018, MNRAS, 476, 3432

Preibisch T., Brown A. G. A., Bridges T., Guenther E., Zinnecker H., 2002, AJ, 124, 404

Salpeter E. E., 1955, ApJ, 121, 161

Schneider F. R. N. et al., 2018, Science, 359, 69

Strong A. W., Moskalenko I. V., 1998, ApJ, 509, 212

Strong A. W., Porter T. A., Digel S. W., Jóhannesson G., Martin P., Moskalenko I. V., Murphy E. J., Orlando E., 2010, ApJ, 722, L58

Sukhbold T., Ertl T., Woosley S. E., Brown J. M., Janka H. T., 2016, ApJ, 821, 38

Tatischeff V., Raymond J. C., Duprat J., Gabici S., Recchia S., 2021, MNRAS, 508, 1321

Tripathi R. K., Cucinotta F. A., Wilson J. W., 1996, Nucl. Instrum. Methods Phys. Res. B, 117, 347

Tripathi R. K., Cucinotta F. A., Wilson J. W., 1999, Nucl. Instrum. Methods Phys. Res. B, 155, 349

Vanbeveren D., 2009, New Astron. Rev., 53, 27

Vasileiadis A., Nordlund Å., Bizzarro M., 2013, ApJ, 769, L8

Vieu T., Gabici S., Tatischeff V., Ravikularaman S., 2022, MNRAS, 512, 1275

Voss R., Diehl R., Hartmann D. H., Cerviño M., Vink J. S., Meynet G., Limongi M., Chieffi A., 2009, A&A, 504, 531

Wallner A. et al., 2016, Nature, 532, 69

Wallner A. et al., 2021, Science, 372, 742

Weidner C., Vink J. S., 2010, A&A, 524, A98

Wiedenbeck M. E. et al., 1999, ApJ, 523, L61

Woosley S. E., Heger A., 2007, Phys. Rep., 442, 269

Wright N. J., 2020, New Astron. Rev., 90, 101549

Yan Z., Jerabkova T., Kroupa P., 2023, A&A, 670, A151

Yanasak N. E. et al., 2001, ApJ, 563, 768

Young E. D., 2014, Earth Planet. Sci. Lett., 392, 16

Zucker C. et al., 2022, Nature, 601, 334

Zucker C., Alves J., Goodman A., Meingast S., Galli P., 2023, in Inutsuka S.-i., Aikawa Y., Muto T., Tomida K., Tamura M., eds, ASP Conf. Ser. Vol. 534, Protostars and Planets VII. Astron. Soc. Pac., San Francisco, p. 43

This paper has been typeset from a TeX/Latex file prepared by the author.
Pantaleoni González M., Maíz Apellániz J., Barbá R. H., Reed B. C., 2021, MNRAS, 504, 2968

Parizot E., Marcowith A., van der Swaluw E., Bykov A. M., Tatischeff V., 2004, A&A, 424, 747

Pellizza L. J., Mignani R. P., Grenier I. A., Mirabel I. F., 2005, A&A, 435, 625

Peron G., Aharonian F., Casanova S., Yang R., Zanin R., 2021, ApJ, 907, L11

Phan V. H. M., Schulze F., Mertsch P., Recchia S., Gabici S., 2021, Phys. Rev. Lett., 127, 141101

Phan V. H. M., Recchia S., Mertsch P., Gabici S., 2023, Phys. Rev. D, 107, 123006

Prantzos N., Abia C., Limongi M., Chieffi A., Cristallo S., 2018, MNRAS, 476, 3432

Preibisch T., Brown A. G. A., Bridges T., Guenther E., Zinnecker H., 2002, AJ, 124, 404

Salpeter E. E., 1955, ApJ, 121, 161

Schneider F. R. N. et al., 2018, Science, 359, 69

Strong A. W., Moskalenko I. V., 1998, ApJ, 509, 212

Strong A. W., Porter T. A., Digel S. W., Jóhannesson G., Martin P., Moskalenko I. V., Murphy E. J., Orlando E., 2010, ApJ, 722, L58

Sukhbold T., Ertl T., Woosley S. E., Brown J. M., Janka H. T., 2016, ApJ, 821, 38

Tatischeff V., Raymond J. C., Duprat J., Gabici S., Recchia S., 2021, MNRAS, 508, 1321

Tripathi R. K., Cucinotta F. A., Wilson J. W., 1996, Nucl. Instrum. Methods Phys. Res. B, 117, 347

Tripathi R. K., Cucinotta F. A., Wilson J. W., 1999, Nucl. Instrum. Methods Phys. Res. B, 155, 349

Vanbeveren D., 2009, New Astron. Rev., 53, 27

Vasileiadis A., Nordlund Å., Bizzarro M., 2013, ApJ, 769, L8

Vieu T., Gabici S., Tatischeff V., Ravikularaman S., 2022, MNRAS, 512, 1275

Voss R., Diehl R., Hartmann D. H., Cerviño M., Vink J. S., Meynet G., Limongi M., Chieffi A., 2009, A&A, 504, 531

Wallner A. et al., 2016, Nature, 532, 69

Wallner A. et al., 2021, Science, 372, 742

Weidner C., Vink J. S., 2010, A&A, 524, A98

Wiedenbeck M. E. et al., 1999, ApJ, 523, L61

Woosley S. E., Heger A., 2007, Phys. Rep., 442, 269

Wright N. J., 2020, New Astron. Rev., 90, 101549

Yan Z., Jerabkova T., Kroupa P., 2023, A&A, 670, A151

Yanasak N. E. et al., 2001, ApJ, 563, 768

Young E. D., 2014, Earth Planet. Sci. Lett., 392, 16

Zucker C. et al., 2022, Nature, 601, 334

Zucker C., Alves J., Goodman A., Meingast S., Galli P., 2023, in Inutsuka S.-i., Aikawa Y., Muto T., Tomida K., Tamura M., eds, ASP Conf. Ser. Vol. 534, Protostars and Planets VII. Astron. Soc. Pac., San Francisco, p. 43

This paper has been typeset from a TeX/LATEX file prepared by the author.