© 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 near-Sun 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 near-Sun 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 near-Sun 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 non-inertial 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 so-called thermo-physical 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 time-centered space scheme following the above-mentioned 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 large-scale 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 (self-heating, 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_eqme^−z∕l_S 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⁻⁵ 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 non-rotating 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 spectrum-averaged 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 so-called 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 well-studied 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: low-Fe olivine (Mg_1.9Fe_0.1SiO₄; Fabian et al. 2001¹) and magnetite (Fe₃O₄; Triaud, A. H. M. J.²; 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, spectrum-averaged 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 spectrum-averaged Q_pr value, b_λ in Eq. (9) is the normalized Planck function, and is the spectrum-averaged 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 ~ 100-1000 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 (face-centered 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 ra hovers over the planar regolith at the height of H(≫ ra). 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 r0(≫ H ≫ ra). No radiation outside this region is considered because the patch right below the particle already fills most of the field of view (r0 ≫ H).

Fig. 2 Calculated Qpr value averaged over blackbody spectra of different temperatures ( in our notation). The values are plotted as a function of the particle size ra 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 anti-solar 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₀ is set as 1% of the asteroid’s radius throughout this work. Here, but is almost unity at the initial condition, as r₀ ≫ 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₀. 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 z-axis components of a_⊙ in Eq. (12) and a_ther in Eq. (15), we re-formulated 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 cross-sectional areas of the particle exposed to the spectrum of temperature T_S and T_⊙. The ratio of the two acceleration mechanisms along the z-axis 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 near-Earth 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 near-Sun 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(log₁₀20)∕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(log₁₀800)∕49−1} km for integer k = 0 to 49). For S, B, and C-type 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.1-0.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³); 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 above-mentioned 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).

