Supernova Dust Evolution Probed by Deep-sea ⁶⁰Fe Time History

The presence of deep-sea pulses of live ⁶⁰Fe is now compelling, providing strong evidence for recent nearby supernovae. Figure 1 compiles the published ⁶⁰Fe data. The figure shows the detected ⁶⁰Fe/Fe isotope fraction versus time, categorized by both sample type and research group; none of the data have been background corrected. The data show two distinct peaks, with all groups agreeing on a ⁶⁰Fe peak around 2–3 Mya and the Wallner et al. (2016, 2021) data indicating another peak around 6–7 Mya. Note that this peak at 7 Mya only appears clearly in the Wallner et al. (2021) data, whose ⁶⁰Fe machine background has finally been lowered enough to show a distinct peak; the background in earlier efforts obscured the signal. Not included in this figure are the recent ⁶⁰Fe flux measurements by Koll et al. (2019) and Wallner et al. (2020), as they cover collectively only the last 30 kyr and would not be visibly distinguishable from the origin (see Section 2.2). The Fimiani et al. (2016) lunar data are included for completeness only: due to micrometeorite gardening effects on the lunar regolith, the data cannot be time resolved to better than ∼8 Myr. ¹⁰ Because of this time range, the lunar data include the signals due to all supernovae in the last 10 Myr; however, due to the half-life of ⁶⁰Fe, the contribution of ⁶⁰Fe from the pulse at 7 Mya is only around 10%.

Zoom In Zoom Out Reset image size

When comparing the ⁶⁰Fe/Fe measurements, one should bear in mind that the geographical distribution of ⁶⁰Fe may not be uniform, and that the uptake, U, of ⁶⁰Fe varies in different materials. Specifically, it was shown in Fry et al. (2016) that the transport of ⁶⁰Fe through the atmosphere should not be isotropic, but rather would favor middle latitudes ∼60°, with minima at the equator and poles (see also Dhomse et al. 2013), yielding a factor of ∼3 in the global difference. This may explain some of the large difference between the flux values reported by Wallner et al. (2016) taken at ∼38°S versus the those in the data of Ludwig et al. (2016) taken at 3°S. In addition, ocean currents may cause variations with longitude in the deposition of ⁶⁰Fe in FeMn crusts and deep-sea sediments at similar latitudes. Differences in analytical technique and sample processing can also affect the final result. We note finally that the uptake is expected to be 100% for sediments and snow, whereas the uptake in FeMn crusts is subject to considerable uncertainty and debate, and may vary depending on location and local conditions (Bishop & Egli 2011).

In this paper, we focus on the sediment data when estimating the timescale of astrophysical ⁶⁰Fe deposition, in view of the availability of multiple samples with time resolution from deep-sea sediments (Ludwig et al. 2016; Wallner et al. 2016). We note that the Wallner et al. (2021) ⁶⁰Fe data from the FeMn crust also has excellent time resolution compared to earlier studies. However, analysis of the deposition timescale requires strict accounting for geophysical processes that might disturb the signal, and this is more straightforward when the data are from the same sample type. By focusing on the unbinned sediment data, we can control most of the variables between the two data sets and therefore make fair comparisons. It should be noted that this approach by necessity deals with very small-count statistics; we address the issue in Section 3.1. We look forward to additional measurements of the 7 Myr peak, ideally in sediments, to allow for a similar analysis of that event.

Koll et al. (2019) and Wallner et al. (2020) have shown independently that there is recent and ongoing infall of extrasolar ⁶⁰Fe onto Earth. Koll et al. (2019) detected an ⁶⁰Fe signal in Antarctic snow deposited over the last 20 yr and use isotopic ratios to show it is not from meteoritic material and must therefore be extrasolar. They suggest that the signal is due to the solar system passing through the Local Interstellar Cloud (a nearby higher-density region of the Local Bubble). Wallner et al. (2020) also detected a fairly steady ⁶⁰Fe signal in deep-sea sediments over the last 33 kyr, in line with the Koll et al. (2019) detection. However, they do not see the sharp increase in the ⁶⁰Fe signal that would be expected from the solar system entering the Local Interstellar Cloud. There are other possibilities for the persistent ⁶⁰Fe flux, including continued delivery of dust following the most recent pulse ∼3 Mya, or flux from the Earth's motion through the local ISM. Further analysis of ⁶⁰Fe in FeMn crusts and deep-sea sediments focused on the age range from 40 Kya to 1 Mya could add significant insight into the origin of the observed recent infall.

That said, the main purpose of this paper is to characterize the peaks that are clearly evident in the data, and we do not attempt to fit the low-level ⁶⁰Fe flux outside the peaks.

The purpose of this work is to determine the deposition or "raindown" time onto Earth of ⁶⁰Fe during the recent pulse ∼3 Mya, and to interpret what this timescale implies about the propagation of supernova-produced material inside the remnant. In order to perform this analysis, we fit the observed ⁶⁰Fe signal with a number of 3-parameter pulse shapes to see which one best described the data, while assuming the errors are dominated by Poisson counting statistics. These shapes included a Gaussian, sawtooth, reverse sawtooth (for comparison, despite our doubt that this is a physically plausible profile), and a symmetric triangle; expressions appear in Appendix A. We also performed one 4-parameter fit with an asymmetric triangle to see if there is a preferred slant to the data.

The general shape of the ⁶⁰Fe data is of particular interest, because the sharpness and slant of the shape can provide insight into the astrophysics of the deposition of the ⁶⁰Fe and its path within the supernova remnant. Predictions in the literature to date have assumed the ⁶⁰Fe traces the gas phase of the blast. Under this assumption, Fry et al. (2015) showed that if the ⁶⁰Fe is well mixed in a Sedov blast, the signal appears discontinuously with the arrival of the forward shock, and decreases thereafter from this maximum. Chaikin et al. (2022) included the effects of incomplete mixing of ⁶⁰Fe in the remnant gas, and also allowed for effects of the Earth's motion. They too found that the ⁶⁰Fe pulse begins abruptly and is concentrated in a pulse, but found that the signal can linger thereafter. Thus, these models would favor the discontinuous profiles we have considered—the sawtooth form, or a cut exponential. As we will see, there are other possibilities if the ⁶⁰Fe is in dust that is decoupled from the gas, so we have chosen a suite of different fitting functions to allow for a range of possible ⁶⁰Fe flux histories.

In the published versions of their work, both Wallner et al. (2016) and Ludwig et al. (2016) bin their data, which serves to demonstrate the strength and overall peak of the signal. For the purposes of this analysis, we use the original, unbinned data, as it removes extraneous smoothing, and we are specifically interested in fitting the unbinned shape. The unbinned ⁶⁰Fe data can be found in Table S.4 of Wallner et al. (2016) and in Tables A.1 and A.2 of Ludwig (2015). Both Wallner et al. (2016) and Ludwig et al. (2016) assume zero background for their analyses. As mentioned above, recent works by Wallner et al. (2020) and Koll et al. (2019) have found a nonzero ⁶⁰Fe background today. This minor discrepancy is discussed in more detail in Section 2.2. For this work, we use the zero background estimate assumed by the original analyses.

