© The Authors 2025

1 Introduction

The physical behavior of asteroids near the Sun is mysterious. Granvik et al. (2016) compared a detected asteroid population with the de-biased near-Earth object population model and suggested a lack of near-Sun asteroid detections, implying that objects with low perihelion distances may be prone to catastrophic disruptions. A near-Sun B-type asteroid, (3200) Phaethon, may help answer this question because its observed mass loss has been reported. It has an orbit with a large eccentricity (e = 0.89) and a short perihelion distance (q = 0.14 au). Although the exact mechanism driving the near-Sun “supercatastrophic” disruption is unknown, thermal processes are suspected (Granvik et al. 2016). In this context, studying Phaethon may provide insight into the effects of the extreme variation in the diurnal and seasonal temperatures on asteroids (Bach & Ishiguro 2021; Holt et al. 2022). Phaethon is also notable for being the parent body of the Geminid meteoroid stream, which is responsible for one of the most prominent annual meteor showers. The proposed mechanisms for the formation of Geminids include cometary activity (Gustafson 1989), thermal fracture (Jewitt & Li 2010; Jewitt 2012), and rotational instability (Nakano & Hirabayashi 2020; Jo & Ishiguro 2024). On the other hand, Phaethon has exhibited near-perihelion activity in multiple perihelion passages. These activities were initially interpreted as dust ejection (Jewitt & Li 2010; Li & Jewitt 2013; Jewitt et al. 2013; Hui & Li 2017). However, recent studies favor the idea that the observed activities may instead be due to the outgassing of sodium or iron (Lisse & Steckloff 2022; Hui 2023; Zhang et al. 2023). In any case, the connection between the present-day activity of Phaethon and the Geminid formation event remains uncertain. Phaethon has been chosen as the flyby target for the JAXA Demonstration and Experiment of Space Technology for INterplanetary voYage with Phaethon fLyby and dUst Science (DESTINY⁺) mission (Arai et al. 2018). This space mission may therefore provide clues regarding the various mysterious aspects of asteroids near the Sun.

Phaethon has a diameter of ∼5 km and a rotation period of 3.603957 hours (Kim et al. 2018; Hanuš et al. 2018; MacLennan et al. 2022). Although this rotation period is above the critical spin barrier for asteroids (Jewitt 2012), the centrifugal force is comparable to its gravity, creating a complex surface dynamical environment. For this reason, fast-rotating asteroids selected as targets for space missions have been subjected to preflight geophysical analysis. For example, Scheeres et al. (2016) provided the pre-encounter surface dynamics analysis of Origins, Spectral Interpretation, Resource Identification, Security, Regolith Explorer (OSIRIS-REx) target (101955) Bennu, which was later followed up by the post-encounter analysis of Scheeres et al. (2019). Similarly, Yu et al. (2019) studied the surface stability of (65803) Didymos prior to the Double Asteroid Redirection Test (DART) mission. The surface dynamical state of (16) Psyche was investigated by Moura et al. (2020) as part of the preparation for the recently launched Psyche mission. It is crucial to perform a geophysical study of Phaethon to anticipate observational features that DESTINY⁺ will detect during the flyby phase and use this geological information to study the link between Phaethon’s current and past activities.

We investigated the geophysical environment of Phaethon using the latest shape model as part of the preflight preparation efforts for the DESTINY⁺ mission. Since DESTINY⁺ is a flyby mission that can only observe about half of Phaethon’s surface, preflight planning and the selection of the closest encounter timing are crucial. Section 2 outlines the methods and parameters we used to calculate Phaethon’s key physical quantity: its geopotential. In Sect. 3 we present the results of these calculations and various derivations related to the dynamical environment on and around Phaethon. The uncertainties and implications of our findings are discussed in Sect. 4. In Sect. 5 we summarize this work.

2 Method

To calculate and interpret the physical quantities related to the surface dynamics of Phaethon, we had to first establish the physical characteristics of the asteroid. In this section we describe our method for deriving the geopotential of Phaethon.

The geopotential is the sum of the gravitational potential and rotational potential, forming the effective potential in the corotating frame of the body. The geopotential equation comprises two parts: the centrifugal and gravitational terms. For the centrifugal term, following observational reports, we assumed that Phaethon is a principal axis rotator with a period of 3.603957 hours (Kim et al. 2018; Hanuš et al. 2018). In this assumption, the reference frame has its origin at the mass center of Phaethon, with the z-axis aligned with the rotation axis. With this setup, the geopotential at a point r is given as follows: $$V(\mathit{r})=-\frac{1}{2}(\mathit{\omega }\times \mathit{r})\cdot (\omega \times \mathit{r})+U(\mathit{r}),$$(1)

where U(r) is the gravitational potential at r and ω = ωẑ is the spin vector of Phaethon with angular velocity ω. Although the centrifugal term can be uniquely computed once the rotational period is given, the calculation of the gravitational potential is more complicated. This is because, like most small bodies, Phaethon has an irregular shape. The analytical formula for the gravitational potential of an irregular shaped polyhedron was established with the key assumption that the density of the polyhedron is homogeneous (Werner & Scheeres 1996). Another method is the mass concentration (mascon) method (Geissler et al. 1996), where the body is filled with point masses of equal mass along a cubic grid space, whose combined gravity represents the gravity of the body. Although this configuration was criticized for its inaccuracy (Werner & Scheeres 1996), the spatial distribution of the point masses can be adjusted to reduce this limitation and emulate potential heterogeneity in asteroids (Yu 2016; Pearl & Hitt 2022).

In this work, we used the soft-sphere discrete element method (SSDEM) in the parallel N-body code PKDGRAV (Richardson et al. 2000; Stadel 2001; Schwartz et al. 2012) to fill Phaethon’s shape model with spherical particles. We divided the asteroid into two regions: the core and the envelope. The core consists of a 1.5 km radius spherical area where the particles, each 80 m in radius, are arranged in a hexagonal close packing (HCP) configuration. After generating the core, we surrounded it with a cloud of randomly distributed particles. By allowing these particles to gravitate naturally toward the core, we covered the core with a random close packing (RCP) envelope. The radii of the envelope particles are within the range of 10–80 m and follow the boulder size-frequency distribution of the asteroid (162173) Ryugu, as investigated by Hayabusa2 (Michikami et al. 2019). This two-layer structure represents a rubble-pile asteroid with a dense core, chosen to account for the possibility that the Geminids formed from Yarkovsky–O’Keefe–Radzievskii–Paddack (YORP) effect-induced rotational instability (Jo & Ishiguro 2024). Previous numerical studies have shown that surface mass shedding due to rotational instability is more likely on bodies with a dense core or particles of a wide size range (Hirabayashi et al. 2015; Zhang et al. 2017; Zhang & Michel 2021).

