Devolatilization of extrasolar planetesimals by 60Fe and 26Al heating

Whilst the formation of Solar system planets is constrained by meteoritic evidence, the geophysical history of low-mass exoplanets is much less clear. The bulk composition and climate states of rocky exoplanets may vary significantly based on the composition and properties of the planetesimals they form from. An important factor influenced by planetesimal composition is water content, where the desiccation of accreting planetesimals impacts the final water content of the resultant planets. While the inner planets of the Solar system are comparatively water-poor, recent observational evidence from exoplanet bulk densities and planetary formation models suggest that rocky exoplanets engulfed by substantial layers of high-pressure ices or massive steam atmospheres could be widespread. Here we quantify variations in planetesimal desiccation due to potential fractionation of the two short-lived radioisotopes 26Al and 60Fe relevant for internal heating on planetary formation timescales. We focus on how order of magnitude variations in 60Fe can affect the water content of planetesimals, and how this may alter the formation of extrasolar ocean worlds. We find that heating by 26Al is the dominant cause of planetesimal heating in any Solar system analogue scenario, thus validating previous works focussing only on this radioisotope. However, 60Fe can become the primary heating source in the case of high levels of supernova enrichment in massive star-forming regions. These diverging scenarios can affect the formation pathways, bulk volatile budget, and climate diversity of low-mass exoplanets.

Substantial internal heating of planetesimals by SLRs, in particular 26 Al, in the Solar System was first theorised by Urey (1955).Evidence of a substantial 26 Al enrichment in the early Solar System in the form of an excess of the 26 Al decay product 26 Mg in calcium-aluminiumrich inclusions (CAIs) was discovered by Gray & Compston (1974) and Lee et al. (1976).With such levels of enrichment, 26 Al was shown to be a dominant heating mechanism for planetesimals and asteroidal ★ E-mail: j.w.eatson@sheffield.ac.uk † Royal Society Dorothy Hodgkin Fellow precursor bodies (Grimm & McSween 1989;Grimm & McSween 1993).The excess levels of 60 Fe have been less clear from meteoritic data (Tang & Dauphas 2012;Telus et al. 2018;Trappitsch et al. 2018;Cook et al. 2021), with estimates varying over two orders of magnitude (from 60 Fe/ 56 Fe = 10 −8 to 10 −6 ) in recent years, with the lower 60 Fe levels favoured in most recent measurements.At the lower end limit, the abundances of 60 Fe are not expected to be sufficient to cause substantial heating and volatile loss of planetesimals, however, substantial variations of the 26 Al/ 60 Fe ratio are expected across starforming regions (Vasileiadis et al. 2013;Boss 2019Boss , 2022;;Fatuzzo & Adams 2022;Parker et al. 2023).Multiple formation mechanisms exist for both 26 Al and 60 Fe (Lugaro et al. 2018;Diehl et al. 2021;Laird et al. 2023).While spallation from stellar cosmic rays has been suggested as a potential mechanism for enrichment of 26 Al (Adams 2021), the low levels of 10 Be (Mc-Keegan et al. 2000;Dunham et al. 2022) and mostly homogenous distribution of decay products throughout the Solar System (Desch et al. 2022) suggest that this mechanism is not the most influential (Parker 2020).Mechanisms external to a stellar system are more probable, such as through supernovae (Chevalier 2000;Adams 2010;Ouellette et al. 2010;Parker et al. 2014;Pfalzner et al. 2015), stellar winds (Gaidos et al. 2009), or AGB stars (Lugaro et al. 2018).Supernovae in particular produce both 60 Fe and 26 Al in abundance, with usually elevated 60 Fe/ 26 Al ratios relative to the Solar System (Kuffmeier et al. 2016;Lichtenberg et al. 2016), while stellar winds from massive evolved stars (such as Wolf-Rayets) can produce significant quantities of 26 Al over their lifetime (Limongi & Chieffi 2006, 2018).As the distribution of almost all SLRs are approximately homogenous, the half life times are short compared to the lifetime of the protoplanetary disks, and 26 Al-devoid CAIs (FUN-CAIs) are rare (Villeneuve et al. 2009;Desch et al. 2023), it is likely that SLRs are deposited immediately prior to or during the formation of the first CAIs.As such, we can infer that the enrichment of SLRs in a stellar system is sensitive to the local stellar environment immediately leading up to the formation of the stellar system (Adams 2010;Portegies Zwart 2019;Parker 2020).
Heating from SLRs has a significant impact on the volatile evolution of initially ice-rich planetesimals (Castillo-Rogez & Young 2017; Golabek & Jutzi 2021;Newcombe et al. 2023), which can have a knock-on effect on the properties of subsequent planets (Grimm & McSween 1993;Lichtenberg et al. 2021).Of importance for this paper, planetesimal internal heating from the decay of SLRs can result in melting, heating, vaporisation and subsequent outgassing of volatiles -particularly water and volatile carbon compounds.In the case of stellar systems with significant SLR enrichment, this can result in a significant reduction in surface and internal water content, which is a prime marker for bulk composition and affects the longterm climatic evolution of rocky planets.Constraining the influence of radiogenic heating for multiple radioisotopes on the volatile loss from planetesimals during accretion would allow for improvements in planetary formation models, better estimations of the ocean world population, and quantify the possible fractionation between major atmophile elements across planetary populations (Wang et al. 2019;Spaargaren et al. 2023).
Previous simulations of planetesimal devolatilization primarily focused on the influence of 26 Al on planetesimal water and carbon compounds (Lichtenberg et al. 2019a;Lichtenberg & Krĳt 2021;Lichtenberg & Clement 2022).This work, in contrast, specifically focusses on the influence of 60 Fe and thus varying 60 Fe/ 26 Al ratios on planetesimal desiccation.A parameter space exploration is performed to determine how 60 Fe heating can affect planetesimals of varying sizes, iron contents and degree of 60 Fe enrichment relative to the Solar System.Internal core-mantle differentiation and thus 60 Fe distribution is simulated by use of two distinct internal structure models, which are discussed in Section 2. As there are multiple parameters to be explored, the number of simulations was extremely large, as such, 2D models with a temperature-based rather than a chemistry-based model were utilised.Additionally, constraining some of these parameters for subsequent simulation sets was the basis for a number of simulation subsets, which are described in Section 3. Finally, the results are discussed and conclusions are made in Section 4 & 5.