In each sediment, AMS measures individual ⁶⁰Fe atoms in each sediment segment corresponding to a time bin. The numbers of counts measured is small: in all of their sediments, Ludwig et al. (2016) detected a total of 89 atoms of ⁶⁰Fe, while Wallner et al. (2016) measured a total of 288 atoms. The n_i number of ⁶⁰Fe counts in each time bin i is therefore also small, with n_i ≤ 23 and often n_i < 10. As a result, the Poisson errors in the counts dominate the uncertainties in the resulting ⁶⁰Fe flux. We therefore tailor our analysis to identify the dependence on count numbers n_i and to treat these Poisson errors appropriately.

The observed ⁶⁰Fe flux values Φ₆₀ (the deposition into the material) are computed by combining the measured ⁶⁰Fe counts with properties of the sediment, as follows. For each time t_i , a number n_i of ⁶⁰Fe atoms are measured. From this, the observed ⁶⁰Fe/Fe ratio is determined:

where the scaling s_i is due to variations in efficiency and other factors, and unique for each measurement. We infer the scalings from the reported n_i and (⁶⁰Fe/Fe)_i . In the cases of null measurements, we follow the same procedure as both experimental groups, using the Feldman & Cousins (1998) prescription for calculating 69.29% CL limits based on n_i = 0 counts, and a background b = 0. This gives an effective limit , which we use to infer the scaling s_i for these measurements.

We use the following equations to correct the data for decays:

Here and throughout, the subscript "c" indicates that the quantity is decay corrected, t_i is the observed time with uncertainty σ_i , τ is the mean lifetime of ⁶⁰Fe, and is the uncertainty in the decay-corrected ⁶⁰Fe/Fe ratio. To calculate the ⁶⁰Fe flux, we use

In Equation (4), Φ is the flux of ⁶⁰Fe, and n₆₀ is the ⁶⁰Fe number density. The sediment mass density is ρ, is the sedimentation rate, and X_Fe = ρ_Fe/ρ is the mass fraction of iron in the sediment. The factor c_Fe = X_Fe/A_Fe m_u measures the concentration of iron in the sediment in atoms per unit mass, with A_Fe as the mean molecular weight of iron and m_u as the atomic mass unit.

Due to the format of the data provided in the relevant papers, the conversion of the ⁶⁰Fe/Fe ratio into a flux was by necessity different for the two groups. Wallner et al. (2016) calculated the flux in their Table S.4, following the formula in Equation (4), and we use those numbers directly. However, Ludwig et al. (2016) used a protocol in their chemical sample treatment (citrate-bicarbonate-dithionite, hereafter CBD) that effectively separates larger iron-bearing grains from smaller grains, which they argue should arise from fossilized magnetotatic bacteria. Their AMS measurements thus give a (⁶⁰Fe/Fe)_CBD ratio for this material, and thus their scaling c_Fe of iron per mass in the bulk unprocessed sediment includes a factor Y_CBD for the fraction of iron selected by the CBD process. This implicitly assumes that the Fe-bearing material excluded from the CBD protocol does not contain ⁶⁰Fe. For the flux calculation in Equation (4), the unbinned ⁶⁰Fe/Fe ratio and extracted iron are from Tables A.1 and A.2 in Ludwig (2015; the binned versions of which appear in Ludwig et al. 2016), the sediment density is taken from Table 6.3 in Ludwig (2015), and the sedimentation rate is from Figure 1 in Ludwig et al. (2016). ¹¹

Combining Equations (1) and (4), we arrive at the relationship between the flux and ⁶⁰Fe counts. For time t_i , we have

Given the counts n_i and the other observables, Equation (7) allows us to determine the flux scalings φ_i . As we now see, these allow us to compare the observed fluxes to model predictions.

Our goal is to fit the observed flux data with several simple models that capture different qualitative trends one might expect in the ⁶⁰Fe flux versus time, Φ_model. The flux profile versus time is described by a set of parameters θ , so that we have Φ(t; θ ). In designing our fitting procedure, we are guided by the fact that the uncertainties in the ⁶⁰Fe flux are dominated by the Poisson errors in the ⁶⁰Fe counts, which is the case for both the Wallner et al. (2016) and Ludwig et al. (2016) data sets. We have designed our fitting procedure to accommodate this situation, closely following the approach laid out by Cash (1979), originally for determining X-ray fluxes from photon count measurements dominated by Poisson uncertainties.

Our analysis requires that at each time t_i we specify the expected number μ_i of events, based on the model flux. To do this, we evaluate Φ_model(t_i ; θ ) and then infer μ_i ( θ ) = Φ_model(t_i ; θ )/φ_i , using the the scaling φ_i found in Equation (7). For each measured number of counts n_i , the fit function with parameters θ gives an expected value μ_i ( θ ).

We now construct a Poisson-based likelihood for the fit given the data. For time t_i , the likelihood for the fit given the data is just the Poisson probability , where P(n∣μ) = μⁿ e^μ /n ! is the Poisson probability of n counts given a mean μ. The total probability for the fit given all of the data is just the product of the likelihoods at each time:

where the fit parameter dependence is through μ_i ( θ ). The negative of the logarithm of the total likelihood in Equation (8) is

Here the constant sums terms with dependencies that do not depend on the fit parameters θ , which means that it does not affect the relative likelihoods of different parameter choices, so we follow the usual practice and neglect it. The function C thus depends on the data and the fit parameters, and determines the goodness of fit in a way closely analogous to the role of a χ² for data that is continuous rather than discrete.

For a given data set and fitting function Φ_model(t; θ ), Equation (9) will give different values to C for different choices of the parameters θ . The likelihood is maximized at the minimum value for C, which we call , and which we will use to assess goodness of fit. The parameters θ giving are the best-fit values. We use Monte Carlo methods to explore the parameter space, find the best-fit parameters, and characterize their uncertainties.

We note that the physical picture of supernova radioisotope deposition could include an abrupt onset to the flux. This could occur if ⁶⁰Fe is entrained in a blast wave, so that the onset of the ⁶⁰Fe coincides with the forward shock's arrival at the solar system. To allow for this, we consider some fitting functions where the flux has an abrupt onset and/or halt. In these cases, there are times before or after the blast passage, for which there is no signal: μ_i = 0, which means that the expected number of counts must be zero. If the measured number of counts is nonzero for these times, then the Poisson probability is zero (C → −∞), and this set of fit parameters is completely ruled out. In other words, our method automatically rejects the models (regions of parameter space) that predict no counts where some are observed. On the other hand, the converse is not true: if the fit has μ_i > 0, Poisson statistics allow for cases where n_i = 0, albeit with a penalty.

