© The Authors 2026

1 Introduction

Recent advances in observational methods allow us to map out the structure in the nearby galactic interstellar medium (ISM) in 3D (e.g., Arenou et al. 1992; Lallement et al. 2019; Edenhofer et al. 2024). One of the most prominent structures, mapped by these techniques is the Local Bubble (hereafter LB Cox & Reynolds 1987; Linsky & Redfield 2021), a large cavity several hundred parsecs across, centered around the Solar System (Zucker et al. 2022, hereafter Z+22), which has been linked to diffuse soft X-ray emission (e.g., Snowden et al. 2000; Yeung et al. 2024). Superbubbles (SBs) such as the LB are believed to be carved out by the various feedback channels of massive stars, such as stellar winds (Tenorio-Tagle et al. 1990), ionizing radiation, (Linsky & Redfield 2021) and supernovae (hereafter SNe Tenorio-Tagle et al. 1990; Breitschwerdt & de Avillez 2006; Wallner et al. 2021).

Over the past few decades, many studies have tried to constrain the origin of the LB. Studies using astrometry, tracing back the positions of nearby star-clusters suggest that ∼ 10−20 SNe originating from the Scorpius-Centaurus OB association (Sco-Cen) might have contributed to the expansion of the LB ∼ 10 Myr ago (Maíz-Apellániz 2001; Fuchs et al. 2006; Zucker et al. 2022). Numerical simulations of SNe expanding into a turbulent, stratified ISM, confirm that ∼ 20 SNe exploding sequentially throughout the last ∼ 14 Myr could explain the observed size of the LB, the observed column densities of O VI in the interior, and the deposition of sedimentary radioisotopes, such as ⁶⁰Fe (Breitschwerdt & de Avillez 2006; Breitschwerdt et al. 2016; Schulreich et al. 2023). Fossil records show that the incorporation rate of ⁶⁰Fe and ²⁴⁴Pu on earth has peaked ≲ 4 Myr and ≲ 7 Myr ago (Wallner et al. 2016, 2021). Combining this with the work of Z+22, who find that the Solar System would have entered into the LB ∼ 5 Myr ago, suggests that at least the more recent peak might coincide with SNe associated with the LB (see also Breitschwerdt et al. 2016).

In (Romano et al. 2025a, hereafter Paper I), we introduce the SISSI (Supernovae In a Stratified, Shearing ISM) zoom-in simulation project, focusing on the properties of highly resolved, simulated supernova remnants (SNRs) expanding into the self-consistently generated ISM of an isolated, Milky-Way-like disk-galaxy. In Paper I, we find that SNRs in a realistic environment, expand faster than previously expected at later times t ≳ 1 Myr, due to galactic shear, the gravitational influence of nearby substructures and the presence of low-density channels. A comparison of our simulated sample of SNRs at t=10 Myr with the LB, suggested that it would be too small for a SB of its age, powered by SNe exploding at the rate suggested by previous studies.

In this work, we follow up on this curiosity, by presenting a new analysis of the LB, derived from the 3D dust maps of Edenhofer et al. (2024) (hereafter E+24) using a method similar to that of O’Neill et al. (2024) (hereafter O+24). As opposed to Z+22, who assume that the stellar kinematics trace the dynamics of the gas through triggered star formation, we directly estimate the size, momentum and ambient density of the LB from the 3D dust maps of E+24 and by utilizing the dynamical evolution of the SNRs in the SISSI simulation (see Appendix A for a brief description of the simulation suite) we constrain its age and the number of SNe powering its evolution.

2 The age of the LB

The properties of a SB are determined by the strength of its energy source, the density of the environment and its age. Previous studies suggest that the radial momentum imparted per SN onto an expanding radiative shell only weakly depends on the density and is largely independent of the age of the SB (Walch & Naab 2015; Oku et al. 2022); namely, by knowing the momentum, we can precisely determine the number of SNe responsible for the expansion of the SB. In Paper I we have found that the average momentum imparted per SN is approximately $${\hat{p}}_{\mathrm{S}\mathrm{N}}\sim (2.6\pm 0.4)\times {10}^5{n}_0^{-0.13}{\mathrm{M}}_{\odot }\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1},$$(1) where n_H=n₀ cm⁻³ is the ambient number density of hydrogen and we estimate the systematic uncertainty based on the range of values found in numerical studies (e.g. Walch & Naab 2015; Kim & Ostriker 2015; Martizzi et al. 2015), dependent on the detailed cooling physics, resolution and other numerical choices. In a purely energy-driven model that can arise in situations in which cooling is inefficient (Kahn 1998), p̂_SN increases slightly with time (El-Badry et al. 2019; Zucker et al. 2022), though for the range of ages considered in the SISSI simulations these differences are minor.

In Appendix B, we find the LB’s average ambient density of ∼ 0.5 cm⁻³, an age (t_age)-dependent momentum of ∼ 3 × 10⁷(Myr/t_age M_⊙ km s⁻¹), a mass of ∼ 6.2 × 10⁵ M_⊙ and a size of ∼ 212.3 pc (see Table B.1 for uncertainties and a comparison with previous estimates). Using these estimates we derive a momentum input per SN of p̂_{SN, LB} ∼(2.9 ± 0.5_sys ± 0.005_stat) × 10⁵ M_⊙ km s⁻¹. We related the number of SNe to its age by utilizing the age-dependent momentum estimate $$\begin{array}{rl}{N}_{\mathrm{S}\mathrm{N},\mathrm{L}\mathrm{B}}\left(t_{\mathrm{a}\mathrm{g}\mathrm{e}}\right)& =p_{\mathrm{L}\mathrm{B}}/{\hat{p}}_{\mathrm{S}\mathrm{N},\mathrm{L}\mathrm{B}}\\ & \sim (104\pm 2_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}}\pm {59}_{\mathrm{s}\mathrm{y}\mathrm{s}})\times {\left(t_{\mathrm{a}\mathrm{g}\mathrm{e}}/\mathrm{M}\mathrm{y}\mathrm{r}\right)}^{-1}\end{array}$$(2)

The corresponding average SN rate was obtained by dividing by the age, Ṅ_{SN, LB}=N_{SN, LB}/t_age. We validated the applicability of this estimate in Appendix C.2.

In contrast to the momentum, the size of a SB depends only weakly on the energy input, and in turn depends most strongly on its age. In Fig. 1, we show the time evolution of the effective size of our simulated sample of SNRs, introduced in Paper I, for SNRs with ambient densities within 0.3 dex of that of the LB. The SBs evolving in a realistic galactic environment (solid lines) grow significantly faster than what is expected from blastwave models in uniform environments (see Appendix D El-Badry et al. 2019; Oku et al. 2022). The size of the SB driven by subsequent SNe (N1x10) closely matches that of the SNR of a single SN (N1) for t ≲ 3 Myr, while it approaches that of the SNR of 10 SNe exploding all at once (N10) after ∼ 10 Myr. For the models with multiple SNe, the size of the LB is reached after ∼ 4.5 Myr (N10) and after 6−7 Myr (N1x10).

At the age at which the simulated SNRs (N1x10) reach the size of the LB, 17 ± 10_sys SNe are required to explain its observed momentum (Eq. (2)), in contrast to the 6-7 SNe that exploded in the simulation. The difference arises from the fact that the simulations were not fine-tuned in order to match the observationally based momentum of the LB. A SB powered by ∼ 17 SNe would have reached the same size even earlier than in the simulation making the difference even large in LB age between our evaluation and previous age estimates (Maíz-Apellániz 2001; Breitschwerdt & de Avillez 2006; Zucker et al. 2022). While the simulations are not fine-tuned to the observation, one can get a rough estimate of the difference using theoretical considerations.