Numerical & heating simulation
Numerical simulations of planetesimals were conducted using the I2ELVIS thermo-mechanical geodynamic code 1 , which uses a combination of the marker-in-cell and conservative finite difference schemes to simulate planetesimal evolution (Gerya & Yuen 2003;Gerya & Yuen 2007), with a similar setup as in Lichtenberg et al. (2019a) and Lichtenberg et al. (2021).Simulations were conducted with a fixed grid resolution of 250 × 250 cells.Higher resolutions 1 https://github.com/FormingWorlds/i2elvis_planetA comparison of H 2 O desiccation for a set of simulations with varying resolution.These simulations were conducted to determine whether simulation resolution impacted the results of the simulations.We found that there was no significant change as resolution varies, as such a value of 250 2 cells was utilised as a good compromise between speed and accuracy.
were tested, but found to produce similar results at the expense of significantly increased processing time (Fig. 1).Simulations evolve using a variable timestep, which is constrained by the following parameters: • maxxystep: The maximum displacement distance for material in a cell (/ • Δ).
• maxtkstep: The maximum temperature change per step.
Additionally, the maximum timestep is set to 2500 years throughout the simulation, in order to prevent numerical instability.
Upon initialisation I2ELVIS assigns properties to cells in order to make up the planetesimal and its surrounding environment.These properties are analogous to geological compositions within the planetesimal: • "Sticky air", the surrounding medium that ensures thermal stability, a proxy for the medium surrounding the planetesimal (Crameri et al. 2012).
• Hydrous silicates which contain water and the simulation's supply of 26 Al.
• Iron, which contains the simulation's supply of 60 Fe.
Hydrous silicates have three important states in a simulation that determine the water content of the planetesimal.Initially water in the wet silicates are assumed to be stored as ice, which can undergo phase changes based on temperature.Beyond a temperature of 273 K ( melt ) the ice melts, which delays further temperature increase as the enthalpy of fusion needs to be overcome.This phase change is reversible, and re-freezing can occur if the cell falls below the threshold temperature.Beyond a temperature of 1223 K, the upper limit of the amphibolite stability field, the entire water content of the cell becomes totally and irreversibly vaporised, and is assumed to leave the planetesimal through out-gassing processes (though these are not simulated); this temperature is referred to as  vap .Water retention across the planetesimal is stored as a fraction of the form: where Heating of the planetesimal is provided by the decay of 26 Al and 60 Fe isotopes.The specific heating of an SLR is given by the formula: where  A is Avogadro's constant 2 ,  SLR is the SLR decay energy,  is the decay constant of the SLR ( = ln(2)/ 1/2 ) and  SLR is the molar mass of the SLR.The specific heating rate for SLRs embedded within a rocky body can be calculated with the equation: where  E is the chrondritic elemental abundance and  SLR is the radioisotopic enrichment (Castillo-Rogez et al. 2009).Enrichment is measured relative to the early Solar System estimates of 60 Fe and 26 Al enrichment; the Solar System enrichment values for 26 Al and 60 Fe are detailed in Table 1.The total heating rate can then be calculated through the formula: where  26Al and  60Fe are the individual heating contributions per unit mass of 26 Al and 60 Fe.Fig. 2 shows decay heating from 26 Al and 60 Fe sources over time, as well as specific heating rates for the SLRs with Solar System abundances.As can be seen, 26 Al specific heating is significantly higher for the first 2 million years after CAI formation, however this is supplanted by 60 Fe after this point.In the case of Solar System abundances, however, 26 Al heating dominates due to the lower quantity of 60 Fe over the entire period of planetesimal formation and cooling.In other star systems that had interactions with supernovae during their formation we may find a greater abundance of 60 Fe, which would impact radiogenic heating and the water budget of subsequently forming planets.Whilst the specific heating rate of 60 Fe remains higher for longer, Solar System abundances result in a minimal impact on the total heating rate.The tables are produced using Eq. 2 and Eq. 3 using parameters detailed in Table 1 Throughout this paper the SLR enrichment is described relative to the estimated isotopic enrichment of the Solar System after CAI formation.This is parametrised as the SLR enrichment ratio, Λ.
where  SLR,★ is the post-CAI formation SLR enrichment of the system described in the simulation and  SLR,⊙ is the equivalent for the Solar System.Table 1 notes the enrichment values for each for Solar System estimates.Simulations are run for a total of 10 Myr, we assumed that the planetesimal being simulated is formed at 1 Myr after formation of the first CAIs.As such 1 Myr of radioactive decay is simulated prior to the start of hydrodynamical simulations to reflect this.

