Issue 
A&A
Volume 654, October 2021



Article Number  A113  
Number of page(s)  19  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202040151  
Published online  19 October 2021 
Thermal radiation pressure as a possible mechanism for losing small particles on asteroids
^{1}
Department of Physics and Astronomy, Seoul National University,
Gwanakro 1,
Gwanakgu,
Seoul
08826,
Republic of Korea
^{2}
SNU Astronomy Research Center, Department of Physics and Astronomy, Seoul National University,
Gwanakro 1,
Gwanakgu,
Seoul
08826,
Republic of Korea
email: ishiguro@astro.snu.ac.kr
Received:
17
December
2020
Accepted:
23
July
2021
Context. Recent observations of dust ejections from active asteroids, including (3200) Phaethon, have drawn considerable interest from planetary astronomers studying the generation and removal of small dust particles on asteroids.
Aims. In this work, we aim to investigate the importance of thermal radiation pressure from asteroid regolith (AR) acting on small dust particles over the surface of the AR. In particular, we aim to understand the role of thermal radiation in the nearSun environment.
Methods. We describe the acceleration of particles over the AR within the radiation fields (direct solar, reflected (scattered) solar, and thermal radiation) in addition to the asteroid’s rotation and gravitational field. Mie theory is used because the particles of interest have sizes comparable to thermal wavelengths (~1–100 μm), and thus the geometric approximation is not applicable. A new set of formalisms is developed for the purpose.
Results. We find that the acceleration of particles with spherical radius ≲1 μm to ~10 μm is dominated by the thermal radiation from the AR when the asteroid is in the nearSun environment (heliocentric distance r_{h} ≲ 0.8 au). Under thermal radiation dominance, the net acceleration is towards space, that is, outwards from the AR. This outward acceleration is the strongest for particles of ~1 μm in radius, regardless of other parameters. A preliminary trajectory integration using the Phaethonlike model shows that such particles escape from the gravitational field within about 10 min. Our results are consistent with the previous observational studies on Phaethon in that the ejected dust particles have a spherical radius of ~1 μm.
Key words: minor planets, asteroids: general / minor planets, asteroids: individual: 3200 Phaethon / meteorites, meteors, meteoroids / interplanetary medium
© ESO 2021
1 Introduction
Recent discoveries of active asteroids (reviewed in, e.g., Jewitt 2012) and the paucity of small dust particles on the Hayabusa2 target asteroid (162173) Ryugu (Grott et al. 2019) have drawn considerable interest from planetary astronomers studying the generation and removal of small (submillimeter; subm(mm) dust particles. Centimeter(cm)sized particles were observed being ejected from the asteroid (101955) Bennu (Lauretta et al. 2019), and later their orbits were physically modeled (McMahon et al. 2020). Also, the nearSun asteroid (3200) Phaethon, the target body of the JAXA DESTINY^{+} mission (Arai et al. 2018), is a symbolic asteroid, as recurrent dust ejections have been witnessed through telescopic observations around its perihelion over the past decade (Jewitt & Li 2010; Jewitt et al. 2013; Hui & Li 2017).
From observations, the effective radius of the particles ejected from Phaethon is estimated to be as small as ~1 μm (Jewitt et al. 2013; Hui & Li 2017). Although the generation mechanism of such small particles has been thoroughly studied in terms of thermal fracturing in an astronomical context (Jewitt & Li 2010; Delbo et al. 2014; Molaro et al. 2017), the mechanism responsible for accelerating and ejecting these dust particles around the perihelion remains unexplored.
In this context, the first question to tackle is that regarding the unknown mechanism responsible for ejecting such particles from asteroids. The mechanisms on Bennu and Phaethon may likely differ because of the large size difference of the particles (cm and μm) and the fact that Phaethon showed recurrent activity near its perihelion. Phaethon is an important object because it is the only object with detectable recurrent dust ejection. There have been questions such as whether Phaethon’s dust ejections near perihelion can explain the total mass of the Geminid meteor stream (Blaauw 2017; Ryabova 2017; Jewitt et al. 2019). Such questions can be discussed more systematically once the dust ejection mechanisms are known. The aim of this work is to find a mechanism able to explain many of the observations of Phaethon’s activity, at least to some degree.
When an asteroid comes close to the Sun, the surface is heated by solar radiation. At a heliocentric distance of 0.2 au, the maximum surface temperature reaches ~900 K. Along the terminator on a nearSun asteroid, the rapid heating or cooling may induce severe stress on the surface material (Molaro et al. 2017), which will lead to thermal fracturing, as demonstrated in laboratory experiments (Delbo et al. 2014). Accordingly,small dust particles may be generated on the asteroid surface (Jewitt & Li 2010). Meteorites are the most direct sources withwhich to infer the size of such small particles. Studies have found meteoritic matrix particles with sizes down to 0.01–5 μm (Scott & Krot 2005; Righter et al. 2006), which suggests that dust particles of ~ 1 μm or smaller could be present on asteroids.
Applications of thermal radiation pressure in similar contexts are presented in McMahon et al. (2020) (see Lauretta et al. 2019 for the activity of Bennu) and the spacecraft control (Hesar et al. 2017). However, the formalisms in those works are applicable only to large particles (geometric regime) at large initial height (typically around 10 m). Using the formalism, McMahon et al. (2020) found that particles with initial speed < 10 cm s^{1} cannot escape from Bennu.
In this work, we consider small stationary particles (radius 0.5 to tens of microns) at a small height of the order of 1 cm above the asteroid regolith (AR) without any initial speed. Our results demonstrate that neither initial thrust nor large initial height are essential for such small particles to be accelerated and escape into space in the nearSun environment.
2 Methods
A particle floating above the AR experiences the following forces: direct solar radiation, indirect solar radiation reflected (scattered) on AR, thermal radiation from AR, and the gravity of the asteroid. The acceleration vectors associated to these are denoted a_{⊙}, a_{ref}, a_{ther}, and a_{grav}, respectively. In addition to gravity, the centrifugal term (a_{cen}) is included if we choose a noninertial rotating frame on the asteroid. The net radiative acceleration is defined such that a_{rad} = a_{⊙} + a_{ther} + a_{ref}. The net acceleration is the sum of all accelerations, a_{total} = a_{rad} + (a_{grav} + a_{cen}). Here, a_{ther} and a_{ref} are directed outward from the AR, while a_{⊙} is usually directed inward. The timescale of interest, for example the time required for a particle to leave the asteroid’s gravitational field, is assumed to be short compared to the asteroid’s orbital timescale, and hence the orbital motion of the asteroid is not considered.
2.1 Problem statement
In this work, we aim to calculate the net acceleration acting on a particle that is placed above the AR at an initial height of H. The way in which we achieve this initial condition is discussed in Sect. 4.3; it is taken as granted for the modeling in this work. To physically handle this problem with manageable complexity, the following assumptions are used.
The particle is a homogeneous sphere of radius r_{a} with a known refractive index. The condition of height, H ≫ r_{a}, should be satisfied all the time. Then any radiations from AR can be approximated as a plane wave (similar to solar irradiation), ignoring the curvature of the particle. Therefore, the Mie theory (Lorenz–Mie theory) is applicable.
Under the Mie scattering theory (see Bohren & Huffman 1998, chapter 4), the effective cross section for radiative pressure is calculated for a given particle size, refractive index, and wavelength (Sect. 2.3). If the spectrum of the incident radiation is known, the radiation pressure and acceleration are calculated by integrating the contribution from all wavelengths.
Throughout this work, AR is assumed as a Lambertian scatterer. Although more comprehensive models are available (e.g., Hapke 2012), they would introduce multiple new free parameters, while the Lambertian scattering is a good approximation to the complicated scattering process.
The calculation of thermal radiation from AR requires the temperature distribution on the AR to be known a priori. For this, a smooth spherical asteroid with a 1D thermal conduction model is adopted to calculate the temperature of the AR (see Sect. 2.2).
By combining the assumptions mentioned above, a simple set of descriptions of forces is established. In the following sections, we describe how the AR surface temperature T_{S} is obtained using thermal modeling, how the Mie theory is applied to calculate the radiative accelerations, and a detailed set of derivations of the formalisms used throughout this work. Justifications of selected parameter values used in this work are then provided. Finally, we outline a simple way to compare the magnitudes of solar and thermal radiations based on the developed formalism.
2.2 Thermal modeling
The surface temperature of an asteroid is necessary for estimating a_{ther}. The surface temperature is calculated by the socalled thermophysical model (TPM). In this section, we describe the outline of the widely used scheme (Spencer et al. 1989; Delbo 2004; Mueller 2007) adopted in this work. A reparameterization of the problem is then summarized; this is used throughout the present work.
In the scheme, the following classic 1D thermal conduction equation (Fourier’s equation) on the surface is solved under proper boundary conditions: (1)
Here, ρ, c_{s}, and κ denote the material mass density, the specific heat at constant pressure, and the thermal conductivity, respectively. T is the temperature of the topmost cell in the simulation (we refer to this as the AR surface temperature), t is the time, and z is the depth (Spencer et al. 1989; Mueller 2007). This conduction equation is solved by two boundary conditions: (i) the energy at the AR surface is balanced between the incident solar radiation energy, the subsurface conduction energy, and the radiative emission energy; and (ii) the temperature gradient vanishes at infinite depth (). The first boundary condition (i) is solved via the forward timecentered space scheme following the abovementioned references. The emissivity used in (i) must match with the emissivity ϵ_{S} (λ) to be used in the a_{ther} calculation below. The “infinite depth” in the second condition (ii) is set to ten times the thermal skin depth , where P_{rot} is the rotational period of the asteroid (following, e.g., Spencer et al. 1989; Mueller 2007). The Sun is assumed as a point source, although it can be as large as ~2.5° at 0.2 au.
The modeled asteroid is assumed to be a smooth sphere; that is, any roughness and largescale geological structures are ignored. There are a few advantages of ignoring roughness. First, the model is mostly described by the thermal parameter (Eq. (6)) except for the heliocentric distance and pole orientation. Thus, the discussions become analytically simple. An example is the effect of the rotational period and thermal inertia, as described below Eq. (6). Second, a direct comparison with other theoretical studies such as Ravaji et al. (2019) is possible as they also used smooth spherical asteroids (see Sect. 4.3). Consideration of roughness (selfheating, shadowing, and/or 3D conduction) is computationally demanding and increases the dimensionality of the parameter space enormously.
The rotational period P_{rot} is assumed to be significantly shorter than the orbital one. Thermal equilibrium is therefore reached at each fixed heliocentric distance, and the seasonal effect is neglected. To calculate the AR surface temperature, an initial condition of T(z) = T_{eqm}e^{−z∕lS} is set for depth z at the local noon for a given latitude, where T_{eqm} is defined below in Eq. (2). This array is then evolved for rotations, and the rotation is further repeated until either the specified temperature stability of ~ 10^{−5} K is reached at the surface or the iteration reaches 5000. It is confirmed that changes the surface temperatureby ≪0.1 K, given . The resulting temperature along the longitude and depth during the last rotation becomes the temperature of the AR. If the Sun does not rise at a given latitude, the surface temperature is fixed to 0 K for all longitudes. The temperature is calculated for patches with a resolution of1° in both latitude and longitude. From these gridded values of the surface temperature, the temperature at an arbitrary position on the asteroid is calculated by linear interpolation (RectBivariateSpline of scipy without smoothing factor).
Equation (1) implies that there exists a set of effective ρ, c_{s}, and κ – and sometimes including the porosity factor – that describes the realistic regolith. It is further assumed that ρ, c_{s}, and κ, and thus thermal inertia (), are independent of temperature (see Sect. 4.2 for further discussion). The thermal inertia has the unit of tiu ≡J m^{2} K^{1} s^{1∕2}. The notation “tiu” stands for “thermal inertia unit”, first coined by Putzig (2006).
All equations used in TPM are then finally solved by specifying the spin vector and two parameters. The first parameter is T_{eqm}, the equilibrium temperature at the subsolar point of a nonrotating asteroid or a plate: (2)
The variables used in this equation are explained as follows. First, (3)
is the Bond albedo of an asteroid. G is the slope parameter under the IAU H, G magnitude system and p_{V} is the geometric albedo in V band (first definition by Bowell et al. 1989, corrected by Myhrvold 2016). Also, (4)
is a constant for an AR model, and is the spectrumaveraged hemispherical emissivity under a blackbody spectrum of temperature T_{S}, defined in Eq. (10). Conventionally, is fixed as a constant. It is fixed to 0.90 in this work. As A_{B} ≪ 1 and , the constant c_{T} ≈ 1. The temperature (5)
is the T_{eqm} of a blackbody ( and A_{B} = 0 from Kirchhoff’s law) at 1 au. The second parameter describing the TPM result is the thermal parameter (Spencer et al. 1989; Mueller 2007): (6)
The TPM calculates the profile of a dimensionless quantity, T∕T_{eqm}, on the AR over time (rotation). This profile is uniquely determined once the spin vector and Θ are specified for the smooth model; converting it to the absolute temperature then requires T_{eqm}. An advantage of assuming a smooth, spherical asteroid is that Θ is proportional to . Therefore, changing P_{rot} by a factor of 0.5, for example, is identical to changing Γ to because both Θ and T_{eqm} are not changedby the modification (this fact is used in Sect. 3.4). Only the centrifugal acceleration must be updated according toP_{rot}.
2.3 Mie calculation
Mie theory requires (A) a homogeneous spherical particle and (B) plane parallel incident electromagnetic waves (see, e.g., Bohren & Huffman 1983). It then enables us to calculate the radiation pressure coefficient once the refractive index is known. Here, stands for the information on the particle, the refractive index. The radiation pressure experienced by the particle is proportional to Q_{pr}σ, as well as the incident irradiance. If the extinction and scattering coefficients (Q_{ext} and Q_{sca}, respectively) are given, then Q_{pr} is calculated as (7)
where ϑ is the scattering angle. The second term (Q_{sca} multiplied by the socalled asymmetry factor, ⟨cos ϑ⟩) takes only the momentum transferred towards the propagation direction (see e.g., Sect 2.3 of van de Hulst 1981 and Sect. 4.5 of Bohren & Huffman 1983). Each term in the equation can be calculated once the particle radius r_{a} and refractive index are found using Mie theory. For large particles (r_{a} ≫ λ), a geometric limit is expected: Q_{pr} → 1.
For the optical indices, two wellstudied materials are considered, namely olivine and magnetite. The former is transparent and the latter is opaque material in optical wavelengths. The optical indices are directly adopted from previous experimental works: lowFe olivine (Mg_{1.9}Fe_{0.1}SiO_{4}; Fabian et al. 2001^{1}) and magnetite (Fe_{3}O_{4}; Triaud, A. H. M. J.^{2}; Mutschke, H., 2021, priv. comm.). Although they are not representative materials in the real AR compared to mixtures like meteoritic samples, our model can easily be extended to any real material with a known refractive index. These materials are chosen because the refractive indices are accurately known. The calculated Q_{pr} (λ, r_{a}) values for these two materials are given in Appendix A.
For the requirement (B) to hold, an assumption H ≫ r_{a} is needed such that the solid angle of the particle viewed from an infinitesimal AR patch () is much smaller than unity. As shown in Sect. 3, the r_{a} of interest does not exceed a few tens of microns, and hence an arbitrary choice of H = 1 cm satisfies the requirement.
2.4 Radiative accelerations
To derive the formalism for radiative accelerations, we consider an arbitrary spectral irradiance (J_{λ}) from a black body of temperatureT, spectrumaveraged emissivity , and solid angle dΩ on the particle from a certain incident direction (Fig. 1). For simplicity, we define a few auxiliary quantities (Appendix B.1): (8) (9) (10)
Here, B_{λ}(λ, T) is the Planck function, σ_{SB} is the Stefan–Boltzmann constant, and ϵ(λ) is the emissivity function of the radiating source (either the Sun or the AR). Then, in Eq. (8) is the spectrumaveraged Q_{pr} value, b_{λ} in Eq. (9) is the normalized Planck function, and is the spectrumaveraged emissivity. Using these parameters, the magnitude of the radiative acceleration on a Mie particle of radius r_{a} and mass density ρ_{m} is (Appendix B.2): (11)
We use the Q_{pr} values calculated in Appendix A to integrate in Eq. (8). Figure 2 shows the calculated for magnetite and olivine whenϵ, and therefore , are constants. The trend that as r_{a} ≫ λ (i.e., the geometric optics regime) is clearly seen in Fig. 2. For olivine, the optical index is measured separately for three polarization directions with respect to the crystalline axis orientations. (Fabian et al. 2001; see Appendix A for details). The calculated values in Fig. 2a are, however, nearly the same for all three cases, at least under the temperatures in which we are interested (AR temperature of T_{S} ~ 1001000 K and T = T_{⊙} = 5777 K). Thus, the average of the three values is adopted throughout this work. Magnetite does not have multiple optical indices, because of its simplesymmetric (facecentered cubic) structure.
The values are calculated in the range of r_{a} ∈ [0.1, 1000] μm (at 0.1 μm intervals from 0.1 to 2.0 μm, 1 μm intervals to 10 μm, 10 μm intervals to 100 μm, and 100 μm intervals to 1000 μm) and T_{S} ∈ [10, 1300] K (with 1 K interval) for both magnetite and olivine. The values are then linearly interpolated when necessary. Figure 2 shows values with respect to r_{a} for fixed T.
Fig. 1 Schematic diagram of a spherical particle hovering over the surface. The particle of radius r_{a} hovers over the planar regolith at the height of H(≫ r_{a}). The solar irradiation (the long red arrow) is originated from the direction . A tiny surface patch with an area of dA reflects the solar radiation and emits the thermal radiation to the particle (the dashed and squiggly arrows, respectively, to the directionof ) with the emission angle of e′. In our model, e″ = e′ holds because of the planar surface assumption. The isothermal regolith region (gray dashed) is assumed to be a circle with a radius of r_{0}(≫ H ≫ r_{a}). No radiation outside this region is considered because the patch right below the particle already fills most of the field of view (r_{0} ≫ H). 
Fig. 2 Calculated Q_{pr} value averaged over blackbody spectra of different temperatures ( in our notation). The values are plotted as a function of the particle size r_{a} and blackbody temperature T of the incident radiation for two minerals: a) olivine and b) magnetite. The emissivity ϵ is assumed to be constant, so cancels out in Eq. (9). For olivine, three polarization directions with respect to the crystalline axes can lead to small differences in the results (three different markers; see Appendix A), and the average value (solid line) is used in this work. converges tounity for a large particle limit, as predicted from the geometrical optics approximation. 
2.4.1 Solar radiation a_{⊙}
The irradiance from the Sun is directed inward to the AR surface with an incident angle i. If it is assumed to be a uniform beam from the solid angle of , Eq. (11) gives (12)
The solar radiation is assumed as a perfect blackbody in this work, so . a_{⊙} decreases as as the solar irradiance decreases. For the Sun, we use the notation .
Because a_{⊙} is oriented in the antisolar direction, a_{⊙} is in the inward direction on the dayside of the asteroid. When the particle goes behind the shadow of the asteroid, we regard a_{⊙} = 0.
2.4.2 Reflected solar radiation a_{ref}
For a reflected (scattered) component of the acceleration, the radiative acceleration contributed from each infinitesimal surface patch in the field of view must be integrated. After proper integration (Appendix B.3), the acceleration from a Lambertian scattering AR is derived as (13)
where μ_{i} := max {cosi, 0} for the incidence angle i, and the height factor (14)
The r_{0} is set as 1% of the asteroid’s radius throughout this work. Here, but is almost unity at the initial condition, as r_{0} ≫ H = 1 cm. We did not assume in this general formula, because it plays an important role in the trajectory integration (Sect. 4.4) when H ≳ r_{0}. Because A_{B} ≪ 1 for many asteroids, one obtains a_{ref} ≪ a_{⊙} as expected.
2.4.3 Thermal radiation a_{ther}
Similar to a_{ref}, an integration must take place over the surface patches in the field of view (Appendix B.4): (15)
From the formula, we see that a_{ther} is very sensitive to the AR temperature T_{S}. Also, it is linearly proportional to , which implies the importance of the adequate Mie calculation, especially for small particles (when deviates significantly from unity).
2.4.4 Comparing a⊙ and a_{ther}
To compare the zaxis components of a_{⊙} in Eq. (12) and a_{ther} in Eq. (15), we reformulated a_{ther} using Eqs. (2) and (12): (16)
where is the relative efficacy of the thermal radiation from the AR compared to the solar radiation. In other words, it can be understood as the ratio of effective crosssectional areas of the particle exposed to the spectrum of temperature T_{S} and T_{⊙}. The ratio of the two acceleration mechanisms along the zaxis is then given by (17)
Here, a_{⊙,z} = a_{⊙}μ_{i}.
If the thermal inertia Γ = 0, the surface temperature is the same as the equilibrium temperature, T_{S} = T_{eqm}. This is the case for the nearEarth asteroid thermal model (NEATM; Harris 1998) for example, with the beaming parameter of unity. Thus, the ratio in Eq. (17) stays constant during rotation of the asteroid for Γ = 0 unless μ_{i} vanishes. In reality, however, this ratio changes over the local time. In particular, it increases towards the evening terminator. This is because T_{S} does not cool immediately like NEATM (Sect. 3.1), and so . Therefore, one would expect a_{ther} to dominate a_{⊙} in the afternoon under TPM, which is a more realistic version of NEATM. This is a key point in the motivation of the present study.
2.5 Gravity and rotation
The gravitational acceleration is calculated by (18)
where ρ_{M} is the bulk mass density of the asteroid (not necessarily equal to the mass density of grains ρ_{m}; likely ρ_{M} < ρ_{m} due to the bulk porosity), and D is the asteroid diameter. The gravitational acceleration is directed inward from the AR. Figure 5 of Jewitt (2012) can be reproduced by comparing a_{grav} in Eq. (18) with a_{⊙} in Eq. (12).
The last remaining acceleration component is a fictitious one, namely the centrifugal acceleration. This is introduced only when describing the net acceleration of the particle from an observer sitting on the AR (i.e., in the rotating frame on the asteroid). The centrifugal acceleration on the equator is calculated from (19)
for the asteroid’s rotational period P_{rot}. In the main part of this work, a_{cen} is taken into account for calculating the effective gravity by . As we are considering an initially stationary particle with respect to the observer, the Coriolis force is zero.
2.6 Parameter selection of model asteroids
For fictitious model asteroids, two aspect angles (the angle between the spin vector and the asteroid–Sun vector) are considered, namely, 90° and 45°, which are called the “perpendicular” and “aspect 45°” models in this work, respectively. It is assumed that the fiducial model asteroids are located at a heliocentric distance of r_{h} = 0.2 au, that is, in the nearSun environment. We further changed it from 0.1 to 2 au for test cases (25 equally spaced r_{h} values in logarithmic scale by r_{h} = 10^{k(log1020)∕24−1} au for integer k = 0 to 24). Fiducial asteroids have smooth spherical shapes with D = 1 km, and for test cases, D ranges from 0.1 to 80 km (50 equally spaced D values in logarithmic scale by D = 10^{k(log10800)∕49−1} km for integer k = 0 to 49). For S, B, and Ctype asteroids, the bulk mass density ρ_{M} ~ 2660 ± 1290, 2190 ± 1000, and 1570 ± 1380 kg m^{3}, respectively (Carry 2012). Therefore, as a simplifying assumption, the model asteroids have a fixed ρ_{M} = 2000 kg m^{3}. The Bond albedo is fixed at A_{B} = 0.05, a typical value for an asteroid with p_{V} ~ 0.10.2 using Eq. (3). The rotation period of the fiducial model is also fixed to P_{rot} = 6 h, a representative value of the asteroids in the size range (Warner et al. 2009, updated in 2020^{3}); the P_{rot} value has almost no effect on these preliminary results when it is changed by a factor of two (see Sect. 3.4). The thermal inertia is taken as 200 tiu (a typical value for asteroids of a similar size range in Delbo et al. 2015).
All the abovementioned parameters are summarized in Table 1 (Model asteroids columns). The surface temperatures of these asteroids are calculated by solving the 1D heat conduction equation using these parameters (Sect. 2.2).
Physical parameters of Phaethon (Hanuš et al. 2018) and the fiducial and tested parameters for fictitious model asteroids are given.
Fig. 3 Equatorial temperature profiles. a, the fictitious asteroid at r_{h} = 0.20 au and b, Phaethon at perihelion (r_{h} = 0.14 au). For the model asteroids, different aspect angles are considered: perpendicular (solid) and 45° (dashed). The different colors of the lines indicate different thermal inertia values: black lines for Γ = 0, light blue lines for fiducial value for model asteroids and Phaethon (Table 1), and red lines for the values considering the temperaturedependent thermal inertia (Γ = 669 tiu and 2622 tiu for model asteroids and Phaethon, respectively; see Sect. 4.2). 
2.7 Selection of particle parameters
The particles have mass densities of ρ_{m} = 3000 kg m^{3} following previous works on dust ejection or electrostatic dust lofting simulations (e.g., Li & Jewitt 2013; Jewitt et al. 2013; Senshu et al. 2015; Hui & Li 2017; Orger et al. 2018). The assumed mass density ρ_{m} is close to the (25143) Itokawa sample (3400 kg m^{3} Tsuchiyama et al. 2011). The particles are assumed to be homogeneous^{4} and spherical, and placed on the surface (at a height of H = 1 cm over the surface). The radiation pressure coefficient (Q_{pr}) is obtained from Mie theory (Sect. 2.3).
We do not consider microporosity (ϕ) in this work.If ϕ is to be introduced, the particle density ρ_{m} must be updated to (1 −ϕ)ρ_{m}, but the Mie theory (Sect. 2.3) still applies if the particle is still assumed to be homogeneous. The radiative accelerations are then affected linearly: a_{⊙}, a_{ref} and a_{ther} (Eqs. (12), (13), (15), respectively) have the simple identical dependency, and hence the net radiative acceleration will be scaled as .
3 Results
3.1 Temperature profile on asteroids
To test the validity of our TPM calculation, a few temperature profiles over time were generated and compared with previous works (e.g., Mueller 2007, Fig. 2.2), and the code passed the tests. Temperature profiles in Fig. 3 resemble some representative ones from fictitious model asteroids and Phaethon used in this work. The thermal inertia difference does not induce large temperature profile changes during the day compared to the night (for the physical reasons, see Sect. 4.2).
The nighttime temperature increases significantly as thermal inertia. Also, the nighttime temperatures for perpendicular and aspect 45° models are very similar, although they differ considerably in the daytime. Therefore, it is expected that would be insensitive to the aspect angle at night. The largest thermal inertias are the fiducial values corrected by the temperature dependency (Sect. 4.2).
Fig. 4 Components of the radiative accelerations acting on a r_{a} = 2 μm particle at the height of H = 1 cm on the equator of D = 1 km model asteroid as a function of longitude (local time), for a) olivine; b) magnetite. The black horizontal lines show the gravity (dashed) and the effective gravity taking the centrifugal acceleration into account (solid). The absolute value of is used as it is negative. Both a_{⊙} and a_{ref} turn off almost immediately after Sunset, which causes the spiky feature of at longitude of 270°. a_{ther} persists throughout the night thanks to nonzero thermal inertia. 
3.2 Components of acceleration
The resultant radiative acceleration components for r_{a} = 2 μm olivine and magnetite particles on the equator of the fiducial model asteroid (D = 1 km) are shown in Fig. 4. Here we see that, especially at the evening side, the thermal component, a_{ther}, dominates the total radiative acceleration, a_{rad}. At night, is larger than the threshold , that is, the particles are accelerated outwards from the AR. The difference between these curves of olivine and magnetite comes from the difference in the refractive indices (Sect. 2.3 and Appendix A), which makes different (Fig. 2).
The sharp spike of at the evening terminator is caused by the sudden quench of the solar radiation near the terminator (when the particle enters the shadow behind the asteroid). In the computation, a_{⊙} is not turned off until it reaches this shadow even though its footprint reached the longitude of 270°, but this lag is invisible due to the small height H = 1 cm. The reflected component a_{ref} turns off right after Sunset, that is, when the footprint of the particle is in shadow, regardless of the particle’s height.
3.3 Acceleration over local time
Similar to Fig. 4, the net radiative accelerations for various particle radius r_{a} on different asteroids are given in Fig. 5. The fiducial model asteroid (D = 1 km) and the Phaethonlike model (Table 1) are considered. These figures indicate that particles with radii from r_{a} < 1 to r_{a} ~ 10 μm tend to have larger accelerations outwards from AR. This holds for both our fiducial model asteroids (D = 1 km) and for Phaethon (D = 5.1 km).
As mentioned in Sect. 3.1, the thermal inertia uncertainty does not induce a large uncertainty in during the day compared to the night. The difference is small between the perpendicular (solid) and aspect 45° models (dashed), especially at night. The cause of the inward acceleration of the 1 μm olivine particle is given in Sect. 4.1.
3.4 Conditions for outward acceleration
In this section, the fictitious model asteroids are used to investigate the effects of some parameters. First, the effect of heliocentric distance r_{h} is vital to testing whether the outward dust acceleration is a phenomenon appearing universally on any asteroid or if it happens only in the nearSun environment. Also, the effect of asteroid size D and rotational period P_{rot} give hints as to whether the likelihood of losing small particles is dependent on the asteroid’s size and rotation. Here, these parameters are investigated by simulating the particles on fictitious asteroids of varying D, aspect angle, and P_{rot}, at several r_{h}. Two particle types (magnetite and olivine) are investigated for varying radius r_{a}.
The total radial acceleration at the initial state, that is, a stationary particle at height H = 1 cm above the AR, is expressed asa function of D and r_{a} at given r_{h} by . The maximum of along the equator is then found and plotted as a function of (D, r_{a}) for a few selected models (aspect angle, rotational period) and heliocentric distances, as in Fig. 6.
From the figure, the maximum of peaks at r_{a} ≈ 1 μm for magnetite and at r_{a} ≈ 23 μm for olivine. The same trend must be true for realistic AR materials if they have and ρ_{m} similar to these particles. There is a general trend that less particles can achieve outward net acceleration if r_{h} or D increases. The former is because T_{S} gets smaller as r_{h} increases (thus reducing a_{ther} in Eq. (15)), and the latter is because a_{grav} (Eq. (18)) dominates the radial total acceleration.
No large change is visible for models with different P_{rot} in Fig. 6. As mentioned in Sect. 2.2, the temperature on AR after halving P_{rot} is identical to that when Γ = 200 tiu is updated to , because both Θ and T_{eqm} are the same for these two cases. However, such a small change in Γ does not result in a large change in the temperature profile at noon in the nearSun environment (Fig. 3). As the maximum appear near noon (Fig. 5), it is expected that the change in P_{rot} by a factor of two will not result in a large change in a_{rad} and thus in a_{total}. Therefore, the change in P_{rot} does not change the results shown in Fig. 6 significantly.
Comparing the rows in Fig. 6, it is often seen that particles on P_{rot} = 6 h have slightly greater outward acceleration than P_{rot} = 3 h. This result may contradict an intuitive prediction that centrifugal acceleration of a fast rotator is larger than that of a slow rotator. A change in Θ when changing P_{rot} (as mentioned above, P_{rot} → 0.5P_{rot} has identical effect as ) slightly changes the maximum temperature as in Fig. 3. Hence, a_{rad}, which is a sensitive function of the temperature (Eq. (15)) is reduced near noon, where the maximum outward acceleration likely happens. As these two effects take place at the same time, the contours in Fig. 6 are not necessarily expected to show a simple trend as a function of P_{rot}.
Also, in Fig. 6, some bumps are seen especially at contours around r_{a} ~ 1 μm. Figure 5a is an example for understanding the olivine (blue dashed) contour in Fig. 6a. In the figure, the location where occurs changes from the evening terminator (when r_{a} ≲ 1 μm) to noon (when r_{a} ≳ 1.5 μm). Therefore, the olivine contour in Fig. 6a will change its trend at around r_{a} ~ 1 μm.
Finally, using the calculated values, we constrain the diameter of asteroids that can accelerate particles. As seen in Fig. 6, the condition for diameter gives only the upper bound under the parameter space considered in this study, and hence the maximum diameter is plotted in Fig. 7. The derived maximum diameters of asteroids are smaller for olivine than magnetite, because olivine is less affected by the thermal radiation (smaller in Fig. 2). Also, the clear trend of decreasing maximum diameter is seen as r_{h} increases, that is, as the thermal radiation weakens. The profiles in Fig. 7 approximately follow a power law , as indicated by the black lines in the figure. At r_{h} = 1 au, particles can be accelerated outward to space only if D ≲ 1 km under the model used in this study.
Summarizing the above results in this section, we find the following key points: (1) the thermal radiation becomes critical as an asteroid approaches the Sun, (2) smaller asteroids tend to lose particles more easily than larger asteroids, (3) particles of r_{a} ~ 1 μm are most efficiently accelerated outwards, and (4) these results are less sensitive to rotational period unless a_{cen}≳ a_{grav}.
Fig. 5 Radial accelerations on the equatorial surface at heights of H = 1 cm, for the corresponding particle radius r_{a} with respect to longitude. The net radial radiative accelerations () are shown.The black horizontal lines show the gravity and the effective gravity as in Fig. 4. a, c) olivine; b, d) magnetite. Panels a and b are calculated for fiducial perpendicular (solid) and aspect 45° (dashed) model asteroids (D = 1 km) at the heliocentric distance r_{h} = 0.2 au for given particle radii r_{a} in the legend. For the perpendicular model, the 25% ambiguity inthe thermal inertia (Γ = 200 ± 50 tiu) is expressed as shading. Panels c and d are calculated for Phaethon using the nominal physical parameter values (Table 1) at the perihelion. The shading signifies the uncertainty range in the thermal inertia (Γ = 600 ± 200 tiu, Hanuš et al. 2018) of Phaethon. From these figures, it is obvious that most of ≲ 10μmsized particles areaccelerated outwards from the AR. 
Fig. 6 Maximum total radial acceleration, [mm s^{2}] of particleson the equator of the fictitious model asteroids. Each row corresponds to a different asteroid model a) to d) perpendicular model with P_{rot} = 6 h, e) to h) perpendicular model with P_{rot} = 3 h, i) to l) aspect 45° model with P_{rot} = 6 h, and m) to p) aspect 45° model with P_{rot} = 3 h). Each column corresponds to different heliocentric distances. The maximum acceleration is calculated by finding the maximum value along the longitude on the equator for the given asteroid diameter D and particle radius r_{a}. The blue dashed lines and red solid lines indicate olivine and magnetite, respectively. The thickness of the contour increases as the acceleration value increases for fixed contour levels . Negative (inward) accelerations are not drawn. 
Fig. 7 Maximum diameters of asteroids (D) that producepositive on the equator as a function of heliocentric distance (r_{h} ∈ [0.1, 2.0] au). The red solid and blue dashed lines are the results for magnetite and olivine, respectively. The results of perpendicular and aspect 45° models are discriminated by the line thicknesses, and the results of different rotational periods (6 h and 3 h) are discriminated by the transparency of the lines. The black lines show the slope according to a power law. 
3.5 Phaethon
Asteroid (3200) Phaethon is the most appropriate object with which to test the validity of our model because it has shown activity near the perihelion (Jewitt & Li 2010; Li & Jewitt 2013; Jewitt et al. 2013; Hui & Li 2017) but has not shown lingering activity at r_{h} ~ 1 au (Green et al. 1985; Luu & Jewitt 1992; Chamberlin et al. 1996; Hsieh & Jewitt 2005; Wiegert et al. 2008; Jewitt et al. 2018, 2019). If the recurrent nearperihelion activities are caused by the same mechanism over time, any model explaining the phenomenon must coincide with the observations, in that: (A) Phaethon’s activity (dust ejection) must be stronger near the perihelion, and (B) the ejected particles must have a spherical equivalent radius of r_{a} ~ 1 μm (Jewitt et al. 2013; Hui & Li 2017).
To test the ability of our model to fulfil the requirements, the maximum was found, as in Sect. 3.4, along Phaethon’s equator. The physical parameters are as given in Table 1. The results are shown in Fig. 8 as a function of heliocentric distance r_{h} and particle radius r_{a} for magnetite and olivine. If the maximum on the equator can be a proxy of the likelihood of dust ejection, this figure tells us that Phaethon should show the strongest dust ejection near the perihelion. Over time, as Phaethon’s r_{h} increases (Sect. 4.3.3), the likelihood will rapidly drop, and eventually converge to zero when r_{h} ≳ 0.8 au. Therefore, our model satisfies the requirement (A).
In addition, Fig. 8 indicates that particles of radius from r_{a} < 1 μm to r_{a} ~ 23 μm are likely to have larger accelerations outwards from the AR, depending on the material. This is similar to the results shown in Sect. 3.3 (Fig. 5) and 3.4 (Fig. 6), which used fictitious model asteroids. Moreover, this is highly consistent with observational studies, as in requirement (B). Therefore, we conclude that our model satisfies the two important requirements outlined above and based on previous observational studies of Phaethon.
The two plots in Fig. 8 (pre and postperihelion) appear similar because, by coincidence, the aspect angle of Phaethon at perihelion is almost 90° (strictly speaking, 101°). This means the geometric configuration of the Sun and AR patch on the equator are nearly timesymmetric with respect to the perihelion.
Fig. 8 Maximum radial acceleration on the Phaethon’s equator. As in Fig. 6, the maximum acceleration on the equator is calculated for a wide range of particle sizes (0.5 40 μm) before (a) and after (b) the perihelion passage. The black vertical line corresponds to thePhaethon’s perihelion. The numbers on the contour plots are in units of mm s^{2}, and the thickness of the contour increases as the acceleration value increases for fixed contour levels . The blue dashed and red solid lines indicate olivine and magnetite, respectively. 
4 Discussion
4.1 The acceleration peak at ra ~ 1 μm
In this section, we provide qualitative discussions regarding the result that r_{a} ~ 1 μm particles have the largest acceleration. Larger particles tend to have negative values because gravity dominates a_{total}, as the areatomassratio decreases. The coefficient plays a key role, especially for small particles. In Fig. 2, it is seen that gets small for smaller particle sizes. This means that solar radiation dominates the thermal radiation, and hence the particle is easily pushed into the AR.
This is why particles of r_{a} ~ 1 μm tend to experience outward acceleration, while larger or smaller particles tend to experience weaker, or even inward acceleration (). For olivine, the ratio R_{Q} decreases steeply so that solar radiation dominates during the daytime. For this reason, the olivine particle of r_{a} = 1 μm in Fig. 5 displays an inward acceleration during the daytime.
Fig. 9 Identical plot to Fig. 5, except the thermal inertia of asteroids is changed according to Eq. (20). (a, c) Olivine; (b, d) magnetite. The shade in each plot signifies the uncertainty on thermal inertia, scaled by the same factor, . 
4.2 Effect of thermal inertia Γ
Delbo’ et al. (2007) and Delbo et al. (2015) discuss the effect of the temperature dependency of Γ. The dependency is caused mostly by the temperature dependency of the thermal conductivity, κ. Although Γ must change at each voxel in TPM (depending on the temperature at each position on the asteroid and depth), a firstorder approximation to cope with this temperature dependence is to roughly scale Γ by or (20)
where Γ_{1} is the thermal inertia at 1 au (Delbo et al. 2015). In the present work, the thermal inertia of the model asteroids is taken to be 200 tiu from the Γ_{1} values of the similarsized asteroids (Delbo et al. 2015, Fig. 9). If the asteroids were placed at r_{h} = 0.2 au, the true effective thermal inertia would be Γ = 669 tiu following this relationship. Phaethon’s thermal inertia is also determined at r_{h} ~ 1 au (Hanuš et al. 2016, 2018), so the value at its perihelion, r_{h} = 0.14 au, is Γ = 2622 tiu. The temperature profiles using these Γ values are plotted in Fig. 3.
Updating the thermal inertia values of the model asteroids and Phaethon to these larger values, we obtain Fig. 9. When the thermal inertia Γ is updated from the fiducial value, the morning and noon temperatures decrease because the heating process takes longer. On the other hand, the afternoon and night temperatures increase because cooling also takes longer, as depicted in Fig. 3. This implies that the a_{ther} for the models with increased Γ will be larger than fiducial models in the afternoon. This effect is confirmed by comparing Figs. 5 and 9.
4.3 Possible mechanisms to generate and lift particles to initial height
The accelerations shown in this work, such as the radiative accelerations, have a maximum order of < 100 mm s^{2}. However, this acceleration is still too small to overcome the cohesion force of the regolith particles. According to Eq. (12) of Hartzell & Scheeres (2011), the acceleration equivalent to the cohesion force is (21)
where S = 0.11 is the ‘cleaness’ parameter and the cohesion force in the original work is divided by particle mass . Thus, an event to break this cohesion force is likely to have happened to free the particle from the regolith and achieve the required initial condition H ≫ r_{a}. Two such mechanisms are briefly discussed, and then the implication to Phaethon’s orbit is discussed.
Fig. 10 Thermal fatigue plot. The thermal parameter (Θ) over heliocentric distance (r_{h}) for Phaethon (solid lines, r_{h} ≥ 0.14 au) and fiducial perpendicular model asteroids with varying r_{h} (dashedlines) are shown. The black star and plus symbols indicate that Γ is fixed, while the green triangle and cross markers show results where Γ is assumed to be changed by Eq. (20). The blue solid lines are the contours of , and the brown dashed lines are that of . The values shown for the contours are the and values for the case when c_{T} = 1 and P_{rot} = 1 h. The contours get thicker for larger values to guide the eyes. 
4.3.1 Thermal fatigue
Ravaji et al. (2019) provided a simple visualization of two (insufficient, as they mention, yet useful) criteria to investigate the probability of thermal fatigue (Fig. 5 of their work). First is where is the daily maximum temperature difference on the asteroid’s surface, that is, the amplitude of the curves in Fig. 3. The second criterion is , where is the maximum timerate of change in surface temperatures, that is, the maximum time derivative of the temperature profile. In the original work, the thermal inertia and Bond albedo of the asteroids were fixed while rotational period and heliocentric distance are changed.
We devise a similar but more general plot as shown in Fig. 10. This figure visualizes the likelihood of thermal fatigue on general small bodies of the Solar System. Taking the Bond albedo, the emissivity, and the rotational period as free parameters, the two criteria are reformulated by (22)
For a detailed derivation, see Appendix C.
The fiducial perpendicular model and Phaethon have (c_{T}, P_{rot}) = (1.014, 6 h) and (1.015, 3.604 h), and so the righthand side of Eq. (22) is 98.7 [K] and 98.6 [K], respectively.This criterion can be checked by tracing blue solid contours in Fig. 10. The righthand side of Eq. (23) becomes 11.8 [K min^{1}] and 7.1 [K min^{1}], respectively (brown dashed contours in Fig. 10). From the figure, the fiducial model at r_{h} ≲ 1.5 au and Phaethon at r_{h} ≲ 0.8 au satisfy both of the criteria, and hence, thermal fatigue is likely to occur in those regions (Region III in the notation of Ravaji et al. 2019). The probability of the occurrence of the fracture and the initial size and velocity distribution of the generated particle are beyond the scope of this work. Previous work (Jewitt & Li 2010) suggested that thermal expansion potential energy is translatedinto kinetic energy and can accelerate the particle to a certain initial velocity immediately after the thermal fatigue.
4.3.2 Electrostatic dust lofting and levitation
The electrostatic lofting mechanism, which is the mechanism that launches the particles from the AR, has recently been studied extensively both theoretically (Zimmerman et al. 2016) and experimentally (Wang et al. 2016; Hood et al. 2018). If small particles are residing on the rocks by cohesion force, a supercharging effect near the terminator or the shadow on asteroids may free the particle from the AR with initial kinetic energy and achieve the condition H ≫ r_{a}.
From Wang et al. (2016) (Fig. 4 of them) and Hood et al. (2018), particles of radius r_{a} = 522 μm (2r_{a} = 1044 μm) could achieve an initial speed of up to v_{z,0} ~ 0.6 m s^{1} due to the supercharging effect at r_{h} = 1 au. This is the aftermath of Coulomb interaction, and therefore the particle size range and initial speed are independent of asteroid size (D) or r_{h}. Thus, electrostatic dust lofting due to the supercharging effect can achieve the initial condition (H ≫ r_{a}).
It is found that the dust lofting mechanism does not repeat forever on the asteroid; rather, the dust lofting rate decreases over time during the experiment (Hood et al. 2018). Therefore, it is expected that the electrostatic dust lofting will occur only immediately after the AR is reset, that is, when fresh small particles are generated by micrometeoritic bombardment or thermal fatigue.
4.3.3 Phaethon’s orbital movement
As described in Sect. 2.2, the TPM calculates the equilibrium temperature after rotating the asteroid at least 50 times and ignoring seasonal variation. Phaethon has a large eccentricity of 0.89, and therefore the true anomaly changes dramatically near its perihelion. Figure 11 shows the change in r_{h} of Phaethon and the aspect angle using the pole orientation (Table 1). During the ten asteroidal days (≈ 1.5 day) near the perihelion, r_{h} changes by ~10% and aspect angle changes >10°. The environment of Phaethon at its perihelion is a chaos of thermal inequilibria under strong solar irradiation.
Although this weakens our assumption that the seasonal effect is negligible, such information can provide important hints as to the probability of thermal fatigue and electrostatic dust lofting. The chaotic inequilibria result in a large and sudden change in temperaturenot only along the longitude but also along the latitude due to the quick change in the aspect angle. As soon as fracture happens and small particles are generated, they may have reached an initial height or be lofted by the supercharging effect. In this work, initial velocity was set to zero with a small initial height above AR (H = 1 cm) to show the validity of our idea even when the initial condition is unfavorable. If the particles have a positive initial speed outwards, more particles could have ejected. Moreover, there is a timesymmetry in the illuminationcondition on Phaethon along the latitude: The thermal (illumination) history of the north during the perihelion passage is nearly the timereversal of that of the south because the aspect angle is nearly 90°. This may have implications on the latitudinal particle size distribution on Phaethon, but investigation of this matter is beyond the scope of this work.
4.4 Particle trajectory integration
One of the most important points to improve in this work is that the large acceleration does not necessarily guarantee the escape of particles. This problem is at least twofold. First, the particle trajectory integration must be performed for varying height of the particle, until at least the Hill radius (r_{H} = 66 km; Jewitt et al. 2019). Second, even though the particle can be ejected by thermal radiation, it must have been generated (by, e.g., micrometeorite bombardment or thermal fatigue), and the initial conditions must be known accurately. This second part is beyond the scope of this work and has already been discussed in previous sections to some extent, but we develop the first point below.
We consider the particle trajectories in a frame in which the observer orbits around the Sun with Phaethon, but does not rotate around the asteroid. Thus, the particles are set to have an initial speed of v_{0} = πD∕P_{rot} along the direction of the asteroid’s rotation. Also, the centrifugal acceleration from the orbital motion and the gravitational acceleration of the Sun are balanced. In the given reference frame, Coriolis force is negligible.
For the particles’ trajectories, the simplex integration scheme (leapfrog kickdriftkick) with a fixed time interval of d t = 0.025 s is used. For selected test cases, longer time intervals (e.g., 0.1 s) sometimes showed artifacts, but time intervals much shorter than d t = 0.025 s (e.g., 0.001 s) did not change the result. In this preliminary work, Eqs. (13) and (15) are used for all the cases. Nevertheless, it is meaningful to conduct the trajectory integration, because the purpose of this preliminary study is to reveal the possibility that small particles can escape from the asteroid.
The initial positions of particles are selected by ±15° from the morning and evening terminators at 5° intervals for r_{a} = 1 and 2 μm (Fig. 12).
A particle is flagged by three conditions:
 (A)
