Evolution of Mercury's Earliest Atmosphere

MESSENGER observations suggest a magma ocean formed on proto-Mercury, during which evaporation of metals and outgassing of C- and H-bearing volatiles produced an early atmosphere. Atmospheric escape subsequently occurred by plasma heating, photoevaporation, Jeans escape, and photoionization. To quantify atmospheric loss, we combine constraints on the lifetime of surficial melt, melt composition, and atmospheric composition. Consideration of two initial Mercury sizes and four magma ocean compositions determine the atmospheric speciation at a given surface temperature. A coupled interior-atmosphere model determines the cooling rate and therefore the lifetime of surficial melt. Combining the melt lifetime and escape flux calculations provide estimates for the total mass loss from early Mercury. Loss rates by Jeans escape are negligible. Plasma heating and photoionization are limited by homopause diffusion rates of $\sim10^{6}$ kg/s. Loss by photoevaporation depends on the timing of Mercury formation and assumed heating efficiency and ranges from $\sim10^{6.6}$ to $\sim10^{9.6}$ kg/s. The material for photoevaporation is sourced from below the homopause and is therefore energy-limited rather than diffusion-limited. The timescale for efficient interior-atmosphere chemical exchange is less than ten thousand years. Therefore, escape processes only account for an equivalent loss of less than 2.3 km of crust ($0.3\%$ of Mercury's mass). Accordingly, $\leq0.02\%$ of the total mass of H$_2$O and Na is lost. Therefore, cumulative loss cannot significantly modify Mercury's bulk mantle composition during the magma ocean stage. Mercury's high core:mantle ratio and volatile-rich surface may instead reflect chemical variations in its building blocks resulting from its solar-proximal accretion environment.


INTRODUCTION
MESSENGER data from X-ray, gamma ray, and neutron spectrometers constrain the composition of Mercury's surface, and motivate theories and models to understand Mercury's bulk composition, formation, and evolution. The surface composition and geology of Mer-cury is compatible with partial melting of cumulates that were originally formed by magma ocean crystallization (McCoy et al. 2018). Subsequent impact excavation exposed the cumulates at the surface (McCoy et al. 2018;Charlier et al. 2013). The low oxygen fugacity (f O 2 ) of the uppermost layer of Mercury's regolith, together with Mercury's large core size, suggest a reduced mantle where nominally lithophile elements such as Ca, Mn, Cr, and Ti are present in sulfides rather than silicates (Vander Kaaden & McCubbin 2016). Relative to basaltic rocks exposed at the surface of other terres-trial planets, a large amount of the moderately volatile element Na (3-5 wt%) is detected on Mercury's surface (Peplowski et al. 2014). Observations of Na variation in Mercury's exosphere may relate to night-side deposit formation and dawn re-emission (e.g., Cassidy et al. 2016). Hence, it remains an open question how moderately volatile elements such as Na may have accumulated on the surface-whether from an extant or now extinct process-and how their abundance compares to Mercury's bulk composition.
Magma oceans are pivotal in determining the initial conditions and subsequent evolution and chemical differentiation of terrestrial planets in the solar system (e.g., Elkins-Tanton 2012; Chao et al. 2021). Radiometric dating reveals that magmatic iron meteorites, which represent planetesimal cores, formed within 2 Myr of solar system formation (Kruijer et al. 2014). The rocky planet whose mass is most similar to that of Mercury, and for which samples are available, Mars, likely accreted, formed an iron core, and underwent complete solidification of its magma ocean ocean within about 20 Myr of solar system formation (Bouvier et al. 2018). Crucially, rapid core formation in terrestrial planets requires a magma ocean to enable efficient metal segregation (Stevenson 1990). By analogy, and given its solar-proximal location, extensive melting is therefore expected to have occurred on proto-Mercury (Vander Kaaden & McCubbin 2016;Brown & Elkins-Tanton 2009). Following its crystallisation, partial melting of magma ocean cumulates have been invoked to explain Mercury's contemporary surface composition (McCoy et al. 2018).
Energy from accretion and radiogenic heat (e.g. from 26 Al) may have driven the differentiation of Mercury if it formed sufficiently early in solar system history (Bhatia & Sahijpal 2017;Siegfried & Solomon 1974). Following a phase of rapid growth, the subsequent reduction of impactor flux would have enabled Mercury's magma ocean to cool and crystallize without additional largescale remelting. During this time, the mantle is expected to have stratified into a basal layer of olivine and a plagioclase and clinopyroxene dominated crust which is now observed on Mercury's surface (Brown & Elkins-Tanton 2009). During the cooling of the magma ocean when the surface remains mostly molten, chemical species readily exchange between the interior, atmosphere, and exosphere-as occurred for other terrestrial planets in the solar system (e.g., Elkins-Tanton 2008). Fegley & Cameron (1987) addressed the hypothesis that the anomalously high bulk density of Mercury (owing to a high core/mantle ratio) is the result of evaporation of silicate melt components from the surface of a Hermean magma ocean. They presumed atmospheric loss was sufficiently slow that the atmosphere remained in equilibrium with the magma ocean. In their model, vapor was removed in a step-wise fashion and the composition of the magma ocean evolved accordingly. In reality, however, evaporated species are transported, mixed, and lost from the atmosphere and exosphere, with the flux at which loss occurs integrated over the magma ocean lifetime ultimately dictating the total mass loss. Therefore, consideration of interior, atmospheric, and exospheric processes are necessary to assess whether significant quantities of rock-derived atmospheres can be lost during the Hermean magma ocean stage.
Based on observations of solar-mass stars, the early solar extreme ultraviolet (EUV) and X-ray fluxes were likely 400 times larger than they are today. This would have made photoionization a highly efficient nonthermal, and photoevaporation a highly efficient thermal, atmospheric escape mechanism Tu et al. 2015). Other loss mechanisms of potential importance include atmospheric sputtering and kinetic escape (e.g. Jeans escape) that occur over the lifetime of the magma ocean. Non-thermal loss rates can be constrained by known plasma pressures at proto-Mercury due to the incoming solar wind, as well as EUV luminosities of the early Sun estimated from population studies of nearby Sun-like stars (Ribas et al. 2014;Tu et al. 2015).
In this paper we establish the extent of element evaporation and loss from Mercury during its early magma ocean phase. Models are constructed of the (1) coupled evolution of the magma ocean and atmosphere, (2) evaporation of metals and metal-oxide species from the Hermean magma ocean, (3) the mixing ratios and abundances of molecular species throughout the atmosphere and at the exobase, and finally (4) loss rates of these species from the upper atmosphere. We discuss these results in the context of the chemical evolution of Mercury's surface environment, bulk composition, and present-day observations.

