Exomoons as sources of white dwarf pollution

Polluted white dwarfs offer a unique way to study the bulk compositions of exoplanetary material, but it is not always clear if this material originates from comets, asteroids, moons, or planets. We combine N-body simulations with an analytical model to assess the prevalence of extrasolar moons as white dwarf (WD) polluters. Using a sample of observed polluted white dwarfs we find that the extrapolated parent body masses of the polluters are often more consistent with those of many solar system moons, rather than solar-like asteroids. We provide a framework for estimating the fraction of white dwarfs currently undergoing observable moon accretion based on results from simulated white dwarf planetary and moon systems. Focusing on a three-planet white dwarf system of Super-Earth to Neptune-mass bodies, we find that we could expect about one percent of such systems to be currently undergoing moon accretions as opposed to asteroid accretion.


INTRODUCTION
White dwarfs (WDs) are the end-states of medium mass stars (M * 8M ), and their high gravity leads to rapid sinking of any elements heavier than helium. This sinking occurs on the timescales of days to millions of years (Koester 2009;Blouin et al. 2018), and generally leaves the spectra of white dwarfs devoid of metal features. Nevertheless, over a thousand white dwarfs have been observed to be 'polluted' with heavy elements (Coutu et al. 2019), with estimates that up to half of all white dwarfs are polluted . The source of the pollution is believed to be the remains of rocky parent bodies, which survived the post-main sequence evolution of the host star. Surviving plan-ets can scatter these bodies onto highly eccentric orbits, such that the objects approach the white dwarf and are tidally disrupted and subsequently accreted by the star (e.g., Debes & Sigurdsson 2002;Jura 2003).
The majority of WD polluters appear rocky (Doyle et al. 2020;Swan et al. 2019) and a few appear icy (Hoskin et al. 2020;Farihi et al. 2013). Many accreted bodies are chondritic in composition, and based on this and their apparent masses, it has been often assumed that these parent bodies were small rocky bodies analogous to the asteroids, or Kuiper Belt objects (Xu et al. 2017), in the solar system. Here we attempt to quantify the fraction of polluting bodies that are exomoons rather than asteroids based on the parent body masses required for observed WDs and numerical simulations of the frequency of each population's accretion. We are motivated by the recent discovery of beryllium in a polluted white dwarf (Klein et al. 2021) and the interpretation that the observed large excess in Be relative to other rock-forming elements is a tell-tale indicator that the parent body accreted by the WD was an icy moon.
In this interpretation, the excess Be is the result of irradiation of the icy moon in the radiation belt of its host giant planet . Payne et al. (2016) and Payne et al. (2017) showed that moons can be liberated by close encounters between planets and that liberated moons could be accreted by a white dwarf. Here we examine this proposal in greater detail using a statistical analysis and N-body simulations of the liberation and accretion of moons.
The probability of observing a moon versus an asteroid depends on the frequency and duration of the respective accretion events relative to the detectable amount of pollution on a typical WD. Based on occurence rates of dust around A-type stars, debris belts are expected to be ubiquitous amongst polluted white dwarf progenitors (Melis 2016). Accordingly, assuming that any planetary systems with moons available for accretion would also have a debris belt available, the probability of observing moon versus asteroid accretion in a given WD system is P moon accretion P asteroid accretion = T moons T asteroids , where T is the fraction of time that the associated population provides observable pollution. Because we are interested in cumulative times, T is dependent on the accretion rate of each population (T = T (accretion rate)). Because asteroids are expected to accrete much more frequently than moons due to their sheer number, P moon accretion /P asteroid accretion will depend on whether the accumulation of successive asteroid accretions is sufficient to sustain high masses of WD pollution, compared to the single accretion events of relatively larger-mass moons.
In this paper we estimate the parameters of Equation 1 for a system of three super-Earth to Neptune mass planets. A number of planetary architectures are capable of becoming unstable and aiding to feed meterial to the white dwarf (e.g., Maldonado et al. 2022;Stephan et al. 2017;Veras & Gänsicke 2015), however, we focus on this particular system because existing studies of asteroid (debris) belts in this architecture provide a baseline against which to compare moon accretions.
Our paper is organized around Equation 1. In §2 we introduce an analytical model for masses of polluting elements in the WD convection zone, which we apply to all accretion events throughout this work. We show that this model implies that observed polluted WDs require masses much larger than one would expect for typical asteroid belt objects, emphasizing the need to understand moon accretions. In §3 we use observations of polluted WDs to put limits on the levels of pollution required for the accreter to be detectable. §4 finds T asteroids , the cu-mulative timescale of detectability of asteroid pollution, using extrapolated asteroid accretion frequencies from previous studies. We present the results from our own N-body simulations for moon accretions in §4.2, and find T moons . Finally, we summarize all quantities and discuss the implications of Equation 1 in §5.