Each of the four sediment cores was fitted separately, as each core is a separate time column, and in many cases, data points from different cores lie on the same time slice, which creates difficulties with Poisson statistics. We were also interested in examining any differences we could find in the pulse arrival time and the peak pulse time, to see if there were differences in the timing calibration between the different sediments. In order to ensure a universal resolution for all the pulse shapes, we enforced the same initial parameters across all of the 3-parameter fits (similar values were used for the 4-parameter fit, although it is not statistically comparable).

To test our methodology, we generated simulated data points, drawn from a predetermined flux history Φ_true(t). We randomly chose sample times t_i , and drew counts n_i from a Poisson distribution appropriate for our flux history. We found that our method generally performed well: the best-fit parameters were close to the parameters of the known input Φ_true. However, we did find that the accuracy and precision of the width and thus timescale parameter depends on the functional form and the time distributions of measurements. In particular, we found that a crucial factor is the number of measured points with zero counts before and after the bulk of the signal. The more leading and trailing null points, the better the width was determined. On the other hand, if there were only one or two null points on either side of the signal, the width was less well constrained, with the true width being at the low end of the range allowed by the fit. This is a manifestation of the physical effect that in noisy data with small mean numbers of counts, one or two bins with zero counts do not strongly constrain the fit; rather, many are needed to exclude the models that span wide timescales. We will see this behavior manifested in the fits below.

To allow for a variety of possible ⁶⁰Fe time histories, we fit the flux data with six possible 3-parameter shapes, chosen for mathematical simplicity and resemblance to physically motivated trends suggested in the literature. Their mathematical expressions are given in Appendix A. Three of these give a signal that has a finite duration: (1) a symmetric triangle with equal duration linear rise and fall around a peak, (2) a sawtooth that begins abruptly at a peak and falls linearly to zero, (3) and a reverse sawtooth that rises linearly from zero to a peak, then drops to zero. We explored three additional profiles that allow comparison to traditional fits and explore a more gradual rise and fall that formally never goes to zero: (1) a Gaussian, (2) a Lorentzian, and (3) a cut exponential that starts at a peak and then drops exponentially.

The data for each sediment core were fit with each specific fit model. The results are summarized in Tables 2 and 3, but for brevity we plot results only for select cases. Figures 2–5 show results for all cores using the sawtooth fit, while Figure 6 shows the results for the Ludwig 848 core using a Gaussian fit. Each fit figure shows the flux versus time data for a single core. The middle part of the figure plots the ⁶⁰Fe flux (in 10⁴ atoms cm⁻² kyr⁻¹) versus time (in Mya), including both detections and nondetections. Overlaid on top of these points is the best-fit curve, with the statistic and the fit's peak time (t_peak), peak flux (Φ_peak), and width (σ_t ) given in the upper left corner. The top three plots in each figure display two-dimensional projections of the three parameters of each fit model, exhibiting the contour confidence levels of ΔC = corresponding to 1σ, 2σ, and 3σ for a three-dimensional Gaussian; the bottom three plots show the marginalized likelihood for each parameter normalized so that the peak is at 1.

Zoom In Zoom Out Reset image size

Table 2. Goodness of Fit: Values for 3-parameter Fits

	Values
	Wallner et al. (2016)		Ludwig et al. (2016)
Model	Core 4521	Core 4953	Core 848	Core 851
Cut Exponential	−311.17	−158.58	63.33	30.32
Gaussian	−312.56	−174.25	62.24	26.55
Lorentz	−312.45	−173.15	65.04	24.62
Sawtooth	−311.02	−158.26	57.73	29.00
Reverse Saw	−311.72	−153.66	80.61	39.01
Triangle	−312.59	−174.34	59.85	25.56

Note. The most negative value for each core gives the best fit, which is shown in boldface. Fits that have nearly identical (for the respective core) are in italics. Best-fit values should be compared between models for the same core (vertical columns), and not between cores (horizontal rows).

Download table as:  ASCIITypeset image

Figure 2 shows the results for a sawtooth fit to the Ludwig core 848, and serves as an exemplar for similar plots of other cores and fits. In the top panels, we see two-dimensional slices of the θ = (t_peak, σ_t , Φ_peak) parameter space. The red dot gives the location of the best fit, where . The surrounding contours correspond to the 68%, 95%, and 99% confidence-level values. We see that the best-fit regions are relatively compact, meaning that the best fit is well determined. The peak flux and peak time values are the best determined, with quite small uncertainties. The width parameter shows a broader distribution and hence larger uncertainty. These trends are reflected in the bottom panel of Figure 2, where we give the one-dimensional marginalized likelihood distributions such as

We see that unlike the well-determined peak time and flux, the sawtooth width distribution shows a rapid rise and a long tail toward long durations.

The middle panel of Figure 2 shows the data for this core, overlaid with the curve for the best-fit parameters. We see that this optimal sawtooth function "turns on" essentially at the earliest nonzero measurement (largest time before present), denoted t_first in Table 1. The sawtooth onset is the peak time t_peak, and so for the best fit, this essentially is set by the first nonzero measurement. The best-fit curve goes to zero soon after the last nonzero point (smallest time before present). This means that the width parameter σ_t is in this case essentially set by the time interval between the first and last nonzero measurements shown in Table 1. We see that the height of the sawtooth at onset is a compromise among the data points so that, unsurprisingly, the measurements with the smallest errors determine the peak height Φ_peak and also influence the slope and hence the width.

We can understand the broad width σ_t distribution for this and other sawtooth fits by considering the interplay between the data and the best-fit curve in the middle panel of Figure 2. The nonzero measurements exact a "cost" in goodness of fit C for models that miss them, with the most extreme case being outright rejection (C → ∞) of models that predict zero flux where counts are nonzero. Conversely, the points with zero counts impose a cost in C for fits that are nonzero, but this penalty is less severe, corresponding to allowing for Poisson fluctuations. This means that sawtooth fits will be highly suppressed if they are narrower than the nonzero count data, but will have some freedom to extend beyond the width of the data, until the available zero-count points suppress fits that are much wider than the data. This is the trend we see in the width distribution.

These insights from Figure 2 elucidate also the trends in sawtooth fits to the other sediment cores. Figure 3 shows sawtooth fits for Ludwig core 851. Here again we see that the peak time and peak flux values are fairly well determined; though the top and bottom panels show that the peak time has a sharp lower limit, but its distribution extends for about 0.4 Myr beyond this. We can understand this feature from the middle panel: there is a lack of data between the earliest measured nonzero point and prior zero-points. This gap is about 0.3 Myr, and the lack of data here means that there is no penalty for fits that have an onset t_peak anywhere in this range. This leads to the width in t_peak. ¹²

Zoom In Zoom Out Reset image size

