X-Ray Energy Deposition Model for Simulating Asteroid Response to a Nuclear Planetary Defense Mitigation Mission

The energy deposition model described in this paper strives to encompass as many realistic scenarios as possible spanned by the relevant parameters listed below.

Materials—Silicon dioxide (SiO₂), forsterite (Mg₂SiO₄), ice (H₂O), and iron (Fe) were chosen to represent the NEO materials. Though a potential impactor's composition would be much more complex, the materials listed are common components of more intricate minerals and expected to behave similarly. Bulk densities of asteroids with diameters less than 50 km (of interest for mitigation missions) are clustered in the 1–3 g cm⁻³ range. These can be reasonably well represented with our set of materials by including either macro- or microporosity (or a combination; Carry 2012). For deflection simulations, previous work by Dearborn et al. (2020) has shown that final velocities are relatively insensitive to minor differences in the type of silicate EOS used (e.g., using hydrated forsterite produced similar outcomes to SiO₂). The largest compositional lever on deflection outcomes is the atomic mass of the elements, which controls the X-ray energy deposition length scale. For this reason, we choose materials that are reasonable low-Z (H₂O, ice) and high-Z (Fe) end-members for materials present at asteroid or comet surfaces, along with two mid-Z rocky materials (SiO₂ and Mg₂SiO₄). Each selected material has well-established EOS data to inform the radiation transport simulations.

Porosity—The energy deposition model was developed using simulations with microporosities (ϕ) ranging from 0% to 80% for all materials, which covers the vast majority of possible NEOs. Extrapolation to higher porosity is possible but should be attempted with caution. Larger macroporosity features are better suited for full-scale hydrodynamic models and are excluded to keep the size and scope of the study manageable.

X-ray source spectrum—The NED thermal radiation can be approximated by a blackbody (BB) source with a peak at ∼4 keV according to Sections 7.73 and 7.79 of Glasstone & Dolan (1977). For simplicity and to provide multiple options, the NED sources we simulate are approximated as 1 and 2 keV BBs with spectra pictured in Figure 1 and peaks at ∼2.8 and ∼5.6 keV, respectively.

Zoom In Zoom Out Reset image size

X-ray source strength—The X-ray fluences (f) that were modeled range from 10⁻⁴ to 1 kt m⁻². The lowest value represents a small surface layer melt, which is well suited for the edges of illumination in a deflection mission. The highest point in the fluence range roughly corresponds to the energy per area seen by the closest point of an NEO when a 1 Mt NED is detonated 10 m from the surface, a realistic disruption scenario.

Source duration—The NED's radiation is released over a length of time (t_end) that depends on the number of neutron generations it takes to deliver the device's yield. The fast neutron generation time in fission is approximately 0.01 μs or 10 ns. Section 1.57 of Glasstone & Dolan (1977) sets a lower limit of the NED source duration as ∼10 ns (the last e-folding of the energy production) and estimates that the vast majority of the yield in a 100 kt device would be released within 70 ns. To encompass a broad range of scenarios, NED sources ranging between 10 and 100 ns were modeled and included in the deposition formulation. All of the sources were approximated as constant/having no time dependency, though uncertainties arising from this assumption are explored at the end of the paper.

To keep the scope of this study manageable and not introduce additional layers of complexity, all test asteroids simulated in this study will be homogeneous spheres. Additional features such as macroscopic variation in properties and irregular shapes will be the subject of future work and are discussed in Section 6.

The radiation-hydrodynamics code selected to simulate the process of X-rays impinging on an NEO surface was the Kull multiphysics code (Rathkopf et al. 2000). Developed at LLNL primarily for simulating inertial confinement fusion processes at the National Ignition Facility, the code is multidimensional, mesh-based, and able to accommodate many types of physics, including multiple radiation transport options: diffusion, discrete ordinates (also known as S_N ), and implicit Monte Carlo (IMC).

Of the three, IMC was selected as the transport package of choice for the simulations described in this study. This was done because the simulation has an angularly homogeneous point source of radiation situated in a vacuum. The diffusion approximation loses all angular information about the direction of radiation propagation and calculates the radiation flux via Fick's law using a diffusivity proportional to the inverse of the opacity. Diffusion thus has particularly large errors when simulating radiation in a vacuum and would deposit energy even in the shadowed part of the asteroid. The S_N method computes the propagation of radiation along specific directions chosen during problem setup; these specific directions are called discrete ordinates in the nomenclature of the method. In nonopaque materials, the radiation propagating in the specific directions will diverge as they move away from a point source, leaving regions of anomalously low radiation intensity between them. This discretization error is known as a "ray effect." The ray effects introduced by S_N transport would have created nonsmooth energy deposition on the asteroid surface because in a vacuum, ray effects are not reduced by absorption or scattering (Castor 2004).

Since it was designed for high-density and high-energy conditions, Kull would occasionally have difficulty simulating expanding, low-density material, such as would be required in a nuclear mitigation mission. In some instances, depending on the EOS data, Kull would pull the material into tension, which decreased the simulated blowoff momentum in an unrealistic manner. For those cases, the Kull simulation could be handed off to a similar mesh-based hydrodynamic code better suited for the conditions, before tension arose, in order to complete the full simulation with integrity. The handoff code, Ares (Darlington & McAbee 2001; Bender et al. 2021), was also used for comparison and cross-checking pure hydrodynamic runs initialized with the energy deposition model.