Overview
Our combined modeling strategy provides insight into the initial composition and evolution of Mercury's exosphere by considering (1-2) energy and mass exchange between the interior and atmosphere, (3) speciation in the atmosphere, (4-5) loss from the exosphere: 1. SPIDER Bower et al. 2018) is a coupled interioratmosphere model used to determine the surface temperature and lifetime of melt at the surface, as well as the pressure-temperature structure of the atmosphere. Magma ocean cooling is regulated by the atmospheric opacity, which depends on the pressure of atmospheric species and hence the solubility of species in melt.
2. VapoRock calculates the equilibrium partial pressures of metal-bearing gas species of the elements Si, Al, Mg, Ca, Na, Fe and K above the magma ocean surface (Wolf et al. 2021). This determines the metal-bearing composition of the atmosphere as a function of temperature and the bulk composition of the magma ocean. It utilises ENKI's Ther-moEngine (http://enki-portal.org) and combines estimates for element activities in silicate melts with thermodynamic data for metal and metal oxide vapor species (Lamoreaux et al. 1987;Lamoreaux & Hildenbrand 1984).
3. VULCAN (Tsai et al. 2017(Tsai et al. , 2021 solves for the equilibrium chemistry of the atmosphere as a function of altitude by using element abundances for metals (output by VapoRock), volatile abundances (output by SPIDER), and the atmospheric pressure-temperature structure (also output from SPIDER). This provides the mixing ratios of atmospheric species, which are required to calculate escape at the exobase.
4. DISHOOM (Oza et al. 2019;Gebek & Oza 2020) is an atmospheric evolution model which computes the total mass loss of gaseous species to space due to ultraviolet (EUV) heating, surface heating, as well as plasma heating. borino et al. 2019;Vorburger et al. 2015;Wurz & Lammer 2003) determines the rate of exospheric escape of particles due to Jeans escape and photoionization. It tracks particle trajectories using a thermal energy distribution that depends on the temperature at the exobase.

Cooling of the magma ocean
Previous thermal modeling of Mercury's interior has focused either on the accretion phase (Bhatia & Sahijpal 2017) or its long-term evolution over billions of years (e.g., Tosi et al. 2013;Grott et al. 2011;Spohn 1991;Stevenson et al. 1983). Here, we model the thermal evolution of Mercury's magma ocean at the end of its accretion phase. At this time, the final magma ocean cools and crystallizes on a timescale short enough such that there is negligible disruption so our results remain independent of its accretion history. We model the thermal evolution of Mercury's magma ocean using SPIDER (Bower et al. 2018; to constrain the duration of melt at the surface as it cools from 2400 K to 1500 K. This is necessary to compute the evaporation of metals and metal oxides at the planetary surface, prior to the formation of a surface lid around 1500 K. Heating by the decay of radiogenic isotopes 26 Al, 40 K, 232 Th, 235 U, 238 U is included in our model and the model starts from solar system time zero to obtain an upper estimate of the surface cooling time. The main parameters are provided in Table A1 (Appendix A)and are guided by the parameters and results from previous models of Mercury (Bhatia & Sahijpal 2017;Tosi et al. 2013).
Cases prefixed by 'S' ('Small Mercury', Table 1) have a planetary radius of 2440 km, which is the present-day radius of Mercury. Cases prefixed by 'L' have a radius of 3290 km, which assumes Mercury was larger than at present day, perhaps due to mantle stripping driven by an impactor (Chau et al. 2018;Asphaug & Reufer 2014;Benz et al. 2008). Cases with 'V' (volatile) consider the partitioning of carbon and hydrogen species, here termed volatiles, between the melt and atmosphere as well as redox reactions . By contrast, cases with 'N' (non-volatile) do not consider volatiles but rather assume SiO is the only IR absorbing species; a suffix of '5' denotes a low SiO opacity (10 −5 m 2 kg −1 ) and a suffix of '3' denotes a large SiO opacity (10 −3 m 2 kg −1 ), both at a reference pressure of 3 × 10 −6 bar ( Fig. 2, Semenov et al. 2003). For non-volatile cases, the surface pressure of SiO is imposed in SPIDER as a function of surface temperature as determined by VapoRock calculations. Then the atmospheric opacity and hence magma ocean cooling rate can be determined.
For cases SV and LV, carbon and hydrogen can either exist as reduced (CO, H 2 ) or oxidised (CO 2 , H 2 O) species, where the f O 2 is constrained to one log unit below the iron-wustite buffer (∆IW = −1, IW-1 hereafter). This is marginally higher than the most recent estimates for the f O 2 of Mercury's mantle (Cartier & Wood 2019). For cases SV and LV, we determine the total H and C abundances by calculating the ppmw necessary for an Earth-size planet to produce a 100 bar CO 2 (i.e. Venus-like atmosphere) and 270 bar H 2 O (i.e. one Earth ocean mass) atmosphere. The abundances of H and C are equivalent to 330 ppmw of H 2 O and 120 ppmw of CO 2 , respectively. The mass of large Mercury's mantle is about a factor of 5 larger than small Mercury's mantle, resulting in a 5 times increase in the total volatile budget.  (2017), respectively. Second letter of the case name denotes with volatiles (V) and no volatiles (N). Non-volatile cases have an additional number of either 5 or 3, to denote small (10 −5 m 2 kg −1 ) and large (10 −3 m 2 kg −1 ) SiO opacity, respectively (Semenov et al. 2003).

Evaporation from the magma ocean
At the high surface temperatures that characterize a magma ocean (>1500 K), the partial pressures of the vapor species of the major rock-forming oxides (e.g. SiO 2 , NaO 0.5 , KO 0.5 ) can become significant (Sossi & Fegley 2018;Sossi et al. 2019;Visscher & Fegley 2013). Gasliquid equilibria for these elements are described by congruent evaporation, generalised as: where M is the metal, x is the oxidation state of the metal in its gaseous state, and n the number of electrons exchanged in the reaction. Both x and n are integer values that may be ≥ 0 or ≤ 0. At equilibrium, the partial pressure of any given metal or metal-oxide species in an ideal gas is given by: (2) where K (1) is the equilibrium constant of the reaction (Eq. 1), X is the mole fraction, and γ is the activity coefficient of the metal oxide melt species, M x+n O (x+n)/2 . Equilibrium constants involving 31 gas species (Table B2, Appendix B) are calculated according to their thermodynamic properties given in Lamoreaux et al. (1987); Lamoreaux & Hildenbrand (1984). Evident from Eq. 2 is that estimates for the composition of the silicate melt in addition to the activity coefficients of its constituent components are required to correctly predict partial pressures. To this end, likely compositions representative of Mercury's crust, mantle and possible precursors are shown in Table 2. The MELTS algorithm is used to estimate activity coefficients of melt oxide species (Ghiorso & Sack 1995). The f O 2 is constrained to lie one log unit below the IW buffer (IW-1), which is calculated according to O'Neill & Eggins (2002): where T is temperature and R the gas constant. These ingredients together comprise the VapoRock code and permit calculation of equilibrium partial pressures over a range of temperatures, f O 2 , and silicate melt compositions (Wolf et al. 2021).

Magma ocean composition
The composition of the early Hermean mantle is uncertain. To address how variability in the surface composition affects the evolved partial pressures of metal and metal oxide gas species, four compositions are investigated (Table 2): (1) Enstatite Chondrites (EH4), (2) Bencubbin Chondrites (CB), (3) Northern Smooth Plain (NSP) lava, and (4) NSP source. The first composition assumes that core-mantle differentiation was sluggish, with the composition of the magma ocean being approximated by Enstatite Chondrites (EH4), often cited as appropriate starting compositions for Mercury due to their high bulk iron content, strongly reduced nature, and the resemblance of partial melts thereof to Hermean surface compositions (Nittler et al. 2011). We assume that all FeO is extracted in the form of metallic iron to form Mercury's core, which results in high SiO 2 and MgO contents in the complementary silicate fraction.
The second composition is based on chondrules found in the Bencubbin-class of Carbonaceous Chondrites (CB), which best reproduce Mercury's surface composition based on MESSENGER data (McCoy et al. 2018;Brown & Elkins-Tanton 2009). Spectrometric measurements of Mercury's surface show a crust rich in Na and S and poor in Fe relative to other basaltic rocks (Nittler & Weider 2019). Presuming these abundances are representative of bulk Mercury, two other compositions are investigated; the Northern Smooth Plain lava (NSP melt) that represents a volatile-rich composition observed on the surface, and its inferred mantle source (NSP source) (Nittler & Weider 2019). These two compositions represent Mercury's crust and mantle, respectively. Although it is not anticipated that magma ocean crystallisation produced the NSP melt composition directly, it is included to define an end-member opposing the CB composition that is Na and K poor, and comparably Fe-rich (Table 2).  (Wiik 1956); CB chondrite chondrule data with bulk CB Na and K mass balanced for chondrules to fit bulk meteorite iron-silicate ratio (Lauretta et al. 2007;Weisberg et al. 2000Weisberg et al. , 1990; and northern smooth plains (NSP) composition for the lava and source (Nittler & Weider 2019;Namur et al. 2016) 2.5. Atmospheric structure The atmospheric structure is constrained by the temperature at the magma ocean-atmosphere interface, the planetary equilibrium temperature, and the atmospheric composition and pressure (Appendix E). For non-volatile cases, the calculated vapor pressures of Si, Na, K, Fe, Mg, Al, and Ca oxide species in equilibrium with the magma ocean ( Fig. 2) are used directly in the exospheric loss model (Section 2.6). SiO vapour pressures reported in Table 1 at T surf = 2000 K are not strongly affected by the magma ocean composition and range from 10 −4±0.1 bar. For cases SV and LV in which Mercury's atmosphere contains outgassed Hand C-bearing gases, the partial pressures of H 2 , H 2 O, CO, and CO 2 are calculated by SPIDER according to volatile solubility and f O 2 buffered by the magma ocean at IW-1. A modified version of the VULCAN code is then used to compute the equilibrium chemical speciation in the atmosphere that contains both the metalbearing gases and H and C volatiles (Tsai et al. 2017(Tsai et al. , 2021. VULCAN computes the atmospheric mixing ra-tios using the pressure-temperature (P -T ) structure of the atmosphere.
VULCAN by default includes about 300 reactions for C, H, O, and N to which we added reactions involving Si, Mg, Ca, Fe, Na and K to obtain their equilibrium speciation (Table C3, Appendix C). Table 3 shows the initial element-to-hydrogen ratios used in the VULCAN calculations. Surface vapor pressures of Ca and Al-bearing species did not exceed 10 −6 bar in the magma ocean temperature range investigated, and are thus excluded. For the remaining species we used a case dependent P -T profile from SPIDER at a surface temperature of 2000 K to determine the mixing ratio of species in the atmosphere as a function of altitude. The P -T profile may imply condensation of certain elements initially present in the vapor, which may rain out of the atmosphere prior to escape. To assess this possibility, Gibbs Free Energy minimisation of the atmospheric composition were preformed throughout the atmospheric column using Fact-Sage 7.3. (Bale et al. 2016).
SPIDER, and for volatile cases also VULCAN, provide descriptions of the atmospheric structure and composition needed to determine the altitude of the homopause and exobase. The homopause is the altitude at which molecular diffusion exceeds diffusion by eddies and thus separates the well-mixed lower atmosphere from the mass-separated upper atmosphere. The exobase is the altitude at which gas is loosely bound to the planet and is collisionless (Knudsen number = 1) resulting in efficient escape.

Homopause level and diffusion
To determine the homopause level z hom we require the particle density at the homopause n hom . For a steadystate homopause height, the molecular coefficient D ik is equal to the eddy diffusion coefficient K zz , allowing us to solve for the particle density n hom . The diffusion coefficient D ik (m 2 /s) within the homosphere is calculated for each major species using the Chapman-Enskog relation (Chapman & Cowling 1970). It determines the binary diffusion rate of a gaseous species i with mass m i within a gas of average mass m k : where, k B is the Boltzmann constant and T skin the absolute temperature at the homopause (skin temperature). Homopause pressures (P hom = n hom k B T skin ) for each species are thereby about 10 −6 bar for all cases. We approximate the intermolecular distance σ ik = σi+σ k 2 with Note-Ratios are based on SPIDER and VapoRock results for magma ocean compositions given in Table 2. We neglect Al and Ca as their vapor pressure do not exceed 10 −9 bar at 2000 K for any composition.
the radius of the species relative to the mean species diameter weighted by the mixing ratio. The dimensionless collision integral Ω ik is assumed to be unity. The eddy velocity is often approximated by the atmospheric species thermal speed v th , and the characteristic eddy length scale L eddy is approximated by the atmospheric scale height H (e.g. Atreya et al. 1986). Values for K(z) = v eddy L eddy that are calculated based on this assumption exceed the suggested eddy diffusion coefficient upper limit of 320 m 2 /s by several orders of magnitude (Vlasov & Kelley 2015). Hence we use this upper limit in the volatile cases to determine n homo , which is based on the energy dissipation rate within the Earth's atmosphere.
To compute the z hom and T hom for the volatile cases we determine the altitude at which the P -T profile reaches a number density n homo (K(z)). As the number density at the homopause only depends on K zz with the same order of magnitude for both volatile and non-volatile cases, n homo (K(z)) is approximately 10 18 − 10 19 at m −3 . For the non-volatile cases we do not obtain a P(z)-T profile from VULCAN for the non-volatile cases due to the absence of H based species. We therefore use the barometric formula with gravity as a function of height to compute z hom .

Exobase level
Due to the large difference in number density between the homopause and the exobase, the barometric formula is not applicable assuming an isothermal upper atmosphere with height-dependent gravity. We approximate the exobase height of early Mercury, which is subject to extensive loss, by finding an exobase height that results in a loss rate that is in equilibrium with the homopause diffusion rate. The loss from the exobase is proportional to the exobase height (increasing surface area) whereas diffusion from the homopause to the exobase is inversely proportional to the exobase height (decreasing density gradient). The exobase altitude z exo of each species is determined for all cases by setting the homopause diffusion rateṀ diff,i equal to the largest, diffusion limited mass loss rate of photoionizationṀ ion (Eq. 9, Section 2.6).
The homopause diffusion rateṀ diff in kg/s of a species i is obtained by multiplying the diffusion coefficient D ik with the species number density gradient, the species mixing ratio at the homopause n i /n hom , and the homopause surface area A hom : The number density of particles at the exobase n exo is necessary to determine the number density gradient between the homopause and exobase and ultimately z exo . For a single species i, the exobase is defined at a altitude at which the particle free path (λ col = 1/n exo σ col ) is equal to the exospheric scale height (H = k B T /mg), therefore (i.e., Gronoff et al. 2020): with the collision cross section σ i , the skin temperature T skin , and the acceleration of gravity g(h) at the exobase altitude z exo . In a multi-species atmosphere, each species has a specific mass and collision cross section, leading to a species-specific scale height and exobase density and altitude. The collision cross sections (CCSs) of each species are approximated as their respective atom or molecule size (Table D4, Appendix D). Typical values for n exo are around 10 12 − 10 13 atom m −3 which coincides with ∼ 10 −12 bar. The skin temperature used for determining z exo and z hom is derived from the atmosphere model (Appendix E): with the magma ocean surface temperature T surf , emissivity (depends on optical depth and hence atmospheric composition and pressure), and equilibrium temperature T ∞ .

Exospheric loss
The E-MC escape model focuses on Jeans escape, and photoionization and photodissociation to investigate the loss of proto-Mercury's exosphere. These mechanisms compete for importance; Jeans escape acts at high exospheric temperatures, whereas photoionization and photodissociation act at large solar EUV and X-ray fluxes present during early times, respectively. The E-MC model simulates escape by tracking ∼ 10 5 exospheric particles, with trajectories initiated at the exobase with a random angle and energy selected from a Maxwellian velocity distribution function (Vorburger et al. 2015). The initial energy of the exospheric particles depends on the exobase temperature, which decreases with time due to magma ocean and atmospheric cooling.
The loss processes in the E-MC model are calculated on a particle by particle basis. As soon as a particle reaches Mercury's Hill radius (Table 4), it is assumed to have escaped Mercury's gravitational attraction and is subsequently removed from the simulation. Another loss process is through interaction with photons. At each altitude step starting from the exobase and moving away from the planet, the E-MC model calculates the probability of a particle being photo-dissociated or photoionized. If the particle is photodissociated, the code calculates the corresponding trajectories of the fragments and assesses the chance of escaping the gravitational well and the potential for subsequent photoionization. Ionized particles are considered lost from the exosphere, assuming they are picked up by the electromagnetic forces of the solar wind plasma or Mercury's magnetospheric plasma. Photoionization and photodissociation rates are scaled for each dominant species using the EUV flux of the early Sun. The EUV flux and mass loss are dependent on the rotational evolution of the Sun, with a fast rotator being much more active than a slow rotator Tu et al. 2015). We consider a moderately fast rotating Sun, where the EUV luminosity L EUV (J/s) is: where t (Ma) is the time since the formation of the solar system. Typical values for the incident EUV fluxes at Mercury at 1 Ma and 5 Ma are thereby 10 3.0 and 10 2.2 J/s m −2 respectively. The loss is equal to the sum of the exospheric particles that have either been photoionized or lost through gravitational escape. The loss of a given species from the exosphere at a given time is then calculated using its mixing ratio at the exobase. For a given exobase temperature, the loss rate (Ṁ = dm/dt) from the exosphere by photoionizationṀ ion is: where A exo is the surface area of the exobase and 0 ≤ ξ i ≤ 1 is the fraction of lost particles of a species i with mass m i , thermal speed of v therm i , and particle flux leaving the exobase of ψ source The area of the exobase is equal to sum of the R P and z exo (Section 2.5). The total loss is determined by integrating the loss flux over the lifetime of surficial melt.

Atmospheric loss and surface evaporation
Atmospheric loss by thermal processes (photoevaporation) and non-thermal processes (plasma-heating) are determined using DISHOOM (Oza et al. 2019), using Eq. 10 and Eq. 11, respectively. Preliminary calculation of Jeans-escape using DISHOOM demonstrated negligible loss due to surface heating compared to all other mechanisms. The escape parameters appropriate to proto-Mercury are summarized in Table 4.
Irradiation from the impinging solar wind plasma and high energy photons may heat the atmosphere and drive escape at a level that is significantly larger than surface heating and photoionization. Upper-atmospheric heating (photoevaporation) is caused by incoming X-ray and EUV photons that deposit heat into a neutral medium via molecular absorption (Watson et al. 1981) or photoelectric heating (Murray-Clay et al. 2009). This expands the atmospheric envelope beyond the gravitational influence of the body (R H in Table 4). The heating can be estimated by energy-limited escape driven by EUV photons (Watson et al. 1981), which is a reasonable ap- proximation to thermally-driven hydrodynamic escape (Krenn et al. 2021;Volkov & Johnson 2013): where, z abs is the absorption altitude, generally taken to be 1.25 planetary radii for an outgassing atmosphere (e.g., Johnson et al. 2015) where the X-ray and EUV photons can absorb and thereby deposit heat into the atmospheric molecules. We note that we use 1.25 Mercury radii R P , as a conservative lower-limit in absorption altitude as the homopause situated at ∼ 1.4 R P represents an upper-limit. The efficiency at which the atmosphere is heated, η EU V , is uncertain so we use 10 −3 as a conservative lower estimate (Ito & Ikoma 2021) and 10 −1 as an upper limit (Mordasini 2020). Both efficiencies that we used were previously applied to atmospheres that use vastly different planet parameters but similar enough as they consider a metal oxide (non-volatile) or a H/He (volatile) atmosphere respectively. Hot Jupiter H/He envelopes as well as volcanic atmospheres suggest η EU V may be as large as 0.35 (Murray-Clay et al. 2009;Lellouch et al. 1992).
Mass loss due to plasma-heating is observed at Jupiter's moon Io and is fundamentally driven by plasma ram pressure and magnetic pressure interacting with the atmosphere (e.g., Johnson 1990). Therefore, we estimate the atmospheric loss from an impinging plasma on Mercury by scaling to the plasma pressures measured at the Galilean satellites (Johnson 2004 (2020), the mass loss of a species i by plasma heating at proto-Mercury is: where x i is the element fraction of the species i in the atmosphere, andṀ i,Io its atmospheric sputtering loss rate at Io. The total plasma pressureP , gravitational binding energyÛ , and ion velocityv ion of Mercury are expressed as non-dimensional values that are scaled to Io's corresponding values. The total plasma pressure is additive where P tot = P mag + P ram . For the calculations we use a magnetic pressure P mag = 1.7 nPa based on an estimation of Mercury's magnetic moment of 2.76 ×10 12 T m 3 at the magnetopause stand-off distance of 1.4 R P . The ram pressure due to the solar wind varies from ∼ 10-30 nPa (Korth et al. 2012), yielding a total pressure of ∼ 12-32 nPa. Evaporation from the magma ocean and atmospheric loss have to be equal to retain the atmospheric pressure and thus a steady state. The evaporation rate of a species is approximated by the Hertz-Knudsen-Langmuir equation. The evaporation rate of a species i with molar mass M i over the surface of Mercury in mol/s is given by: with the evaporation and condensation coefficients γ (set to unity for a liquid), surface pressure P s , equilibrium pressure P eq , and the Mercury radius R P . By setting the homopause diffusion rateṀ diff,i equal to the evaporation rateṀ evap,i , the equation is solved for the ratio of surface to equilibrium pressure, p i,s /p i,eq for each species at each temperature step. For homopause diffusion rates, the ratio lies > 0.99, therefore the atmosphere up to the homopause is considered to be in equilibrium.

Surficial melt lifetime and atmospheric structure
The surface temperature of the Hermean magma ocean cools from 2400 to 1500 K in around 400 -9000 years, depending on the planetary size (i.e. mantle mass) and efficiency of radiative energy loss to space (Fig. 1). The cooling rate is inversely proportional to R P since it depends on the ratio of the planetary surface area to mantle mass. Hence, a large Mercury takes longer to cool than a small Mercury for otherwise identical parameters.

Non-volatile cases
The cooling trajectory of non-volatile cases is characterised by two episodes. First, when the surface temperature is high (early time), the pressure and hence opacity of SiO is large and therefore cooling is slow. The second episode of cooling is rapid, since even a small decrease in surface temperature produces a drastic fall in both SiO pressure and opacity, driving the planet towards cooling like an ideal black body. Therefore, the cooling timescale for N5 cases is only marginally greater than for an ideal black body, which bounds the minimum cooling time to 470 years. Increasing SiO opacity by 2 orders of magnitude increases the minimum cooling time to around 1000 years (N3 cases), marginally affected by the chosen magma ocean composition.
Non-volatile cases at 2400 K have thin atmospheres of < 0.1 bar and comprise > 99.8% gaseous Na, SiO, Fe, K, and Mg. The major constituents are Na and SiO at high temperatures whereas the mixing ratio of SiO rapidly decreases below 2400 K (Fig. 2). The partial pressure of Mg behaves similar to SiO but does not exceed 1% of the atmospheric mixing ratio for any composition. The mixing ratios of Fe and K, however, are more variable because they depend on the assumed magma ocean composition and reach their highest mixing ratios of 9% and 2%, respectively, for CB and EH4 compositions with high FeO and K 2 O (Table 2, Fig. 2). Refractory components-AlO and Ca-have negligible partial pressures (< 10 −6 bar at 2400 K) and are thus ignored in further calculations. The highest total surface pressure of metal-bearing species at low temperatures is obtained with the NSP melt composition. At 2400 K, the total surface pressure is 6.16 × 10 −2 bar, which decreases to 8.52 × 10 −6 bar at 1500 K. Figure 3 shows the homopause levels of Na for the non-volatile cases as a function of time for CB and NSP melt compositions. Exobase levels lie within a few 100s of kilometers of the homopause and are omitted in the log-log plot as a result. During the magma ocean phase, the levels evolve within 470-660 years for N5 cases, and 1100-1480 years for N3 cases. The homopause and exobase locations are only weakly sensitive to the planet size and gravity. The magma ocean composition, however, exerts a strong influence on the atmospheric structure. For the CB case, the homopause lies at 685 km whereas for the NSP melt composition it lies around 1258 km at a magma ocean surface temperature of 2400 K. The early inflation of an atmosphere above a cooling magma ocean is due to increasing T skin caused by decreasing IR-opacity as the partial pressure of SiO decreases (Eq. 7). Following this stage, the homopause altitude falls to 83 km (CB) and 439 km (NSP melt) at 1500 K. The exobase density and location is further dependent on the mean cross section of atmospheric species, which is tied to the composition-dependent vapor pressures (Eq. 6). The high vapor pressure of Na in the NSP melt relative to the CB composition lowers the mean molecular weight and the mean collision cross section of the atmosphere, which both increase its extent.  Table 1 for case parameters.

Volatile cases
Volatile-bearing cases result in cooling times of 3400 years (Case SV) and 8900 years (Case LV). Both small and large proto-Mercury have the same initial volatile abundances of C and H by ppmw, but this manifests in a larger total reservoir size of volatiles for a large proto-Mercury compared to a small one. The mass of volatiles in the atmosphere defines the surface atmospheric pressure, which in turn determines the optical thickness of the atmosphere and hence the efficiency of radiative cooling.
Atmospheres of volatile cases around a small and large proto-Mercury reach surface pressures of about 5 and 12 bar at a magma ocean surface temperature T surf of 2000 K, respectively (Fig. 4). This result is independent of the partial pressures of the metal-bearing species, as their contribution is ≤ 0.1 bar at T surf = 2000 K. Thus, it is the outgassed hydrogen and carbon species (which depends on their solubilities) that dictates the surface pressure. VULCAN is then used to compute the equilibrium chemistry of the atmosphere accounting for the   (Table 4) Sodium hydride (NaH) is also present at the surface at pressures one order of magnitude lower than NaOH and remains about constant throughout the atmosphere, reaching similar mixing ratios to K. Potassium hydride (KH) is ignored as no rate constant exists in the NIST kinetics database (kinetics.nist.gov). Exobase levels are situated up to 2910 and 2590 km for small and large Mercury with homopause levels down to 2360 and 2160 km, respectively. The P -T structure of the atmosphere gives T skin values for the homopause and exobase of about 1021 and 893 K (Table 4), respectively, for small and large cases, which lies well below the T skin calculated for the non-volatile cases (T skin = 1615 K). The exobase levels of the volatile cases are thus comparable to the non-volatile cases. Unlike the non-volatile cases, the planet size has a large impact on the atmospheric structure, which is solely due to the difference in the total volatile reservoir (Section 3.1).

Atmospheric loss
We find that the major atmospheric escape mechanism is photoevaporation (Eq. 10) with a lower limit of photoevaporation constrained to 10 6.6 kg/s and an upper-limit of 10 9.6 kg/s for both non-volatile (Na) and volatile cases (H, C, and O). The upper limit is thereby roughly three orders of magnitude larger than the photoionization of the major atmospheric species. For the high heating efficiency (η EUV = 10 −1 ) case, loss rates become evaporation limited when reaching 1600 K as the surface to equilibrium pressure ratio approaches zero. For low heating efficiencies (η EUV = 10 −3 ), the ratio of surface to equilibrium pressure p i,s /p i,eq remains at >0.93 (Eq. 12). Photoevaporation as an approximation of thermally-driven hydrodynamic escape (Eq. 10) therefore expresses the highest uncertainty on the stability of the atmosphere. Figure 5 shows the integrated mass loss over the most extensive surficial melt lifetime of 8900 years. The photoevaporative erosion dR of the surface can be estimated by assuming mass conservation where: Using a mantle density of ρ mantle = 3.5 g/cm 3 and assuming a high EUV heating efficiency of 10 −1 at a large Mercury size allows for ∼ 2.3 km loss of crust over the 8900 years of the volatile case surficial melt lifetime. We have shown using Eq. 12 that this case becomes evaporation limited due to the high photoevaporation rates. The total integrated loss when assuming an EUV heating efficiency of 10 −1 is therefore lower than shown in Figure 5 but not significantly. This is due to most loss occurring during the early magma ocean stage, when the surface temperatures are high and evaporation is not limiting the potentially high photoevaporation rates.
Plasma-driven escape is diffusion-limited as a supply is required at the exobase (Eq. 5). TheṀ P calculated lie orders of magnitudes belowṀ ion as given in (Table  5) and therefore do not affect the atmospheric structure. Using an intermediate total pressure of 25 nPa, we find aboutṀ P = 10 3.4 kg/s for small and large, non-volatile Mercury cases and aboutṀ P = 10 3.1 kg/s  for volatile Mercury cases, respectively. The loss rates are thereby comparable to the plasma-driven escape observed on Jupiter's moon Io in SO 2 ∼ 10 3 kg/s (Thomas et al. 2004).
Time-averaged mass loss rates by photoionization are given in (Table 5). Like non-thermal plasma-driven escape, non-thermal escape due to photoionization at the exobase (Eq. 9) is diffusion-limited. As photoionization rates of non-volatile cases follow the same trends independent of planet size, we report results focusing on a small proto-Mercury only, omitting the large Mercury photoionization rates, which are mostly within a factor two for the dominant species (Table 5). The results of non-volatile N5 cases are also not reported, as an almost isothermal atmosphere results in less than a factor two larger loss rates at high temperatures.
Regarding the presence of metal oxide derived gaseous species, Gibbs Free Energy minimisation of the vapor phase (using FactSage) along a case-dependent atmospheric P -T profile ( Figure E1, Appendix E) indicates that Mg and SiO condense into clinopyroxene (1900 K) and then into olivine (≈ 1700 K) during cooling, by which temperature their fraction remaining in the gas is negligible. Iron persists in the vapor to lower temperatures, condensing partially into olivine before iron metal condenses at 1350 K. Therefore, while Mg, Si (and Ca and Al) all condense prior to reaching the exobase (T ≈ 1680 K), Fe is likely to partially reside in the vapor phase. Sodium never fully condenses (Nepheline, its major host mineral, condenses in very minor proportions below 1500 K), while K remains entirely in the vapor phase down to at least 950 K.

Loss from non-volatile atmosphereṡ
M ion in non-volatile cases is sensitive to the chosen initial composition (Fig. 6). The loss fluxes of the nonvolatile species of interest -SiO, Na and K-are proportional to their mixing ratios in the atmosphere (Fig. 2). In cases with high initial SiO partial pressures the mixing ratios and hence the diffusion-limited loss rates of Na and K increase during initial cooling as SiO becomes less abundant. This is most evident in the CB loss flux with an initially increasing loss rate despite decreasing temperatures (Fig. 6).
In the low-Na and low-K composition CB, where SiO is the dominant metal oxide at high tempera-

Photoevaporation -integrated loss
Case LV LN3 LN5 Figure 5. Integrated mass loss rates of the bulk atmosphere through photoevaporation from 2400 K to 1500 K for large Mercury cases. The shaded areas represent the uncertainty of loss rates based on initial conditions. For the N3 case only the upper and lower limit are shown as dashed lines as it lies within both V and N5 areas. The highest loss rates assume upper limits for EUV luminosity only 1 Ma after Sun formation, LEUV(1 Ma), and heating efficiency of 10 −1 , and the lowest assume LEUV(5 Ma), and heating efficiency of 10 −3 . Using a mantle density of ρ mantle = 3.5 g/cm 3 results in a maximum of crust material being lost to space ranging between 2.3 km and 16 cm, depending on the degree of XUV intensity and heating efficiency.
tures, loss fluxes of SiO reach up to 9.6 × 10 5 kg/s at T surface = 2400 K. The vapor pressure of SiO declines with respect to other dominant gas species (Na and K), thereby reducing its mixing ratio rapidly with decreasing temperature. The lower mixing ratio of SiO causes loss rates to drop to 9.4 × 10 2 kg/s at 1500 K. For the same composition, loss rates for Na and K are around 2.0 × 10 5 kg/s and 3.1×10 3 kg/s at T surf = 2400 K, increasing to 3.2 × 10 5 kg/s and 4.2 × 10 3 kg/s at T surf = 1500 K, respectively. For the high-Na end-member composition of NSP melt, the diffusion-limited loss rates for Na and K are 5.9×10 5 kg/s and 1.2×10 4 kg/s respectively when at T surf = 2400 K, and decrease to 3.9 × 10 5 kg/s and 4.3 × 10 3 kg/s respectively when T surf = 1500 K. The ratio of Na to K diffusion rates is ∼ 100 and hence about one order of magnitude higher than their mixing ratios in the atmosphere. Integrated diffusion-limited losses over the surficial melt lifetimes are shown in Figure 7. If we assume based on FactSage results, that SiO, and Mg are absent and only Na, K and Fe remain in the atmosphere, the diffusion-limited loss rates for all cases peak in the small Mercury NSPm case with a fairly temperature independent Na loss rate of about 10 6 kg/s. Loss rates for K and Fe are thereby about three orders of magnitude lower than that of Na and do not significantly contribute to the total loss. The difference in change ofṀ diff and thusṀ ion between species with continued magma ocean cooling results from the constantly dropping total surface pressure simultaneously to shifting partial pressures. At low magma ocean surface temperatures, Na exerts most of the pressure, whereas at high temperatures condensing species, such as SiO, are the predominant contributors to the total pressure and mean molecular mass of the atmosphere (Fig. 2). A lower surface pressure results in a lower homopause level and leads to a smaller diffusion area and rate, however the rapid decline of the SiO partial pressure in the atmosphere (n SiO /n hom ) results in the sharp drop of SiO homopause diffusion rates (i.e. Eq 5) but increases the T skin and therefore the extent of the atmosphere (Fig. 3). Similarly, the loss rates through ionization at the exobase for a relatively low mean molecular weight ('light') Na, K, and Fe based atmosphere are up to a factor five higher than the 2 × 10 5 kg/s loss of Na from a 'heavy' atmosphere, which includes Mg and notably SiO.
These limits are significant for atmospheric escape estimates by plasma heating, photoionization or Jeans escape. All three of those processes are calculated from the exobase, which is for plasma heating and photoionization where ions and photons can access a rarefied neutral atmosphere. This is therefore diffusion-limited as a high flux is required to source the neutral species experiencing a momentum transfer from the plasma. Photoevaporation is not necessarily diffusion-limited so long as a sufficiently large column exists at the altitude where EUV photons are able to absorb on to infrared emitting molecules (z abs , Fig. 8). In N 2 /CH 4 atmospheres (e.g. Kuiper Belt Objects) the critical column density is estimated to be 10 18 /cm 2 (Johnson et al. 2015) which is easily achieved at a fiducial absorption altitude of z abs situated at ∼ 1.25 R P , where the column density is equivalent to ∼10 21 /cm 2 for an isothermal scale height of H ≈ 150 km at a T surf = 2000 K. For a magma-silicate atmosphere, as studied here, SiO or a similar species would be able to re-emit in the infrared resulting in upper atmospheric expansion, and Roche-lobe overflow to space.

Loss from volatile atmosphere
In the volatile cases, assuming a speciation as encountered at the H homopause, the diffusion limited loss fluxes of primary species lie between 10 4 kg/s and 10 5 kg/s for all major species (H 2 , CO, H 2 O and CO 2 ) in the small and large proto-Mercury cases.
The loss flux of minor species Na, and K are several orders of magnitude lower than those of the major species, at 10 1 , and 10 0 kg/s respectively. In the non-volatile cases, loss fluxes of Na and K are directly proportional to their thermodynamic activities in the melt. Sodium activity increases by about a factor 4.5 from the CB to NSP melt composition and K by a factor two between NSP source and EH4, respectively. Relative to non-volatile cases, loss fluxes for Na and K are several orders of magnitude lower in the high pressure, volatilerich atmosphere at T surf = 2000 K. If we assume the loss fluxes of the dominant H, C, and O based species at T surf = 2000 K to be constant over the lifetime of the molten surface (Section 3.1) and integrate them for small and large Mercury volatile cases we obtain a total mass loss by photoionization of 4.1 × 10 16 kg and 1.8 × 10 17 kg, respectively. This exceeds the total photoionization mass loss from the low absorbing N5 non-volatile case by only about one order of magnitude (Fig. 7). The mass loss of Na in the volatile cases, however, only contributes about 10 12 kg of the total, which is about four orders of magnitude below the total mass loss of the non-volatile cases. Again, this assumes that the Na loss flux is constant in the volatile case. This is deemed appropriate, because Na is only a minor component of such atmospheres, is lost at slow rates that represent an insignificant fraction of its total budget, and does not condense before reaching T surf = 1500 K.  Table 5 tabulates the total atmospheric loss rates due to the following escape mechanisms: ionizationṀ ion , photoevaporationṀ U , and plasma-heatingṀ P . Figure 8 illustrates the atmospheric level from where the degassed magma ocean atmosphere is escaping.
We find that the loss fluxes from the exobase caused by photoionization and atmospheric sputtering are supplylimited (in all volatile and non-volatile cases) by homopause diffusionṀ diff , which dictates the exobase elevation in order to remain in steady state. The maximuṁ M ion of volatile and non-volatile cases are comparable, even though their atmospheres are comprised of different major species (non-volatile case: Na, and SiO, volatile case: H 2 , H 2 O, CO and CO 2 ). This similarity is attributed to the higher T skin of the non-volatile cases at T surf = 2000 K, caused by IR opacity which is tied to the mixing ratio of SiO. SiO pressure rapidly decreases with decreasing temperatures, which is contrary to CO 2 and H 2 O in the volatile cases. A higher skin temperature in the non-volatile cases hence compensates for the lower atmospheric pressures.
The photoevaporation rate is limited by the degree of upper atmospheric heating efficiency, η EU V , rather than  Figure 7. Integrated mass loss rates of exospheric oxides and elements through photoionization as the magma ocean cools from a surface temperature of 2400 K to 1500 K. SiO loss is shown although it is unlikely to persist in the upper atmosphere as it is consistently below its highest condensation temperature of T = 1900.
byṀ diff , as well as by the supply of gases from surface evaporation (eq. 12). The assumption of a constant photoevaporation rate is thus only valid for evaporation at high temperatures above 1600 K at which evaporation rates are fast enough for supply to be sustained, or for moderate mass loss rates of about 10 7 kg/s. The majority of mass loss occurs at high temperatures ( Fig.  5) when surface evaporation rates are high compared to photoevaporative loss rates, the later of which are independent of temperature and depend instead on the EUV flux.
The diffusion-limited loss rates by photoionization of the four major volatile species: H 2 , H 2 O, CO and CO 2 from a volatile-rich atmosphere total ∼ 10 5 kg/s. The loss of Na from a thick, volatile rich atmosphere is inhibited by its low mixing ratios at the homopause and exobase. Therefore, diffusion-limited loss of Na is most efficient when the atmosphere is thin, reaching a few 10 5 kg/s, which coincides with the total mass loss rates from volatile cases. The total integrated mass loss by photoionization from Mercury's exosphere is low for small volatile and non-volatile cases with ≤ 4.1 × 10 16 kg and ≤ 3.0 × 10 16 kg, respectively.
The mass loss of single species is negligible compared to the total inventory of the magma ocean reservoir. For example, 0.033 wt% H 2 O and a low estimate of 0.1 wt% Na in a total mass of M MO ≈ ×10 23 kg corresponds to a reduction of the total H 2 O reservoir mass (volatile cases) and Na (non-volatile cases) by ≤ 0.02%. Assuming a well-mixed mantle reservoir, the bulk composition of Mercury would not significantly change even for species with low abundance in the reservoir and large loss rates such as H 2 O, CO 2 , and Na. Energy-limited escape via photoevaporation, however, can erode up to 2.3 km of Mercury's crust which coincides with 0.3% of small Mercury ( Figure 5). Assuming small heating efficiencies as well as a lower EUV flux, leads to integrated photoevaporation losses and eroded crust thicknesses that are reduced by four orders of magnitude.
Physical segregation between crystal and liquid during magma ocean cooling will induce chemical fractionation of element abundances with respect to those of the bulk mantle. Namely, the incompatible lithophile elements (Na, K, Al, and Ca) become enriched in latestage liquids of a Hermean magma ocean. This effect is simulated by considering the composition of the NSP melt as a surface magma ocean analogue relative to that of its inferred source. These differences notwithstanding, the partial pressures of metal-bearing gas species vary only marginally among EH4, CB and NSP compositions. This is due to two factors: 1) vapor pressures of different elements vary by orders of magnitude among one another (e.g. compare Na with AlO) whereas abundances of these major elements vary only by a factor of 2-3 in most cases; and 2) higher mole fractions of Na and K in the NSP melt are partially compensated  Figure 8. Mass loss processes and their rates demonstrate the coupling between various atmospheric layers. RP is the planet radius, zexo is the exobase altitude, and z hom is the homopause altitude which governs exospheric loss processes of Jeans escape, plasma heatingṀP , and photo-ionizationṀion. z abs is the absorption altitude where upper-atmospheric heating (photoevaporation)ṀU commences. The absorption altitude is assumed to lie below the homopause and therefore photoevaporation is not limited by homopause diffusion. by their lower activity coefficients relative to the NSP source or EH4 composition. As we show in Section 3, all elements other than Na and K, and potentially Fe condense before reaching the exobase. We can therefore conclude that the atmospheric pressure and speciation around proto-Mercury only depends on the abundance of volatile-and moderately volatile elements.

Early origin of surface Na
In order to determine the potential impact of a magma ocean-generated atmosphere on the surface composition of a small proto-Mercury, we calculate the total mass of Na in the atmosphere. We consider a hypothetical scenario in which the atmosphere collapses as soon as the first crust forms, coinciding with the termination of the magma ocean stage when a surface temperature of 1500 K is reached. To obtain a result which is consistent with the notion of Na-poor building blocks (e.g., Humayun & Cassen 2000), we use the low-Na CB composition (Table 2) and its H-, and C-absent, pure Na atmosphere composition. The resulting total amount of Na integrated over the whole atmosphere yields about 10 11 kg.
Whether it precipitates as Na metal or as another compound depends on the composition of the atmosphere that exists. Although not considered in our model, such an atmosphere would contain significant quantities of other moderately volatile elements that could combine with Na to form complex molecules, ηEUV(10 −1 ) S 9.5 8.6 L 9.6 8.7 Note-Loss rates of plasma heatingṀP , pho-toionizationṀion, and photoevaporationṀU . Photoevaporation is insensitive to the atmosphere's emissivity and composition but depends on the EUV flux and the EUV heating efficiency (end members of 10 −3 and 10 −1 ). The EUV flux is a function of the age of the solar system (Eq. 8, after Johnstone et al. 2015;Tu et al. 2015).
namely, F, Cl and S. The species NaCl is inferred to be stable among volcanic gases (e.g., Aiuppa et al. 2003) and has been directly observed in Io's atmosphere (Lellouch et al. 2003;Moullet et al. 2010), and may therefore be a potential candidate to form surficial deposits. Sodium chloride is also observed as a stable precipitate from experimentally-generated volcanic gas analogues (Renggli & Klemme 2020), and is therefore likely to occur as a Na-bearing phase on the Mercurian surface. This is supported by the coinciding distribution pattern of Na and Cl from more recent volcanic deposits found in gamma ray spectrometer data (Evans et al. 2015).
Here we consider a simplified case, for which the mass of sodium in the atmosphere is uniformly distributed over the surface of Mercury as pure, low density, Na metal, results in a layer less than 1 mm for the CB case. Using a more Na-rich composition like EH4 combined with an increased atmospheric reservoir size of large proto-Mercury would lead to a factor four thicker Na layer, but still less than 1 mm. In the volatile cases (i.e., with CO 2 and H 2 O), the amount of Na in the at-mosphere is identical to the volatile-free cases, as, in our model, the partial pressure of Na is independent of the presence of volatiles. The small dissolved quantities of CO 2 and H 2 O in the silicate melt (≤ 1000 ppm) should thereby not influence the activity coefficients of the major rock-forming species. This hypothetical Na metal layer would not outlast meteorite impacts, which are assumed to have removed 50 m to 10 km of early crust (Hyodo et al. 2021). For the enrichment to be preserved, the atmospheric sodium would have to be incorporated into a layer with a thickness exceeding the removed crust. However, for a minimum layer of 50 m we obtain a total Na wt% increase of merely ∼1 ppm and ∼10 ppb for small and large proto-Mercury cases, respectively. We thus conclude, that the collapse of an early Na-rich atmosphere would not contribute to a notable increase of Na in the surface.

Controls on mass loss
The mean column density at the exobase depends on the weighted average of the dominant species' collision cross sections (CCSs). Loss rates are directly related to the exobase density. However, using CCSs from Kim & Desclaux (2002) that are about one order of magnitude smaller would reduce the homopause and exobase levels by a few tens of kilometers and decreases the homopause diffusion limited loss by ≤ 2%. The sensitivity of mass loss to the chosen CCSs is therefore weak. In non-volatile cases, if we consider that all species except Na, K, and Fe condense (FactSage in Section 3), then mean molecular mass and the collision cross section of the atmosphere decreases which enhances molecular diffusion (Eq. 4). This pushes the homopause and therefore the exobase further from the planet surface, increasing the atmospheric surface area and therefore loss. Furthermore, the absence of SiO leads to a hotter skin temperature as the atmosphere becomes IR transparent, further enhancing loss. The difference of the ionization mass loss rate at the exobase between a Na, K, and Fe atmosphere and an atmosphere where SiO is a major component at high temperature is thereby about a factor three larger for all cases.
For rocky exoplanets on short orbits, the atmospheric temperature around our homopause levels (10 −7 bar) can be as high as 3800 K for a surface temperature of 2400 K (Ito et al. 2015). Mercury posesses different planet parameters (1 M Earth and 0.02 AU vs. 0.055 M Earth and 0.3 AU for Mercury), however, the more intense early UV flux experienced by Mercury could similarly boost the temperature at the homopause. Calculations with an increased skin temper-ature of 3800 K at the homopause resulted in about a factor two higher photoionization loss rates for all cases.
We used photoevaporation as a proxy for thermally driven hydrodynamic escape. Krenn et al. (2021) has shown for a large range of parameters that photoevaporation can underestimate hydrodynamic escape especially at low EUV fluxes. Given our large incident EUV fluxes of about 10 2 − 10 3 Js −1 m −2 and our escape parameters (Table 4) we expect to be within one order of magnitude of thermally-driven hydrodynamic escape rates (compare EUV fluxes and escape parameters to Fig. 4 in Krenn et al. 2021, although our escape parameters for small Mercury lie just below the shown range). Figure 8 illustrates the transport of mass away from different levels in the atmosphere. In our model, atmospheric escape can be either energy-limited (e.g.Ṁ U ) or diffusion-limited (Jeans escape,Ṁ P ,Ṁ ion ). Below we describe the role of enhanced atmospheric heating or cooling on diffusion and energy-limited escape.

Atmospheric evolution and structure
Diffusion rates are tied to the homopause density, which determines the homopause altitude. The eddy diffusion coefficient (K zz ) needed to determine n hom bears large uncertainties, however. For all cases, a K zz larger than the Earth-derived upper limit of 3.2 × 10 6 cm 2 /s would most likely be adequate to accommodate proto-Mercury's increased atmospheric temperature, increasing z hom and lowering n hom , leading to a slightly larger diffusion and therefore loss rate. Even if we assume a larger K zz , however, homopause diffusion will remain the limiting factor for mass loss. We find for volatile cases, that even if K zz is three orders of magnitude larger, the total loss for volatile cases increases by a factor of less than two. The sensitivity ofṀ diff and thereforeṀ ion to the eddy diffusion coefficient is therefore weak.
Ionization could further increase the exobase temperature, and hence the reported diffusion-limited loss fluxes could be a lower limit. Whether mass loss occurs from the exobase surface, or whether it is the result of an advective outflow, is canonically assessed by the escape parameter (e.g. Genda & Abe 2003). If the escape parameter λ 0 ≤ 3 the atmosphere experiences mass outflow due to its non-zero net velocity, and escape occurs inward of the exobase at the sonic point (where the thermal velocity exceeds the sound speed). If λ 0 3 the atmosphere escapes because the mean free path is longer than the scale height, and Jeans escape prevails. Table  4 shows how the exospheric escape parameters are all 4 ≤ λ 0 ≤ 15, a near-transitional escape regime between Jeans-and hydrodynamic end members, which was re-cently determined to be relevant for the putative magma ocean on the Moon (Tucker et al. 2021). These authors demonstrated via DSMC simulations (Bird 1994) that cooling due to escape is important for λ 0 15. Therefore in Table 4 based on our escape parameters, it would appear that although ionization may further enhance escape, cooling may temper this loss. In addition, the significant ionization rates ofṀ ion 10 6 kg/s at the semi-major axis of proto-Mercury promotes the generation of an ionosphere that is modulated by the planetary magnetic field. Simulations on an early Mars analogue have demonstrated that ion escape is efficient at removing material (Egan et al. 2019). Therefore, it is possible that we are underestimating escape by not considering magnetic interactions.
In the concurrent 'energy-limited' regime it would appear that if EUV photons are able to absorb on to a sufficiently high flux of molecules (Section 3.2.1), heating would overwhelm cooling. However, based on the escape parameters in Table 4 it appears that cooling associated with escape may be important, arresting loss. At the same time, the study of low-mass, close-in exoplanets orbiting Sun-like stars has posited the idea that low-mass planets are nevertheless born with hydrogen/helium (H/He) envelopes although these are rapidly lost due to photoevaporation (Mordasini 2020). For a H/He envelope equivalent to 1% the mass of proto-Mercury, we find that our upper-limit on photoevaporation results in the dissipation of a H/He envelope in 10 4.4 years, which is larger than the lifetime of the molten surface. The possibility of a H/He envelope to persist during the molten surface lifetime is therefore non-trivial, and could result in significant heating which could not only enhance escape but also elongate the melt lifetime past the 10 4 years we study here.

Origin and evolution of Mercury
The elevated core:mantle ratio, coupled with an Naand S-rich surface, distinguish Mercury from the other terrestrial planets. Two key hypotheses exist to account for these characteristics; (1) the preferential loss of silicate material, either by evaporation (Fegley & Cameron 1987) or by collisional stripping (Benz et al. 1988) and (2) equilibrium condensation and sorting of metal from silicate in the solar nebula (Lewis 1972;Weidenschilling 1978). In evaluating hypothesis (1), Fegley & Cameron (1987) concluded that ∼ 75-79% of silicate material would need to be lost during a fractional vaporisation hypothesis to reproduce the core:mantle ratio of present-day Mercury. In this work, we show that such high fractions of loss of silicate material are untenable, be it from a small or a large proto-Mercury (total mass losses are below 0.3%). The principal reason is that atmospheric cooling timescales are too rapid with respect to evaporation and escape timescales, meaning that integrated loss rates over ∼ 10 4 yr are small with respect to the mass of proto-Mercury. Moreover, substantial amounts of atmospheric-or collisional escape of Mercury's crust is not represented in the high K/U ratio of its surface (McCubbin et al. 2012), as preferential loss of silicate material will predominantly deplete its incompatible lithophile element budget (O'Neill & Palme 2008).
There are several caveats to our conclusions, namely, that our results are valid for dry or C-, and H-bearing atmospheres, but do not consider the effect of other minor volatiles (Cl, S, F) on the volatility behaviour of metals. Metal chlorides and metal sulfides may be important gaseous species under moderate temperatures (∼ 1000 K, Renggli et al. 2017), increasing their volatility. Secondly, conditions on the surface of Mercury may have been considerably more reduced than modelled herein (IW-5; Cartier & Wood 2019). Because the partial pressures of most metal-bearing species increase with decreasing f O 2 (Eq. 1) vaporisation rates for alkali metals may be an order of magnitude higher (considering that their exponent, n = 1/4; Sossi et al. 2019).
However, these faster evaporation rates may be offset by the presence of a surficial graphite layer on the magma ocean (Keppler & Golabek 2019). Such a layer is promoted under reducing conditions as the solubility of C in silicate melt decreases from ∼ 360 ppm at the IW buffer to 1 ppm at IW-4 (Duncan et al. 2017;Keppler & Golabek 2019). The extent of a graphite layer therefore depends on the C content of Mercury and its f O 2 , both of which are poorly known. A surficial lid would additionally delay cooling of the mantle, unless the lid is regularly broken as possibly occurred for the flotation crust on the moon (Perera et al. 2018). Nevertheless, the net effect of a graphite lid on Mercury's magma ocean would be to reduce the extent of degassing calculated herein. Therefore, we conclude that the physico-chemical characteristics of Mercury cannot have been produced during a magma ocean stage on a near fully-grown planet.
These obstacles are ameliorated when considering vapour loss from planetary building blocks. Should Mercury have accreted from smaller, km-size planetesimals, then melting and vaporisation on the precursor bodies would have led to more efficient mass loss (e.g., Hin et al. 2017). Thus, vaporisation may still be a physically viable mechanism to explain Mercury's composition, provided it occurred on its precursor bodies. However, another problem arises because moderately volatile elements, such as Na, S and K, are always more volatile (i.e., their partial pressures are higher for a given activity) than the major mantle components, such as Mg and Si (Sossi et al. 2019). Moreover, as demonstrated herein, Na is more easily lost with respect to Mg and Si due to its lower molar mass and higher tendency to remain in the gas phase in an adiabatically-expanding atmosphere (3.2.1). As such, appealing to evaporative loss of Mg and Si to increase the core:mantle ratio while retaining Na and K is inconsistent with evaporation from a silicate melt on small planetary bodies. Therefore, other hypotheses should be considered.

