Post-main sequence thermal evolution of planetesimals

White dwarfs that have accreted planetary materials provide a powerful tool to probe the interiors and formation of exoplanets. In particular, the high Fe/Si ratio of some white dwarf pollutants suggests that they are fragments of bodies that were heated enough to undergo large-scale melting and iron core formation. In the solar system, this phenomenon is associated with bodies that formed early and so had short-lived radionuclides to power their melting, and/or grew large. However, if the planetary bodies accreted by white dwarfs formed during the (pre)-main sequence lifetime of the host star, they will have potentially been exposed to a second era of heating during the star's giant branches. This work aims to quantify the effect of stellar irradiation during the giant branches on planetary bodies by coupling stellar evolution to thermal and orbital evolution of planetesimals. We find that large-scale melting, sufficient to form an iron core, can be induced by stellar irradiation, but only in close-in small bodies: planetesimals with radii $\lesssim$ 30 km originally within $\sim$ 2 AU orbiting a 1$-$3$\,M_{\odot}$ host star with solar metallicity. Most of the observed white dwarf pollutants are too massive to be explained by the accretion of these small planetesimals that are melted during the giant branches. Therefore, we conclude that those white dwarfs that have accreted large masses of materials with enhanced or reduced Fe/Si remain an indicator of planetesimal's differentiation shortly after formation, potentially linked to radiogenic heating.