Fig. 3 Equatorial temperature profiles. a, the fictitious asteroid at rh = 0.20 au and b, Phaethon at perihelion (rh = 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 temperature-dependent 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⁴ 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 ra = 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 aref turn off almost immediately after Sunset, which causes the spiky feature of at longitude of 270°. ather 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 Phaethon-like 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 near-Sun 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 ≈ 2-3 μ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 near-Sun 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 ra 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 rh = 0.2 au for given particle radii ra 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μm-sized 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 Prot = 6 h, e) to h) perpendicular model with Prot = 3 h, i) to l) aspect 45° model with Prot = 6 h, and m) to p) aspect 45° model with Prot = 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 ra. 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 (rh ∈ [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 near-perihelion 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 ~ 2-3 μ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 post-perihelion) 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 time-symmetric 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 area-to-mass-ratio 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 first-order approximation to cope with this temperature dependence is to roughly scale Γ by or (20)

where Γ₁ 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 Γ₁ values of the similar-sized 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.1-1 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 (rh) for Phaethon (solid lines, rh ≥ 0.14 au) and fiducial perpendicular model asteroids with varying rh (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 cT = 1 and Prot = 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 time-rate 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 re-formulated by (22)

and (23)

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 right-hand 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 right-hand 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 = 5-22 μm (2r_a = 10-44 μ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 time-symmetry in the illuminationcondition on Phaethon along the latitude: The thermal (illumination) history of the north during the perihelion passage is nearly the time-reversal 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 two-fold. 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₀ = π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 kick-drift-kick) 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 asteroid-centric 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 first-order 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 one-tenth of the rotational period of Phaethon. Only within about one day will the particle reach the distance of 2.5 × 10⁵ 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 anti-solar 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 post-perihelion, respectively, where perihelion is indicated by the thick red dashed vertical line. The thin blue dashed vertical line in (a) indicates ± 10Prot of Phaethon from its perihelion. The markers in (a) are separated by Prot 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 same-colored trajectory lines: ra = 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 anti-clockwise as indicated by the arrow, and the Sun is at x = − ∞ as shown by the hallow arrows, but below the xy-plane since the aspect angle is 101° > 90°, so atotal has non-zeroz component in this figure. This is why the orange solid line (ra = 1 μm at longitude 80°) seems to penetrate into the asteroid, as the trajectory is projected onto the xy-plane. (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., escape-guaranteed. (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² and M ∝ ρ_MD³, the Yarkovsky acceleration a_Y ∝ 1∕ρ_MD. As a_Y is the observable, it remains constant, that is, ρ_MD = const. The gravitational acceleration a_grav ∝ ρ_MD (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 near-Sun 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 asteroid-centric 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. NRF-2018R1D1A1A09084105, 2019R1A6A1A10073437). Y.P.B. wants to thank Hye-Eun 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 (low-Fe) 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.0-4097 μ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.2-2.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, Qpr(λ, ra) of magnetite (top) and the function bλ(λ, T) (bottom; Eq. 9) for selected temperatures.

Fig. A.2 Radiation pressure coefficient, Qpr(λ, ra) of olivine (first three rows) and the function bλ(λ, T) (bottom; Eq. 9) for selected temperatures as in Fig. A.1. The three different Qpr 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 spectrum-averaged 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 cross-sectional 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 z-axis will cancel out. The acceleration along the z-axis 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 z-direction. 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₀ 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 re-formulating 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 right-hand 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 right-hand 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 uniformly-spaced 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 re-formulated 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⁵. 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 (04-accplot.ipynb), which is used for generating Fig. 6, 7, and 8, is accessible via a separate figshare service⁶, due to its large size. The package miepython is also available via the GitHub service⁷.

References

All Tables

All Figures

Fig. 1 Schematic diagram of a spherical particle hovering over the surface. The particle of radius ra hovers over the planar regolith at the height of H(≫ ra). 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 r0(≫ H ≫ ra). No radiation outside this region is considered because the patch right below the particle already fills most of the field of view (r0 ≫ H).

In the text

Fig. 2 Calculated Qpr value averaged over blackbody spectra of different temperatures ( in our notation). The values are plotted as a function of the particle size ra 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 rh = 0.20 au and b, Phaethon at perihelion (rh = 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 temperature-dependent 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 ra = 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 aref turn off almost immediately after Sunset, which causes the spiky feature of at longitude of 270°. ather 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 ra 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 rh = 0.2 au for given particle radii ra 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μm-sized 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 Prot = 6 h, e) to h) perpendicular model with Prot = 3 h, i) to l) aspect 45° model with Prot = 6 h, and m) to p) aspect 45° model with Prot = 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 ra. 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 (rh ∈ [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 (rh) for Phaethon (solid lines, rh ≥ 0.14 au) and fiducial perpendicular model asteroids with varying rh (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 cT = 1 and Prot = 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 post-perihelion, respectively, where perihelion is indicated by the thick red dashed vertical line. The thin blue dashed vertical line in (a) indicates ± 10Prot of Phaethon from its perihelion. The markers in (a) are separated by Prot 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 same-colored trajectory lines: ra = 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 anti-clockwise as indicated by the arrow, and the Sun is at x = − ∞ as shown by the hallow arrows, but below the xy-plane since the aspect angle is 101° > 90°, so atotal has non-zeroz component in this figure. This is why the orange solid line (ra = 1 μm at longitude 80°) seems to penetrate into the asteroid, as the trajectory is projected onto the xy-plane. (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., escape-guaranteed. (a, c) Olivine; (b, d) magnetite particles. (a, b) morning terminator; (c, d) evening terminator.
In the text

	Fig. A.1 Radiation pressure coefficient, Qpr(λ, ra) of magnetite (top) and the function bλ(λ, T) (bottom; Eq. 9) for selected temperatures.
In the text

	Fig. A.2 Radiation pressure coefficient, Qpr(λ, ra) of olivine (first three rows) and the function bλ(λ, T) (bottom; Eq. 9) for selected temperatures as in Fig. A.1. The three different Qpr 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