Even more striking is the width parameter in Figure 3: we see in the bottom panel that the width distribution comes to a peak at 1.84 Myr, but the distribution is highly skewed. The σ_t likelihood cuts off rapidly below the peak, but shows a long tail beyond the peak that extends out to the longest values that were allowed. Again, the data in the middle panel show the reason: after the last nonzero point (earliest time before present), there is only a single point with zero counts. Thus there is little penalty for fits with a width extending far beyond the interval between the nonzero point. Furthermore, the scatter of the nonzero points allows for a wide range of slopes, which also permits large widths. The lesson here is clear and reasonable: to get a strong constraint on the width is it essential to measure multiple points with zero⁶⁰Fe counts before and after the points with nonzero counts. For Ludwig core 848, this is the case, and the width is better constrained than that of Ludwig core 851.

This lesson is underscored when we consider Figures 4 and 5, which show sawtooth fits to Wallner cores 4521 and 4953 respectively. For both cores, we see that effectively we can only set a lower limit to the width parameter. That is, the width likelihood in the bottom center panel begins to rise after some minimum σ_t , and continues to increase up to the highest allowed value. Looking at the data, we see that both cores have no points with zero counts before or after the nonzero counts. Thus, the nonzero count duration sets a lower limit to the width, corresponding to the onset of the likelihood rise, about (1.6, 2.5) Myr for cores (4521, 4593). However, the available data essentially set no upper limit to the width for these cores. We also see that the peak time is poorly constrained, again due to the lack of zero-count data at times earlier than the first nonzero point.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have produced plots in the style of Figures 2–5 for the other fitting functions. For brevity we show here only the Gaussian fit to Ludwig core 848, which appears in Figure 6. The main trends are similar to those we saw for the sawtooth fit to this core in Figure 2, and the peak width and peak flux are quite well determined. For the Gaussian case, we also find that the width is well determined, better than the sawtooth case. For the cores not shown, it is illuminating to compare the trends in the Gaussian fit to those found in the sawtooth fits shown in Figures 3–5. For the Ludwig core 851, the peak time and width are broader than those in core 848, but are still well determined. On the other hand, for Wallner core 4521, the Gaussian fits also give only a lower limit to the width, and the peak time likelihood does have a clear maximum but broad tails on either side. Interestingly, the Wallner core 4953 fits are better determined, with the width and other likelihoods showing clear peaks and wings that go to zero on both sides. We believe this is due to the effect of a few nonzero points with high-precision fluxes, which anchor the Gaussian fit and prevent large excursions away from the best-fit region.

Zoom In Zoom Out Reset image size

Among the other fits we tried, the reverse sawtooth case also deserves mention. Here the fit has a linear rise at early times, ending with an abrupt cutoff at late times; this was chosen to contrast with the more physically motivated sawtooth case. For these fits, we find that the width parameter likelihoods set only lower limits for all cases. This includes the Ludwig cores for which it had been possible to determine the width in the ordinary sawtooth case.

Figure 7 summarizes the best-fit curves for the six different fitting functions. Several trends emerge. Regarding the pulse widths, we see that, for the symmetric functions (Gaussian, Lorentzian, triangle), the best-fit curves all span at least 1 Myr, with Wallner cores spanning >2 Myr. For the asymmetric functions (sawtooth, reverse sawtooth, cut exponential), the best-fit curves are generally even wider. The initial signal arrival is sharply defined in the sawtooth and cut exponential cases, where it exceeds 3 Mya for at least one Wallner and one Ludwig core. In the other cases, the onset is gradual but also begins no later than 3 Mya.

Zoom In Zoom Out Reset image size

Turning to the signal amplitude, we see that in all cases the Wallner peak fluxes are higher than those of Ludwig. These differences were discussed in Section 2, and may point to geophysical differences in ⁶⁰Fe fallout, and could also reflect differences in ⁶⁰Fe extraction techniques. We note that the ⁶⁰Fe/Fe ratios in the two Ludwig sediments are similar, as are the isotope ratios in the two Wallner sediments. The differences between the fluxes in the two Ludwig (and two Wallner) cores are thus largely due to the differences in sedimentation rates. For example, the sedimentation rate for L-848 is a factor ∼3 lower than that of L-851. These differences propagate to our fits summarized in Figure 7, with the Gaussian and Lorenzian profiles showing similar ratios. But the fit shape has an influence, as the triangle fits give lower ratios of peak flux, particularly for the Wallner cores, while the cut exponential and sawtooth fits give larger ratios. ¹³

To provide a basis for comparing quantitatively the different fits for each core, we compiled all of the best-fit values for the six fits shown in Figure 8 and Table 2; for each data set, C measures the negative of the log of the likelihood, in a close analogy to a χ² for continuous data. The minimum value thus identifies the point of maximum likelihood, playing the role of the minimum χ². We note that, for a given core, quantifies goodness of fit, in that the more negative the value, the better the fit. As with χ², variations of around the minimum are not statistically significant, while larger variations indicate a preference for the model with lower . We also note that, while the can be compared between fits for the same core, it is not appropriate to compare the relative values for fits between different cores (for example, the fact that the values for the Ludwig cores are above zero and the values for the Wallner cores are below zero says nothing about the relative goodness of fits to the Wallner and Ludwig cores). In Table 2, the most negative value for each core has been shown in boldface. We show in italics values that are nearly identical for the respective core have been shown in italics, indicating that the respective fit functions for the bold and italic values are all of similar quality.

Zoom In Zoom Out Reset image size

We find that, for Wallner Core 4521, the triangle, Gaussian, and Lorentz fits are all equally good. As seen in Figure 4, the data for Core 4521 are quite irregular, which could account for the multiple favorite fit shapes. Wallner Core 4953 also has a preference for the triangle and Gaussian fits. Meanwhile, the Ludwig Cores 848 and 851, which are more heavily sampled than the Wallner cores, have a clear preference for the sawtooth and Lorentz fits, respectively. The main takeaway from these results is that the ⁶⁰Fe pulse does not have a preferred shape, even across the same core. Although we frequently describe the widths of the six fitting functions as the width timescales, it should be noted that these are not actually good measures for comparing different fit shapes. For example, σ_t for the sawtooth fit is the full width (FW) of the fitting function and therefore is the actual timescale for that curve shape. However, σ_t for the Gaussian is not given directly by the actual start and stop times of the function, since the Gaussian never reaches zero flux. Therefore, in order to compare the timescales for each function (and find a preferred time width for the supernova pulse), we have examined the traditional FW at half maximum time for each fit for each core. Since our primary interest is the maximum width of the function, we have also plotted the FW at 0.1 the maximum (which is closer to the true timescale). Table 3 and Figure 9 show differences in the width determination between the Wallner and Ludwig samples. The values corresponding to the best-fit values from Table 2 are also given in boldface and italics where relevant. As we have noted, these reflect differences in sampling rather than a true discrepancy; hence the robust conclusion is the lower limit this places on the signal duration.