For Kull simulations where porosity was nonzero, a simple snowplow model ⁴ was used to simulate the material's behavior under compression. Ares utilized a more sophisticated model that lowers porosity when the pressure is above a corresponding crush curve pressure, but a similar snowplow model was also implemented for cross-code comparison. There were no distinguishable differences in the Ares results between the two porosity models, which lends confidence to the effectiveness of the snowplow model in this scenario.

2.2.1. EOSs and Initial Conditions

The full, 2D geometrical setup of our simulations with all relevant variables can be found in Figure 2. Because the problem is cylindrically symmetric around the axis connecting the NED source and the center of the NEO, a 3D approximation can be achieved by rotating around the same axis. The NEO 2D mesh was constructed using a number of equally spaced wedges fanning out from the θ = 0 axis. The number of required wedges would vary based on the HOB and the simulation being used. The wedges would then be radially segmented into ratioed zones based on the specified (very small) thickness of the surface zone and the zone count. The number of required zones depended on the code used, the material-opacity combination, and the flux of energy arriving from the source.

Zoom In Zoom Out Reset image size

A 1D simulation setup is similar to a single wedge in the 2D geometrical setup, except the wedge shape is exchanged for a 1 cm² column. The ratioed radial zoning is kept the same, and the end of the mesh is cut off between 10 and 100 cm, depending on how far the energy propagates in the given simulation time, rather than a column length equal to R. The β between the NED radiation and NEO surface is also preserved.

For every mesh parameter, a 1D resolution study was performed based on material and NED fluence level to ensure the accuracy of the results to within a few percent. Each study involves running a suite of simulations with a varying mesh parameter to determine a resolution value where the result converges. Though the number of zones in each direction matters, the most sensitive parameter was the thickness of the surface zone. Results from a suite of 1D IMC simulations at normal incidence for each material at a range of fluences can be seen in Figure 3, which convey when the surface zone thickness converges in each case. The results were drawn from cases using a 1 keV BB spectrum, a source duration of 10 ns, and zero porosity. For a 2 keV BB or nonzero porosity, the required zone thickness would be slightly coarser, while a longer source duration would require slightly finer resolution. All simulations were run with a surface thickness resolution smaller than would be required for convergence.

Zoom In Zoom Out Reset image size

As previously discussed, the primary metric used in this study for comparison between simulations is the amount of momentum carried by the mass ejected from the surface of the asteroid: the blowoff momentum, or p_b . By Newton's third law of motion, the blowoff momentum is closely coupled to the momentum imparted to the full asteroid by the NED and would ultimately determine the change in trajectory, Δv, which is what matters most in a mitigation mission. For our purposes, the blowoff momentum is a variable local to the surface of the asteroid, which allows for smaller simulations.

A zone would be considered "blown off" if its material is both melted and moving faster, on average, than the NEO's escape velocity (v_esc) in the radial direction. To be labeled melted, a zone's temperature would be higher than the material's melt temperature at the local density, which was given by a melt curve taken from the EOS data. As the material in a blowoff momentum zone is often moving much faster than an NEO's low v_esc, the melting criterion was the dominant requirement. All 1D simulations were run with the v_esc criterion determined using a 300 m NEO. Sensitivity studies indicated that the v_esc threshold would begin to compete with the temperature requirement on NEO size scales similar to the Moon (R ∼ 1700 km).

The suites of 1D simulations were all run to a completion time of 10 μs (at least 100 times longer than the source duration), with the final p_b value in each used for comparison and analysis. In many cases, particularly at high fluence, p_b will still be increasing at 10 μs and could take several orders of magnitude more time to reach convergence. Thus, this study does not attempt to reach convergence with its simulations and makes no concrete claims concerning the "final" blowoff momentum or Δv value for any given asteroid. All p_b and Δv values that are listed are taken at t = 10 μs unless otherwise specified. Comparing values between simulations is done with the assumption that blowoff momentum trajectories of IMC and model-initiated hydro simulations will continue to behave similarly at longer timescales, and the difference between them will not noticeably increase. A detailed analysis supporting this assumption can be found in Section 4.2. Figure 4 depicts a series of example blowoff momentum versus time curves for the range of fluences explored in this problem using 1D simulations with both logarithmic and linear y-axis scaling. Though the logarithmic curves appear to be semiconverged, the linear data suggest otherwise.

Zoom In Zoom Out Reset image size

The basis of this energy deposition model utilizes energy density (per volume) profiles taken from a 1D IMC simulation at a specific time slice, shortly after the source is finished illuminating the surface of the material. The energy profile consists only of the internal energy imparted by the NED radiation. Any displaced meshing is returned to its original position before it is utilized for model construction. The optimal choice of time slice was selected by initializing separate, purely hydrodynamic 1D Kull simulations with the same mesh/material construction using energy density profile splines taken from IMC simulations and comparing blowoff momentum values between the two. Figure 5 shows examples of energy density profiles for silicon dioxide at various points after the source's illumination is finished. How well each time slice was able to recreate the blowoff momentum of the original IMC simulation with various materials at a range of fluences is shown on the right.

Zoom In Zoom Out Reset image size

Though the accuracy over the times is quite similar, 5 ns was selected for subsequent energy profile harvesting. At 5 ns, the vast majority (>90%) of reradiation has already occurred at fluence levels where it is a dominating effect, but fairly little energy has been turned into motion.

To ensure this method would eventually scale up to 2D, similar energy density profiles were collected from 2D IMC simulations and compared to the profiles collected in 1D. Figure 6 shows a comparison between the two sets of profiles, each taken at the previously selected t = 5 ns.