Previous models of radiative blastwaves in uniform ambient media (El-Badry et al. 2019; Oku et al. 2022; Lancaster et al. 2024) suggest that the effective size is roughly $\propto {\dot{N}}_{\mathrm{SN}}^{0.15-0.25}$. A blastwave powered by more frequent SNe would thus only lead to a mild increase in the size and would be largely consistent with an age of ≲ 6 Myr. In Fig. 2 we show that for R ∝〈Ṅ_SN〉^α with α ≲ 0.225 the size matches that of the LB for t_age ∼ 3.5−5.5 Myr, while for α>0.225, there are no solutions that can simultaneously explain the LB’s size and momentum. Not shown are scenarios with extreme Ṅ_SN, but we confirm that they reproduce a similar range of ages (∼ 2.8−5.8 Myr) and correspondingly lower and higher numbers of SNe ranging from ∼ 7 to ∼ 59 SNe. A more detailed evaluation is outside of the scope of this work and will be subject to future investigation (Romano et al., in prep.).

By invoking standard assumptions about the stellar initial mass function, the average time between SNe in a star cluster can be linked to its mass (Kim et al. 2017): $$\Delta t_{\mathrm{SN}}\sim 0.4{\left(\frac{M_{\mathrm{cl}}}{{10}^4{\mathrm{M}}_{\odot }}\right)}^{-1}\mathrm{Myr},$$(3) which for the range of SN rates found above would correspond to a mass of M_{cl, LB} ∼ 0.5−8.4 × 10⁴ M_⊙ for the progenitor cluster of the LB.

We compared this number to the swept-up mass, and assuming that the swept up mass roughly corresponds to the mass of the birth cloud of the cluster driving the expansion of the LB (i.e. neglecting the background ISM) we obtained an upper limit for the cloud-scale star formation efficiency, $${ϵ}_{\star }=\frac{M_{\mathrm{cl},\mathrm{LB}}}{M_{\mathrm{sw},\mathrm{LB}}}\lesssim 0.8-14\%,$$(4) in agreement with the range of values found in observational (Lee et al. 2016; Chevance et al. 2020) and numerical studies (Grudić et al. 2022; Farias et al. 2024).

We summarize our results and compare them to those of Z+22 in Table 1. The obtained cluster mass is somewhat higher than the stellar mass of Sco-Cen (Krause et al. 2018), which is oftentimes assumed to be the progenitor cluster of the LB (Maíz-Apellániz 2001; Zucker et al. 2022). Similarly, our result that the LB should be younger than was previously believed is in tension with the older age (∼ 15−17 Myr) of the stars in Sco-Cen, indicating that Sco-Cen stars alone might not alone power the LB.

An additional constraint for the SN origin of the LB comes from the fossil records of sedimentary radionuclides (Wallner et al. 2021; Ertel et al. 2023; Ertel & Fields 2024), which requires that the Solar System to have entered the LB at least ≳ 4 Myr ago. Both the scenario proposed in this work (Appendix E) and that of Z+22 satisfy this constraint. A more detailed analysis of the dynamical incorporation of the dusty, radioactive material may provide additional evidence in favor of either scenario.

According to Swiggum et al. (2024) (hereafter, S+24), Sco-Cen is part of a larger family of stars known as the α-Persei family (α Per), which reportedly had an average time between SNe of $\Delta t_{\alpha \mathrm{Per}}=0.{32}_{-0.23}^{+0.50}\mathrm{Myr}$. Moreover, the star formation history (SFH) of αPer indicates episodic star formation over the past ≳ 60 Myr, with peaks every ≳ 12.5 Myr and the latest peak within the last ∼ 12.5 Myr.

Our results are in agreement with the SFH in αPer if we assume that the LB was powered by the latest peak of star formation rather than the trough associated with the formation of Sco-Cen. In Appendix F we constructed a more detailed star formation and SN-rate history over the past 40 Myr to substantiate this claim. The SN rate and SFH resulting from these considerations are shown in Fig. F.2. The upper panel shows that the peak of star formation in the LB area is before the onset of SNe that fuel its expansion, followed by a sudden distinct downturn, coincident with the range of ages found here. This suggests that the expansion of the LB may have predominantly quenched star formation in the solar neighborhood, rather than triggering it. While LB-triggered star formation is not entirely ruled out, especially when one considers young star-forming regions on the LB’s edge, it appears to be subdominant.

In Fig. F.3, we show the orbits of the clusters potentially contributing SNe to the LB. We find that while the majority of the active clusters are coincident with Sco-Cen, in line with the expectation that the LB likely originated from Sco-Cen (Fuchs et al. 2006), there were likely significant contributions from other stellar associations in the solar neighborhood.

The high expansion speed v_LB ∼ 0.5 R_LB/t_age ≳ 10 km s⁻¹ (Table 1), suggests that the LB might be just above the regime where it is starting to be affected by its Galactic environment, which typically dominates as v ≲ σ_turb ∼ 10 km s⁻¹ (Paper I). Thus, models for SBs expanding into a uniform medium (Appendix D) might provide a reasonable age estimate after all. For the slowly-cooling wind model used by Z+22 with an updated ambient density and the SN rate following S+24, we obtained $t_{\text{age}}^{\mathrm{SCW}}\sim 6.6_{-2.3}^{+2.4}\mathrm{Myr}$, while for a rapidly-cooling wind model (Oku et al. 2022) we obtained $t_{\text{age}}^{\mathrm{RCW}}\sim 9_{-4}^{+5}\mathrm{Myr}$, in rough agreement with both scenarios.

de Avillez & Breitschwerdt (2012) compare the local ISM simulations of the LB and the adjacent Loop I SB to the observed ion abundances of C IV, N V, and O VI to constrain the expansion history of the LB. They find that, within the scope of their model of the local ISM, the last SN explosion within the LB must have occurred 0.5−0.8 Myr ago. These numbers are in rough agreement with the lower end of SN rates found here (low N_SN column of Table 1). They are somewhat larger than the SN rates, based on recent Gaia data ($\sim 0.{32}_{-0.23}^{+0.50}\mathrm{Myr},\mathrm{S}+24)$. One reason for this discrepancy could be the slightly higher average ambient density of ∼ 1 cm⁻³ in the simulations of de Avillez & Breitschwerdt (2012), leading to shorter overall chemistry timescales. In a lower-density environment, the same constraint might thus be in better agreement with the higher SN rates found here. A detailed comparison is outside of the scope of this work and will be subject to future investigation.

The LB appears to be at the tipping point between microscopic SNRs, whose dynamics are dominated by local shock physics, and mesoscopic SBs that are strongly affected by galactic-scale processes. More detailed models of such large SBs as well as more precise stellar ages (e.g., Miret-Roig et al. 2022; Ratzenböck et al. 2023; Swiggum et al. 2024) are needed to shed light on the progenitor cluster and expansion history of the LB.

Fig. 1 Time evolution of effective size of the simulated sample of SNRs of the SISSI sample between 2 and 10 Myr, for SNRs with ambient densities within 0.3 dex of the LB (Table B.1). Gray, red and blue lines correspond to different explosion models. Solid lines correspond to the median size of the simulated bubbles, with shaded areas corresponding to the range between the 30th and 70th percentiles. Dotted lines and the hatched blue contour correspond to theoretical models based on radiative blastwaves in uniform media (Appendix D). For model N10, the size of the LB corresponds to an age of ∼ 4.5 Myr, while for N1x10 it corresponds to ∼ 6−7 Myr.

Fig. 2 Effective size as a function of age or number of SNe, respectively, extrapolated from the model N 1 x 10 and the mean of the momentum constraint Eq. (2) for various values of α; namely, Rα=RN 1 × 10(t=tage) ×[ṄSN, LB(tage)/1 Myr−1]α, where ṄSN, LB= NSN, LB/tage. For α ≲ 0.225, the extrapolated sizes match the LB for ages in the range 3.5−5.5 Myr, while for α ≳ 0.225 the two constraints on the size and momentum cannot be satisfied simultaneously.