Escape from the gravity field: The radial speed
 (B)
Landing on the surface: The particle height gets smaller than 0.1 cm.
 (C)
Oscillation over the surface: The particle starts to fall back to the surface.
The condition (A) is applied when the particles are guaranteed to escape from the gravitational field of Phaethon. In the equation of (A), M is the Phaethon’s mass and r is the asteroidcentric distance to the particle. v_{esc} on the equator in our model simulation is 2.5 m s^{1}. There was no particle that once reached and hit the ground in our test cases. This situation may happen due to the solar radiation when the particle is ejected at high speed towards the Sun. The condition (B) is applied when the particles reach the asteroid surface. The condition (C) captures the artifact caused by our assumptions that underestimate a_{ther} and a_{ref} as the particle height increases. Hence, the particle may fail to further accelerate outwards, even though it could have escaped from the asteroid in reality. There are two ways to flag (C). First, a particle’s height is divided by an arbitrary scale and if the quotient (i.e., the integer part of ) is decreased during the trajectory integration, it is flagged (C). We set scale = 5 m. Another case is when the particle’s longitude changes from lon ∈ [0°, 180°) (morning) to lon ∈ [180°, 360°) (evening) orvice versa.
In Fig. 12, the trajectories of particles (r_{a} = 1 and 2 μm) are drawn until one of these three flags is applied. The flag is also indicated at the initial position of each particle. If a particle was initially near noon, the trajectory first goes up (increased H) but then down (decreased H) and repeats this oscillatory behavior until it escapes or hits the ground. This is because of the underestimation of a_{ther} + a_{ref}, while a_{⊙} is estimated correctly. This is why the initial positions were chosen only near the terminators in Fig. 12. All the particles with the flag (A) in Fig. 12 reach this condition, i.e., escape, within 8 min and many reach this condition within 5 min, except for olivine initially located at a longitude of 75°, which took 14 and 19 min to reach this condition for r_{a} = 1 and 2 μm, respectively.
Equation (12) gives a_{⊙} ~ 60 mm s^{2}, which is ~50a_{grav}, for fiducial particles with r_{a} = 1 μm. Thus, ignoring any thermal radiation and assuming a_{total} ≈ a_{⊙} for a firstorder approximation (as in Jewitt et al. 2013), the particle reaches the Hill sphere radius (r_{H} = 66 km; Jewitt et al. 2019) in 25 min or about onetenth of the rotational period of Phaethon. Only within about one day will the particle reach the distance of 2.5 × 10^{5} km, the dust tail length of Phaethon (Jewitt et al. 2013). Therefore, the activity should be visible very quickly (< 1 day) after the particles start to be generated.
The fact that most particles could not escape in Fig. 12a implies the importance of thermal radiation: If it were not for a_{ther}, particles could not escape easily. For example, let us take a closer look at trajectories of olivine particles in Fig. 12a. Many olivine particles could not escape from Phaethon. Among them, the olivine particles initially at longitude of 75° (15° before the morning terminator) could nevertheless escape and were flagged as A. What mechanism removed them from the AR? These particles have marginally positive net acceleration, (a_{ther} + a_{ref} > a_{grav} − a_{cen} in Fig. 5). Therefore, due to a_{ther}, it is possiblefor these particles to achieve nonzero outward speeds while rotating above the AR before sunrise. When the Sun rises and a_{⊙} comes into play, the particles gain high speed along the antisolar direction within a few minutes. Then, as the footprints of particles move towards midnight (longitude = 0 in Fig. 5), the contribution of a_{ther} increases again because of the higher temperature (Fig. 5). Meanwhile, v_{esc} (r) decreases as . Eventually, the particles are flagged as (A) at x ~ 320 m and 615 m for r_{a} = 1 and 2 μm, respectively. A similar explanation is applicable for magnetite in Fig. 12b. Magnetite particles in Fig. 12b are more strongly affected by a_{ther} than olivine at a longitude of < 90° (compare Figs. 5c and d), and hence escape from Phaethon more easily than olivine. For all particles initially at 90° ≲longitude ≲ 105°, the net acceleration is always inwards (toward the AR) in our Phaethon model (Figs. 5c and d), and therefore they cannot escape and are always flagged as (B).
If the thermal inertia Γ is updated as discussed in Sect. 4.2, one would expect particles to become more likely to escape from Phaethon if the initial position is near the terminator. The results are shown in Fig. 13, and as expected, particles are more easily ejected from the surface.
Fig. 11 Heliocentric distance (red plus; left ordinate) and the aspect angle (green dot; right ordinate) of asteroid (3200) Phaethon over time. (a) 50 rotations before and after the perihelion; (b) over one full orbit. Negative and positive times denote pre and postperihelion, respectively, where perihelion is indicated by the thick red dashed vertical line. The thin blue dashed vertical line in (a) indicates ± 10P_{rot} of Phaethon from its perihelion. The markers in (a) are separated by P_{rot} of Phaethon, while those in (b) are separated by 1 day. We note that the heliocentric distance changes by ~ 10% over ten rotations, and the aspect angle changes by ~1.5° per rotation near the perihelion. 
Fig. 12 Trajectory of particles projected onto Phaethon’s equatorial plane at the perihelion, initially at H = 1 cm and rotating with Phaethon’s rotational speed. Hollow circles, separated by 5° in longitude, show the initial position of the particles of the samecolored trajectory lines: r_{a} = 1 μm (solid) and 2 μm (dotted). Letters show the flag (see text) of 1 and 2 μm particles for the first and second line, respectively. The rotation is anticlockwise as indicated by the arrow, and the Sun is at x = − ∞ as shown by the hallow arrows, but below the xyplane since the aspect angle is 101° > 90°, so a_{total} has nonzeroz component in this figure. This is why the orange solid line (r_{a} = 1 μm at longitude 80°) seems to penetrate into the asteroid, as the trajectory is projected onto the xyplane. (a, c) olivine; (b, d) magnetite particles. (a, b) Morning terminator; (c, d) evening terminator. 
Fig. 13 Similar figures to Fig. 12, but the thermal inertia Γ is corrected based on Eq. (20). We note that many more cases are flagged as (A), i.e., escapeguaranteed. (a, c) Olivine; (b, d) magnetite particles. (a, b) morning terminator; (c, d) evening terminator. 
4.5 Phaethon’s diameter
It is important to mention that the diameter of Phaethon is not constrained well among previous publications. In this work, D = 5.1 km is adopted, because it is derived by the work that determined the thermophysical parameters and spin vector in Table 1 (Hanuš et al. 2016, 2018). A different thermal modeling study suggested with geometric albedo p_{V} = 0.16 ± 0.02 and thermal inertia, (Masiero et al. 2019), whereas a radar study suggested D > 6 km (Taylor et al. 2019).
Changes in diameter may affect thermal inertia and/or geometric albedo determination because the albedo is linked with the diameter and the absolute magnitude, and thermal inertia is a parameter used in TPM along with the albedo. It is beyond the scope of this work to discuss the detailed effects of the uncertainties on these parameters. However, the change in the bulk mass density ρ_{M} (Hanuš et al. 2018) does not affect our results. ρ_{M} is determined by Yarkovsky drift, that is, the change in the semimajor axis of the asteroid. The Yarkovsky acceleration of asteroid a_{Y} will be proportional to the reaction force from thermal radiation from the asteroid, F_{tr}, divided by its mass, M. Because F_{tr} ∝ D^{2} and M ∝ ρ_{M}D^{3}, the Yarkovsky acceleration a_{Y} ∝ 1∕ρ_{M}D. As a_{Y} is the observable, it remains constant, that is, ρ_{M}D = const. The gravitational acceleration a_{grav} ∝ ρ_{M}D (Eq. (18)), and therefore it must remain constant even if D is updated, because ρ_{M} is determined from the Yarkovsky drift. Hence, we argue that the update in D will only lead to a minor correction (e.g., initial speed of particles on the surface) to our results on Phaethon.
5 Conclusions and future work
In this work, the importance of thermal radiation from asteroidal regolith surface acting on particles of radii of a few microns is investigated for the first time. Our results show that particles of r_{a} ~ 1 μm are likely to have the largest outward net acceleration on asteroids in the nearSun environment (heliocentric distance ≲ 0.8 au). The full range of “ejectable” particle size spans from ~0.5 μm to ≲10 μm, although the exact range depends on the physical properties of the parent asteroid and particle, as well as initial conditions. This is consistent with observational studies of Phaethon’s activity at its perihelion that suggested the ejected particles have r_{a} ~ 1 μm.
However, this work should be taken with care because it is agnostic to the mechanisms that generate the particles and achieve the initial height above the regolith, similarly to other works like studies of electrostatic dust levitation. Therefore, future works may include investigations of those mechanisms, and could also include the likelihood of the occurrence of thermal fatigue as functions of heliocentric distance and location on the asteroid and distributions of the generated particles in the size, velocity, and asteroidcentric location domain. Other areas that remain to be investigated are the effects of simplifying assumptions used in thermal modeling. Most importantly, the problem of underestimation of the thermal radiation for heights greater than the isothermal length scale on the asteroid surface must be solved.
Acknowledgements
This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MEST) (No. NRF2018R1D1A1A09084105, 2019R1A6A1A10073437). Y.P.B. wants to thank HyeEun Jang for the discussions on the radiation pressure coefficient and Mie theory, Hojin Cho for providing discussions and priceless comments on many of the details of the manuscript, Harim Jin for discussing the physics of thermal fracturing, Sanghyuk Moon for discussing the numerical calculation of dust particle trajectory, Myung Gyoon Lee for commenting on observational results (all at the same institute as YPB), Scott Prahl at the Oregon Institute of Technology for sharing and revising the Mie theory code,Xu Wang at the University of Colorado, Boulder, and Xiao Ping Zhang at Macau University of Science and Technology for their kind discussions and explanations about the electrostatic dust lofting mechanism, Josef Hanuš at Charles University in Prague for sharing information on how the widely used thermal modeling code works, and Harald Mutschke at University of Jena for kind discussions on the optical index experiments of olivine and magnetite. We thank the reviewer, Oleksiy Golubov, for thoroughly reviewing the manuscript and greatly improved its quality in all aspects.
Appendix A Radiation pressure coefficients (Q_{pr}(λ, r_{a}))
As described in Sect. 2.3, the radiation pressure coefficient, Q_{pr} (λ, r_{a}), is calculated for two materials, olivine and magnetite. The calculated Q_{pr} values are shown as functions of wavelength λ for selected particle radii in Fig. A.2 (olivine) and Fig. A.1 (magnetite).
The optical index is available at λ > 0.1 μm for magnetite (Sect. 2.3) The derived Q_{pr} from Mie theory is shown in Fig. A.1.
For (lowFe) olivine, two datasets of optical indices are combined (Fabian et al. 2001; Sect. 2.3). The first dataset is composed of three indices as functions of wavelength λ = 2.04097 μm. The three indices are obtained from incident light polarized parallel to each of the three different crystalline axes of olivine (orthorhombic). These are denoted E∥x, E∥y, or E∥z following the notationin the original publication. The second dataset is obtained based on a measurement for a randomly oriented olivine sample in λ = 0.22.0 μm. After combining the two datasets, the optical index is identical for λ < 2.0 μm, but there are three separate cases for λ > 2.0 μm. These three different optical indices result in three separate Q_{pr} curves in Fig. A.2 and curves in Fig 2.
We further assumed Q_{pr} = 1 for 0.1 μm ≤ λ ≤ 0.2 μm for olivine to match the wavelength range with magnetite. For small particles of r_{a} ≪ 1 μm, the value of Q_{pr} (λ ∈ [0.1, 0.2] μm, r_{a}) may be different from unity (as in the r_{a} = 0.1 μm curve in Fig. A.2). However, this deviation is not critical because the contri bution of the radiation flux at the short wavelengths to the total acceleration is small (for blackbody radiation of temperature T ≲ T_{⊙}). Furthermore, as described in Sect. 2.4 and Fig. 2, the resulting difference in is negligible. Hence, for E∥x, E∥y, and E∥z are averaged to get olivine’s value used throughout this study.
Fig. A.1 Radiation pressure coefficient, Q_{pr}(λ, r_{a}) of magnetite (top) and the function b_{λ}(λ, T) (bottom; Eq. 9) for selected temperatures. 
Fig. A.2 Radiation pressure coefficient, Q_{pr}(λ, r_{a}) of olivine (first three rows) and the function b_{λ}(λ, T) (bottom; Eq. 9) for selected temperatures as in Fig. A.1. The three different Q_{pr} values are derived from three different optical indices depending on the polarization of the incident light with respectto the crystalline structure (see text). 
Appendix B Radiative accelerations derivations
In this section, we derive the equations we used in Sect. 2.4.
B.1 Auxiliary notations
The blackbody spectrumaveraged emissivity can be defined as (B.1)
We derived Eq. (10) using the Stefan–Boltzmann law, (B.2)
Next, the radiance (Planck function multiplied by the emissivity) can be normalized by the total flux: (B.3)
Subsequently, substituting Eq. (B.1) and (B.2), we derive Eq. (9). Similarly, Q_{pr} can be averaged over the radiance, i.e., (B.4)
We then useEq. (9) to derive Eq. (8). Using the above notations, a simplified notation without integral is developed: (B.5)
This relation is used to derive Eq. (11).
B.2 Derivation of Eq. (11)
When an object with a geometric crosssectional area of is illuminated by a spectral irradiance of J_{λ}(λ), the magnitude of the radiative force, i.e., the momentum flux per time, is (B.6)
The magnitude of the radiative acceleration that the particle experiences is a = ∥F∥∕m, where is the particle mass. If the source is a blackbody with a solid angle of dΩ and emissivity ϵ, the irradiance J_{λ} = ϵB_{λ}dΩ. Substituting m, J_{λ}, and Eq. (B.5), a can be rewritten as (B.7)
which is Eq.(11).
B.3 Derivation of Eq. (13)
For the a_{ref} calculation, the irradiance must be integrated over the surface elements within the field of view. As we assume the AR is a Lambertian scatterer, it leads to (B.8)
where is the bidirectional reflectance, and A_{B} is the Bond albedo of the Lambertian patch (for a derivation, see, e.g., Hapke 2012, Sect. 8.5.1). Therefore, the irradiance from the infinitesimal patch of solid angle dΩ is (B.9)
where μ_{i} = max {cosi, 0}. By using μ_{i}, the reflected radiation is “turned off” immediately after the sunset when calculated by a code.
The Lambertian scattering leads to an azimuthal symmetry of the scattered radiation from the local surface, so the acceleration of the particle perpendicular to the zaxis will cancel out. The acceleration along the zaxis from the infinitesimal patch is therefore (see Fig. 1) (B.10)
The infinitesimal annulus at angle e″ and width de″ has the solid angle of , and any patch within this annulus will give an identical acceleration contribution to the particle along the zdirection. Thus, the contribution of the annulus is obtained by replacing dΩ with in Eq (B.10). The integration of the annulus of radius from 0 to r_{0} is identical to the integration of e″ from 0 to . If the height factor is defined as Eq. (14), (B.11)
Thus, after performing the integration, Eq. (13) is obtained.
B.4 Derivation of Eq. (15)
For a_{ther}, similar to a_{ref}, the azimuthal symmetry is guaranteed. Meanwhile, the acceleration (da_{ther}) due to the radiation originating from an infinitesimal area on the AR with solid angle is given by Eq. (B.7), with , T = T_{S}, and . Then, multiplying cose″ and integrating over , in the same way as done for a_{ref}, leads to Eq. (15).
Appendix C Thermal fatigue reparameterization
Here, we derive the two criteria we used in Sect. 4.3.1. The two criteria are from Ravaji et al. (2019), but we devise a more general view using the (r_{h}, Θ) coordinate by reformulating them.
For the first one, , we rewrite
for a dimensionless temperature used in TPM, u = T∕T_{eqm} (and ). The u profile, and thus the Δû_{S}, depend only on the thermal parameter Θ (and the spin vector; see Sect. 2.2). The term on the righthand side, , is the ΔT_{S} when the asteroid is at r_{h} but with c_{T} = 1 (e.g., A_{B} = 0 and ; see Eq. 4). All values other than c_{T} are calculable once the coordinate (r_{h}, Θ) is given, and so it is mathematically convenient to use (C.1)
rather than to represent general asteroids.
Similarly, the second criterion is reformulated using
where is the time bin and nlon is the number of longitude bin (nlon is fixed in the TPM simulation). Δu(Θ) on the righthand side is calculated numerically by taking the difference of the consecutive u array elements in the code (u[1:] – u[: −1]). Therefore, is obtained for the coordinate grid (r_{h}, Θ). is a constant for the conducted TPM. The second criterion is obtained by taking the maximum of the absolute values of the Δu array in the code, and multiplying it by the constants. Thus, it is more convenient to use (C.2)
rather than to represent general asteroids.
To generate the plot in Fig. 10, we set the grid of Θ on a logarithmic scale: Θ = 10^{x} for 50 uniformlyspaced x ∈ [−2, +2]. The aspect angle is fixed to 90°. First the equatorial temperature profile for these 50 models, u(Θ) = T∕T_{eqm}, are calculated (see Sect. 2.2). Then the temperature is restored by multiplying T_{eqm} at each r_{h} with c_{T} = 1. The contours are calculated for the reformulated formulae, i.e., the contours are and for a fictitious object of c_{T} = 1 and P_{rot} = 1 h.
For c_{T} ≪ 1 (A_{B} → 1), the first condition pushes the region for thermal fatigue to the lower left corner in Fig. 10. This is expected because, if most sunlight is reflected, we need small thermal conduction (Γ → 0) and higher insolation (r_{h} → 0). Meanwhile, for a fast rotator (P_{rot} → 0), the second criterion requires a similar condition. This is also expected, because the temperature profile will be nearly constant at a given latitude for such an asteroid, and the temperature change rate will be very small unless it is close to the Sun and there is small thermal conduction.
Appendix D Code and data availability
All the scripts, source codes and developed packages used to generate datasets and plots in this work are available via the GitHub service^{5}. As the model calculation may take a long time, the tables of Q_{pr} and , as well as the refractive indices of olivine and magnetite, are available from the GitHub repository. Also, the dataset generated by one of the notebooks in the GitHub repository (04accplot.ipynb), which is used for generating Fig. 6, 7, and 8, is accessible via a separate figshare service^{6}, due to its large size. The package miepython is also available via the GitHub service^{7}.
References
 Arai, T., Kobayashi, M., Ishibashi, K., et al. 2018, in Lunar and Planetary Science Conference, 2570 [Google Scholar]
 Blaauw, R. C. 2017, Planet. Space Sci., 143, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Bohren, C. F., & Huffman, D. R. 1983, Absorption and Scattering of Light by Small Particles (New York: Wiley) [Google Scholar]
 Bohren, C. F., & Huffman, D. R. 1998, Absorption and Scattering of Light by Small Particles (New York: Wiley) [CrossRef] [Google Scholar]
 Bowell, E., Hapke, B., Domingue, D., et al. 1989, in Asteroids II, eds. R. P. Binzel, T. Gehrels, & M. S. Matthews, 524 [Google Scholar]
 Carry, B. 2012, Planet. Space Sci., 73, 98 [CrossRef] [Google Scholar]
 Chamberlin, A. B., McFadden, L.A., Schulz, R., Schleicher, D. G., & Bus, S. J. 1996, Icarus, 119, 173 [NASA ADS] [CrossRef] [Google Scholar]
 Delbo, M. 2004, PhD thesis, Free University of Berlin [Google Scholar]
 Delbo’, M., dell’Oro, A., Harris, A. W., Mottola, S., & Mueller, M. 2007, Icarus, 190, 236 [CrossRef] [Google Scholar]
 Delbo, M., Libourel, G., Wilkerson, J., et al. 2014, Nature, 508, 233 [Google Scholar]
 Delbo, M., Mueller, M., Emery, J. P., Rozitis, B., & Capria, M. T. 2015, Asteroid Thermophysical Modeling, 107 [Google Scholar]
 Fabian, D., Henning, T., Jäger, C., et al. 2001, A&A, 378, 228 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Green, S. F., Meadows, A. J., & Davies, J. K. 1985, MNRAS, 214, 29P [NASA ADS] [CrossRef] [Google Scholar]
 Grott, M., Knollenberg, J., Hamm, M., et al. 2019, Nat. Astron., 3, 971 [Google Scholar]
 Hanuš, J., Delbo’, M., Vokrouhlický, D., et al. 2016, A&A, 592, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hanuš, J., Vokrouhlický, D., Delbo’, M., et al. 2018, A&A, 620, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hapke, B. 2012, Theory of Reflectance and Emittance Spectroscopy (Cambridge University Press) [Google Scholar]
 Harris, A. W. 1998, Icarus, 131, 291 [Google Scholar]
 Hartzell, C. M., & Scheeres, D. J. 2011, Planet. Space Sci., 59, 1758 [NASA ADS] [CrossRef] [Google Scholar]
 Hesar, S. G., Scheeres, D. J., McMahon, J. W., & Rozitis, B. 2017, J. Guidance Control Dyn., 40, 2432 [NASA ADS] [CrossRef] [Google Scholar]
 Hood, N., Carroll, A., Mike, R., et al. 2018, Geophys. Res. Lett., 45, 13, 206 [Google Scholar]
 Hsieh, H. H., & Jewitt, D. 2005, ApJ, 624, 1093 [NASA ADS] [CrossRef] [Google Scholar]
 Hui, M.T., & Li, J. 2017, AJ, 153, 23 [NASA ADS] [Google Scholar]
 Jewitt, D. 2012, AJ, 143, 66 [CrossRef] [Google Scholar]
 Jewitt, D., & Li, J. 2010, AJ, 140, 1519 [NASA ADS] [CrossRef] [Google Scholar]
 Jewitt, D., Li, J., & Agarwal, J. 2013, ApJ, 771, L36 [NASA ADS] [CrossRef] [Google Scholar]
 Jewitt, D., Mutchler, M., Agarwal, J., & Li, J. 2018, AJ, 156, 238 [NASA ADS] [CrossRef] [Google Scholar]
 Jewitt, D., Asmus, D., Yang, B., & Li, J. 2019, AJ, 157, 193 [NASA ADS] [CrossRef] [Google Scholar]
 Lauretta, D. S., Hergenrother, C. W., Chesley, S. R., et al. 2019, Science, 366, 3544 [NASA ADS] [CrossRef] [Google Scholar]
 Li, J., & Jewitt, D. 2013, AJ, 145, 154 [NASA ADS] [CrossRef] [Google Scholar]
 Luu, J. X., & Jewitt, D. C. 1992, Icarus, 97, 276 [NASA ADS] [CrossRef] [Google Scholar]
 Masiero, J. R., Wright, E. L., & Mainzer, A. K. 2019, AJ, 158, 97 [NASA ADS] [CrossRef] [Google Scholar]
 McMahon, J. W., Scheeres, D. J., Chesley, S. R., et al. 2020, J. Geophys. Res.: Planets [Google Scholar]
 Molaro, J. L., Byrne, S., & Le, J. L. 2017, Icarus, 294, 247 [CrossRef] [Google Scholar]
 Mueller, M. 2007, PhD thesis, Free University of Berlin [Google Scholar]
 Myhrvold, N. 2016, PASP, 128, 045004 [NASA ADS] [CrossRef] [Google Scholar]
 Orger, N. C., Cordova Alarcon, J. R., Toyoda, K., & Cho, M. 2018, Adv. Space Res., 62, 896 [NASA ADS] [CrossRef] [Google Scholar]
 Putzig, N. E. 2006, PhD thesis, University of Colorado at Boulder [Google Scholar]
 Ravaji, B., AlíLagoa, V., Delbo, M., & Wilkerson, J. W. 2019, J. Geophys. Res. (Planets), 124, 3304 [NASA ADS] [CrossRef] [Google Scholar]
 Righter, K., Drake, M. J., & Scott, E. R. D. 2006, Compositional Relationships Between Meteorites and Terrestrial Planets, eds. D. S. Lauretta, & H. Y. McSween, 803 [Google Scholar]
 Robie, R. A., & Bethke, P. M. 1962, Molar Volumes and Densities of Minerals (Unites States (Unites States Department of the Interior Geological Survey) [Google Scholar]
 Ryabova, G. O. 2017, Planet. Space Sci., 143, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Scott, E. R. D., & Krot, A. N. 2005, ApJ, 623, 571 [NASA ADS] [CrossRef] [Google Scholar]
 Senshu, H., Kimura, H., Yamamoto, T., et al. 2015, Planet. Space Sci., 116, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Spencer, J. R., Lebofsky, L. A., & Sykes, M. V. 1989, Icarus, 78, 337 [Google Scholar]
 Taylor, P. A., RiveraValentín, E. G., Benner, L. A. M., et al. 2019, Planet. Space Sci., 167, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Tsuchiyama, A., Uesugi, M., Matsushima, T., et al. 2011, Science, 333, 1125 [NASA ADS] [CrossRef] [Google Scholar]
 van de Hulst, H. C. 1981, Light Scattering by Small Particles [Google Scholar]
 Wang, X., Schwan, J., Hsu, H. W., Grün, E., & Horányi, M. 2016, Geophys. Res. Lett., 43, 6103 [NASA ADS] [CrossRef] [Google Scholar]
 Warner, B. D., Harris, A. W., & Pravec, P. 2009, Icarus, 202, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Wiegert, P. A., Houde, M., & Peng, R. 2008, Icarus, 194, 843 [NASA ADS] [CrossRef] [Google Scholar]
 Zimmerman, M. I., Farrell, W. M., Hartzell, C. M., et al. 2016, J. Geophys. Res. (Planets), 121, 2150 [NASA ADS] [CrossRef] [Google Scholar]
The compact material density (the density when microporosity is zero) is ~3200 kg m^{3} for forsterite (nonFe olivine; Mg_{2}SiO_{4}), ~4400 kg m^{3} for fayalite (nonMg olivine; Fe_{2}SiO_{4}), and ~5200 kg m^{3} for iron (II, III) oxide (pure Fe_{3}O_{4}; all densities from Robie & Bethke 1962). The adopted density ρ_{m} = 3000 kg m^{3} can therefore be understood as a homogeneously porous sphere with certain microporosity in the Mie theory.
https://github.com/scottprahl/miepython; the specific version we used is https://github.com/scottprahl/miepython/tree/0129cedb231c5e60170860b69cdb06eb26be5669
All Tables
Physical parameters of Phaethon (Hanuš et al. 2018) and the fiducial and tested parameters for fictitious model asteroids are given.
All Figures
Fig. 1 Schematic diagram of a spherical particle hovering over the surface. The particle of radius r_{a} hovers over the planar regolith at the height of H(≫ r_{a}). The solar irradiation (the long red arrow) is originated from the direction . A tiny surface patch with an area of dA reflects the solar radiation and emits the thermal radiation to the particle (the dashed and squiggly arrows, respectively, to the directionof ) with the emission angle of e′. In our model, e″ = e′ holds because of the planar surface assumption. The isothermal regolith region (gray dashed) is assumed to be a circle with a radius of r_{0}(≫ H ≫ r_{a}). No radiation outside this region is considered because the patch right below the particle already fills most of the field of view (r_{0} ≫ H). 

In the text 
Fig. 2 Calculated Q_{pr} value averaged over blackbody spectra of different temperatures ( in our notation). The values are plotted as a function of the particle size r_{a} and blackbody temperature T of the incident radiation for two minerals: a) olivine and b) magnetite. The emissivity ϵ is assumed to be constant, so cancels out in Eq. (9). For olivine, three polarization directions with respect to the crystalline axes can lead to small differences in the results (three different markers; see Appendix A), and the average value (solid line) is used in this work. converges tounity for a large particle limit, as predicted from the geometrical optics approximation. 