Zoom In Zoom Out Reset image size

2D simulations cannot match the level of IMC particle statistics per zone of a 1D simulation, so the profiles are substantially more noisy, especially at large angles of incidence. However, the results for each dimension still matched closely enough to lend confidence in the model's eventual implementation in two dimensions.

Furthermore, the symmetry of a homogeneous, spherical asteroid allows for 1D simulations to successfully reproduce the blowoff momentum of a full 2D simulation if integrated over the illuminated angles. As mentioned previously, the 2D rad-hydro simulations performed for this study are computationally expensive and often plagued by various mesh-tangling issues, which are absent in 1D simulations. In cases where an asteroid's Δv is needed more quickly than a proper spherical 2D model would take to run to completion (or if the simulation encounters mesh tangling), a representative set of 1D simulations spanning θ = 0​​​​​–${\theta }_{\mathrm{Max}}$ can be substituted instead. The set would require a sufficient number of simulations to properly resolve the blowoff momentum versus θ dependency in order for the calculation to be accurate. Each 1D simulation must also be set up to accurately match its 2D analog while correctly accounting for the change in angle of incidence ($A=\cos \beta$). With variables referenced from Figure 2, and where subscripts i and j refer to the surface zone's inner and outer edge compared to θ = 0, respectively, the source arrival time should be adjusted according to t = (S_i + S_j )/(2c). ⁵ If directly comparing between IMC simulations, the source duration should be increased as θ increases by adding the time the source takes to transverse the zone to the original source duration: t_zone = (S_i − S_j )/c. The fluence each zone sees should also be calculated based on ${YA}/(4\pi {S}_{\mathrm{avg}}^{2})$, with S_avg = (S_i + S_j )/2 and Y representing the NED yield. Once these adjustments are made, the blowoff momentum for each simulation can be constructed into a p_b (θ) curve and integrated over θ using the formula below to produce an aggregate 2D blowoff momentum:

$$\begin{eqnarray}&&{p}_{b,\mathrm{int}}=2\pi {R}^{2}{\int }_{0}^{{\theta }_{\mathrm{Max}}}{p}_{b}(\theta )\sin (\theta )\cos (\theta )d\theta .\end{eqnarray} \tag{ 1 }$$

Based on over 100 2D and 1D sweep simulation pairs, the discrepancy distribution between the blowoff momentum produced by each was centered at 1.2% with a standard deviation of σ = 4.8%. This technique is deployed in Section 5 to replace failed 2D simulations and serve as a quality check.

With a time slice selected and 2D implementation verified, construction of the energy deposition model could begin. The baseline energy density profiles used for assembling the model would be harvested from a suite of 504 1D IMC simulations, hereafter called the "original" suite. The simulations would cover the four different materials, both source spectra, nine different fluences spanning f = 10⁻⁴–1 kt m⁻² at regular logarithmic intervals (depicted in Figure 4), and seven source lengths (t_end = 10, 20, 30, 40, 60, 80, and 100 ns). All simulations would be run at a normal (A = 1) angle of incidence and with 0% porosity. Adjustments to the model to incorporate varying A and porosity (ϕ) are discussed in later sections. The energy density profiles collected from the original 504 simulations were carefully reviewed to determine an optimal way to fit their shape across fluences and source lengths such that any potential scenario could eventually be generated. After an extensive search through many different functional forms, a piecewise function for each of the different sections was eventually selected as the best-fitting option. A picture of a standard deposition profile can be seen in Figure 7 with each of the function's components.

Zoom In Zoom Out Reset image size

At shallow depths, the first section of the piecewise function encompasses the "thermal wave" portion of the energy profile, which is dominated by the heated and expanding material. Depicted in red, this portion was best fit by a semiellipse. The next two sections, in blue and yellow, depict a shock wave, where the ejection of expanding material creates a pressure pulse in the opposite direction that propagates deeper into the material. Despite appearances and after extensive testing, the blue tail's behavior was not steep enough to be accurately fit with an exponential function. Instead, a 1/x function was found to better describe the profile's shape. The shock-wave portion ends with a sharp cutoff for all materials and source spectra except for one case: ice with a 2 keV BB source. For this combination, at higher fluences, the thermal wave was strong enough to begin overtaking the shock wave, which required the varying of the length of the cutoff function to accommodate the longer tail. Aside from this exception, after accounting for continuity constraints, the profile only required six parameters to fully describe an energy profile's shape. The final parameters with explanations are found in Table 2.

Table 2. The Parameters Needed to Define the Energy Deposition Profile Piecewise Functional Fit (with Notes on Each)

Name (${{ \mathcal P }}_{i}$)	Formula	Notes
X-Scale (${ \mathcal X }$)	c1	Horizontal radius of the ellipse
Y-Scale (${ \mathcal Y }$)	c2	Scaling factor used to determine the height of the ellipse
Ex-Frac (${ \mathcal E }$)	$1-\tfrac{{c}_{4}}{2{c}_{1}}$	Linked to the decay rate of the tail; configured as a fraction to prevent the fitting routine from converging on unrealistic values
Square (${ \mathcal R }$)	$\tfrac{{c}_{3}}{2{c}_{1}-{c}_{4}}+{c}_{5}$	Height of the "square" beneath the half-ellipse; must equal the 1/x function at the horizontal ellipse diameter for continuity
Offset (${ \mathcal B }$)	c5	Represents the constant height of the shock wave; taken from the ${ \mathcal M }x+{ \mathcal B }$ linear function before ${ \mathcal M }$ was assessed to be unnecessary for the profile fits
Cutoff (${ \mathcal C }$)	c6	Depth where the shock-wave energy density sharply drops; selected prior to fitting and used as fixed parameter to stabilize the fitting results
Expit (${ \mathcal I }$)	c7	Defines how quickly the shock wave drops to 0 using the expit function in Python, $1/(1+{e}^{-{ \mathcal I }x})$; set to 104 in all cases except Ice-2 keV, where the parameter is allowed to float

