© ESO, 2017

1. Introduction

Our solar system resides in a low-density region of the interstellar medium (ISM), partially filled with hot, soft X-ray emitting plasma, known as the Local Bubble (LB). While the existence of this superbubble (SB) is widely accepted today, not least because it was found to be responsible for about 60% of the 1 / 4-keV flux in the Galactic plane (Galeazzi et al. 2014), its exact origin is still a matter of ongoing debate. The LB extends roughly 200 pc in the Galactic plane, and 600 pc perpendicular to it, with an inclination of about 20° relative to the axis of Galactic rotation (cf. Lallement et al. 2003). If we assume that it was created by several nearby supernova (SN) explosions in the last ≃10 Myr (e.g. Smith & Cox 2001), we are faced with the problem that no young stellar cluster could be found inside its boundaries.

One obvious way out of this dilemma was to scan the solar neighbourhood for moving groups of stars that might have once crossed the boundary of the LB and, additionally, provided an adequate number of SNe in order to explain both the present LB extent and the detected soft X-ray emission. For that purpose, Berghöfer & Breitschwerdt (2002) calculated stellar trajectories of the nearby Pleiades subgroup B1 backwards in time. Applying an integer-binned initial mass function (IMF) they inferred the number of SNe in the past, determined their explosion sites from the trajectories, and their explosion times from their main-sequence life times (assuming that all stars in the cluster are born coevally). Similar conclusions were drawn from an independent study by Maíz-Apellániz (2001). Later, Fuchs et al. (2006) performed an ab initio search scrutinising all stars with sufficient kinematical data available within a heliocentric volume of 400 pc diameter. Using the Hipparcos catalogue and a compilation of radial velocities, they came up with a sample of B stars concentrated in both real and velocity space, which were then identified as members of the Sco-Cen OB association. Tracing back in time this sample, Fuchs et al. found that the stars had entered the present LB volume rather off-centre 10 to 15 Myr ago, and that since then, 14 to 20 SN explosions should have occurred, which delivered enough energy to explain the current size of the LB and its ionisation state.

A deeper understanding of the nonlinear evolution and structure of the ISM in general and the LB in particular can only be gained from solving the full-blown set of hydrodynamical (HD) and magnetohydrodynamical (MHD) equations through numerical simulations. While early studies of the global ISM could only consider very few physical processes on small, two-dimensional domains at low numerical resolution (Bania & Lyon 1980; Chiang & Prendergast 1985; Chiang & Bregman 1988; Rosen et al. 1993), the advance in computational capabilities not only allowed for the inclusion of additional physics (like magnetic fields) and larger-scale domains (Passot et al. 1996; Rosen & Bregman 1995), but eventually also the leap into the third dimension (Mac Low et al. 1998; Stone et al. 1998; Padoan & Nordlund 1999). More recent works can be roughly divided into three classes according to their box size (and the highest resolution associated with it): small-scale (Balsara et al. 2004; Piontek & Ostriker 2007), meso-scale (de Avillez & Breitschwerdt 2004, 2005; Slyz et al. 2005; Gressel et al. 2008; Joung et al. 2009; Girichidis et al. 2015), and large-scale (Hanasz et al. 2009; Dobbs et al. 2011). A physically realistic three-dimensional model that covers all dynamical ranges of an entire galaxy is still out of reach, though there are efforts to achieve this to some degree by coupling the adaptive mesh refinement (AMR) technique with zoom-ins to specific regions (see e.g. Renaud et al. 2013, and follow-up studies).

Taking their modelled ISM as a typical background medium, Breitschwerdt & de Avillez (2006) further tested the aforementioned LB formation scenario under more realistic conditions. They found that the bubble’s extension, as well as its ion column density ratios (see de Avillez & Breitschwerdt 2009, 2012), matched ultraviolet data, as obtained e.g. with FUSE (Oegerle et al. 2005; Savage & Lehner 2006), remarkably well.

Another interesting way to validate evolution models for the LB is to search for traces of the involved SNe on Earth, manifesting themselves as geological isotope anomalies. Ellis et al. (1996) identified several possible target isotopes for such a study, including ⁶⁰Fe (see also Korschinek et al. 1996), which is believed to be produced primarily in core-collapse and electron-capture SNe of massive stars, with mass M ≳ 8 M_⊙, via successive neutron capture on preexisting Fe isotopes, both before and during the explosion (Timmes et al. 1995; Wanajo et al. 2013). Other suggested sources are asymptotic giant branch (AGB; Busso et al. 2003) and super-AGB (SAGB; Lugaro et al. 2012) stars, as well as a rare class of high-density thermonuclear SNe (Woosley 1997).

Owing to its comparatively long half-life of t_{1 / 2} ≃ 2.6 Myr (Rugel et al. 2009; Wallner et al. 2015), ejected ⁶⁰Fe survives long enough not only to travel significant distances (a few hundred pc) through the ISM, but also to be detectable by its β⁻ decay via 60Co and gamma-ray emission at 1173 and 1333 keV. Measurements with the INTEGRAL satellite (Wang et al. 2007) revealed an average photon flux per gamma-ray line of about 4.4 × 10^-5 cm^-2 s^-1 rad^-1 for the inner Galaxy region, which might be translated into a steady-state Galactic ⁶⁰Fe mass of 1.53 M_⊙ (Diehl 2016). Very recently, Binns et al. (2016) succeeded in directly detecting ⁶⁰Fe nuclei in Galactic cosmic rays, showing that ⁶⁰Fe sources must have been within the distance the high-energy particles can travel for the duration of t_{1 / 2}, which is typically less than about 1 kpc.

However, the half-life of ⁶⁰Fe is very short when compared to the age of the solar system (≃4.5 Gyr), implying that the ⁶⁰Fe, which, as we know from the analysis of primitive meteorites (Tachibana & Huss 2003; Mostefaoui et al. 2005; Bizzarro et al. 2007), was present at about that time, has completely decayed away. This adds to the fact that the isotope’s natural terrestrial background, arising mainly from the influx of pre-enriched extraterrestrial material (e.g. dust and meteorites) onto the Earth, and to a lesser extent from cosmogenic spallative production due to the steady cosmic-ray flux into the Earth’s atmosphere, as well as in situ production via the penetrating muon and neutron flux, is low (Fields & Ellis 1999).

These considerations were motivated by the work of Knie et al. (1999) who had, for exactly these reasons, used accelerator mass spectrometry to search for an ⁶⁰Fe excess in a ferromanganese (FeMn) crust from the South Pacific Ocean floor and indeed found one in two layers dating back about 6 Myr. In a more detailed analysis of another FeMn crust sample (termed 237KD), stemming from the equatorial Pacific, Knie et al. (2004) narrowed down the ⁶⁰Fe anomaly to the time interval 1.7–2.6 Myr before present (cf. Feige et al. 2012). This definitive signal was corroborated in another study of 237KD conducted by Fitoussi et al. (2008), who also looked into a different geological reservoir: marine sediments. These grow about a thousand times faster than FeMn crusts and thus have a far better time resolution, albeit the isotopic concentration of each layer is lower. Unfortunately the sediment core analysed by Fitoussi et al. (2008), which originates from the North Atlantic Ocean and spans the same time interval as the crust 237KD, did not show a compatible ⁶⁰Fe signal.

In contrast, Wallner et al. (2016) found a signal consistent with the expected time from the FeMn crust study in four sediment cores from the Indian Ocean. Remarkably, the signal is broader, rather pointing towards an accumulation of SNe than a single event. Apart from the sediment cores, Wallner et al. also analysed two FeMn crusts (237KD not included) and two FeMn nodules stemming from the Pacific and Atlantic oceans, respectively. Since each of these deep-sea archives contain elevated levels of ⁶⁰Fe for the time range 1.5–3.2 Myr ago, the authors further concluded that the ⁶⁰Fe signal must be global. For one of the two FeMn crusts they also detected an additional smaller ⁶⁰Fe peak 6.5–8.7 Myr ago, whose exact origin is yet rather elusive. Remarkably, a collision of asteroids in the main belt 8.3 Myr ago, connected to a boosted, possibly ⁶⁰Fe enriched, influx of interplanetary dust particles and micrometeorites (Farley et al. 2006) falls within this particular time range.