In the text 
Fig. 3 Equatorial temperature profiles. a, the fictitious asteroid at r_{h} = 0.20 au and b, Phaethon at perihelion (r_{h} = 0.14 au). For the model asteroids, different aspect angles are considered: perpendicular (solid) and 45° (dashed). The different colors of the lines indicate different thermal inertia values: black lines for Γ = 0, light blue lines for fiducial value for model asteroids and Phaethon (Table 1), and red lines for the values considering the temperaturedependent thermal inertia (Γ = 669 tiu and 2622 tiu for model asteroids and Phaethon, respectively; see Sect. 4.2). 

In the text 
Fig. 4 Components of the radiative accelerations acting on a r_{a} = 2 μm particle at the height of H = 1 cm on the equator of D = 1 km model asteroid as a function of longitude (local time), for a) olivine; b) magnetite. The black horizontal lines show the gravity (dashed) and the effective gravity taking the centrifugal acceleration into account (solid). The absolute value of is used as it is negative. Both a_{⊙} and a_{ref} turn off almost immediately after Sunset, which causes the spiky feature of at longitude of 270°. a_{ther} persists throughout the night thanks to nonzero thermal inertia. 

In the text 
Fig. 5 Radial accelerations on the equatorial surface at heights of H = 1 cm, for the corresponding particle radius r_{a} with respect to longitude. The net radial radiative accelerations () are shown.The black horizontal lines show the gravity and the effective gravity as in Fig. 4. a, c) olivine; b, d) magnetite. Panels a and b are calculated for fiducial perpendicular (solid) and aspect 45° (dashed) model asteroids (D = 1 km) at the heliocentric distance r_{h} = 0.2 au for given particle radii r_{a} in the legend. For the perpendicular model, the 25% ambiguity inthe thermal inertia (Γ = 200 ± 50 tiu) is expressed as shading. Panels c and d are calculated for Phaethon using the nominal physical parameter values (Table 1) at the perihelion. The shading signifies the uncertainty range in the thermal inertia (Γ = 600 ± 200 tiu, Hanuš et al. 2018) of Phaethon. From these figures, it is obvious that most of ≲ 10μmsized particles areaccelerated outwards from the AR. 