Zoom In Zoom Out Reset image size

Thus the presently available sediment data carry enough uncertainty (particularly in sampling) to prevent an unambiguous measurement of the preferred shape or timescale for all cores. As future measurements are made, these questions will become all the more important, because a demonstrated consistency among the measurements will allow us to probe the underlying astrophysics.

Our key result is that the FW at 0.1 maximum height, for all functions and cores, shows that the width of the deposition timescale is at least 1 Myr. This is significantly longer than the prediction of the traditional Sedov model.

We now discuss how terrestrial and geophysical effects might smear the timescale of the ⁶⁰Fe pulse.

The atmosphere. When the dust, which is traveling at up to 100 km s⁻¹, hits the Earth's atmosphere, it is vaporized. The iron atoms then combine in the upper atmosphere with molecules such as ozone and hydroxyl radicals, and are eventually deposited on land and in the ocean. Fry et al. (2016) provides excellent approximations for the residence times of the iron in the atmosphere and demonstrates that the iron settles out of the atmosphere in less than 10 yr.

Deposition on land. There is no currently known way to detect the ⁶⁰Fe signal on land, unless it is deposited on an ice sheet, e.g., in Antarctica, where ongoing deposition has recently been measured (Koll et al. 2019).

Deposition in the ocean. The residence, i.e., removal, time of iron in the ocean is relatively short, on the order of ∼500 yr at most (Bruland et al. 1994; Boyle 1997; Resing et al. 2015). When the dust settles on the ocean floor, it may be taken up by FeMn crusts or deposited as sediment.

FeMn crusts. Once iron is absorbed by a FeMn crust, it remains "locked in" and unchanged until analysis, and the time signal is accurately preserved. However, due to the slow growth rate, it is very difficult to time resolve FeMn crusts on the order of kyr, and only the recent work of Wallner et al. (2021) has done so. There is also the possibility of crust porosity, which would enable the iron to attach below the surface, thus smearing the signal.

Sediments. It is possible to resolve the time structure of deposits in sediments on the order of kiloyears, but there are a number of effects that can smear the time signature. The two most important are bioturbation (the churning of the upper few layers of sediment by macroscopic organisms) and chemical reducing environments (created by bacteria and leading to the movement of iron within the sediment column). However, all the sediment samples we consider are from the deep sea where bioturbation effects are negligible, and were carefully selected to ensure that a reducing environment did not occur (Fitoussi et al. 2008; Ludwig et al. 2016).

In summary, there are a number of geochemical, geophysical, oceanic, and biological effects that can in principle distort the ⁶⁰Fe timescale. However, the combination of these effects is less than Δt ∼ 10³ yr, which is far shorter than the Δt ∼ 10⁶ yr timescale found in the data. We conclude, therefore, that the measured timescale must be astrophysical in origin. Accordingly, we need a model for the origin and transport to Earth of the dust that can accommodate the signal width and might also be able to predict the line shape, which could be constrained by future data with higher precision.

There is significant interest in analyzing the possibility of multiple supernovae to account for the signal ∼3 Myr ago, e.g., Breitschwerdt et al. (2016) and Schulreich et al. (2018) propose that some 16 to 19 supernovae have exploded in the Local Bubble, and have contributed to the extended ⁶⁰Fe signal. Unfortunately, the data are too noisy to cleanly distinguish multiple superposed supernovae peaks. We have seen that the sediment measurements cannot unambiguously distinguish the relatively simple pulse shapes we have tried, which have shapes as different as possible for a singly peaked structure. Thus, the data in hand cannot exclude a more complex pulse shape that would superpose multiple supernova pulses.

Despite these limitations, the ⁶⁰Fe data carries substantial information bearing on the question of multiple supernovae creating the observed 3 Myr pulse. In particular (a) singly peaked pulse shapes provide adequate descriptions of the measurements, and (b) none of the data show clear groupings of points within the time range of detected points. The available data are thus consistent with a single peak, and do not require multiple events.

Furthermore, an accounting for the ⁶⁰Fe data must explain not only the long signal width seen in the samples but also the discovery of distinct pulses at 3 and 7 Myr, with no apparent signal in between. If there are multiple events in the 3 Myr peak, then there would likely need to be a similar set of events for the 7 Myr peak, but then the gap between the two would need explanation. These considerations will inform models for multiple supernovae, but do not rule them out.

We do not consider it a productive exercise to fit for multiple supernovae given the limitations of the available data, so we restrict our analysis to a single supernova. We discuss models for the pulse width below in Section 4.

Additional time-resolved measurements of the 3 Myr ago signal (such as could be provided by more dense sampling) could help enormously, but would require significant extra efforts beyond those already made to gather the current data set. Another possibility would be to measure the fluxes of additional isotopes across the ⁶⁰Fe signal region. If there were significant variations in the isotope ratios, these could constitute evidence for multiple supernovae with different combinations of nucleosynthesis mechanisms.

In order to examine the ambiguity of the preferred fit shape for the data, we include a 4-parameter "sharktooth" fit in our analysis. The purpose of this fit is to analyze whether the data is better described by an asymmetric shape, and at what level of preference. We chose to work with a 2-width triangle shape, of which the sawtooth, reverse sawtooth, and symmetric triangle are the three-parameter extremes. The function is defined below in Appendix A. The shape was given the full range from almost full sawtooth to almost full reverse sawtooth, with only a minor initial minimum of 0.05 Myr to prevent either width from being zero. As with the 3-parameter fits, we also enforced a maximum total width of 3.5 Myr, which is clearly seen to be relevant for Wallner Core 4521 in Figure 11. This width cap is to prevent the fit from stretching all the way to t = 0 Myr for the noisier data, and was picked as the minimum width needed before the triangle shapes stopped changing dramatically.

Results for the Ludwig 851 core appear in Figure 10. We see in the bottom panels that the two widths have well-defined peak likelihoods, but broad distributions. In the upper left panel, we see that the two widths have some anticorrelation, as one might expect if they have to at least sum to the spread in the nonzero data points.

Zoom In Zoom Out Reset image size

Examining the best-fit shapes for the four cores in Figure 11, we can see that most of the data are best fit with a sharktooth shape that is intermediate between being fully symmetric or asymmetric. Only Ludwig Core 848 prefers a perfect sharp sawtooth shape, with the fit function capped by the minimum width initial parameter. The two Wallner cores prefer a slightly sawtooth shape, while Ludwig Core 851 actually has a slight preference toward a reverse sawtooth.

Zoom In Zoom Out Reset image size

It should be noted once again that the sharktooth fits are not statistically comparable to the 3-parameter fits in Section 3.2. However, in view of the lack of strong preference for a specific 3-parameter fit, it is an interesting possibility to explore further. Given the variations in the current data, we cannot meaningfully pick a preferred shape for the ⁶⁰Fe pulse, so the pulse shape does not provide significant information on the underlying astrophysics. However, the different pulse widths extracted using the different preferred fit shapes can be used to make some inferences on the astrophysical processes inside the supernova remnant.