Returning briefly to the subject of marine sediments it should be added that these also harbour so-called magnetotactic bacteria, which gather iron from their environment to create nano-sized crystals of magnetite (Fe₃O₄). The latter are used by the microbes to orient themselves within Earth’s magnetic field, and thus to navigate to their preferred conditions – a behaviour known as magnetotaxis. Bishop & Egli (2011) speculated that nearby SN acitivity should have left its mark also in the fossilised remains of such bacteria. In order to prove this hypothesis, they analysed magnetofossils of a sediment core from the eastern equatorial Pacific. The results (Ludwig et al. 2016) indeed show a comparatively weaker ⁶⁰Fe signal in layers spanning a time range of about 1 Myr, with the maximum again located at ≃2.2 Myr ago.

In view of all these discoveries it seems little surprising that several of the lunar rocks that were recovered in the course of the Apollo missions also showed an enhanced concentration of ⁶⁰Fe (Fimiani et al. 2016). Unfortunately, time-resolved studies are impossible as our cosmic neighbour is regularly hit by a broad range of impactors (Langevin & Arnold 1977), which continuously churn up its surface and mix different layers. Besides this so-called gardening, there is also almost not sedimentation on the Moon. Lunar samples can hence only provide a first hint for the existence of a signature (cf. Feige et al. 2013).

Theoretical attempts to link the FeMn crust measurements to the formation of the LB, but also to a single SN event, were made by Feige (2010) and Fry et al. (2015), respectively, both on the basis of analytical calculations. Feige particularly used an SN model developed by Kahn (1998), which does not assume a homogeneous density distribution for the blast wave to expand in, but rather a medium that had already been modified by a previous SN event. More specifically it was assumed that the first SN runs into a homogeneous ISM, whereas all subsequent SNe occur within the spherical region excavated by the first SN, while being exposed to a radial density profile of the form ρ(R) ∝ R^{9 / 2}. Following Fuchs et al. (2006), an IMF with one star per mass bin was used, and the most probable explosion sites were calculated. Including isotopic yields from stellar evolution calculations, the measured ⁶⁰Fe profile could be reproduced with surprisingly high accuracy.

Results of such analytical models should however be treated with some caution, since the similarity solutions applied have a number of serious deficiencies. First, the external pressure is required to be small compared to the pressure in the expanding bubble. This poses no strong limitations for the modelling of a second or subsequent SN event, whose blast wave takes only little time to move through the hot and thin pre-stratified medium. However, the situation is completely different for a single SN, or the very first one in a series, whose blast wave must penetrate through a totally inhomogeneous ISM and thus would soon undergo strong distortions. By the same token, the density distribution of the ambient medium has to be either homogeneous or obey a power law, implying that possible dense obstacles are not taken into account. This includes the LB’s outer shell, which was neglected in the preliminarily work by Feige (2010), but was later implemented as a Weaver et al. (1977) wind bubble for the improved analytical study contained in Breitschwerdt et al. (2016, hereafter BFS+). Second, the external medium is taken to be constant over time. In reality, however, the medium changes after each SN. Feige (and also BFS+) at least switched from a homogeneous to a power-law density distribution after the first SN went off. Nevertheless, they assumed that all subsequent explosions occur in the same environment. Last but not least, turbulent mixing and mass loading, as present in real SBs, are not taken into account in those models.

A definitive answer to the question whether the observed ⁶⁰Fe excess is explainable in the context of the LB formation can therefore only be given through three-dimensional high-resolution numerical simulations. This paper details on our first step working towards that goal (see also Schulreich 2015). For this purpose we built upon the numerical simulations of BFS+, presenting an additional homogeneous model (see Table 1, model A). As it will be shown this sets a lower limit for possible homogeneous media. We particularly used the masses of the same 16 SN progenitor stars, together with the explosion times, total ejected masses, and ⁶⁰Fe mass fractions derived from there. The applied SN explosion sites, on the other hand, result from the most probable stellar trajectories (see also BFS+). Like Breitschwerdt & de Avillez (2006), we studied the formation and evolution of the LB not in isolation, but in concert with the neighbouring SB Loop I. For the sake of conciseness, we currently neglected the complex multiphase nature of the (local) ISM and considered instead only homogeneous background media, including the one already used in BFS+. Deeper conclusions can be drawn now, however, as we not only traced the global ⁶⁰Fe dynamics in our solar neighborhood but also the contributions of the individual SNe. This was achieved by tagging each of them with their own passive scalar (or tracer) – a quantity that behaves like a drop of dye when dispersed in a liquid. The influence of the inhomogeneous SN-driven environment, already teased in BFS+, will be extensively analysed in the forthcoming Paper II, which will also contain a discussion of the properties of turbulence in such a self-consistently evolved medium.

The article is organised as follows. Section 2 describes the model setup, in Sect. 3 our results with respect to SB evolution and comparison to the FeMn crust measurements are discussed, and we close with a general discussion and conclusions in Sect. 4.

2. Model description

Our HD simulations were carried out with the tree-based AMR code Ramses (Teyssier 2002), which allows for solving the discretised Euler equations in their conservative form by means of a second-order unsplit Godunov method for perfect gases. In particular, we employed the MUSCL-Hancock scheme (van Leer 1979) together with the Local Lax-Friedrichs (or Rusanov) approximate Riemann solver (Lax 1954; Rusanov 1961) and the MinMod slope limiter (Roe 1986). The computational domain is cubic with 3 kpc side length and covered by a base Cartesian grid of 64³ uniform cells. In regions of steep density and pressure gradients, the grid was recursively refined up to six additional levels, resulting into a peak spatial resolution of 0.7 pc. This high resolution was also kept at the coordinate origin, which marked the location of Earth (or Sun; for a distinction one would require several additional levels of grid refinement, and thus much more computing time). The grid was initially filled with a uniform and static medium. The adopted values for the gas temperature T, particle density n, and metallicity (in solar units) Z/Z_⊙ in our two best-fit models are given in Table 1. They actually mimic two specific “phases” of the classical three-phase ISM model by McKee & Ostriker (1977), namely the warm ionised medium (model A) and the warm neutral medium (model B). The SB-generating SNe were initialised as blast waves right at the start of their Sedov-Taylor phase, each liberating a canonical explosion energy of E_SN = 10⁵¹ erg. The self-gravitating gas was exposed to an isotropic radiation field (that includes absorption through the Galactic plane as well as the cosmic microwave background), but was also allowed to loose energy via various interstellar cooling processes. This energy gain and loss was modelled under the assumption of collisional ionisation equilibrium (CIE) using the temperature-dependent net cooling functions generated for several gas densities with the spectral synthesis code Cloudy (Ferland et al. 1998). The boundaries of our computational domain were formally periodic, but were actually never reached by expanding shells during the time frame considered. As a further “boundary condition”, we set up a neighbouring SB, intended to represent Loop I. This is required since ROSAT observations have shown that soft X-rays are absorbed by a nearby neutral “wall” enclosed by an even denser ring-like structure, which is most likely the result of an interaction between the LB and the Loop I SB (Egger & Aschenbach 1995).

2.1. Setting up the Loop I superbubble

For identifying possible SN progenitor stars of Loop I, we performed a search for B stars analogous to the one of the LB (Fuchs et al. 2006), but for a spherical search volume of twice the size. Since less kinematical data is available for stars at larger distances, it becomes increasingly difficult to find bright stars. It turns out that Tr 10 and the association Vel OB2 have recently passed through the present volume of Loop I rather off-centre. Like for the LB this is not a problem per se, as the extension and morphology of evolved SBs are primarily determined by the ambient density and pressure distribution.