Assumptions & limitations
Our simulations are based on static planetesimals that do not accrete further, meaning that ongoing accretion of, for example, varying composition of pebbles that grow planetesimals is not taken into account (Lichtenberg et al. 2021;Johansen et al. 2023).The addition of accreting pebble layers on top of planetesimals has been shown to alter their thermal evolution (Sturtz et al. 2022).Furthermore, we do not model geochemical reactions operating in the host rock, such as the oxidation of iron metal by water flow.This is an important aspect, as the diversity of FeO and hydrogen-bearing phases in inner Solar System asteroids bear evidence of rapid water-rock reactions during planetary formation (Sarafian et al. 2017;Monteux et al. 2018;McCubbin & Barnes 2019;Grant et al. 2023), and may influence the final volatile fractionation on planetesimals.Following up from this work, we will thus investigate the more detailed multi-phase dynamics of water-rock aggregates (Travis & Schubert 2005;Bland & Travis 2017).This is a substantial simplification as mobilization and transport of volatile ices, fluids, and gases strongly influences their redistribution between forming metal core, planetary mantle, and outgassed reservoir (Suer et al. 2023).Additionally, the physical dislocation in multi-phase fluid treatments can differ substantially if the initial aggregate consists of a mixture that becomes comparable in ice and rock composition (Gerya 2019).Zhang ( 2023) recently suggested sublimation of volatile ice phases as an additional devolatilization mechanism for ice-rich planetesimals that cross to the inside of the water snow line (due to the snow line migrating outwards during the Class I stage of the solar protoplanetary disk, Lichtenberg et al. 2021), further highlighting the importance of multiphasic treatments of planetesimals.In a similar vein, ongoing metal-silicate differentiation is not treated directly, but bracketed by our assumptions on Fe distribution, which we discuss in the next subsection.

60 Fe distribution models
Shortly after the formation of the stellar system, subsequent rapid formation of planetesimals, and initial radogenic heating, the hot, undifferentiated liquidus phase material of the planetesimal begins to undergo differentiation (Nimmo & Kleine 2015).Segregation of silicates and iron occur on an approximately similar timescale, which occur over a period of 0.4 Myr to 10 Myr (Lichtenberg et al. 2019b), with the upper limit being the approximate timescale of our simulation.The segregation time of the planetesimal material is dependent on the initial formation time, degree of radiogenic heating (Dodds et al. 2021), and mode of core formation (Walte et al. 2020).As 26 Al heating drives melting of silicates, a strongly convective mantle forms, which can inhibit the growth of the core (Neumann et al. 2012).As the planetesimal cools this convective flow decreases and eventually the mantle solidifies, this occurs on a timescale of 10 Myr to 50 Myr (Neumann et al. 2018).There are two qualitative types of core formation in planetesimals, which differ in timescale: in percolative core formation, S-rich metals may form an interconnected vain network, which can remove mantle metals quickyl via gravitational drainage before the appearance of the first silicate melts (Yoshino et al. 2003;Ghanbarzadeh et al. 2017).However, laboratory experiments find that typically some metals remain stranded in the silicate matrix (Bagdassarov et al. 2009;Cerantola et al. 2015;Walte et al. 2023), which are then removed from the mantle once the silicates reach the rheological transition.As the range of core formation times across these two major modes intersects the epoch our simulations occur in, and a full model of core formation and material segrega- tion are beyond the scope of this work, two distinct planetesimal structures are simulated in this paper: (i) The core model: Iron is contained in a large metal core surrounded by silicates.This represents an idealised end-state for a pre-differentiated planetesimal, in a scenario where metal-silicate separation operates rapidly and early.
(ii) The grain model: Iron is randomly distributed throughout the planetesimal.This represents the end-member scenario of metal core formation completing slowly, after the main stage of internal devolatilization.
Fig. 3 illustrates the physical differences of these models, as well as their model-dependent parameters.An additional model with a mantle containing iron grains as well as a small, still forming core was considered, as it represented the mid-point between the two models.However, simulations with these models were not conducted as the results from the first two models were very similar, as discussed in Section 3.2.In addition to the parameters dictating iron abundance in each model, Ψ and Φ, the planetesimal radius and isotopic enrichment are varied, in order to explore a broad parameter space.