INTRODUCTION
In an era of exoplanet detection, there is a growing interest in the interior dynamics and volatile content of exoplanets, two properties that are crucial for their habitability.However, our understanding of exoplanet interiors are limited by detection techniques which only reveal their bulk properties, and poorly constrained interior modeling (Dorn et al. 2015;Wang et al. 2018Wang et al. , 2022)).Fortunately, white dwarfs (WDs) that have accreted debris of tidally disrupted planetesimals (planetary building blocks) provide a potentially powerful tool to probe the interiors and formation processes of exoplanets.
White dwarfs, remnants of degenerate stellar cores, should originally preserve atmospheres predominantly composed of hydrogen/helium due to the rapid gravitational settling of heavier elements after radiative levitation becomes negligible (effective temperature ≲ 20000 K).Therefore, the metallic absorption features in the spectra of ∼ 20%-50% of the WDs possibly originate from recent/ongoing accretions (Zuckerman et al. 2003(Zuckerman et al. , 2010;;Koester et al. 2014;Wilson et al. 2019).
The observed WD pollutants reveal diverse compositions for planetary bodies (Putirka & Xu 2021).A large fraction of pollutants resemble the bulk Earth, consistent with tidal disruption and accretion of thermally processed rocky planetesimals (Harrison et al. 2018;Doyle et al. 2019;Harrison et al. 2021a;Trierweiler et al. 2023).
★ E-mail: yl817@cam.ac.ukExceptions with excess oxygen compared to that hosted in metal oxides indicate the likely accretion of water (Farihi et al. 2013;Hoskin et al. 2020).High Ca/Na, Ca/Mn relative to stellar abundances, corresponding to depletion of moderately volatiles is typically assumed to originate from the formation processes of the planetary system, for instance, incomplete condensation and devolatilisation during secondary melting (Lodders 2003;O'Neill & Palme 2008;Pringle et al. 2014;Siebert et al. 2018;Harrison et al. 2021a).A significant dispersion in Fe/Si is usually interpreted as a natural consequence of asynchronously accreted core/mantle-rich fragments from bodies differentiated due to decay of short-lived radioactive elements (e.g., 26 Al) (Jura & Young 2014;Bonsor et al. 2020;Buchan et al. 2022;Curry et al. 2022;Brouwers et al. 2023).
In the solar system, differentiation of planetesimals: large-scale melting and formation of global magma ocean, followed by gravitational segregation of iron from silicates and the formation of an iron core, is thought to be powered by the decay of short-lived radionuclides (e.g., 26 Al) and/or violent impacts (Keil 2000;Chambers 2004).This magma ocean phase, accompanied with further (moderately) volatile depletion distinct from that of incomplete condensation (Schaefer & Fegley 2008;Vollstaedt et al. 2020), shapes the interiors of terrestrial planets.
However, there is no evidence that these early-stage thermal processes inside the solar system planetary bodies are common in the exoplanetary systems.For materials accreted by white dwarfs, latestage thermal processes induced by stellar evolution, for instance,

Observed pollutants
Early AGB Thermally pulsing AGB Thermal history of white dwarf pollutants: formation from protoplanetary disks, main sequence and post-main sequence life of the star, scattering into the Roche limit, circumstellar disks formation, accretion.The color bars show the general luminosity and radius variations of the host star, and the thermally pulsing asymptotic giant branch stellar wind density.
heating due to stellar irradiation and tidal dissipation, may also result in (moderately) volatile depletion (Jura & Xu 2010;Malamud & Perets 2016, 2017a,b), large-scale melting and iron core formation, thereby mimicking their early-stage counterparts taking place around planetary formation.To probe the formation stages of a planetary system via WD pollutants, it is thus essential to distinguish between the early-stage and late-stage thermal processes.
This paper models the thermal evolution of planetesimals induced by stellar irradiation and quantifies in what locations and under what conditions planetesimals are sufficiently heated to undergo differentiation (large-scale melting and formation of an iron core) during the giant branches of the host star.As a result, we comment on the possibility that a high level of Fe/Si in the atmosphere of a white dwarf originates from planetary bodies differentiated under stellar irradiation instead of radiogenic heating/impacts.

Journey of planetesimals from formation to accretion onto the white dwarf
The thermal and orbital evolution of planetesimals accreted by the white dwarf is coupled to stellar evolution (Figure 1).After the formation of the planetary system, remaining planetesimals experience a long period of steady heating from the main-sequence host star, with its interior temperature approaching a constant value.Afterwards, the interior temperature of planetesimals rises significantly when the host star climbs towards the tip of the red giant branch (RGB, shell hydrogen fusion).Low-mass stars go through helium flash (runaway nuclear fusion) in their degenerate core, while intermediate-mass stars can ignite core helium burning quietly, with the former reaching much higher tip RGB luminosity and radius.Then, stellar luminosity drops, until rising up again near the end of core helium burning (CHB) phase, before entering asymptotic giant branch (AGB, shell helium fusion), when the (average) luminosity keeps rising.Towards the end of AGB, these stars undergo thermal pulsations (TPAGB, cyclical thin shell fusion), accompanied with a rapid increase in average luminosity (Hansen & Kawaler 1994;Cristallo et al. 2015).At this time, planetesimals experience a short period of intense heating, predominantly affecting a depth equal to the corresponding diffusion length scale.Close-in bodies also experience stronger drag force and/or tidal decay when interacting with the expanding (circumstellar) stellar envelope, spiralling inwards, potentially being disrupted and engulfed by the host star (Maloney & Gallagher 2010;Mustill & Villaver 2012;Villaver et al. 2014;Jia & Spruit 2018).
When the host star enters its WD phase after losing a significant fraction of its mass, gravitational perturbation from the planetary system becomes stronger.Planetesimals may be scattered onto highlyeccentric orbits, which in turn decay under tidal interactions with the orbital energy dissipating as heat inside these planetesimals.Finally, once scattered into the Roche limit, planetesimals become tidally disrupted, forming a circumstellar disk, circularized and acrreted onto the WD via a combination of the Poynting-Robertson (PR) effect, gas drag, and the Yarkovsky effect (Rafikov 2011b,a;Bochkarev & Rafikov 2011;Metzger et al. 2012;Veras et al. 2015;Malamud et al. 2020;Veras 2020;Veras et al. 2022).

Paper layout
This paper is organized as follows.In Section 2 we summarize the methodology, a synthesis of stellar (2.1), planetesimal's orbital (2.2) and thermal evolution (2.3) codes to obtain the temperature profile evolution of planetesimals, from which we quantify the degree of differentiation (2.4) in a given system.In Section 3, we present theoretical stellar evolutionary tracks (3.1), and the resultant temperature profile evolution of sample planetesimals (3.2).We also present the simulation results of parameter space study in Section 3.4: the maximum central temperature (3.4.1) and the degree of differentiation (3.4.2) throughout the thermal history of individual planetesimals, with various sizes and orbits.Then, we quantify the mass fraction of planetesimals in a given population that are differentiated as a consequence of stellar evolution until asymptotic giant branch in Section 3.4.3.In Section 4, we discuss the limitations, their possible alterations to the results (4.1), and the implications (4.2) for observations of white dwarf pollutants.In Section 5, we summarize this study.

METHOD
Figure 2 summarises our method.The thermal evolution of planetesimals crucially depends on their sizes (  ) and the stellar irradiation.The strength of this irradiation is quantified by the equilibrium temperature at the surface of the planetesimal, a function of stellar luminosity and distance of the body from the host star.First, we model the evolutionary track of the host star utilizing MIST (Section 2.1), from which we deduce the orbital evolution ( ()) and survival of a planetesimal for a given initial position (2.2).Second, based on  (), and the stellar luminosity variation ( * ()) from the stellar evolutionary track, we compute the equilibrium surface temperature of the planetesimal ( ( = ±  , )), serving as the boundary condition of the thermal evolution code (2.3.1).Third, we solve the time-depend temperature profile ( (, )) inside the planetesimal numerically (2.3).Finally, by choosing a critical point   such that any region of a planetesimal heated above which is assumed to melt, we quantify the degree of stellar-induced differentiation in planetesimal size (  )-initial position ( 0 ) space (3.4).

Stellar evolution
Evolution of the host star plays an important role in the thermal processing and survival of close-in rocky bodies.Stellar luminosity determines the irradiation on planetesimals, whilst stellar mass loss triggers orbital expansion of planetesimals, potentially weakening the irradiation.Additionally, the radius of the star is the key to the survival of planetesimals, as close-in bodies may be engulfed by the stellar envelope.
The set of evolutionary tracks MIST (MESA isochrones and stellar tracks, Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015;;Dotter 2016;Choi et al. 2016) are used to follow the stellar parameters with a range of initial conditions: initial mass of 1, 2 and 3  ⊙ with an initial metallicity of 0.014.MIST interpolates a grid of single stellar evolutionary tracks from the one dimensional MESA (Modules for Experiments in Stellar Astrophysics) code (Jermyn et al. 2023), without the need to specify overshooting parameters and fine-tuning of resolutions in order to resolve thermal pulsations of the thermally pulsating AGB (TPAGB) star.MESA solves stellar evolutionary equations utilizing implicit Newton-Raphson method, additionally accounting for timedependent convection in non-steady state.MESA includes a variety of equation of states, e.g., Skye for fully ionized matter (Jermyn et al. 2021), FreeEOS for partially ionized matter (Irwin 2012).In terms of opacity, MESA accounts for molecular opacity at low temperature region (e.g., close to the surface of TPAGB star), Compton opacity at extremely high temperature and conductive opacity in degenerate region.As a result, MESA is able to model complex post-main sequence evolution of stars, including thermal pulsations driven by nuclear burning in a thin shell.

Orbital evolution
The orbital evolution of a planetesimal, together with the stellar luminosity, determines the strength of irradiation on the planetesimal.As the star loses mass, planetary orbits expand.This orbital expansion conserves the specific angular momentum of the planetesimal and reduces the stellar irradiation received (2.2.1).With the expansion of stellar radius, close-in planetesimals approach and/or enter the stellar envelope.Survival and orbital evolution of these planetesimals depend on their interactions with the environment, via friction, ablation and disruption (2.2.2).We neglect orbital tidal decay (Appendix B2, Veras & Fuller 2019;Mustill & Villaver 2012) since it is usually negligible for planetesimal-sized bodies.

Stellar mass loss
The gravitational force on a close-in planetesimal is usually dominated by the host star during the majority of its (post) main sequence life.The direction of gravitational attraction is parallel to the planetesimal-star separation vector.As a result, there is no net torque on the planetesimal ( r × F = 0), and its specific angular momentum  is conserved: where  0 ,  0 and  * ,0 are the initial orbital semi-major axis, eccentricity and stellar mass, respectively, and (), () and  * () are the corresponding values at a later time .
For close-in rocky bodies, the mass loss timescale far exceeds the orbital timescale of planetesimals.In this case, mass loss is adiabatic (Appendix A), without inducing any eccentricity variation (Veras et al. 2011).Consequently, conservation of  can be simplified to: , which together provide the stellar irradiation, whose strength is quantified as the planetesimal's surface temperature ( (  = ±  ,  )).The thermal evolution of the planetesimal is then simulated to predict its internal temperature evolution ( ( ,  )), which reveals the fraction of the body that undergoes large-scale melting and differentiation.The roundrects and rectangles representing numerical simulations and outputs, respectively. is the radial coordinate inside the planetesimal,  is the time coordinate. (2)

Expanding envelope
During the thermally pulsing asymptotic giant branch (TPAGB), the stellar envelope of a low/intermediate-mass star expands to several AUs.Meanwhile, drastic mass loss of the TPAGB star leads to the formation of a circumstellar envelope.As a result, close-in planetesimals directly interact with stellar materials, spiralling inwards and losing masses.First, environmental gases/dust exert ram pressure on the moving planetesimal.If ram pressure exceeds the binding energy of the body per unit volume, disruption (break up of planetesimal) would occur.The disruption threshold can be approximated as (Jia & Spruit 2018): where v is the velocity of the planetesimal relative to the environmental materials, and  ext is the density of the environment (stellar/circumstellar envelope): where  sw is the radial speed of isotropic stellar wind driven by radiation pressure,  sw ∼ −  *   * ,  is the speed of light.The left hand side of Equation 3 represents ram pressure ( ram ) exerted on the planetesimal and the right hand side is the gravitational binding energy (  ) of the planetesimal per unit volume.For primitive (uniform composition) bodies with   ∝  3  , as   is independent of   and   scales as  2  , disruption of the planetesimal is possibly catastrophic.
Second, motion relative to environmental gases/dust leads to friction and loss of planetesimal's angular momentum, causing inward spiraling (Villaver & Livio 2009).The orbital speed of a planetesimal that may enter the stellar envelope during pulsations (∼ 10 km/s) far exceeds the planetesimal's escape velocity (∼ 100 m/s), indicating that the drag acting on the planetesimal is always in the hydrodynamic regime, and that gravitational focusing and accretion of environmental gases/dust onto the planetesimal is negligible.The hydrodynamic drag force can be expressed as (Staff et al. 2016;Jia & Spruit 2018;MacLeod et al. 2018;O'Connor et al. 2023): where   is the drag coefficient of the planetesimal, which is a function of Mach number and Reynolds number, and approaches a constant for large relative speed  (appropriate for our case) (O'Connor et al. 2023).We choose   = 1 and 0.5 to investigate the qualitative behaviour of the planetesimal's inward spiralling.We track the orbital evolution of the planetesimal throughout TPAGB by solving the equation of motion numerically with the additional friction force term in polar coordinates (, ), and consider that the body is lost once its disruption threshold (Equation 3) is met: If the planetesimal's equation of motion can be decoupled to radial and tangential components, assuming that the motion of the envelope materials during thermal pulsations is dominated by the radial component ( v = (  −  ext ) r +   φ), then: We additionally include the mass loss/size shrinking of the planetesimal in collisions with gases/dust (non-thermal ablation), which can be approximated as (Jia & Spruit 2018): We neglect thermal ablation (sublimation) as each thermal pulse only lasts for ∼ 100 yr, corresponding to a diffusion length scale of ≲ 100 m.

Thermal evolution
The key to assessing the regions where a thermal process occurs is the interior temperature profile of the planetesimal.This is tracked during the evolution of the star by the radial thermal diffusion equation (Narasimhan 1999;Jura & Xu 2010), where we assume that spherical symmetry is preserved.When the temperature at a nodal point reaches the critical point,  crit , the corresponding thermal process is assumed to occur instantaneously (see Appendix C for differentiation timescale).We make the simplification that thermal and physical properties of the planetesimal remain constant.Whilst a convenient simplification for the modeling it is also the case that heat only transports via conduction as convection in 1-D is suppressed by thermal inversion (strong irradiation on the surface).
Based on these assumptions, the radial thermal diffusion equation can be expressed as: where  is the radial coordinate, thermal diffusivity   =    is a function of thermal conductivity, , bulk density  and specific heat capacity   of the body.We choose similar thermal/physical properties as Ricard et al. 2009 andLichtenberg et al. 2021 (Table 1).For planetesimals with identical radius, the dominating factor distinguishing their thermal evolution is   .The sensitivity of our results to   is discussed in Section 4.1.3.
We solve Equation 9 numerically by implicit difference method to ensure the stability and accuracy of the long-timescale numerical integration (Gerya 2019): where  represents the number of radial steps and  represents the number of time steps.We assign 10 4 nodal points to a planetesimal's interior, and 10 5 time steps to a planetesimal's thermal evolution during each stellar evolutionary phase (main sequence, red giant branch, core helium burning and asymptotic giant branch).
In practice, thermal properties are pressure and temperature dependent.Furthermore, heating in the body can be asymmetric, leading to convection in 3D.These factors may be important in thermal processing of a planetesimal and we discuss further the corresponding consequences in Section 4.1.3.

Surface equilibrium temperature
The boundary condition of Equation 9is the surface equilibrium temperature of the planetesimal, a quantification for the strength of irradiation.The surface temperature acts as the driving force of heat penetration into the interior of the planetesimal.We simplify the problem by assuming circular orbits for all planetesimals, but will discuss the effect of eccentricity by introducing time-averaged equilibrium temperature over one orbital period in Section 4.1.4(Mé ndez & Rivera-Valentín 2017).The surface equilibrium temperature of a planetesimal on a circular orbit has the form: where  is the distance of the planetesimal from the host star, equivalent to semi-major axis  for circular orbits,  * is the stellar luminosity,   is the bond albedo of the planetesimal,  is infrared emissivity ( ≈ 1),  is Stefan-Boltzmann constant,   is the fraction of re-radiation (  ranges from 0.5 for a tidally locked body, to 1 for a fast rotator).For simplicity, we choose   = 1 to preserve spherical symmetry.

Differentiation fraction
Based on the temperature profile evolution, We quantify the degree of differentiation in each planetesimal by the fraction of volume (   ) that could have melted (had its temperature increase above  crit ) during the giant branches.We then repeat this process for a population of bodies with different sizes and initial semi-major axis.Noting that planetesimals are not uniformly distributed in planetesimal size (  )-initial semi-major axis ( 0 ) space, we further introduce the distribution of planetesimals   - 0 space, in order to quantify the degree of differentiation of the planetesimal population.
We consider that the white dwarfs accrete rocky planetesimals that started with semi-major axes between   (critical orbit within which planetesimals end up being engulfed by the host star) and 10 AU, with the total mass of planetesimals lying between  0 and  0 +  0 to be given by  total  0 ∝  − 0 , where  is chosen to be 1 2 based on Minimum Mass Solar Nebula model (MMSN) (Hayashi 1981;Crida 2009).
We further assume that observable white dwarf pollutants mainly come from rocky planetesimals with radii between 10 and 100 km, with the number of planetesimals between   and   +   to be given by    ∝  −   , where  is chosen to be 4 based on the distribution of minor objects in the Solar system (Ivezić et al. 2001;Schlichting et al. 2013).The size distribution can be transformed to We combine spatial and size distributions of planetesimals and express the joint distribution as: The choice of power law indexes above may not apply generally in other systems.We discuss the effect of these indexes in Section 4.1.4.
To assess the degree of differentiation in a system based on the joint distribution of planetesimals, we consider two scenarios.In the first,  Mm , the ratio of total mass of material whose temperature exceed  crit to the total mass of planetesimals accreted onto the white dwarf is considered, regardless of whether planetesimals fully differentiate: where   (  ,  0 ) is the volume fraction of differentiation for a planetesimal of radius   with initial semi-major axis  0 .
In the second,  Mp , only those planetesimals that melt > 95% of their volume (defined as fully differentiated/iron core formation) are considered, under the assumption that the processes leading to the identification of core or mantle-rich material in white dwarf atmospheres require the formation of an iron core in the parent body: where   [   ] is of the form: The applied integration limits to Equation 13 and 14 are summarized in Table 2.

Stellar evolutionary tracks
In this section we summarise and compare the MIST evolutionary tracks (Section 2.1) of two samples, I: 1  ⊙ low-mass star and II: 2  ⊙ intermediate-mass star (Figure 3), both of which go through AGB and end in C-O white dwarfs, while undertaking distinct evolutionary tracks near their RGB tips.The zoom-in plots correspond to thermally pulsing asymptotic giant branch (TPAGB) with the origin of the time axes marks the start of TPAGB phase.The age of the star at the start of of its TPAGB phase ( 0 ) is added.
Figure 3 and Table 3 illustrate that sample I is more enhanced in size and luminosity during its RGB tip (approaching half of its tip AGB values) compared to sample II, whose tip RGB luminosity (radius) is around 1% of its tip AGB value.This occurs because low-mass (∼ 0.8-2  ⊙ ) stars, unlike their intermediate-mass (∼ 2-8  ⊙ ) counterparts, lack gravitational pressure in their cores and start helium burning with helium flash (runaway fusion in the contracted degenerate cores), until thermal pressure dominates, boosting their sizes and luminosity significantly at the RGB tips.
Mass loss patterns for both samples are similar: rapid mass loss near the end of TPAGB.In this case, Equation 2 indicates that outward migration of planetesimals mainly occurs at the end of thermal pulsations, suppressing the growth in planetesimal surface temperature (  , Equation 11) with stellar luminosity ( * ).The   of a planetesimal initially orbiting sample II at 1.5 AU shows clear decreasing trend near the end of TPAGB, despite the generally increasing  * (Figure 3, lower panel).For the same reason, planetesimals are much more likely to escape from the envelope via stellar mass loss during late pulsations.

Temperature profile evolution
As is described in Section 2.3, the time evolution of planetesimal's interior temperature profile is calculated based on the stellar irradiation, quantified by   (Equation 11).We show the thermal evolution of 3 sample systems (listed in Table 4) in Figure 4.In all 3 samples, the temperature evolution in the interior of planetesimals lags behind that of its surface temperature (  ).The time delay of a layer's temperature relative to   can be estimated by the diffusion time scale ( diff =  2 diff   ) corresponding to the depth of this layer and is especially relevant for rapid   variations, e.g., during the tip RGB and TPAGB.For instance, thermal pulsations lasting ∼ 100-1000 yr correspond to a diffusion length scale  diff ∼ 100 m, thus only affecting a thin surface layer of this size.
The interior temperature of the planetesimals in sample IA and IIA rises to the corresponding equilibrium temperature within a small fraction of stellar main sequence life, since the diffusion timescale ( diff ) of 30 km is around 30 Myr, much shorter compared to the main-sequence lifetime of the host stars.The initial temperature profile of planetesimals are unimportant unless the host star is more massive than 3  ⊙ and   exceeds 100 km (at which point, the diffusion length scale corresponds to the main sequence lifetime of these shorter-lived stars satisfies:  diff ( main sequence ) ≲   ).This is supported by the uniform temperature profile of these 3 samples at the start of RGB (Figure 4, lower panel).
Most of the thermal processes in planetesimals of sample IA and IB occur at the tip RGB, while AGB heating in IA and IB only penetrates layers of depths ∼ 10 km because of its short duration.On the contrary in sample IIA, AGB heating is much more significant than its RGB counterpart.This is consistent with the results in Section 3.1: the enhanced RGB luminosity of low-mass stars due to helium flash increases their ability to heat planetesimals during the RGB stage, compared with their AGB stage, whereas intermediate-mass stars do not go through the He flash, and therefore their planetesimal heating is dominated by the AGB phase.
When exposed to intense stellar irradiation (as in the case during the giant branches), a planetesimal's interior temperature decreases inwards (thermal inversion).Thermal processing requiring strong irradiation, e.g., differentiation, starts at the exterior of planetesimals.This thermal inversion ceases at the start of CHB phase, when the stellar luminosity drops before rising up again.At this point, the interior temperature of a planetesimal decreases outwards at its outer layer, preserving longer if the body is larger (sample IB).

Envelope entry
In this section we track the orbital evolution of planetesimals that eventually enter the stellar envelope during its TPAGB, based on planetesimal's equation of motion (Section 2.2.2).Our simulation verifies the planetesimal's orbital evolution is still dominated by outward migration from stellar mass loss, rather than stellar wind drag before envelope entry (see Appendix B1).We present the orbital evolution and mass loss for 4 sample planetesimals of   = 100 km after entering the envelope of a 2  ⊙ star, with  0 = 1.3 and 1.4 AU, and   = 1 and 0.5, respectively, in Figure 5.
Planetesimals in all four systems considered spiral inwards to their disruption limits (Equation 3) within 50 yr.According to Equation 5, deceleration due to drag force scales as , indicating that smaller planetesimals experience a stronger deceleration from drag and spiral in faster.Therefore, planetesimals with   ≲ 100 km cannot survive in the stellar envelope when ignoring other effects (e.g., interactions with the planetary system).Meanwhile, much larger planetesimals/planets, for instance, with   ∼ 1000 km feels less drag effect (≲ 10% decay in  during the first envelope entry) and may have enough energy to (partially) unbind the envelope.
Planetesimal mass loss due to ablation is negligible in all 4 samples.Ablation mainly occurs when the planetesimal is close to its disruption limit, where both environmental density and planetesimal's speed relative to the stellar envelope increase rapidly.

Thermal processing inside planetesimals
In this section we present the degree of thermal processing/differentiation in planetesimal radius (  )-initial semi-major axis ( 0 ) space based on thermal equations described in Section 2.3.Section 3.4.1 shows the maximum temperature raised at the centre ( = 0) and half radius ( =   2 ) of the sample planetesimals in our parameter space.Section 3.4.2demonstrates the volume fraction of melting (   ) in this space for two temperature thresholds (  ), 1400 K and 1800 K. Based on   (  ,  0 ), we present, in Section 3.4.3, the estimated mass fraction of differentiated samples according to two definitions,  Mm in Equation 13and  Mp in 14.

Maximum central temperature
We present the maximum temperature at  = 0 and   2 in planetesimals of two sample systems, I: 1  ⊙ and II: 2  ⊙ host star, in Figure 6.This figure shows that any thermal process triggered by stellar evolution take place preferentially in smaller-sized bodies closer to the host star, with its   and  0 dependence determined by the corresponding temperature threshold.
As is shown in the left panels of Figure 6, all completely melted planetesimals (which satisfy  ( = 0) ≥   = 1800 K) later enter the stellar envelope during the thermally pulsing AGB (fall to the left of the red dotted line).Considering the uncertainties in modeling thermal pulsations of AGB stars, we smooth out the oscillations of stellar radius and compute the corresponding engulfment limit, represented by the blue dashed lines in Figure 6 (see Section 4.1.1).This new engulfment limit leaves a narrow triangle-like parameter space of completely melted planetesimals that survive the AGB phase of the host star (sample I: 0.8 AU ≲  0 ≲ 1 AU,   ≲ 28 km, II: 1 AU ≲  0 ≲ 1.3 AU,   ≲ 18 km).This result emphasizes the importance of thermal pulsations on the survival of the most thermally processed planetesimals.
On the other hand, if  crit is lowered to 1400 K (see Section 4.1.2for the reason of this choice), some completely melted planetesimals survive the host star's thermal pulsations (sample I: 1.3 AU ≲  0 ≲ 1.7 AU,   ≲ 28 km, II: 1.5 AU ≲  0 ≲ 2.2 AU,   ≲ 22 km).Meanwhile, there is an increasingly larger area under a lower temperature contour in Figure 6.In other words, a lower  crit corresponds to a much larger   - 0 space where the planetesimal is heated above this threshold over a given fraction of its volume.This  4).The period of each evolutionary phase is normalised by its overall timescale, with the absolute time along the horizontal axis of the plot.The lower panel shows the temperature profile evolution of planetesimals in sample system IA, IB and IIA.The contour plots start from the red giant branch (RGB) of the host star, followed by the core helium burning phase (CHB) and the asymptotic giant branch (AGB), with the age of the star increasing clock-wisely.Same normalisation as the upper panel is applied.
result stresses the strong dependence of the degree of differentiation on  crit .
The thermal processing in sample I is mainly limited by the strength of heating, whereas in sample II the thermal processing is mainly limited by the timescale of heating.Compared to sample I (upper panels of Figure 6) that undergoes helium flash during its RGB, a given temperature contour in sample II (lower panels, quiet ignition of core helium burning) intersects the axes at smaller   but larger  0 (for instance, see the 1800 K contour lines).Furthermore, the increase in   value of each contour line-the   axis intersections from  = 0 (left panels of Figure 6) to  =   2 (right panels) is larger for sample II (lower panels), indicating that compared to sample I, the temperature gradient in sample II is steeper.These features are consistent with the facts that 1. the maximum thermal processing in sample I occurs at the RGB phase, which is much longer (lasting around 1.5 Gyr, corresponding to  diff ∼ 100 km) with much smaller average luminosity gradient (1.7  ⊙ /Myr) compared to the AGB life of sample II (20 Myr,  diff ∼ 10 km, 272.3  ⊙ /Myr, when planetesimals undergo the maximum heating), and 2. the tip AGB of sample II is around five times as luminous as the tip RGB of sample I.

Volume fraction of melting
We present the fraction of volume that reaches sufficient temperatures for melting at any point during its evolution of each planetesimal in   - 0 space (   (  ,  0 )) for two sample systems, I: 1  ⊙ and II: 2  ⊙ host star considering two  crit , 1400 K (right panels) and 1800 K (left panels) in Figure 7.
The similarities between Figure 6 and Figure 7 are summarized below: • Smaller planetesimals closer to the star undergo larger-scale melting • All planetesimals where more than 10% by volume reaches temperatures above 1800 K are engulfed • Compared to 1800 K, the fraction of planetesimals where a portion of their volume reaches above 1400 K increases significantly With the decrease in   , the contours of Figure 7 converge to the critical value of  0 (black dashed vertical lines in the left panels), beyond which the surface temperature of planetesimals can never reach   .This convergence results from the decrease of the duration when   >   with the increase in  0 .

The degree of differentiation
Due to the fact that observations of white dwarf pollutants may occur at any point of a white dwarf's accretion history, such that the sizes and initial positions of the accreted bodies are unconstrained, it is useful to consider the average effect of stellar irradiation on the interiors of a population of planetesimals.We quantify this average effect by the mass fraction of the planetesimal population where large-scale melting is induced by the host star's giant branches.In order to calculate this fraction, we consider a nominal population of planetesimals, orbiting within 10 AU surviving the giant branches, with sizes between 10 and 100 km, as described in Section 2.4 with . The volume fraction of melting (   (  ,  0 )), as described in Section 3.4.2 is integrated over planetesimal size (  )-initial semi-major axis ( 0 ) space to get the mass fraction of melted regions (  Mm , Equation 13), and the mass fraction of planetesimals melted over 95% of the volume (  Mp with  = 95, Equation 14), relative to the total mass of the planetesimal population.In this section we present  Mm and  Mp at different  crit in 3 sample systems, I: 1  ⊙ and II: 2  ⊙ and III: 3  ⊙ host star in Figure 8.
Both  Mm and  Mp in all 3 samples increase rapidly with the decrease in the critical temperature of melting,  crit .There are 2 orders of magnitude difference from 800 K to 1800 K in   for all 3 samples. Mp always lie below the corresponding  Mm , converging rapidly to 0 when  crit approaches the maximum surface temperature of planetesimals within 100 K, indicating its much stricter requirement of small   and  0 (long-lasting intense heating).For  Mm (  Mp ) to exceed 10%,  crit requires lowering to ≲ 1300 K (1200 K) for sample I, ≲ 1400 K (1200 K) for sample II and ≲ 1400 K (1000 K) for sample III.
At  crit ≳ 1100 K,  Mm of identical  crit decreases with stellar mass.In this region, the rapid increase in stellar luminosity with stellar mass is the dominant factor affecting the degree of differentiation.The contrary occurs at  crit ≲ 900 K, where the duration of intense irradiation becomes more important.A similar feature presents in  Mp , at  crit ≳ 1700 K and ≲ 1200 K, stressing the growing importance of the timescale of intense heating for large-scale melting.

DISCUSSION
Our model predicts that stellar evolution may lead to differentiation in small planetesimals close to the host star.However, the size of the parameter space where large-scale melting occurs, as well as its contributions to the white dwarf pollutants will largely depend on how thermal pulsations occur and the critical temperature of melting: for differentiated planetesimals to occupy 10% of the mass of the population described in Section 3.4.3, crit must be lowered to ∼ 1300 K.
In this section we discuss the main limitations in this study: stellar/planetesimal orbital evolution during TPAGB (Section 4.1.1),thermal evolution of planetesimals (Section 4.1.2and 4.1.3),and planetesimal distributions (Section 4.1.4),together with their possible impacts on the results.Based on the white dwarf atmospheric model, we discuss the observability of accretion of planetesimals differentiated under stellar irradiation (Section 4.1.5).
By comparing our predictions to observational data of polluted white dwarfs, we further discuss the necessity of early thermal processes such as radiogenic heating, impacts and incomplete condensation (Section 4.2).

Thermal pulsing asymptotic giant branch
Thermal pulsations are accompanied with unstable cyclical nuclear burning at the base of He and H shells.No steady state is available because of the significantly different reaction rates of H and He burning.Powered by the intense helium flash, a convective zone forms in the He intershell, affecting stellar pulsations and dredge-ups (Hansen & Kawaler 1994;Cristallo et al. 2015).As a result, thermal pulsations are extremely sensitive to specific assumptions of convection, for instance, convective overshooting.In reality, convection is a 3-D process, which cannot be fully simulated by MESA, a 1-D stellar evolution code (Paxton 2021;Jermyn et al. 2023).
On the other hand, with the rapid mass loss of the host star during thermal pulsations, the gravitational attraction from the host star weakens, while the interactions of planetesimals with the planetary system becomes stronger.The orbital evolution of planetesimals may be dominated by resonances, and planetesimals may be scattered inwards/outwards.A massive planet can also unbind the loosely bound pulsating stellar envelope, altering the condition of engulfment and the stellar evolutionary track.
However, we have shown that the expansion of the stellar envelope during thermal pulsations, without any planetesimal-planet interaction, leads to engulfment of the most thermally processed planetesimals.For instance, Figure 7 shows that almost all planetesimals with melting covering > 15% of their volumes (   ≳ 15%) with  crit = 1800 K are engulfed during the TPAGB.
To account for the uncertainties in modelling thermal pulsations and the additional effect of the planetary system, we repeat the analysis in Section 3.4.3while smoothing out pulsations (similar to SSE code, Hurley et al. 2013) for two sample systems, I: 1  ⊙ and II: 2  ⊙ host star, shown in Figure 9.In both samples, there is ∼ 5%-10% increase in  Mm and  Mp .This increase is especially significant at high  crit , e.g., 1600 K-1800 K, where orders of magnitude difference are present.
Furthermore, the luminous AGB phase may lead to non-negligible Yarkovsky effect and YORP effect arising from asymmetric thermal emission of planetesimals, altering their orbital (semi-major axis and eccentricity) and spin parameters (spin angular velocity and obliquity).We simulate the maximum Yarkovsky effect and the maximum YORP effect acting on a planetesimal (see Appendix D).Our sim- ulations illustrate that for our smallest planetesimals (  = 10 km ), the outward/inward Yarkovsky drift can reach ∼ 0.4 AU at the end of TPAGB life of the system.However, compared to the uncertainty in the planetesimal's critical semi-major axis at the end of the TPAGB phase (∼ 1 AU) arising from modelling of thermal pulsations, and accounting for the fact that the maximum Yarkovsky drift is only reached under extreme circumstances (e.g., tidal locking, maximum/minimum entries of Yarkovsky matrix), the Yarkovsky effect may only play a minor role in the survival of close-in planetesimals.On the other hand, the maximum YORP spin up has the possibility to break up the smallest planetesimals whose angular velocity already exceeds ∼ 0.5  break (= √︃ Another uncertainty of TPAGB is the mass loss pattern.Although an adiabatic approximation is still valid in our parameter space, mass loss could be axial instead of isotropic, such that only a small proportion of stellar luminosity contributes to the momentum of ejecting masses, leading to much slower stellar wind and denser circumstellar envelope.When exposed to the dense axial stellar wind, ablation of planetesimal becomes faster, accompanied with faster inward spiraling.As a result, although the orbital evolution of planetesimals unaffected by the axial stellar wind remain dominated by stellar mass loss, the engulfment regions, and thus the critical orbits, extend outwards for planetesimals inside the overdense stellar wind.More close-in (the most thermally processed) planetesimals are engulfed during the TPAGB, reducing the degree of thermal processing of the whole planetesimal population ending up accreted by the white dwarf.
In this study we assume a constant metallicity, 0.014 for all sample Comparison between  Mm , mass fraction of melted regions in single planetesimals, and  Mp , mass fraction of planetesimals melted over 95% of the volume, relative to the total mass of the planetesimal population, in 1 (blue, upper panel) and 2 (red, lower panel) solar mass star systems, before and after smoothing out pulsations (green), with the same assumptions in Figure 8. host stars, which is unable to represent the whole population of white dwarfs.For stars undergoing identical evolutionary phases, a higher metallicity corresponds to larger opacity, and hence a less compact star with more energy dissipated in the interior, resulting in stronger thermal expansion during pulsations.Meanwhile, stars with considerably distinct metallicty may undergo distinct evolutionary phases.For instance, when the metallicity increases from 0.014 to ∼ 0.05 for a 1 solar mass star, the star no longer undergoes an AGB phase, corresponding to a scenario that more close-in planetesimals avoid engulfment with the degree of thermal processing (dominated by the RGB phase) in individual bodies remained.The degree of thermal processing averaging over all planetesimals accreted by the white dwarf hence increases.

Critical point of melting
In this work we consider a simple model in which any portion of a planetesimal heated above a critical temperature is considered to have melted.The key question is whether a single critical point is valid and if so, what the appropriate critical temperature would be.Planetesimals have mixed compositions and hence melting ranges, whose boundaries correspond to the highest temperature for solidus ( solidus , no melting) and the lowest temperature for liquidus ( liquidus complete melting), instead of a single melting point marking the critical transition from solid to liquid.However, we can define a critical point of melting, as the transition from solid-like to liquidlike behaviours, which occurs typically around a liquid fraction of 40% (Costa 2005;Caricchi et al. 2007).This critical liquid fraction may be reached at different temperatures (within the melting range) in different regions depending on the exact composition, local gravity, etc.However, given the other uncertainties in this model, using a single temperature provides a reasonable qualitative behaviour.
The uncertainty in the critical point ( crit ) lies in both the poorly constrained melting range ( solidus - liquidus ) and the position of  crit within this melting range: the degree of melting required for efficient melt migration.
According to the model in Johansen et al. 2023, the pressure dependent melting range boundaries,  solidus and  liquidus for planetary silicates can be estimated as: where  0 ∼ 1661 K,  0 ∼ 1.336 GPa,  ∼ 0.134 for  solidus and  0 ∼ 1982 K,  0 ∼ 6.594 GPa,  ∼ 0.186 for  liquidus .For the largest primitive planetesimal in our parameter space (  =100 km), there is only an increase of ≈ 3 (1) K in  solidus ( liquidus ) from the surface to the centre of the body because of pressure variation On the other hand, the melting range covers ∼ 300 K (1661 K-1982 K at 0 pressure), which contributes to the main uncertainty in   .Furthermore, Scheinberg et al. 2015 suggests a slightly different melting range of ∼ 1400 K-1800 K, and McCoy et al. 1999 concludes that   ∼ 1600 K with   ∼ 1700 K based on experiments.Based on these studies,  crit may vary from 1400 K to 2000 K by maximum, corresponding to around 10% (5%) difference in  Mm (  Mp ) as is shown in Figure 8.

Thermal penetration
Equally important to the critical temperature is the thermal diffusivity   , which has a strong influence on the ability of heat to penetrate into the planetesimals ( diff ∝ √   ), and which we approximate as a constant.In reality,  is a potentially strong function of both the composition of the planetesimals and local temperature.As planetesimals heat up,  likely decreases and thus, could decrease the ability of heat to penetrate during the AGB/RGB.A quick calculation is presented based on the model in Miao et al. 2014: where   has the unit of mm 2 /s,  ranges from ∼ 0.13-0.20 and  ranges from ∼ 210-330 for different types of rocks.There is a factor of ∼ 5 decrease in   from our initial condition,  = 150 K, to  solidus ∼ 1400 K, equivalent to a factor of 2 decrease in the penetration depth during AGB/RGB.This means that the radius of the largest planetesimal that undergoes large-scale melting in planetesimal may be reduced by maximum to 15 km from 30 km.Meanwhile, a planetesimal's instantaneous surface temperature is a function of latitude instead of a constant, resulting in asymmetric heating and a loss of spherical symmetry.This temperature contrast can lead to 3-D convection in melted regions and possibly faster thermal penetration.
In reality, a thermal evolution model involving not only conduction, but also convection, phase transition, and variations in the local thermal and physical properties may be required to model the evolving interior of the planetesimal, at the cost of longer computational time and additional free parameters, for instance, initial composition of the body (Malamud & Perets 2016).

Initial distributions of planetesimals accreted by WDs
Whilst the limits on the maximum size of a planetesimal where melting occurs as a function of distance to the star is robust (given the above assumptions), the fraction of planetesimals accreted by white dwarfs that are core-mantle differentiated is much more uncertain, due to planetesimal distributions and the position/size of the accretion zone.
We apply ) in the simulations, assuming   and  0 are independent variables and the distributions are unaltered by scattering.In terms of spatial distribution (), we utilize the result of Minimum Mass Solar Nebula model (MMSN, Hayashi 1981;Crida 2009), not necessarily applicable to other systems.In terms of size distribution, , we apply the assumption of collisional evolution (Dohnanyi 1969;Gáspár et al. 2012;Pan & Schlichting 2012), which is steeper than that of planetesimals formed in the streaming instability (Simon et al. 2016(Simon et al. , 2017)), with  ≲ 3.
There is no guarantee that these values apply to all planetary systems and to all size ranges.Furthermore, the planetary system can alter the distributions via resonances.In Figure 10, we show how  and  may change our results,  Mm ( crit ) and  Mp ( crit ), in a system with a 1 solar mass host star.With the increase in  (), more masses of the planetesimal population concentrate towards the host star (are occupied by small-sized bodies), leading to a higher degree of differentiation of the whole planetesimal population.Scattering and accretion of planetesimals towards the WD rely on the interactions with the planetary system, which can be nonuniform in space and size.For instance, Li et al. 2022 simulates the accretion of solar system planetesimals onto the Sun when it becomes a white dwarf and concludes that the geometry of the solar system leads to inside-out accretion (the inner bodies are accreted earlier than their outer counterparts).In Figure 11, we present  Mm ( crit ) and  Mp ( crit ) for three subsets of planetesimal population described in Section 3.4.3,initially lying between a: 1.5-3 AU, b: 1.5-5 AU, and c: 3-5 AU around a 1 solar mass host star.In general,  Mm and  Mp decrease rapidly with the outward migration or expansion of the scattering zone.
In this study, we only consider spherically symmetric planetesimals and neglect their shape distribution, as well as the possible alterations to their shapes after melting.In reality, the shape of a planetesimal may affect its orbital and thermal evolution, as well as the scattering and disruption processes, for instance, disruption outside the Roche limit (Katz 2018;Veras et al. 2020;McDonald & Veras 2021), which are beyond the scope of this study.
In the solar system, it is common for asteroids to reside in an elliptical orbit (Malhotra & Wang 2017).In this case, Equation 11 only represents instantaneous equilibrium temperature which varies with the true anomaly.We introduce the time-averaged equilibrium temperature (<  eq >) for one orbital period () (Mé ndez & Rivera-Valentín 2017): (M *, Ca ) lim (kg) Figure 12.The minimum size of the planetesimal ( lim,current ) and the corresponding mass of calcium in the atmospheres of the white dwarf ((  * ,Ca ) lim ), that leads to observable photospheric calcium.The mass fraction of calcium is assumed to be 1.5% in the sample planetesimals.We fit a linear relationship between white dwarf's effective temperature ( eff ) and the detection limit of log 10 (Ca) (Hx) , to the transformed Ca/H and Ca/He (original data in Zuckerman et al. 2003Zuckerman et al. , 2010) ) at a limiting equivalent width of 14 mÅ for a S/N of 30 and spectrograph resolution of 40000 at given  eff .The kink at  eff ∼ 12000 K and 22000 K for white dwarfs with H and He-dominated atmospheres marks the transition from radiative to convective behaviour (Bergeron et al. 2011;Caron et al. 2023).
where  = (1 −  cos ), with  the eccentric anomaly,  2 is the complete elliptic integral of the second kind,  eq () = 11) .Compared to  eq , <  eq > decreases by 10% for identical semi-major axis and  ≳ 0.1.Meanwhile, as the pericentre distance of elliptical orbits ( peri ) satisfies  peri = (1 − ) < , critical  0 increases compared to that for circular orbits.For identical initial semi-major axis distribution, a population of planetesimals on circular orbits are of the highest degree of differentiation with other conditions kept constant.