Given a sample of 80 stars that might have entered this volume about 12.3 Myr ago (main sequence lifetime of a still existing 8.2 M_⊙ star), we first needed to estimate how many have already exploded. For that purpose, we fitted an IMF that is typical for young massive stars with initial masses M, $\frac{\mathrm{d𝒩}}{\mathrm{d}\mathit{M}}\mathrm{=}{\frac{\mathrm{d𝒩}}{\mathrm{d}\mathit{M}}|}_{\mathrm{0}}{\mathit{M}}^{\mathrm{\Gamma }\mathrm{-}\mathrm{1}}\mathit{,}$(1)where Γ = −1.1 ± 0.1 (Massey et al. 1995; Berghöfer & Breitschwerdt 2002; Fuchs et al. 2006). With the lower and upper limit of the relevant mass range being defined by A0 and B0 stars, respectively (i.e. M_l = 2.6 M_⊙ and M_u = 8.2 M_⊙), we obtained the normalisation constant via integration of Eq. (1); $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\frac{\mathrm{d𝒩}}{\mathrm{d}\mathit{M}}|}_{\mathrm{0}}\mathrm{=}\mathrm{𝒩}{\left[{\mathrm{\int }}_{\mathrm{2.6}}^{\mathrm{8.2}}{\mathit{M}}^{-2.1}\hspace{0.17em}\mathrm{d}\mathit{M}\right]}^{-1}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{=}\mathrm{80}{\left[\frac{{\mathrm{2.6}}^{-1.1}\mathrm{-}{\mathrm{8.2}}^{-1.1}}{\mathrm{1.1}}\right]}^{-1}\mathrm{=}\mathrm{351.}\end{array}& & \end{array}$(2)The lifespan of the first and most massive star that exploded in Loop I might be expressed as ${\mathit{\tau }}_{{\mathit{M}}_{{\mathit{t}}^{\mathrm{\prime }}}}\mathrm{=}{\mathit{\tau }}_{\mathrm{u}}\mathrm{-}{\mathit{t}}^{\mathrm{\prime }}\mathit{,}$(3)where τ_u is the lifetime of a star of mass M_u and t′ is the time at which the association entered the present volume of Loop I. When quantifying the main sequence lifetime of stars in the mass range 2 ≤ M ≤ 67 M_⊙ by a fit to the isochrone data of Schaller et al. (1992), $\mathit{\tau }\mathrm{=}{\mathit{\tau }}_{\mathrm{0}}{\mathit{M}}^{\mathrm{-}\mathit{\beta }}\mathit{,}$(4)where τ₀ = 1.6 × 10⁸ yr and β = 0.932 (Fuchs et al. 2006), Eq. (3) translates into ${\mathit{M}}_{{\mathit{t}}^{\mathrm{\prime }}}\mathrm{=}{\left({\mathit{M}}_{\mathrm{u}}^{\mathrm{-}\mathit{\beta }}\mathrm{-}\frac{{\mathit{t}}^{\mathrm{\prime }}}{{\mathit{\tau }}_{\mathrm{0}}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathit{\beta }}\mathrm{·}$(5)As t′ = 12.3 Myr, the mass of the most massive star that has already exploded is M_t′ = 19.2 M_⊙. One therefore finds for the number of missing stars $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathrm{𝒩}}_{\mathrm{SN}}\mathrm{=}{\mathrm{\int }}_{\mathrm{8.2}}^{\mathrm{19.2}}{\frac{\mathrm{d𝒩}}{\mathrm{d}\mathit{M}}|}_{\mathrm{0}}{\mathit{M}}^{-2.1}\hspace{0.17em}\mathrm{d}\mathit{M}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{=}\mathrm{351}\left[\frac{{\mathrm{8.2}}^{-1.1}\mathrm{-}{\mathrm{19.2}}^{-1.1}}{\mathrm{1.1}}\right]\mathrm{=}\mathrm{19.}\end{array}& & \end{array}$(6)Obviously the masses of these SN progenitors should be based on the IMF introduced above. In order to obtain the statistically most probable distribution (Maíz Apellániz & Úbeda 2005), we performed the mass binning by assigning exactly one star to each mass interval of the IMF, where the actual stellar masses are taken to be the average mass in each bin. Analogously, Feige (2010) (and also BFS+) derived the initial masses of the LB stars, which we also employed in this work. Using Eq. (1), we obtained for the mass of a particular star lying in the mass bin between M₁ and M₂ (>M₁), $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{M}}_{\mathrm{\ast }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{M}}_{\mathrm{1}}\mathrm{+}{\mathit{M}}_{\mathrm{2}}\mathrm{)}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\begin{array}{c}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}{\mathit{M}}_{\mathrm{1}}\mathrm{+}{\left[{\mathit{M}}_{\mathrm{1}}^{-1.1}\mathrm{-}\mathrm{1.1}{\left({\frac{\mathrm{d𝒩}}{\mathrm{d}\mathit{M}}|}_{\mathrm{0}}\right)}^{-1}\right]}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{1.1}}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\mathit{.}\end{array}& & \end{array}$(7)These initial masses served as a basis for determining the total ejected masses (using stellar evolution data from Woosley & Heger 2007) as well as the ⁶⁰Fe mass fractions (see Sect. 2.2). The explosion times, on the other hand, are defined by t_exp = τ−τ_c, where τ is as given in Eq. (4) and τ_c is the age of the cluster, which should be comparable to the lifespan of a 8.2 M_⊙ star, i.e. τ_c = 23 Myr. For arranging the stellar explosion centres in space, we calculated the trajectories of all known B stars of Tr 10 and Vel OB2. Since not all of the member stars pass through the Loop I volume we selected for the 19 SN progenitors only the closest to the centre. A full list of the Loop I input parameters is given in Table A.1.

2.2. Calculating the amount of SN-released ⁶⁰Fe that arrives on Earth

According to stellar evolution models (e.g. Limongi & Chieffi 2006; Woosley & Heger 2007), there is a close link between the ⁶⁰Fe content of a star and its mass. For the considered stellar mass range of ≃10–30 M_⊙, the calculated ⁶⁰Fe yields however exhibit a large scatter between ≃10^-6 and 10^-4M_⊙. This is due to the sensitivity of nucleosynthesis to a variety of factors, including cross-sections, mass loss, and rotation. For our model, we neglected these subtleties and instead used the exponential fit to recent stellar evolution data provided in BFS+.

Owing to its low concentration, ⁶⁰Fe should have no dynamical influence on the (turbulent) flow. We therefore treated it as a passive scalar, obeying an advection-diffusion equation of the form (see e.g. Davidson 2004) $\frac{\mathit{\partial Z}}{\mathit{\partial t}}\mathrm{+}\mathrm{(}{u}\mathrm{·}{\nabla }\mathrm{)}\mathit{Z}\mathrm{=}\mathit{\alpha }{\mathrm{\nabla }}^{\mathrm{2}}\mathit{Z,}$(8)where Z is the ⁶⁰Fe mass fraction, u the fluid velocity, and α the diffusivity of the contaminant (which is assumed here to be isotropic). Since for the case of the ISM the Peclet number, Pe = UL/α, is large (U and L denote characteristic velocity and length scales), diffusive effects are restricted to the very small, so-called microscale of turbulence. At larger scales, the right-hand side of Eq. (8) can hence be neglected. The scalar Z then acts like a marker that tags the chemically enriched fluid parcels. However, there is still a need for some diffusive process ultimately operating at small scales, even though its exact value is quite insignificant for the turbulent mixing at the larger scales, as demonstrated by de Avillez & Mac Low (2002) in ISM simulations with different resolutions. Like them we did not implement any physical diffusion term, but rather left the mixing to numerical diffusion, which is generally faster and operates on a larger scale than its physical counterpart. As a consequence, the timescale of mixing in our simulations rather poses a lower limit to the mixing time resulting from physical diffusion.