Core model
The first model considered for this research consists of an iron core surrounded by hydrous silicates.This model is comparatively easy to implement, and required no modification of the underlying I2ELVIS code.However, this model is not a realistic case for planetesimal formation while still retaining water, as iron core formation by liquid iron alloy percolation requires heating of a planetesimal above the dehydration temperature of 1223 K to 1700 K (Neumann et al. 2012).Whilst this case therefore does not represent and accurate account of how metal and volatile parts of a planetesimal would behave, it can approximate SLR repartitioning during core formation, which occurs in tandem with dehydration.As such, it represents an idealised endmember case, and can be compared with the successor grain model to determine how different distributions of iron and its related 60 Fe heating throughout the planetesimal body may affect the final water fraction.
The size of the core is controlled by the core radius ratio parameter, given by the equation: where  c is the core radius and  pl is the total planetesimal radius.
Throughout these simulations this parameter is varied from 0.0 up to 0.99; values beyond this were considered redundant as this is already an extreme value, and that there would be no hydrous silicate cells left to measure desiccation from.Whilst simple to execute, this model is less physically accurate than the grain model also described in this paper.Primarily, the degree of differentiation for even large planetesimals is comparatively slow, and occurs over a timescale of ∼ 1 Myr to ∼ 10 Myr after CAI formation.As such, a differentiated body is unlikely for a large planetesimal so early after its formation.However, with more extreme degrees of radiogenic heating from 60 Fe the process of iron melt segregation may be accelerated due to faster melting of the body.

Grain model
The grain model provides a more realistic simulation of the planetesimal shortly after formation.Upon initialisation a percentage of the cells containing hydrous silicates are converted at random to have an iron composition instead.The model therefore produces planetesimals where core formation has not occurred.A pseudorandom number generator is used to determine which cells are converted, as true random numbers were not deemed necessary.
As the percentage of cells affected by the change increases the amount of iron in the planetesimal increases, we define the Fe grain volume ratio, Φ, as the ratio of the number of cells changed at initialisation to the number of grains that were unaffected: where  Fe is the number of cells with the iron marker and  S is the number of cells with the hydrous silicate marker upon initialisation.This parameter is used as a stand-in for the iron content in the planetesimal; similarly to the core model, 60 Fe enrichment is a separate parameter.This model offers significant improvements over the iron core model in terms of physical accuracy, as it offers a closer analogue to the post-formation conditions and properties of a typical planetesimal.Additionally, planetesimals under these conditions should evenly heat faster than planetesimals with an iron core, due to the even distribution of 60 Fe throughout the rocky body of the planetesimal.However, whilst more accurate than the iron core model, the grain model still has some physical inaccuracies.In particular, iron and silicates cannot segregate due to limitations in the one-phase fluid approximation, though below the iron alloy melting temperature this is not a concern (Keller & Suckale 2019;Zhang et al. 2021).

Data recording & processing
2D outputs are performed every 50 timesteps, and contain the primitive variables (density, velocity & pressure), temperature and radiogenic heating rate.These outputs can then be used to calculate other values and can be averaged for specific regions of the planetesimal.Additionally, hydrous fraction, ice fraction and water fraction are dumped for the entire planetesimal at every time step.

RESULTS
For simulations where desiccation occurs we found that the evolution of the planetesimals was broadly the same.As the planetesimal begins to heat up over the first 10 5 yr of the simulation some cells  exceed the melting temperature, before subsequently exceeding the vaporisation temperature.Desiccation occurs rapidly at this point, before the planetesimal begins to thermalise, then cool, preventing any further desiccation from occurring.This progression can be seen in Fig. 4, where this radpid desiccation can be observed.
The parameter space exploration of this paper was broken into a series of simulations.Firstly, the core model was utilised in order to get a baseline of how 60 Fe heating affects water content, and its relative heating impact compared to 26 Al.This model was also used to constrain the ideal size of planetesimal for subsequent simulations, and also to determine how iron content can affect the efficiency of 26 Al as a heating element.Subsequently, the grain model was used for a separate set of simulations, in order to determine how iron distribution affects heating, in order to determine if heating is the only major variable in desiccation.

Core size ratio
First, we performed a set of simulations in order to explore how core size -and therefore iron abundance -affects desiccation of planetesimals.The core-to-radius ratio, Ψ, was varied between 0.05 and 0.95 in steps of 0.05, with additional simulations of 0.01 and 0.99 for completeness.Simulations for Ψ = 0 and Ψ = 1 were not conducted as these were determined to be redundant (Section 2.3.1).The common parameters for these simulations were a radius of 100 km with 60 Fe enrichment a factor of 10 3 greater than Solar System estimates.No 26 Al enrichment was included in these simulations, in order to focus solely on desiccation due to 60 Fe heating.The results of this set of simulations were then used to constrain the core size parameter space, in order to reduce the number of required simulations.
Fig. 5 shows that initially desiccation is fairly limited, but rapidly increases as the core-to-planetesimal ratio increases.Eventually a point is reached where increasing the core size does not lead to any appreciable increases in desiccation amount.Desiccation increases rapidly when the planetesimal is ∼ 1% iron by volume, plateaus at ∼ 10%, and decreases again at ≳ 50% iron by volume.As a result, two values of Ψ were chosen for further simulations: • Ψ = 0.25, which is in the range of the initial increase in desiccation, and corresponds to a Fe volume fraction of 1.56%.
, peak mean planetesimal temperature peak mean mantle temperature for simulations varying core size, Ψ, with 60 Fe enrichment 10 3 times greater than Solar System estimates.Desiccation becomes significant after the core radius,  c , is greater than 20% of the radius of the planetesimal,  p , quickly reaches a maxima and decreases rapidly beyond 90% due to only crustal deposits remaining.
• Ψ = 0.50, which is in the trough of maximum desiccation, and corresponds to a Fe volume fraction of 12.5%.

