Dust Accretion onto Exoplanets

Accretion of interplanetary dust onto gas giant exoplanets is considered. Poynting-Robertson drag causes dust particles from distant reservoirs to slowly inspiral toward the star. Orbital simulations for the three-body system of the star, planet, and dust particle show that a significant fraction of the dust may accrete onto massive planets in close orbits. The deceleration of the supersonic dust in the planet's atmosphere is modeled, including ablation by thermal evaporation and sputtering. The fraction of the accreted dust mass deposited as gas-phase atoms is found to be large for close-in orbits and massive planets. If mass outflow and vertical mixing are sufficiently weak, the accreted dust produces a constant mixing ratio of atoms and remnant dust grains below the stopping layer. When vertical mixing is included along with settling, the solutions interpolate between the mixing ratio due to the meteoric source above the homopause, and that of the well-mixed deeper atmosphere below the homopause. The line opacity from atoms and continuum opacity from remnant dust may be observable in transmission spectra for sufficiently large dust accretion rates, a grain size distribution tilted toward the blowout size, and sufficiently weak vertical mixing. If mixing is strong, the meteoric source may still act to augment heavy elements mixed up from the deep atmosphere as well as provide nucleation sites for the formation of larger particles. The possible role of the Lorentz drag force in limiting the flow speeds and mixing coefficient for pressures $\la 1\, \rm mbar$ is discussed.