When generalising the analysis to the spatiotemporal evolution of m passive scalars with concentrations Z_i(i = 1,...,m), the time-dependent Euler equations in conservation-law form read as (see e.g. Toro 2009) $\frac{\mathit{\partial }{U}}{\mathit{\partial t}}\mathrm{+}\frac{\mathit{\partial }{F}\mathrm{\left(}{U}\mathrm{\right)}}{\mathit{\partial x}}\mathrm{+}\frac{\mathit{\partial }{G}\mathrm{\left(}{U}\mathrm{\right)}}{\mathit{\partial y}}\mathrm{+}\frac{\mathit{\partial }{H}\mathrm{\left(}{U}\mathrm{\right)}}{\mathit{\partial z}}\mathrm{=}{S}\mathrm{\left(}{U}\mathrm{\right)}\mathit{,}$(9)with the state vector ${U}\mathrm{=}{\left[\mathit{\rho ,\rho u,\rho v,\rho w,E,\rho }{\mathit{Z}}_{\mathrm{1}}\mathit{,}\mathit{...}\mathit{,\rho }{\mathit{Z}}_{\mathit{m}}\right]}^{\mathrm{T}}\mathit{,}$(10)the vectors of fluxes in the corresponding coordinate directions $\begin{array}{ccc}{F}& \mathrm{=}& {\mathrm{[}\mathit{\rho u,\rho }{\mathit{u}}^{\mathrm{2}}\mathrm{+}\mathit{P,\rho uv,\rho uw,u}\mathrm{(}\mathit{E}\mathrm{+}\mathit{P}\mathrm{)}\mathit{,\rho u}{\mathit{Z}}_{\mathrm{1}}\mathit{,}\mathit{...}\mathit{,\rho u}{\mathit{Z}}_{\mathit{m}}\mathrm{]}}^{\mathrm{T}}\mathit{,}\\ {G}& \mathrm{=}& {\mathrm{[}\mathit{\rho v,\rho uv,\rho }{\mathit{v}}^{\mathrm{2}}\mathrm{+}\mathit{P,\rho vw,v}\mathrm{(}\mathit{E}\mathrm{+}\mathit{P}\mathrm{)}\mathit{,\rho v}{\mathit{Z}}_{\mathrm{1}}\mathit{,}\mathit{...}\mathit{,\rho v}{\mathit{Z}}_{\mathit{m}}\mathrm{]}}^{\mathrm{T}}\mathit{,}\\ {H}& \mathrm{=}& {\mathrm{[}\mathit{\rho w,\rho uw,\rho vw,\rho }{\mathit{w}}^{\mathrm{2}}\mathrm{+}\mathit{P,w}\mathrm{(}\mathit{E}\mathrm{+}\mathit{P}\mathrm{)}\mathit{,\rho w}{\mathit{Z}}_{\mathrm{1}}\mathit{,}\mathit{...}\mathit{,\rho w}{\mathit{Z}}_{\mathit{m}}\mathrm{]}}^{\mathrm{T}}\mathit{,}\end{array}$and S(U) denoting an appropriate source term vector. Here, ρ is the fluid density, u, v, w the components of the velocity vector u, E the total energy, and P the (thermal) pressure. Imprinted in this evolution is the radioactive decay of ⁶⁰Fe, for which we also accounted in our simulations. The quantitative comparison between our model results and the FeMn crust measurements was then accomplished via the following steps.

First, we determined the local interstellar fluence, ℱ, which is the flux of the isotope at the Earth’s position integrated over time. Since this position coincides with the point of contact of eight, highest resolved grid cells, with the HD variables being always defined at the cell centres, we had to average over the ⁶⁰Fe mass fluxes, ρ | u | Z, in those cells. Hence, $\mathrm{ℱ}\mathrm{=}\frac{\mathrm{\left(}\mathit{\rho }\mathrm{\right|}{u}\mathrm{|}\mathit{Z}{\mathrm{)}}_{\mathrm{VA}}}{\mathit{A}{\mathit{m}}_{\mathrm{u}}}\mathrm{\Delta }\mathit{t,}$(14)where (... )_VA denotes the volume average, A = 60 the mass number of ⁶⁰Fe, m_u ≃ 1.66 × 10^-24 g the atomic mass unit, and Δt the last simulation time step.

For the remaining steps we took the fall-out on Earth to be isotropic, which was demonstrated by Fry et al. (2016) to be a suitable assumption. The surface density of ⁶⁰Fe atoms deposited at time t before present in the FeMn crust can then be expressed as¹ (cf. Fields & Ellis 1999) $\mathrm{\Sigma }\mathrm{=}\frac{\mathit{fU}}{\mathrm{4}}\mathrm{ℱ}\exp \mathrm{(}\mathrm{-}\mathit{\lambda t}\mathrm{)}\mathit{,}$(15)where the factor ¹/ 4 stems from relating Earth’s cross-section (${\mathit{R}}_{\mathrm{\oplus }}^{\mathrm{2}}\mathit{\pi }$) to its surface area ($\mathrm{4}{\mathit{R}}_{\mathrm{\oplus }}^{\mathrm{2}}\mathit{\pi }$), the exponential term describes the decay of the radionuclide during its final disposal in the FeMn crust (λ ≡ ln(2) /t_{1 / 2}), and the product fU represents the ⁶⁰Fe survival fraction.

Specifically, f denotes the fraction of ⁶⁰Fe atoms that reaches the Earth’s orbit after overcoming a variety of filtering processes, while being condensed into dust grains. After estimating the dust condensation of ⁶⁰Fe at departure from source, as well as its destruction through interstellar passage, heliosphere, and solar radiation pressure, Fry et al. (2015) concluded that f ≃ 0.01.

The other, so-called uptake factor, U, accounts for the fact that only a certain fraction of ⁶⁰Fe spread over Earth’s surface actually finds its way into the FeMn crust. The underlying chemical selection process is only poorly understood, implying that the value of U is, as of f, very uncertain. Knie et al. (2004) quantified U using the relative concentrations of iron and manganese in water and the FeMn crust, and the uptake of manganese (4%), while implicitly assuming that f = 1. They found fU ≡ U ≃ 0.006. More recent studies, however, suggest that U = 0.5–1 (Bishop & Egli 2011; Feige et al. 2012). Hence, by taking U = 0.6 and f as derived by Fry et al. (2015), one retains the numerical value claimed by Knie et al. (2004), but still adds some (very simplistic) sub-grid dust physics to the problem.

We finally obtained the number density of ⁶⁰Fe for each crust layer by summing Σ in the time bins lying within the corresponding time intervals and then dividing by the thickness of the layer. Measurements demand to relate these densities to the one of stable iron, i.e. n_⁶⁰Fe/n_Fe, or, in short, ⁶⁰Fe/Fe. The required normalisation constant follows from (cf. Feige 2014) ${\mathit{n}}_{\mathrm{Fe}}\mathrm{=}\frac{{\mathit{w}}_{\mathrm{Fe}}\mathit{\rho }{\mathit{N}}_{\mathrm{A}}}{\mathit{A}}\mathit{,}$(16)where w_Fe ≃ 0.1527 is the weight fraction of iron in the FeMn crust sample whose mass density ρ is 1.5 g cm^-3, N_A ≃ 6.022 × 10²³ mol^-1 is the Avogadro constant, and A ≃ 55.845 g mol^-1 is the molar mass for iron. We thus have n_Fe ≃ 2.470 × 10²¹ cm^-3.

2.3. Overview of input data