ANALYTICAL MODEL FOR CONVECTION ZONE MASSES
To determine the duration and pollution levels associated with an accretion event, we make use of the model by Jura et al. (2009) for the buildup of accreted material in a white dwarf atmosphere. Throughout this work we will refer to the model as J09. The model describes the time-dependent mass of the polluter currently observable in the convection zone of a white dwarf as a function of polluting parent body mass, element settling times through the WD atmosphere, and the duration of the accretion event. The model assumes the accretion disk, and therefore accretion rate, decays exponentially as one might expect from dissipative forces that depend on mass. Under this assumption, the mass of element Z that is observed to be in the convection zone of the white dwarf at the time of observation t after the start of the accretion event, M CV (Z, t), is: (2) where M PB (Z) is the mass of element Z in the parent body, τ set (Z) is the e-folding settling time of element Z, and τ disk is the characteristic lifetime of the accretion disk.
The observable pollution mass M CV (Z, t) in J09 depends on the settling timescale, which in turn depends on the properties of the host star. Therefore, variations in white dwarf temperature and composition have significant impacts on the maximum accumulations of pollution in the WD atmosphere and the timescales during which pollution levels are sufficiently high to be observable. In Figure 1 we show the pollution masses for a representative element and parent body in a DA (top) and a DB (bottom) white dwarf in order to illustrate these differences. Defined spectroscopically, DAs have atmospheres dominated by hydrogen, while those of DB white dwarfs are helium-dominated. While white dwarf classifications extend well beyond these two categories, for simplicity we will restrict our discussion to these two broad categories, using the terms DA and DB to mean hydrogen and helium-dominated in what follows. The primary difference between the two types of WDs in the present context is that DAs generally have settling timescales for the heavy elements of days to thousands of years while the DBs have settling timescales of 10 5 to 10 6 yr for these elements. Additionally, because hot DAs will have minimal, if any, convection zones, we consider M CV to more generally represent the mass of observable pollution in the WD atmosphere for these cases.
The maximum heavy element mass in the convection zone for each WD type occurs where settling times are approached. In other words, the maximum is the amount of accreted material that can build up in the atmosphere before sinking exerts an influence on the abundances. In the example to follow, we assume a settling time of 10 2 yr for the DA and 10 6 yr for the DB. Both cases are assigned an accretion disk e-folding lifetime of 10 5 yr.
In Figure 1, we show the fraction of the parent body mass that is currently observable in the convection zone as a function of time since the start of accretion, for an exemplary heavy element. Both cases illustrate the three phases of pollution: while accretion is ongoing and before settling begins, the mass of the element builds up in the atmosphere (increasing phase). When the settling and accretion rates equalize (the blue point in Figure 1), the mass stays relatively constant (steady state). Once settling dominates as accretion wanes, the mass of the pollution decreases rapidly (decreasing phase). Note that, in this model, the majority of a single accretion event is in a decreasing phase.
Because settling in the DBs begins after the majority of the parent body has been accreted onto the white dwarf, DBs can exhibit much larger fractions of the polluting metals than DAs, for the same parent body mass. The DA steady state phase occurs earlier and therefore at relatively higher settling and accretion rates than the DB phase, and the two DA rates conspire to cause very small fractions of the parent body to be observable at any given time.
Estimates for τ disk generally range from 10 4 to 10 6 yr (e.g., Girven et al. 2012). In Figure 2 we show how varying the disk timescale changes the pollution curves for the theoretical DA and DB stars shown in Figure  1. In particular, note that the maximum of the DA pollution decreases much more rapidly with increasing disk lifetimes than in the case of DBs.
We employ J09 in two ways. First, we estimate the parent body masses responsible for observed pollution in a sample of 21 white dwarfs and compare these masses with the distribution of asteroid and moon masses in our own solar system ( §2.1). Secondly, we use the model in conjunction with the average masses of asteroids and moons to determine the timescales of observability for both populations ( §3).  Figure 1. Examples of the fraction of parent body masses of a particular element that are predicted to be in the convection zone of a WD as a function of time, according to the J09 model (Equation 2). We assume an accretion disk lifetime of 10 5 yr and settling times of 10 2 yr for the DA (top) and 10 6 yr for the DB (bottom). Note the three phases of accretion: increasing, steady state, and decreasing. While analogous phases occur both for the DA and DB, the timescales defining the boundaries of each phase are swapped due to the longer DB settling times.

Steady-State Parent Body Masses
In order to calculate parent body masses for polluted white dwarfs, observed pollution masses (M CV (Z)) and settling times were collected from the references in Table 1. The total masses of heavy elements in each white dwarfs are shown in Figure 3. For an assumed disk timescale we then solve Equation 2 for the parent body mass at a range of possible elapsed accretion times for each observed element in the white dwarf. This gives an expression for the mass of element Z in the parent DB, T set = 10 6 yrs disk = 10 4 yrs disk = 10 5 yrs disk = 5 × 10 5 yrs Figure 2. Fraction of parent body mass for a typical rockforming element in the convection zone of a WD as a function of time, for three different accretion disk lifetimes. The DA pollution curves decrease by nearly a factor of ten with increasing disk times, while the DB curves decrease less significantly. The peaks of both curves shift to later times as the disk accretion time increases, reflecting changes in the time where steady state between accretion and settling is achieved.
body for an assumed elapsed accretion time t elapse and an observed metal mass of M CV (Z): (3) Note that we change the time variable to t elapse to emphasize the difference in the meaning of time between Equations 2 and 3. Equation 2 solved for the variation in polluting element mass with time since accretion for a single parent body mass. Equation 3 instead takes an observed heavy metal mass and provides a range of parent body solutions that depend on the time at which one assumes the observation was taken. To obtain the total parent body mass solution, we sum over all observed elements at a given elapsed accretion time, such that M PB (t elapse ) = M PB (Z, t elapse ). Observed M CV, Z , g DA DB Figure 3. Total masses of heavy metals observed for the white dwarfs in Table 1, sorted by type. The DAs tend to have lower observed masses of polluting elements compared to the DBs.