3 Concluding remarks

We compared the momentum and the effective size of the LB, derived from publicly available, high-quality 3D dust maps, to our sample of simulated SNRs expanding into the shearing, stratified ISM of the isolated Milky-Way-like SISSI galaxy. We find that in order to match constraints on both the radial momentum and the effective size, the LB would have to be driven by a larger number of SNe and be significantly younger than was previously reported.

In particular, a rough agreement requires an age since the onset of SNe of t_{age, LB} ∼ 2.8−5.8 Myr for N_{SN, LB} ∼ 7–59, in agreement with estimates of the average time between SNe in αPer (S+24). Investigating the SFH of the solar neighborhood (Fig. F.2), such a young age indicates that SNe-feedback-induced formation of the LB overall quenched local star formation although new clouds forming at its shell might now trigger new generations (Z+22). Another possibility is that the nearby star-forming regions (e.g., Sco-Cen, Taurus) were already formed before the formation of the LB. Indeed, these regions are all closer (Zucker et al. 2021) than the effective radius of the LB (∼ 210 pc), indicating that they might be entrained by it (see e.g., Fig. 13 of Romano et al. 2024). This might trigger further star formation. While our current methodology, which defines the bubble’s boundary as the first prominent density peak along each line of sight, cannot disentangle these details, a more sophisticated method (e.g., O’Neill et al., EAS 2025) might provide more insight.

We stress that the conclusions of this work that star formation near the sun might be quenched and driven (Z+22) by the expansion of the LB are not mutually exclusive. Star formation on large scales can be quenched, while locally the conditions for new star formation are triggered. The degree to which SN-driven SBs such as the LB, might trigger and quench star formation on large (≳ 100 pc) scales is a highly interesting question that is not well understood. The SISSI simulations are well suited for this task and an investigation of this pertinent question is currently underway (Romano et al., in prep.).

We also not that while our analysis focuses on SNe only, a fraction of the energy is also provided from other sources such as stellar winds and H II regions. Though, it is worth stressing that Schulreich et al. (2023) consider the role of stellar winds and find that before the onset of SNe the LB maintains an approximately constant size of ∼ 30 pc for ∼ 10 Myr, suggesting that the contribution from winds might be minor. In this context the age of the LB could be inflated by prolonged periods of negligible growth. We do not count such prolonged periods of halted growth toward the age of the LB.

While our results are in tension with previous estimates, which assumed that the LB was exclusively powered by SNe in Sco-Cen, they are in better agreement with more recent observational data on the more global star formation in the solar neighborhood. We have used scaling relations in order to extrapolate our simulations to the unique conditions of the LB. More detailed models for old (t ≳ 1 Myr) SBs expanding into a realistic galactic ISM as well as additional observational constraints are now required to resolve remaining uncertainties and obtain a clearer picture for the origin and history of the LB.

Data availability

The Julia source-code for our analysis is made available at https://doi.org/10.5281/zenodo.17054923.

Acknowledgements

We thank Manuel Behrendt for detailed and very fruitful discussions and comments that substantially improved the quality of the manuscript and for providing an early version of the AVALON galaxy simulation code which is the basis for the SISSI simulations. We also thank the anonymous referee for their insightful comments and suggestions that helped to improve the quality of this work. This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2094 – 390783311. LR thanks the developers of the following software and packages that were used in this work: Julia v1.10.0 (Bezanson et al. 2017), Matplotlib v3.5.1 (Hunter 2007), Healpix v2.3.0 (Tomasi & Li 2021), and Galpy v.1.10.2 (Bovy 2015).

References

Appendix A The SISSI simulations

The SISSI simulations are a suite of highly resolved (Δ x ≲ 1 pc) hydrodynamical simulations of SNe exploding in thirty distinct regions of the multi-phase ISM of an isolated Milky-Way-like galaxy. The ISM of the SISSI galaxy is shaped by the complex interplay of stellar feedback, realistic radiative cooling and heating (M. Behrendt et al. in prep.), secular evolution in an axisymmetric background gravitational potential, and local gravitational collapse and star formation.

As such the SISSI simulations represent a natural next step in terms of physical realism for models of SNRs in environments of increasing realism and complexity (e.g. Kim & Ostriker 2015; Fielding et al. 2018; Kim & Ostriker 2017). Due to the large range of probed physical conditions they complement simulations of specific SNRs such as the LB (e.g. Breitschwerdt & de Avillez 2006; de Avillez & Breitschwerdt 2012; Schulreich et al. 2023) by drawing general conclusions about SNR evolution in a realistic ISM.

In SISSI we considered three separate explosion scenarios: (i) N1 corresponds to a single SN explosion at each site at the beginning of the zoom-in period (t=0), (ii) N10 corresponds to 10 simultaneous SN explosions at each site at t=0, and (iii) N 1 × 10 corresponds to a single SN at each site, every 1 Myr starting from t=0. The zoom-in period begins after a prolonged ISM-generation period at lower resolution, where stellar feedback is modeled in an effective way. During the zoom-in period lasting 10 Myr, the effective stellar feedback is shut-off in order to not interfere with the zoomed-in SNe.

Some key features of the SISSI simulations include

1. Ambient densities in the range 10⁻³ ≲ n_H[cm⁻³] ≲ 10³.

2. A multi-phase ISM with a typical gas velocity dispersion of ∼ 10 km/s with significant scatter.