INTRODUCTION
In our solar system, dust is constantly produced by both asteroidal and cometary sources (see Moro-Martin 2013 for a review). Under the action of Poynting-Robertson (PR) drag, this dust will slowly spiral inward toward the Sun. When the dust approaches a planet, there are several possible outcomes: a physical collision with the planet; ejection from the solar system; or the dust will pass by the planet and again inspiral until it reaches another planet or eventually the Sun. Dust which approaches close enough to the Sun may sublimate completely or be ejected due to enhanced radiation pressure after partial sublimation (the β-meteoroids).
Interplanetary dust may be directly observed as zodiacal light through the scattering of sunlight, and also through thermal emission (Krüger & Grün 2014). The inner part of the zodiacal light is called the F-corona (Szalay et al. 2020), and it dominates the reflected light down to a few solar radii, inside of which the dust is destroyed (Mann et al. 2004;Howard et al. 2019). Upon entering the Earth's atmosphere, the stopping of the fast-moving dust particles is observable through radar reflections off their ionized trails, or as visible meteoroids (Whipple 1950;Opik 1976;Plane 2012). The gas-phase atoms released by ablation give rise to an enhancement of heavy atoms, which is utilized by astronomers to make laser guide-star observations with adaptive optics systems (Olivier & Max 1994). Some gas-phase atoms condense back into small dust grains called meteoric smoke particles (Rosinski & Snow 1961;Hunten et al. 1980). They may act as nucleation sites around which an ice mantle is accreted, generating significant scattering opacity. These high altitude noctilucent clouds are coincident with the dust-stopping layer, and their composition and influx rate have been constrained by extinction observations (Hervig et al. 2017).
This paper presents a study of the accretion of dust and its possible observable consequences for exoplanets orbiting very near their star. Infrared excess due to dust thermal emission has been observed from other arXiv:2206.09093v1 [astro-ph.EP] 18 Jun 2022 stars, with the far, mid, and near infrared arising from cold, warm, and hot dust at different distances from the star (Trilling et al. 2008;Absil et al. 2013;Eiroa et al. 2013;Ertel et al. 2014Ertel et al. , 2020Absil et al. 2021). It is the hot dust component that would directly interact with close-in planets. In our solar system, the inner solar system sources of dust may be dominated by Jupiterfamily comets with a smaller component from the asteroid belt (Nesvorný et al. 2010(Nesvorný et al. , 2011. The dust mass input rate (Ṁ d,iss ) to the inner solar system, which is required to explain the zodiacal light and meteoroids impacting the Earth's atmosphere, is roughly constrained to beṀ d,iss ∼ 10 7 g s −1 (Nesvorný et al. 2011).
There are two reasons why higher dust accretion rates may be possible for extrasolar planets as compared to the Earth. First, in order for a star to have a detectable excess from hot dust requires vastly more dust than found in our own solar system (the zodi; Roberge et al. 2012). The near-infrared interferometric detection of a H-or K-band excess requires > ∼ 3600 zodi's (Absil et al. 2013). Detection is more favorable for warm dust at ∼ au separations, e.g., LBTI can detect warm dust down to tens of zodi for some systems (Ertel et al. 2020). While ∼ 15 − 20% of stars have detectable hot dust excess (Absil et al. 2013;Ertel et al. 2014;Absil et al. 2021), undetected systems with far largerṀ d could be more common and perhaps ubiquitous.
Second, Bonsor et al. (2018) have shown that massive, close-in planets are surprisingly effective at accreting dust as compared to smaller and more distant planets. The Earth is estimated to accrete only a tiny fraction of the dust passing through its orbital radius (Ṁ Earth ∼ 10 2 − 10 3 g s −1 , Nesvorný et al. 2011); hence, some exoplanets may experience much larger dust accretion rates much larger than Earth's, with correspondingly higher gas-phase atomic and dust densities in the atmosphere.
The goal of this paper is to estimate, in the context of simplified models, if these densities can be large enough to affect the atmospheric abundances, or to give rise to significantly larger transit radii. It will be shown that this may be accomplished for a range of ratesṀ d = 10 7 − 10 9 g s −1 for systems with favorable parameters.
One motivation for this work came from the puzzling transmission spectra of extrasolar planets (see the detailed review in Madhusudhan 2019). For the case of a clear atmosphere, atomic lines such as the Na D doublet may be expected to dominate the opacity over a wide band in the optical, matching onto Rayleigh scattering by H 2 toward shorter wavelengths. While such a feature has been observed (Nikolov et al. 2018), it is much more common that only a narrow range of wavelengths will show the expected Na D feature, and outside the center of the line the transit radius varies slowly with wavelength. At longer wavelengths, absorption features from molecules such as H 2 O and CH 4 are expected in a clear atmosphere, but again these features are often muted. An interpretation of these observations is that the continuum absorption corresponds to aerosols made of condensates or due to photochemical reactions high in the atmosphere. The estimated pressure range for these high altitude clouds must be µbar − mbar, well above the ∼ 0.1 bar level where the optical continuum is absorbed. Lofting large particles to such high altitudes requires strong vertical mixing (e.g. Powell et al. 2019).
In addition, species like Mg, Si, and Fe, which readily condense to form clouds near the P = 1 bar level of the atmosphere (e.g. Burrows & Sharp 1999), are observed with transit radii well above that of the optical continuum (e.g. Madhusudhan 2019). Since the abundances of these elements should have been greatly reduced above the cloud deck (Visscher et al. 2010), the origin of these elements in atomic form high in the atmosphere is difficult to explain. While large amounts of vertical mixing are required to mix heavy atoms and dust grains to high in the atmosphere, the meteoric source discussed in this paper would provide an alternative explanation that does not require strong mixing.
Earlier work by Lavvas & Koskinen (2017) has already raised the possibility that meteoroid ablation may seed the upper atmosphere with heavy elements, which may then form the aerosols and give rise to the observed continuum absorption in the transit spectrum. However, that study was mainly concerned with the formation of photochemical hazes through mixing of heavy molecular species upward from the deep atmosphere to the ∼ 1 µbar level, where they would undergo photolysis. The authors concluded that while meteoroid accretion would deposit heavy elements in the upper atmosphere, the rate and hence the importance of the meteoric source was uncertain. This was a motivation for the present work in which a plausible range for the dust accretion rate is discussed, even if the value for individual systems is still uncertain. The more detailed calculations in Lavvas & Koskinen (2017) provide context for the present study in that once a dust accretion rate is assumed, the coagulation of the small particles and effect on the transit spectrum may be estimated.
A possible source of heavy elements that may be accreted into the planet's atmosphere, aside from the case of interplanetary dust studied here, is gas loss from hypothetical satellites orbiting the planet (Johnson & Huggins 2006;Cassidy et al. 2009;Oza et al. 2019), similar to the Io plasma torus around Jupiter.
A somewhat analogous situation to the exoplanets occurs for white dwarfs (Gänsicke et al. 2006;Farihi et al. 2009Farihi et al. , 2010. As settling times on white dwarfs can be a small fraction of their age, one might expect pristine photospheres showing little evidence of heavy elements. However, evidence of heavy elements are indeed found in their spectra for a significant fraction of white dwarfs, implying recent accretion. The observed abundance in the atmosphere is then a competition between accretion and settling (e.g. Bauer & Bildsten 2019), similar to the case studied here for accretion of dust onto planets. The existence of accreted heavy elements is perhaps more surprising in the white dwarf context as the post-mainsequence evolution of the star would effectively destroy their inner solar systems, and the orbits of any remaining bodies would expand significantly as the star loses mass (e.g. Mustill & Villaver 2012;Sun et al. 2018).
The outline of the paper is as follows. In Section 2, the dust orbital dynamics calculation of Bonsor et al. (2018) is revisited and their central result affirmed. Section 3 discusses the entry of the dust into the planet's atmosphere, the deceleration to subsonic speeds, and subsequent ablation as a function of altitude. Section 4 shows that the mixing ratios below the source, if gravitational settling is dominant, may approach solar abundance for a dust accretion rate comparable toṀ d,iss . Section 5 contains a discussion of the expected amount of mixing by fluid circulation in the upper atmosphere. The dust source layer is matched onto the rest of the atmosphere in a more complete manner in Section 6, including the effects of turbulent mixing, gravitational settling, and upward escape of hydrogen. Section 7 uses the profiles of mixing ratio to compute the dependence of the transit radius on vertical mixing coefficient and dust accretion rate. The size distribution and number density of the remnant dust is discussed in Section 8. Discussion, comparison to previous work, and estimates for Lorentz drag are given in Section 9. Section 10 contains a summary of the calculations and main results. The Appendix contains a calculation of the accretion fraction onto the planet using theÖpik approximation, and comparison to the Rebound simulations.
Calculations throughout this paper will be presented for three different planet sizes, a small planet with M p = 0.1 M J and R p = 0.7 R J , a medium-size planet with M p = M J and R p = R J , and a large planet with with M p = 2.0 M J and R p = 2.0 R J . A Sun-like star has parameters M s = M , luminosity L s = L , radius R s = R and dust sublimation radius R sub = 4 R .

DUST ORBITAL DYNAMICS
Solutions of the three-body problem, including the orbital motion of the star, planet, and dust particle are considered in this section. The goal is to calculate the fraction of the dust that is accreted onto the planet, sublimated near the star, or ejected from the system. For the dust that hits the planet, the impact speed and angle of incidence will be computed.
Two solution methods are used: numerical integration of the equations of motion (Section 2.1), which makes few approximations but is time-consuming; and theÖpik approximation (Appendix A), which greatly speeds up the calculation at the expense of approximations for the orbital dynamics.
The collisional dust source and any distant planets that scatter the dust are assumed to be distant so that PR drag has nearly circularized the dust orbit by the time it reaches the inner planet of interest. This planet is assumed to be on a circular orbit around the star. The dust feels the gravitational forces from both the star and planet, as well as radiation pressure and PR drag forces from the star. There are several phases for the orbital dynamics, as illustrated in Figure 1.
First, as PR drag removes energy and angular momentum from the orbit, the dust will slowly spiral inward on nearly circular orbits. Second, when it gets close enough to the planet, and for assumed small initial eccentricity and inclination, it will enter into a mean-motion resonance with the planet, temporarily halting the inspiral. Continued PR drag causes the eccentricity to rise to the point that the dust orbit begins to cross the planet's orbit. Due to the resonance, the planet is initially far away from the dust near its pericenter and close encounters are avoided. However, as PR drag pushes the orbit further into the resonance, the libration amplitude of the resonant angle increases and eventually the orbit leaves resonance, ending the second phase.
In the third phase, the dust orbits are no longer in resonance, and begin with semi-major axis outside that of the planet, but pericenter inside the planet's orbit, allowing close encounters between the two. There are three fates for the dust particle in this phase: physical collision with the planet, orbits that intersect the dust sublimation zone of the star, or ejection of the dust to infinity. For planets that are sufficiently distant from the star, dust particles that avoid both ejection and collision with the planet may have time to circularize before they encounter the dust sublimation zone. For close orbits, the eccentricity may be raised to sufficiently large values to enter the sublimation zone directly.  Direct numerical integration of the equations of motion is performed with the Rebound code (version 2.20.5, Rein & Liu 2012). The star, of mass M s and dust sublimation radius R sub , and planet, of mass M p and radius R p , are on a circular orbit of radius a p . The 15th order Gauss-Radau integrator with adaptive stepsize (IAS15) is used. The radiation pressure and PR drag forces due to the stellar flux are (Burns et al. 1979) and are included as additional, velocity dependent forces. The dust position is x and the positions of the star and planet are x s and x p , respectively. The stardust separation vector is then x ds = x − x s = r ds n ds , and the relative velocity is is the ratio of radiation to gravitational force on the dust, L s is the stellar luminosity, ρ d is the dust particle density, s is the dust radius, and m d = 4πρ d s 3 /3 is the dust mass. The quantity Q pr is a dimensionless coefficient depending on the dust size, composition and the stellar spectrum (Burns et al. 1979). For large dust, which obeys the geometrical optics limit, Q pr 1 and β 2.9 × 10 −3 (2 g cm −3 /ρ d )(s/100µm) −1 (L s /L )(M s /M ) −1 . The radiation pressure term effectively decreases the mass of the star to M s (1 − β) so that the dust would orbit more slowly than the planet for a circular orbit at the same separation. The inspiral timescale, for a circular orbit of radius a, is much longer than the orbital timescale. The simulations are run until a physical collision or ejection occurs, or until a maximum run time of 90 t pr (a p ), which is sufficient to ensure that most dust particles have been accreted or ejected. A by-product of using fixed integration times is that some runs in which particles have diffused to a a p orbits will not have had sufficient time to finish. However, these runs are small in number and do not affect the outcome probabilities significantly.
Rebound's direct collision scheme is used to detect collisions of the dust with either planet or the sublimation zone around the star. The relative velocity and angle of the velocity with respect to the radial direction are tabulated at the time when the dust particle hits the surface.
Calculations were carried out for a Sun-like star, each of the three planet sizes, and a range of semi-major axis a p = (4 − 20) × R . Results are presented for a dust orbital inclination I = 10 • relative to the planet's orbit. The initial dust semi-major axis is a = 3a p , initial eccentricity is zero, and for each orbital separation, 1000 runs were carried out distributed uniformly over initial orbital phase.
The dust does not evolve during the simulation in that the value of β is fixed and the dust is assumed to be instantly destroyed if it enters the sublimation radius R sub around the star. For a Sun-like star this corresponds to a radius where the dust equilibrium temperature would be higher than 2000 K. Depending on the sublimation timescale as compared to the orbital period, dust mass loss may lead to an increase in the importance of the radiation pressure force. In the case without a planet, this has been shown to lead to ejection from the system when β > ∼ 1 (the β-metereoroids; e.g. Kobayashi et al. 2009). However, when a planet is included collisions with the planet will compete with ejection. This effect is not included in the present study.
Selected results are presented in Figures 1 -6. An example trajectory for β = 0.01 is shown in Figure 1, and can be interpreted in terms of the analytic calculations in Weidenschilling & Jackson (1993). The top panel shows the semi-major axis versus eccentricity for the orbit of the dust around the star. Initially the orbit has small eccentricity and is decaying due to PR drag. It then enters a 2:1 mean motion resonance with the planet (Sicardy et al. 1993;Weidenschilling & Jackson 1993;Quillen 2006), after which the semi-major axis becomes constant at a = 2 2/3 (1−β) 1/3 a p and the eccentricity rises. The trajectory crosses the dashed line for a(1 − e) = a p and the dust and planet orbits now cross. The eccentricity rises to a maximum e 0.5 (Weidenschilling & Jackson 1993), and then decreases again as oscillations in a increase. When the resonance is broken, subsequent encounters with the planet's Hill sphere give rise to large changes in a and e (cometary diffusion, see e.g. Yabushita 1980, and the Appendix). This partic-ular trajectory ends with the dust particle hitting the planet.
The middle panel of Figure 1 shows the resonant angle ϕ = λ − λ − versus time, where λ and are the mean longitude and longitude of pericenter for the dust particle, and λ is the mean longitude of the planet. Initially this angle is circulating, but begins to librate in resonance at t 600 yr. Thereafter, the libration amplitude slowly increases over ∼ 15 t pr (a p ). Near t 1800 yr, the orbit falls out of resonance and ϕ begins to circulate. At this point close encounters may occur, and a short time later the dust particle hits the planet.
The lower panel of Figure 1 shows the Jacobi constant as calculated in the inertial frame. In the absence of PR drag, T is a constant of the motion. The value of T as the resonance is exited is crucial for determining the fate of the dust and the speed at which it will hit the planet. Here r dp is the dust-planet separation, v p = (GM s /a p ) 1/2 is the planet's orbital speed and n p = (GM s /a 3 p ) 1/2 is the planet's mean motion. The figure shows that T initially decreases due to PR drag during the inspiral phase. The planet enters the resonance at T 3.1 and continues to decrease, albeit more slowly, and exits the resonance with T 2.7. Subsequent decay of T on a timescale ∼ t pr (a p ) can then occur due to PR drag outside of resonance, but it does not change due to encounters inside the Hill sphere. Figure 2 shows the fraction for each of the four possible outcomes as a function of a p for the medium-size planet. These fractions will be called f pl , f st and f ej to hit the planet, hit the dust sublimation zone, or to be ejected, respectively. Each line is labeled with the dust fate as well as the value of β.
Immediately apparent is that physical collisions with the planet (blue lines) are a significant fraction of the outcomes for massive planets near the star. This important result was first shown by Bonsor et al. (2018). It motivates the following sections of this work on the effect of dust accretion on the planetary upper atmosphere.
For dust sublimation radius R sub = 4 R , the a p = 4.2 R runs have f st 100%. At larger a p , f st is still significant as the random walk in a and e allows some particles to lower their pericenter while holding T fixed. Interestingly, f st is sensitive to PR drag since the β = 0.01 run has much smaller f st than the β = 0.1 and 0.5 runs.
The ejection fraction f ej is the opposite, with ejections more common for smaller β. The fraction hitting the planet is larger for smaller β dust due to the longer inspiral timescale. The ejected fraction grows with a p and Outcome probability versus planet semi-major axis ap for a Sun-like star and medium-size planet. The three physical outcomes are colliding with the planet (pl, blue lines), dust sublimation near the star (st, red lines), and ejection from the system (ej, cyan lines). Since the simulations are run for a fixed time, 90tpr(ap), a small number had not yet achieved one of the three outcomes and are incomplete (inc, magenta lines). Each line is labeled by β = 0.01, 0.1, and 0.5. An initial inclination I = 10 • has been used for all runs.
is expected to dominate at sufficiently large a p (see the Appendix). Figure 3 shows the results for the small planet with β = 0.1, 0.5. The main difference from Figure 2 is now the complete absence of ejected particles. Impacts on the dust sublimation zone dominate, and the fraction hitting the planet has values of f pl = 10 − 20%, similar to the medium-planet case for the same β. Additional runs for an Earth-mass planet and β = 0.01 and 0.1 find maximum f pl = 1−10% and a similar outward decrease.
The opposite happens for the large planet case shown in Figure 4. Now ejections dominate and the sublimation zone only dominates at very small a p . The fraction hitting the planet is larger than the medium-size case for the same β.
While not shown here, runs were also performed for inclination I = 0 • and they showed significant differences. In particular, hitting the planet was even more dominant, with f pl > ∼ 0.9 for the medium-size planet over nearly the whole semi-major axis range, and f st much smaller. A few runs were then done at I = 1 • − 2 • and it was found that the f pl rapidly decreased to near the I = 10 • value over that small range. This may be due to a critical inclination ∼ R p /a p ∼ 1 • for which z < ∼ R p and collisions are enhanced. Hence, unless the dust accretion disk is very thin, the more conservative I = 10 • results may be more appropriate, but in any event low inclination dust sources favor accretion onto the planet. The value I = 10 • is comparable to the Hill inclination r Hill /a p (M p /3M s ) 1/3 ∼ 5 • for which z < ∼ r Hill . For larger inclinations it is expected that the collision rate decreases. It is also expected that if the planet were allowed finite eccentricity that the results would be similar for e p < ∼ r Hill /a p (M p /3M s ) 1/3 ∼ 0.07. The impact speed at which the dust enters the planet's atmosphere is closely related to the Jacobi constant. Let u = v−v p be the relative velocity of the dust and planet. Using the rotating-frame expression for the Jacobi constant evaluated at the position of the planet, the impact speed can be written as where u esc = 2GM p /R p is the escape speed from the surface of the planet and the second term u ∞ v p √ 3 − 2β − T ∞ is the relative velocity at infinity as the particle enters the planet's Hill radius. For closein orbits the two terms in Equation 5 are comparable   speed varies significantly between the small, medium, and large planets. In this context, a Keplerian orbit around the star may be used to approximate 1/2 cos(I)(6) to compare to Figure 1. Figure 5 shows a histogram of u ∞ /v p for four runs with Jupiter-like planet, a Sun-like star, and (a p /R , β) = (10, 0.01), (20, 0.01), (10, 0.1), (20, 0.1). The normalizations of each curve are different due to different f pl . Each distribution is peaked near u ∞ /v p = 0.4−0.6, with a tail extending to lower and higher values. For small β this corresponds to T ∞ = 2.65 − 2.85 at the time of impact. As in Figure 1, these values are much smaller than when the particle enters resonance, and are determined by the decrease in T during resonance as well as the subsequent decrease after the resonance is broken. Figure 6 shows the cosine of the angle between the inward radial direction and the dust velocity when it impacts the planet, for the same parameters as used in Figure 5. Aside from the difference in normalization due to different f pl in each run, each distribution is broad but dominated by vertical entry.

STOPPING AND ABLATION
The goal of this section is to study the stopping and ablation of the dust accreted by the planet. The distribution of ablated mass with altitude forms the meteoric source for gas-phase atoms discussed in Sections 4 -7. The size distribution of the partially ablated remnant dust particles will be used in Section 8 to compute dust abundance profiles and continuum opacity.

background
As dust grains enter the atmosphere at highly supersonic speeds, they are decelerated by collisions with gas particles. The collisions occur in the kinetic regime as the mean free paths of gas particles are much larger than the dust grain radius in the upper atmosphere. For the size ranges of interest here, s ∼ 0.1 − 100 µm, two main mechanisms contribute to dust ablation, sputtering, and thermal evaporation.
In thermal evaporation, frictional heating raises the dust temperature, leading to melting, vaporization, and mass loss from the dust particle. This is most important for larger dust particles that take longer to decelerate.
Sputtering refers to the process where an incident nucleus penetrates the surface layers of the solid grain, inducing a cascade of collisions among a large number of atoms in the solid, some producing a recoil sufficiently near the surface such that a target atom exits the solid (Johnson 1990). Sputtering occurs for all grain sizes, but is most important for small grains. For both processes, the ablated gas-phase atoms emerge at high speeds, and will be collisionally ionized. However, since the subsequent settling times of the ions (see Section 4 below) are longer than ionization and recombination timescales, they will soon come to the ionization state for that layer in the atmosphere.

radiative heating
While the dust particle is orbiting the star, the stellar radiation field sets the dust temperature. However, as it approaches the planet and is accreted, it also encounters the planet's radiation field, which has comparable flux to that of the star when near the planet. The equilibrium temperature of the planet's optically thick lower atmosphere is which can approach the dust grain sublimation temperature for close-in orbits. Here T s is the stellar effective temperature. As grains enter the atmosphere, they may be heated by both the stellar and planetary radiation fields, as well as frictional heating with the atmosphere.
After stopping, dust grain temperatures in the upper atmosphere are set by a balance of radiative heating and cooling, and are decoupled from the gas temperature. This may allow grains to survive even when the gas temperature is well above the melting temperature of the dust, as for H II regions in the interstellar medium. For gas temperature T ∼ T eq , number density n = P/k B T and thermal velocity v th = (8k B T /πµm u ) 1/2 , the heating rate due to collisions with the gas is nv th πs 2 2k B (T −T d ) (Draine 2011). The largegrain radiative heating rate due to a blackbody radiation field at temperature T eq is 4πs 2 σT 4 eq . Collisions with gas atoms are then only important (Chiang & Goldreich 1997;Draine 2011) at pressure levels deeper than well below the stopping layer. At P < P d,th , the dust temperature is set by the stellar and planetary radiation, while below this layer the dust temperature is regulated to be near the gas temperature.
To compute the radiative heating of the grain, the stellar radiation field is taken to be a blackbody with tem-perature T s and solid angle on the sky Ω = π(R s /a p ) 2 . The optically thick lower atmosphere of the planet is assumed to produce radiation with temperature T atm T eq and solid angle Ω = 2π. The gas in the upper atmosphere is optically thin to the bulk of the radiation field. For dust temperature T d , the net radiative heating rate is thenĖ where the step function Θ(day) is 1 on the dayside and 0 on the nightside. In the numerical calculations, the stellar irradiation will be included for simplicity, Θ(day) = 1, as if all the dust was entering on the dayside of the planet. These changes in numerical factors will affect the critical dust sizes for which emissivity is inefficient and also for which the temperature may be raised above the melting temperature by frictional heating. The dust grain absorption efficiency for a blackbody radiation field of temperature T is taken to be ε(s, T ) = min 1, (s/0.1 µm)(T/10 3 K) 2 , where the numerical coefficient in the second term can vary by 30% depending on composition (Draine 2011). This formula interpolates between the case of large dust for which ε = 1, and small dust with inefficient emission and absorption with ε ∝ s. The critical dust size below which the efficiency is less than unity is s cr (T ) = 0.1 µm (T/10 3 K) −2 . Large grains for which s > s cr (T eq ) will have dayside temperature T d = (3/2) 1/4 T eq 1.1 T eq independent of grain size. Small grains with s < s cr (T s ) can be superheated with dayside temperature T d 1.6 T eq (a p /10R s ) 1/6 , again independent of grain size 1 . Medium-size grains have a temperature dependent on size, interpolating between these two regimes.

stopping calculations
Numerical solutions are computed here using the equations from Plane (2004), which are used to model stopping in the Earth's atmosphere, as well as alterations described below. While this model only includes thermal evaporation, sputtering will be included here as well (see also Vondrak et al. 2008 for a more comprehensive treatment of both processes). The equations for the altitude z, vertical velocity u z and horizontal velocity u 1 Assuming that they are not so small that heating by collisions with gas particles is dominant over radiative heating and cooling.
The stopping layer is idealized as being geometrically thin with constant gravity g, as well as isothermal with gas temperature T = 3000 K and mean molecular weight µ = 1.23, appropriate for a mixture of atomic H and He at a pressure level near P = 1 µbar (Huang et al. 2017).
Here ρ(z) = ρ 0 exp(−z/H) is the atmospheric density, ρ 0 is the density at the 1 µbar level and H = k B T /µm p g is the scale height, ρ d = 2 g cm −3 is the assumed dust grain density, u = u 2 x + u 2 z and Γ = 1.2 is the drag coefficient. Below the melting temperature, or above the melting temperature but for net negative heating, the equations for the dust temperature T d and radius s are while above the melting temperature and for net positive heating they become Here Λ = 0.5 is the heat transfer coefficient, C d = 10 7 erg g −1 K −1 is the specific heat, σ is the Stefan-Boltzmann constant, L = 3 × 10 10 erg g −1 K −1 is the latent heat, and the melting temperature is taken to be T melt = 1800 K. The choices for the coefficients are discussed in Plane (2004) and Vondrak et al. (2008). We include a prescription for the emissivity of small dust, which is more important for the high temperatures near the star than for the Earth.
In the frame of the dust grain, atmospheric atoms of mass m approach the grain with energies E = 50 eV(m/m u )(u/100 km s −1 ) 2 , where m u is the atomic mass unit. The sputtering rate is proportional to the number of impacts with energy E above the threshold E thr required for sputtering, which is typically much larger than the lattice energy. While H atoms are the most abundant atmospheric particles, their energies E are significantly lower, and He and O impactors will also be included as they are more robustly above the sputtering threshold energies.
The low-energy approximation for sputtering yields from Draine & Salpeter (1979) is used (see their Section 4.b), in which the dust mass-loss rate is computed aṡ and the sputtering yield has the form Y Here A is a normalization constant and the sputtering threshold energy is E thr,i . This form for the mass-loss rate has been used to derive the change in radius in Equations 14 and 16. A MgSiO 3 composition is assumed for the dust, for which the average ejected atomic mass from the target ism t = 20 m u . The sum is over the number density n i and sputtering yield Y i of each of the three impactors, with the formulas given in Equations 31-33 of Draine & Salpeter (1979), and shown in their Figure 4. The low energy formula increases as E 2 well above threshold, and is used until the broad plateau seen in their Figure 4 is reached. A maximum yield Y i,max is then enforced. For the number densities n i of the impactors, solar abundance ratios are used for H, He, and O 2 . The threshold speeds for Y i to be nonzero are 76, 33, and 17 km s −1 for H, He and O, respectively, and they approach the maximum yield at significantly higher speeds. Hence, H can only sputter at fairly high incident speeds, while O's threshold speed is below the escape speed for the planets considered here. While O has a low threshold speed, it's low abundance (here taken to be solar) gives a maximum yield at high energy of (n O /n)Y O < ∼ 5×10 −4 . The high energy yield of H and He are much higher, at (n H /n)Y H ∼ (n He /n)Y He ∼ 3×10 −3 .
Sputtering typically leads to a tail of mass deposition higher in the atmosphere than that for thermal evaporation, since sputtering shuts off after deceleration. For ablation solely by sputtering, and assuming the yield can be expressed in the power-law form i n i Y i ≡ nY 0 (u/u 0 ) α , then for vertical incidence the ds/dt and du z /dt equations can be combined to yield a solution for the final grain mass m fin = m init exp(−γ), where the exponent is Hence, for α = 4 above threshold and below the maximum yield, these fiducial parameters imply that each dust grain will undergo mass loss by the amount of m fin /m init exp(−0.2) = 0.8. This base level of mass loss is roughly what is found for dust which does not undergo thermal evaporation. The exact values can vary depending on if and by how much each impactor is above threshold. For instance, if all species are at such high speeds than the yields are constant (α = 0), a power-law solution m fin /m init (u fin /u init ) γ results, where now γ in Equation 18 is redefined without the 1/α factor and Y 0 may be a factor of 2 higher if the contribution from H is included.
Ignoring ablation, deceleration requires the dust particle of mass m d = 4πρ d s 3 /3 to collide with its own mass in atmospheric particles, πs 2 Hρ m d . This allows an estimate of the pressure level for stopping to be ρ d gs 0.2 µbar ρ d /2 g cm −3 g/10 3 cm s −2 (s/1 µm). This estimate is applicable for smaller dust that does not undergo thermal ablation. For larger dust, a better approximation is that most of the mass is deposited where the particle reaches the melting temperature. Equating the frictional heating per unit area, Λρu 3 /2, to the thermal emission per area, 4σT 4 d , gives an estimate for the stopping pressure to be in agreement with the numerical results. Note that this estimate is independent of the dust size. The stopping pressure depends strongly on both the melting temperature and the entry speed. The value T melt = 1800 K is adopted here, motivated by the more detailed study in (Vondrak et al. 2008). Models in which elements can have different melting temperatures have also been used to better model the deposition of Na and K for the Earth's atmosphere (Vondrak et al. 2008). With T melt = 1200 K for Na, they will be deposited at a factor of 10 lower pressure. Similarly the higher melting temperatures for Ca and Al minerals would lead to deposition at slightly larger pressures.
Such sublimation could also occur before the dust accretes onto the planet for orbits close enough to the star. The case of Na and K are particularly interesting. It is expected that planets interior to the sublimation zone for Na and K would not accrete these elements as they would be lost from the dust particles before they are accreted onto the planet. If dust accretion can affect the transmission spectra of exoplanets, it may be expected to see a difference for planets interior to the sublimation zone as compared to exterior planets. Figure 7 shows trajectories for three different size dust particles. The parameters used are a medium-size planet, a Sun-like star, an orbital radius of a p = 10 R , and an initial speed u = 90 km s −1 at vertical incidence. The equilibrium temperature is T eq = 1300 K, and the expected dust temperature set by radiation is T d,eq 1.5 1/4 T eq = 1440 K. The three dust sizes are s 3, 30, and 300 µm.

results
The s = 3 µm dust (solid line) rapidly adjusts to T d,eq and frictional heating leads to only a small temperature increase. It decelerates by pressure levels P ∼ 1 µbar and suffers only minor mass loss due to sputtering. The s = 30 µm dust (dotted line) must travel deeper into the atmosphere before it begins to decelerate. Frictional heating causes significant decrease in dust size over a fraction of a scale height. Rapid deceleration then leads to a decrease in frictional heating and temperature, and it cools back to T d,eq . The s = 300 µm dust (dashed line) also undergoes strong frictional heating, but does not decelerate in time to avoid catastrophic mass loss. The trajectory was stopped when the dust reached a size s = 0.01 µm.
While this example used a fixed entry speed, Equation 19 implies that a range of speeds will lead to a range of stopping altitudes and distribution of mass deposited with altitude. Since large dust undergoes more mass loss, grain size distributions with mass dominated by large grains may deposit more mass as gas-phase atoms, while small grain size distributions may retain more remnant grains that have undergone little ablation. Lastly, planets closer to the star will have T d,eq closer to T melt and the smaller required heating to melt implies mass deposition higher in the atmosphere as compared to planets further from the star.
Further progress requires a choice of grain size distribution that depends on the size distribution and number density at the source, PR drag and radiation pressure, collisional agglomeration, and erosion, encounters with planets, and the stellar wind environment near the star (Grun et al. 1985;Ishimoto 2000;Wyatt 2005;Thébault & Augereau 2007;Vondrak et al. 2008;Kobayashi et al. 2009;Nesvorný et al. 2010Nesvorný et al. , 2011Plane 2012;van Lieshout et al. 2014;Carrillo-Sánchez et al. 2015;Rieke et al. 2016). Two extreme limits are considered here. Grun et al. (1985), G85 from here on, gave a model for the solar system dust size distribution, including both PR drag and grain-grain collisions. Due to the relatively low grain density inside 1au, the collision rate is small for this model and it has both mass and area dominated by larger grains s ∼ 100 µm. We augment their result by imposing a cutoff at the blowout radius. By contrast, the calculations in van Lieshout et al. (2014) are for a much higher grain density in which erosive collisions lead to a distribution in which number, mass, and area are dominated by grains with β ∼ 0.5 near the blowout size.
Motivated by Figure 5, the distribution of u ∞ /v p is chosen to be a Gaussian with mean 0.5 and standard deviation 0.1, and the escape speed contribution is then included as u = u 2 ∞ + 2GM p /R p . The angle distribution is chosen to be P (µ)dµ = 2µdµ to represent the curves in Figure 6. Figures 8 and 9 show the results for dust ablation for the G85 size distribution. For comparison, the G85 distribution is shown in Figure 19 as the dashed line. Figure 8 shows the ablated mass per volume per time, ρ, as a function of pressure. Three different planet sizes are shown as well as four different orbital separations. For all planet masses, at orbital separation a p = 6 R the equilibrium temperature is high enough that only a small amount of frictional heating leads to thermal evaporation immediately near the top of the atmosphere, which is set to be at P = 1 nbar. For the other three orbital separations, the mass is distributed much lower in the atmosphere. All three planet masses show a peak higher up due to sputtering, while the M p = 1 and 10 M J runs show a secondary peak lower in the atmosphere due to thermal evaporation. The fraction of the incident dust mass that is ablated is larger for the high mass planet as compared to the two lower mass cases. Figure 9 shows the fraction ( abl ) of incident dust mass which is ablated as gas-phase atoms and which remains as remnant dust particles ( rem = 1 − abl ). Again, the large speeds at small a p imply all the dust is ablated. For larger separations, the escape speed has a significant effect, with the massive planet having a much larger amount of ablation. Beyond a p ∼ 10 R the escape speed contribution dominates over the orbital relative velocity. Figure 10 is the same as Figure 8, but for a deltafunction size distribution with all initial grains with β = 0.5, or s = 0.57 µm. The contribution from sputtering is still present in this case, but the deeper thermal ablation peaks are absent as the grains are able to decelerate before rising to the melting temperature. The fraction of mass which is ablated and remaining in the remnant grains is similar to that found for the G85 distribution in Figure 9.
The scenario studied here, of a distant dust source and nearly circular orbits, will produce small collision speeds. The collision speeds for highly eccentric orbits will tend to be larger. A dust particle with a a p and a(1 − e) a p will have relative speed u ∞ 3 1/2 v p due to radial motion near the escape speed. A dust particle with a a p but a(1 − e) a p will have u ∞ v p 2(1 − β) ± 1 for prograde and retrograde dust orbits. This formula gives 12 − 72 km s −1 for the Earth. Significantly higher dust speeds will result in more of the dust mass being ablated into gas-phase atoms than the cases considered here.

GAS-PHASE ATOMS
In this section the abundance profile of a trace species deposited into the atmosphere by dust accretion, and slowly drifting downward due to gravitational settling, will be considered. It will be shown that a constant mixing ratio results (see e.g., Hunten et al. 1980), and that the resulting abundances are comparable to solar composition for a dust accretion rate comparable toṀ d,iss .
The effects of vertical mixing, molecular diffusion, and upward motion of the background species due to atmospheric escape will also be included in Section 6.
Let the trace species have mass m, number density n(z), vertical speed u(z) and flux F (z) = n(z)u(z). Similarly, the background species, which we will consider to be H 2 , has mass m b , density n b (z) velocity u b (z) and flux F b (z) = n b (z)u b (z). The gas is assumed isothermal with temperature T . The collision frequency of the trace with the background will be written as ν = σv n b .
In the hard-sphere model, the transport collision rate coefficient between a trace particle and the background particles is (Chamberlain & Hunten 1987 .7 × 10 −9 cm 3 s −1 T 10 3 K 1/2 Figure 8. Deposition rate of ablated dust mass per volume as a function of pressure for the G85 size distribution. The black, red, and blue lines correspond to small, medium, and large planets. The solid lines that form spikes near P ∼ 1 nbar are for close-in planets at ap = 6 R . The dotted, short-dashed, and long-dashed lines are then for ap/R = 10, 14, and 18. The assumed accretion rate into the atmosphere is f plṀd = 10 8 g s −1 . Here m b is the mass of the background species and R coll,b and R coll are the hard-sphere radii for the background and trace. Atomic radii are only measured for a few species, as discussed in García Muñoz (2007) and Banks & Kockarts (1973). Following Banks & Kockarts (1973), from here on we will assume a radius R coll,t = R coll,b = 1.5 × 10 −8 cm for all species. A heavy atom accelerated downward by gravity g, and colliding with the lighter background atoms with timescale ν −1 will have a downward drift velocity of where the gravity of Jupiter is g J = 2700 cm s −2 . The settling time over a scale height is t settle = Hν g = 1.4 × 10 8 s P 1µbar Figure 9. Fraction of decelerated dust mass in gas-phase atoms and remnant dust particles. The G85 dust distribution has been used, with a cutoff at β = 1. Solid lines are for dust and dashed for gas. The black, red, and blue lines represent small, medium, and large planets.
This timescale is much longer than photoionization and recombination times, and the atoms quickly come into local ionization equilibrium 3 . The settling speed can be used to find the mixing ratio where F < 0 for downward motion. Since the background density has canceled out of the right-hand side, the result will vary slowly with altitude.
3 Here the atoms are treated as neutral. Including the upward charge-separation electric field will effectively decrease the downward gravity, and using ion-neutral cross sections will give larger collision rates than the neutral-neutral rates considered here. Hence, the results underestimate the densities of ionized species. Neutral atoms are considered for simplicity.
The particle flux below the source can be written HereṀ d is the accretion rate of the dust population as it moves toward the star under PR drag, f is the abundance by number of this element in the dust particles, andm 15 m u is the mean mass of the elements making up the dust particles.
The results are now assembled to compute the mixing ratio below the meteoric source using Equation 23. The result is n n b = 4.1 × 10 −3 cm −3 f f pl ablṀd 10 8 g s −1 Figure 10. Same as Figure 8 but for all dust with the same initial size, set by β = 0.5.
The relative number of different elements scales as f /m for m m b , and favor small planets due to the smaller gravity.
The mixing ratios below the meteoric source from Equation 25 are shown in Table 1 and Figure 11. If f pl ablṀd 10 8 − 10 9 g s −1 the abundances of C, O, Na, Mg, Al, Si, S, Ca, Cr, Fe, and Ni are comparable to their solar abundance. Given the previous successes in modeling transmission spectra for several lines using solar abundances (e.g. Huang et al. 2017), the meteoric dust model provides a natural explanation for these abundances. However, some elements are notably underabundant, the main example being He, which has tiny abundance in meteoritic dust.
Given the rough agreement between the abundances due to the meteoric source and solar abundances, and the difficulty in measuring absolute abundances for a number of elements, it may be difficult to tell the differ-ence between heavy elements mixed up from below and the effect of dust accretion.

MIXING IN THE UPPER ATMOSPHERE
A crucial consideration for high altitude absorbers is the size of the vertical mixing coefficient K (short for K zz ). Opposite to the case of heavy elements mixed up from below, where mixing must be strong to loft them to high altitudes, mixing must be sufficiently weak for the meteoric source to dominate and set a constant mixing ratio 4 . Turbulent mixing, as parameterized by K, arises naturally in a convective atmosphere, where K ∼ u eddy eddy is related to the size ( eddy ) and speed (u eddy ) of the energy-bearing eddies. However, for strongly irradiated extrasolar planets, the entire atmosphere at P < ∼ 10 2 − 10 3 bar may be stably stratified (e.g. Arras & Bildsten 2006) and mixing is instead due to winds driven by day-night temperature differences. An analytic theory for K due to winds in stably stratified atmospheres has been developed by Zhang & Showman (2018), and a detailed comparison to global circulation models including tracer particles was carried out by Komacek et al. (2019). The latter paper includes radiation transfer in the two-stream and twofrequency band (optical and infrared) approximation, with the upper boundary near the infrared photosphere at P ∼ 1 mbar and the lower boundary well below the optical photosphere at P ∼ 0.1 bar. A drag force is included near the lower boundary, to represent friction with the interior planet where winds are negligible. Rayleigh drag is also included with constant drag timescale throughout the simulation domain. Excluding Rayleigh drag, and for long chemical equilibration times, their numerical results show turbulent mixing as large as K ∼ 0.01 c s H b ∼ 10 10 cm 2 s −1 for the case of strongly irradiated atmospheres with high wind speeds, where c s ∼ 1 km s −1 is the sound speed. Since the wind speeds in this simulation increased upward, K showed a strong vertical increase. When Rayleigh drag is in-cluded, wind speeds can be strongly suppressed for large drag, thereby decreasing K.
A key question is how to extend these results for K to lower pressures in the P ∼ mbar − µbar range appropriate for transit observations.
Since vertical ballistic motion of a fluid element at speed ∼ c s can only result in 1 scale height of displacement, the deposition of the stellar optical continuum flux at the optical photosphere, and subsequent transport of heat up to the infrared photosphere by radiative diffusion and fluid motion, may not be able to coherently drive fluid motions in the P ∼ mbar − µbar region, well into the optically thin portion of the atmosphere (along a vertical path) 5 . Rather, the local imbalance of optically thin heating and cooling, which taps into only a small fraction of the total stellar radiation, may drive the motions in the upper atmosphere. Simulations includ- Note-Columns 2 and 3 are the solar photospheric and CI chondrite abundances taken from Lodders (2003). All elements with abundances 1% that of Si or larger are shown. The fourth column shows the mixing ratio below the meteoric dust layer from Equation 25 for f pl ablṀd = 10 8 g s −1 , Mp = MJ and T = 1000 K.
ing appropriate input physics (optically thin, non-LTE, magnetic field effects, etc.) are required to address this question. A second effect that becomes of primary importance in the upper atmosphere is the Lorentz drag. Effectively, due to the finite magnetic diffusivity, gas moving perpendicular to the magnetic field will experience a drag force attempting to damp out that motion 6 . This is analogous to the friction with the deep planetary interior. The weather-layer, with strong winds relative to the mean planetary rotation, should be confined between these upper and lower layers where gas is forced to corotate.
To briefly review the origin of this force, stellar heating-driven flows at speed u attempt to advect the magnetic field (B), creating currents (J). Ohmic diffu-6 For a Jupiter-size magnetic field, magnetic pressure dominates gas pressure for P < ∼ 1 µbar and the gas will be incapable of inducing large perturbations to the field in the high conductivity limit.
sion 7 attempts to dissipate these currents with a diffusion coefficient η. When advection balances diffusion, currents of size J ∼ (c/4π)(uB/η) give rise to a drag force per volume JB/c ∼ uB 2 /(4πη) and Lorentz drag timescale t drag = 4πρη/B 2 . When the drag time becomes smaller than the characteristic flow timescale R p /u ∼ (10 10 cm/1 km s −1 ) ∼ 1 day, the drag forces can greatly decrease the flow speeds (Zhang & Showman 2018) and hence the mixing coefficient K.
In particular the Lorentz drag timescale becomes short at lower pressures, due to the decrease in both density ρ and also diffusivity η. The magnetic diffusivity is given by η 230 cm 2 s −1 T/K x −1 e (Perna et al. 2010), where x e is the mixing ratio of electrons. While it is appropriate to determine x e using collisional ionization of Na and K deep in the atmosphere, above the P ∼ 1 mbar level photoionization dominates (Fortney et al. 2003;Lavvas et al. 2014;Huang et al. 2017;Lavvas & Koskinen 2017). If a solar abundance of Na were fully ionized, the resulting x e 4.3 × 10 −6 would give a diffusivity η = 2 × 10 9 cm 2 s −1 T/10 3 K and a drag timescale where Jupiter's equatorial magnetic field has been used. In a detailed calculation for HD 189733b, Lavvas & Koskinen (2017) find x e ∼ 10 −8 (10 −6 ) at P = 1 mbar (1 µbar), giving the range t drag = 10 4 (10 −1 ) s. Even the longer drag time t drag = 10 4 s was found to greatly decrease the flow speeds and K by Komacek et al. (2019).
To estimate the upper atmosphere flow speeds and mixing coefficient, we follow the one-zone model of Zhang & Showman (2017), now applied locally to the upper atmosphere instead of the lower atmosphere. Using dimensional analysis on the horizontal Euler equation gives where the terms represent fluid acceleration, the Coriolis force, the day-night pressure gradient force, and the drag force. Here we have assumed large day-night temperature difference and ignored the factors of order unity in the more detailed formula in Zhang & Showman (2017). For weak drag and Coriolis terms, the flow speed is u ∼ c s . However, when drag forces decrease the flow speed, the two terms on the left-hand side become negligible and a balance of the two terms on the right-hand side gives which is smaller than c s when t drag is smaller than the fiducial (drag-free) horizontal advection timescale t adv ≡ R p /c s ∼ 1 day. As a numerical example, for the range t drag = 10 4 − 10 −1 s, and advection time is t adv = 10 5 s, the flow speed would be reduced a factor of t drag /t adv = 10 −1 − 10 −6 over the range P = mbar − µbar. Clearly, in the upper atmosphere of planets receiving strong ionizing stellar flux and which creates high ionization levels, the Lorentz drag timescale may be orders of magnitude smaller than t adv , and significantly decrease the flow speeds. Komacek et al. (2019) find a day-night morphology for the flow when drag is strong, as opposed to an equatorial zonal wind when drag is weak.
The characteristic vertical speed w ∼ u(H b /R p ) gives rise to a vertical mixing coefficient (Zhang & Showman 2018;Komacek et al. 2019) where t chem is the timescale for source terms to significantly change the tracer density due to e.g. photoionization and recombination times for atoms, and growth times for dust. Hence, K ∝ u 2 (u) for short (long) t chem . For long t chem and weak drag, u c s gives K c s H b (H b /R p ). This is the expected value for the lower atmosphere when drag is weak. In the strong drag case, plugging in the flow speed from Equation 27 gives which is smaller by a factor of t drag /t adv than the weak drag estimate. If the horizontal flow timescale t adv (t drag /t adv ) −1 becomes longer than the chemical timescale, the expression becomes which is now smaller by both the factors t chem /t adv 1 and (t drag /t adv ) 2 1. For the fiducial range t drag /t adv = 10 −1 − 10 −6 over the range P = mbar − µbar, this would decrease K by 10 −1 − 10 −6 for the linear scaling and at least 10 −2 − 10 −12 for the quadratic scaling.
While a more detailed calculation is warranted, it's clear that Lorentz drag is an important effect over the pressure range P ∼ mbar − µbar and that most likely K is much smaller than the fiducial value c s H 2 b /R p ∼ 10 10 cm 2 s −1 applicable in the lower atmosphere. If K is indeed reduced by many orders of magnitude, then the homopause at D = K may occur near near the base of the upper atmosphere, where Lorentz drag first begins to suppress the day-night flows. In this case, since heavy elements may not be mixed up from below, the meteoric source may be required for heavy elements to be present in the upper atmosphere probed by optical and infrared transmission spectra.

DIFFUSION, MIXING, AND ESCAPE
In Section 4, only the effect of gravitational settling was included. Here vertical mixing, molecular diffusion, and upward drag due to atmospheric escape of the background species are included as well. To goal is to understand the conditions under which the meteoric dust abundance estimates from Section 4 are applicable.
For a heavy trace species, the abundance profiles can be affected by their much smaller scale height. Let a = m/m b be the mass ratio of the heavy trace particle to the background particle. The scale height of the background is H b = k B T /(m b g) and that of the trace is H = H b /a. A second important consideration is the size of the molecular diffusion coefficient D as compared to the vertical mixing coefficient K. Typically D will increase upward faster than K and above the homopause altitude z h , at which D(z h ) = K(z h ), the atmosphere will no longer be well mixed and diffusion effects will become important.
The vertical momentum equation of the trace is which can be rewritten as where the diffusion coefficient is D = k B T /(mν). For the case of zero fluxes, F b = F = 0, the solution is called diffusive equilibrium where the background density profile n b (z) = n b (z 0 ) exp(−(z − z 0 )/H b ). Species heavier than the background (a 1) will decrease upward much faster than the background and their mixing ratios will become small. Hence, in diffusive equilibrium, the optical depth of the trace will decrease upward strongly and the transit radius will be limited to be near the homopause.
The transition between the lower, well-mixed layer and the upper layer where diffusion is important can be evaluated by adding in the phenomenological mixing term where K is the vertical mixing coefficient. The term with K tries to keep the scale height of the trace to be the same as that of the background. Following Chamberlain & Hunten (1987), the effects of a meteoric source layer and drag forces from atmospheric escape can be included through solution of the two equations for n(z) and F (z) While the detailed profiles of dF/dz = fρ/m from Section 3 could be used as the source here, for simplicity, the atomic source has been given by a Gaussian profile with total flux F s > 0, centered at altitude z s , and with width σ s . The density of the background will be referenced to the homopause as n Assuming a constant K everywhere, the diffusion coefficient can be written D(z) = K exp(Z). The boundary conditions are the mixing ratio ξ ≡ n(z 0 )/n b (z 0 ) at a reference altitude z 0 deep in the well-mixed layer, and F (z 0 ) = F esc − F s , where F esc is only nonzero when blowoff is occurring for large F b , as discussed below. This choice for the flux gives F = F esc at the top of the atmosphere and the source gives a downward flux −F s . The left panel of Figure 12 shows examples of abundances profiles including a meteoric source 5 scale heights above the homopause, and for three different abundances deep in the atmosphere. The mixing ratio due to the meteoric source is close to the ξ = 10 −6 line and so it looks almost straight. However, when the deep and meteoric mixing ratios do not match, there is a region a few scale heights thick over which the ratio transitions. For this case of a meteoric source well above the homopause, the solutions including the source (solid lines) retain large mixing ratio while the diffusive equilibrium solutions (dashed lines) decay sharply above the homopause for this case of a = 10.
The right panel of Figure 12 uses the same source parameters and a deep mixing ratio ξ = 10 −5 . The different lines represent different escape flux F b of the background. For small F b , the solution agrees with the rightmost solid line in the left panel, including the diffusive equilibrium above the source, giving rise to rapid falloff. As F b increases, the mixing ratio below the source increases and the upward decrease above the source is slower. For the value F b = F b,cr ≡ Kn b (z h )(a − 1)/H b the abundance actually increases below the source layer and is constant above. Figure 13 shows the case of a dust source two scale heights below the homopause. For the case of meteoric source layer mixing ratio smaller or comparable to the deep mixing ratio the effect of the dust source is small. However, when the dust source mixing ratio is larger than that of the deep atmosphere the abundances can have a significant impact on the abundance.
These numerical results can be recovered with an analytic solution below the meteoric source (Chamberlain & Hunten 1987) where the effective particle mass is Z 0 = (z 0 − z h )/H b and the mixing ratio below the meteoric source is This solution interpolates between the value deep in the well-mixed layer, ξ, and the value between the homopause and the meteoric source, Ξ. Above the meteoric source, the solution again decays as the diffusive equilibrium form but now with the slower decay a eff instead of a. When a eff ≤ 1, the drag force from the upward-moving light background particles overcomes gravity. This occurs for particle masses below a critical mass or equivalently, above the critical flux Above this critical flux, the particles with mass m may be carried along with the outflow, and their upward flux F esc will be proportional to the background flux F b (Chamberlain & Hunten 1987). The critical mass loss rate of the background species to give a blowoff state for the trace iṡ The numerical estimate is appropriate for an atomic hydrogen background and 28 Si atoms and for a mediumsize planet. Due to the m 2 scaling, the required mass-loss rate to cause a blowoff state ranges from ∼ 10 11 g s −1 for He to ∼ 10 13 g s −1 for Fe. The latter mass-loss rate is large enough to evaporate the entire planet away in ∼ 6 Gyr. Large transit depths of heavy elements require their presence high in the atmo-sphere, and a blowoff state may be invoked to accomplish this, requiring high mass-loss rates. The alternative suggested here is that the meteoric source may supply heavy elements to the upper atmosphere without the need for such high mass-loss rates.

TRANSIT RADIUS FOR ATOMIC LINES
Atomic resonance lines interact with the gas much more strongly than the neighboring continuum, lead-ing to larger transit depth in the neighborhood of the lines. If the τ = 1 radius for the continuum is R p and that of the line is R p + z(λ), then this annulus of area 2πR p z(λ) gives a contribution to the transit depth of size ∆F λ /F λ 2R p z(λ)/R 2 s over that of the continuum. The goal of this section is to give examples of how z(λ) depends on dust accretion rate and vertical mixing coefficient K.
For definiteness, the Na D doublet will be used for the line. Different values of vertical mixing coefficient will be used to show the effect of the homopause. The source is taken to be a delta function at P = 10 −2 µbar. Deep in the well-mixed layer the abundance is assumed to be solar with mixing ratio ξ relative to the background species. The abundance below the meteoric source will depend on the ablation rate abl f plṀd .
Constant temperature is not a good approximation over the entire atmosphere, as it is expected to vary from T ∼ T eq ∼ 1000 K near the P = 1 mbar level to T ∼ 3000 − 10, 000 K above the P = 1 µbar level due to thermospheric heating. While the wings of the line may have their τ = 1 altitude where T ∼ T eq , wavelengths near line center may have τ = 1 where T T eq . To take these two cases into account, but still using the isothermal assumption for the atmosphere, two cases will be shown. In the first case, more appropriate for the wings, the temperature is chosen to be T = 1500 K atmosphere with H 2 as the background. In the second case, more appropriate for line center, a T = 5000 K atmosphere with H atoms is used as the background. In both cases, the P = 1 mbar level with be chosen as the z = 0 point, but keep in mind that the scale height is very different in the two cases.
For absorber number density that increases strongly downward, the optical depth can be evaluated as the effective path length 2πR p H b times the density at the altitude of closest approach, giving an optical depth Here the cross section is where f os is the oscillator strength, H(x, a v ) is the Voigt function, x = (ν − ν 0 )/∆, a v = Γ/(4π∆), and ∆ = ν 0 2k B T /mc 2 . The number density n(z) of the absorber is the product of the background density and the mixing ratio in Equation 37. Plugging into Equation 44, root solving can then be used to find the value of z(λ) where τ (z(λ), λ) = 1. Figure 14 shows z(λ) as a function of wavelength λ for a T = 1500 K atmosphere with H 2 background for a medium-size planet. Four different profiles are shown, labeled by dust accretion rate and vertical mixing coefficient. When the meteoric source is absent (dashed lines), the density decreases strongly above the homopause since the Na scale height is 11.5 times smaller than the H 2 scale height. Even at the center of the line, where the cross section is largest, z(λ) extends above the homopause by less than a scale height. Hence, in the absence of dust accretion, the transit radius is limited by the homopause. Only in the case that the homopause extends to very high altitudes, where τ < 1 even with constant mixing ratio ξ, is the effect of the meteoric source on the line profile ignorable.
By contrast, the meteoric source, here with f pl ablṀd = 10 8 g s −1 , allows the transit altitude to be set by the source, rather than the homopause, for both cases where the homopause is well below the source. The solid red and blue lines agree near the center of the line, and only begin to disagree below the respective homopauses where the mixing ratio becomes ξ rather than η. Figures 15 -17 show the transit altitude for the hotter atmosphere with T = 5000 K and H atoms as the background, over a large range of f pl ablṀd in each plot and three different values of the mixing coefficient K.
For the large mixing ratio K = 10 10 cm 2 s −1 used in Figure 17, the transit radius is nearly the same for all the accretion rates, since the homopause is near the meteoric source, and the mixing ratio is ξ over nearly the entire atmosphere. By contrast, for the smaller value K = 10 6 cm 2 s −1 used in Figure 15, there is a significant difference in the transit altitude even for very small accretion rates. At f pl ablṀd > ∼ 10 2 − 10 3 g s −1 , the transit altitude of the line center begins to extend above the homopause. By 10 5 g s −1 , well belowṀ d,iss , the line core is approaching the source altitude. For the 10 7 − 10 8 g s −1 cases, the line core extends above the source altitude, and the line wings are also significantly affected. The intermediate case of K = 10 8 cm 2 s −1 shown in Figure 16, still shows a difference on the line wings for 10 9 g s −1 , and the line core is near the source for f pl ablṀd > ∼ 10 7 g s −1 .
These results shown that it is possible for the meteoric source to affect the transmission spectrum of strong lines if the vertical mixing is sufficiently small. The importance will be limited if very strong mixing extends all the way up to the meteoric source region. source (0, 10 6 ) (0, 10 8 ) (10 8 , 10 6 ) (10 8 , 10 8 ) Figure 14. Transit radius versus wavelength for the Na doublet for an isothermal atmosphere with T = 1500 K and H2 background. The medium-mass planet is used. The z = 0 altitude is set at P = 1 mbar and the source is taken to be a delta function at P = 10 −2 µbar, shown as the black line labeled "source". The four lines represent two different vertical mixing coefficients and two dust accretion rates, and are labeled by (f pl ablṀd (g s −1 ), K(cm 2 s −1 )). The red and blue horizontal lines are the homopause altitudes for K = 10 6 and 10 8 cm 2 s −1 , respectively.

size distribution
Let the size distribution of the remnant dust be P(s), which is normalized to unity, dsP(s) = 1, and has mean massm d = dsP(s)(4π/3)ρ d s 3 . Below the meteoric source, there is a flux of dust particles per size interval ds which is denoted Φ(s). This flux can be written as where the negative sign denotes a downward flux. The size distribution for the remnant dust was computed as part of the stopping calculations in Section 3. Figure 18 uses the G85 size distribution and parameters for the small planet. Shown are the size and mass distributions for the incident dust before it hits the atmosphere (dashed line) and the remnant dust after it has stopped (solid lines), for three different orbital separations. Large millimeter-size dust is depleted by as much as an order of magnitude in mass; however, the smaller sizes are less affected. Ablation is larger for closer-in planets.
Results are not shown for the a p = 6 R case since all remnant dust was in the lowest size bin (s = 0.01 µm), implying that it would have been further ablated if the calculation had not been stopped. All runs at this orbital separation have such large entry speeds that little mass is left in remnant dust. Figure 19 shows results for the medium-size planet. Now there is significant ablation at the peak of the mass distribution at all orbital separations shown, due to the higher escape speed. Ablation is larger for closer planets. Also apparent is small remnant dust, below the blowout size at s = 0.3 µm, created by ablation of larger dust. The amount of small remnant dust is larger for closer orbits. Figure 20 shows results for the large-size planet. In addition to significant depletion at large masses, micronsize dust is also noticeably depleted. The abundance of remnant dust below the blowout size is now significantly larger, although still smaller than the peak of the initial distribution.
The above cases were for the G85 dust-size model. For the single-dust-size model with β = 0.5, or s = 0.58 µm, the dust size is small enough that thermal evaporation does not occur and there is only mass loss due to sputtering.  mixing is negligible, the dust density is determined solely by the vertical settling velocity u, as in the atomic case discussed in Section 4. The number density of dust per size interval ds is a function of pressure and size, and will be denoted N (P, s). It is normalized to the total number density n(P ) = dsN (P, s), and is given by The dust velocity is computed as u(P, s) = −u set (P, s) where the (positive) settling speed is computed as u set = max(u ep , u st ), where the larger of the Epstein and Stokes drift speeds is used. At high altitudes, the gas mean free path is larger than the dust radius s, and the Epstein settling velocity applies with value 8 where ρ = µm u P/(k B T ) is the mass density, µ = 2.3 is the mean molecular weight in the lower atmosphere discussed here. For a delta-function size distribution P(s ) = δ(s − s), the number density is Since the dust density is proportional to the gas density, the mixing ratio is fairly constant with altitude. This result can be re-expressed as the ratio of dust and gas densities This estimate shows that for dust accretion rates near the solar system value, and micronsize dust, the dust is settling sufficiently fast that its mass fraction is small compared to the solar metallicity value 0.01.
8 This formula assumes that the downward gravity force is balanced by the upward drag force. However, above the optical and infrared photospheres, where the radiation field is highly anisotropic, both the stellar and planetary radiation fields imply a radiation force on the dust grain. In the case of the upward planetary radiation field, the dust will feel an upward radiative acceleration βpg so that the effective gravity of the planet becomes (1 − βp)g, and the setting speed is reduced by a factor 1 − βp. This may increase the grain density for small grains with βp < ∼ 1. Further, grains with the appropriate size and composition to give βp > ∼ 1 would experience a net upward radiative levitation force. While these effects may be important, they are ignored here for simplicity.
Along a vertical path, the optical depth due to remnant dust settling in the atmosphere is τ vert ∼ f d (P/g)πs 2 /m d 2.9 × 10 3 f d (P/1 mbar)(g/g J ) −1 (s/1 µm) −1 . Hence, at the P ∼ 1 mbar level a mass fraction f d > ∼ 3.4 × 10 −4 is needed for vertical optical depth unity, which would require dust accretion rates far larger thanṀ d,iss .
While the Epstein formula applies, the main contribution to the optical depth along a straight-line path through the atmosphere during transit occurs at the point of closest approach, where the effective path length for an isothermal atmosphere is 2πH b R p . The resultant extinction optical depth along a path with maximum pressure P is In the last step the dust size was expressed in terms of β using Equation 2. An extinction cross section 2πs 2 has been used, which assumes that 2πs > ∼ λ.
Equation 51 shows that the optical depth during transit may be order unity at ∼ mbar pressures. By using β = 0.5, the dust mass is efficiently used to give a large area to block starlight. Grain size distributions with most of the mass at sizes much larger than a micron will give smaller optical depths by comparison. For the small, medium and large planets, the ratio of (GM p /R p ) −3/2 is 20 : 1 : 0.1, which favors lowmass planets. This is offset partially by the decreased fraction f pl accreted onto smaller mass and radius planets, but reinforced by the larger rem for small planets. Lastly, for all other parameters fixed, smaller stellar mass may lead to significantly higher optical depth, at constant β, as the blowout size is smaller. For instance, the value of τ for a K star such as HD , with M s = 0.85 M and L s = 0.33 L , will be larger by a factor (0.84/0.33) 2 = 6.6.
At lower altitudes, the higher gas density may lead to < ∼ s, in which case Stokes drag applies with where η dyn 1.3 × 10 −4 g cm −1 s −1 T/10 3 K 1/2 is the dynamic viscosity for a H 2 dominated atmosphere using the hard sphere cross section. In this case the settling speed varies slowly with altitude, giving rise to a roughly constant number density and hence decreasing mixing ratio with depth. The pressure level at which the settling speed switches from the Epstein to Stokes limit is P sw (s) = 0.7 bar 1 µm s T 10 3 K .
The Epstein formula applies for small dust down to deep in the atmosphere, while larger dust switches to the Stokes formula higher in the atmosphere.

dust distribution including settling and vertical mixing
In Section ??, the dust was assumed to settle below the source. When vertical mixing is included, the dust abundance may be significantly altered below the dust homopause where K > ∼ u set H b . The dust abundance will be determined from solutions of the equation with constants K, H b and Φ(s). The settling speed u set interpolates between the Stokes and Epstein limits. Equation 54 is numerically integrated for each dust size s from a base pressure of P = 10 bar up to P = 1 µbar. The dust flux is taken from Equation 46, and self-consistently includes the incident size distribution of the remnant dust, G85 or single-size, as well as rem for the given planet size and distance from the star. The accretion rate onto the planet is taken to be f plṀd = 10 8 g s −1 in all results shown. Since Equation 54 is linear, the optical depth can be scaled up or down for different accretion rates.
The boundary condition N = 0 is used at the base. This represents destruction of dust in deeper layers where the temperature goes above the dust melting temperature. With this boundary condition, the effect of mixing is to decrease the dust density below the K = u set H b layer, as dust is mixed down to the base where it is destroyed. This boundary condition neglects dust production through upward mixing and condensation, and focuses solely on meteoric source and the remnant dust settling down from above.
In order to include the reduced extinction cross section when 2πs < λ, the cross section used is σ(s, λ) = 2πs 2 min(1, (2πs/λ) 4 ). This will lead to decreasing τ when the bulk of the dust is smaller than the wavelength. Figure 21 shows optical depth versus wavelength, with the different lines representing pressure levels from P = 1 µbar − 1 bar. The G85 dust size distribution is used, as well as a medium-planet size. The distance from the star is a p = 10 R . The three panels show results for mixing coefficients K = 10 2 , 10 6 and 10 10 cm 2 s −1 .
In the top panel of Figure 21, the small values of K lead to mixing being unimportant even for small dust, and optical depth increases proportional to pressure as in Equation 51. The maximum optical depth achieved is τ ∼ 0.1 at the P = 1 bar level. In the middle and lower panels, increased mixing decreases the dust abundance at deeper levels in the atmosphere leading to a smaller maximum optical depth. Figure 22 shows the results for the single-size distribution with β = 0.5. The change from small to large dust compared to the wavelength is now more apparent.
The key difference is in the maximum optical depth. As the same dust accretion rate is now in smaller dust, the number density is higher, leading to significantly larger τ ∼ 10 2 at the 1 bar pressure level, for the case of small mixing. Again for larger mixing the maximum optical depth achieved is much smaller at the deeper levels.

DISCUSSION
The calculations presented in this paper motivate that accretion of interplanetary dust may occur at much larger rates for close-in, massive exoplanets than for the Earth, and that this meteoric source of heavy elements in the upper atmosphere could in principle contribute significantly to atmospheric opacity. One would like to be able to estimate the contribution from both the deep atmosphere and the meteoric source to the opacity for any specific planet. There are, however, two main sources of uncertainty that prevent firm conclusions about the importance of dust accretion at this time.
The first uncertainty concerns mixing. As discussed in Section 5, the level of mixing must be sufficiently weak in order for the meteoric source to be important, while for strong mixing the meteoric source is overwhelmed by the reservoir of heavy elements deep in the planet. Given that the region of interest is not undergoing thermal convection with the associated turbulent mixing, mixing occurs by the atmospheric circulation patterns and estimates of K are more complicated and uncertain (e.g. Zhang & Showman 2018). More accurate calculations appropriate for the region of interest likely require resistive MHD simulations including a conductivity profile with dayside photoionization in the upper atmosphere, and are beyond the scope of this paper.
The second uncertainty is observational. The accretion rate onto the planet has been written as a product of three factors, with f pl ablṀd for gas and f pl remṀd for remnant dust. For distant debris disk sources and for dust densities low enough that dust collisions are infrequent on the inspiral time, the efficiencies of accre-tion onto the planet f pl and amount ablated into gas or remaining as dust can be estimated for any particular planetary system. The dust size distribution and the accretion rate of dust toward the star,Ṁ d , are however poorly constrained. One might have hoped that emission from the outer debris disk could be observed, and henceṀ d near the star estimated. However, Absil et al. (2021) find no strong correlation between the presence of a near infrared excess for the hot dust and a mid or far infrared excess for an outer debris disk. Hot dust excesses are not necessarily accompanied by warm or cold dust excesses. This raises the issue of the origin of the hot dust. Presumably the debris disk cannot be in the inner solar system as its collisional lifetime is then much shorter than the age of the star (Wyatt et al. 2007). But if the debris disk is very far from the star, a mechanism must deliver the dust to the small radii where it is observed (e.g. Bonsor et al. 2012), and the dust level in the outer disk must remain small enough to remain undetected. How many zodis of dust are required to affect the transmission spectrum for gas-phase elements? If heavy element abundances near solar are required, then Equation 25, Table 1 and Figure 11 show that f pl ablṀd ∼ 10 8 g s −1 is required for a medium-size planet. Orbital radius a = 8 R and dust with β = 0.1 would have f pl abl 0.1. Hence the required dust accretion rate to achieve solar abundance would beṀ d = 10 9 g s −1 ∼ 100Ṁ d,iss (100 zodi's). While significantly more than what is found in our own solar system, it would likely be undetectable as a H-or K-band excess by an order of magnitude. Implicit in the above argument is that dust density and thermal emission is proportional to accretion ratė M d . If the collision time is shorter than the PR inspiral time, then destructive collisions may then limiṫ M d < ∼ a few × 10 8 g s −1 for a solar type star and grains near the blowout size (Wyatt 2005;van Lieshout et al. 2014). If true, then this maximum accretion rate may limit the abundances due to accreted dust to be subsolar by a factor of 10 or more. Such small accretion rates would also likely not deliver dust to near the sublimation zone at a high enough rate to explain the observed near-infrared excesses (Absil et al. 2021).
While there are examples of systems with both a nearinfrared excess and also a close-in planet (e.g. HD 4113 and HD 20794), we are not aware of any that are transiting. Hence, at present there is no system where transmission spectroscopy could be used to see if dust excess directly translates into an affect on the line or continuum spectra. Bonsor et al. (2018) calculated the fraction of dust hitting the planet and showed that it can be order unity for massive planets near the star. The results presented in Section 2.1 confirm this central result. As Bonsor et al. (2018) did not include a finite size star, significant differences are expected for planets in close orbits as the dust pericenter may lower to the point that it hits the sublimation zone around the star, giving large values of f st . In addition to the orbit integrations, modeling of each scattering event by theÖpik approximation (see the Appendix) shows qualitative and rough quantitative agreement with the orbit integrations, and is a much faster method. Lastly, the present study also computed the impact speed and angle as the dust collides with the planet's atmosphere.
While the orbit simulations in Section 2.1 included the radiation pressure and PR drag due to the star's light, the contribution from the planet's light was neglected for simplicity. This is a good approximation beyond a few planetary radii away from the planet, which is adequate for most of the dust particle's evolution. But since the outgoing flux from the planet is comparable to the incoming flux from the star, when very near the planet, its flux should be included. Effectively the planet will have a mass M p (1 − β p ), where β p is the ratio of radiation pressure to gravitational force for the planet, and the results from a lower-planet mass model could be used as a first approximation. To model this well, the position-dependent scattered optical light and thermal infrared emission on the dayside and nightside should be included, as well as the difference in β for the different wavelengths of the stellar and planetary flux. This is beyond the scope of the present paper.
Likewise, for grains in the upper atmosphere either due to dust accretion or upward mixing, the slower settling speeds implied by radiation pressure β p < ∼ 1, or even radiative levitation for small grains with β p > ∼ 1, has been ignored. This issue deserves further study. Lavvas & Koskinen (2017) considered stopping and ablation of dust in exoplanet atmospheres using the equations from Moses (1992). The latter study included thermal evaporation but not sputtering (see Section 3), which was found here to be a significant source of ablation for the large impact speeds implied for massive planets and close orbits. Sputtering gives rise to a distribution of ablated dust mass extending to lower pressures than the P ∼ 1 µbar found for thermal evaporation.
Another key aspect of Lavvas & Koskinen (2017) was to include radiative heating and cooling for grains. Their numerical results show that grain temperatures are close to the gas temperature at P > ∼ 10 mbar and are determined by the radiation field at lower pressures, where differences are found between the temperatures of large and small grains in that study. The analytic description in Section 3.2 uses the radiative absorption efficiency for blackbody radiation from Draine (2011) to show the effect of super-heating of small grains. When the grain size is smaller than λ/2π for its emitted thermal radiation (s = 0.1 µm for T d = 1000 K), the dust temperature rises, as compared to T eq , due to inefficient cooling. This rise continues until the dust becomes so small (s = 30Å for T s = 5800 K) that absorption of stellar radiation becomes inefficient as well.
Although not discussed in detail here, the increased temperature for small grains may lead to significantly increased detachment rates, due to the exponential dependence on temperature, as compared to the case where gas and grain temperatures are equal (Draine 1981;Lazzati 2008). The resulting smaller nucleation rates, which would occur at low pressures when the grain temperature is set by radiation, may make nucleation at such low pressures more difficult as compared to higher gas pressure regions where the grains are not superheated.
While the focus here has been on the effect of dust on the transmission spectrum, it may also affect the dayside emission spectrum. If the vertical optical depth of grains at P < ∼ 1 mbar becomes order unity, they may mask the expected dayside molecular absorption features (e.g. H 2 O rovibrational transitions, Sing et al. 2016). Further, if this dust is superheated, its emission temperature may be significantly different than would occur from gas at the that pressure level, similar to the mid-infrared spectrum in the protoplanetary disk case (Chiang & Goldreich 1997). This paper has focused on the gas-phase atoms and remnant dust directly produced by the incident dust particles, and has ignored subsequent molecule formation and dust growth. However, the gas-phase atoms and remnant dust may augment heavy elements mixed up from below by providing seeds for heterogeneous nucleation, in addition to those produced by homogeneous nucleation (Lee et al. 2018;Powell et al. 2018), as well as heavy-element monomers that may stick to the seeds. A parameter study for heterogeneous growth of dust grains through an assumed downward flux of seeds at the upper boundary was carried out in Gao et al. (2018).
In addition, as described Lavvas & Koskinen (2017), the gas-phase atoms may form molecules that then undergo photolysis to generate soot layers up near ∼ µbar levels in the atmosphere. The mass fluxes reported there were in the range of 10 −14 − 10 −12 g cm −2 s −1 , which are comparable to the accretion rates f pl abl fṀ d /(4πR 2 p ) = 10 −13 g cm −2 s −1 (f pl abl fṀ d /10 8 g s −1 )(R J /R p ) 2 from dust accretion at rates > ∼Ṁ d,iss . The results are similar as these dust accretion rates give rise to nearly solar composition abundances of species contained in the dust, while in Lavvas & Koskinen (2017) they assume solar abundances deep in the atmosphere and strong vertical mixing.

SUMMARY
Observed near-infrared excesses around nearby stars imply large amounts of hot dust orbiting near the sublimation radius around the star. The goal of this paper was to understand if a close-in, massive planet could accrete a significant fraction of the hot dust, and whether or not this could affect atmospheric abundances of gasphase atoms or dust enough to be detected through the transmission spectrum. The calculations and results are as follows.
1. Three-body simulations of the star, planet, and dust particle including PR drag were carried out to compute the fraction of inspiraling dust that is accreted by the planet. Under the assumption that the dust orbits have circularized by the time they reach the planet, and that the dust and planet orbits are nearly coplanar, it is found that Jupitermass planets near the star may accrete a significant fraction, up to 10-50%, of the dust for favorable parameters.
2. Calculations for the stopping and ablation of the dust particles in the planet's atmosphere were carried out to determine the distribution of deposited gas-phase atoms and remnant dust particles with altitude. Both thermal evaporation and sputtering are included, and two extreme dust size distributions were examined. It is found that the large impact speeds for planets close to the star or for massive planets tend to cause dust to be efficiently ablated into gas-phase atoms, while small planets further from the star tend to lead to less ablation and more mass in remnant dust. The duststopping layer is at low pressures P < ∼ µbar, where the gas temperature may be above the dust melting temperature. However, the dust temperature in these regions is set by radiation, and is regulated to be near the equilibrium temperature for large dust.
3. The abundance of gas-phase, neutral atoms below the meteoric source is computed assuming that the particles descend at their terminal velocity. For dust accretion rates roughly 10-100 larger than that in our inner solar system, the heavy element abundances in the upper atmosphere can approach the solar value. This is a central result of the paper.
4. It is argued that the increase of electrical conductivity and Lorentz drag with altitude may lead to significantly smaller flow speeds and mixing in the upper atmosphere.
5. For atomic resonance lines, the dependence of the transit radius on the level of mixing and dust accretion rate is exhibited for a range of each parameter. In the absence of the meteoric source, the transit radius decreases as the level of mixing decreases. Including the meteoric source, even for no mixing, the transit radius can approach the meteoric source altitude even for dust levels of order one zodi. As the dust accretion rate increases, the transit radius is increased further from line center. These results motivate that, for sufficiently large dust accretion rate and sufficiently weak mixing, that large transit radii may reflect the altitude of the dust source.
6. The continuum opacity of the remnant dust particles is maximized when the size distribution is peaked in near the blowout radius and mixing is weak. Optical depths of order 0.1 at mbar pressure levels may then be achieved for dust accretion rates a factor of 100 higher than in our solar system.
this time. The relative velocity vector u entering the Hill sphere can be constructed from the orbital parameters up to the signs of the radial and vertical component, which are chosen randomly. The impact parameter b ⊥ u is chosen uniformly in area inside the Hill radius and uniformly in angle. If the minimum separation with the planet r min (b) = (b 2 /b 0 )/(1 + 1 + (b/b 0 ) 2 ) is smaller than R p then a physical collision has occurred and the simulation ends. If there is no physical collision, then in the planet's reference frame, the velocity kick given to the dust particle is where the impact parameter for 90 • scattering is b 0 = GM p /u 2 . The first term is the parallel kick, slowing the particle down relative to its incoming trajectory, and the second term is the perpendicular kick along the impact parameter b. This velocity kick rotates u while keeping the same speed before and after the collision. Given the post-encounter velocity vector u + ∆u, it is checked if the particle is ejected to infinity, and if so the simulation is stopped. If the particle is not ejected, the new values of a, e and cos(I) are reconstructed from the position, which is assumed not to have changed, and the post-encounter velocity. After each PR step and each kick, it is checked if a(1 − e) < R sub , and if so the particle has hit the dust sublimation radius and the simulation is stopped. Further it is checked if a(1 + e) < a p , so that no further encounters with the Hill sphere will occur, and the simulation is stopped as the particle will eventually hit the star due to PR drag. With the focus on orbits close to the star, only a small number of encounters may occur before the dust fate is sealed, and hence each encounter is simulated. For more distant orbits for which many encounters happen before ejection or collision, the Fokker-Planck approximation may be more efficient (Yabushita 1980;Duncan et al. 1987;Tremaine 1993;Brasser & Duncan 2008). In the present case though, even when cometary diffusion of a and e occurs, diffusion toward ejection must compete with diffusion toward hitting the star, which would present an additional boundary condition that is complicated to enforce. Simulating each Hill sphere encounter is efficient for the cases considered here.   Figure 2. A single initial value u ∞ /v p = 2.5 was used forÖpik's method. The range of a p has been extended out to 40 R to see trends more clearly. There is both qualitative and rough quantitative agreement with the full numerical simulations, at the ∼ 10% level for this particular u ∞ . While not shown here, other values of u ∞ = 2.4 − 2.8 are qualitatively similar and show agreement to 10-20%. The true answer, as shown in Figure 5, is that particle leave the resonance with a range of u ∞ and so an average should be performed over u ∞ to try to obtain a more precise answer. That is beyond the scope of the present paper, where rough quantitative agreement is sufficient.