In the text 
Fig. 6 Maximum total radial acceleration, [mm s^{2}] of particleson the equator of the fictitious model asteroids. Each row corresponds to a different asteroid model a) to d) perpendicular model with P_{rot} = 6 h, e) to h) perpendicular model with P_{rot} = 3 h, i) to l) aspect 45° model with P_{rot} = 6 h, and m) to p) aspect 45° model with P_{rot} = 3 h). Each column corresponds to different heliocentric distances. The maximum acceleration is calculated by finding the maximum value along the longitude on the equator for the given asteroid diameter D and particle radius r_{a}. The blue dashed lines and red solid lines indicate olivine and magnetite, respectively. The thickness of the contour increases as the acceleration value increases for fixed contour levels . Negative (inward) accelerations are not drawn. 

In the text 
Fig. 7 Maximum diameters of asteroids (D) that producepositive on the equator as a function of heliocentric distance (r_{h} ∈ [0.1, 2.0] au). The red solid and blue dashed lines are the results for magnetite and olivine, respectively. The results of perpendicular and aspect 45° models are discriminated by the line thicknesses, and the results of different rotational periods (6 h and 3 h) are discriminated by the transparency of the lines. The black lines show the slope according to a power law. 

In the text 
Fig. 8 Maximum radial acceleration on the Phaethon’s equator. As in Fig. 6, the maximum acceleration on the equator is calculated for a wide range of particle sizes (0.5 40 μm) before (a) and after (b) the perihelion passage. The black vertical line corresponds to thePhaethon’s perihelion. The numbers on the contour plots are in units of mm s^{2}, and the thickness of the contour increases as the acceleration value increases for fixed contour levels . The blue dashed and red solid lines indicate olivine and magnetite, respectively. 