The material parameters used in this study follow conventional values used in previous studies (e.g., Zhang & Michel 2020). We summarize these values, along with their references, in Table 2. For the static friction coefficient, μ_n, and the shape parameter, β, we selected values that align with the approximate range of friction angles observed in terrestrial materials (Hirabayashi et al. 2015; Zhang et al. 2017). The normal spring constant k_n and the simulation timestep were chosen to ensure that the maximum velocity, equivalent to the freefall speed, would limit inter-particle overlaps to no more than 1% of the minimum particle radius during the simulation. We note that this k_n value is roughly equivalent to the elastic Young’s modulus of ∼1 MPa, which is approximately analogous to porous terrestrial gravels (DeMartini et al. 2019). Studies of Ryugu by the Hayabusa2 mission and (101955) Bennu by the OSIRIS-REx mission suggest that the regoliths of both asteroids have nearly zero cohesion (Arakawa et al. 2020; Walsh et al. 2022). Roberts et al. (2021) and Zhang et al. (2022) proposed that the interiors of these asteroids also have low cohesion. In accordance with these findings, we assumed that the particles in our model are cohesionless.

Once particle aggregation was complete, the modeled asteroid was rotated at Phaethon’s rotation period until the particles stabilized in the rotating environment. Next, we carved the aggregate into the shape model of Phaethon. The Phaethon shape model used in this work is based on 12 days of radar observations (from two apparitions), 195 light curves (from 18 oppositions), and eight occultations (from three apparitions). Its volumetric mean diameter is 5.0 km, with the x, y, and z diameters of the dynamically equivalent, equal volume ellipsoid measuring 5.88, 5.43, and 4.01 km, respectively (Marshall et al., in prep.).

Finally, the carved model was rotated again at Phaethon’s rotation period to complete the model preparation process. A cross-sectional view of the modeled asteroid is shown in Fig. 1. The envelope has a packing efficiency of 60.2%. For the core, the theoretical HCP packing efficiency is $\frac{\pi }{3\sqrt{2}}\approx 74\mathrm{\%}$, leading to a total packing efficiency of about 63.1%. This corresponds to a bulk porosity of 36.9% for the model. We adjusted the mass density of each particle to match Phaethon’s bulk density of 1.58 g/cm³ (MacLennan et al. 2022). In summary, we created a Phaethon discrete element model with a two-layer structure consisting of a robust HCP core and a polydisperse RCP envelope. Moving forward, we applied this model in the mascon method, approximating the total gravitational potential of the body by summing the gravitational potential of each particle.

Fig. 1 Cross section of the mascon model built using PKDGRAV. The white grid overlay shows the original 3D shape model of Phaethon. The core (blue) and envelope (yellow) particles are distinguished.

3 Results

In this section we present our results on the geopotential in Sect. 3.1, the dynamical state on the surface in Sect. 3.2 and above the surface in Sect. 3.3.

3.1 Geopotential

Using the mascon model, we calculated the geopotential at the centroid of each triangular face of the shape model, representing the geopotential of the face itself. In some faces, the absolute geopotential values were unusually high due to certain particle mass centers being too close to the surface. In these cases, we adjusted the distances from the problematic particles to the surface to be equal to their radii. Figure 2 shows the resulting geopotential map of Phaethon.

To confirm the accuracy of the mascon method, we calculated the geopotential using the polyhedron method for comparison and present the results in the appendix. Figure A.1 shows the relative difference between the two methods, assuming a uniform density that is equal to the bulk density of the mascon model for the polyhedron method. To account for the different structure assumptions of the two methods, we also calculated the polyhedron method by setting the bulk density to that of the mascon envelope (1.40 g/cm³), then replacing the density of the central 1.5 km radius spherical volume with the core bulk density (1.77 g/cm³). The result is shown in Fig. A.2, where one can see that the relative difference between the two methods is reduced.

In both cases, while some of the relative error can be attributed to the mascon resolution, the mascon geopotential tends to be higher in the northern mid-to-high-latitude region compared to the polyhedron method. It is important to note that the mascon model was carved from a larger spherical aggregate, where the particles that are closer to the core would have been more compressed, causing them to be more tightly packed than those near the surface. This feature may have remained after being carved into Phaethon’s shape model, as the higherlatitude regions in the northern hemisphere are radially closer to the core, they correspond to the deeper and denser part of the original sphere model. Nonetheless, our “realistic” mascon model produces geopotential values consistent with the polyhedron method’s analytical solution, with only minor differences likely due to local heterogeneity and model resolution.

From the geopotential, we derived several quantities to understand the possible dynamical state of Phaethon’s surface and near-surface materials, and eventually to consider science cases before the DESTINY⁺ flyby. The rest of this section mainly employs the methods in Scheeres (2012) and Scheeres et al. (2016). Unless otherwise specified, the following results are based on the mascon method calculation.

Fig. 2 Geopotential maps calculated using the mascon method. From left to right and from top to bottom, each panel shows the geopotential map oriented in the −z, −x, −y, +x, +y, and +z directions.

3.2 Surface dynamical state

Based on our model, we calculated the acceleration on the surface, the tilt, and the slope to characterize the dynamical environment on the Phaethon surface. Each result is reported in the following subsections.

3.2.1 Surface acceleration

The gradient of the geopotential ∇V at the position r corresponds to the acceleration due to both gravity and the rotation of the body. By extracting the component orthogonal to the normal surface vector, we calculated the surface acceleration, a_t: $$a_t=\mathrm{\nabla }V-(\mathrm{\nabla }V\cdot \hat{n})\hat{n},$$(2)

where $\hat{n}$ is the surface normal unit vector. Figure A.3 illustrates the direction and magnitude of Phaethon’s surface acceleration, which is higher in the high-latitude regions and lower in the equatorial region, implying that, at the current state, surface materials tend to migrate from polar to equatorial areas.

In realistic conditions, the acceleration of the surface by the geopotential is countered by other forces, particularly friction. For static friction, it is given by $$f=({\mu }_s\mathrm{\nabla }V\cdot \hat{n})\frac{a_t}{\left|a_t\right|},$$(3)