3. Explosion sites at three distinct galactocentric radii (R ∼ 2, 4.5, and 8 kpc, corresponding to orbital timescales of ∼ 50−200 Myr.

Here we recapitulate some of the key differences between SNRs exploding in the ISM of the SISSI galaxy versus SNRs exploding in uniform media, found in Paper I:

1. Strong deviations from spherical symmetry (Minor-to-major-ratios ≲ 0.2, semi-major-to-major-ratios ≲ 0.4).

2. A systematically larger size.

These differences generally manifest after a few Myr or equivalently a few percent of an orbit, likely due galactic shear, the galactic density structure and the coupling to pre-existing blow-out regions. At the low number of SNe considered in SISSI, SBs are not expected to break out of the disk significantly, unless they expand through already carved out low-density channels (Schulreich et al. 2023; Romano et al. 2025b).

The parameter space covers the physical conditions of the LB (e.g., ambient density, geometry), allowing us to draw tentative conclusions about its origin, under the assumption that it is not too different from a typical SNR in a similar environment from the SISSI simulation suite.

Appendix A.1 Comparison to simulations of the local ISM

There are a number of studies studying the properties and history of the Local Bubble, using dedicated simulations of SNe in the local ISM (Breitschwerdt & de Avillez 2006; de Avillez & Breitschwerdt 2012; Schulreich et al. 2023). In these studies, a stratified ∼ 1 kpc × 1 kpc ISM patch, resolved with a maximum resolution of ∼ 1 pc, and physical conditions roughly matching those in the local ISM is set up before the stellar feedback input (SNe and in more recent studies stellar winds) from a sample of nearby stellar groups, associated with the present-day Sco-Cen association are allowed to perturb the medium and carve out a SB identified with the LB.

These studies are in excellent agreement with various observables in the local ISM, such as the column-density of O VI and other ions (de Avillez & Breitschwerdt 2012) and the deposition of radio-isotopes in the Solar System (Breitschwerdt et al. 2016). These simulations suggest an age of the LB of 10−15 Myr since the onset of SN-explosions in Sco-Cen.

Despite these successes it should be noted, that several simplifications have been made that might affect their results.

1. The simulations neglect the self-gravity of the gas, which might affect the structure of the ISM, the orbits of the stellar-groups depositing their feedback and the expansion of the LB once it has sufficiently slowed down, especially in those regions approaching denser regions of the local ISM.

2. Galactic rotation and radial gradients are neglected possibly affecting the geometry of the stars and gas of the ISM into which the simulated LB expands. At the suggested age of ∼ 14 Myr these effects will likely already be relevant as shown in Paper I.

3. A slightly higher average density and in Schulreich et al. (2023) a uniform midplane density field likely affect the typical size of the simulated LB in these simulations. The absence of a medium, structured by turbulence and self-gravity, leads to the absence of a volume-filling warm low-density phase and low-density channels through which SNRs can expand more rapidly than in a uniform medium (Ohlin et al. 2019).

4. There are significant uncertainties in the ages and explosion times of the progenitor clusters of the LB (Hunt & Reffert 2023; Ratzenböck et al. 2023; Swiggum et al. 2024). By drawing only a single sample these studies neglect these studies likely underestimate the uncertainty arising from this simplification.

In SISSI these simplifications are avoided by placing the SNe in the self-consistently generated ISM of an entire Milky-Way-like galaxy, at the cost of the anyway quite uncertain spatial and temporal sequence of SN events in the local ISM.

Appendix B Analysis of the LB using 3D dust maps

We recovered the properties of the LB by applying a similar analysis as O+24, based on the 3D dust maps of E+24, to obtain estimates for the LB’s geometry, through the means of the shape tensor, defined in Paper I (See also Appendix C.3 for the definition). Some steps of our analysis required the original data products of E+24, so we opted for an independent analysis, instead of directly using the data products of O+24.

We followed the same basic steps of smoothing, peak finding and finally an analysis of the mass distribution enclosed by the peaks, used by O+24, with only minor differences. In contrast to O+24 we directly worked with the logarithmically spaced grid of the 12 sample dust maps rather than sampling the 3D dust maps on a linearly spaced grid of the mean of the sample maps. Moreover, we weighted the smoothing kernel to explicitly account for the non-uniform volume-elements dV ∝ r² dr (Appendix B.1) and we linearly extrapolated the density field beyond the grid boundaries to better capture the profile-shape in their vicinity (Appendix B.2). In Appendix B.3 we show that our results largely agree with those of O+24, demonstrating the robustness of the method. The following steps of peak identification and conversion between differential extinction and hydrogen number density are identical to O+24, who estimate a systematic uncertainty of ∼ 10% for the conversion factor, which we propagate into our mass and momentum estimates.

To obtain the average density of the LB, we divided the total mass, corresponding to the volume integral of the density from r = 0 to the outer edge of the shell, by the volume. The outer edge is defined as the first point past the peak with half the peak’s prominence. We accounted for the mass contained within the central ∼ 69 pc, where the differential extinction is missing from the dust maps, by assuming a constant hydrogen number density along each line-of-sight within this central volume, proportional to the integrated extinction at ∼ 69 pc. The mass in this region accounts for ∼ 7 × 10³ M_⊙, providing a conservative upper limit on the systematic uncertainty due to this assumption.

We estimated the momentum as a function of age by assuming a rapidly cooling wind expansion, R ∝ t^1/2 (Oku et al. 2022; Lancaster et al. 2024), and homologous expansion v ∝ r, giving $$p_{\mathrm{LB}}\left(t_{\mathrm{age}}\right)={\int }_{\mathrm{LB}}v\mathrm{d}M\approx \frac{1}{2t_{\mathrm{age}}}{\int \int }_0^{R_{\mathrm{out}}\left(\Omega \right)}\mu n_{\mathrm{H}}(r,\Omega )r^3\mathrm{d}r\mathrm{d}\Omega ,$$(B.1) where radial integrals sum over radial bins of the unsmoothed profile within the outer shell radius, angular integrals sum over all Healpix lines-of-sight with constant d Ω ≈ 1.6 × 10⁻⁵ rad², and we used v ≈ r/(2 t_age), where t_age is the age of the LB. As shown in Appendix C.1, Eq. B.1 overestimates the momentum by a factor of ∼ 2, which we accordingly corrected for. The resulting momentum as a function of time for fixed size R is consistent with SB expansion under a constant momentum injection rate¹ ṗ_SN, giving p_LB ∝ R⁴/t ∝ ṗ_SN t.

The velocity estimate for momentum-driven expansion provides a conservative lower limit for the momentum. In the case of purely energy-driven expansion, the pre-factor for the velocity would be 3/5, and lead to a 20% higher momentum estimate. The LB might operate in between these regimes.

The assumption of homologous expansion leads to further systematic uncertainties, which we estimate by considering the momentum in the case that all of the gas along each line-of-sight is expanding at the same speed. In this approximation, the momentum is ∼ 12% higher than our estimate from homologous expansion. Conversely a more concentrated velocity distribution within the LB will lead to a slightly lower momentum. We thus estimate the uncertainty due to the assumed velocity structure to be on the order of ∼ 10%, though this is likely highly correlated with the factor-of-two correction factor that we employed.

O+24 have identified an open “chimney”-like topology towards the nothern galactic pole, suggesting that it might be part of a larger blow-out structure that is not covered by the 3D dust maps. While our shell-identification method recovers a well-defined shell position for each line-of-sight, the possibility of such a chimney is not entirely ruled out, since the 3D dust maps are increasingly uncertain towards low-density regions such as the galactic poles. Thus in principle, our derived estimates for momentum, mass and size are to be understood as lower limits.

We summarize the thus obtained properties of the LB in Tab. B.1. For comparison, we also show the values obtained by Z+22 and O+24 where available. Overall our results are similar to the more recent investigation of O+24, however we also calculated the momentum, which we can use to estimate the age of the LB. Uncertainties are dominated by systematic uncertainties at the 60 percent level.

Appendix B.1 Modified smoothing Kernel

O+24 used a Gaussian smoothing-kernel W to obtain smooth density profiles before applying their peak-finding algorithm. This procedure is not manifestly mass-conserving and might lead to numerical artifacts that can bias the results.

In order to account for this potential bias, we used an explicitly mass-conserving approach, by introducing weights that account for the non-uniformity of the volume elements dV = r² dr dΩ. In particular, the condition that the mass is identical between the unsmoothed and smoothed profiles can be written as $$M={\sum }_i\mathrm{d}{V}_i{\rho }_i={\sum }_i\mathrm{d}{V}_i{\overline{\rho }}_i,$$(B.2) where the smoothed density $${\overline{\rho }}_i={\sum }_j{W}_{ij}{\rho }_j.$$(B.3) Condition B.2, leads to the normalization condition for the smoothing matrix W $$\mathrm{d}{V}_j={\sum }_i\mathrm{d}{V}_i{W}_{ij}.$$(B.4) We find that with these modifications a smoothing length of σ_smth=9 pc, slightly above the value used by O+24 yields robust results. However, we find that unless one uses the smoothed density profile for integrals – which like O+24, we do not – the differences due to the mass-conserving Gaussian kernel are only minor, reinforcing the robustness of both methods.

Fig. B.1 Smoothed density profiles using different boundary treatments for two characteristic lines-of-sight. The boundary is indicated by the dashed vertical line. The original data is shown as a gray line, and the extrapolations beyond the boundary are shown at decreased opacity. The top panel shows a line-of-sight with a peak near the inner boundary, while the bottom panel shows a line of sight with a steep decline near the boundary. With reflective boundary conditions, the peak is missed entirely, while simple normalization adds a spurious peak even in the case of the steep decline.

Appendix B.2 The role of the boundary treatment

The choice of the boundary conditions can qualitatively affect the results, by introducing spurious peaks or erasing physical peaks near the boundary.

The simplest way to deal with the boundaries is not to deal with them at all. In particular, consider the smoothing matrix W and assume that it is properly normalized: $${\sum }_i{w}_i{W}_{ij}=w_j,$$(B.5) where w_i are some weights (e.g., volume, a constant, etc.), then weights of the points near the boundary are systematically boosted. This results in smoothed data that will systematically lie below (above) the original data values, if they exhibit a negative (positive) gradient near the boundary. On the other hand, if proper normalization is omitted, points near the boundaries lose weight and the smoothed profile is always systematically lower than the data.

A similar way to deal with boundaries are so-called reflective boundary conditions, where the values of the data as well as their positions are mirrored beyond the boundary, namely, x_−i = x₀ − (x_i−x₀) and y_−i = y_i. This method increases the weight of points near the boundary, but in a more controlled way.

Another common method is to set the data outside of the boundary region to a constant value, for instance the value of the data at the boundary. This leads to a slight over- (under-) estimation near the boundary, in the case of a positive (negative) slope in the original data.

Finally, linear extrapolation of the data beyond the boundary promises to capture steep gradients near the boundaries well. In particular, here we consider a linear extrapolation centered around the boundary point: $$y_{-i}=y_0+\frac{y_i-y_0}{x_i-x_0}(x_{-i}-x_0).$$(B.6)

Since all of these methods are linear, they can be directly incorporated into the smoothing kernel matrix W, which can be precomputed.

In Fig. B.1 we show how these various methods compare for two directions, exhibiting some of the features where differences between the methods are likely to arise, namely, a peak near the boundary, and a steep negative gradient.

In the case of the peak near the boundary it is critical, that the smoother does not remove the peak. However, we find that this is the case for reflective boundary conditions. In the other cases with extrapolation, a shallow peak remains, indicating that in certain cases the peak might still disappear. Only in the case without extrapolation, we recover a strong peak. On the other hand, in the case of a steep negative gradient, it is important that no spurious peaks emerge, a condition that is not satisfied for the smoother without extrapolation alone.

In summary, reflective boundary conditions as well as no extrapolation are not well suited for the task at hand, while linear extrapolation and constant boundaries appear to capture the most important features. Due to the slightly better performance at capturing peaks near the boundary, we opted for linear extrapolation.

Appendix B.3 Comparison to previous work

In order to ensure that our results are robust, we tried to reproduce the shell properties reported by O+24 using the method described in their paper as well as the method described in Sec. B, both applied to the sample mean 3D dust map.

In particular, we sampled the dust map between 69 and 1244 pc at uniform 1 pc intervals and smoothed it with a Gaussian kernel with a smoothing length of σ_smth = 7 pc with reflective boundary conditions. Peak finding and the location of the inner and outer peak boundaries are done identically between the two methods.

In Tab. B.2 we compare the thus obtained results. For the distances to the peak (boundaries) and the width of the shell we show the median, including ± 2 σ uncertainties, corresponding to the 2.28^th and 97.72^th percentiles, while for the shell mass, namely, the mass between the inner and outer edge of the shell, we simply show the value obtained from the analysis of the sample mean 3D dust map. While we could not reproduce the results of O+24 exactly, our results are reasonably close and small differences are likely due to minor details omitted in their description of their analysis.

Curiously, in the reconstruction of the shell by O+24 there are small holes in a few directions where their method cannot identify any prominent peaks, while our method fails to reproduce them. We inspected the profiles and peak identification for a handful of these directions and do indeed see prominent peaks, indicating that these holes might have arisen from the initial sampling step. Due to their small number, they do not meaningfully affect our conclusion about the average geometry.

Appendix C Evaluation of analysis methods

Appendix C.1 Momentum estimate

In order to assess the accuracy of the momentum estimate in Eq. B.1, we computed such geometrical momentum estimates for the SISSI sample of simulated SNRs at different points in time and compared them to the true momentum obtained from the velocity field. The result of this comparison is shown in Fig. C.1, which shows that the geometrical momentum estimate is slightly biased towards greater values, by a factor of ∼ 2, indicative of the (expected) departure from homologous expansion.

The choice of the reference central point, as well as other geometric effects, for instance due to the definition of the SNRs volume may affect this estimate and introduce additional systematic uncertainties reflected by the representative error bar shown in Fig. C.1. The error bar, corresponding to a range of extreme scenarios provides a conservative upper limit for the systematic uncertainty induces by these effects of ∼ 0.24 dex.

Appendix C.2 Number of SN explosions

We assessed the accuracy of the SN count estimate in Eq. 2, using our sample of simulated SBs, by comparing the thus estimated SN count with the actual number of injected SNe. The results are shown in Fig. C.2 which suggests that for moderate densities comparable to that of the LB, the estimate is fairly accurate, while for higher and lower densities it increasingly overestimates the number of SNe as time passes.

This overestimate might be driven by mergers of the simulated, controlled SBs with preexisting SBs in their environment, for which we have no handle to accurately count the number of SNe that they would contribute. At low densities, SBs grow much larger and are thus more likely to merge with other large SBs, while at high densities, the star formation rate is elevated, leading to a more active environment, with a high density of SNRs.

Fig. C.1 Comparison of the geometric momentum estimate Eq. B.1 to the true momentum for the simulated sample of SBs from the SISSI simulation at different points in time, as indicated by the color scheme. A representative error bar corresponding to uncertainties due to different ways of defining the center (geometric center vs. center of mass) of the SBs and their volume (threshold value for passive-scalar tracer-variable) is shown in the upper-left corner. The geometrical momentum estimate is slightly biased towards larger values, by a factor of ∼ 2.

Appendix C.3 Shape Tensor

In Paper I the shape tensor is defined as $$S_{ij}=V_{\mathrm{SNR}}^{-1}{\int }_{\mathrm{SNR}}({\parallel \mathbf{x}\parallel }^2{\delta }_{ij}-x_i{x}_j){\mathrm{d}}^3\mathbf{x}$$(C.1) By assuming an approximately ellipsoidal shape, the three ellipsoidal radii, are defined by $$r_i=\sqrt{2.5(tr\left(S\right)-2S_i)},$$(C.2) where S_i are the eigenvalues of S_ij and tr(S) is the trace. The smallest, intermediate and largest eigenvalues correspond the minor a, semi-major b and major c axis, respectively. The effective size of an SNR is the geometric mean of the three eigenvalues $$r_{\mathrm{eff}}={\left(abc\right)}^{1/3}.$$(C.3)

To determine the alignment of the LB within the Galaxy, we measured the pitch angle α and polar direction cos (θ) for both the major and minor axes. The pitch angle is defined relative to the direction of Galactic rotation, with α=90^° and α=−90^° corresponding to the Galactic center and anti-center, respectively. The magnitude of the polar direction is 0 (1) for directions parallel (perpendicular) to the Galactic plane.

Fig. C.2 Time-evolution of the SN count estimate Eq. 2 for the simulated sample of SBs in different density ranges. The true SN count is shown as a black dashed line. At extremely high and low ambient densities the estimate appears to overestimate the number of SNe, while for moderate densities, around the density of the LB, the estimate is quite accurate, within about ≲ 50%.

Appendix D Models for radiative blastwaves in uniform media

We distinguish between blastwave models driven by a single explosion at t = 0 and blastwaves that are driven by continuous energy- and momentum-injection.

The dynamics of a radiative blastwave driven by a single explosion are determined the conservation of the radial momentum that was acquired before the onset of radiative cooling. This phase is also known as the momentum-conserving snowplow phase (e.g., Oku et al. 2022). The shock radius of the momentum-conserving snowplow as a function of expansion time t is given by (e.g., Oku et al. 2022) $$R_{\mathrm{MCS}}={\left(\frac{3N_{\mathrm{SN}}{\widehat{p}}_{\mathrm{SN}}t}{\pi {\rho }_{\mathrm{ISM}}}\right)}^{1/4},$$(D.1) where N_SN is the number of SNe exploding at t = 0, p̂_SN is the momentum injected per SN, given by Eq. 1 and ρ_ISM is the ambient density.

The dynamics of continuously driven blastwaves are more uncertain and seem to depend on the detailed balance of cooling and energy injection (El-Badry et al. 2019; Oku et al. 2022; Lancaster et al. 2024). If the cooling is relatively gentle and a large fraction of the injected energy can reach the radiative shell, the dynamics of the blastwave match those of an energy-driven wind with a scaled down energy injection rate (El-Badry et al. 2019) $$R_{\mathrm{EDW}}={\xi }_{\mathrm{EDW}}{\left(\frac{(1-\theta )E_{\mathrm{SN}}{t}^3}{\Delta t_{\mathrm{SN}}{\rho }_{\mathrm{ISM}}}\right)}^{1/5},$$(D.2) where ξ_W ∼ 0.88 is determined by considering the internal structure of the blastwave, E_SN=10⁵¹ erg is the total injected per SN, (1−θ) ∼ 0.3 (Zucker et al. 2022) is a reduction factor, accounting for the effect of radiative cooling and Δt_SN is the average time between SN explosions. On the other hand, if cooling is sufficiently rapid each SN contributes the same momentum p̂_SN and the dynamics are determined by the conservation of the momentum, which is imparted onto the radiative shell at a constant rate (Oku et al. 2022; Lancaster et al. 2024) $$R_{\mathrm{MDW}}={\left(\frac{3{\widehat{p}}_{\mathrm{SN}}{t}^2}{2\pi \Delta t_{\mathrm{SN}}{\rho }_{\mathrm{ISM}}}\right)}^{1/4},$$(D.3) which may be compared to Eq. D.1 upon substituting N_SN(t) = t/Δt_SN. This model is often referred to as momentum-driven or rapidly cooling wind (Lancaster et al. 2024).

Fig. E.1 Crossing of the expanding LB’s shell and the Solar System. Gray, red, blue and purple lines are the same as in Fig. 1. Black dashed lines correspond to the distance of the Solar System from the center of the LB (Eq. E.1) for different ages of the LB. The Solar System passes earlier through the LB’s shell for more powerful explosions, but even for a single SN explosion it would have crossed the shell within the first ∼ 1 Myr of its expansion.

Appendix E Passage of the Solar System

Fossil records of sedimentary ⁶⁰Fe and ²⁴⁴Pu on earth suggest an enrichment with SN ejecta over the past 4 Myr (Wallner et al. 2021). If this enrichment was due to the SNe powering the LB, the Solar System must have entered the LB shortly before the enrichment started (Breitschwerdt et al. 2016; Ertel et al. 2023; Schulreich et al. 2023). Making use of this insight, Z+22 estimated that the Solar System has entered the LB ∼ 5 Myr ago.

For a given expansion model for the LB’s shell, we estimated when the Solar System entered the LB, by comparing the traceback distance of the Solar System from the center of the LB to the radius of the expanding shell. We adopted the same magnitude of v_⊙ ∼ 20 km s⁻¹ for the speed of the Solar System as Z+22 and assumed that it moves radially², with a present-day distance from the center of d_⊙(t_age)=0.

Under these assumptions, the distance of the Solar System from the center of the LB is $$d_{\odot }\sim v_{\odot }(t_{\text{age\hspace{0.17em}}}-t).$$(E.1) In Fig. E.1 we show the passage of the Solar System through the LB’s shell. For values of t_age ≳ 4 Myr, the Solar System would have crossed the LB’s shell during the first ≲1 Myr of its expansion. Incorporation of sediments into earth’s crust over the past ∼ 4 Myr could have commenced shortly after the passage (Ertel et al. 2023; Ertel & Fields 2024) and is therefore in agreement with the value t_age ≳ 4 Myr obtained above.

Appendix F Recent star formation and SNe in the solar neighborhood

S+24 have used the star cluster catalogue of (Hunt & Reffert 2023) to group the clusters in the solar neighborhood into three distinct families of star clusters that likely share a common origin. One of these cluster families, named αPer, after its prominent member Alpha Persei can be both spatially and temporally associated with the LB.

Previous analyses of the origin of the LB usually focused on the role of Sco-Cen in driving its expansion (Maíz-Apellániz 2001; Breitschwerdt et al. 2016; Zucker et al. 2022). While the findings of S+24 agree with these studies if only the SNe in Sco-Cen are considered, the potential contribution from other members of αPer remains unclear. S+24 report the SFH in αPer which suggests the existence of a peak in the star formation activity associated with young clusters in Sco-Cen ∼ 10 Myr ago, which might be linked to a more recent burst in SN activity than previous studies considered. Unfortunately, the low time resolution of the reported SFH does not allow to draw any further conclusions.

In order to investigate the role of this recent increase in star formation, we used the data products of S+24 to compute the recent SFH at higher temporal resolution and use it to estimate the SN-rate history for αPer. To this end, we used the 16%, 50% and 84% quantiles of the cluster ages, masses and the number of SNe that exploded since their formation for all clusters associated with αPer.

The SFH is the mass-weighted sum of the age distributions of each cluster 𝒫_i(t), convolved with a window function W(t; Δt): $$SFH\left(t\right)={\int }_{-\infty }^{\infty }W(t-t_{\text{form}};\Delta t){\sum }_{i\in \text{family}}{m}_i{\mathcal{P}}_i\left(t_{\text{form}}\right)\mathrm{d}{t}_{\text{form}}.$$(F.1) Here m_i is the mass of each cluster, which we sample from its posterior distribution.

We can reproduce the SFH reported by S+24 if we use a top-hat filter with Δt=12.5 Myr as the window function, use the median cluster mass for m_i and use delta-peaks, peaked at the median age for the age distributions of the clusters 𝒫_i(t) = δ(t − median (t_form,i)). However, this treatment does not necessarily capture the sizeable uncertainties in the cluster ages and masses, encoded in their 16% and 84% quantiles.

To account for the uncertainty in the age and simultaneously push the time resolution of the SFH to the limit, we reconstructed plausible age distributions (Appendix F.1) from the available data and set W=δ(t−t_form), namely, we took the SFH as the mass-weighted sum of the age distributions of the clusters. We accounted for the uncertainty in the mass by sampling from the posterior distribution, leaving us with a number of sample SFHs.

We estimated the SN-rate history, by assuming that the SN rate of each cluster remains at an approximately constant value between the onset of SN explosions t_delay after the clusters formation and the life-time t_active of the least massive stars undergoing type-II SN explosions (Leitherer et al. 1999; Kim et al. 2017) and is zero otherwise: $${\mathcal{R}}_{\mathrm{SN},i}\left(t\right)=\frac{N_{\mathrm{SN},i}}{min(t_{\text{form},i},t_{\text{active}})-t_{\text{delay}}}{\mathbf{1}}_{(t_{\text{delay}},t_{\text{active}})}(t_{\text{form},i}-t),$$(F.2) where 1_X(x) is the indicator function, which evaluates to unity if x in X and zero otherwise and N_SN,i is the number of SNe that exploded since the formation of the cluster. We sampled both N_SN,i and t_form,i from their respective posterior distributions.

Fig. F.1 Reconstructed posterior age-pdf for the cluster ADS_16795 for various values of β and γ. Dotted, blue lines correspond to the 16%, 50% and 84% quantiles. For smaller values of β the mode of the distribution is closest to the 16% quantile, while for β ∼ 1 it is closer to the median. For β ≥ 1 the pdf approaches a constant non-zero value for x → 0. While the choice of γ has little influence on the shape of the pdf, for γ ≳ 1.1 the pdf is not longer positive-definite. The pdf corresponding to our fiducial choice of parameters is depicted as a solid, black line.

Appendix F.1 Probability Distributions

Given only a number of quantiles x₀ ≤ … ≤ x_i ≤ … ≤ x_n corresponding to percentiles 0 = p₀<…<p_i<…<p_n = 1 and a prior probability density function (pdf) q(x), the most plausible posterior pdf is a locally reweighted version of the prior, where the weights are chosen such that the cumulative probabilities match those of the quantile constraints (Jaynes 1957): $$\mathcal{P}(x\mid q,{\left(x_i\right)}_{i=0,...,n},{\left(p_i\right)}_{i=0,...,n})=q\left(x\right){\sum }_{i=1}^n{w}_i{\mathbf{1}}_{(x_{i-1},x_i)}\left(x\right),$$(F.3) with weights $$w_i=\frac{p_i-p_{i-1}}{\mathbb{P}(x_{i-1}<x<x_i\mid q)}.$$(F.4) With most smooth priors, the posterior distribution will exhibit large jumps, a feature that might be undesirable for our reconstruction of the SFH.

Generating smooth pdfs from quantile constraints is a common data-science application to which a number of creative approaches exist. One particularly simple and flexible method that makes inverse-sampling trivial are so-called quantile-parametrized distributions (Keelin & Powley 2011).

The method uses a parameterization of the quantile function that is linear in the quantile data constraints, using a number of arbitrary basis functions. The corresponding pdf can then be simply derived using the inverse function rule of calculus, since the quantile function is the inverse of the cumulative distribution function. However, care needs to be taken to ensure that the resulting pdf is positive definite.

We adopted this method to reconstruct the mass- and age posterior-pdfs using the following parameterization $$x\left(p\right)=\frac{p^{\beta }}{{(1-p)}^{\gamma }}(a_1+a_{2}p+a_{3}sin\left(\pi p\right)),$$(F.5) where β>0, γ>0, p in [0, 1] is the cumulative probability, and a_i are the parameters determined by the quantile constraints x(p_i)=x_i. For β<1, the pdf tends towards zero for y → 0, while for β>1 it tends toward a constant value >0. We adopted γ=1 and β=0.7 for which we confirmed that all pdfs are positive definite. In Fig. F.1 we show how the pdf depends on the values of β and γ for one of the clusters in αPer.

Since the number of SNe is limited to positive integers we did not need to worry about smoothness. However, there are other complications. The data provided by S+24 correspond to the quantiles of the marginal distribution. Yet, it is clear that the number of SNe has to correlated with the age of the cluster, at the very least, because it can only be non-zero for clusters older than t_delay.

To account for this correlation, we set the probability of N_SN>0 SNe to zero for t_form<t_delay, but otherwise kept the pdf for N_SN>0 unchanged with respect to the marginal distribution. The marginal and joint probabilities of N_SN=0 SNe are then linked by the condition $$\mathcal{P}(N_{\mathrm{SN}}=0)=p_0+(1-p_0)\mathcal{P}(N_{\mathrm{SN}}=0\mid t_{\text{form}}>t_{\text{delay}}),$$(F.6) where p₀=ℙ(t_form<t_delay).

To define a sensible marginal pdf for the number of SNe, we assumed a constant prior for N_SN<N_SN, 84%, an exponential tail for N_SN>N_SN, 84% and we set $$\mathcal{P}(N_{\mathrm{SN}}=0)=max(p_0,{max}_{\{i\mid Q_i=0\}}{p}_i,{min}_{\{i\mid Q_i>0\}}\frac{p_i-p_{i-1}}{1+Q_i}),$$(F.7) where Q_i are the quantiles of the marginal distribution. The joint pdf of t_form and N_SN is then fully specified by applying Eqs. F.3 and F.6.

Appendix F.2 SFH and SN-rate in αPer

For each cluster in αPer we drew a million samples of the cluster mass, age and number of SNe from their respective pdfs and used them to derive a million realizations of the SFH (Eq. F.1) and the SN-rates (Eq. F.2). For each realization, we summed up the SN-rates of the clusters to obtain the SN-rate history of αPer. We also derived the average time between SNe, by taking the reciprocal.

In Fig. F.2 we show the 16%, 50% and 84% quantile histories, obtained by drawing computing the quantiles at each point in time. We mark the onset of the expansion of the LB according to our model and that of Z+22 as well as the SN rates, required by the respective models.

We find that there is a broad peak in the SFH ∼ 4−10 Myr ago, which according to S+24 is associated with the formation of Upper Scorpius, Corona Australis, ρ Ophiuchus, which are all part of Sco-Cen, and various other clusters in Taurus. While Ratzenböck et al. (2023) report a dominant peak in the SFH of Sco-Cen ∼ 15 Myr ago, the SFH in Fig. F.2 associated with the LB is quite low ≲ 100 M_⊙ Myr⁻¹ before the onset of its peak ∼ 10 Myr. These differences can likely be explained due to the different normalization (mass vs. star count), contributions from other clusters in our sample, and prominently the more accurate age estimates used by (Ratzenböck et al. 2023), further highlighting the need for accurate stellar ages. Curiously, the range of ages for the LB found by our analysis above coincides with the downturn in the peak of star formation activity, which is consistent with a scenario where the birth of the LB, powered by the delayed feedback from first stars formed in the peak, is quenching any further star formation in the region, which is in stark contrast to the triggered star formation scenario proposed by Z+22. This difference in interpretation results from the fact that the formation of the LB in the model of Z+22 starts before the peak of the star formation activity while in our scenario it starts after.

Fig. F.2 Top panel: Star formation History of αPer over the past 40 Myr. The SFR is peaked around 8 Myr ago. The peak is associated with recent star formation in Sco-Cen S+24. However, the dominant peak corresponds to the formation of the ∼ 1000 M⊙ cluster Theia 38, at a distance of ∼ 500 pc from the Solar System. Middle and bottom panel: SN rate and average time between SNe in αPer over the past 40 Myr for different values of tdelay. The SN rate continuously grows (the time between SNe shortens) as the stellar mass of αPer builds up, with a sharp increase in SN activity in the last few Myr, associated with the peak in the SFR, 8 Myr ago. The steepness of the increase slightly depends on the choice of tdelay, with a more gradual increase for tdelay=0 and a steeper increase for 0<tdelay<8 Myr. The SFH, SN Rate and average time between SNe of the clusters spatially associated with the LB are shown in purple, where we assume a fiducial value of tdelay=3 Myr. Also shown as shaded regions are the age of the LB as well as the average time between SNe found by Z+22 (orange) and this work (blue).

Fig. F.3 Trajectories of the member clusters of αPer in the local standard of rest, within the vicinity of the LB, from left to right, over the past 5, 5–10, 10–15 and 15–20 Myr. Trajectories are plotted as lines of increasing width, where increasingly wider lines correspond to later times. Trajectories are colored based on the clusters’ instantaneous SN Rate and the opacity is set to the instantaneous probability that the cluster is between tdelay and tactive after its formation and is therefore potentially contributing SNe at a given time. While it is unlikely that a particular cluster colored black has contributed a SN, due to the large stochasticity in low mass clusters, the numerous black clusters might still collectively contribute a significant number of SNe. The projected gas density as well as a projection of the LB’s shell at the present time are also shown to provide context. We also show the number of SNe contributed within each time frame by the clusters, that are co-spatial with the current extent of the LB, using the purple dashed line to decide whether a cluster belongs to the LB or not. With the exception of about 10 clusters in the top right (a few outside the frame) the majority of the clusters has been co-spatial with the current extent of the LB (in projection) for the past few Myr and therefore could have contributed SNe. The number of SNe contributed by nearby clusters in the past 5 Myr is consistent with our estimate in Tab. 1.

The SFH can be compared with the SN rate, which has been steadily increasing from ∼ 1 SN/Myr 20 Myr ago to ≲ 10 SNe/Myr at the present day. The SN rate grows increasingly larger than the value of $\Delta t_{\mathrm{SN}}^{-1}\sim 1{\mathrm{Myr}}^{-1}$ required by the model of Z+22, while it remains well within the range of values required by our analysis above.

To verify, whether this SN activity is actually spatially associated with the LB, in Fig. F.3 we show the trajectories of the clusters in αPer in the vicinity of the LB, colored by their instantaneous SN Rate, over the past 5, 5−10, 10−15, and 15−20 Myr. Trajectories were calculated using the galactic dynamics package galpy (Bovy 2015), following the parameter choices outlined in S+24. We find that the majority of the clusters has been co-spatial with the current extent of the LB over the past few Myr, with just over 10 of the 66 cluster members being too distant to have contributed.

In order to separate out the contributions from these distant clusters, we performed a simple spatial cut, indicated by the purple-dashed line in Fig. F.3. The SN Rate of the remaining clusters is shown as a purple line in Fig. F.2, which shows that the SN Rate of the clusters co-spatial with the LB is still consistent with our findings, highlighted in Tab. 1.

In Fig. F.3 we also show the number of SNe contributed by the clusters associated with the LB during each time frame. Even though there were on average 9 SNe in the time window from 10−20 Myr ago, the fact that most of the contributing clusters have low SN rates suggests, that these SNe merged with the ISM in isolation before they could combine to form a coherent SB. On the other hand, the more frequent SNe in the past 10 Myr might be more clustered, especially towards the Sco-Cen region, where 4–5 exceptionally active clusters reside. These clustered SNe likely overlapped, forming the coherent LB that suppressed further star formation in the local ISM and reached a large enough size to be coherently powered by most subsequent SNe in the solar neighborhood.

All Tables

All Figures

Fig. 1 Time evolution of effective size of the simulated sample of SNRs of the SISSI sample between 2 and 10 Myr, for SNRs with ambient densities within 0.3 dex of the LB (Table B.1). Gray, red and blue lines correspond to different explosion models. Solid lines correspond to the median size of the simulated bubbles, with shaded areas corresponding to the range between the 30th and 70th percentiles. Dotted lines and the hatched blue contour correspond to theoretical models based on radiative blastwaves in uniform media (Appendix D). For model N10, the size of the LB corresponds to an age of ∼ 4.5 Myr, while for N1x10 it corresponds to ∼ 6−7 Myr.

In the text

Fig. 2 Effective size as a function of age or number of SNe, respectively, extrapolated from the model N 1 x 10 and the mean of the momentum constraint Eq. (2) for various values of α; namely, Rα=RN 1 × 10(t=tage) ×[ṄSN, LB(tage)/1 Myr−1]α, where ṄSN, LB= NSN, LB/tage. For α ≲ 0.225, the extrapolated sizes match the LB for ages in the range 3.5−5.5 Myr, while for α ≳ 0.225 the two constraints on the size and momentum cannot be satisfied simultaneously.

In the text

Fig. B.1 Smoothed density profiles using different boundary treatments for two characteristic lines-of-sight. The boundary is indicated by the dashed vertical line. The original data is shown as a gray line, and the extrapolations beyond the boundary are shown at decreased opacity. The top panel shows a line-of-sight with a peak near the inner boundary, while the bottom panel shows a line of sight with a steep decline near the boundary. With reflective boundary conditions, the peak is missed entirely, while simple normalization adds a spurious peak even in the case of the steep decline.

In the text

Fig. C.1 Comparison of the geometric momentum estimate Eq. B.1 to the true momentum for the simulated sample of SBs from the SISSI simulation at different points in time, as indicated by the color scheme. A representative error bar corresponding to uncertainties due to different ways of defining the center (geometric center vs. center of mass) of the SBs and their volume (threshold value for passive-scalar tracer-variable) is shown in the upper-left corner. The geometrical momentum estimate is slightly biased towards larger values, by a factor of ∼ 2.

In the text

	Fig. C.2 Time-evolution of the SN count estimate Eq. 2 for the simulated sample of SBs in different density ranges. The true SN count is shown as a black dashed line. At extremely high and low ambient densities the estimate appears to overestimate the number of SNe, while for moderate densities, around the density of the LB, the estimate is quite accurate, within about ≲ 50%.
In the text

Fig. E.1 Crossing of the expanding LB’s shell and the Solar System. Gray, red, blue and purple lines are the same as in Fig. 1. Black dashed lines correspond to the distance of the Solar System from the center of the LB (Eq. E.1) for different ages of the LB. The Solar System passes earlier through the LB’s shell for more powerful explosions, but even for a single SN explosion it would have crossed the shell within the first ∼ 1 Myr of its expansion.

In the text

Fig. F.1 Reconstructed posterior age-pdf for the cluster ADS_16795 for various values of β and γ. Dotted, blue lines correspond to the 16%, 50% and 84% quantiles. For smaller values of β the mode of the distribution is closest to the 16% quantile, while for β ∼ 1 it is closer to the median. For β ≥ 1 the pdf approaches a constant non-zero value for x → 0. While the choice of γ has little influence on the shape of the pdf, for γ ≳ 1.1 the pdf is not longer positive-definite. The pdf corresponding to our fiducial choice of parameters is depicted as a solid, black line.

In the text

Fig. F.2 Top panel: Star formation History of αPer over the past 40 Myr. The SFR is peaked around 8 Myr ago. The peak is associated with recent star formation in Sco-Cen S+24. However, the dominant peak corresponds to the formation of the ∼ 1000 M⊙ cluster Theia 38, at a distance of ∼ 500 pc from the Solar System. Middle and bottom panel: SN rate and average time between SNe in αPer over the past 40 Myr for different values of tdelay. The SN rate continuously grows (the time between SNe shortens) as the stellar mass of αPer builds up, with a sharp increase in SN activity in the last few Myr, associated with the peak in the SFR, 8 Myr ago. The steepness of the increase slightly depends on the choice of tdelay, with a more gradual increase for tdelay=0 and a steeper increase for 0<tdelay<8 Myr. The SFH, SN Rate and average time between SNe of the clusters spatially associated with the LB are shown in purple, where we assume a fiducial value of tdelay=3 Myr. Also shown as shaded regions are the age of the LB as well as the average time between SNe found by Z+22 (orange) and this work (blue).

In the text

Fig. F.3 Trajectories of the member clusters of αPer in the local standard of rest, within the vicinity of the LB, from left to right, over the past 5, 5–10, 10–15 and 15–20 Myr. Trajectories are plotted as lines of increasing width, where increasingly wider lines correspond to later times. Trajectories are colored based on the clusters’ instantaneous SN Rate and the opacity is set to the instantaneous probability that the cluster is between tdelay and tactive after its formation and is therefore potentially contributing SNe at a given time. While it is unlikely that a particular cluster colored black has contributed a SN, due to the large stochasticity in low mass clusters, the numerous black clusters might still collectively contribute a significant number of SNe. The projected gas density as well as a projection of the LB’s shell at the present time are also shown to provide context. We also show the number of SNe contributed within each time frame by the clusters, that are co-spatial with the current extent of the LB, using the purple dashed line to decide whether a cluster belongs to the LB or not. With the exception of about 10 clusters in the top right (a few outside the frame) the majority of the clusters has been co-spatial with the current extent of the LB (in projection) for the past few Myr and therefore could have contributed SNe. The number of SNe contributed by nearby clusters in the past 5 Myr is consistent with our estimate in Tab. 1.

In the text