NLTT
As an example, Figure 4 shows the application of Equation 3 to the DA G149-28, and DB WD 2207+121. The vertical lines show the settling times for each element associated with each WD and the assumed disk lifetimes of 10 5 yr. Note that the shape of the parent body solution is roughly the inverse of the pollution mass curve, reaching a minimum during the steady state phase of accretion, when the pollution mass is at a maximum. Equation 3 shows that the steady state point, dM PB (Z, t elapse )/dt elapse = 0, coinciding with the minimum estimate for the parent body mass, will occur at time: In practice, when summing over multiple elements as in Figure 4, each element would reach the steady state point at slightly different times, due to the variations in settling times. Nonetheless, as long as the range of settling times are well above or below the disk lifetime, solving for elemental abundances at the time corresponding to steady state will give a minimum estimate for the total parent body mass. Thus, for the remainder of this work one can think of the 'minimum parent body mass' to be analogous to the 'parent body solution assuming steady state.' Note that this does not necessarily mean we are assuming all white dwarfs are in steady state, but rather that any other phase of accretion would re- quire a more massive parent body than the steady state solution to explain the observed metal pollution. We calculate parent body masses for the white dwarfs in Table 1 using the effective temperature and log g values reported in the references to obtain the WD convection zone masses and elemental settling times from the Montreal White Dwarf Database (MWDD) . We then derive M CV (Z) values from the relative abundances reported in the references, using the MWDD settling times and white dwarf envelope mass fractions. For comparison, we also considered the models provided by Koester et al. (2020), which calculate settling times that include the lack of convection zones in hot, hydrogen dominated white dwarfs. We find that the settling timescales derived in the Koester models are generally comparable to those reported by the MWDD, and therefore do not significantly change the resulting parent body masses.
Additionally, while we use all stellar parameters from the MWDD instead of those reported in each reference, we find that in most cases the values are in agreement, to within a factor of a few.
Minimum parent body solutions for the observed white dwarfs are shown in Figure 5, for a range of disk timescales. Because of the marked dependence on disk lifetime ( Figure 2), we expect variations in assumed disk lifetime to change the minimum parent body estimates for DAs much more than for DBs. This expectation is realized, as shown in Figure 5.
We find that the majority of parent body estimates for DB WDs are between roughly 10 23 -10 24 g regardless of the disk timescale chosen. DA estimates are generally lower unless we increase the disk timescale to 10 7 yr. This is beyond current estimates for disk lifetimes given in the literature, but does provide a more satisfactory agreement between parent body masses for DA and DB white dwarfs. Despite that effect, we will adopt a disk timescale of 10 5 yr for the remainder of our analysis as that value is more generally accepted in the literature and results in lower parent body masses that we can take as minimum estimates. Figure 6 shows the elapsed accretion times corresponding to the minimum parent body solution for each WD in the sample, using the stellar parameters from  MWDD, and with a disk lifetime of 10 5 years. The upper and lower limits show the range of accretion times for which the parent body solution is within a factor of two of the minimum. In Figure 7 we show how the distribution of calculated parent body masses (assuming a disk lifetime of 10 5 yr) compares to the distributions of moon masses in DA DB disk = 10 4 yrs disk = 10 5 yrs disk = 10 6 yrs disk = 10 7 yrs Figure 5. Parent body masses calculated assuming the white dwarf is in steady state, with four different assumed accretion disk timescales. DA parent body solutions tend to vary more dramatically with disk timescale than the DBs. Most DB parent body mass solutions are close to ∼ 10 23 ,and while the DA solutions are somewhat lower, the two begin to converge as disk timescales increase. Note that the 10 7 yr timescale is longer than most estimates for disk timescales.
the solar system as well as to the approximate distribution of asteroid masses. Moon masses and radii are from the JPL Solar System Dynamics group 1 . We calculate the asteroid masses assuming the distribution of asteroid radii is dN ∝ r −3.5 dr (Dohnanyi 1969), and that all bodies have the same density of 3 g/cm 3 . In reality, the majority of asteroids have lower densities (∼ 2.5g/cm 3 ), so we can consider the asteroid masses as upper limits. We assume that the range of asteroid diameters is 1-1000 km, corresponding with a range of masses of approximately 10 15 − 10 24 g. The WD pollution parent body masses are generally far larger than those defined by the asteroid distribution. While the DA parent body masses and the majority of the DB masses reside in the range of the largest bodies in our asteroid belt, such as Ceres (∼ 10 24 g) or Vesta (∼ 10 23 g), they are well above the bulk of the asteroid mass distribution. From our sample, for a disk timescale of 10 5 yr, we consider a mean DA parent body to be ∼ 10 21 g and a mean DB to be ∼ 10 23 g. In our solar system, there are about 30 and 15 moons that fall above these DA and DB masses, respectively. The large calculated parent body masses compared to minor solar system bodies implicates an observational bias. This may be partially due to the large masses required to detect pollution, which we outline in the next section.