Note. The formulae reference variables defined in Figure 7.

Download table as:  ASCIITypeset image

After settling on a functional form, the 504 profiles were fit using a nonlinear least-squares fitting routine in Python to obtain values for each parameter. These parameters were then collected and separated based on each material–source spectrum combination. Each of the 49 individual parameter sets ⁶ were then themselves fit to a 2D polynomial to determine the dependency on fluence (f in log₁₀(kt m⁻²) units) and source duration (t_end in ns). The 2D polynomial had 16 available parameters with the format

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal P }}_{i} & = & {a}_{0}+{a}_{1}{t}_{\mathrm{end}}+{a}_{2}{t}_{\mathrm{end}}^{2}+{a}_{3}{t}_{\mathrm{end}}^{3}+{a}_{4}f\\ & & +\ {a}_{5}{t}_{\mathrm{end}}f+{a}_{6}{{ft}}_{\mathrm{end}}^{2}+{a}_{7}{{ft}}_{\mathrm{end}}^{3}+{a}_{8}{f}^{2}\\ & & +\ {a}_{9}{t}_{\mathrm{end}}{f}^{2}+{a}_{10}{t}_{\mathrm{end}}^{2}{f}^{2}+{a}_{11}{f}^{2}{t}_{\mathrm{end}}^{3}+{a}_{12}{f}^{3}\\ & & +\ {a}_{13}{t}_{\mathrm{end}}{f}^{3}+{a}_{14}{t}_{\mathrm{end}}^{2}{f}^{3}+{a}_{15}{t}_{\mathrm{end}}^{3}{f}^{3}.\end{array}\end{eqnarray} \tag{ 2 }$$

In an attempt to scale back the potential 49 × 16 numbers generated from all the fitting, each of the profile fitting parameters (${{ \mathcal P }}_{i}$) in Table 2 were first run through a minimization algorithm before the final fitting. This algorithm found the best fit for each ${{ \mathcal P }}_{i}$ across all eight material–source temperature combinations for all possible numbers and combinations of a_i , such that a subset of the 16 polynomial parameters could be selected for fitting without significantly impacting the final accuracy. The final formats of the truncated 2D polynomials can be found in the appendices, and the final fitting parameters can be found in the python energy deposition script on GitHub at Burkey & Managan (2023).

With the profile fitting complete, a piecewise energy density function can be evaluated for any combination of material, source spectrum, source duration, and fluence. The accuracy of this initial capability was tested by regenerating functional profiles for the same suite of 504 1D simulations that informed the model's construction and using them to initialize a complementary suite of purely hydrodynamic simulations. Utilizing the generated profiles, the internal energy of the relevant 1D simulation surface zones was set at t = 0 by taking the average energy per volume from the piecewise function over the two zone edge depths for each zone with increasing depth from the surface until the function yielded 0. When the suite of function-initiated hydrodynamic (FXN) simulations was run to 10 μs and the resulting blowoff momentum values were compared with their IMC counterparts, the discrepancy percentage distribution was defined by an average and standard deviation of 8.1% ± 25.7%. A full discussion of the discrepancies, their potential causes, and various physics-based solution attempts can be found in Appendix B. The aforementioned study indicated that the time dependence of the material opacities was the dominant source of discrepancy, which could not be corrected by modifying the deposition function directly. Instead, a simple scaling factor was introduced to the functional fit, which would be used to match the FXN simulation's final p_b to the IMC result, indirectly correcting all different discrepancy sources simultaneously.

To develop a detailed dependency of p_b versus the new scaling factor, an additional suite of 1D FXN simulations was run, where each function of the 504 cases was scaled ⁷ from 0.3 to 1.6 in increments of 0.1. Using the blowoff momentum values from the scaled FXN simulation suite, p_b,Sc/p_b,1 versus the scaling factor versus fluence was fit to another 2D polynomial (Equation (2)), this time with all 16 terms. The source duration was found to have a negligible impact on the scaling factor's results and was not included in the fit. Utilizing the new 2D fitted function, all 504 scaling factors needed to adjust the FXN suite's p_b s to match the IMC simulations were extracted. In the same manner as the previous set of parameters in Table 2, the scaling factors were divided according to each material/BB combination and fit to another optimized, truncated 2D polynomial (${ \mathcal S }$, also found in Appendix A). With the new combination of profile-generating parameter functions multiplied by the scaling polynomial, the updated model was used to regenerate/initialize another set of 504 1D hydrodynamic simulations. A histogram depicting the distribution of percent differences between the IMC simulations' p_b and the scaled FXN simulations' p_b is shown in Figure 10 with the label "Original." The standard deviation of the collected original 504 residuals is σ = 3.9%, and the individual σ for each material/BB combination can be found in Table 3.