Having examined the cores individually, we now turn to global fits that use all the cores to analyze jointly the ⁶⁰Fe deposition history. To do this, we first note that the ⁶⁰Fe flux time profile (i.e., its shape) should be common to all of the sediments regardless of their location on the Earth. We therefore fit all of the samples using a single pulse shape and thus the same width parameter σ_t .

We account for potential systematic differences between the samples by allowing the other fit parameters to vary independently for each core. To include possible systematic offsets in the absolute dating of the different samples, we fit different peak times for each core; this would correspond to shifts in the inferred time history between the samples. Small differences in the peak times of individual fits would suggest that differences in the dating are small, but large differences would point to the presences of unaccounted systematics, or the use of a bad fitting function. We also allow for differences in the peak heights, which could arise due to different infall and uptake at different locations. To perform these fits, we use a Markov Chain Monte Carlo approach, which enables us to search efficiently the nine-dimensional parameter space using the emcee package (Foreman-Mackey et al. 2013). ¹⁴

Figure 12 shows our results for the case of a sawtooth time profile; the best-fit parameters and other statistics are given in Table 4. Turning first to the width, we see that the likelihood is zero until about 2 Myr, then rises until the highest allowed value. Thus, as we have seen with the individual fits in Figures 2–5, the global sawtooth form only sets a lower limit to the signal width. We find that the 95% confidence lower limit to the FW at 0.1 maximum is 2.7 Myr. This limit is significantly larger than the time span of nonzero data points shown in Table 1, showing that the requirement of a sawtooth form leads to wider deposition time.

Zoom In Zoom Out Reset image size

Table 4. Global Fit Parameters

		Fit Type
Parameter		Gaussian	Sawtooth	Triangle
Width parameter σt	[Myr]	0.47 ± 0.05	3.18 ± 0.26	0.99 ± 0.10
Full width at 0.1 max Δt	[Myr]	2.00 ± 0.21	2.86 ± 0.23	1.8 ± 0.2
Δt 95% CL lower limit	[Myr]	1.7	2.7	1.6
Peak time tpeak [Myr]	L-848	2.16 ± 0.10	2.7 ± 0.05	2.1 ± 0.09
	L-851	2.44 ± 0.11	3.2 ± 0.2	2.4 ± 0.9
	W-4521	2.12 ± 0.11	3.3 ± 0.5	2.2 ± 0.9
	W-4953	2.44 ± 0.05	3.5 ± 0.2	2.5 ± 0.4
Peak flux Φpeak [atoms cm−2 kyr−1]	L-848	0.037 ± 0.006	0.039 ± 0.021	0.042 ± 0.007
	L-851	0.116 ± 0.019	0.13 ± 0.04	0.12 ± 0.02
	W-4521	3.89 ± 0.41	6.0 ± 1.9	4.5 ± 0.5
	W-4953	3.19 ± 0.42	2.7 ± 0.5	3.6 ± 0.5
Fluence [107 atoms cm−2]	L-848	0.043 ± 0.07	0.060 ± 0.023	0.041 ± 0.007
	L-851	0.14 ± 0.03	0.19 ± 0.05	0.12 ± 0.02
	W-4521	4.4 ± 0.6	9.5 ± 3.2	4.4 ± 0.6
	W-4953	3.7 ± 0.4	4.4 ± 0.8	3.6 ± 0.3
Goodness of fit		−195.5	−185.6	−196.8

Note. Fits with a single width parameter but varying peak times and peak flux values for different cores. L denotes Ludwig et al. (2016), W denotes Wallner et al. (2016). We list below the best-fit values and 1σ uncertainties. The peak flux is given in units of 10⁴ atoms cm⁻² kyr⁻¹, and the fluence in units of 10⁷ atoms cm⁻². For the goodness of fit , bold indicates the highest likelihood, and italics the close second.

Download table as:  ASCIITypeset image

The shapes of the peak times are similar to those in the individual core fits in Figures 2–5, and show the same asymmetries. We find that the most probable peak times differ significantly among the samples. This reflects the fact that, for a sawtooth, the peak is also the point with the earliest nonzero ⁶⁰Fe counts, itself an artifact of the sampling. These differences between the peak times (upper right panel of Figure 12) are also similar to those in the individual fits.

The peak ⁶⁰Fe flux values (lower left panel) and the ⁶⁰Fe fluences (lower right panel) measured in the two Wallner samples are similar, and also those of the two Ludwig samples. However, we do see significant differences in the peak ⁶⁰Fe flux values (lower left panel) and the ⁶⁰Fe fluences (lower right panel) between the Wallner and Ludwig measurements. As noted above, this may be due to differences in the analysis techniques, and also perhaps due to differences in the deposition rates caused by varying infall or uptake factors. We note that the global sawtooth fits for the Wallner cores give values that are slightly higher than in the individual fits, but well within errors.

We have also performed global fits for the Gaussian and triangle fit functions. The results are summarized in Table 4, and the Gaussian fit is shown as Figure 14. For the Gaussian case, we see that the FW at 0.1 maximum of 2.00 ± 0.21 Myr (upper left panel) is a compromise between the values of the Gaussian widths found in individual sample fits as seen in Table 3. On the other hand, the triangle case gave a well-determined width, similar to the global Gaussian fit and the individual triangle fits. We find that, for the Gaussian fit, the 95% confidence-level lower limit to the FW at 0.1 maximum is 1.7 Myr, while for the triangle fit the same limit is 1.6 Myr. The upshot is that, for these other functional forms, the global fits once again give a long timescale for the ⁶⁰Fe deposition. Thus, using all of the data, we find that the pulse timescale is at least

which is also close to the interval between the earliest and latest nonzero ⁶⁰Fe counts across all crusts shown in Table 1. ¹⁵ Any model for ⁶⁰Fe delivery must account for this timescale.

For the Gaussian and triangle global fits, we also find that the peak time and peak flux values for each sediment are similar to those for the corresponding single-sediment fits seen in Figures 2–6. The differences among the peak times are relatively small, with significant overlaps between the different fits. This suggests that the absolute dating of the various samples does not suffer from significant unknown systematic uncertainties. This also supports our conclusion that the differences in peak times in the sawtooth case (Figure 12) are due to the abrupt onset of that fitting form, and the sampling.

Turning to the fluence, Table 4 shows that the Gaussian and triangle results are very similar, indicating that the integral nature of the fluence is not very sensitive to the differences between these fitting functions. On the other hand, the sawtooth case gives substantially higher fluences for all sediments. Since the sawtooth fluence is given by Φ_peak σ_t /2, and the peak flux values are consistent across different fitting function, it is clear that the timescale σ_t is the source of the discrepancy. Indeed, we have seen that the sawtooth timescale is poorly determined by individual fits in Figures 2–5. For this reason, we see that determinations of the fluence are model-dependent given the current data.