In the text 
Fig. 9 Identical plot to Fig. 5, except the thermal inertia of asteroids is changed according to Eq. (20). (a, c) Olivine; (b, d) magnetite. The shade in each plot signifies the uncertainty on thermal inertia, scaled by the same factor, . 

In the text 
Fig. 10 Thermal fatigue plot. The thermal parameter (Θ) over heliocentric distance (r_{h}) for Phaethon (solid lines, r_{h} ≥ 0.14 au) and fiducial perpendicular model asteroids with varying r_{h} (dashedlines) are shown. The black star and plus symbols indicate that Γ is fixed, while the green triangle and cross markers show results where Γ is assumed to be changed by Eq. (20). The blue solid lines are the contours of , and the brown dashed lines are that of . The values shown for the contours are the and values for the case when c_{T} = 1 and P_{rot} = 1 h. The contours get thicker for larger values to guide the eyes. 

In the text 
Fig. 11 Heliocentric distance (red plus; left ordinate) and the aspect angle (green dot; right ordinate) of asteroid (3200) Phaethon over time. (a) 50 rotations before and after the perihelion; (b) over one full orbit. Negative and positive times denote pre and postperihelion, respectively, where perihelion is indicated by the thick red dashed vertical line. The thin blue dashed vertical line in (a) indicates ± 10P_{rot} of Phaethon from its perihelion. The markers in (a) are separated by P_{rot} of Phaethon, while those in (b) are separated by 1 day. We note that the heliocentric distance changes by ~ 10% over ten rotations, and the aspect angle changes by ~1.5° per rotation near the perihelion. 