To conclude this section we give a full list of the input data and recapitulate how it is obtained:

• sample of cluster stars: selected according to compactnessin real and velocity space from astrometric catalogue data(such as Hipparcos) for a heliocentric sphere with diameterD ~ 400 pc (LB) and 800 pc (Loop I);

• cluster age (τ_c): derived by comparison of the cluster turn-off point with isochrones of Schaller et al. (1992);

• number and initial masses of SN progenitors (, M_∗): estimated via IMF fitting of sample stars (Massey et al. 1995), assuming one star per mass bin with the bin’s average mass for the perished stars (Eqs. (6) and (7));

• main sequence lifetimes of SN progenitors (τ): calculated as a function of mass only via Eq. (4);

• SN explosion times (t_exp): calculated from subtracting τ_c from τ, thus assuming coeval star formation in the cluster;

• SN explosion centres (x,y,z): derived from the center-of-mass stellar trajectories, as determined via the epicyclic equations of motion (Lindblad 1959; Wielen 1982); for the LB, the most probable ones are calculated for a Gaussian error distribution in the Hipparcos positions and proper motions;

• total ejecta mass and ⁶⁰Fe mass fractions (M_ej, Z): taken, corresponding to M_∗, from stellar evolution models (Woosley & Weaver 1995; Rauscher et al. 2002; Limongi & Chieffi 2006; Woosley & Heger 2007);

• ⁶⁰Fe survival fraction (fU): quantified by combining estimates from a dust survival model (Fry et al. 2015, and references therein) and measurements in terrestrial archives (Knie et al. 2004; Bishop & Egli 2011; Feige et al. 2012); the dynamics of dust and gas is assumed to be perfectly coupled, where the latter is traced via passive scalars.

3. Results

3.1. Evolution and properties of the Local and Loop I superbubbles

Fig. 1 Colour-coded maps of the logarithmic total gas column density in model A for three different times before present. Integrations along the y- and z-axes are shown in the upper and lower panels, respectively. Earth’s projected position is marked by a cross. Circles indicate the projected centres of SN explosions that occurred in the time frame before the snapshot was taken. The sizes of the circles correspond to the initial masses of the progenitor stars (see legend in the upper left panel).

Fig. 2 As for Fig. 1, but for model B.

Fig. 3 Colour-coded maps of the logarithmic 60Fe mass density in model A at three different times before present. Vertical (y = 0) and horizontal (z = 0) cuts through the computational box are shown in the upper and lower panels, respectively. Symbols are as in Figs. 1 and 2.

Fig. 4 As for Fig. 3, but for model B.

The SN explosions generate a coherent bubble structure whose evolution is followed over 12.6 Myr (age of the LB) until the present time. In Figs. 1 and 2 we present snapshots of the total gas column density for the models A and B, respectively. Note that in these, as well as in all remaining maps, the x, y, and z-axes point towards the Galactic centre, into the direction of the Galactic rotation, and towards the Galactic north pole, respectively.

After their almost coeval formation, the two SBs evolve at first completely independently from each other. The shocked ambient medium cools very quickly and gets compressed into thin shells. These supershells are subject to thermal and fluid dynamical instabilities, which give rise to the formation of cold and dense clumps. In particular, during the “breaks” between two consecutive SN explosions, when the supershell is decelerating, Rayleigh-Taylor instabilities develop at the contact discontinuity that separates the shocked SN ejecta from the cooler and denser swept-up ambient medium. The filaments of dense gas produced thereby grow into the cavity where they can mix with the hot bubble gas. In addition, an overstability sets in that arises from a mismatch between thermal and ram pressure at the shell surface. This “Vishniac instability” (see Ryu & Vishniac 1991, and references therein) can cause growing, oscillating ripples on that surface. Although such structures are particularly prominent for extinct SBs, whose shell motion has stalled due to the lack of energy sources (i.e. stellar winds and SNe; cf. Krause et al. 2013), they may already be visible at the late stages of the simulations presented in this work; see e.g. the spike protruding from the shell of Loop I at (x,y) = (240,−600) pc in the lower right panel of Fig. 2. It is also interesting to note a criss-cross network of shocks as a result of sequential SN explosions. This may lead to stronger coalescing shock waves and thus to local pressure enhancements, which may also be responsible for local deviations from sphericity of the LB and Loop I despite the homogeneous ambient medium.

Certainly an effect of the latter is that the supershells retain their spherical shape on large scales. After about 3 Myr (t ≃ 9.6 Myr before present; model A) and 4.6 Myr (t ≃ 8 Myr before present; model B) the shells collide, giving rise to an even denser and cooler interaction layer. Also this interface becomes eventually Rayleigh-Taylor-unstable, since the pressure in the two SB volumes is not equal. This leads to the formation of cold and dense cloudlets that travel into the direction of lower pressure (i.e. the LB interior) as soon as the instability becomes fully non-linear (cf. Breitschwerdt et al. 2000). In model A, the interaction shell breaks up after an evolution time of about 6.5 Myr (t ≃ 6.1 Myr before present), allowing for immediate exchange of hot gas between the SBs. In model B, on the other hand, no break-up occurs during the whole evolution time considered. The final size of the LB is about 800 by 600 pc (model A) and 580 by 480 pc (model B) in the (x,y)-plane, and about 760 pc (model A) and 540 pc (model B) perpendicular to it. The atomic hydrogen number density and temperature of the LB interior, particularly in the solar neighbourhood, is about 10^-4.2–10^-3.9 cm^-3 and 10^6.9–10^7.1 K in model A, and about 10^-4.2–10^-3 cm^-3 and 10^5.8–10⁷ K in model B. If the ambient medium is homogeneous, it is generally the first SN that determines the final appearance of the SB most. This is because within a few sound-crossing times the bubble has a homogeneous pressure acting on the shell of the first explosion. The centre of the final SB should consequently lie very close to the explosion centre of the earliest SN. That this is indeed the case, also in our models, can be seen best in the bottom panels of Figs. 1 and 2 (the largest black circle marks the locations of the earliest explosion). Furthermore, since we used the same SN data for both models, the relative arrangement of the two SBs does not change either.

It should be noted that the calculated final x- and y-extensions only poorly match the sizes derived from observations. Also the elongation in z-direction is highly underestimated. This is because the background medium is assumed to be homogeneous, which is certainly an oversimplification of the actual conditions in the solar neighbourhood. As seen from the extended, three-dimensional differential opacity maps of the local ISM by Lallement et al. (2014), the LB is bound today by a series of dense clouds located in the first and fourth quadrants, as well as in the anti-centre region of the Galactic plane. These clouds, which also possess non-negligible vertical extensions, are known to be part of Gould’s belt. Moreover, there is a 500–1000 pc wide cavity in the third quadrant, as well as a small, elongated cavity in the opposite direction. The larger cavity is centred below the Galactic plane and was identified as the counterpart of the so-called GSH 238+00+09 supershell, which is visible in radio wavelengths (see Puspitarini et al. 2014, and references therein). Assuming that these structures were present in a similar fashion as observed today throughout the entire evolution time, the LB’s outer shell would probably have experienced early on a greater “resistance” in the second, fourth, and (partly) the first quadrant, and would thus not have grown to such a great extent in those directions. This should particularly hold for model A, in which the SBs expand more rapidly. The evolution in the third quadrant is more difficult to predict since in our models, this region is partly occupied by Loop I. The presence of the dense interstellar cloud at a distance of ≃200 pc and longitude 240° might nevertheless restrict the LB growth in a similar fashion as the developing interaction layer in our simulations. The combined effect of all these structures might also “shift” the solar system closer to the centre of the then smaller LB. Further knowledge can only be gained from additional simulations whose initial conditions are tailored to match these recent observations. These will be the subject of a forthcoming paper.