In order to initialize either a 2D or 3D asteroid, the energy deposition model from the previous section would require another functional dependency on the angle of incidence. Fortunately, on a fundamental scale, the only change between material receiving X-rays at normal incidence and at an angle is the percentage of reflected photons. Therefore, no changes to any of the functions from the previous section would be required to incorporate the angle of incidence. Instead, the impinging fluence used to generate the energy density profile would be scaled down by the fraction of reflected photons. Thus, the angle of incidence dependency would add a separate function to the model. However, as the angle of incidence increases, a high enough percentage of photons are reflected such that the fraction of photons absorbed by the material and then reradiated is also significantly impacted. The IMC simulations were unable to distinguish between reflected and reradiated photons, so another, indirect method of disentangling the change in fluence was pursued.

An additional product of the original model's development was enough data to inform a detailed functional dependency between fluence/source duration and blowoff momentum for each material/BB combination Once this dependency was fit (again using Equation (2)) for a specific source duration, the new 2D function could use the p_b from an A ≠ 1 incident-angle simulation to determine what fraction of fluence the surface "saw" compared to it is A = 1 counterpart, thus determining the effective reflected photon fraction. To quantify the reflected photon fraction versus A behavior, another suite of 768 1D IMC simulations was run covering four non-normal angles of incidence (A = 0.8, 0.6, 0.4, and 0.2), eight fluences (all from the original suite except 10⁻⁴ kt m⁻²), three source durations (10, 40, and 100 ns), and the eight material/BB combinations. The lowest fluence was left out of the suite due to the function-generating fluence being pushed below the range used to build the energy profile generation by the angle of incidence. The p_b values at 10 μs for each of the new suite's simulations were input into the aforementioned f versus p_b and t_end 2D polynomial function to return the fluence absorbed by the surface without the reflected photons. The new f values were divided by f_A=1 to determine the fluence scaling factors based on angle of incidence. As with all of the previous parameters, the 768 fluence scaling factors were separated into the eight material/BB combinations for fitting. For reference, Figure 8 contains several sample surface plots of the fluence scaling factors (each material represented) based on f and A. To account for the now three dependent variables, Equation (2) was utilized for fitting again, with t_end swapped out for A and expanded by an additional 10 terms to account for t_end. The final functional dependency was again minimized to find the most important of the 26 parameters. The resulting truncated, 12-term 3D polynomial (Θ; format found in Appendix A, Equation (A1)) was used for fitting the eight sets of fluence scaling factors. With the angle-dependency function in place and using the same methodology as the original suite, another new set of FXN simulations were run, with p_b values collected at 10 μs. The resulting distribution of percent differences between the IMC and FXN simulations' p_b has a standard deviation of σ = 5.9% and is also shown in Figure 10 with the label "Angle."

Zoom In Zoom Out Reset image size

Though not included in the above analysis, for cases at high angles of incidence, where the scaled fluence falls below this model's range, the piecewise function at 10⁻⁴ kt m⁻² scaled down by f/10⁻⁴ kt m⁻² would be utilized as an approximation (Equation (A4)). In subsequent sections, the product of the original incident fluence and Θ will be called the angle-scaled fluence. It should also be noted that when constructing a 2D setup, the incident fluence is scaled down by $A=\cos \beta$ according to the Lambert cosine law. Θ does not incorporate this additional scaling, and it should be taken into account before generating energy deposition profiles.

Bulk density measurements from asteroids and comets are typically well below the grain densities of the materials that comprise these small bodies, which indicates that there is empty space enclosed in the material. The voids could come from either macroporosity, such as space between boulders, or microporosity, such as space between grains. In order to accurately model the X-rays' penetration into an asteroid, nonzero microporosity must also be folded into the energy deposition model, while macroporosity would be handled with a large-scale hydrodynamics code. The details of our microporosity simulation setup are discussed in Section 2.2. Based on preliminary IMC energy density profiles with nonzero porosity, scaling the depth and amplitude using the porosity fraction would return a close match to the original zero-porosity profile. At low fluences, the porosity-scaling method was extremely successful and only required adjusting the shock-wave depth to match the scaled original piecewise energy profile to the nonzero-porosity IMC data. However, with increasing porosity and fluence, some modifications to other parameters of the piecewise function were required to reproduce the IMC energy density profiles. Such discrepancies can be seen in the example set of original and scaled energy density profiles in Figure 9.

Zoom In Zoom Out Reset image size

The energy density profiles needed to find, fit, and construct the energy deposition model's dependency on porosity were taken from a third suite of 864 1D IMC simulations covering four nonzero-porosity values (20%, 40%, 60%, and 80%), nine fluences (values from the original suite), three source durations (10, 40, and 100 ns), and the eight material/BB combinations. Of this set, 786 ran to completion, with the failed simulations randomly scattered throughout the initial parameters. After studying the new suite of energy density profiles, it was determined that adding porosity dependence to the piecewise function parameters ${ \mathcal S }$, ${ \mathcal X }$, ${ \mathcal B }$, ${ \mathcal C }$, and ${ \mathcal I }$ would be sufficient to account for changes to the piecewise function and scaling. The 786 new energy density profiles were then fit using the same piecewise functional fitting routine from the original suite, this time with the unmodified parameters ${ \mathcal Y }$, ${ \mathcal E }$, and ${ \mathcal R }$ fixed to the original 0% porosity values. The fitted values for the remaining floating parameters were then divided by their corresponding original 0% porosity values to form more sets of scaling ratios. When combined with the 0% porosity data points (ratios = 1.0), each parameter's data set (separated by material/BB combination) could be fit to another functional form that would scale the original parameter function up or down based on the porosity value while still yielding the original value with zero porosity. Similar to the construction of the original ${{ \mathcal P }}_{i}$ functions, the porosity-adjusted parameters were fit to a 3D function dependent on source duration (t_end), fluence (f), and porosity (ϕ). This time, a set of artificial uncertainty values were input into the fit, with nonunity values having 1% uncertainty and 1.0 values having 0.1% uncertainty, which would ensure that the integrity of the original model remained intact. In all cases except for ${ \mathcal X }$, a 3D polynomial minimized for the least possible number of terms and with a₀ = 1.0 was sufficient to capture the behavior and fit to the ratio data. For ${ \mathcal X }$, a function with the form $1+\exp (\phi )$ was chosen instead.

