The Solar wind prevents re-accretion of debris after Mercury's giant impact

The planet Mercury possesses an anomalously large iron core, and a correspondingly high bulk density. Numerous hypotheses have been proposed in order to explain such a large iron content. A long-standing idea holds that Mercury once possessed a larger silicate mantle which was removed by a giant impact early in the the Solar system's history. A central problem with this idea has been that material ejected from Mercury is typically re-accreted onto the planet after a short (~Myr) timescale. Here, we show that the primordial Solar wind would have provided sufficient drag upon ejected debris to remove them from Mercury-crossing trajectories before re-impacting the planet's surface. Specifically, the young Sun likely possessed a stronger wind, fast rotation and strong magnetic field. Depending upon the time of the giant impact, the ram pressure associated with this wind would push particles outward into the Solar system, or inward toward the Sun, on sub-Myr timescales, depending upon the size of ejected debris. Accordingly, the giant impact hypothesis remains a viable pathway toward the removal of planetary mantles, both on Mercury and extrasolar planets, particularly those close to young stars with strong winds.


Introduction
The compositions of Mercury, Venus, Earth, and Mars, serve as direct windows into the conditions that persisted during the opening 100 million years of our solar system's formation. Inferred mass ratios between iron and silicates for Venus, Earth, and Mars are all roughly consistent with a common, chondritic abundance (Righter et al. 2006;Ebel & Stewart 2018). In contrast, Mercury exhibits an anomalously high iron content. Mercury's iron core possesses a radius exceeding 80% of the planet's radius, as compared to Earth's more modest 50%, alongside a correspondingly high bulk density that rivals Earth's, despite the order of magnitude difference in planetary masses (Ash et al. 1971). The explanation for such a disparity in iron content remains mysterious.
Over the decades, many hypotheses have been brought forward to understand Mercury's large core (Solomon 2003;Ebel & Stewart 2018). Broadly, these ideas separate into those that suppose Mercury's high iron-to-silicate ratio is primordial, versus those that propose Mercury formed with a chondritic composition, but that the silicates were subsequently removed, usually by collisions. Primordial scenarios, in contrast, draw upon mechanisms that separate iron and silicates within the protoplanetary material before becoming incorporated into a planet.
Perhaps the earliest class of ideas supposed that nebular conditions were at such a temperature where silicates were gaseous, but iron had condensed (Urey 1951;Lewis 1973). A more recent proposal, still emphasizing Mercury's close-in, high-temperature position in the disk, suggested that the solar nebular conditions became sufficiently hot at Mercury's orbital distance that Mercury's mantle was vaporized and removed along with the nebular gas (Cameron 1985). In general, temperature-driven mechanisms are inconsistent with modern views of nebular temperature distributions during planet formation (Armitage 2011;Hartmann et al. 2016), and appear to contradict the widespread existence of extrasolar planets with periods shorter than Mercury's (Borucki 2016).
Alternative fractionation pathways exist that do not rely upon temperature, such as photophoresis (an effect related to conductivity; Wurm et al. 2013), gas drag (separation of densities; Weidenschilling 1978), or magnetic erosion (magnetic attraction of Fe-rich material; Hubbard 2014). Whereas these mechanisms may reproduce Mercury's Fe/Si ratio in principle, they typically require special nebular conditions that are either not expected, or contradictory to other findings. Research continues to develop in the field of circumstellar disk chemistry such that primordial fractionation of Si and Fe is not necessarily ruled out.
In this work, we focus on a separate class of post-formational models that suppose Mercury to have initially conglomerated with a chondritic elemental composition, similar to that of the other terrestrial planets. Subsequently, the silicate mantle was stripped by way of a giant impact (Wetherill 1985;Benz et al. 1988Benz et al. , 2007Chapman 1988). This model is appealing because an epoch of giant impacts is expected during the latter stages of most envisioned planet formation scenarios (Raymond et al. 2004(Raymond et al. , 2018. Indeed, giant collisions continue to serve as viable explanations for many other solar system mysteries, including the formation of the Earth's Moon (Hartmann & Davis 1975;Canup & Asphaug 2001;Ćuk & Stewart 2012;Lock et al. 2016), the Martian moons (Citron et al. 2015;Hyodo et al. 2017) the Pluto-Charon system (Canup 2005), and even the tilt of Uranus (Safronov 1966;Morbidelli et al. 2012).
Generally speaking, simulations of Mercury's giant impact have shown that a high-energy impact onto a proto-Mercury of about 2.25 times Mercury's current mass is capable of removing the required mass of silicates to match Mercury's bulk density (Benz et al. 1988(Benz et al. , 2007Chau et al. 2018). However, a major drawback of the giant impact theory is that much of the material launched from Mercury as vapor tends to condense into solid spherules (Johnson & Melosh 2012) and remain on Mercury-crossing trajectories in the aftermath of the collision. Consequently, the ejected material simply reaccretes over ∼10 Myr timescales (Benz et al. 2007;Gladman & Coffey 2009) and the final planetary density is not sufficiently enhanced.
Proposed solutions to the problem of reaccretion have called upon Poynting-Robertson drag, which acts on a similar timescale to reaccretion for centimeter-sized particles, the size of debris expected from the condensing ejecta cloud (Burns et al. 1979;Benz et al. 2007). However, the high optical depth of the heliocentric ejecta ring can significantly self-shield (Gladman & Coffey 2009) leaving it unclear whether Poynting-Robertson drag alone is enough to prevent reaccretion. This problem has motivated several modifications of the initial impact theory. In one scenario, a larger impactor survives the collision intact, having stripped Mercury's mantle (Asphaug et al. 2006;Asphaug & Reufer 2014). In another, numerous smaller impacts erode Mercury's surface (Svetsov 2011).
Here, we present the hypothesis that the solar wind constitutes a robust source of orbital decay for material ejected from Mercury. The solar wind consists of a nearly radial outflow of ionized particles, traveling outward into the solar system at hundreds of km s −1 (Phillips et al. 1995). Any objects in orbit must pass through a nonzero interplanetary density of plasma, inducing a drag force upon their orbits (Burns et al. 1979;Mukai & Yamamoto 1982).
Today, the solar wind-induced drag felt by particles larger than a few microns is roughly one-third of that from Poynting-Robertson drag. However, if the young Sun possessed a substantially stronger wind than today, the solar wind drag would increase in turn. It remains difficult to determine the wind properties of stars during the planet-forming period because few direct measurements exist of such early winds (Jardine & Collier Cameron 2018). Nevertheless winds from Sun-like stars are expected to decay with time, with measurements suggesting 100 times the modern mass-loss rate in younger solar analogs (Wood et al. 2014). In contrast, the magnitude of Poynting-Robertson drag may only have varied by a factor of a few in the early solar system, as the Sun's luminosity dropped during the Hayashi track, arriving onto the pre-main sequence slightly dimmer than today (Maeder 2008).
In this work, we show that increasing the solar wind by a factor of 10-100 compared to modern leads to removal of ejecta from Mercury within ∼1 million year timescales, safely within the 10 Myr reaccretion time. In Sections 2 and 3, we construct a simple, axisymmetric solar wind model, and compute its impact upon ejected debris in Section 4. In Section 5, we demonstrate the effect of the solar wind drag within N-body simulations, and in Section 6 we discuss the implications of the solar wind hypothesis.

Physical Setup
In this section, we outline the physical scenario describing the trajectories of particles ejected from Mercury (see Figure 1 for a schematic of the problem). We suppose that the precursor of Mercury is impacted by another object, ejecting a large fraction of Mercury's silicate-rich mantle into heliocentric orbit with some velocity v ej (Benz et al. 1988;Chapman 1988;Gladman & Coffey 2009). Once in orbit, the material cools and condenses into particles, assumed non-interacting spheres of radius s and density ρ s =3 g cm −3 . Each particle follows a Keplerian orbit, of semimajor axis a, eccentricity e, and velocity (Murray & Dermott 1999) = - a  e  e  f  e  f  e  f  1  1  cos  sin  1 cos , where we have defined the true anomaly f, stellar mass M å , gravitational constant G, and f e r, are unit vectors in the plane of Mercury's orbit (Figure 1). Subsequent to ejection, the azimuthal ram pressure from the solar wind plasma leads to orbital evolution of the particles. We construct a model for the young Sun's wind in order to compute whether the solar wind can remove debris from Mercury-crossing trajectories before reimpacting Mercury's surface (Gladman & Coffey 2009). For the purposes of the analysis that follows, we assume that this ejecta remains in the plane of Mercury's orbit, thus simplifying the interaction between ejecta and solar wind plasma. This coplanar assumption will be lifted when we turn to N-body simulations in Section 5.
The modern solar wind consists of an ionized plasma expanding into the interplanetary medium with a predominantly radial velocity, which ranges from ∼400 to 600 km s −1 (the slow and fast wind; Phillips et al. 1995). The magnitude of mass loss in the modern solar system is Mass-loss rates of Sun-like stars appear to decrease with time, with measured values of about 100 times values in solar analogs at ages of 100 s of millions of years (Wood et al. 2014). However, few measurements exist of stellar winds during the first 100 million years subsequent to disk-dispersal. During the disk-hosting stage itself, star-disk interactions may increase winds to upwards of ∼10 5 times modern (White & Hillenbrand 2004;Matt & Pudritz 2008;Plavchan et al. 2009), but it is unlikely that such strong wind persists long after the nebular gas has been removed. Accordingly, stellar wind magnitudes during the first 100 million years of planet formation remain mysterious, with estimates ranging from modern levels up to 10,000 times modern (Cohen & Drake 2014;Ó Fionnagáin & Vidotto 2018).
In addition to the magnitude of mass loss, young solar analogs generally differ from the Sun by possessing shorter spin periods (ranging between 1 and 10 days; Bouvier et al. 2014) and stronger magnetic field strengths (0.1-1 kGauss; Donati & Landstreet 2009;Vidotto et al. 2014). The rapidly rotating magnetosphere would accelerate wind plasma to drive a significant azimuthal component to the wind velocity that we compute below (Lovelace et al. 2008;Spalding 2018;Carolan et al. 2019).
A small azimuthal component of the wind's velocity would reduce the net drag upon orbiting debris. However, if the stellar magnetosphere was strong enough, its wind may reach super-Keplerian velocities, causing an acceleration of orbiting material instead of a drag. We suppose that the solar wind plasma possesses a velocity v sw , this would result in a ram pressure of r | upon orbiting debris, where C D is a drag coefficient. We set C D =1 (Mukai & Yamamoto 1982), but its exact value depends upon the particle's composition and size.
The above ram pressure leads to a force of where the wind-facing surface area of the particles is taken as πs 2 . We may obtain a first order approximation to F by assuming null eccentricity, together with considering the limits where the azimuthal component of the wind is negligible (as it is today) and the wind's radial velocity v r greatly exceeds the Keplerian velocity. Finally, we approximate the wind as spherical to write the plasma density as With the above simplifications, the drag given by Equation (2) may be used to compute a timescale of orbital decay (Spalding 2018): where M  is the mass-loss rate of the star. Notice that under the current level of approximation, the drift timescale is independent of the wind's velocity. The modern Sun's mass loss is = »´-M M 2 10 14    M e yr −1 (Phillips et al. 1995), but as mentioned above was likely exceeded by orders of magnitude during terrestrial planet formation. Drift timescales are proportional to particle radius s, which is likely to vary with impactor parameters and the composition of debris (Melosh & Vickery 1991;Benz et al. 2007;Johnson & Melosh 2012;Chau et al. 2018). By combining Smoothed-Particle Hydrodynamics (SPH) simulations with thermodynamic computations, Benz et al. (2007) calculated that ejected particles possess typical sizes of s∼1 cm and it is this value that we assume throughout. However, we return to this problem in Section 6.2 to discuss avenues toward updated estimates.
Given the parameters above, we find that the drift timescale takes the magnitude which is plotted for various orbital radii in Figure 2. As can be seen, solar winds exceeding roughly 10-20 times the present magnitude are sufficient to cause outward migration of centimeter-sized particles over million year timescales. The reaccretion timescale of ejecta found previously is on the order of 1-10 Myr (Gladman & Coffey 2009). Thus, from the simple treatment above, the solar wind constituted a significant component to the orbital evolution of material ejected from the primordial Mercury subsequent to a purported giant collision.
Despite the promising results of timescale considerations, ejecta from Mercury are unlikely to be launched onto circular, heliocentric orbits. Moreover, we have not yet included the Figure 1. Schematic of the physical scenario considered in this work. Material is launched from Mercury (1) subsequent to a giant impact. This material initially follows the outer, red orbit that crosses the orbit of Mercury (gray ellipse). A drag from solar wind plasma draws angular momentum from the particle's orbit (2), causing the orbit to spiral inward toward the Sun over time (3), and eventually detach from Mercury's orbit (4), inhibiting reaccretion. The zoom-in on the right denotes the coordinate frame used to prescribe the initial trajectories of collisional debris in our N-body simulations (Section 5). The timescale of in-spiral grows linearly with particle radius, where a radius of 1 cm is assumed throughout (see Section 6.2 for a discussion of particle sizes). effect of the azimuthal component of the wind velocity, which may increase or decrease drift timescales depending upon the magnitude and sign of the quantity v sw -v K .

Stellar Wind Model
One of the earliest quantitative models of the solar wind considered a one-dimensional, steady-state approximation of an isothermal gas moving radially outward from an unmagnetized, nonrotating star (Parker 1965). Such a model is remarkably efficient at recreating the radial velocity of the solar wind, but is incapable of deriving an azimuthal velocity. A two-dimensional, axisymmetric model that takes rotation and magnetism into account was developed soon after (Weber & Davis 1967). Such a "Weber-Davis" model assumed that the wind consisted of an infinitely conductive plasma that interacted with the star's rotating magnetosphere through the principles of ideal magnetohydrodynamics.
In the intervening decades, solar wind models have continuously grown in sophistication, including fully threedimensional magnetohydrodynamics (MHD) simulations (Cranmer & Saar 2011;Cohen & Drake 2014;van der Holst et al. 2014). Our primary goal for this work is to compute an approximate radial profile of the wind's azimuthal velocity, given a stellar mass-loss rate, magnetic field, and spin rate. Thus, the Weber-Davis model will suffice for our purposes, but we draw the evolution of M  , magnetic field and spin rate from published observations (Vidotto et al. 2014) and three-dimensional MHD simulations (Ó Fionnagáin & Vidotto 2018).

Weber-Davis Model
The picture of the solar wind we employ considers the wind plasma to expand radially outwards along the stellar equator. Its outward motion is governed by both pressure and magnetic stresses, but its azimuthal acceleration is dominated by magnetic torques. A pathway toward a solution of this "Weber-Davis" model has been presented in numerous works (Weber & Davis 1967;Hartmann & MacGregor 1982;Lamers et al. 1999) and so we only provide a brief outline of the solution here. Closely following Hartmann & MacGregor (1982), we assume that both the magnetic field B and wind velocity are axisymmetric, and possess only radial (v r , B r ) and where B r,0 is the magnetic field strength at the stellar photosphere, taken to be when r=R å , the stellar radius. Note that all calculations are carried out in the equatorial plane. If Faraday's law is combined with the assumptions of ideal MHD (i.e., infinite conductivity), we obtain a steady-state relationship of where the stellar angular velocity W  has been introduced. For simplicity, we assume that the wind is isothermal, with isothermal sound speed given in terms of the photosphere temperature T=10 6 K, the mean molecular number μ=0.6, the mass of a proton m H , and the Boltzmann constant k (Lamers et al. 1999).
Each parcel of plasma carries with it a constant specific angular momentum L and energy E, which are derived from the steady-state, azimuthal and radial components of the momentum equation respectively: The conserved angular momentum is extracted via integration as where the Alfvén radius r A is defined as the point where the radial wind velocity is equal to the Alfvén velocity, defined as m r º v B r r A, 0 sw (Lamers et al. 1999;Bellan 2008). Likewise, integration yields the conserved energy Finally, rearranging the momentum equation yields to the relationship 0 sw . By inspection of Equation (12), any solution passing through a point where the denominator vanishes must also cause the numerator to vanish. This situation occurs at 2 radii, which are r f (the fast point) and r s (the slow point), and physically correspond to the transition of the radial wind velocity through the fast and slow magnetosonic points (Weber & Davis 1967;Hartmann & MacGregor 1982;Lamers et al. 1999;Bellan 2008).
The requirement of vanishing numerator and denominator at both points constitutes a system of four simultaneous equations. Two further equations arise from the constancy of the energy E at the stellar surface r=R å and at the Alfvén radius r=r A . Altogether, these six equations must be solved for unknowns, which include the three radii r r r , , , the values of radial wind velocity at r f and r s , and finally the velocity of the wind at the stellar photosphere. We obtain these six parameters by way of a Newton-Rhapson solver, from which the solar wind properties, and thus the ejecta drift timescales, may be obtained as a function of radius.
Illustrative radial profiles of v r and v f are presented in Figure 3 at three epochs-3, 10, and 30 Myr. The radial velocity increases monotonically with distance, reaching a constant as  ¥ r . In contrast, the azimuthal velocity displays an extremum in velocity slightly interior to the Alfvén radius (vertical dotted lines). Physically, this peak arises because close to the star the magnetic field is strong enough to accelerate an outwardly expanding ring of plasma over a timescale that is short compared to the radial expansion timescale of the plasma. Approximately, this strong coupling locks plasma to solid-body rotation close to the star's surface (denoted by the slanting, dashed red line). At larger distances (r > r A ), magnetic coupling weakens, such that the expanding ring of plasma no longer receives significant acceleration from the weakening magnetic stresses. In order to approximately conserve its angular momentum v f approximately falls as 1/r (Carolan et al. 2019).
At times earlier than about 10 Myr, a substantial range of radii exist where v f exceeds the Keplerian velocity (the black line in the bottom panel). We discuss the significance of this in greater detail below, but qualitatively it suggests that debris occupying Mercury's orbit at earlier times are more likely to experience orbital expansion than at later epochs, when orbital contraction is favored.
Up until recently, in situ measurements of the solar wind's azimuthal velocity were unavailable at distances comparable to Mercury's orbit, and so the theoretically expected peak in v f could not tested. However, the Parker Solar Probe recently detected azimuthal winds, increasing toward the Sun at distances inward of 0.2 au (Kasper et al. 2019). Surprisingly, these winds possessed azimuthal components of up to 40 km s −1 , at least an order of magnitude above those predicted for the Sun at its current age. Indeed, these speeds are more consistent with the winds we compute at an age of 30 Myr (bottom panel of Figure 3). It is impossible to tell as of yet how these theoretical underestimations propagate to younger, more active stars, and thus whether or not our calculations serve as lower limits of the early wind strength. Undoubtedly, the upcoming closer passes of the Solar Probe to the Sun will help shed light on this unfolding mystery.

Early Stellar Properties
The epoch of giant impacts persisted throughout approximately the first 100 million years of solar system history (Raymond et al. 2018), though the exact time of Mercury's hypothesized impact is uncertain (see Section 6.1). As mentioned above, young stars typically possess spin periods between 1 and 10 days and magnetic fields between 0.1 and 1 kGauss (Bouvier et al. 2014;Vidotto et al. 2014). A less wellconstrained aspect lies in choosing the mass-loss rate M  , alongside its time-dependence. Given the lack of direct measurements, previous work has typically resorted to 3D magnetohydrodynamic models from which the mass loss emerges as an output of the model (Cohen & Drake 2014;Ó Fionnagáin & Vidotto 2018).
Unfortunately, these computational approaches do not approach a consensus, with estimates varying widely, though generally the periods P å , M  , and B å all decrease with time. We with magnitudes that are chosen to match observed T Tauri magnetic fields and spin periods, as well as allowing the wind to pass through M 100   during the period of terrestrial planet formation: Though the uncertainty in the above parameters is large, it is important to point out that the time at which Mercury's purported giant impact occurred is currently unconstrained. Thus, even if an exact knowledge existed of the young Sun's stellar parameters, the conditions persisting during Mercury's giant impact would remain with large uncertainties. Accordingly, our goal is to explore whether a broad range of stellar parameters and impact times are consistent with removal of ejected material. In addition to the above parameters, the stellar radius itself would be changing throughout the terrestrial planet formation period. Specifically, the Sun contracts along the fully convective Hayashi track, entering the Henyey track after several tens of Myr (Shu et al. 1987;Maeder 2008). Using stellar models from the online database of Siess et al. (2000), we computed the time evolution of a Sun-like star's radius. We find that this contraction is well-approximated by the equation where R e is the present-day Sun's radius. The decay of the stellar radius provides a significant reduction in the effect of the wind, because of the  R r 2 2 dependence of B r . With the above prescriptions, the stellar wind properties as a function of r and t are fully specified. We now turn to the influence this wind has upon material ejected from Mercury's surface.

Orbital Evolution of Ejecta
In the previous section, we obtained a prescription for the solar wind as a function of heliocentric distance and time. Here we use these results to compute the orbital evolution of ejecta launched from Mercury. The drag F experienced by ejecta is given in vector form by Equation (2). There exist analytic expressions that relate the radial F r and azimuthal F f components of this drag to the time evolution of a particle's eccentricity e and semimajor axis a, which take the form (Burns 1976;Hedman 2018 . Unfortunately, the above is not a trivial computation when particles possess significant eccentricities, which is likely to be the case appropriate for ejecta flung from Mercury into heliocentric orbit (Gladman & Coffey 2009). To illustrate this, suppose that a particle is ejected with velocity v ej relative to Mercury, directed in Mercury's orbital plane but at an angle α to Mercury's orbital velocity. For simplicity, we assume Mercury's orbit to be circular, such that the eccentricity and semimajor axis of the particle are given by where b º v v ej Kep is the normalized relative ejecta velocity. Using the relationships above, we may write down a formula for the orbital and eccentricity decay time as a function of the particle's ejected velocity and direction: These timescales are illustrated in Figure 4 in order to depict the dependence of drift and eccentricity decay timescale upon the velocity of ejection. As illustrated in Figure 4, the eccentricity-damping timescale is similar in both sign and magnitude to τ a . Interestingly, this means that during the first few million years, when azimuthal wind velocities are high, semimajor axes and eccentricities of debris both grow at similar rates. At later times, they both damp at similar rates. We do not explore the eccentricity evolution in great detail here, but it serves as a potentially important dynamical influence upon debris disks around young stars.
In Figure 5, we set e=0 and plot the orbital evolution timescale owing to wind-induced drag at Mercury's orbit as a function of the age of the system (beginning at 3 Myr). For comparison, we include the timescale associated with Poynting-Robertson drag, which takes the form (Burns et al. 1979 in terms of the Sun's luminosity L=3.8×10 26 W and the speed of light c=3×10 8 m s −1 . As mentioned above, the Sun's luminosity will slightly exceed the modern value at 3 Myr and drop slightly below it at 100 Myr (Siess et al. 2000;Maeder 2008). Given the relatively small magnitude of such a variation we simply use the Sun's current luminosity for comparison. The drift timescale associated with the solar wind is significantly shorter than that associated with Poynting-Robertson drag over most of the wind's evolution. However, it is worth noting that between about 10 and 20 Myr, the windinduced timescale grows substantially. This occurs because the wind transitions from a super-Keplerian to a sub-Keplerian configuration, which passes briefly through a phase where the relative velocity between debris and wind is small.
Most importantly for our work, the magnitude of τ a is consistently shorter than 1 Myr across a wide range of launch parameters, dropping below 0.2 Myr for low launch velocities. Accordingly, over the entire feasible range of launch velocities, ejecta is removed over a timescale shorter than the 10 Myr reaccretion time presented previously (Gladman & Coffey 2009).

Numerical Simulations
In the quantitative treatment above, we concluded that the orbital evolution timescales for centimeter-sized particles ejected from the surface of Mercury typically lie in the range of 0.1-1 Myr. This semianalytic treatment suggests that the ejecta would be removed from Mercury-crossing trajectories before reaccreting. However, perturbations from Mercury itself, as well as from the other planets complicate this picture beyond the above scenario of a single particle orbiting the Sun. Secular and mean motion resonances permeate the solar system, potentially affecting the evolution of debris (Gladman et al. 1996).
In order to test our analytic results within a more realistic early solar system scenario, we perform N-body simulations of the ejecta orbits, subject to the gravitational forces of Mercury, Venus, Earth, Mars, Jupiter, and the Sun. We carry out simulations, both with and without the solar wind-induced drag in order to compute the efficiency with which the solar wind may remove debris from Mercury's vicinity.

Computational Set Up
All simulations were carried out using the REBOUND symplectic integration package, where the hybrid integration scheme "Mercurius" was employed (Rein & Liu 2012;Rein & Tamayo 2015). This integration scheme conducts a WHFast algorithm (Wisdom & Holman 1991), where our chosen timestep was 0.2 days, transitioning to an IAS15 algorithm to deal with close encounters (Rein & Spiegel 2014). Multiple timesteps were checked, ranging from 2 to 0.2 days, but there was no systematic dependence observed between collision rate and step size. In order to model solar windinduced drag, we prescribed exponential evolution of both the eccentricities and semimajor axes of ejecta particles by way of the additional forces options within the REBOUNDx package (Tamayo et al. 2020).  We began with a control simulation that did not include the solar wind, similarly to Gladman & Coffey (2009). These integrations were performed for 10 million years each, and 2 particle ejection velocities (the velocity of the ejecta relative to Mercury) were chosen, v ej ={4, 10} km s −1 . In parallel, we simulated these two launch velocities under the action of the solar wind. In these simulations, eccentricities and semimajor axes ejecta orbits were forced to exponentially decay over respective timescales of τ e and τ a , computed according to Equation (19).
Given that eccentricity-damping timescales computed above are relatively fast, and for the sake of computational speed, we employed Equation (19) to lowest order in eccentricity, such that ( )  and assumed that s=1 cm, the typical size of debris computed by Benz et al. (2007). Accordingly, the drift times are and were updated with the particles' new semimajor axes after each 1000 yr interval throughout the simulation. We simulated two cases for each of v ej =4 km s −1 and v ej =10 km s −1 , with = M M 10, 100 { }   . Cumulatively, we ran three cases for each of the two launch velocities, with these three cases denoted as "weak wind," "strong wind," and "wind-free" in Figure 6.

Initial Conditions
Each simulation began with 110 particles, ejected isotropically from Mercury. During the 10 Myr simulation these particles were removed if either (1) they collided with one of the six massive bodies (including the Sun), or (2) their apocenter fell below 0.2 au (which we considered as safely removed from Mercury's vicinity). Our earlier, analytic treatment assumed coplanarity between the ejecta orbits and Mercury's orbit. This assumption is easily lifted within N-body simulations. If the planets and ejecta were assumed coplanar, reaccretion timescales would be deceptively short, owing to the enhanced probability of collision between coplanar orbits as compared to mutually inclined orbits (Farinella & Donald 1992).
We initialize the planetary orbits with similar eccentricities, inclinations, and semimajor axes to their modern values, but with uniformly randomized mean anomalies, arguments of pericenter, and longitudes of ascending node. The exception is Mercury itself, which we began with mean anomaly, argument of pericenter, and longitude of ascending node all set to zero, i.e., the planet began at y=0 on our coordinate grid (Figure 1). This set-up was simply to aid in prescribing the initial particle orbits relative to Mercury, which we now describe.
With respect to the ejected particles, we adopted a uniform distribution of ejecta angles over 4π steradians. Specifically, we chose N=11 angles, uniformly spaced in ranges 0<f i 2π and N−1 angles from q -< < 1 cos 1 j ( ) , where θ i and f j take their usual definitions within a spherical coordinate system centered upon Mercury 4 and are defined as yielding N(N-1)=110 ejected particles. Given the above choices of angles and ejecta velocities, we set up the particles on heliocentric orbits. Having chosen the angles θ i and f j , we initialize each particle with a position x sin cos , sin sin , cos 24 relative to Mercury, such that the initial heliocentric position is + x x ij M where x M is Mercury's initial position. Likewise, the initial velocity of the particle relative to Mercury is given by To ensure that the particle is well clear of Mercury at the beginning of the integration, we set x 0 =0.006 au, which is approximately 4 Hill radii from Mercury. The value of v 0 took either of 4, 10 { } km s −1 .

N-body Results
The cumulative number of collisions onto Mercury within our numerical simulations is presented in Figure 6. The windfree examples agree reasonably well with Gladman & Coffey (2009), in that typical reaccretion timescales are 10 million years. Note that we underestimate the number of collisions arising from 4 km s −1 ejections (∼20 % versus ∼40% after 10 Myr). The reason for this discrepancy likely arises owing to the stochasticity inherent to the problem, together with the smaller number of ejected particles (110) used in the present study, as compared to the 7000 used in Gladman & Coffey (2009). For the purposes of the present study, a comparison between the cases with and without wind-induced drag is the primary focus. Figure 6. Number of particles, from an initial population of 110, that reimpact Mercury subsequent to being launched. Red and blue correspond to launch velocities of 4 km s −1 and 10 km s −1 respectively. Short-dashed lines correspond to drag associated with a wind of M 10   ("weak wind") and longdashed lines to winds of M 100   ("strong wind"). Even moderate winds significantly reduce the number of particles that reimpact Mercury's surface, allowing them to remain unbound from the planet.
The reaccretion rates occurring under the influence of windinduced drag were also presented in Figure 6. Both = M M 100    (strong wind) and even = M M 10    (weak wind) were sufficient to prevent most reaccretion, with the former allowing only one or two collisions upon Mercury, as opposed to ∼20−30 with no drag. These numerical results confirm the expectation from our analytics that drag timescales corresponding to the early solar wind prevent most of the reaccretion of material ejected from Mercury's surface.
In conclusion, within the scope of our analysis, the giant impact origin for Mercury's large core is consistent with the dynamical trajectory taken by material ejected from the primordial Mercury within the early solar system, permeated by the young Sun's enhanced wind.

Discussion and Conclusions
By way of a combination of semianalytical theory and numerical simulations, we have shown that the drag associated with the early Sun's wind was likely sufficient to cause orbital migration of centimeter-sized particles over <1 Myr timescales. Previous work has suggested that a giant impact upon Mercury's progenitor would expel its mantle into heliocentric orbit, where the material condenses into ∼ centimeter-sized objects. Accordingly, we demonstrated that the addition of the solar wind, which would have been significantly larger than today at the time of impact, reduces the likelihood that much of this ejected material would reaccrete onto Mercury.
Alternative solutions to the problem of reaccretion have been proposed. The Poynting-Robertson effect, for example, is most analogous to the influence of the solar wind (Burns et al. 1979). Previous work has shown that the timescale associated with Poynting-Robertson drag is indeed shorter than 10 Myr for nominal parameters. Criticism of the Poynting-Robertson drag as a solution has included the argument that the optical depth of the purported ring of debris may exceed unity (Gladman & Coffey 2009). This would increase the timescale of removal, allowing more time for reaccretion.
The benefit of the solar wind in this context is that it is potentially 10-100 times stronger, so even if the penetration of solar wind ions is somewhat reduced by the presence of a dense ring of material, there exist multiple orbital decay times within which to remove the ring. In addition, differential apsidal precession of ring particles would give the ring a nonzero scale height that would reduce the degree of self-shielding. The issue of self-shielding would be worth exploring in follow-up work.

Timing of Impact
Age constraints upon the final formation times of the inner terrestrial planets are largely based upon isotopic data of both Earth rocks and Martian meteorites (Raymond et al. 2018). In particular, the radioactive decay system of W-Hf suggests that Earth's formation time spanned 10 s of millions of years, with its final giant impact perhaps occurring upwards of 100 Myr after the beginning of planet formation (Kleine et al. 2009;Fischer & Nimmo 2018). Mars, in contrast, likely formed substantially earlier, within a few million years and before the dissipation of the protoplanetary nebula gas (Dauphas & Pourmand 2011).
Unfortunately, similar direct geochemical constraints are absent for Mercury (and Venus) owing to a lack of Mercurian meteoritic samples. The likely time of Mercury's mantle-stripping impact is therefore poorly constrained. Indirect inferences may be made using the Earth-Moon system, which likely formed from a similar giant impact (Canup & Asphaug 2001), hinting that Mercury's final giant impact may have occurred within a similar timeframe; within ∼100 Myr. In principle, a lower limit on the age of Mercury's surface might be deduced from crater size-frequency analysis (Strom et al. 2011). However, an absolute chronology for Mercury's cratering record is largely absent, such that it is difficult to rule out impacts that occurred within the first few million years after disk dissipation.
Uncertainty in the timing of the impact is not critical to our hypothesis in a qualitative sense; ejecta launched within the first ∼100 Myr is uniformly removed from Mercury's orbit. However, the decaying strength of the wind suggests that if the impact occurred earlier than ∼10 Myr, debris would be pushed outward to wider orbits than Mercury, potentially becoming part of Venus or Earth (see Figure 5). This contrasts with debris launched later, which is more likely to simply decay onto the Sun's surface. If the final impact occurred earlier than 10 Myr, it would suggest that Mercury's formation time is more in-line with Mars (Dauphas & Pourmand 2011). Mars is consistent with a planetary embryo, which had its growth halted, possibly owing to the inward migration of Jupiter (Walsh et al. 2011;Batygin & Laughlin 2015). Interestingly, Mercury's suspected initial mass, at 2.25 times its current value (Benz et al. 2007) is only about 10% larger than Mars' current mass. This raises the possibility that Mercury, too, might have been a stranded embryo, later altered during the phase of giant impacts (Raymond et al. 2006(Raymond et al. , 2018.
Ultimately, there is no specific reason to favor late or early impact scenarios for Mercury, owing to incomplete ageconstraints upon its surface. However, large uncertainties likewise exist in the early wind's properties, and so even with accurate ages for Mercury's impact it may still be difficult to say with great confidence whether the debris ended up in the Sun or in the other terrestrial planets. Nevertheless, our modeling here indicates that at all times earlier than 100 Myr, and subsequent to disk dissipation, the Sun's wind was likely strong enough to remove ejected material from Mercury's orbital vicinity and prevent substantial reaccretion.

Size Distribution of Ejecta
Throughout this work, we implicitly assumed that the majority of the mass of ejected particles existed as centimeter-sized spheres. This assumption was chosen in light of thermodynamic arguments presented in Benz et al. (2007). All timescales computed in our work depend linearly upon particle sizes, and thus it is important to briefly return to the robustness of our assumption of 1 cm particles. In particular, Benz et al. (2007) followed the thermodynamic evolution of material that is initially shocked to a supercritical fluid state, before subsequently cooling along an isentrope. As it cools, the ejecta eventually becomes subcritical, with the most energetic material initially condensing to a vapor, while less energetic material first reaches a liquid state (Johnson & Melosh 2012). A full treatment of the cooling process must apply different approaches for each of these thermodynamic outcomes. The work of Benz et al. (2007) used classical nucleation theory (Raizer 1960;Zel'Dovich & Raizer 2002) to model the highestenergy material, but the lower energy material required a different approach (Grady 1982;Melosh & Vickery 1991).
A more recent treatment of the condensation of ejected material was presented in Johnson & Melosh (2012), but within a different impact regime-that of the 10 km impactor that marked the end of the Cretaceous Period on Earth (Alvarez et al. 1980). Their model consisted of a hemispherical, expanding plume of vapor, out of which solid spherules condensed. The size distribution of condensed material was strongly peaked, with a mean value that scales approximately linearly with the impactor radius. Additionally, the average size of debris exhibited strong variation with respect to the impactor's velocity. Specifically, for an impactor of 1000 km in diameter, the size of debris ranges from roughly 0.01 cm at impact velocities of 15 km s −1 , rising to ∼10 cm for 25 km s −1 , but falling gradually again to 0.01 cm at speeds of 50 km s −1 .
In the case of Mercury, the impactor's diameter (assuming the single impact case Benz et al. 2007) likely exceeded 1000 km, perhaps approaching ∼6000 km (Chau et al. 2018). Furthermore, the most likely velocity of the impact has been suggested to lie within the range of 20-30 km s −1 (Benz et al. 2007;Chau et al. 2018). Given that a 10 km impactor produces debris with characteristic size s=350 μm, the larger impactor for Mercury will create correspondingly larger debris. For the expected range of impactor sizes, 1000-6000 km, and for an assumed impact speed of 30 km s −1 , we can scale up the results from Johnson & Melosh (2012) to estimate that the mean radius of debris would lie in the range s≈3−20 cm. These values are roughly consistent with our assumed size (1 cm), but larger debris sizes are possible, and will depend upon the impact velocity and impactor size. Despite the uncertainty in particle sizes, we may reframe the problem in terms of the largest particles s max that survive after a time τ subsequent to the giant impact. Using Equation (4) which suggests that even if the typical size of particles was closer to 10 cm, as our extrapolations from Johnson & Melosh (2012) suggest, wind strengths of M 100   remain sufficient for removal. It is unclear what fraction of the total mass of ejected material would reside in the largest particles, and thus it is worth noting that all ejecta smaller than s max would also be removed.
Observational constraints upon the remnants of giant impacts may become available in the form of debris disks around other stars , some of which are already interpreted as indicative of past giant impacts, for example, BD +20 307 (Weinberger et al. 2010) and HD 172555 (Lisse et al. 2009;). This latter system possesses a population of fine dust (s100 μm) in addition to a population of larger dust inferred to match Spitzer data. There are few direct probes available of centimeter-sized populations around analogs of the young solar system, and thus it remains difficult to validate the assumption of centimeter sizes using existing data.
In addition to particle size, it is essential to consider whether an appropriate mass of silicates is removed from the system to account for Mercury's composition. While a full treatment of this problem is beyond the scope of the present investigation, we note that at the high energies required to match present-day Mercury the majority of escaping material is likely to become vaporized (Benz et al. 2007;Chau et al. 2018). The debris considered in this study condense from this expanding vapor cloud (Johnson & Melosh 2012), but eventually the condensation sequence is quenched. Accordingly, an uncertain fraction of the vapor column is converted into solid debris-the rest remains in vapor form and is likely removed by photodissociation on a short timescale (potentially 1 yr . Accordingly, the influence of solar wind-induced loss of solids upon Mercury's composition depends critically upon the ratio of vapor relative to condensates. At least within the regime of a 10 km impactor, collision velocities of 30 km s −1 appear to lead to a roughly equal split in mass between condensates and vapor (Johnson & Melosh 2012), but it is unclear whether a 50% vapor fraction is appropriate for Mercury-like giant collisions. In order to determine this fraction with greater precision, it is necessary in the future to undertake simulations of the giant impact itself, while tracking the debris with an appropriate equation of state for shocked silica (Pierazzo et al. 1997;Kraus et al. 2012). In tandem, imposing thermodynamical evolution models such as those of Johnson & Melosh (2012) on top of SPH simulations (Chau et al. 2018), promises to improve estimates of the particle size distribution. If particle sizes derived from these more involved treatments exceed several meters in radius, then unrealistically large solar winds may be required to remove the debris, and alternative pathways toward removal of material would be required.

Implications for Exoplanets
Mercury's high density is anomalous in the solar system, but a growing number of exoplanetary density measurements are facilitating a comparison to exoplanetary analogs. While Mercury-sized planets have been detected around other stars (Barclay et al. 2013), few density estimates exist of this size class (Jontof-Hutter et al. 2016). The only current example, Kepler-138b, is a Mars-sized body with both mass and radius measurements, but does not show enhanced density and furthermore exists around an M-dwarf star (Jontof-Hutter et al. 2015).
If the stellar wind is instrumental in removing debris from giant impacts, we would make the qualitative prediction that closer-in rocky planets are capable of reaching densities that exceed those of more distant planets. This prediction follows from the expectation that subsequent giant impacts selectively remove silicates (Marcus et al. 2009;Bonomo et al. 2019), but only if the planet is close enough to the host star for the material to be removed. Such a trend has recently been revealed (Swain et al. 2019), where hotter terrestrial planets reach higher densities than cooler examples. Nevertheless, such a general prediction is confounded by the unknown degree to which stellar winds vary from system to system, and the enhanced impact speeds and reduced reaccretion times expected at shorter orbital periods (Volk & Gladman 2015).
Extrasolar material in the more exotic form of polluted white dwarfs (Jura & Young 2014) provide additional constraints upon the occurrence of Mercury-like planetary compositions. Specifically, white dwarfs present in their spectra evidence of rocky material recently and continuously falling onto their surfaces. It recently became possible to analyze the composition of this material and determine the Fe:Si ratios and oxygen fugacities they possess (Doyle et al. 2019). Although the parent bodies and dynamical histories of the white dwarf pollution remain unidentified, even among this set of cosmochemical measurements Mercury's oxygen fugacity and Fe:Si ratio stand as anomalies.

The Status of the Giant Impact Hypothesis
The giant impact scenario for Mercury's formation has existed for decades and taken numerous forms (Chapman 1988;Benz et al. 2007;Asphaug & Reufer 2014;Chau et al. 2018). In this paper, we have focused upon one specific concern with this model, and that is the likelihood that any material blasted off Mercury is likely to reaccrete over a short timescale. We have shown that the solar wind constitutes a viable pathway to removing this ejecta, and thereby alleviating this problem with the giant impact hypothesis.
The MESSENGER mission has afforded a higher level view of the evidence than merely the Si:Fe mass ratio (Solomon et al. 2018). In addition to Mercury's large iron content, the planet's surface was revealed to possess an abnormally low oxygen fugacity with respect to the other terrestrial planets (Nittler et al. 2011;McCubbin et al. 2012). These chemical constraints are difficult to compare with the giant impact framework, in part owing to computational limitations related to following individual chemical species within the SPH simulations typically applied to giant impact investigations (Lock & Stewart 2017;Chau et al. 2018), and furthermore due to the difficulty in addressing the degree of equilibration between silicates and iron within each body (Dahl & Stevenson 2010).
An additional revelation, that the surface of Mercury is volatile-rich, was once considered as inconsistent with a giant impact model (Peplowski & Evans 2011). The idea that giant impacts preferentially remove volatiles has typically been motivated by low volatile contents of the Moon , itself an outcome of a giant impact. However, the volatile depletion of the Moon remains difficult to replicate within giant impact simulations, requiring fractionation between a Moon-forming disk and the proto-Earth (Canup et al. 2015). In other words, the volatile depletion seen on the Moon occurred as a result of it being the remnant of the impact, forming out of circumplanetary material that experienced volatile fractionation (Lock et al. , 2018Lock & Stewart 2017;Nakajima & Stevenson 2018). With respect to Mercury's high volatile content, then, it is more appropriate to compare Mercury to the Earth and not to the Moon .
A final point regarding volatiles is that giant impacts constitute a robust expectation during the latter stages of terrestrial planet formation generally, across a range of models Raymond et al. 2018). Therefore, it is likely that Venus, Earth, and even Mars experienced giant collisions and yet exhibit similar volatile contents to Mercury (Ebel & Stewart 2018). Accordingly, the high volatile content of Mercury no longer appears inconsistent with the planet having experiencing a giant impact at the end of its formation epoch.
The feasibility of the giant impact hypothesis for the origin of Mercury's anomalously large core remains contentious. Nevertheless, the problem of reaccretion of ejected material becomes significantly reduced when the influence of the solar wind is taken into account. Thus, in early planetary systems, the stellar wind may play an instrumental, yet underappreciated role in the final architectures of planetary systems. Within this context, Mercury's anomalous density stands as the nearest, and most familiar example of the role of the Sun's wind in the early solar system. C.S. thanks Gregory Laughlin, Ravit Helled, Brett Gladman, and Sarah Stewart for useful discussions, along with the 51 Pegasi b Heising-Simons Foundation grant for their generous support. F.C.A. acknowledges support from the NASA Exoplanets Research Program. We acknowledge the thoughtful comments of two anonymous reviewers, whose input substantially improved the quality of the manuscript.