POLLUTION DETECTION LIMITS
In §2.1, we showed that in the context of the J09 model, most observed WD pollution requires parent body masses consistent with the more massive moons of the solar system, and the extremely rare most massive asteroids. We now turn to observations to determine a lower limit of observable pollution in a WD atmosphere to place further constraints on the differences between moon and asteroid pollution.
Lower limits for observable masses of pollution in white dwarf atmospheres are obtained from measurements of calcium masses in polluted WDs. We choose calcium because there is a large sample of observations for Ca available in the literature with which we can assess minimum masses. We collect Ca masses for DA WDs from the SPY Survey (Koester et al. 2005) and . Fraction of bodies in each population with a given mass. Densities of 2.5 and 2 g/cm 3 are assumed to calculate masses for asteroids and moons, respectively. The asteroid curve is derived from Dohnanyi (1969) and the moon data are from the JPL Solar System Dynamics group. The parent body masses are split into DA and DB populations and assume a disk lifetime of 10 5 yr. The DBs tend to have larger parent body mass solutions, though the masses for both the DAs and DBs are situated towards the tail end of asteroid masses and mid-to-high moon masses. Note that the parent body masses shown are lower limit solutions.
masses for the DBs from Zuckerman et al. (2010). We first approximate a line to the lowest calcium masses in the SPY survey to derive an expression for minimum mass as a function of effective temperature (equation 5).
We then use the DB data to renormalize the line for the DB WDs. Because the data are expressed by number relative to hydrogen/helium, the relation between the minimum masses and temperature is dependent on the mass of the convection zone (mass of hydrogen for DAs or helium for DBs). Our expressions for the minimum masses of Ca that are observable in the polluted WDs are M min Ca, DA = 10 where m Ca is the mass of calcium, m H is the mass of hydrogen, m He is the mass of helium, and M CV is the mass of the WD convection zone or WD atmosphere. As the atmospheres of DAs tend to be relatively small (see Table 1 orders of magnitude lower. The lower limit for calcium detection in the DAs is consistent with their lower observed M CV,Z values, as described above. Taking the minimum observable calcium masses for each of the observed WDs, we now use the J09 model to calculate the parent body mass associated with each calcium mass limit. For each WD, we first apply the J09 model to the minimum detectable calcium mass to find the minimum mass of calcium in the parent body. We then calculate the total parent body mass by assuming chondritic composition (∼ 1% calcium by mass), an accretion disk e-folding time of 10 5 years, and settling times provided by the MWDD based on the effective temperature of each WD. Figure 8 shows the resulting parent body masses associated with the extrapolated minimum observed calcium mass for each WD in our sample. These bodies would provide just enough calcium pollution to be detected with current technology.
While the instantaneous minimum observable pollution masses derived from equation 5 require much higher masses for the DBs, we find that DAs and DBs require similar overall parent body masses to produce observ-able pollution. This is due to the difference in fractions of parent body that can build up in each type of atmosphere (Figure 1). Assuming each parent body is accreting as a single event, the resulting minimum parent body masses required for observable pollution for both types of WDs are generally larger than the mean for solar system asteroids, and closer to the masses of solar system moons, suggesting an observational bias against "typical" asteroids.
In the following sections we will examine how moons and asteroids compare when we allow for continuous accretion of material. This approach is particularly necessary for asteroids, which are thought to reach accretion rates that require material from multiple objects to be present in the WD atmosphere at any one time. Going forward, we will adopt the instantaneous minimum observable mass limits of 5.3 × 10 16 and 1.3 × 10 20 g of total heavy elements for a typical DA (T eff = 10000 K) and DB (T eff = 14000 K), respectively. For single chondritic accretions, these limits correspond to minimum observable total parent body masses of 1.3 × 10 19 g for the DA and 1.7 × 10 20 g for the DB. Note that because DAs accumulate much smaller fractions of parent body in their atmospheres at steady state, the minimum observable parent body masses for DAs and DBs are similar, despite orders of magnitude differences in their instantaneous limits.

CONTINUOUS ACCRETION MODELS FOR ASTEROIDS AND MOONS
While describing the J09 model in §2, we considered the increasing, steady state, and decreasing phases for a single accreting body. However, it is possible that pollution in the WD atmosphere could be from multiple parent bodies. In particular, Mustill et al. (2018) find that asteroids can accrete onto WDs continuously for up to billions of years. To assess how much of the pollution from continuous accretion could fall within observable limits, we consider populations of potential polluters based on mass frequency distributions of asteroids and moons in the solar system as guides. We then use total accretion rates determined previously from N-body simulations to construct synthetic pollution curves by applying the J09 model (Equation 2) to each event.
We currently restrict the comparison of moon and asteroid accretion to a three-planet system of super-Earths and Neptunes, a system that has previously been shown to result in high rates of asteroid accretion by Mustill et al. (2018). Additionally, we focus on the first 200 Myr past WD formation, the time frame in which asteroid pollution levels are expected to be at their maximum.

Asteroid Accretion
Mustill's simulations show that in the first few million years after WD formation, a debris belt can reach a peak accretion rate of ∼ 10 −4 N AB /Myr, where N AB is the number of asteroid belt objects. Because the 200 Myr time span is relatively short compared to the full length of time considered in the Mustill simulations, we consider the accretion rate to be approximately constant in our calculations. We therefore use the accretion rate to constrain the total number of accretions that can occur within the 200 Myr period, and then pick the specific accretion event times randomly from a uniform distribution, such that each point in time is equally likely to be the start of an accretion event.
The masses of accreting bodies for each event are chosen at random, without replacement, from the distribution of masses representing the solar system asteroid belt. We assume a range of radii of 0.5 − 500 km and a total mass of the asteroid belt of 3 × 10 24 g. Dohnanyi (1969) found the radius distribution for a collisionally-generated debris belt to be dN ∝ r −3.5 dr. Translating this distribution into masses, we have an asteroid mass range of 1.6 × 10 15 − 1.6 × 10 24 g and dN ∝ ρ 5/6 m −11/6 dm. The total mass of the belt is then M belt ∝ mdN , which we constrain to be the mass of the asteroid belt. Assuming that all bodies are spherical, we obtain M belt = A ρ 4π 3 r 3 r −3.5 dr, where A is a constant, and ρ is the density of the asteroids. Assuming a constant density of 3 g/cm 3 and our adopted radius range of 0.5 − 500 km, we find dN = 1.7 × 10 19 r −3.5 dr. This gives a total of about 12 million objects in the debris belt, for an accretion rate of approximately 1200 events/Myr. We therefore will only accrete a fraction of the total number of bodies in the full 200 Myr time period.
To sample this distribution, we split the range of asteroid masses into 12 bins, spaced logarithmically in mass, a choice which leaves one object in the Ceres-mass bin. We assign each object in the belt a mass bin according to the distribution described previously, again assuming a constant asteroid density of 3 g/cm 3 . For each accretion event, we pick an object at random, identify its mass bin, and select a mass at random from the range of masses associated with that bin. We obtain masses of individual elements in the parent body by assuming chondritic composition, and then evolve these masses through the J09 model (Equation 2) to track the total masses of polluting elements in the WD atmosphere. The number of bodies in the selected bin is then decreased by one.
The upper panels of Figure 9 show the results of this calculation assuming a disk e-folding time of 10 5 yr and settling times for a 10000 K hydrogen-dominated WD and 14000 K helium-dominated WD. Each peak in the figure corresponds with an accretion event, and is followed by a tail during which mass sinks out of the atmosphere. The colored curves show the mass of individual elements comprising the polluting debris in the convection zone as a function of time and the black dashed curve shows the total mass of heavy elements in the convection zone. Note that because asteroid accretions are very frequent compared to the settling timescales, material from at least one body can be found in the mixing layer for the majority of the 200 Myr time interval.
In §3 we found that the minimum convection zone pollution mass that is observable is about 5.3 × 10 16 g for a typical DA WD and 1.3 × 10 20 g for a typical DB, as derived from observed calcium masses. The horizontal lines in Figure 9 show the detectability threshold, and any peaks in pollution that exceed these limits are considered detectable periods of accretion.
In both the DA and DB cases, a long-term mass of pollution is sustained in the atmosphere, with peaks when more massive asteroids accrete. The sustained background pollution is generally not enough to exceed detection limits, but the more massive asteroid accretions are observable. This suggests that while multiple debris belt objects may be contributing to observed pollution, it is likely that the majority of the observed mass is due to a single, more massive body. We find that asteroid accretions onto the DA can exceed the pollution detection limit for a total of ∼ 29 Myr, resulting in a cumulative fraction of observable time of T asteroids,DA = 0.145 All else equal, pollution levels are higher for the DB case than the DA case as more mass can build up in the DB WD atmosphere due to slow settling times. Overall, in the DB simulation, the pollution is at an observable level for a total of ∼ 72 Myr out of the 200 Myr interval, such that T asteroids,DB = 0.360.
We note that in the asteroid accretion scenarios, the maximum masses of heavy elements that accumulate in the WD atmosphere are ∼ 5 × 10 16 g for the DA and ∼ 10 20 g for the DB. From Figure 3, we see that the WDs in the observed sample have total masses of heavy metals that can exceed the maxima reached by our continuous asteroid accretion simulation by factors of up to ∼ 10 3 . In addition to the expected vicissitudes of extrasolar asteroid belt masses, we consider that the higher observed masses could be due accretion of moons.

Moon Accretion Simulations
We carried out the same calculations for moons around WDs as we did for the asteroids in order to compare the expected detectability of accretion events for the two sources. For this purpose we require an accretion rate for moons. This rate comes from the efficacy  Mustill et al. (2018) and assuming the distribution of masses follows that of the solar system asteroid belt. Bottom row: Accretion of a single moon of mass ∼ 10 20 g onto a typical DA (left) and a moon of mass ∼ 10 21 g onto a typical DB (right). These masses represent the median mass of the subset of solar system moons that satisfy the criteria for being liberated from their host planets and providing observable levels of metal pollution. of liberating moons from host planets, and the accretion rate of the liberated moons onto the WD.
We simulate the separation of moons from their host planets during the stellar mass loss event that produces the WDs in order to estimate f moons , the frequency of moon accretions by WDs. Because moon orbital periods require very short time steps for integration, we break down f moons into two parts to accommodate computational limits. We define f moons as N moons × f accrete , where N moons is the number of moons in the system that can be liberated from their host planets and are able to provide detectable levels of pollution on a WD and f accrete is how frequently a moon from this population reaches the white dwarf.
We use the N-body code REBOUND (Rein & Liu 2012) to model three planets each with two moons. The IAS15 adaptive time-step integrator (Rein & Spiegel 2015) allows time-steps to be shortened or lengthened according to the occurrence of close encounters between particles. Due to computational limits, we set the smallest allowable timestep to be 0.1 days. Following the approach described in Payne et al. (2017) for moon liberations, we start each simulation with a three-planet system (no moons) and integrate the planet orbits during stellar mass loss. The mass loss excites the planetary orbits, eventually leading to orbit crossings. We halt this portion of the simulation at the first orbit crossing, when the periapsis of any planet falls below the apoapsis of the adjacent interior planet, or, alternatively, when the apoapsis exceeds the periapsis of the adjacent exterior planet. At this point we insert two test particle moons around each planet (6 moons in total). The simulation is then resumed, including stellar mass loss if it is still ongoing.
During the moon and planet portion of the simulation, we consider a moon to be accreted by the white dwarf if it passes within the Roche limit of the white dwarf (∼ 0.005 AU). However, we note that it is possible that bodies farther away than 0.005 AU could still be accreted by the white dwarf, for example through Alfven wave drag (Zhang et al. 2021). Averaging the frequency of moon accretions across all simulation trials gives f accrete .

Initial Conditions for Moon Simulations
While Payne et al. use planetary architectures as simulated by Veras & Gänsicke (2015) and Veras et al. (2016) to assess moon accretion, we carry out our tests using the same three-planet system simulated by Mustill et al. (2018) ( §4) in order to compare our results for moon accretions to those of debris belt accretions.
Each of our simulations begins with three planets around a 3M main sequence host star that evolves into a 0.75M white dwarf. The planets have masses of 1.3, 30.6, and 7.8 M ⊕ and initial semi-major axes of 10, 11.6, and 13.07 AU, respectively, as prescribed by the Mustill simulations. We focus on this particular set of planets as they resulted in the largest fractions of asteroids engulfed by the white dwarf in the Mustill study. We choose random initial inclinations in the range [0 • , 1 • ], initial eccentricities of zero, and random values between 0 • and 360 • for all other orbital angles.
Stellar mass loss is incorporated by updating the stellar mass according to the analytical formulae for singlestar evolution by Hurley et al. (2000). This code calculates stellar properties over 1 Gyr. For simplicity each of our simulations begins at the start of stellar mass loss such that white dwarf formation occurs ∼ 100 Myr after the start of the simulation.
For the three-planet systems, the first orbit crossing usually occurs before the end of the stellar mass loss event. Because mass loss speeds up significantly towards the end of the 100 Myr interval, the stellar mass is usually still close to ∼ 3M , and the semi-major axes of the planets have generally not increased dramatically, when the simulation is paused for moon insertion.
We follow the prescription for moon insertion described by Payne et al. (2017). The semi-major axis of each moon relative to the planet is chosen randomly in the range r min < a moon /R H < 0.4, where R H is the instantaneous Hill radius of the planet at the moment of the first orbit crossing. We considered two r min values of 0.004 and 0.04. This range results in numerically manageable, stable orbits. Payne et al. found that once moons are liberated from their host planet, the initial conditions of the moons cease to matter due to the intense scattering each moon experiences. The initial inclination of each moon is randomly chosen to be within 1 • of the plane of planetary orbits, the eccentricity is set to zero, and all other angles are randomly chosen.
For comparison with our initial conditions, Figure  10 shows the distributions of semi-major axes relative to the host planet Hill radii for solar system moons. Data are gathered from the JPL Solar System Dynamics group, and we assume a uniform density of 2 g/cm 3 for all bodies. The top panel shows the distribution of semi-major axes of moons by planet, and the lower panel shows the masses for the solar system moons. The range in a moon /R H in our simulations is seen to coincide with the upper third of the values exhibited by solar system moons; semi-major axes of a/R H = [0.04, 0.4] includes most of Jupiter's and Saturn's moons, but not the majority of Uranus' or Neptune's moons (Figure 10).  Figure 11 shows an example of the results from one simulation. The lower plot shows the instantaneous periapse of each moon, while the upper figures show snapshots of the orbital configurations. In this simulation, the closest approach of a moon to the WD is ∼ 0.007 AU. Table 2 lists the closest approaches of each of the moon test particles across all simulations, as well as the initial conditions for each moon's orbit, relative to its host planet.

Accretion Rates of Liberated Moons
Of the 60 moons comprising 10 simulations, one enters the Roche limit ( 0.005 AU) of the white dwarf within the first 200 Myr after white dwarf formation. We therefore set the accretion rate of liberated moons at 1/60 per 200 Myr after the formation of the host WD. Note that like the asteroid N-body results, the moon accretion rates are dependent on the number of bodies available in the system. Our result should be taken as the fraction of liberated moons that accrete.
In our simulations, all moons are liberated from their host planets. Consistent with Payne et al. (2017), as a Figure 11. Top: Snapshots of orbits from simulation, for planets (thin black lines) and moons (colors). The times indicate the time since the WD formed, and the WD is shown as a black star at 0 AU. Bottom: Periapse versus time for the moons, in the same simulation as show in the top plot. Each curve corresponds with the orbits of the same color in the top plot. Two moons are scattered onto hyperbolic orbits following a close approach to the WD.
conservative estimate, we assume that any objects outside of 0.04R H of their planet will be liberated. We therefore take the population of solar system moons with semi-major axes in this range as the population of potential moon parent bodies that could pollute a WD. In our solar system, out of a total of about 200 moons, 145 are situated at more than 0.04 R H , including the irregular satellites.
Because moon accretions are single events, we can use the J09 accretion model to set a lower limit on the mass of a moon that can provide the minimum observable mass. For the DA of 10000 K, the least massive observable parent body is 1.3 × 10 19 g and for the 14000 K DB the limit is 1.7 × 10 20 g. For the population of solar system moons exterior to our limit for liberation of 0.04R H , this gives a total of 22 moons that can provide observable pollution on a DA and 10 that can do so on a DB.
We now return to our expression for the frequency of moon accretions, f moons as N moons × f accrete . We found that 1/60 of moons liberated from their host planets are expected to find their way to the WD. Inserting the accretion rate from the N-body simulations, and considering the number of moons that can both be accreted and observed on the surface of the WD based on the solar system population of moons, gives

J09 Model for Moons
Applying the accretion rates of f moons,DA = 0.0019Myr −1 and f moons,DB = 0.0009Myr −1 results in 0.38 and 0.17 accretions in the first 200 Myr past WD formation, for the DA and DB WDs, respectively. One concludes that moons are expected to be visible as single accretion events, with no build up of pollution from multiple objects, unlike the case for asteroids. We use solar system moons as a guide for computing the median mass expected for accretion events with moons as parent bodies. In the previous section, we showed that only moons with total masses greater than 1.3 × 10 19 g and 1.7 × 10 20 g can be detected on the surface of a DA and DB, respectively. In order to obtain the population of observable moons, we therefore select the solar system moons that are above these mass limits, and are situated outside of our assumed liberation limit of 0.04R H . We assume the resulting moons represent the population of bodies that could both be liberated from their host planets and provide a detectable amount of pollution if they were to be accreted by their host WD. Recall that we found that 22 solar system moons meet these requirements for the DA accretion events and 10 Table 2. Minimum periapse (q) reached by each moon, in AU, and the parameters used to initialize the moon orbits. amoon is given relative to the Hill radius of the planet. Two moons were inserted around each planet. The WD Roche limit is at 0.005 AU. solar system moons for the DB accretion events. Assuming any of these moons would be equally likely to accrete, we used the median mass moon of each population to represent the median mass moon expected to be observed accreting onto each WD type. This results in a median mass for an observed accreted moon on a DA of 1.2 × 10 20 g and 3.7 × 10 21 g for a DB.
In the lower panels of Figure 9, we show the pollution curve due to accretions of these median mass objects on the same time axis as the asteroid accretions to illustrate how the single events compare with the continuous accretions. We assume that each moon has chondritic composition. Because these are single events, they follow the same phases of accretion as shown in Figure 1, however we now also show the variation in the abundances of individual elements. The heavy metal limit outlined in §3 is shown as a horizontal line in each plot.
For a single moon accretion, the DA case exceeds the mass limit for 0.57Myr while the DB is observable for 3.89Myr. Note that these timescales far exceed the duration of observability for a single asteroid on either WD type. However, because asteroid pollution accumulates from multiple bodies, the continuous asteroid calculations result in greater cumulative timescales of observability.
For a mean fraction of observable time for the moons we multiply these numbers by the ex- We can now use our results to evaluate the relative probabilities of detecting asteroids from a debris belt and moons in polluted WDs.
Returning to Equation 1, we now fill in the values derived for the DA and DB cases to find the relative probability of observing accretion of moons and aster-oids, yielding P moon accretion P asteroid accretion DA = 0.001 0.145 = 0.007 (6) and P moon accretion P asteroid accretion DB = 0.003 0.360 = 0.008. (7) Therefore, for three-planet Super-Earth/Neptune systems with both moons and asteroids available for accretion, we would only expect up to 1% of polluted DAs and DBs to currently have observable amounts of pollution due to moon accretion.

DISCUSSION
Over 1000 white dwarfs have observations of at least one polluting element in the atmosphere (Coutu et al. 2019). Around 20 are considered strongly polluted, with multiple rock-forming elements detected. Therefore, ∼ 2% of all polluted WDs can be considered 'highly polluted'. The large, moon-like mass solutions calculated in §2 represent the parent bodies associated with this highly polluted sample. Assuming that all of these high-mass parent bodies are indeed moons, this suggests that ∼ 2% of polluted white dwarfs are accreting moons. Of course, as shown by the pollution masses in Figure  9, accretion of the most massive asteroids may also result in high masses of pollution, so the population of the most highly polluted white dwarfs may also include the most massive debris belt members. Nonetheless, this statistic is consistent with our results derived from Nbody accretion rates, that ∼ 1% of overall pollution is expected to come from moons as opposed to less massive debris belt objects.
Uncertainties in this study include the timescales assumed in the J09 model, observational constraints be-yond what has been considered in our minimum detectable mass method, exoplanet moon and asteroid populations, and the effects of planetary architectures beyond the three-planet system considered as our test case. We now explore these various possibilities and their effects on the parent body solutions or numerical accretion models.

Disk e-folding timescale
The e-folding lifetime of the debris disks around the WDs, τ disk , determines how quickly pollution accretes onto the WD and, in conjunction with the settling times, limits the maximum mass of any given element that can accumulate in the WD atmosphere ( Figure 2). Throughout the asteroid and moon comparison in this paper we assume τ disk = 10 5 yr, as an accommodation for estimates that span from 10 4 to 10 6 yr.
In §2 we derived t min , the assumed elapsed accretion time which recovers a minimum solution for the parent body mass, or equivalently the steady state point in the J09 model. Plugging in t min to the J09 model, we can therefore write the general solution for the minimum parent body mass solution relative to the observed mass in the atmosphere as a function of the ratio of the disk and settling timescales, for an arbitrary element: where M PB (Z, t min ) is the parent body solution assuming steady state for the element Z and M CV (Z) is the observed mass of the heavy metal. Figure 12 shows Equation 8 applied to a range of disk-to-settling time ratios.
Writing the parent body mass expression in terms of the ratio of characteristic timescales shows why DAs and DBs have different sensitivities to changes in τ disk (outlined in §2) by illustrating the two limits of the parent body mass calculation. If disk timescales far exceed settling times, as they would for DAs, the minimum parent body solution will increase rapidly. For τ disk = 10 5 yr and a DA settling time of days, parent body solutions can reach factors of 10 7 times the observed metal mass. On the other hand, if settling times exceed disk timescales, as they may for DBs, τ disk /τ set approaches 1, so that the minimum parent body mass is equal to current mass of metal in the WD atmosphere.
Increasing the assumed disk timescale would shift parent body mass solutions to the right in Figure 12. While DB WDs with very long settling times would not be strongly affected (parent body solutions would still be roughly equivalent to the mass of metal in the atmosphere), minimum parent body solutions trend approximately linearly with τ disk /τ set when the disk timescale exceeds the settling time (DAs). Decreasing the assumed disk timescale would similarly not strongly affect the DBs but would decrease the parent body solutions of the DAs proportionally. Note that WDs with settling times within the range of disk timescale estimates will have non-linear dependencies on disk timescale. These effects can be seen in the parent body solutions calculated for the observed WDs ( Figure 5).
We now return to our parameters T moons and T asteroids , the timescales during which pollution exceeds detectable levels. For simplicity, we start by considering a single accretion with a generic observability timescale T . For a single event, T depends on the peak mass of pollution (M CV ) that can build up in the atmosphere, the duration for which a high mass of pollution can be sustained, and the detection limit associated with the WD.
The maximum mass of pollution deposited by a given parent body mass is the inverse of Equation 8. Therefore, pollution accumulation for DAs varies inversely with changing disk timescales while peak masses of pollution in DBs remain roughly constant. The amount of time that relatively large masses of pollution can be sustained in the WD atmosphere can be approximated as the difference between τ disk and τ set . As seen in Figure 4, these timescales are on either side of t min , where the maximum M CV is reached. For simplicity, we will consider DA settling times well below, and DB settling times well above, the range of possible disk timescales, so that |τ disk − τ set | will be approximately equal to the longest of the two timescales.
Based on this approximation, for DBs |τ disk − τ set | ∼ τ set , and to a reasonable approximation, changing τ disk should not affect T for the DBs. Therefore, and we anticipate T moons and T asteroids to be robust against disk timescales for DBs of sufficiently long settling times.
DA settling times are short, so |τ disk −τ set | ∼ τ disk , and increasing the disk timescale will lengthen the amount of time that peak pollution levels can be sustained. However, increasing τ disk decreases the maximum mass of pollution that can be accumulated. If the decreased pollution masses still exceed the detection limit, then T would increase, but if the new pollution masses do not exceed detection limits, T would fall to zero. The overall effect on T moons and T asteroids will therefore depend on the detection limits and distribution of debris belt masses considered. For our distributions of moons and asteroids, generally only the most massive bodies are contributing towards observable pollution, so assuming that most pollution would remain above the detection threshold, we would expect T moons and T asteroids to increase with increases in disk timescale.

Asteroid belt mass and size distribution
The accretion rate extrapolated from N-body simulations depends on the number of objects in the asteroid belt, which in turn depends on the total mass and assumed distribution of radii for the population of asteroids. Furthermore, the number of bodies contributing to heavy element pollution at any given time, and therefore the median mass of heavy elements in the WD convection zone, varies directly with the accretion rate. In this work, we assumed a solar system mass asteroid belt, with radii spanning 0.5 − 500km following a distribution of dN ∝ r −3.5 dr.
Due to the continuous nature of asteroid accretion, pollution remains at a relatively stable minimum for the duration of accretion (∼ 10 16 g total in the DA, ∼ 10 20 g for the DB), with short spikes to greater, observable masses when a particularly massive asteroid accretes. Whether asteroid accretion is observable for long time periods is therefore an 'all or nothing' issue, and is very sensitive to how a typical amount of accumulated heavy elements compares to the detection limit. If the limit is just above the typical mass that can build up from multiple accretions, we will observe only the peaks of the most massive asteroid accretions. However, if the detectability limit is just below the typical mass of accumulated metals, accretion is observable for the entire time period.
Given that many polluted white dwarf progenitors are estimated to be more massive than the Sun (Coutu et al. 2019), it is reasonable to assume that many of these systems may have had debris belts more massive than the asteroid belt. If the debris in these more massive belts follows the same collisionally-produced power-law distribution as described for the asteroid belt, there would be correspondingly more objects of any given mass. Assuming the accretion rates per number of available bodies is unchanged, a larger debris belt would increase the total number of accretion events, and therefore increase the typical pollution mass at every point in time, resulting in a larger T asteroids . If the moon accretions remain unchanged, this increase in T asteroids would decrease the fraction of pollution that would be due to moons.

The impact of planet spacing on moon liberation
Because moon liberations depend on close encounters between planets, and planet separations determine how quickly systems can become unpacked during stellar evolution, we expect the frequency of moon liberations to vary with planetary system architectures. The threeplanet system used throughout §4 has spacings of 5-7 mutual Hill radii, with the innermost planet situated at 10 AU. This arrangement of planets is somewhat more tightly packed than most observed systems. Separations for systems detected by the Kepler satellite generally peak around 14-20 mutual Hill radii (Pu & Wu 2015;Weiss et al. 2018), however these observed planets are all on orbits interior to ∼ 2 AU, and it is unclear if this trend would directly apply to outer planets.
One constraint on outer planet spacings is HR 8799 (Marois et al. 2008(Marois et al. , 2010, which hosts four giant planets (> 5M Jup ) exterior to 10 AU, and is a likely candidate for a future polluted white dwarf system (Veras & Hinkley 2021). Considering the most stable configuration for these planets (Goździewski & Migaszewski 2020), separations are approximately 2-3 mutual R H . From the sample of directly-imaged exoplanets, Nielsen et al. (2019) found that approximately 9% of stars more massive than 1.5M could host such massive planets outside of 10 AU.
In Figure 13 we show the orbital crossings that result in simulations with initial planet spacings of 5-10, 10-20 and 20-30 mutual Hill radii. Planet masses for all three simulations are 1.3, 30.6, and 7.8 M ⊕ , the same as used in §4. These simulations result in 502, 148, and 173 crossings for the 5-10, 10-20, and 20-30 cases, respectively. We find that while the number of crossings varies with planet spacings, crossings do still occur even at spacings more consistent with the majority of observed systems. From this simple comparison we do not 10-20 mutual R H 5-10 mutual R H 20-30 mutual R H Figure 13. The periapse evolution of three-planet systems for three different ranges of planet spacings in terms of mutual Hill radii: 5-10 (top), 10-20 (middle), and 20-30 (bottom). In each system, the innermost planet is at 10 AU and the spacings between each consecutive planet are randomly chosen from the stated ranges. Planet masses for all three simulations are 1.3, 30.6, and 7.8 M⊕. While the most closely packed system experiences the most orbital crossings between planets, crossings still occur in the more widely spaced systems.
anticipate that orbital crossings, and consequently moon liberations, would be entirely eliminated for more widely separated systems. Further study is necessary for more detailed connections between liberations and accretions and planetary architectures.

Populations of exomoons
In this paper, we assumed that moon populations would be similar to those in the solar system. It is possible that moon populations in exoplanet systems do not resemble the solar system. However, as there are no confirmed detections of rocky or icy exomoons, it is difficult to determine how many such moons a typical exoplanet system might have. Additionally, if exomoons in general are not found in all planetary systems, the probabilities of moon accretions derived in this work should be multiplied by the fraction of polluted white dwarf systems that do host exomoons.
A first-order limit on the number of moons in a given WD system should be the placement of moon-hosting planets. Dobos et al. (2021) found that whether exomoons can survive in a stable orbit around a given planet depends on the proximity of the planet to the host star. They conclude that planets on very short period orbits (< 100 days) are most limited in the fraction of moons they can retain, while planets with longer periods could retain at least 60% of their original moons. For the purposes of determining the populations of moons available to accrete onto a WD, the longer-period planets are likely more relevant, as the inner planets risk being engulfed by the star in its red giant phase.
In our three-planet system, planets were originally located 10-13 AU around a 3M star, with periods of 18−27 years. According to the Dobos et al. study, these planets could retain about 70% of their moons. The giant planets of the solar system have orbital periods of about 11 to 160 years, and have a similar moon retention fraction of 60 − 80%. This suggests that whether or not the distributions in mass and semi-major axis of moons in our theoretical three-planet system match those for our solar system's moons, systems resembling our threeplanet system should have a large portion of their moons intact and bound to planets when the white dwarf forms. Further constraints on exomoon populations could be made as Transit Timing Variation searches for exomoons progress (e.g., Teachey & Kipping 2021;Kipping 2021).

CONCLUSIONS
Motivated by the large parent body masses required to explain observed levels of white dwarf pollution, and the recent discovery of beryllium in a white dwarf atmosphere, we have used N-body simulations and the Jura et al. (2009) accretion model to assess the likelihood that a WD will be polluted by a moon. We focus this study on the first 200 Myr past WD formation for a planetary system containing three Super-Earth/Neptune-class planets. Extrapolating from asteroid N-body simulations, we find that such a planetary system could sustain an asteroid accretion rate of approximately 1200 objects per Myr. Assuming that the system had a population of moons similar to the regular solar system moons, we find from N-body simulations that we could expect up to about 0.4 moon accretions per 200 Myr.
Using the population of observed white dwarfs with calcium detections, we find that the pollution must have a calcium mass component of at least 5.8 × 10 14 g to be observable in a DA atmosphere, and 1.4 × 10 18 g to be detected in a DB. We use these limits to determine the cumulative fractions of time that moons and aster-oids can produce observable levels of pollution. Based on our numerical accretion model, we expect ∼ 1% of white dwarf pollution to come from moons as opposed to asteroids. If we consider, as a first-order approximation, that all of the most highly polluted WDs (requiring the most massive parent bodies) are polluted by moons, our parent body mass approach returns a similar statistic, of moons making up about 2% of polluters.