CONCLUSIONS
We combined chemical and thermodynamic equilibrium models of the thermal evolution of Mercury's magma ocean and gaseous species derived thereof, to model the thermochemical evolution of an early atmosphere on Mercury. For an initially large Hermean mantle with initial C and H budgets comparable to those of other rocky planets, namely Earth ('volatile cases'), the lifetime of surficial melt may have reached almost 10 4 years. Compared to a present-day sized proto-Mercury without a greenhouse atmosphere, this lifetime is an order of magnitude larger and thereby may enable early atmospheric mass loss to occur over an extended duration. Cases with C and H show that Mercury could have started with a 5-12 bar atmosphere. By contrast, excluding the presence of C and H species results in a thin, short-lived metal-and metal oxide-bearing atmosphere. The upper atmospheres of volatile cases are dominated by H 2 , and CO whereas non-volatile cases are mostly Na and SiO.
Photoionization is a minor exospheric loss mechanism, limited by homopause diffusion (Ṁ diff ) up to a maximum of a few 10 5 kg/s. If C and H volatiles are absent from the atmosphere, theṀ diff limit applies to SiO and Na. Mass loss rates via photoevaporation, M U ≤10 9.5 kg/s, exceed those from all other known mechanisms due to the high EUV luminosity of the early Sun. This could in the best case scenario erode an equivalent thickness of up to 2.3 km of proto-Mercury's crust when assuming high EUV heating efficiencies of 10 −1 . Atmospheric sputteringṀ U ∼ 10 3.4 kg/s (also limited byṀ diff ) occurs at the exobase, knocking-off neutral gas molecules due to the ram pressure of the solar wind.
By integrating atmospheric loss rates over surficial melt lifetimes, we bracket the expected total mass loss from Mercury's early atmosphere. Based on photoionization, Jeans escape, and plasma heating, the evaporation and loss of the magma ocean of proto-Mercury did not significantly modify its bulk composition. This is because magma ocean cooling times are too short to drive substantial total loss for the determined atmospheric loss fluxes. Photoevaporation can remove an equivalent crustal thickness of up to 2.3 km in about 10,000 years, which is approximately ∼ 10 20 kg of material. Integrated losses of even the most volatile elements considered here, Na and K, are insignificant with respect to their total budgets when escape is diffusionlimited (≤ 0.02% decrease of the initial Na composition which would be a difference of 3×10 −4 wt%). Hence the present Na-rich surface composition may indicate that catastrophic volatile loss during the magma ocean stage did not occur, and that Mercury's peculiar composition is inherited from that of the solar-proximal region of the nebula from which it accreted. The evolving surface temperature of the Hermean magma ocean is calculated using the SPIDER code, which is described in detail in Bower et al. (2018); ; Bower et al. (2021). Table A1 shows the parameters used to model proto-Mercury. The mass absorption coefficients of H and C volatile species are determined at 1.01 bar, and the coefficients of SiO at 3 × 10 −6 bar. 0.993 -H2 mass absorption (CIA) 5 × 10 −5 m 2 kg −1 H2 solubility law * H2O mass absorption 10 −2 m 2 kg −1 H2O solubility law * * CO mass absorption 10 −5 m 2 kg −1 CO solubility law * CO2 mass absorption 10 −4 m 2 kg −1 CO2 solubility law * * SiO mass absorption (large) 10 −3 m 2 kg −1 SiO mass absorption (small) 10 −5 m 2 kg −1 Initial surface temp 2400 ‡ K † † Average Al abundance based on the composition of EH4 and NSP source (Table 2). † Average current estimates for bulk heat source from Tosi et al. (2013) and natural abundances from Ruedas (2017). ‡ Similar to maximum temperature estimate of Mercury's surface during accretion and differentiation (Bhatia & Sahijpal 2017).