Current results however suggest that the exact size and shape of the LB does not play a great role for modelling the ⁶⁰Fe amount that arrives on Earth, except for the case in which ⁶⁰Fe is carried to the solar system by a reflected shock (see Sect. 3.2). Generally it is necessary that the solar system lies within the LB so that it can be reached by the expanding SN remnants. For this, the gas-dynamical properties of the medium between the SN explosion centres and the locus of the Earth are crucial, as they determine the arrival time of the LB’s outer shell and hence set the date for the earliest possible ⁶⁰Fe signal in the terrestrial record.

We further note that the relative locations of the Local and the Loop I SB are not the same as in the joint-evolution model by Breitschwerdt & de Avillez (2006). This is because they used different stellar trajectories for the LB progenitors, and, more importantly, a different formation scenario for Loop I, based on other moving group stars from the Sco-Cen cluster for which they assumed a single explosion centre that is somewhat arbitrarily shifted 200 pc relative to the centre of the LB in positive x-direction. Therefore, and contrary to our model, their LB is generated by 19 SNe, and their Loop I SB by 39 SNe. Apart from that there is some discussion in the literature whether the large circular feature in the radio continuum sky, known as the North Polar Spur, which is associated with Loop I, is part of the local ISM at all or rather belongs to the Galactic centre region (see Puspitarini et al. 2014, and references therein). The presence of nearby young stars in the direction of the Galactic centre, however, strongly argues for SN explosions that generated the Loop I SB.

Using a passive scalar as discussed in Sect. 2.2, we followed the spatiotemporal evolution of the SN-released ⁶⁰Fe. Corresponding simulation snapshots are shown in Figs. 3 and 4 for models A and B, respectively. From these maps it is seen that inhomogeneities in the ⁶⁰Fe distribution, arising from recent SNe, are smoothed out over time. Responsible for this mixing process are turbulent shear flows resulting from SN blast waves that run into the surrounding supershell and thus generate asymmetric reflected shock waves. The zones of underpressure left behind by these reflected discontinuities manifest themselves as the blue “lakes” particularly visible in the lower panels of Figs. 3 and 4. With a characteristic large-eddy size of ℓ ≃ 100 pc and an average sound speed after an SN of the order of a ≃ 10² km s^-1, we obtained for the mixing timescale within the SB cavity, τ_m ~ ℓ/a, a value of about 1 Myr. Considering that the last SN occurred roughly 1.5 Myr ago, it is therefore not surprising that the ⁶⁰Fe mass is distributed quite uniformly in the LB cavity at present time.

3.2. Comparison with crust measurements

At each simulation time step, the total local interstellar fluence of ⁶⁰Fe atoms was calculated according to Eq. (14). Figure 5 shows the resulting profiles for both models. Upon closer inspection (see inlay), the signals come in three different types, which are all embedded in a “background noise” of ≃10⁵ cm^-2 amplitude.

Fig. 5 Local interstellar fluence of 60Fe atoms due to all SNe in the LB volume as a function of time before present in model A (grey line) and B (black line). The inlay magnifies the three different enhancement types, being SN blast waves (high and narrow signals), reflected shocks (weaker and broader signals), and the supershell (very broad rightmost signal).

The first type are sharp and intense (ℱ ≃ 0.3–2.4 × 10⁶ cm^-2) sawtooth waves that occur seven times in model A, but only once in model B. Such a wave form is characteristic of the density profile of an SN remnant in the Sedov-Taylor phase. These signals must hence originate from the shells of individual remnants that cross the Earth’s orbit while expanding into the LB volume. The time of direct deposition of the SN ejecta on Earth, given by the (scale) width of the signal, is a measure of the remnant’s shell thickness. The time range of direct exposure in our models, Δt ≃ 70–130 kyr, roughly lies between the values proposed by Fields et al. (2005) and Bishop & Egli (2011), but almost perfectly matches the predictions given in Feige (2014), which are based on the SN model by Kahn (1998).

The second type, which appears to be unique to model A, are signals that occur after almost every sawtooth wave (note that the time axis is from right to left), while being always more extended and weaker than the waves they succeed, almost as if they were the SN explosions’ “echoes”. Remarkably, the time lag between the two signals increases the closer the present time is approached.

The third type is a joint feature of both profiles, occurring only once and right at the beginning (t ≃ 6.1 and 2.2 Myr in model A and B, respectively). Here, the signals are quite broad (Δt ≃ 300 and 450 kyr in model A and B, respectively) and possess intermediate amplitudes of about 2.9 × 10⁵ (model A) and 9.6 × 10⁵ cm^-2 (model B).

Fig. 6 Comparison of the measured (symbols with error bars) and simulated (histograms) 60Fe/Fe ratio as a function of the terrestrial archive’s layer age for model A (upper panel) and B (lower panel). Simulation data is set upon the instrumental background derived specifically for 237KD, which is indicated by the dashed line. The grey levels of the histograms refer to the 60Fe survival fractions employed.