In sum, we see that when all of the sediment data are combined; the global fits find that the ⁶⁰Fe fallout time is long, ≈2 Myr. We recall that, as discussed in Section 3.3, this duration is not a geophysical artifact, but reflects the underlying astrophysical fallout time. We now turn to the interpretation of this result.

In this picture, the dust grains move with the gas, so that the grain velocity is the same in magnitude and direction as the plasma bulk velocity v_gr = v_gas, i.e., the motion is radial, overlaid with perturbations due to turbulence. In addition, we assume the mass density ρ_dust of dust grains is always proportional to the gas density ρ_gas ≈ m_p n_gas. This model thus posits a direct proportionality between the mass fluxes of the grain particles and the gas: J_gr = ρ_dust v_gr ∝ ρ_gas v_gas. Thus, the time history of the grain flux is determined by that of the gas flux.

This model is the one adopted by Fry et al. (2015), whose key results we summarize here. The model focuses on the Sedov phase of the remnant, in which the blast radius evolves as follows as a function of time t: , where E is the kinetic energy of the ejecta into a uniform-density medium with , and β = 1.1517 for monatomic gas, with adiabatic index γ = 5/3. Inverting this relation, we derive the following estimate of the blast arrival time at radius r:

where we scale using benchmark estimates of the Local Bubble density and total blast energy (see, e.g., Fry et al. 2020). The corresponding speed of expansion is

Just behind the shock, the gas density is always the same constant multiple of the ISM density, namely ρ_gas = 4ρ_ism, for γ = 5/3. However, because the mass of the supernova ejecta is fixed initially, the ejecta density must drop as the blast volume grows. To fix numbers, we assume that the ejecta are entrained with the gas, and approximate the blast as a thin uniform-density shell of fractional width x = Δr_shell/r_blast ≪ 1, in which case mass conservation implies x ≈ 1/12.

The duration of the flux is the timescale for the blast shell to pass by. At a fixed distance r, and using self-similarity, the shell crossing time is

This is nearly 2 orders of magnitude smaller than the observed ⁶⁰Fe pulse width.

The ISM density would have to be ≳100 cm⁻³ in order to overcome this discrepancy. Such a density is characteristic of a giant molecular cloud complex, which could have been a plausible density if our ⁶⁰Fe-depositing supernova were the first to explode in a nearby star-forming region. However, modeling of the Local Bubble indicates that it has hosted multiple supernovae over timescales of 3 Myr or longer (Smith & Cox 2001; Breitschwerdt et al. 2016), and there is now clear evidence for ⁶⁰Fe deposition by an explosion ∼7 Mya (Wallner et al. 2021). After the first explosion, the local interstellar region would have a much smaller density, so we are driven to consider other explanations for the long ⁶⁰Fe timescale.

It is also worth highlighting that terrestrial (anthropogenic) explosions do not exhibit efficient mixing of ejecta with the larger blast wave. Although both move outward rapidly, the ejecta (i.e., the material responsible for the explosion) remains confined relatively near the center of the explosion while the forward shock travels much greater distances and without carrying ejecta material with it. For example, in the case of the Chelyabinsk bolide, the larger meteorite fragments were found in a 32 km long, 10 km wide region along the original trajectory of the meteor (Popova et al. 2013). Meanwhile, the smaller aerosols and dust particles formed a cloud that rose vertically 11 km over 80 s then stabilized at that altitude for the next 120 s (Gorkavyi et al. 2013). In contrast, the shock wave produced by the bolide propagated radially outward traveling 23 km in 76 s and 52 km in 173 s (Popova et al. 2013) showing a definitive decoupling between the shock wave and ejected material. Such observations of decoupling further constrain the entrainment model's validity.

Several possible explanations of the long ⁶⁰Fe deposition timescale have been discussed elsewhere. (1) The prolongation of the signal over time could reflect multiple supernovae, each with a narrow pulse width (Breitschwerdt et al. 2016). As we have noted, the data does not demand this, but also cannot exclude multiple supernovae in the 3 Myr peak. However, because we observe a broad ⁶⁰Fe pulse ∼3 Mya, preceded by a gap and and another pulse ∼7 Mya that also seems broader than suggested by the entrainment model (12), this scenario would require two bursts of supernovae in rapid sequence, with a well-defined lull in between. (2) The ⁶⁰Fe flux could reflect the motion of the Sun through the supernova material (Chaikin et al. 2022), with a complex time history needed to accommodate the two broad pulses. (3) Another possibility is that some ⁶⁰Fe-bearing dust was trapped by the Local Interstellar Cloud, an ∼5 pc feature that envelops the solar system (Koll et al. 2019; Linsky et al. 2019). One can quantify this by considering the stopping power of the cloud. For dust grains the size a_gr ≳ 0.3 μm needed to overcome solar radiation pressure, the stopping distance due to drag is 60 pc for n_LIC = 0.2 cm⁻³. Thus we do not expect drag to efficiently stop the grains unless they are much smaller. (4) Opher & Loeb (2022) model the effects of a neutral cloud (seeded with ⁶⁰Fe) passing through the solar system, and show that, if the cloud density is very high, it can compress the heliosphere within 1 au. If the cloud is also large enough, the passage can last for the required time.

We consider these scenarios to be worthy of further investigation, but also not without their challenges. Here we propose another solution motivated by our recent work (Fry et al. 2020), assess its merits and drawbacks, and offer observational tests.

We propose that the ⁶⁰Fe arrives at the Earth (and Moon) as part of charged dust grains that were created in the supernova, whose propagation was largely determined by the magnetic structure of the remnant impact into the surrounding medium. Fry et al. (2020) performed detailed calculations of the propagation of charged dust in a supernova remnant, motivated by the ⁶⁰Fe data. Here we summarize the rich physics influencing dust grain evolution and propagation.

Dust formation in supernova remnants is a subject of intense ongoing research, but it is clear that a substantial amount of dust is formed very soon after the explosion, e.g., SN 1987A shows infrared emission consistent with all of the supernova-produced iron being locked into grains within tens of years after the explosion (Matsuura et al. 2011, 2017, 2019), while recently Niculescu-Duvaz et al. (2022) examine a large sample of supernovae finding that, after ∼30 yr, on average of material is condensed into dust. We therefore follow Fry et al. (2020) in assuming that dust is present, likely with a range of grain sizes, and initially entrained with the gas from which it formed. Thus a range of dust compositions, sizes, and velocities is present.

After the dust is created, it suffers collisions with gas particles that lead to drag, sputtering, and charging. However, these effects are minimal because the dust is comoving with the gas. But as the supernova remnant evolves, a reverse shock propagates inward. This shocks and slows the gas, thereby decoupling the dust from the gas. Grains suffer some damage at the shock, but those large enough to survive crossing the reverse shocks will then be subjected to increased drag, sputtering, and charging.