Observability of planetesimals differentiated under stellar irradiation
The smallest planetesimals, which the model predicts to undergo the largest degree of melting, also produce the weakest signal when accreted, and are not necessarily observable throughout the cooling age of the white dwarf.Our ability to detect calcium (Ca) in a white dwarf atmosphere depends crucially on the signal to noise (S/N) achieved in the observed spectra, alongside the spectral resolution.In order to investigate whether the planetesimals differentiated under stellar irradiation are detectable, we compute the minimum detectable atmospheric Ca/H and Ca/He of the polluted white dwarfs in Zuckerman et al. 2003Zuckerman et al. , 2010 12 and that the observed metal pollution originates from a single body that accreted onto the white dwarf over 1 Myr (upper panel) and 1 Kyr (lower panel) assuming  =  acc in Equation 22, respectively.The mass fraction of calcium is assumed to be 1.5% in the sample planetesimals.
dwarf. ≈ 3.8 × 10 −4 (4.2 × 10 −4 ),  ≈ −14.2(−16.6)for hydrogen (helium) atmospheres.The minimum mass of calcium in the atmosphere of the white dwarf (( * ,Ca ) lim ) that is observable is of the form: where  is the ratio of white dwarf's atmosphere mass to the white dwarf mass, obtained from Montreal White Dwarf Database (MWDD, Dufour et al. 2017;Bédard et al. 2020),  Ca and  Hx are atomic masses of calcium and hydrogen/helium, respectively.We can estimate the minimum size of the individually accreted body ( lim ) that can explain the current atmospheric Ca abundance of the white dwarf by: where  * ,Ca is the mass of Ca in the atmosphere of the white dwarf.We assume that sample planetesimals have a bulk-Earth like calcium mass fraction  ,Ca ≈ 1.5%.The observational limit (minimum size of the accreted body  lim,current , and minimum calcium mass in the atmosphere of the white dwarf, ( * ,Ca ) lim that can lead to observable signals) as a function of WD's effective temperature calculated from Equation 19, 20 and 21 is shown in Figure 12 for WDs with  * = 0.5, 0.6, 0.7  ⊙ and H, He-dominated atmospheres, respectively.Observations of instantaneous accretion are more sensitive to WDs with H-dominated atmosphere until  eff ≲ 8000 K ( lim ∼ 1 km), after which  lim of WDs with He-dominated atmosphere is lower.
In reality, planetesimals are unlikely to accrete instantaneously, rather spreading their accretion over a longer time period that we call  acc (we define the start of accretion to be  = 0, such that  acc also marks the time when the accretion stops).This is important to consider, as the abundance of Ca in the atmosphere of the white dwarf ( * ,Ca ) has a different relationship to the accretion rate (   , assumed to be constant) before and after  =  acc (Koester 2009;Jura & Young 2014): where  Ca is the sinking timescale of Ca obtained from MWDD (Paquette et al. 1986).We consider a special case where  =  acc (the maximum of  * ,Ca ()) in this study and estimate the minimum size of a planetesimal ( lim ) that produces observable Ca lines: where  ,lim = .
By including the sinking timescale of Ca and the accretion timescale of planetesimal debris, we show the corresponding observational limit ( lim ) in Figure 13 with 2 distinct  acc , i: 1 Myr (upper panel) and ii: 1 Kyr (lower panel) in the same parameter space as Figure 12.
In white dwarfs where the pollutants accrete over 1 Myr, the minimum size of a planetesimal that leads to observable Ca lines,   decreases by ≈ 2 orders of magnitude from  eff = 25000 K to 7000 K for white dwarfs with both H and He-dominated atmospheres.If the pollutants come from a single body, accretion of a body with   ≲ 10 km is only detectable in a white dwarf with He-dominated atmosphere and     ≲ 10000 K.
If, instead, pollutants accrete over 1000 yr, the predicted   approaches that of instantaneous accretion (Figure 12) at  eff ≲ 7000 K (20000 K) for white dwarfs with a H (He)-dominated atmospheres.These transitions occur at    ≫   , such that the mass of white dwarf pollutants at the end of accretion is representative of the total mass of the accreted body.The dependence of  lim on  acc is stronger for white dwarfs with H-dominated atmosphere because of their much shorter  Ca compared to their counterparts with He-dominated atmosphere at identical  eff (about 2-3 orders of magnitude).
Large-scale melting (≳ 95% of the volume) triggered by stellar evolution, as is shown in Figure 7, can only occur in planetesimals with   ≲ 30 km even when  crit is lowered to 1400 K. Figure 12 and 13 illustrate that those small (  ≲ 30 km) planetesimals potentially undergo large-scale melting induced by giant branch evolution are only detectable (assuming single body accretion) in the scenario that the body accretes rapidly ( acc ≲ 1 Kyr) or the white dwarf cools down to     ≲ 16000 K.
For He-white dwarfs cooler than ∼ 16000 K, it is possible to detect planetesimals as small as ∼ 30 km, with this limit decreasing to ∼ 1 km when the white dwarf cools to     = 7000 K.For these white dwarfs, more smaller-sized bodies, potentially being more thermally processed become observable, when accreted individually.However, the coolest He-white dwarfs also have extremely long Ca sinking timescale, which adds the possibility that the pollutants consist of several accreted bodies.
Compared to the preferred parameter space where planetesimals differentiate under radiogenic heating (  ≳ 50 km (Curry et al. 2022)), stellar irradiation preferentially differentiate smaller-sized bodies.Therefore, irradiation-differentiated bodies become observable later in the cooling age of the white dwarf, with a definitive upper bound of the absolute pollutant abundances equal to the largest body that may differentiate under stellar irradiation (30 km according our model).