Once the adjustments to the piecewise profile function were complete, they were used to initiate another set of FXN simulations, from which the blowoff momentum values were then used to fit/construct a similar functional porosity dependence for the scaling factor ${ \mathcal S }$. The full set of final parameter functions, with the porosity scaling included, can be found in Appendix A. With all of the adjusted parameters, including ${ \mathcal S }$, in place, a final suite of 786 FXN simulations was run to quality-check p_b with the IMC results. Figure 10 shows the distribution of percent differences in p_b between the two simulation sets with the label "Porosity" with a standard deviation of 6.6%.

Zoom In Zoom Out Reset image size

Though this energy deposition model covers a vast span of parameter space, the adventurous modeler may need to go beyond its bounds. For the occasions where extrapolation is necessary, the following paragraphs contain guidance as to how much caution should be taken.

The most likely extrapolation was already included in the final angle × porosity suite: a high angle of incidence, leading to a lower fluence than this model was constructed for. Especially for 2D simulations with angular zones at A < 0.05, some degree of accommodation is required. To address this need, for cases with angle-adjusted fluence below 10⁻⁴ kt m⁻², the energy deposition model will return the relevant profile generated at 10⁻⁴ kt m⁻² scaled down by f/10⁻⁴, as discussed in Section 3.5. The FXN simulations with an original fluence of 10⁻⁴ kt m⁻² but scaled down based on A in the angle × porosity suite were initiated in this way, and the final σ includes their discrepancy contributions. When only the f = 10⁻⁴ kt m⁻² cases are selected (n = 287), the FXN/IMC simulation p_b discrepancy distribution is averaged at −5.7%, with a standard deviation of 12.7%. These numbers are not substantially different from the full distribution's average and standard deviation, so extrapolation into this regime can be performed, when appropriate, with confidence that the overall result will not be significantly impacted. In addition, contributions to the total p_b are dominated by the region where A is near 1, as the fluence arriving at the asteroid falls with decreasing A.

On the opposite end of the fluence range, extrapolation to f > 1 kt m⁻² should be attempted with caution, and the results may significantly impact the total p_b of a 2D or 3D NEO simulation. However, to inform any modelers that wish to attempt an ultrahigh-fluence simulation, a small set of 1D simulations was run at 2 and 5 kt m⁻² and A = 1.0, with the remaining parameters spanning the breadth of the angle–porosity suite. The 2 kt m⁻² suite of simulations yielded fairly successful matching between the IMC and FXN runs, with the difference in p_b distribution centered at −1.8% with a standard deviation of σ = 11.3%. Within the 2 kt m⁻² suite, only the Ice-2 keV case would give cause for caution, with the model underestimating the IMC blowoff momentum by an average of −24%. Meanwhile, the 5 kt m⁻² suite was clearly beyond the model's capability for extrapolation, with the majority of material/source cases yielding an average overestimation of >20%. Only the cases using ice yielded average discrepancies that underestimated the blowoff momentum, though with magnitudes still large enough to cause concern. The authors recommend staying below 2 kt m⁻² if extrapolation cannot be avoided.

Available data suggest that NEOs with porosity of >80% are unlikely. However, to test the model's capabilities in this regime, 257 IMC/FXN simulation pairs with 90% porosity, which were an extension of the angle × porosity suite, were run to completion. The discrepancy distribution was characterized by μ = −7.37% and σ = 26.7%, which were collected by excluding approximately 20 egregious (>200%) outliers, likely due to a simulation failure. As a result, it is recommended to examine initial energy depositions carefully before proceeding with limit-exceeding porosity. A number of cases, particularly with either ice at any fluence or the other materials at high fluence, overestimated p_b by >50%.

The source duration explored in this study already encompasses the very shortest possibility, a single neutron generation (10 ns), and extends well beyond the 70 ns upper range recommended by Section 1.57 of Glasstone & Dolan (1977). No extrapolation beyond these bounds is recommended.

This model depends heavily on the assumption that if the blowoff momenta for an IMC and FXN simulation are matched at 10 μs, they will continue to match at much later times. To verify and quantify the validity of this argument, the angle × porosity suite of 1D simulations from Section 4 was extended to 50 μs. 70% of the simulation pairs reached the new completion time, with the blowoff momentum recorded at 10 μs intervals. Comparing the percentage of residual distributions, such as those in Figure 10 between the extended simulation set at the original matching point of 10 μs and separately at 50 μs, yielded means and standard deviations of μ_1k = −3.19%/σ_1k = 8.46% and μ_5k = −3.44%/σ_5k = 8.52%. The small changes between both values indicate that to 50 μs, the assumption for long-time matching appears solid.