Given the negligible magnetic field in the inner portion of the supernova remnant, the dust still travels radially. However, when the dust grains subsequently encounter the magnetized ISM, it acts as a mirror and reflects the dust. The dust grains then pass back through the remnant until they encounter the ISM material once more, and are again reflected by its magnetic field (at each ISM encounter, there is also some probability that the grains become trapped). The resulting dynamics is that the dust is confined to the ejecta region, with repeated bouncing motion, "pinball" style, in the ejecta interior. As they propagate, the dust particles are slowed by drag and are sputtered, becoming smaller and losing mass.

An order of magnitude calculation illustrates the key features of dust evolution found in the Fry et al. (2020) simulations. We model a dust grain as a sphere of radius a, density ρ_dust, and mean atomic mass m_gr. As the dust particle moves through the gas, it suffers collisions at a rate

where in the second expression we approximate the collision cross section with the geometric cross section, and assume that the dust is moving much faster than the gas, so that the relative speed moreover, we assume that the gas is dominated by hydrogen.

Drag has collisional and Coulomb components. In the assumed limit of fast dust particles with , the collisional term is

where v _rel = v _gr − v _gas is the grain speed relative to the local gas speed. This gives a collisional stopping time of order

where M_gr is the grain mass. The Coulomb drag force for high speeds larger than the plasma thermal speeds v² ≫ kT/m_p is independent of grain size: , in contrast to the collisional drag; here is the so-called Coulomb logarithm. The Coulomb stopping time therefore exhibits a different and stronger dependence on grain size:

If the grain is indeed moving much faster than the gas, then sputtering decreases the grain radius at a rate

where Y is the energy- or velocity-averaged yield of atoms liberated per collision, and m_gr is the mean mass of a grain atom. The resulting timescale for grain erosion due to sputtering is

We see that, in the high-velocity regime, the sputtering and drag timescales are related by a factor that depends only on the grain atomic mass and yields. Estimating these values for iron-bearing grains shows the timescales to be comparable, with drag somewhat faster. We thus expect significant drag and sputtering both to occur, and that the grain lifetime is similar to the drag timescale.

Figure 13 shows sample simulations of the trajectories of 0.1 μm dust grains (blue) encountering the magnetic field in the ISM (red). The left panel shows an example where the dust grain bounces (is reflected) straight back from the ISM, and the central panel shows an example where the dust grain is trapped temporarily by the ISM magnetic field, but eventually escapes. These exemplify the "pinball" feature discussed in Fry et al. (2020), whereby the directions of the dust grains' trajectories can be changed, losing memory of the location of their progenitor. Finally, the third panel shows an example where the dust grain is trapped for the duration of the simulation. All three panels illustrate how the timescale for the deposition of live radioisotopes can be extended due to interactions with the ISM.

Zoom In Zoom Out Reset image size

Determining quantitatively the conditions when a grain will be reflected versus trapped (and for how long) is a subject requiring further examination. A grain's charge depends on its composition, speed, and the surrounding environment; the dynamics of the magnetic field is determined by the supernova remnant's plasma dynamics. While the grain's possible motions are understood (i.e., curvature drift, gradient drift, reflections, etc.), when those occur is not fully characterized. The time the grain first encounters the magnetic field, the pitch angle of the grain's velocity with the local magnetic field, the dynamical timescale of the magnetic field, and the size of any turbulent eddies and density perturbations can influence the confinement of the dust grains. Fry et al. (2020) provided an initial statistical result for metallic iron grains, but a more detailed examination is beyond the scope of this work.

To summarize, our model builds upon the results of Fry et al. (2020) to predict the following supernova dust grain history and dynamics as the supernova remnant evolves.

• 1.  
  Free expansion phase: coupled. Dust grains are nucleated in dense ejecta knots, comoving with gas (Fry et al. 2020).
• 2.  
  Start of Sedov phase: decoupling. As the ejecta (composed of ejected gas and dust grains) experiences the reverse shock and decelerates, the surviving dust grains are dynamically decoupled from the gas. Their subsequent propagation is dominated by the Lorentz force and drag. In weakly or nonmagnetized supernova material, their motion is radial until they encounter the magnetized ISM (Fry et al. 2020).
• 3.  
  Sedov and snowplow phases: reflection and trapping. In encounters with the magnetized ISM, grains may be either reflected or trapped. The reflected grains traverse the inner supernova remnant until their next encounter with the ISM. Grain size plays a crucial role here: the smallest grains (a ∼ 0.005 μm) are trapped at the first encounter, sputtered, and lost, whereas the largest grains (a ∼ 1 μm) travel past the supernova–ISM boundary and then are trapped. Intermediate-sized grains (a ∼ 0.05–0.1 μm) have non-negligible probabilities both to be trapped and to escape. Thus over time the grain number density in the supernova material will drop, in favor of a buildup of trapped grains at the supernova–ISM boundaries. The trapped dust motion depends on the grain size and potential (i.e., the charge-to-mass ratio).
• 4.  
  Fadeaway phase: shell buildup and release to ISM. Trapped grain motion in turbulent field lines will be approximately diffusive, leading to a buildup in a shell at the supernova–ISM boundary. Spatial and time changes in the magnetic fields can lead to grain deceleration and acceleration, in addition to the action of drag. The grains are stopped over the drag timescale. The larger grains survive and remain as the forward shock slows to a sound wave, and the supernova remnant fades.

In this picture, the spatial distribution of dust grain changes over time. While the dust grains are predominantly reflected, the grains should roughly uniformly populate the inner supernova remnant. As they become trapped in the surrounding magnetized ISM, the grains should move diffusively in a shell of increasing thickness around the inner remnant. Thus, when the grain-bearing material arrives at the solar system, we expect the particle density to be a mix of a uniform and shell profile, with the shell profile more favored at late times and large distances. Detailed modeling of this distribution, and the resulting ⁶⁰Fe time profile at Earth, is beyond the scope of this paper, but is a problem we intend to revisit in future work.

The three parameters for all fits are θ = (t_peak, σ_t , Φ_peak); these are the time of the peak, a measure of the width, and peak flux, respectively:

The FW0.1M measures, Δt, of the FW are given by

where ε = 0.1 is a somewhat arbitrary choice but corresponds roughly to the level of the background to the measurements.

Finally, for each functional form, we calculate the corresponding fluence, :

We also use an asymmetric triangle or sharktooth shape, with separate width parameters and before and after the peak, respectively. This function allows us to probe the degree of asymmetry or "lean" in the signal:

This encompasses both our sawtooth fit () as well as the triangle fit , as well as a reverse sawtooth shape. The FW0.1M, Δt, and the fluence, , are defined by analogy with those for the 3-parameter triangle.