Applying the rebinning procedure described in Sect. 2.2 allowed us to do a direct comparison with the FeMn crust measurements of Knie et al. (2004) and Fitoussi et al. (2008), as shown in Fig. 6. The ⁶⁰Fe survival fraction, fU, was either set to 0.006 (light grey histograms), as originally suggested by Knie et al., or 0.005 (dark grey histograms), corresponding to a combination of a dust fraction of f = 0.01 with the newly derived lower limit for U (cf. Sect. 2.2). It is seen that for both models, measurements are best matched if the slightly lower uptake factor is employed. This is particularly true for the (global) maximum in the crust layer corresponding to t ≃ 2.2 Myr, whose presence was lately confirmed by the study of Wallner et al. (2016, see grey symbols in Fig. 6, and which is perfectly reproduced by both models (max(⁶⁰Fe / Fe) ≃ 2.3 × 10^-15). It is also interesting to note that model A suggests contributions in five more layers than model B, with two of those, namely for t ≃ 3.1 and 5.7 Myr, lying within the error bars of the measurements. The local maximum calculated for t ≃ 3.9 Myr is however absent in the FeMn crust data available.

Fig. 7 Local interstellar fluence of 60Fe atoms due to each individual SN in the LB volume (see legends) as a function of time before present in models A (upper panel) and B (lower panel). The inlays are zoom-ins on the last group of peaks, respectively.

To explore the exact origin of each signal in the fluence profiles, we rerun the simulations with all parameters unchanged, except that we assigned a passive scalar to each individual SN event. Formally, this implies adding an index i with 1 ≤ i ≤ 16 to the variables ℱ, Z, and Σ in Eqs. (14) and (15). Results are plotted in Fig. 7, where identification numbers correspond to the explosion order. Inspecting these profiles, together with the left column panels of Figs. 3 and 4, reveals that the earliest signals of both models, i.e. the quite broad peaks of intermediate amplitude already known from Fig. 5, are produced when the LB’s supershell runs over the solar system. Separated from the shocked ambient medium by a contact discontinuity, the shell harbours a large amount of ⁶⁰Fe atoms that have by then been expelled in previous SN events (SNe #01–08 and #01–15 in model A and B, respectively) and that have not yet decayed (see also the inlays of Fig. 7). Since the supershell of model B arrives later at the Earth’s orbit, it could already accumulate more SN ejecta (and thus enriched material), which is the reason why the profile of the total local interstellar fluence (Fig. 5) has a first signal that is both broader and higher than in model A. Also note that the first signals do not set in abruptly but in a smooth, ramp-like manner. This is due to the combined effects of numerical diffusion present in all Godunov-like advection schemes and the Rayleigh-Taylor instability operating at the contact discontinuity inside the supershell (see Sect. 3.1), allowing for mixing across the discontinuity.

Even the pulses generated by single SN remnants actually entrain fractions of previously released ⁶⁰Fe, since the average time between successive explosions is short when compared to the half-life of ⁶⁰Fe. This is nicely confirmed by Fig. 7, showing that the high peaks associated with recent SN events indeed superpose a series of peaks with lower amplitudes. The same can be observed for the weaker follow-up signals, which are produced when the SN blast waves return after being reflected from the surrounding supershell; so in this sense they are indeed SN “echoes”. Finally, since the supershell moves away from Earth, the time lag between an explosion and its “echo” becomes greater for later times.

Putting all this information together, we arrive at the following interpretation of Fig. 6. Starting with the most prominent feature, the (global) maximum in the layer for t ≃ 2.2 Myr, we can now safely state that, in model A, it is produced by SN #14 and 15 and their reflected shocks, whereas in model B, it is due to yet undecayed material from SNe #01–15 reaching Earth almost simultaneously while being gathered within the LB’s supershell (see also the middle and left column panels of Figs. 3 and 4, respectively). The additional local maximum at t ≃ 3.9 Myr, which is unique to model A, is mainly due to SN #12 and 13. The upper crust layers (corresponding to t ≲ 1.7 Myr) obtain their ⁶⁰Fe content almost exclusively from SN #16. There is at first the arrival of the SN blast wave itself (seen for model B in the middle column panel of Fig. 4) that generates the enhancement in the layer for t ≃ 1.3 Myr². Later on, turbulent motions induced by the blast wave of the last SN lead to the enhancement in the layer for t ≃ 0.4 Myr. Particularly, all the “background noise” in the modelled fluence profiles is due to turbulence inside the LB cavity.

We hence conclude that the measured ⁶⁰Fe excess in the FeMn crust could be either due to the passages of dust-loaded shells from nearby SN remnants over the Earth’s orbit, sweeping up also enriched material, which was present at that time in the LB volume (model A), or, alternatively, due to the arrival of the supershell of the LB that incorporates the material from several previous SN events (model B). We consider the LB formation scenario of model B to be physically more probable than that of A. Not only are the derived and measured ⁶⁰Fe / Fe ratios in better agreement (there is particularly no second maximum), but also the size and (more importantly for simulations that are based on homogeneous background media) the gas-dynamical properties of the fully developed SB are a better match to observations.

We further emphasise that the additional smaller peak at t ≃ 6.5–8.7 Myr, partly present in the data of Wallner et al. (2016, filled triangles in Fig. 6, does not occur in our profiles. If this enhancement is not of meteoritic origin, as mentioned earlier, but actually connected to the formation of the LB, it would be an indication for a clustering of stars between 12 and 15 M_⊙, thus differing from the most probable binning used for the present study. However, this hypothesis needs further testing and simulations, and is beyond the scope of this work, but will be addressed in a forthcoming paper.

4. Discussion and conclusions

In this paper we presented three-dimensional HD simulations (with resolutions adaptively adjusted down to subparsec scale) aiming to reproduce an ⁶⁰Fe excess measured in a deep-sea FeMn crust in the context of the formation of the LB. Parameters for the LB-generating SNe are based on an IMF normalised to still existing low mass stars of the parent stellar moving group identified by Fuchs et al. (2006). Like in the analytical model by Feige (2010) (see also BFS+), ejected masses, lifetimes and thus explosion times were derived from the initial masses of the progenitor stars, whereas computations of the most probable paths for the individual perished members of the moving group provided most probable explosion sites. For estimating the particular ⁶⁰Fe mass ejecta, we resorted to predictions from stellar evolution models. The formation of the LB, as well as of the neighbouring Loop I SB, was studied in two homogeneous, self-gravitating environments, sharing similarities with the classical warm ionised (model A) and warm neutral medium (model B). Interstellar heating and cooling processes were included via standard CIE net cooling functions. By utilising passive scalars for tagging the ⁶⁰Fe enriched gas released into the SB cavities, we recorded the flux of ⁶⁰Fe atoms at the Earth’s orbit (usually termed local interstellar fluence) over the whole simulation time, and then calculated time averages over a range matching the crust layers in order to do comparisons with the measurements of Knie et al. (2004) and Fitoussi et al. (2008). Our models are able to reproduce both the observed timing and intensity of the ⁶⁰Fe excess with rather high precision. However, the underlying physical processes are different. In model A the signal arises due to the fast-paced blast waves of two individual SN remnants (#14 and 15 out of 16 in total), which cross the Earth’s orbit twice as a result of reflection from the LB’s outer shell. In model B, on the other hand, it is the supershell of the LB itself that injects the yet undecayed ⁶⁰Fe content of all previous SNe (#01–15 in this case) at once, but over a longer time range. For both scenarios it is important to note that ⁶⁰Fe does not reach Earth in gaseous form, but condensed onto dust grains that have to survive a variety of filtering processes throughout their journey. These were summarised by us into the combined multiplicative factor fU, for which we found a better agreement with the crust data if fU = 0.005 rather than 0.006. Both values occur in the current literature. The observed extension of the LB as well as its gas-dynamical properties were better matched by model B. The age of the LB has been estimated from the main sequence lifetime of a still existing 8.2 M_⊙ star (cf. Fuchs et al. 2006).

Primarily for the sake of better comparison with analytical calculations (Feige 2010, BFS+) our model does currently not include the stellar winds blown by the moving group members. For a locally homogeneous ISM, as considered in this work, these winds would have carved an extended, egg-like shaped cavity around the cluster’s center-of-mass trajectory, already before the first stars explode (cf. Weaver et al. 1977, and references therein). As soon as this happens, however, the dynamics of the wind-blown bubble is quickly overtaken by the SNe, as is clear from comparing the mechanical luminosities of the winds and SNe (cf. Tomisaka et al. 1981). The former cannot be larger than $\frac{\mathrm{1}}{\mathrm{2}}\mathit{Ṁ}{}_{\mathrm{tot}}\mathit{v}_{\mathrm{\infty }}^{\mathrm{2}}\mathrm{\simeq }\mathrm{2.5}\mathrm{\times }{\mathrm{10}}^{\mathrm{36}}\mathrm{erg}\hspace{0.17em}{\mathrm{s}}^{-1}$, where Ṁ_tot ≃ 2 × 10^-6M_⊙yr^-1 is the maximum stellar mass loss rate (only reached if the most massive star is still alive) and v_∞ ≃ 2000 km s^-1 is the terminal wind velocity. The latter, on the other hand, is E_SN/ Δτ ≃ 4.2 × 10³⁷ erg s^-1, with Δτ ≃ 7.4 × 10⁵ yr being the average time interval between two consecutive SNe. Note that we only consider cumulative mechanical energy output rates here, which should be justified because the mean separation between the massive stars (and hence the SN explosion centres) of ≃44 pc is most of the time smaller than the size of the SB. So the expansion of mature SBs is determined by SNe, rather than by stellar winds (and also ionising photons), which dominate the early phase of evolution (see also Rogers & Pittard 2013; Krause et al. 2013; Yadav et al. 2016; Kim & Ostriker 2017; Kim et al. 2017).

Although our simulations draw a quite conclusive picture of the LB formation scenario, there are nevertheless some caveats that have to be addressed.

First, the production and expulsion of ⁶⁰Fe is a complicated process that strongly depends on the mixing processes in the interiors of massive stars, and the depth of their convection zones. Predictions of the ⁶⁰Fe yields can thus vary by a factor of a few, even for the same stellar masses, depending on which stellar evolution model is used (cf. Tur et al. 2010).

Second, the uptake factor, quantifying the sedimentation through the Earth’s atmosphere onto the ocean floor, is only loosely constrained. More precise estimates, either from measuring ⁶⁰Fe at other sites or from the analysis of other radioisotopes, are therefore highly welcome, since its value has a profound effect on our quantative results.

The third concerns heliospheric filtering. ⁶⁰Fe-enriched plasma could only reach the solar system via the poles because of strong ram pressure of the solar wind in the equatorial plane. Due to its higher mass this does not affect ⁶⁰Fe-dust. A comprehensive description of the formation and survival of SN dust at different compositions and sizes is therefore desirable, as it is particularly questionable whether the single factor applied in this work can fully account for all processes involved. In principle, similar arguments apply for the previous two points. An obvious advantage of employing simple input parameters however lies in the fact that they can be easily adapted when the underlying physics is better understood or measurements have become more precise. Still, it should be emphasised that these uncertainties affect the absolute values of ⁶⁰Fe, but not relative ones, giving the model a solid foundation.

Last but not least, there is the dependence on the background medium, as well as on the explosion sites in the stellar moving group. For the former we can say from our current studies with a homogeneous background medium that the average external density may not exceed ≃0.3 cm^-3. Otherwise, the LB supershell would arrive later and the ⁶⁰Fe overabundance would no longer appear in the crust layer in which it was discovered. Yet strictly speaking, this is only true for the SN explosion centres applied in this work, which we left unchanged as these are the most probable ones for the considered stellar moving group (see also BFS+). One can however check whether any relevant changes in the calculated ⁶⁰Fe/Fe profiles (Fig. 6) would result if an SN would occur by a maximum of 1σ away from its most probable explosion site lying at distance d. In order to estimate this spread in distance, δR, we traced back the paths of 1000 hypothetical stars that would explode on exactly the same location 10 Myr into the past, while taking into account the error in their mean velocities (≃3 km s^-1 for each component). We obtained δR ≃ 20 pc. As the blast wave velocity in the LB volume is of the order of 10⁴ km s^-1, this translates to a spread in the arrival time of a thus generated ⁶⁰Fe signal of ≃2 kyr, which lies far below the average time resolution of the FeMn crust sample (≃0.6 Myr, see the bin sizes in Fig. 6). Relevant for the height of the signal is the size of the SN shells when they are crossing our solar system. When assuming for the sake of simplicity perfect sphericity and neglecting the clumpiness of these shells as a result of instabilities, the factor by which the height of the signal would change lies between (1 ± δR/d)². For the last three SNe in the LB volume d ≃ 100 pc, yielding for the factor’s lower and upper limit 0.64 and 1.44, respectively.

Directly related to the background medium is the question about the amount of ⁶⁰Fe in the region where the LB began to form. The steady-state Galactic ⁶⁰Fe mass of 1.53 M_⊙ given by Diehl (2016) was derived for the diffuse, extended ISM and thus does not allow conclusions to be drawn on local values in the solar neighborhood today, and even less on those of more than 10 Myr ago, where also radioactive decay should play a crucial role. The only serious way to obtain such information would be long-term chemical evolution simulations of the (local) Galaxy that particularly include chemical mixing. As such models are currently unavailable, we set in our simulations the initial amount of ⁶⁰Fe equal to zero everywhere, thus always representing the lower limit case. We note in passing that a completely analogous problem arises in the context of explaining the origin of ⁶⁰Fe in the early solar system (see e.g. Wasserburg et al. 1998; Huss et al. 2009; Boss & Keiser 2013; Vasileiadis et al. 2013), demonstrating that such studies are strongly needed.

Also worth mentioning is the alternative explanation for the 2.2 Myr-old ⁶⁰Fe signal given by Basu et al. (2007) and Stuart & Lee (2012). Their idea is based on the fact that ⁶⁰Fe is not only produced in SNe but also by spallation reactions between Galactic cosmic rays and nickel in (micro-) meteorites and interplanetary dust particles. Via collisions of asteroids, e.g. in the main belt, these small objects could have been released in large amounts, thus exposing Earth to a boosted isotopic influx. However, the concentration of nickel measured within the FeMn crust is too low to explain the ⁶⁰Fe signal within the framework of this scenario (Fitoussi et al. 2008; Wallner et al. 2016). On the other hand, the cosmic-ray protons, antiprotons, and positrons that would have been copiously released by a nearby SN would certainly have led to a temporal increase of radioisotopes in meteoritic particles. As a matter of fact, there are evidences for one (or more) such events about 2 Myr ago in the cosmic ray energy spectra (e.g. Kachelrieß et al. 2015), fitting nicely into the picture drawn in this work.

Another strong hint for the SN origin of the ⁶⁰Fe signal are the experimental ⁶⁰Fe and ⁵³Mn activities recently obtained from lunar samples (Fimiani et al. 2016). Like ⁶⁰Fe, the long-lived radionuclide ⁵³Mn is generated from cosmic-ray spallation reactions on iron in meteorites. Hence, if the ⁶⁰Fe signal was produced from an enhanced influx of micrometeorites, one would expect the lunar activity ratios of ⁶⁰Fe/⁵³Mn to be consistent with meteoritic values. Yet the lunar soil samples carry an enhanced ⁶⁰Fe signature that exceeds the activity created by cosmic rays.

We can actually test our numerical calculations against these lunar soil measurements, stating an observed fluence in the range between 10⁷ and 6 × 10⁷ cm^-2. As their data is not time-resolved, we have to sum the value of Σ (cf. Eq. (15)) in all time bins available, spanning a total range of 12.6 Myr. When setting U = f = 1, we obtain fluences of 6.18 × 10⁸ cm^-2 for model A and 3.36 × 10⁸ cm^-2 for model B. While it is almost certain that f< 1 (cf. Sect. 2.2), the Moon’s uptake factor should be, due to the lack of an atmosphere, equal or very close to unity. Taking our simulations as a basis, we estimate a range of 0.016 and 0.179 for the ⁶⁰Fe survival fraction fU. The lower limit would be actually compatible to an uptake of 100%, if the value of f were only slightly higher than estimated by Fry et al. (2015), namely 0.016 instead of 0.010. Therefore also the measurements of ⁶⁰Fe in lunar samples are in agreement with our previously performed simulations.

Acknowledgments

M.M.S. and D.B. acknowledge funding by the DFG priority program 1573 “Physics of the Interstellar Medium”. We thank the anonymous referee for her/his valuable comments which helped to improve the manuscript.

References

Appendix A: Additional table

All Tables

All Figures

Fig. 1 Colour-coded maps of the logarithmic total gas column density in model A for three different times before present. Integrations along the y- and z-axes are shown in the upper and lower panels, respectively. Earth’s projected position is marked by a cross. Circles indicate the projected centres of SN explosions that occurred in the time frame before the snapshot was taken. The sizes of the circles correspond to the initial masses of the progenitor stars (see legend in the upper left panel).

In the text

	Fig. 2 As for Fig. 1, but for model B.
In the text

	Fig. 3 Colour-coded maps of the logarithmic 60Fe mass density in model A at three different times before present. Vertical (y = 0) and horizontal (z = 0) cuts through the computational box are shown in the upper and lower panels, respectively. Symbols are as in Figs. 1 and 2.
In the text

	Fig. 4 As for Fig. 3, but for model B.
In the text

	Fig. 5 Local interstellar fluence of 60Fe atoms due to all SNe in the LB volume as a function of time before present in model A (grey line) and B (black line). The inlay magnifies the three different enhancement types, being SN blast waves (high and narrow signals), reflected shocks (weaker and broader signals), and the supershell (very broad rightmost signal).
In the text

Fig. 6 Comparison of the measured (symbols with error bars) and simulated (histograms) 60Fe/Fe ratio as a function of the terrestrial archive’s layer age for model A (upper panel) and B (lower panel). Simulation data is set upon the instrumental background derived specifically for 237KD, which is indicated by the dashed line. The grey levels of the histograms refer to the 60Fe survival fractions employed.

In the text

	Fig. 7 Local interstellar fluence of 60Fe atoms due to each individual SN in the LB volume (see legends) as a function of time before present in models A (upper panel) and B (lower panel). The inlays are zoom-ins on the last group of peaks, respectively.
In the text