Radius comparison
After determining how desiccation amount varies through modification of the core size, we then progressed to determining how desiccation amount varies through the size of the planetesimal itself.Lichtenberg et al. (2019a) notes that there is a strong dependence on planetesimal radius and desiccation, with larger planetesimals undergoing greater desiccation due to the significantly greater mass of 26 Al.Planetesimals were varied from 1 km to 100 km with 1 dex spacing.The core-to-planetesimal radius ratio was maintained at Ψ = 0.5 and 60 Fe enrichment was varied between 1 ≤ Λ 60Fe ≤ 10 4 with 1/3 dex spacing between each simulation. 26Al enrichment was not included, in order to focus on the influence of 60 Fe on water content.
Fig. 6 shows that desiccation is strongly dependent on planetesimal radius.Small planetesimals < 10 km need an extremely high degree of 60 Fe enrichment in order to undergo significant desiccation, while above 10 km there is a rapid increase in desiccation amount, which tapers out by 100 km as water becomes increasingly rarefied.Even in the case of large planetesimals high 60 Fe enrichment is still required for desiccation; this is a recurring theme in later simulation subsets.Lines for the Λ 60Fe = 1 and Λ 60Fe = 10 simulation subsets were not included, as no desiccation was observed in any of these simulations.
In order to constrain the parameter space of our simulations, we use a common planetesimal sizes of 50 km and 100 km for all subsequent simulations, which is approximately at the peak of the birth size distribution of planetesimals generated by the streaming instability mechanism (Li et al. 2019;Simon et al. 2022).

Isotopic enrichment
Once planetesimal geometry had been constrained, we conducted an additional set of simulations in order to detail the influence of isotopic enrichment on desiccation.Fig. 7 shows the results of simulations where only 60 Fe enrichment occurs, Λ 60Fe is varied from 1 to 10 4 and Ψ is varied between 0.25 and 0.5.The results are consistent with the previous simulations, and shows that desiccation becomes significant above 200 times the lower limit for Solar System 60 Fe enrichment (recall that the lower limit is Λ 60Fe = 1.15 × 10 −8 ), with a larger core resulting in more desiccation.
To compare the influence of 60 Fe and 26 Al on water content, a number of simulations with both SLRs at varying degrees of enrichment were conducted.Similar to the previous simulations, two subsets were run, one where Ψ = 0.25 and another where Ψ = 0.50, 60 Fe abundance was varied from 1 to 10 3 times Solar System enrichment in a sequence of the form {10  }, planetesimal radius was not varied, and set as 50 km.Meanwhile, 26 Al enrichment was restricted to 0, 1, 1.78, 5.62 and 10 times Solar System enrichment, a sequence of the form {10 /4 }.
Fig. 8 is a contour plot of the final retained water fraction, H  , for the Ψ = 0.25 simulation set.These plots use a continuous gradient interpolated between these simulations using Delauney triangulation, which can be used to generate contours between sparse, non-linearly spaced data.We find that 26 Al is very effective at reducing the final water content of a planetesimal, with near total desiccation happening in planetesimals with only marginally greater than Solar System 26 Al abundance.Conversely, 60 Fe does have an impact, but only for simu- lation with little to no 26 Al.In the case of a simulation with total 26 Al depletion and Λ 60Fe = 1000 a final water fraction of H  = 0.796 is found, which is higher than the Solar System enrichment estimation simulation, which reports a value of H  = 0.509.Fig. 9 details the temperature change due to 26 Al and 60 Fe enrichment.An identical change is observed, with 60 Fe enrichment having a minimum impact on the peak mean planetesimal temperature, while 26 Al enrichment rapidly pushes the planetesimal temperature through the vaporisation point.
Fig. 10 is similar to Fig. 8 for the Ψ = 0.5 simulation set, a similar lack of influence due to 60 Fe is observed.However, low-26 Al-enrichment simulations do show markedly less desiccation, primarily due to having ∼ 30% of the hydrous silicate mass of the Ψ = 0.25 planetesimals.Conversely, as there is significantly more 60 Fe, we observe a final water abundance fraction for the  26Al = 0, Λ 60Fe = 1000 simulation of H  = 0.153, far below the maximum highest 60 Fe-only desiccation amount of the Ψ = 0.25 simulations.
To summarise, 26 Al is typically the primary radiogenic heating source for planetesimals with a clearly differentiated core, requiring 60 Fe enrichment far greater than the lower Solar System estimates to result in significant H 2 O outgassing.However, in planetary systems that are highly enriched in 60 Fe, for example by supernovae enrichment, 60 Fe can start to significantly contribute to planetesimal devolatilization.