However, if "convergence" is defined as a <1% increase in blowoff momentum over 10 μs, our estimations indicate that high-fluence cases reach that point at ∼1 ms. Meanwhile, the lowest-fluence cases have usually already reached convergence at the endpoint of this study: 50 μs. For the long-time assumption to hold for all cases, the rate of change between the IMC and FXN p_b should be low enough to not cause significant deviation between 50 μs and whenever convergence is reached. The change in percent p_b discrepancy between IMC and FXN simulations between the 40 and 50 μs data points is shown based on the angle-scaled fluence level in Figure 11. Fortunately, the rates of change converge closer to 0 at higher fluences, which require longer run times, compared to the lower-fluence cases, which have larger rates of change but likely do not require run times beyond 50 μs. The average absolute values of the change rate for each of the five fluences (starting with the lowest) are 2.48 × 10⁻⁵, 8.68 × 10⁻⁶, 5.10 × 10⁻⁶, 3.37 × 10⁻⁶, and 3.74 × 10⁻⁶ % ns^–1. Therefore, assuming that these rates remain constant and using the highest-fluence average, a function-initiated simulation may have diverged ∼3.5% from its IMC counterpart when completed at 1 ms. In addition, the projected convergence time of the middle- and lower-fluence cases would be short enough compared to the highest-fluence cases to negate any marginal increase in average percent change rate. Finally, when the absolute values of the percent change rates for all cases were collected at various time intervals, the means/standard deviations were 10–30 μs, μ = 2.9 × 10⁻⁵/σ = 1.2 × 10⁻⁴; 30–40 μs, μ = 1.3 × 10⁻⁵/σ = 3.9 × 10⁻⁵; and 40–50 μs; μ = 0.9 × 10⁻⁵/σ = 2.7 × 10⁻⁵ % ns^–1. These results clearly indicate that the percent change rate is decreasing over time and may eventually converge to 0. Therefore, 3.5% conservatively estimates the error caused by change over time; in reality, it would likely be much lower. Even without the decreasing change rates, 3.5% is well within this model's margin of uncertainty. Based on these arguments, the long-time convergence assumption appears valid.

Zoom In Zoom Out Reset image size

This section aims to summarize the various potential sources of error within this model and condense them into a "total" uncertainty.

The first possible way of introducing error into the model is via insufficient resolution in simulations. Though both sets of 1D IMC and FXN simulations were run with either exactly matching or slightly modified meshes, having an insufficiently resolved mesh will cause the results to deviate from the "true" result. As shown in Figure 3, the p_b values for an underresolved simulation would likely be artificially higher than a fully resolved simulation. Overestimating the effect of a mitigation mission is much more dangerous than underestimating, so a very thorough resolution study was run before model development began (see Section 2.3 for details). However, as simulations with non-normal angles of incidence and nonzero porosity were folded into the model, the expected fully resolved resolutions were extrapolated from the A ≠ 1 scaled fluence and the ϕ > 0 scaled depth rather than performing another full-resolution study. Though this extrapolation is reasonable, a possibility remains that some calculations are underresolved. Also, for simulations that failed to reach 10 μs for any code-related reason, slightly reducing the number of zones would often resolve the issue, introducing further potential for slightly insufficient resolution. Despite making every effort to run simulations at a resolution over what would be required for convergence to the fully resolved p_b , we conservatively estimate that the potential error introduced based on insufficient resolution in any of the 1D simulations used for model development is at most 5%.

Even with a fully resolved simulation, some uncertainty is also introduced via the way the physics packages of the simulation are set up, most notably, the material models and melt definitions. In particular, the regime this study samples (isochoric heating to vaporization and expansion into vacuum) is not well represented in experiments that inform the EOS, thus leaving the thermodynamic data poorly constrained and likely to vary depending on how the EOS is constructed. To come up with a rough estimate of how implementing a different EOS model might affect both the IMC and FXN results, another small suite of ∼300 simulations was run but using a different EOS model for each material. For SiO₂, rather than using an ANEOS model adapted to accommodate geologic materials in the high-temperature regime (Melosh 2007), the stock LEOS model was implemented. The alternate forsterite EOS was one developed for a separate Planetary Defense study of a related material, dunite, which has a similar composition. Alternate EOS files for ice and iron came from an older LEOS model created before data from NIST was added and from another recent study done by Sarah Stewart (Stewart 2020), respectively. The average percent differences between the original and alternate EOS IMC simulations by material were SiO₂: 12.7% ± 13.3%, forsterite: 9.7% ± 6.8%, ice: 0.6% ± 1.1%, and iron: −2.1% ± 10.3%, with an aggregate average of 5.2% ± 11.0%. The aggregate percent difference between the alternate EOS IMC and FXN simulations was −1.4% ± 6.7%. As the alternate EOS IMC versus FXN σ is nearly identical to the percent standard deviation of the suite of IMC versus FXN simulations it mirrored (porosity suite: σ = 6.6%), it can be inferred that FXN fitting weaknesses are the dominant source of the EOS IMC versus FXN discrepancies as well. For a best estimate, the total EOS uncertainty for this study was taken to be the total standard deviation of the percent difference in p_b between the original and alternate EOS IMC: 11.0%. It should be emphasized that this uncertainty is a best guess, and results could vary depending on what EOS is used. The original EOS tables for the materials in this study are included in the GitHub repository containing the python scripts (Burkey & Managan 2023) for reference.

The intrinsic discrepancies between the IMC and FXN simulation p_b values that result from imperfect profile and scaling fits were discussed extensively in the previous sections. Because the suite calls upon all three "layers" of the model and includes incident-angle scaling at low fluence, we recommend using the combined error from the angle × porosity simulation set: 8.2% for the aggregate fitting uncertainty. For cases where a more specialized estimation of error is required, modelers can choose from the eight individual material/source combinations in Table 3 under "Angle × Porosity."

Another discrepancy between reality and the model is the time-varying fluence rate of a proper 2D IMC simulation. Rather than having a constant fluence rate during the exposure time, a surface zone in a 2D IMC simulation would see an increasing fluence rate, while the X-ray photons reached the edge closer to the detonation point. Depending on the angle of the zone and the source duration, the zone's fluence rate would then decrease as the NED radiation crosses the far edge of the zone and into the next. This time dependency was not taken into account when building the energy deposition model in 1D due to the overall effect being quite small compared to the complication of attempting to implement it. However, the error introduced from neglecting this effect is on the order of others listed and should be included. To quantify the time-dependent source illumination uncertainty, 20 random time-dependent rate functions delivering the same total fluence and with beginning and end points at 0 were generated for t_end = 20 and 80 ns. These rates were then incorporated into 1D simulations for each material with a 1 keV BB source and f = 10⁻¹ and 10⁻³ kt m⁻². The p_b standard deviations for each set of 20 simulations were collected to quantify the uncertainty. A conservative error estimate of 5% will be listed here, which represents the average σ from the t_end = 80 ns cases. If they choose, readers could subtract 1%–2% from this value if they require low-density materials (ice or something with high porosity, ρ < 2 g cm⁻²) or by choosing a short t_end (∼20 ns). In addition, this uncertainty would encompass the discrepancy introduced by using a realistic NED output, which would also be nonconstant.

Finally, the issue of long-time behavior discussed in the previous section also introduces some additional potential error. For low-fluence cases, where the blowoff momentum levels out quickly, this source of error could likely be ignored outright. The largest errors likely occur at high fluence, where the previously defined loose convergence requirement (<1%/10 μs) is roughly met at 1 ms. For this reason, we take the 3.5% value calculated in Section 4.2 as a rough but conservative estimate of the error contributed from behavior beyond 10 μs. Alternatively, "convergence" can be defined in many ways, and the percent rate of change estimations from Section 4.2 should be used accordingly if the previous criteria are insufficient. However, the percent change rate's convergence toward zero would likely offset any increase in error from extending the projected end time.

With the assumption that each of these sources of uncertainty is independent from the others, we calculate the combined energy deposition model uncertainty by squaring in quadrature. The final combined energy deposition model uncertainty is 15.8%. This uncertainty is sufficiently small to use the model with confidence while not significantly inflating the overall uncertainty of a mission simulation. The individual sources of error and total uncertainty can be found in Table 4.

This error estimation only represents the uncertainty associated with using the energy deposition model and should not be used as the full uncertainty of a deflection velocity resulting from a mitigation mission. The real uncertainty of any nuclear planetary defense mission would be dominated by the NEO's properties, many of which will be poorly constrained unless a reconnaissance mission is also attempted beforehand. This uncertainty calculation also excludes the effect from exchanging the source spectrum approximations deployed in the model's construction for a real NED output, which will be briefly discussed in Section 5.1.

A final argument in favor of this model's validity comes from comparing these results to previously published work. Section 5 of Dearborn et al. (2020) covers a set of simulations that were run in a third mesh-based radiation-hydrodynamics code not utilized in this work. Figure 11 of Dearborn et al. (2020) presents results that are most relevant to this work. This figure contains deflection velocities from a series of 2D simulations with a 560 m, ϕ = 62% SiO₂ asteroid using 1, 1.5, and 2 keV BB sources, a 1 Mt yield, and standoff distances at 100 and 500 m. The final deflection velocity values from the simulations were taken at 0.25 ms, substantially longer than 10 μs. The deflection velocity values for each of the relevant BB source points in Dearborn et al.'s (2020) Figure 11 were 1.7, 3.1, 3.8, and 7.6 cm s⁻¹ for the 1 keV/500 m, 1 keV/100 m, 2 keV/500 m, and 2 keV/100 m cases, respectively. Drawing from Table 5, all Δvs from this work are between 80% and 87% of the Dearborn values and still increasing. Therefore, based on these sparse data points, no egregious discrepancy can be inferred.

For further comparison, the scripts used to generate the data points for Dearborn et al. (2020) were procured and used to rerun another set of 48 simulations matching the suite from Section 5 while collecting p_b versus time data. Figure 14 represents the p_b values at 10 μs collected from the Dearborn suite and compared to the Kull IMC suite using the same methods from Figure 13.

Zoom In Zoom Out Reset image size

When compared, the percent discrepancy distribution between the Dearborn and Kull IMC suites yielded a standard deviation of 14.5%. It should be emphasized that this study is meant to compare both the Kull IMC results and the closely matched FXN results to another radiation-hydrodynamics code. Though Kull and Ares produce very similar results, there are many likely reasons why the difference between Kull and the Dearborn code is larger. Some potential culprits include utilizing slightly different material models and the lack of a thorough resolution study with the Dearborn code. With more effort, the gap between the two radiation-hydrodynamics codes would mostly likely be bridged with sufficient explanation. However, even the larger percent discrepancy is still much smaller than the uncertainties introduced to the system when all of an asteroid's unknown properties are brought into play.

Finally, Figure 11 of Dearborn et al. (2020) reveals another key detail worthy of examination. Alongside the Δvs for the BB source spectra with a 1 Mt device, the figure also depicts a sweep of simulations at various standoff distances using a "realistic" NED output, which falls between the 1 and 2 keV data points. The writers recommend utilizing this plot when determining which BB source spectrum is best suited for a particular scenario.