In the text 
Fig. 12 Trajectory of particles projected onto Phaethon’s equatorial plane at the perihelion, initially at H = 1 cm and rotating with Phaethon’s rotational speed. Hollow circles, separated by 5° in longitude, show the initial position of the particles of the samecolored trajectory lines: r_{a} = 1 μm (solid) and 2 μm (dotted). Letters show the flag (see text) of 1 and 2 μm particles for the first and second line, respectively. The rotation is anticlockwise as indicated by the arrow, and the Sun is at x = − ∞ as shown by the hallow arrows, but below the xyplane since the aspect angle is 101° > 90°, so a_{total} has nonzeroz component in this figure. This is why the orange solid line (r_{a} = 1 μm at longitude 80°) seems to penetrate into the asteroid, as the trajectory is projected onto the xyplane. (a, c) olivine; (b, d) magnetite particles. (a, b) Morning terminator; (c, d) evening terminator. 

In the text 
Fig. 13 Similar figures to Fig. 12, but the thermal inertia Γ is corrected based on Eq. (20). We note that many more cases are flagged as (A), i.e., escapeguaranteed. (a, c) Olivine; (b, d) magnetite particles. (a, b) morning terminator; (c, d) evening terminator. 

In the text 
Fig. A.1 Radiation pressure coefficient, Q_{pr}(λ, r_{a}) of magnetite (top) and the function b_{λ}(λ, T) (bottom; Eq. 9) for selected temperatures. 

In the text 
Fig. A.2 Radiation pressure coefficient, Q_{pr}(λ, r_{a}) of olivine (first three rows) and the function b_{λ}(λ, T) (bottom; Eq. 9) for selected temperatures as in Fig. A.1. The three different Q_{pr} values are derived from three different optical indices depending on the polarization of the incident light with respectto the crystalline structure (see text). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.