Do we observe post-main sequence differentiation?
Our results show that over-abundances of siderophile/lithophile elements in the atmosphere of a cool white dwarf can result from a parent body differentiated by stellar irradiation, but only in the case that the absolute abundances of WD pollutants drops below that of a 30 km planetesimal.
In this section we further discuss if any observed white dwarfs, whose pollutants are identified as core/mantle-rich, could result from accretion of small-sized planetesimals differentiated during stellar evolution.In order to investigate this we make the broad assumption that the observed atmospheric abundances are dominated by a single body, noting that this may not be the case (Turner & Wyatt 2019).
The majority of white dwarfs with helium-dominated atmospheres identified in the literature as accreting core or mantle-rich (or even crust-rich) material are highly polluted (e.g., SDSSJ0736+4118, SDSSJ0744+4649, Hollands et al. 2017;Harrison et al. 2021a).Such white dwarfs have likely accreted a body of ≳ 50 km in radius to match the observed Ca abundances.Optimistically, a 50 km planetesimal can melt ∼ 40% of its volume under stellar irradiation (Figure 7, right panels).This low degree of differentiation is unlikely to explain the visible siderophile/lithophile character of these white dwarf pollutants.In this case, radiogenic heating, which preferentially differentiates planetesimals with   ≳ 50 km (Curry et al. 2022) may be a better explanation.A notable exception among these He-white dwarfs with core/mantle-rich pollutants is WD0122-227 whose abundances were interpreted as core-rich by Swan et al. 2019;Buchan et al. 2022.The log 10 (Ca) (He) of -10.1 corresponds to the accretion of a ∼ 30 km planetesimal, which if it resided interior to ∼ 1.5 AU of a 1-3  ⊙ star and melted at 1400 K, would have undergone large-scale melting on the giant branches.
Among the white dwarfs with H-dominated atmospheres, WD2105-820, PG0843+516 Swan et al. 2019;Kilic et al. 2020 have been identified as core/mantle-rich and have low metal abundances in their atmospheres, equivalent to accretion of a ≲ 1 km planetesimal.These weakly polluted systems are consistent with the instantaneous accretion of planetesimals sufficiently small to have been differentiated by giant branch heating.However, the uncertainty lies in the short sinking timescale of Ca (ranging from ∼ 10 −2 yr to 10 yr) in the H-dominated atmospheres of white dwarfs, adding the probability that current Ca mass is well below that contained in the parent body whose accretion continues over many sinking timescales.
In general, we do not expect the core/mantle-rich pollutants of white dwarfs with He-dominated atmospheres presented in the literature to date to result from planetesimals differentiated during giant branches.Considering the much longer sinking timescale compared to H-white dwarfs, future observations targeting at pollutants in the He-dominated atmospheres of white dwarfs with over-abundances in siderophile/lithophile elements but low absolute abundances of Ca may provide more evidence for stellar-induced differentiation.