Grain model
For the second set of simulations we implemented the grain model instead of the core model.Two major parameters are varied over these simulations: 60 Fe enrichment Λ 60Fe and iron volume fraction Φ. 26 Al enrichment Λ 26Al is also varied, though only between 0, 1 and 10.The parameter space of 26 Al enrichment was significantly compressed compared to the core model simulation set, in order to accommodate the expanded 60 Fe enrichment space, and introduce the Φ parameter without a drastic increase in simulation count.Λ 60Fe was varied from 1 to 10 4 -with 1 dex steps -while Φ was varied between 0.01 and 0.9.Planetesimal radius was increased from the core model simulation set to 100 km.The results of these simulations , for core model planetesimals where Ψ = 0.50, as there is a significantly lower mass of Al, there is a lower dependence on 26 Al enrichment and a greater 60 Fe enrichment dependence compared to Fig. 8.However, it is still not significant, as 60 Fe enrichment must be 2 orders of magnitude higher than Solar System estimates to have any reasonable impact.are very similar to the results of the core model, 60 Fe enrichment has to be significantly enriched compared to the Solar System in order for desiccation to occur in the case of 26 Al depletion.Meanwhile, 26 Al enrichment relative to canonical Solar System levels results in a greater deal of desiccation, and enrichment beyond that results in near-total desiccation.A divergence from the previous simulation set arises from varying Φ, where high values impede desiccation from 26 Al, as there is significantly less 26 Al available to heat the planetesimal.This is similar to the results observed for the core model Ψ = 0.5 subset as shown in Fig. 10.
Fig. 11 shows H  for varying values of enrichment parameters Λ 26Al , Λ 60Fe and grain volume fraction, Φ.Data is interpolated between simulations in the same manner as Section 3.1.3.Similarly to the core model simulations, 60 Fe enrichment has a significantly lower impact on final water content than 26 Al enrichment.In fact, the link between desiccation and 60 Fe enrichment is less pronounced than with the core model.In the case of fully 26 Al-depleted simulations, we find that no appreciable desiccation occurs until Λ 60Fe > 1000 even for simulations where Φ = 0.9.Fig. 12 shows the associated temperature of these simulations, here we can see similar dependencies as with Fig. 11, and that Tpeak is not exceeded outside of simulations with 26 Al enrichment.In summary, the grain model produces similar results to the core model, as such we can infer that bulk heating is more important than heating specific areas of the planetesimal.Furthermore, this lends further evidence that 26 Al enrichment is significantly more influential with 60 Fe/ 26 Al fractionation factors close to the Solar System, as 60 Fe enrichment must be multiple orders of magnitude higher to produce any significant effect, while a similar effect occurs with a single-digit factor increase in 26 Al enrichment.

60 Fe versus 26 Al
Our simulation results suggest that radiogenic heating from 26 Al is significantly more important compared to 60 Fe for 60 Fe/ 26 Al fractions close to the Solar System value.In certain cases, with 60 Fe enrichment more than 10 2 to 10 3 relative to the Solar System, 60 Fe heating can result in significant desiccation of a planetesimal, though these values may only be realised in massive star-forming regions, or with distinct supernovae enrichment in individual systems.Future considerations for determining which SLRs are important for the process of planetary formation would need to narrow down a typical galactic enrichment range, and whether such 60 Fe enrichment levels are common on a planetary system level.Additional simulations can narrow down this parameter space further by incorporating more comprehensive models of stellar feedback on local scales (Nicholson et al. 2019;Parker et al. 2023;Patel et al. 2023).

Other SLRs
Whilst other SLRs formed by stellar nucleosynthesis could be considered for simulation, such as 41 Ca and 53 Mn (Russell et al. 2001), these would have a lower abundance and enrichment to 60 Fe, resulting in an even lower influence on desiccation.As such, these were not included in our simulations.

60
Fe as a temperature sustainer?
Whilst our initial results show that there is a minimal impact of 60 Fe on planetesimal heating outside highly enriched systems, the longer half-life of 60 Fe could potentially sustain internal temperatures if the planetesimal was sufficiently enriched.This can be seen in Fig. 2 where the specific heating of 60 Fe is greater than the specific heating rate of 26 Al ∼ 2 Myr from CAI formation.In the case of Solar System abundances, however, this does not produce a significant amount of heating, and is still almost 2 orders of magnitude less than the specific heating rate of 26 Al at the specific heating crossing point.For more 60 Fe-enriched planetesimals this could result in sustained heating.
A continued heating source would prevent the planetesimal from cooling and continue the vaporisation and out-gassing processes.In order to infer the influence of 60 Fe over longer time scales from our simulations, we calculated the mean averaged radiogenic heating power density: over the planetesimal.Fig. 13 shows a comparison of  over time between simulations utilising the grain model where Φ = 0.25 with varying isotopic enrichment.After 2 Myr the average power density falls in the case of all simulations, becoming similar for simulations with matching 60 Fe enrichment amounts.It is clear that in all cases with significant 60 Fe enrichment that 60 Fe heating becomes the strongest heating mechanism near the end of the simulations, and should retain adequate heating for some time after the simulations have concluded.Whilst this enduring heating would not change the results with our desiccation model, more complex thermochemical models with out-gassing and without a discrete evaporation temperature would lead to greater desiccation rates.The sustenance of high temperatures is relevant for late-formed planetesimals and cometesimals in extrasolar systems.In the Solar System, there is evidence for prolonged planetesimal formation, essentially until the very end of the disk phase (Kleine et al. 2020;Simon et al. 2022;Lichtenberg et al. 2023).Recent evidence from JWST for water-enriched inner disks in both low-mass and high-mass star-forming regions (Perotti et al. 2023;Ramirez-Tannus et al. 2023) illustrates that desiccation of planetesimals until very late stages of the disk could impact the final composition of rocky and terrestrial exoplanets.How effective the sustenance by 60 Fe is will depend on the late-stage pebble flux and dislocation of vapour from ice-rock mixtures.Our simulations show that a very high 60 Fe enrichment relative to the solar system produces temperatures amenable for devolatilization, which would be sustained until the very late stages of disk evolution.During the late stages of the disk the varied composition of comet-like objects can impact the thermal evolution of still forming bodies within the disk (see Section 2.2) (Golabek & Jutzi 2021;Arakawa & Wakita 2023).This effect is not taken into account in the simulations discussed in this paper.As such, further investigation of the conditions where 60 Fe can still maintain devolatilization in planetary systems should be performed.Such investigations should ideally have a focus on smaller objects, high ice fractions, and varying composition (such as a varied C/O ratio) (Davidsson 2021;Lichtenberg & Krĳt 2021).

Possible influence on exoplanet populations
Gaining information on the SLR distribution across planet-forming systems has been a continuous challenge, which complicates assessing their influence on exoplanet formation (Lugaro et al. 2018;Parker 2020).However, it has been suggested that planetary debris in polluted white dwarf systems may serve as additional constraint on the SLR distribution across planetary systems if the effects of accretion energy and SLR heating can be distinguished (Jura et al. 2013;Jura & Young 2014).With increasing measurement precision and numerical modelling, research has shown that polluted white dwarfs indeed seem to conserve the footprint of planetesimal differentiation across a significant fraction of planetary systems (Bonsor et al. 2020(Bonsor et al. , 2023b;;Curry et al. 2022).
From a star formation perspective, for a typical Initial Mass Function (e.g.Maschberger 2013), massive stars that would produce 60 Fe are produced in significant numbers (i.e.>5) in regions that also form more than 1000 low-mass stars.Notable examples of such regions in our Galaxy are Cyg OB2 (e.g.Wright et al. 2015), Westerlund An alternative enrichment scenario, which avoids the destructive photoevaporation from massive stars, is from Asymptotic Giant Branch stars (e.g.Karakas & Lugaro 2016;Lugaro et al. 2018;Parker & Schoettler 2023).AGB stars are a stellar evolutionary phase that all 1 -8 M ⊙ stars undergo (Herwig 2005), and these stars produce significant yields of 60 Fe and 26 Al via their winds (Ventura et al. 2018).AGB stars were previously discounted as a viable source of SLR enrichment due to the supposed low probability of an evolved star encountering a young star as it forms its planetary system (Kastner & Myers 1994).However, Parker & Schoettler (2023) report the serendipitous Gaia DR3 discovery of an AGB star interloping through the young star-forming region NGC 2264, and furthermore, show that under reasonable assumptions for the initial conditions of the star-forming region, an AGB star could enrich the Solar System in both 60 Fe and 26 Al, with some stars attaining higher 60 Fe ratios than in the Solar system (e.g.Λ 60Fe ∼ 10 −4 ).
In addition, studies of planet-forming disks indicate that planetesimal formation is rapid, typically operating on a timescale of 10 5 yr (Drazkowska et al. 2022;Manara et al. 2022).Combined, this suggests that a significant fraction of exoplanetary systems are expected to undergo significant SLR-driven heating during the formation of their planets, with substantial fractionation possible between 60 Fe and 26 Al.
Exoplanet surveys have just started to reach into the low-mass regime, where it is possible to distinguish bulk-level under-and overdensities in the size regime of super-Earths and below (Wordsworth & Kreidberg 2022;Piette et al. 2023).Further investigation of these population trends will constrain the effects of accretion environment and geophysical evolution on bulk composition.While typically astronomical studies tend to focus on atmospheric escape and disk migration as the main effects on volatile bulk composition (Venturini et al. 2020a;Bean et al. 2021), we have shown here that variations in the main SLRs can drive substantial devolatilization of the planetes- The grain model is utilised where Φ = 0.25, while Λ 60Fe and Λ 26Al are varied.In all cases specific power drops significantly between 2 and 10 Myr after the start of the simulation as 26 Al decays, though simulations with a higher content of 60 Fe still have a significant power density.
imal building blocks.From a Solar System perspective, meteoritic evidence of initial incorporation of water, and other volatiles and rapid devolatilization of the building blocks of Earth and its planetesimal precursors is accumulating (Sarafian et al. 2017;McCubbin & Barnes 2019;Lichtenberg et al. 2021;Hirschmann et al. 2021;Grewal et al. 2021;Peterson et al. 2023;Lewis et al. 2022;Newcombe et al. 2023;Stephant et al. 2023;Grant et al. 2023).This presents a challenge for a wide spread of potentially Earth-like rocky planets with a similar geodynamic an climatic regime, as the mass fraction of surface water to allow for the coexistence of oceans and exposed land with a geodynamic regime as the Earth is limited to within a relatively narrow range (Cowan & Abbot 2014;Schaefer & Sasselov 2015;Noack et al. 2017).Observational insights on the water content and devolatilization trends of planetary debris around polluted white dwarfs (Farihi et al. 2013;Curry et al. 2022) and in debris disks (Lichtenberg & Krĳt 2021;Marino 2022;Bonsor et al. 2023a) will be crucial to narrow down the potential spread in these effects.
Ultimately, increased understanding of the influence of planetesimal degassing is required to understand the information received from exoplanet population analyses, planetary debris, and detailed characterization of individual low-mass exoplanets.

CONCLUSIONS
Overall we conclude that the influence of radiogenic heating on water content in protoplanetary systems due to 60 Fe is weaker compared to heating from the typically more abundant and energetic 26 Al.However, the contribution of 60 Fe to the devolatilization of earlyformed planetesimals becomes substantial for enrichment levels of 60 Fe/ 56 Fe ≳ 10 −6 , which is the upper end range of enrichment values that has previously been inferred for the Solar System.In exoplanet systems that formed in high-mass star-forming regions with substantial supernovae feedback, 60 Fe may thus be an important contributor to SLR-driven heating.Our parameter space explored 26 Al enrichment over an order of magnitude, while 60 Fe was explored over four orders of magnitude, we found that 60 Fe did not become dominant unless the system was 26 Al depleted, 60 Fe enriched, or possessed a high iron content.Each of these characteristics may be fullfilled in exoplanetary systems forming in diverse star-forming regions, however, we anticipate 26 Al to be the main driving factor of planetesimal differentiation and devolatilization.
Fractionation in 60 Fe and 26 Al as a function of the birth starforming region is expected due to intrinsic stochasticity of star formation (where variations in the IMF lead to significantly different numbers of massive stars, and therefore enrichment, e.g.Nicholson & Parker 2017), the stochasticity of dynamical evolution of young star-forming regions, which affects the amount of SLR ejecta that may be captured by an individual planetary system (Patel et al., submitted) and the short half-lives of the main SLRs.Our results thus suggest that SLR-driven internal desiccation of planetesimals may contribute to compositional scatter across exoplanetary systems, following the initial SLR enrichment of exoplanetary systems at their birth.Future work shall explore the geochemical effect of water flow and iron oxidation inside volatile-rich planetesimals, and the compositional effects of ongoing accretion on the final composition and atmospheric diversity of low-mass exoplanets.
A comparison of final water fraction over simulation cell count.The final result is insensitive to simulation resolution.

Figure 1 .
Figure1.A comparison of H 2 O desiccation for a set of simulations with varying resolution.These simulations were conducted to determine whether simulation resolution impacted the results of the simulations.We found that there was no significant change as resolution varies, as such a value of 250 2 cells was utilised as a good compromise between speed and accuracy.
Specific heating, based on Equation 3.

Figure 2 .
Figure2.A comparison of the specific heating,  SLR , and specific heating rates,  SLR , of the SLRs 26 Al and 60 Fe at Solar System abundances and enrichments.Whilst the specific heating rate of 60 Fe remains higher for longer, Solar System abundances result in a minimal impact on the total heating rate.The tables are produced using Eq. 2 and Eq. 3 using parameters detailed in Table1.

Figure 3 .
Figure 3.Comparison of the initial conditions of each model with an approximately equivalent iron fraction of  Fe / T = 0.25.The core model distributes all iron (yellow) into a central core surrounded by silicates (grey), with a planetesimal radius ratio of Ψ (Eq.6).The grain model randomly distributes iron throughout the planetesimal with a cell fraction of Φ (Eq.7).

Figure 4 .
Figure4.Core model comparison of water fraction, H, as a function of simulation time for planetesimals of radius 50 km and an un-enriched iron core with a radius ratio of 0.25.Simulations follow a similar structure, with most desiccation occurring over a very short period, before reaching a minimum as radiogenic heating reduces.

Figure 5 .
Figure 5. Core model simulation of final retained H 2 O fraction,

Figure 6 .Figure 7 .
Figure 6.Final retained H O fraction, H  for simulations with varying planetesimal radius,  pl .The core model is utilised with a core-to-planetesimal radius ratio of Ψ = 0.5.There is a strong inverse correlation between  pl and H  , similar to what was observed in Lichtenberg et al. (2019a) with water fraction, radius and 26 Al enrichment, though with significantly less pronounced desiccation.

Figure 8 .
Figure 8.Comparison of final H 2 O fraction, H  , for simulations with a core ratio, Ψ, of 0.25 in the core model.Desiccation increases dramatically as 26 Al enrichment increases; however 60 Fe enrichment only becomes significant with greater than approximately 3 × 10 2  60Fe,ss,low or higher.

Figure 9 .Figure 10 .
Figure9.Comparison of peak planetesimal mean temperature, Tpeak,pl , for simulations where Ψ = 0.25 in the core model.Similarly to Fig.8we see that there is a strong dependence on planetesimal temperature and 26 Al enrichment but an extremely weak dependence on 60 Fe enrichment and planetesimal temperature.This shows that in the case of the core model 60 Fe has a weak dependence.

Figure 13 .
Figure13.A comparison of mean power density, , as a function of simulation time for various simulations.The grain model is utilised where Φ = 0.25, while Λ 60Fe and Λ 26Al are varied.In all cases specific power drops significantly between 2 and 10 Myr after the start of the simulation as 26 Al decays, though simulations with a higher content of 60 Fe still have a significant power density.
is the volume of hydrous silicates that has never been heated above  vap ,  H 2 O, is the initial volume of hydrous silicates,  wet is the number of hydrous silicate cells still containing water in the liquid or ice phases, and  dry is the number of hydrous silicate cells that have been removed of water.This ratio is calculated and stored for every timestep.An important variable discussed is the final retained water fraction, H  , which is the value of H at the end of the simulation (H  =  H 2 O,  /  H 2 O, ).