B. VAPOROCK SPECIES
The species included in VapoRock are given in Table  B2. We incorporated Na, Si, Mg, K, Fe and their derivatives into VULCAN by adding 12 reactions (from kinetics.nist.gov) to the preexisting chemistry network based on C, H and O (Table C3). We initially added more reactions, but removed those that had a negligible impact on the resulting atmospheric speciation when omitted.

D. COLLISION CROSS SECTIONS
The CCSs are shown in Table D4, which were approximated by the circular area of radius equal to the atom or bond length. Furthermore, all bonds were approximated to be covalent.
E. ATMOSPHERIC P-T PROFILE SPIDER determines an atmospheric pressuretemperature profile through an analytical solution to the radiative transfer equations (see Appendix in Abe & Matsui (1985) and Sect. 3.7.2 in Andrews (2010)). The solution gives rise to the skin temperature equation (Eq. 7). Fig. E1 shows the volatile (V) and non-volatile (N3, N5) atmospheric pressure-temperature profiles that are used for FactSage and VULCAN calculations. An unphysical outcome of assuming only radiative equilibrium (no convection) is a temperature discontinuity between the base of the atmosphere and the surface of Note-The reactions given affect the speciation of Si, Mg, Fe, Na and K and Si between T = 2000 -873 K and P = 11.7 -10 −7 bar. Note-CS values are based on sizes of atomic, single, double, and triple bond data (Clementi et al. 1967;Pyykkö & Atsumi 2009a,b;Pyykkö et al. 2005).
the magma ocean, which is visually more evident for the non-volatile cases that have a small optical depth. Nevertheless, for all cases the surface temperature is 2000 K.