where μ_s denotes the static friction coefficient, which we assumed to be μ_s = 1.0 (see Table 2). Taking friction into account, the net surface acceleration is then given by $$a_{\text{net }}=\{\begin{array}{ll}{a}_t-f,& \text{if}\left|a_t\right|\ge |f|\\ 0,& \text{otherwise}.\end{array}$$(4)

Figure A.4 shows the revised acceleration maps after accounting for friction. On most faces, static friction offsets surface acceleration. In areas where surface acceleration is strong enough to overcome friction, the resulting net surface acceleration reach up to ∼1 × 10⁻⁴ m/s².

In addition to friction, cohesion can help counter surface acceleration. The relation between cohesion c and the cohesive force F_c for a spherical particle of radius r can be approximated as $c\sim \frac{F_c}{2\pi r^2}$ (Sánchez & Scheeres 2014). Therefore, for cohesion to fully counteract the net surface acceleration, $$c\gtrsim \frac{2}{3}\rho s{a}_{\text{net}},$$(5)

where ρ and s are the mass density and radius of a spherical surface object, respectively. Polarimetric observations of Phaethon estimate a characteristic grain diameter of ∼0.36 mm (Ito et al. 2018), while thermal conductivity modeling suggests a range from tens of microns to centimeters (MacLennan et al. 2022). This grain size range is comparable to that observed on Ryugu through polarimetry (Kuroda et al. 2021), as well as the measurements from returned samples of Ryugu and Bennu (Miyazaki et al. 2023; Lauretta et al. 2024). As such, it is not unreasonable to assume that the grain size range on Phaethon is similar to those of Ryugu and Bennu. Based on the returned samples from Ryugu and Bennu, we consider a representative spherical grain particle of 1 mm radius and 1.8 g/cm³ (Miyazaki et al. 2023; Lauretta et al. 2024). In this case, the cohesion needed to counter the net surface acceleration in Fig. A.4 is c ≳ 10⁻⁴ Pa. The estimated cohesion on Ryugu and Bennu is reported to be ≲1 Pa (Arakawa et al. 2020; Walsh et al. 2022). Thus, if Phaethon’s surface exhibits a similar cohesion level, even minimal cohesion of ∼10⁻⁴ Pa should be sufficient to counteract the net surface acceleration.

3.2.2 Slope

On a massive, slowly rotating body like Earth, the geopotential is primarily governed by the gravitational term. As such, the geometric steepness of a surface strongly correlates with the geopotential gradient. However, this is not the case for fast-rotating asteroids, where the centrifugal term of the geopotential is comparable to the gravitational component in the equation. For such bodies, the steepness of a surface does not necessarily indicate larger changes in geopotential. To address this, Scheeres et al. (2016) proposed distinct definitions for the geometric and geopotential orientations of a surface, referring to the former as “tilt” and the latter as “slope.”

Tilt is defined as the angle between the normal vector of the surface and the radial vector from the center of mass to the surface. Meanwhile, slope is the angle between the normal vector of the surface and the geopotential gradient (or geopotential acceleration) vector. Essentially, the concept of slope aligns with the everyday terrestrial intuition that objects tend to move downhill more easily on steeper surfaces.

Figure A.5 displays the slope of each face in Phaethon’s shape model. Slope stability can be evaluated with the “factor of safety” (Iverson 1997; Barnouin et al. 2022). For slope angle θ, a surface can be considered unstable if $$\tan \theta >\tan {\theta }_r{(1-\frac{c}{{\rho }_e|\mathrm{\nabla }V|T\sin {\theta }_r})}^{-1},$$(6)

where θ_r is the angle of friction (or angle of repose), ρ_e is the bulk density of the envelope, and T is the depth of the surface layer that is prone to failure. Using our PKDGRAV model parameters, we find θ_r ∼ 40^∘ (Agrusa et al. 2024). Assuming c = 0, Eq. (6) simply tests whether the slope exceeds the angle of friction; if θ > 40^∘, the surface material on that face is susceptible to movement. For nonzero cohesion, we assume T to be 10 m, similar to Bennu and Ryugu (Jawin et al. 2020; Perry et al. 2022; Barnouin et al. 2022; Arakawa et al. 2020). However, given Phaethon’s larger diameter (∼5 km) compared to Bennu and Ryugu (< 1 km), it is unclear whether this assumption holds, but we adopted this value as an illustrative example. At the upper limit of c = 1 Pa, slopes become unstable when θ ≳ 44^∘. Thus, there is an uncertainty of ∼4^∘ in the slope stability criterion due to cohesion. Additionally, in reality, the angle of friction is influenced by particle geometry, roughness, and size, among other factors (Beakawi Al-Hashemi & Baghabra Al-Amoudi 2018). In fact, the particles in our PKDGRAV model remain still throughout several rotation periods, indicating stability even with zero cohesion. Therefore, we view this analysis not as conclusive evidence of current surface instability but rather as a theoretical framework to identify regions potentially susceptible to mass movement if triggered (Kim et al. 2023).

3.3 Dynamics off the surface

3.3.1 Equilibrium points

We can infer from Fig. 2 that in the tug-of-war between gravity and centrifugal force in the rotating frame, the former is dominant on the Phaethon’s surface. This suggests the existence of a point above the surface where these forces balance each other. Such points, where $\frac{\mathrm{\partial }V}{\mathrm{\partial }r}=0$ are called equilibrium points (EPs; Scheeres 2012; Scheeres et al. 2016). Performing accurate calculations and characterization of EPs using our mascon model was computationally challenging. Therefore, for this section, we used the Minor-Equilibria-NR code (Amarante et al. 2021) and the homogeneous polyhedron method for this section’s EP calculations. Figure 3 summarizes the locations and types of EPs. Their detailed coordinates can be found in Table 3, along with their V values and stability types. The stability type is based on the classification of the topological structure by Jiang et al. (2014) and Amarante et al. (2021). As shown in Sect. 3.1, we recall that the geopotential used in this code would have a <2% deviation from our mascon model on the surface. Nonetheless, we deem the results to be close approximations of the EPs. In summary, there are seven EPs found around Phaethon: one stable point inside Phaethon, three saddles, and three unstable center points outside Phaethon.

3.3.2 Escape speed

If the velocity of a particle is sufficiently high, combined with the rotational velocity of the parent asteroid, the particle can overcome the asteroid’s gravity. Following the definition by Scheeres et al. (1998) and Scheeres (2012), the escape speed is the minimum velocity in the surface-normal direction within the corotating frame required for an object to be immediately ejected on an escape trajectory from the parent body. After converting this definition to the inertial velocity, $v_I=v_{\mathrm{e}\mathrm{s}\mathrm{c}}\hat{n}+(\omega \times \mathit{r})$, and combined with the escape condition $\frac{1}{2}{v}_I^2\ge U$, the assumption that the ejection direction aligns with the surface normal simplifies the equations to a quadratic form. Solving the quadratic equation provides $$v_{\mathrm{e}\mathrm{s}\mathrm{c}}(\mathit{r})=-\hat{n}\cdot (\mathit{\omega }\times \mathit{r})+\sqrt{[\hat{n}\cdot (\mathit{\omega }\times \mathit{r}){]}^2+2U_{max}(\mathit{r})-(\mathit{\omega }\times \mathit{r}{)}^2},$$(7)

where U_max = max [U(r), GM/r], with G being the gravitational constant and M the total mass of the asteroid. The map of the escape speed is displayed in Fig. 4.

Fig. 3 Equilibrium and zero-velocity curves (ZVCs) around Phaethon. Phaethon’s outline is shown with the yellow lines. The ZVCs are represented by gray lines, with larger V values shown with darker hues. The ZVC interval for the left panel was adjusted to enhance the geopotential around the EPs The dotted red lines are the ZVCs of the minimum geopotential outside of Phaethon, also known as the “rotational Roche lobe” (Scheeres 2012).

3.3.3 Rotational Roche lobe and return speed

The 3D surface corresponding to the minimum geopotential EP outside of Phaethon is referred to as the “rotational Roche lobe” (Scheeres 2012). The Roche lobe intersection marks the surface regions where the potential energy is theoretically sufficient for an object to be launched into space. Conversely, the regions enclosed by the Roche lobe indicate where the surface object requires additional energy to achieve escape. This extra energy corresponds to the return speed.

If an object’s initial velocity is below the return speed, it is guaranteed to fall back to the surface. The return speed is given as $$v_{\text{ret}}(\mathit{r})=\{\begin{array}{ll}\sqrt{2[V(\mathit{r})-C^{\ast }]},& \text{ if }V\ge C^{\ast }\\ 0,& \text{ otherwise, }\end{array}$$(8)

where C^* corresponds to the minimum geopotential EP outside of Phaethon (Scheeres 2012; Scheeres et al. 2016). According to Table 3, this value is −3.58 m²/s² at point E5. The results are illustrated in Fig. 5. If the initial velocity is above the return speed but less than the escape speed, the fate of the object depends on its interaction with the EPs. This interaction can result in outcomes ranging from fallback to long-term orbits, or even escape.

3.3.4 Jacobi speed

When an object located at the highest geopotential point on the surface begins to move downslope, it will naturally gain kinetic energy. Assuming no energy was lost due to resisting forces like friction, this velocity corresponds to the Jacobi speed (Scheeres et al. 2016). In other words, the Jacobi speed can be interpreted as the freefall speed from the highest geopotential point to a reference point. It is defined as $$v_J(\mathit{r})=\sqrt{2[V_{max}-V(\mathit{r})]},$$(9)

where V_max is the maximum geopotential on the surface. The result is shown in Fig. A.6. In realistic conditions, the kinetic energy of an object moving across an asteroid’s surface would be affected by factors such as friction and collisions. Nevertheless, the Jacobi speed can act as a useful basis for analyzing the particles’ capability to be lifted off the surface. For instance, the Jacobi speed calculation indicates that an object migrating from the pole to the equator without energy loss can achieve escape speed. Combined with the return speed map in Fig. 5, the Jacobi speed suggests that the most realistic and likely scenario would be for the particles to exceed the return speed and be lifted off the surface after a landslide event. As noted in Sect. 3.3.3, most of Phaethon’s surface exhibits a return speed of zero, meaning that a particle displaced from a high geopotential point could theoretically escape shortly thereafter.

3.3.5 Trajectory of particles outside Phaethon

To summarize, any object with a velocity below the return speed is guaranteed to fall back to the surface, while if the normal velocity is above the escape speed, it is guaranteed to escape the asteroid. However, for velocities that fall between these two thresholds, the object’s trajectory can vary widely, depending on its initial conditions and interactions with EPs (Scheeres et al. 2016). To investigate the dynamical properties of particles once they are lifted off the surface of Phaethon, we conducted a simple simulation in PKDGRAV.

First, we selected 130 outermost surface particles in the equatorial regions, which served as potential launch sites. We then placed 20-meter-radius particles with a one-meter radial offset from each launch site. To minimize their gravitational influence, these particles were given a small density of 0.1 g/cm³. Lastly, each particle was given the same initial velocity in the radially outward direction in addition to the rotational velocity of the corresponding launch site particles. We tested three initial radial velocities: 0.1, 0.5, and 1 m/s. The simulation timestep and parameters are identical to those listed in Table 2. The locations and velocities of the particles were recorded every 1000 timesteps (∼50.2 s). The simulation ran for 26 million timesteps, which corresponds to just over 10 Phaethon rotations.

Although the initial speed of 0.1 m/s is an order of magnitude lower than the escape speed, most areas of Phaethon have zero return speed (see Fig. 5), which is theoretically enough to reach the lowest geopotential EP (E5). However, in our simulation result, all particles could not reach E5 before falling back to the surface. In other words, despite most of Phaethon’s surface having a zero return speed, the particles from these regions can only escape if they can reach E5, which was difficult for most particles in this simulation, as the initial launch direction was radial. To effectively utilize E5 as an escape route even with low velocity, particles far from E5 must be launched with specific initial trajectories aiming toward E5 following the path of least resistance. In the case of particles whose initial launch locations were close to E5, their return speeds exceeded 0.1 m/s. As a result, all simulated particles fell back shortly after launch and only showed brief movements before settling into a stable position.

The initial speed of 0.5 m/s is below the escape speed but exceeds the maximum return speed of Phaethon. At this velocity, the particles have more opportunities to interact with the EPs. The top panel of Fig. 6 exemplifies this interaction; the particle exhibits a complex trajectory after being launched near E6 and appears to engage with E2 and E3 after rebounding. Additionally, the particles gain latitudinal (z-component) motion as well. When the initial velocity increases to 1 m/s, which is just below the escape speed in the equatorial region, this latitudinal motion becomes more pronounced, with some particles reaching beyond Phaethon’s pole, as shown in the case of the bottom panel of Fig. 6. However, based on Fig. 7, when the particles eventually fall back to the surface, they are still concentrated in the low-latitude region within 20^∘. However, we note that this is the result of a simplified simulation setup with limited initial locations and directions. Future studies could explore more detailed scenarios by varying the launch locations and directions, as well as incorporating the effects of solar radiation pressure on small dust particles (Amarante et al. 2021).

Although our simulations did not yield any particles that reach escape velocity through EP interactions, they demonstrated the influences of EPs on the particle trajectories. Figure 8 summarizes the lifetime of the particles, defined as the time taken for each particle to fall back and rest on the surface. Since EPs are products of Phaethon’s irregular shape, we also analyzed a spherical body with the same mass and equivalent radius as Phaethon for comparison. The expected lifetime of this spherical Phaethon is marked with the blue dotted line in Fig. 8. The lifetime around a spherical body is noticeably shorter than that of most particles in our calculation for Phaethon. Moreover, a few particles survive beyond the simulation duration, although they are expected to eventually fall back to the surface, based on their velocities.

Fig. 4 Map of the escape speeds on Phaethon’s surface.

Fig. 5 Map of the return speeds on Phaethon’s surface.

Fig. 6 Trajectory of a particle launched from the surface with an initial radial velocity of 0.5 m/s (a) and 1 m/s (b). The darkness of the line shows the passage of time, from lighter to darker. The left and right panels are the projections of the x − y and x − z planes, respectively.

Fig. 7 Final positions of the launched particles after they return to the surface. The different colors denote the different initial velocities of the particles. Particles that did not return to the surface by the end of the simulation are not shown.

4 Discussion

In this section, we first outline the uncertainties involved in our method. We then review potential activity mechanisms on Phaethon and their possible effects on surface material in light of our findings. Next, we highlight Phaethon’s unique characteristics compared to previously explored near-Earth asteroids (NEAs). Finally, we discuss the implications of this study for DESTINY⁺'s upcoming in situ observations.

4.1 Model uncertainties

4.1.1 Mass density

Although this study uses the nominal mass density 1.58 g/cm³ by MacLennan et al. (2022), the study also reports an uncertainty of 0.45 g/cm³, or about 28%. Figure 9 shows how much the EPs can change within the 1σ uncertainty range (1.13−2.03 g/cm³). First, we highlight the number of EPs decreasing for higher density.

At the lowest density, we found nine EPs, while at the highest density, there are only three. A similar trend was found on Bennu (Scheeres et al. 2016). However, Yu & Hexi (2018) found that increasing the rotation speed resulted in the annihilation of EPs for eight of the ten small bodies they analyzed, including Bennu. Meanwhile, Amarante & Winter (2020) found that lower density leads to fewer EPs for contact binary object Arrokoth. Within this context, we conclude that the number of EPs is sensitive to mass density and shape model.

The change of location for different density assumptions is also an important trend. As the density increases, so does the gravitational attraction of the body. Since the EP is where the centrifugal and gravitational acceleration is balanced, if the gravitational acceleration increases with density, the EPs require stronger centrifugal accelerations. Assuming that the angular velocity is fixed, this leads to the EPs being pushed further outward. Here, we point out that for the lowest density, all the unstable points are near or inside the surface of Phaethon. Hirabayashi & Scheeres (2014) defined this condition as the “first surface shedding.” In other words, at the lower end of the 1σ density range, we can expect Phaethon to be undergoing global deformation and mass shedding. To summarize, within the currently estimated density range of Phaethon, the number and locations of the EPs can vary widely and have substantial consequences for its surface and near-surface environment.

Our assumption of the core-envelope structure is another possible source of uncertainty. As mentioned in Sect. 2, this structure was used to account for the possibility of the Geminids being formed by a mass shedding event caused by rotational instability (Jo & Ishiguro 2024). However, a denser core is not required for mass shedding to occur if the interior is mechanically stronger than its envelope (see Sect. 1.2 of Agrusa et al. 2024, and references therein). A prime example is Bennu, whose surface features suggest a high-cohesion interior (Zhang et al. 2022), despite also being found to have an under-dense core (Scheeres et al. 2020). Nevertheless, given that Phaethon is a larger and denser body than Bennu, we cannot rule out the possibility that their internal structures are different. Furthermore, Ferrari & Tanga (2022) investigated various cases of core-envelope configurations and concluded that density inhomogeneity between the core and envelope has little effect on structural and surface stability, provided the bulk density remains unchanged.

Fig. 8 Histograms of the launched particles’ lifetimes for initial velocities of 0.5 m/s (upper) and 1 m/s (lower). The dotted blue line marks the expected lifetime, assuming Phaethon is a spherical body, while the dotted red line is the mean lifetime of the particles.

4.1.2 Shape model

The shape model has two types of uncertainty that are relevant to this work: size and smoothness. The relative uncertainties of Phaethon’s lengths along its x-, y-, and z-axes are 3%, 3%, and 6%, respectively. Additionally, for radar-based shape modeling, there is a subjective weight that controls how smooth the model will be (Magri et al. 2007). Increasing this weight discourages small-scale topographical roughness at the expense of worse fits to the observational data. The nominal model used in this study is the “rough” model.

4.1.3 Combined effect of the uncertainties

To evaluate the impact of uncertainties in the mass density and shape model on the surface analysis of Phaethon, we conducted additional calculations by varying three parameters: density, model size, and model smoothness. For density, we adjusted the PKDGRAV particle masses to achieve bulk densities of 1.13 and 2.03 g/cm³, corresponding to ±1σ deviations from the nominal density. Similarly, we varied the model size by adjusting the x- and y-axes together by ±3% and the z-axis by ±6% independently. Three levels of smoothness were also tested. Among these parameters, the mass density had the most significant impact, likely due to its larger uncertainty. When the mass density was increased, the unstable areas were considerably reduced, as increased gravity contributed to greater stability. We exhibit the slopes calculated after varying the model parameters in Appendix B. As our nominal results indicate that unstable regions exist on Phaethon, we focused on the most “stable” scenario to determine whether these regions keep their instability amid these uncertainties. The most “stable” case, with the fewest regions showing nonzero net surface acceleration and slopes exceeding the angle of friction, occurred with high mass density, larger size, and minimal smoothness. In this scenario, the net surface acceleration was zero across all facets. Thus, within the current uncertainty range of the parameters, we cannot rule out the possibility that Phaethon may actually be quite stable.

Fig. 9 The equilibrium points of Phaethon assuming different densities, marked with different markers.

4.2 Activity and ejection mechanisms

4.2.1 Possible present-day activity

From Sect. 3.3.2, we note that the escape speed in the equatorial region is ∼1 m/s. This is comparable to the velocity of ejected particles detected on Bennu by OSIRIS-REx (Lauretta et al. 2019). Thermal fatigue is a possible cause of Bennu ejection events (Molaro et al. 2020). While the diurnal and annual temperature variation on Bennu is ∼100 K (Rozitis et al. 2020), Phaethon is expected to experience variations of several factors higher (Yu et al. 2019; Lee et al. 2019; Bach & Ishiguro 2021; MacLennan & Granvik 2024). Hence, it is not unreasonable to expect Phaethonian particles to receive at the very least a similar amount of kinetic energy from such an event. Additionally, the high albedo and thermal inertia of Phaethon (Hanuš et al. 2018) could indicate a history of thermal fractures.

Another candidate for Bennu’s ejection is meteoroid impact (Bottke et al. 2020). Compared to Bennu, due to the small perihelion distance, Phaethon ventures deeper into the inner Solar System, where the meteoroid population is denser. This suggests that meteoroid impacts on Phaethon could be more frequent. Additionally, such impacts would be more energetic on Phaethon as its highly eccentric orbit leads to high relative velocities of interplanetary dust particles (Szalay et al. 2019). Especially near the perihelion, an impact by centimeter-sized particles can lead to an ejection of thousands of kilograms of material (Zubko et al. 2022). In short, if the same events that cause particle ejections on Bennu occur on Phaethon, it is possible for particles to reach escape speed, particularly at lower latitudes.

Zhang et al. (2023) suggested that the present-day brightening of Phaethon in the perihelion passage is due to sodium volatilization. They estimated that only small dust particles (micron-scale or less) would be able to escape from Phaethon from this mechanism. Additionally, recent thermal modeling suggests that thermal decomposition can also lift dust particles from Phaethon (MacLennan & Granvik 2024).

In short, present-day dust activity on Phaethon can be driven by thermal fracture, meteoroid impact, and/or sodium volatilization. While it would be challenging for DESTINY⁺ to directly observe such small particles, the takeaway should be that regolith particles are disturbed by such activity every perihelion passage. If the force exerted by any of these mechanisms is enough to overcome friction and displace the regolith material, it could trigger downslope mass movements, such as landslides where boulder-sized objects can be affected as well. From Fig. A.6, one can infer that if sodium activity displaces some mass at mid to high latitude (high geopotential region), by the time this mass reaches the equatorial region, it may reach a velocity comparable to the escape speed.

4.2.2 Geminid formation event

Whether these mechanisms caused the Geminid meteor stream in the past is a separate question. Key constraints from the Geminid stream are its estimated mass range (∼10¹²−10¹⁷ g; Hughes 1983; Jenniskens 1994; Blaauw 2017; Ryabova 2017; Jo & Ishiguro 2024), age range (∼2−18 kyr; Ryabova 1999, 2007; Jo & Ishiguro 2024), and the existence of large particles of ∼1−10 cm in size (Blaauw 2017). Based on Bennu’s ejection, thermal fracture and meteoroid impacts can produce cm-sized ejecta. However, for all the mechanisms mentioned above, it is unclear if they are able to produce the Geminid stream mass. It is also possible that the Geminid formation mechanism is currently inactive. Historically, comet-like ejection by ice sublimation was the first proposed activity mechanism (Gustafson 1989). However, during its millions of years in the inner Solar system, it is unlikely for water ice in Phaethon to survive long enough to satisfy the Geminid age range (Jewitt & Li 2010; Yu et al. 2019; MacLennan et al. 2022). On the other hand, Geminid formation by YORP-induced rotational instability was suggested by Nakano & Hirabayashi (2020). Jo & Ishiguro (2024) performed dynamical simulations and found that particle ejection by rotational instability that occurred ∼18 kyr ago can create the Geminid meteor shower.

The results of our study indicate that rotational instability is not likely to be occurring at present on Phaethon. However, we cannot rule out the possibility that Phaethon experienced rotational instability in the past, undergoing global mass migration and even mass ejection events before decelerating to its current rotation (Nakano & Hirabayashi 2020). If this event is the origin of the Geminid stream, most of Phaethon’s surface would consist of freshly exposed material. We also note that three decades of light curve data indicate that Phaethon’s rotation is currently accelerating (Marshall et al. 2022). The YORP effect predicts a slower acceleration, but the observed acceleration could be explained by tangential YORP, which depends on small-scale morphology and cannot be calculated accurately from the current shape model (Senshu et al. 2025). Phaethon may undergo cycles of spin-up, spin-down, and spin-up again on timescales shorter than the YORP timescale.

Fig. 10 Stability of Phaethon’s surface. The conditions of instability are (i) a nonzero net surface acceleration and (ii) a slope steeper than 40∘. Purple, green, and yellow represent none, one, and both conditions satisfied, respectively, for each face.

4.3 Comparison with other asteroids

While the recent explorations of NEAs Bennu, Ryugu, and Didymos have provided various new insights into the evolution and structure of asteroids, our results show that Phaethon is unique in this context. Firstly, Phaethon’s large size yields high surface acceleration. The maximum surface acceleration found on Ryugu and Didymos is ∼1.5 × 10⁻⁴ m/s, and ∼8 × 10⁻⁵ m/s for Bennu (Watanabe et al. 2019; Yu et al. 2019; Scheeres et al. 2019). On Phaethon, we estimate the maximum surface acceleration to be 6.4 × 10⁻⁴ m/s, which is up to several factors higher than the aforementioned asteroids. In other words, the Phaethon’s surface should theoretically be more mobile.

In particular, there is one feature on Phaethon where the surface acceleration tends to be highest. Figure 10 summarizes the results in Figs. A.4 and A.5 and identifies the most unstable regions. We note, as explained in Sect. 4.1, that the results are subject to uncertainties. Therefore, we emphasize that Fig. 10 serves as a visual summary of the nominal results and is intended to highlight regional variations in instability rather than provide definitive indicators of instability. Due to the low resolution of the model, regions where only one or two faces are marked unstable should be considered unreliable. However, the large depression at the mid-latitude of the northern hemisphere (marked with the white arrows in Fig. 10) shows a concentration of unstable faces. This region occupies ∼2.3 km², with ∼0.8 km² of it being possibly unstable (yellow or green in Fig. 10), and correspond to the feature (d) identified in the radar observations by Taylor et al. (2019). The slopes in these regions can reach 51.4^∘ and have some of the highest net surface accelerations as well. Although these calculations strongly indicate that this region is unstable, this does not guarantee mass movement in this region. In fact, our PKDGRAV model remained stable for multiple Phaethon rotations, suggesting that despite the steep slope, it is possible for Phaethon to maintain this cliff-like structure. However, we remind that our PKDGRAV particles are more than 10 m in radius. Thus, while boulder-sized objects can sustain this cliff, that might not be the case for grains and dust particles. We conjecture that the unstable structure consists of exposed boulders, while regolith particles moved to lower geopotential areas just south of this structure, forming a deposit. If this deposit was formed by a local landslide event, its shape would follow the angle of repose (Clark et al. 2023). Considering that the angle of repose is typically ∼30−40^∘, the slope at the southern rim of the depression feature is around this range, implying that this is the deposit section in this region. In short, we argue that Phaethon has a large depression consisting of cliffs and deposits, which is a unique geological feature that was not found on recently explored NEAs. The cliffs may have exposed subsurface boulders, which would potentially give us more insight into the interiors of rubble-pile asteroids, and from the shape of the deposits, it would be possible to constrain the regolith material properties to some degree. As DESTINY⁺ is a flyby mission, it will only be able to observe about half of Phaethon’s surface and has to select the observation region by adjusting the timing of the closest approach. We propose that the northern hemisphere depression be considered as a potential focus during this selection process.

4.4 Implications and expectations for DESTINY⁺

When an airless body is exposed to solar winds and micrometeorite bombardments, its surface can undergo physical and/or chemical modifications that result in optical changes over time a process known as space weathering (Clark et al. 2002). Two recently visited top-shaped asteroids, Bennu and Ryugu, both displayed latitudinal color variation, indicating different degrees of space weathering due to global mass movement (Yumoto et al. 2024). However, space weathering on C-complex asteroids is particularly complicated.

On Bennu, smaller (and presumably younger) craters appear redder than its surroundings, implying that Bennu’s surface becomes bluer over time (Deshapriya et al. 2021; Lauretta et al. 2022; Clark et al. 2023). However, Bennu’s equatorial region, which is thought to be the oldest regions on Bennu due to its low geopotential (Scheeres et al. 2019) and high crater density (Bierhaus et al. 2022), exhibits redder and darker characteristics (Li et al. 2021). Based on these observations, Clark et al. (2023) concluded that space weathering on Bennu leads to nonlinear color variations, from red to blue and back to red. However, DellaGiustina et al. (2020) and Li et al. (2021) noted that the equatorial region appears bluer in the blue-to-visible wavelengths, potentially due to 0.55 μm enhancement from magnetite production. On Ryugu, latitudinal variations and the colors of young craters indicate that space weathering results in the reddening and darkening of surface materials (Sugita et al. 2019; Morota et al. 2020; Tatsumi et al. 2021). Yumoto et al. (2024) cross-calibrated the images from Hayabusa2 and OSIRIS-REx to compare the space weathering trends of the two asteroids. They found that Ryugu and Bennu may have started with similar spectra but underwent different weathering processes, possibly due to variations in surface dust size and coating thickness.

Based on the space weathering trends found by these previous missions, we can conjecture what DESTINY⁺ may find on Phaethon. Figures A.3 and A.5 suggest that mass movement from mid-to-low latitudes should be expected on Phaethon, similar to what was observed on Bennu. Using the Telescopic CAmera for Phaethon (TCAP) on board DESTINY⁺, boulder orientations from surface images will be a strong indicator of such a global mass movement trend on Phaethon (Jawin et al. 2020). Additionally, it is reasonable to infer that the equatorial region is the most static and, therefore, likely the most mature and space-weathered. In contrast, mid-latitude regions would be less space-weathered and thus fresher. This latitudinal difference in space weathering can be revealed through multiband observations with the Multiband CAmera for Phaethon (MCAP) on DESTINY⁺. The space weathering timescale is typically thought to be ≲ 10⁶ years for main-belt S-type asteroids (e.g., Vernazza et al. 2009), and as short as 10³ years based on a recent analysis on the Hayabusa images of NEA Itokawa (Jin & Ishiguro 2022). For C-complex asteroids, Hasegawa et al. (2024) estimated the space weathering timescale for main-belt Ch and Cgh asteroids to be ∼10^3.5−10^4.5 years. If we scale this value to Phaethon’s orbit, following $T_{\mathrm{s}\mathrm{w}}\propto a^2\sqrt{1-e}$ (Hapke 2001), where a and e are the semimajor axis and eccentricity, the timescale decreases to ∼10^2.3−10^3.3 years, shorter than the estimated age of the Geminid stream, ∼10³−10⁴ years (Ryabova 2007, 2017; Jo & Ishiguro 2024).

On the opposite side of space weathering is resurfacing, which refers to the various processes, such as tidal interaction, impact-induced seismic shaking, thermal fatigue and YORP spin-up, that expose fresher material from below the surface (Jin & Ishiguro 2022). The potential past and present activity mechanisms, discussed in Sect. 4.2, can all be considered as Phaethon’s resurfacing mechanisms. Due to Phaethon’s fast rotation, these mechanisms may trigger landslides and launch materials into space. The equatorial region is prone to potential fallback impacts from objects that were briefly launched into space (Fig. 7). If Phaethon’s equator is highly populated with boulders, it could indicate that such launch-and-fallback events, and large landslides that caused these events, happen frequently.

Phaethon’s unique status as the only kilometer-sized asteroid with a very short perihelion distance implies that the tension between space weathering and resurfacing will be vigorous near the perihelion. The near-perihelion phase is where Phaethon is subjected to intense solar wind and denser regions of the zodiacal cloud, leading to stronger space weathering. Conversely, it is also where Phaethon experiences extreme temperature variations and more frequent meteoroid impacts, potentially making the resurfacing process more efficient on Phaethon than on Bennu. A further complication at perihelion is sintering. Polarimetric studies by Ito et al. (2018) and Geem et al. (2022) suggest that Phaethon is populated with large grains, possibly produced by sintering near the perihelion, although its role in space weathering and resurfacing remains unclear.

In summary, on a global scale, since the space weathering timescale is shorter than the average age of the surface, which we assume to be equal to the Geminid age, Phaethon’s surface should be at or near saturation in terms of weathering, in the absence of resurfacing events. However, the near-Sun, high-eccentricity orbit of Phaethon makes it susceptible to multiple resurfacing mechanisms. Additionally, Fig. A.3 indicates that the regolith movement might be more efficient on Phaethon compared to other visited NEAs due to its larger size, supporting the likelihood of latitude-dependent variations in space weathering. On a local scale, as discussed in Sect. 4.3, Phaethon features a unique depression geology at the northern hemisphere. The instability of this area implies that it is among the freshest regions on Phaethon, providing an additional reason for it to be strongly considered as part of the DESTINY⁺ observation region.

5 Summary and conclusions

We have presented the first dynamical analysis of the DESTINY⁺ mission target (3200) Phaethon. Using the latest shape model generated from various observations, we computed Phaethon’s gravitational field using the mascon method, filling the model with particles via the SSDEM N-body code PKDGRAV. Next, we calculated the geopotential and derived surface properties such as the surface acceleration, net surface acceleration, and slope. Additionally, we determined EPs, escape speeds, return speeds, and Jacobi speeds to evaluate conditions under which material would escape from Phaethon’s surface. Our findings indicate that even small triggers could induce landslides on Phaethon, potentially launching surface material. According to our numerical simulations, particles launched with velocities between the return speed and escape speed could remain suspended above the surface for extended periods due to interactions with EPs before eventually returning to the surface. We also highlighted a prominent depression in Phaethon’s northern hemisphere; we hypothesize that it could be a cliff deposit structure and suggest it is one of the freshest regions. However, our results are influenced by uncertainties, particularly that of the mass density. We hope that the DESTINY⁺ mission will help reduce these uncertainties and offer further insight into the space weathering trends of carbonaceous asteroids, as well as the resurfacing and activity mechanisms on near-Sun asteroids.

Acknowledgements

This research at Seoul National University was supported by a grant from the Korean National Research Foundation (NRF) (MEST) funded by the Korean government (No. 2023R1A2C1006180). SEM was supported by NASA grant number 80NSSC19K0523 awarded to University of Central Florida. HJ and MI thank Yoonsoo Bach, Sunho Jin, Jooyeon Geem, Bumhoo Lim and Seungyoo Lee for their helpful comments and discussions, as well as Prof. Hyung Mok Lee, Prof. Woong-Tae Kim and Dr. Hee Il Kim for allowing us the usage of the computational facility gmunu.

Appendix A Additional figures

Here we present the figures that show the relative differences of geopotential between the mascon and polyhedron models (Figs. A.1 and A.2) and the map of the surface accelerations (Figs. A.3 and A.4), slopes (Fig. A.5), and Jacobi speeds (Fig. A.6).

Fig. A.1 Maps of relative difference in the surface geopotential calculated using the mascon method in comparison with the polyhedron method. We assume a uniform density for the polyhedron method calculation. The apparent directions in each panel are the same as in Fig. 2.

Fig. A.2 Maps of relative difference in the surface geopotential calculated using the mascon method in comparison with the polyhedron method. Here, the polyhedron method with the core-envelope structure is used. The apparent directions in each panel are the same as in Fig. 2.

Fig. A.3 Maps of the surface acceleration on the Phaethon surface.

Fig. A.4 Maps of the net acceleration on the Phaethon surface.

Fig. A.5 Slope maps of the Phaethon surface

Fig. A.6 Map of the Jacobi speeds on Phaethon’s surface. The value corresponds to the expected velocity of an object when migrated from the highest geopotential point without energy loss.

Appendix B Model parameter influence on stability

Table B.1 lists the minimum, mean and maximum slopes of Phaethon after varying the bulk density, size and smoothness of the model. We note that a slope value of more than 90^∘ indicates that cohesionless material on the region will be detached from the asteroid (Yu et al. 2019).

References

All Tables

All Figures

	Fig. 1 Cross section of the mascon model built using PKDGRAV. The white grid overlay shows the original 3D shape model of Phaethon. The core (blue) and envelope (yellow) particles are distinguished.
In the text

	Fig. 2 Geopotential maps calculated using the mascon method. From left to right and from top to bottom, each panel shows the geopotential map oriented in the −z, −x, −y, +x, +y, and +z directions.
In the text

Fig. 3 Equilibrium and zero-velocity curves (ZVCs) around Phaethon. Phaethon’s outline is shown with the yellow lines. The ZVCs are represented by gray lines, with larger V values shown with darker hues. The ZVC interval for the left panel was adjusted to enhance the geopotential around the EPs The dotted red lines are the ZVCs of the minimum geopotential outside of Phaethon, also known as the “rotational Roche lobe” (Scheeres 2012).

In the text

	Fig. 4 Map of the escape speeds on Phaethon’s surface.
In the text

	Fig. 5 Map of the return speeds on Phaethon’s surface.
In the text

	Fig. 6 Trajectory of a particle launched from the surface with an initial radial velocity of 0.5 m/s (a) and 1 m/s (b). The darkness of the line shows the passage of time, from lighter to darker. The left and right panels are the projections of the x − y and x − z planes, respectively.
In the text

	Fig. 7 Final positions of the launched particles after they return to the surface. The different colors denote the different initial velocities of the particles. Particles that did not return to the surface by the end of the simulation are not shown.
In the text

	Fig. 8 Histograms of the launched particles’ lifetimes for initial velocities of 0.5 m/s (upper) and 1 m/s (lower). The dotted blue line marks the expected lifetime, assuming Phaethon is a spherical body, while the dotted red line is the mean lifetime of the particles.
In the text

	Fig. 9 The equilibrium points of Phaethon assuming different densities, marked with different markers.
In the text

	Fig. 10 Stability of Phaethon’s surface. The conditions of instability are (i) a nonzero net surface acceleration and (ii) a slope steeper than 40∘. Purple, green, and yellow represent none, one, and both conditions satisfied, respectively, for each face.
In the text

	Fig. A.1 Maps of relative difference in the surface geopotential calculated using the mascon method in comparison with the polyhedron method. We assume a uniform density for the polyhedron method calculation. The apparent directions in each panel are the same as in Fig. 2.
In the text

	Fig. A.2 Maps of relative difference in the surface geopotential calculated using the mascon method in comparison with the polyhedron method. Here, the polyhedron method with the core-envelope structure is used. The apparent directions in each panel are the same as in Fig. 2.
In the text

	Fig. A.3 Maps of the surface acceleration on the Phaethon surface.
In the text

	Fig. A.4 Maps of the net acceleration on the Phaethon surface.
In the text

	Fig. A.5 Slope maps of the Phaethon surface
In the text

	Fig. A.6 Map of the Jacobi speeds on Phaethon’s surface. The value corresponds to the expected velocity of an object when migrated from the highest geopotential point without energy loss.
In the text