Depletion of moderately volatiles
The degree of depletion in moderately volatiles, including elements such as Mn, Na, is considered as a powerful tool to probe the early thermal history of a planetesimal (Jura & Young 2014).For instance, bulk Earth is depleted in these species relative to chondritic meteorites (Harrison et al. 2018), most likely explained by the incomplete condensation of Na-bearing minerals out of the nebula gas early in planet formation (O'Neill & Palme 2008; Pringle et al. 2014;Siebert et al. 2018).Planetary materials accreted by white dwarfs are often Na-poor compared to solar compositions.(Doyle et al. 2020;Harrison et al. 2021a) and Harrison et al. 2021b argue that depletion in moderately volatiles could occur early in planet formation, due to incomplete condensation of the nebula gas, or after the nebula gas has dissipated due to large-scale melting in magma oceans following impacts or radiogenic heating.
For those planetesimals that undergo large-scale melting to form a magma ocean induced by stellar irradiation on the giant branches (Section 3.4.2),additional loss of moderately volatiles such as Na and Mn is anticipated, as well as segregation of the iron from the silicates.For those planetesimals heated insufficiently for the silicates to melt and iron to mobilise, depletion of moderately volatiles such as Na largely depends on the thermal properties of the host minerals.Hence, it is unclear whether and how much the moderately volatiles will be lost.If Na is hosted in refractory silicates, whose melting is necessary for Na to be movable, the critical temperature of Na depletion may be coincide with/exceed the melting temperature of silicates.In this case, the observed white dwarf pollutants depleted in moderately volatiles is a signature of heating during the planetary formation stage (Harrison et al. 2021b).Meanwhile, if Na (partially) exists in the minerals of lower thermal stability, as proposed by Masiero et al. 2021 who suggest Na hosted in sodalite or nepheline could degas above ∼ 1000 K, although a high level of moderately volatile-element depletion in massive white dwarf pollutants almost certainly originate from heating around planetary formation, it is unclear whether stellar irradiation could be the cause of moderately volatile-element depletion in the less massive (≲ mass of a 50 km planetesimal) white dwarf pollutants.

CONCLUSION
Many white dwarfs have accreted planetary materials that are rich in either core or mantle material.The formation of iron cores requires a period of large-scale melting.This work shows that very few planetesimals orbiting white dwarfs underwent large-scale melting triggered by stellar irradiation on the giant branches.For a solar like host star, a planetesimal must initially reside between ∼ 1.3 AU and 1.7 AU and be smaller than ∼ 30 km in radius in order to be melted over 95% of the volume and survive the asymptotic giant branch of the host star, even when the critical temperature of melting is lowered to 1400 K.
This work highlights planetesimal size as a key differentiator between large-scale melting that occurs due to heating on the giant branches and that due to the decay of 26 Al soon after planet formation.For Solar System abundances of 26 Al, the latter prefers bodies ≳ 50 km in radius, whilst the former occurs only in bodies where   ≲ 30 km.Thus, accretion of these small planetesimals that underwent large-scale melting due to giant branch heating are most likely to be seen in cool white dwarfs, especially those with relatively low absolute abundances of pollutants.They do not represent the current observed population of core/mantle-rich white dwarf pollutants.(B5) The evolution of semi-major axis accounting for both tidal decay and stellar mass loss can be generalised as :  =  ,tide +  * ,tide −  *  *

𝑎. (B6)
The overlap of four curves in Figure B2 verifies that tidal interactions have negligible effect on planetesimal-sized bodies on eccentric orbits avoiding engulfment, and the orbital evolution of these planetesimals are dominated by stellar mass loss.

APPENDIX C: GRAVITATIONAL SEGREGATION
The equation of motion of over-dense fragments (over-density Δ, radius   ) is of the form (Stokes 2009;Falkovich 2018

APPENDIX D: THE YARKOVSKY EFFECT AND THE YORP EFFECT
The maximum effect of the Yarkovsky and the YORP effect can be approximated as (Veras et al. 2014(Veras et al. , 2019(Veras et al. , 2022;;Ferich et al. 2022 asymmetry of the planetesimal.We adopt  = 10 −3 ,  * = 2  ⊙ ,  0 = 1.5 AU and  = 0 and switch on the Yarkovsky effect and the YORP effect from the AGB life of the system.As is shown in Figure D1, for the smallest planetesimals in our sample (  = 10 km ), the outward/inward Yarkovsky drift can reach ∼ 0.4 AU at the end of TPAGB life of the system, and the maximum YORP spin up is around 0.5  break (= √︃ Figure1.Thermal history of white dwarf pollutants: formation from protoplanetary disks, main sequence and post-main sequence life of the star, scattering into the Roche limit, circumstellar disks formation, accretion.The color bars show the general luminosity and radius variations of the host star, and the thermally pulsing asymptotic giant branch stellar wind density.

Figure 2 .
Figure2.Key steps of methodology: a synthesis of stellar and planetesimal's orbital evolution giving the stellar luminosity ( * ( )) and planetesimal's position ( ( )), which together provide the stellar irradiation, whose strength is quantified as the planetesimal's surface temperature ( (  = ±  ,  )).The thermal evolution of the planetesimal is then simulated to predict its internal temperature evolution ( ( ,  )), which reveals the fraction of the body that undergoes large-scale melting and differentiation.The roundrects and rectangles representing numerical simulations and outputs, respectively. is the radial coordinate inside the planetesimal,  is the time coordinate.

Figure 3 .
Figure 3. Evolution of stars with  * = 1  ⊙ and 2  ⊙ and the resultant equilibrium temperature for planetesimals with  0 = 1.5 AU and  = 0.The zoom-in plots correspond to thermally pulsing asymptotic giant branch (TPAGB) with the origin of the time axes marks the start of TPAGB phase.The age of the star at the start of of its TPAGB phase ( 0 ) is added.

Figure 4 .
Figure 4.The upper panel shows the temperature evolution at the centre ( = 0), half of the radius ( = 2 ) and the surface ( = ) of a planetesimal in sample system IA and IIA (Table4).The period of each evolutionary phase is normalised by its overall timescale, with the absolute time along the horizontal axis of the plot.The lower panel shows the temperature profile evolution of planetesimals in sample system IA, IB and IIA.The contour plots start from the red giant branch (RGB) of the host star, followed by the core helium burning phase (CHB) and the asymptotic giant branch (AGB), with the age of the star increasing clock-wisely.Same normalisation as the upper panel is applied.

Figure 6 .
Figure 6.The maximum temperature at the centre ( = 0) and half of the radius ( =  2 ) of the planetesimal, throughout its thermal history under various circumstances.The blue (red) dashed (dotted) line corresponds to the critical orbit after (before) smoothing pulsations.The positions of sample IA (black circle) IB (triangle) and IIA (square) in Figure 4 are added.

Figure 7 .
Figure 7. Volume fraction that is once melted (   ) for planetesimals of different sizes and initial positions.The blue (red) dashed (dotted) lines correspond to the critical orbit after (before) smoothing stellar thermal pulsations.The black dashed vertical line corresponds to the semi-major axis beyond which the maximum surface temperature of planetesimals drop below the given melting threshold.The positions of sample IA (black circle) IB (triangle) and IIA (square) in Figure 4 are added.

Figure A1 .Figure B1 .
Figure A1.Parameter  in semi-major axis-time space of a 3 solar mass star in its thermal pulsing asymptotic giant branch.The black solid curve corresponds to the radius of the host star.

Figure B2 .Figure C1 .
Figure B2.Comparison between the orbital evolution of a planetesimal including and excluding tidal effect from zero-age main sequence to the end of AGB life of the system.Four curves overlap with one another and cannot be visually distinguished on the plot.  = 300 km, the pericentre is fixed at 1.33 AU (critical distance of engulfment),  * = 1  ⊙ .The zoom-in plot corresponds to the asymptotic giant branch (AGB) with the origin of time marks the start of this phase.0 is the age of the star at the start of its AGB phase.

=,√ 1 −Figure D1 .
Figure D1.The maximum increase/decrease in the planetesimal's semi-major axis due to Yarkovsky effect (solid lines) and the maximum spin up of the planetesimal (dotted lines) due to YORP effect during the TPAGB phase of the host star. = 10 −3 ,  * = 2  ⊙ ,  0 = 1.5 AU and  = 0.  break = √︃ 4  3 .The inward drift case experiences slightly stronger Yarkovsky effect and YORP effect compared to its outward counterpart because of the stronger stellar flux ∝  *  2 .

Table 1 .
Constant thermal properties of sample planetesimals.

Table 3 .
Comparison of the maximum stellar luminosity and radius during the main sequence, red giant branch and asymptotic giant branch between the 1 and 2 solar mass star.

Table 4 .
List of the parameters of 3 sample systems in Figure4.
The minimum size of the planetesimal resulting in detected photospheric calcium assuming the same detection limits as Figure env is the mass of the convective envelope and  conv is the corresponding convective time-scale:  env being the size of the convective envelope (∼  * ), and   ∼ 3.   is given by:  =  ′  1, 2      conv   ∼ 2,  ′ ∼ 9 2 ,   ∼ 1, and   =  √︃  (   + * )  3. Meanwhile, orbital decay due to planetesimal tide is approximated by constant time lag (CTL) model.The orbital evolution of a planetesimal in 1:1 spin-orbit resonance is of the form:  is the modified planetesimal tidal quality factor,  is eccentricity functions of the form: