The Secular Dynamical Evolution of Binary Asteroid System (65803) Didymos Post-DART

The successful impact of NASA’s DART mission with Dimorphos, the secondary body of binary asteroid system (65803) Didymos, altered the attitude, shape, and orbit of Dimorphos. In addition to perturbing the immediate short-term dynamics of the system, these changes have major implications for the binary Yarkovsky–O’Keefe–Radzievskii–Paddack (BYORP) effect and resulting long-term secular evolution. In this study, we assess the range of possible reshaping-induced changes in BYORP. We produce high-fidelity numerical simulations of the orbit-attitude coupled tidal-BYORP dynamical evolution of the Didymos system to constrain how the secular evolution changed from its preimpact behavior. We find that the nature of the dynamics is highly dependent on a variety of initial conditions and assumptions, and it is difficult to fully predict how the system will secularly evolve following the impact. Rather, we provide a range of feasible possibilities within the bounds of observations and current best estimates of the Didymos system parameters. ESA’s Hera mission will help reduce uncertainties surrounding the postimpact state and shape of Dimorphos in 2027, allowing our predictions of the secular evolutionary effects and long-term fate of the Didymos system to be better refined.


Introduction
On 2022 September 26, the first-ever planetary defense demonstration of a kinetic impactor-the NASA Double Asteroid Redirection Test (DART) mission-successfully impacted the smaller, secondary body of the near-Earth binary asteroid system (65803) Didymos, Dimorphos (Daly et al. 2023a).The spacecraft struck Dimorphos with a relative speed of about 6.15 km s −1 (Daly et al. 2023a), resulting in an estimated along-track orbital velocity component reduction (Δv T ) in the mutual orbit of approximately 2.7 mm s −1 (Cheng et al. 2023).DART impacted the leading hemisphere of Dimorphos (roughly along its intermediate axis), reducing the mutual orbit period by approximately 33 minutes and shrinking the semimajor axis of the mutual orbit (Thomas et al. 2023).The main measure of DART's effectiveness in deflecting an asteroid is the momentum enhancement factor, β, which quantifies how much momentum is transferred from the incident momentum of DART to Dimorphos through the impact.The Δv T estimate from Cheng et al. (2023) indicates a β in the range of 2.2-4.9 depending on the unknown mass of Dimorphos, with a nominal value of 3.61 0.25 0.19 -+ (1σ) assuming Didymos and Dimorphos have equal densities of 2400 kg m −3 .This higher β value suggests that more momentum was transferred to Dimorphos from the ejecta induced by the impact, rather than the impact itself (Cheng et al. 2023).The impact also altered the rotation state of Dimorphos.Prior to impact, it was assumed to be synchronous with the mutual orbit (i.e., its spin period equals the mutual orbit period; Richardson et al. 2022).Depending on the exact postimpact shape of Dimorphos, the impact either excited libration on the order of several tens of degrees or triggered attitude unstable tumbling (Nakano et al. 2023a;Richardson et al. 2023).Relevant preand postimpact parameters of the Didymos system are listed in Table 1.
One of the most curious results surrounding the DART mission is that of both the pre-and postimpact shape of Dimorphos.Prior to impact, little was known regarding its shape or surface morphology (Daly et al. 2023a).It was initially believed to be an elongated ellipsoid, as that is the case for the majority of observed near-Earth binary systems (Pravec et al. 2016).However, images taken by the onboard narrowangle DRACO camera moments before impact unveiled an oblate spheroid shape, differing from observations of known prolate rubble-pile secondaries (Daly et al. 2023a).The preliminary oblate spheroid shape model of Dimorphos was originally reported in Daly et al. (2023a), with the most up-todate shape model presented in Daly et al. (2023b).The updated shape model has a volume-equivalent diameter of 150 ± 2.5 m, with extents of 179 ± 1, 169 ± 4, and 115 ± 1 m along the x-, y-, and z-axes, respectively (Table 1).This corresponds to preimpact major-to-intermediate and intermediate-to-minor axis elongation ratios of a/b ≈ 1.06 and b/c ≈ 1.47, respectively.
Initial hypotheses predicted the formation of a crater as a result of the DART impact (Stickle et al. 2017;Rainey et al. 2020;Stickle et al. 2020Stickle et al. , 2022)).However, recent light-curve observations from Pravec et al. (2023) show a postimpact elongation of a/b = 1.34 ± 0.1, suggesting that global deformation due to the impact could have taken place.(Note that these light-curve results could also suggest that Dimorphos entered a tumbling state or is both tumbling and reshaped.)This type of reshaping mechanism for Dimorphos has been predicted by Raducan & Jutzi (2022), falling in the realm of subcatastrophic collisions, a regime that exists between impact cratering and catastrophic disruption of the asteroid (i.e., where more than half of the original mass of the body is lost; Jutzi 2019).More recent work by Raducan et al. (2023) investigates this reshaping mechanism using best guesses of rubble-pile material properties and boulder size-frequency distributions of Dimorphos to match DART observations.They find in their impact simulations that global deformation of Dimorphos is the most likely scenario, with an elongation ratio upper limit of a/b ≈ 1.2.
It has been shown that variations in the shape of Dimorphos can have significant effects on the short-term dynamics of the Didymos system immediately following the impact.Agrusa et al. (2021) shows that the stability of Dimorphos's postimpact attitude is dependent on its elongation ratios, and Nakano et al. (2022) illustrates how reshaping of Dimorphos can cause additional changes in the orbit period postimpact.In addition to these short-term dynamical effects, the DART impact itself and any shape change induced by the impact will also greatly influence the long-term secular effects of the Didymos system.
One of the most important phenomena in understanding the long-term dynamical evolution of rubble-pile binary asteroids is the binary Yarkovsky-O'Keefe-Radzievskii-Paddack (BYORP) effect.This effect occurs when solar radiation pressure (SRP) interacts with a secondary body in synchronous rotation.Asymmetries in the secondary's shape lead to a timevarying SRP force imparted on the body, resulting in secular changes to the binary mutual orbit on timescales of 10 4 -10 5 yr (Ćuk & Burns 2005).The BYORP effect can either expand or contract the mutual orbit depending on the exact geometry of the system.Measurements of the mean anomaly drift for the Didymos system report a value of 0.15°± 0.14°yr −2 (3σ), implying an inward orbital drift (Scheirich & Pravec 2022).This suggests that BYORP is already active on the system preimpact, adding to the justification for this study.While mainly studied as a secular effect, BYORP is in fact highly coupled with the attitude of the secondary (Ćuk & Burns 2004;McMahon & Scheeres 2010a, 2010b;Ćuk & Nesvorný 2010;Quillen et al. 2022).The current BYORP theory outlined later in this paper also depends on the shape of the secondary.Thus, the attitude and shape changes caused by the DART impact will inevitably affect these long-term secular trends.
In this paper, we aim to constrain how the long-term dynamical evolution of the Didymos system was affected by the DART impact.This is achieved by exploring the parameter space within the range of possible postimpact shape changes to bound our predictions.End cases within this parameter space are driven by the strength-dominated cratering and gravitydominated global deformation regimes.In the strengthdominated regime, crater size is driven by Dimorphos's material strength (Stickle et al. 2022), while in the gravitydominated regime, it is driven by its weak gravity, resulting in more displacement and deformation of material (Hirabayashi et al. 2022;Raducan & Jutzi 2022).This work considers a range of possibilities because there is still large uncertainty surrounding the material properties of Dimorphos, so it is not possible to know exactly how Dimorphos's shape changed from the impact (Raducan et al. 2023).This will be refined later when the European Space Agency's Hera spacecraft visits Didymos in 2027 and fully characterizes the postimpact shape and dynamical state of Dimorphos (Michel et al. 2022).This paper first reviews the current BYORP theory and its interaction with tidal dissipation in Section 2. We introduce the dynamical model used in this work in Section 3 and our approach for generating shape models in Section 4. The results showing both shape-change-induced effects on BYORP and high-fidelity simulations of long-term tidal-BYORP evolution are presented in Section 5. A discussion of the implications of these results is found in Section 6, and we summarize the main takeaways from this work in Section 7.

BYORP Theory
The BYORP effect is driven by thermal reradiation of SRP from an asymmetric secondary of a binary system (Ćuk & Burns 2005).Generally for near-circular orbits, BYORP will be stronger for secondaries that are either perfectly or nearly synchronous with the mutual orbit (also referred to as 1:1 spinorbit resonance, or "tidally locked").This is because the fixed relative geometry of the body causes any imbalances in the reradiated SRP force to become pronounced over time and result in a net torque on the orbit (Ćuk & Nesvorný 2010).Small, nonzero librations in the secondary's attitude dynamics can still cause BYORP to affect the secular trends of the mutual orbit, but nonprincipal axis rotation can significantly decrease the BYORP torque (Quillen et al. 2022).BYORP stops acting on the system once the secondary starts circulating, rolling, or tumbling (Ćuk et al. 2021).This is due to the fact that nonsynchronous rotation results in a random orientation of the secondary at any point in time, thus averaging any torques due to BYORP to zero (Ćuk & Burns 2005).Note that here we refer to "libration" simply as the angular deviation of the body's long axis from the primary-secondary line of centers.
The current BYORP theory is fully described in McMahon & Scheeres (2010a) and summarized in this section.The theory models SRP force as a Fourier series over the secondary's rotation because of the periodic nature of SRP from synchronous rotation.Originally derived in Scheeres (2007), the Fourier series is expressed as where N is the number of surface elements used to describe the secondary shape model and n is the order.The quantity P(R) is the SRP defined as where R is the distance from the Sun and F e = 1 × 10 14 kg km s −2 is the solar radiation constant.The quantity f i describes the interaction between SRP and each of the individual surface elements, i, of the shape model.It is expressed as where n i ˆis the surface normal unit vector for each facet; A i is the facet surface area; s i and ρ i are the specular and diffuse reflectivity coefficients, respectively; B L is the Lambertian scattering coefficient of the surface; U is the identity dyad; and ˆ) is the visibility function that accounts for self-shadowing (it equals 1 when the Sun is above the local horizon of a given facet and 0 otherwise).
The quantity u ˆis the position vector of the Sun in the bodyfixed frame of the secondary and is written in terms of the solar latitude, δ s , and solar longitude, λ s : { ˆˆˆ} is the body-fixed frame of the secondary.The solar position vector, u ˆ, is illustrated in Figure 1.From the perspective of the secondary's body-fixed frame, the solar position u ˆwill trace out the surface of a cone, as the Sun will appear to move in a circle as the body rotates.The upper and lower limits of δ s are driven by the solar inclination of the heliocentric orbit, so δ s will vary as the year progresses.
The quantities A n and B n are the Fourier coefficients derived from the SRP expression.They are computed as More details regarding the computation of these coefficients are provided in Scheeres (2007).These coefficients have units of area and are functions of the secondary's shape, reflective properties, and solar latitude in the body-fixed frame.Thus, a set of the Fourier coefficients can be computed provided a shape model of the secondary.The coefficients are vector quantities described in the body-fixed frame, so the along-track A 0 (2) vector component is the dominating term driving secular evolution.The sign of this term provides insight into the direction of evolution driven by BYORP.A positive A 0 (2) coefficient indicates expansive BYORP, while a negative A 0 (2) coefficient indicates contractive BYORP.SRP can also induce torques on the secondary.Similarly to the SRP force, this torque can also be represented by a Fourier where m i = r i × f i is the SRP moment generated by each facet in the shape model, and C n and D n are the no-prime Fourier moment coefficients (Scheeres 2007).They are computed as The term r i is the position vector from the asteroid's center of mass to the surface element, i.
The Fourier coefficients can be normalized for easier comparison of the influence of BYORP across different asteroid shapes.Since the coefficients have units of area, they can be normalized by the mean radius of the secondary body (r m ) squared.The mean radius is defined as the radius of a volume (V ) equivalent sphere, and thus the normalized coefficients can be computed as Normalized 2 0  ( ) magnitudes-often denoted as simply Btypically range from 0 (for a perfectly smooth, symmetric body) to 10 −3 -10 −2 with a higher value corresponding to a rougher and/or more asymmetric body (Jacobson & Scheeres 2011a).For example, Squannit, the secondary of the Moshup near-Earth asteroid (NEA) binary system, has a BYORP coefficient of 2.082 × 10 −2 (McMahon & Scheeres 2010b; Scheirich et al. 2021).

Tidal-BYORP Equilibrium
In addition to the BYORP effect, the other key phenomenon influencing rubble-pile binary evolution is tidal dissipation.In a two-body system, tides are raised on both bodies due to mutual gravitational interactions.The gravity gradient across one body results in stronger forces acting on the side of the body facing the other (tide-raising) body.This results in distortion of the first body, i.e., a tidal bulge (Murray & Dermott 1999).Nonelastic deformation caused by the tidal bulge leads to energy dissipation in that body.While debated in the literature, it is often assumed that this deformation is driven by friction from mechanical stressing (Goldreich & Sari 2009;Nimmo & Matsuyama 2019).In a system where both the primary and secondary asteroids are initially nonsynchronous, tides dissipate rotational energy from both bodies.The first result from tidal dissipation is that the tides raised on the secondary by the primary drive the secondary's rotation rate toward synchronicity (Goldreich & Sari 2009).Once the secondary is tidally locked, tides raised on the primary by the secondary will continue to dissipate energy and work to make the primary synchronous.However, the primary body often rotates much faster than the secondary in these systems, as they are typically close to their spin disruption limit, so this synchronization process takes longer (Goldreich & Sari 2009).While tidal dissipation is spinning down the primary, angular momentum is transferred from the primary to the mutual orbit in order to conserve the total angular momentum of the system, resulting in expansion of the mutual orbit semimajor axis (Goldreich & Sari 2009;Jacobson & Scheeres 2011a).The tidal dissipation process ends when the primary becomes synchronous with the mutual orbit period and the system has reached the doubly synchronous end state (Jacobson & Scheeres 2011a).The two main tidal parameters driving evolution rates are k 2 and Q (Goldreich & Sari 2009).The tidal Love number, k 2 , is the ratio of additional potential due to tidal deformation to the initial tide-raising potential.The tidal quality factor, Q, is essentially a measure of how effective a body is at dissipating energy through tides (Jacobson & Scheeres 2011a).However, there is a large amount of uncertainty regarding the values of these tidal parameters for rubble piles.
Studying how BYORP interacts with tidal effects provides valuable insight into the overall rubble-pile life cycle and possible end states of NEA binaries (Jacobson 2014).It has been hypothesized in Jacobson & Scheeres (2011a) that a tidal-BYORP equilibrium can exist for rubble-pile binaries where tidal expansion is counteracted by contractive BYORP, such that the mutual orbit is not secularly evolving.This is referred to as "joint opposing evolution" (Jacobson & Scheeres 2011a).At distances inside this equilibrium point, tides dominate and expand the orbit to the equilibrium point.At distances outside this equilibrium point, BYORP dominates and shrinks the orbit to the equilibrium point.If a system experiencing joint opposing evolution is at this equilibrium point, the following relationship between tides and BYORP holds: where ρ is density, G 4 3 d w p r = is the tidal disruption limit assuming a sphere, R p is the radius of the primary, q is the system mass ratio, a is the semimajor axis in terms of primary radii, and H e is a solar parameter expressed as with F e being the solar radiation constant and a e and e e being the heliocentric semimajor axis and eccentricity, respectively.The term k p indicates the tidal Love number for the primary.If a system is in such equilibrium, then the unknown parameters B, Q, and k p in Equation (13) can be constrained.Note that B here is referring to the 2 0  ( ) coefficient.If the tidal parameters were known, then B could be determined; alternatively, if we were confident in our theoretical value of B for a given shape model, then k p /Q could be constrained.
The mean anomaly drift of the Didymos system has been measured with a 3σ uncertainty to be ΔM d = 0.15°± 0.14°yr −2 , corresponding to a semimajor axis drift of a 0.09 0.08 cm = -  yr −1 (Scheirich & Pravec 2022).The negative drift in semimajor axis indicates that BYORP must be acting on the system, as this cannot be attributed to tidal expansion.The small value of ΔM d with a near-zero 3σ uncertainty suggests that the system is either in or very close to this tidal-BYORP equilibrium prior to impact.

Full Two-body Problem
We leverage the full two-body problem (F2BP) in order to numerically simulate the long-term dynamical evolution of the Didymos system due to BYORP and tidal effects.The F2BP is able to capture the unique dynamics of the mutual orbit of a binary asteroid system.There is strong coupling between orbital and attitude dynamics in these systems due to the irregular shapes of the bodies and the close proximity between the bodies, leading to non-Keplerian motion (Scheeres 2006).The F2BP has been used extensively in the literature in studies such as Agrusa et al. (2020Agrusa et al. ( , 2021)), Davis & Scheeres (2020), and Meyer et al. (2023b, 2023a).
The dynamical model used in this work is illustrated in Figure 2. Since we are interested in the long-term evolution of the Didymos system, i.e., tens to hundreds of thousands of years, several assumptions and simplifications are made in order to improve computational efficiency.The first simplification is the use of a spherical primary, allowing its equations of motion to be decoupled from the dynamics of the mutual orbit and the secondary's attitude.We also assume that the primary is experiencing principal axis rotation throughout the entirety of its evolution, which is realistic due to the rapidly rotating nature of typical primaries in near-Earth binary asteroids (Ćuk & Nesvorný 2010).Thus, the general forms of the F2BP equations of motion are where r is the position vector of the primary relative to the secondary, ω B/N is the angular velocity of the secondary's body-fixed frame with respect to an arbitrary inertial frame, U is the mutual gravitational potential, [I B ] is the inertia tensor of the secondary, and m is the reduced mass, defined as Note that the subscript A refers to the primary and the subscript B refers to the secondary.The quantities F and M represent any generic forces and moments acting on the system, respectively.Note that these equations are defined in the secondary's bodyfixed frame, so dot derivative quantities represent time derivatives as seen by the body frame.This also means that [I B ] is a simple diagonal matrix, as it is also defined in the secondary's body-fixed frame.
The generic form of the mutual gravitational potential between two bodies in the F2BP, U, is described using a double integral taken over the mass distributions of both the primary and secondary (Davis & Scheeres 2020), where d is the distance between a given finite mass element on the primary (dm A ) and a given finite mass element on the secondary (dm B ).In our simulations, we model U as a seconddegree approximation of the gravitational potential about an ellipsoidal secondary using MacCullagh's formula.This approximation utilizes a body's principal moments of inertia to describe the potential around a deformed body at any given distance from its center of mass (Murray & Dermott 1999).For a triaxial ellipsoid, it is expressed as where A B , B B , and C B represent the minimum, intermediate, and maximum principal moments of inertia, respectively, and  is defined as with (x, y, z) denoting the location of the primary in the secondary's body-fixed frame.This simplified sphere-ellipsoid model has also been used and validated against higher-fidelity models in Meyer et al. (2023b).Ultimately, it allows for efficient computation of the full range of possible 3D attitude and orbital motion while still capturing the short-term librational dynamics of the secondary.The attitude of the secondary's body-fixed frame relative to an inertial frame is tracked using a set of Euler parameters, with the attitude vector defined as , , , ].This attitude set was chosen in order to avoid singularities in the attitude during numerical integration.However, it will be useful to analyze the attitude of the secondary with respect to the rotating Hill frame of the orbit.Thus, the instantaneous rotating frame is built at each time frame of the numerical simulations, allowing for a direction cosine matrix (DCM) between the body frame and rotating frame to be constructed.A set of 1-2-3 Euler angles (θ 1 , θ 2 , θ 3 ) can be extracted from this DCM to evaluate the secondary's attitude relative to an initial singly synchronous configuration.Selecting the 1-2-3 Euler angle set specifically allows for intuitive understanding of the attitude, as (θ 1 , θ 2 , θ 3 ) represent roll, pitch, and yaw, respectively.Full justification for this choice of angles is { ˆˆˆ} is the body-fixed frame of the secondary, x y z , , { ˆˆˆ} is the body-fixed frame of the primary, x y z , ,

I I I
{ ˆˆˆ} is an arbitrary inertial frame, and r h , , q {ˆˆˆ} is the rotating Hill frame.
described in Agrusa et al. (2021) and Meyer & Scheeres (2021), but θ 2 and θ 3 can be thought of as the out-of-plane and planar libration angles of the secondary when θ 1 = 0. Detailed derivations of this process are omitted here, and general details regarding DCM construction can be found in Schaub & Junkins (2018).

SRP
The BYORP effect can be incorporated into the F2BP equations of motion through SRP forces.The location of the Sun in the secondary's body-fixed frame needs to be tracked throughout the simulation in order to compute the SRP force.We define the inertial position of the Sun relative to the secondary using the heliocentric orbital elements of the binary system and use a simple Kepler solver to compute the current location of the binary system along the heliocentric orbit.A DCM to transform from the inertial to the body-fixed frame, [BN], is constructed at each time step using the Euler parameters of the secondary: This DCM is then used to transform the inertial solar position vector into the body-fixed frame.The solar latitude and longitude can then be backed out from the definition of the Sun's position vector in the body-fixed frame in Equation (5): The Fourier force coefficients, A n and B n , are a function of solar latitude.Thus, tracking the solar latitude throughout the simulation allows us to extract the correct coefficients corresponding to the current solar latitude for a given day.
The SRP force can then be computed using Equation (1).SRP torques are included in the equations of motion following the same process.The Fourier moment coefficients, C n and D n , are computed at every solar latitude, and the SRP moment is calculated with Equation (9).

Tidal Torque
We incorporate tides into our simulations through the classical tidal torque model in Murray & Dermott (1999).This tidal model is driven by the tidal parameters k 2 and Q and has been used in the literature (Jacobson & Scheeres 2011b;Meyer et al. 2023b).This tidal torque model for a point source j acting on a spherical body i is expressed as where ω is the angular velocity of the body, ω orb is the angular velocity of the orbit, ρ is the density, M is the mass, R is the reference radius of the body, and r is the separation distance between the two bodies.However, this planar model is not sufficient for our full analysis of 3D motion.Meyer et al. (2023b) derives a formulation to account for full 3D motion, where is the relative spin rate of a body and rˆis the orbital position unit vector.This G ˆterm replaces the −sign(ω − ω orb ) quantity in the classical tidal torque model in Equation (24).Meyer et al. (2023b) also derives an effective tidal acceleration term to account for the total torque on the orbit.The reader is referred to Meyer et al. (2023b) for the full derivation of this term, but the acceleration term due to the torque on the orbit is where Γ orb is equal and opposite to the sum of the torques felt on both bodies: Γ orb = −(Γ A + Γ B ).The tidal acceleration term γ orb is derived from conservation of angular momentum, so it must be computed in an inertial frame and rotated back into the secondary's body-fixed frame for use in the F2BP equations of motion.
Tidal dissipation acts on both bodies of a binary system.Since primaries in these near-Earth systems are rapid rotators near their spin breakup limit, the primary is always rotating faster than the mean motion of the orbit.Thus, the tidal bulge leads the primary-secondary line of centers and will cause spin-down of the primary and expansion of the mutual orbit (Murray & Dermott 1999).However, the librations experienced by the secondary cause the tidal bulge to oscillate between leading and lagging the tide-raising primary.Since the tidal bulge is a physical phenomenon that results from gravity gradients and nonelastic deformation of the bodies, it is nonphysical for the tidal torque to switch signs instantaneously.Thus, we linearize the tidal torque on the secondary about 0 f »  to allow for a finite crossing time of the tidal bulge.The details regarding this linearization are found in Appendix C of Jacobson & Scheeres (2011b), and this linearization is also utilized in Meyer et al. (2023b).

Solar Third-body Effects
The last effect included in our numerical simulations are third-body perturbations from the Sun.Solar third-body effects are important to include in our simulations, as they have been shown to cause the solar node (i.e., the R.A. of the ascending node of the Sun in the secondary's body-fixed frame) to circulate and possibly trigger solar resonances (McMahon & Scheeres 2010a, 2010b;Ćuk & Nesvorný 2010).Full derivation of the acceleration term for this perturbation will be omitted, as it follows the general approach for deriving additional gravitational perturbations from a simplified twobody model.
where r Sun,s is the position vector pointing from the Sun to the secondary, r Sun,p is the position vector pointing from the Sun to the primary, and μ Sun is the gravitational parameter of the Sun.
Note that this acceleration term must also be computed in an inertial frame and then transformed into the body-fixed frame.

Equations of Motion
The final equations of motion used for numerical integration include the mutual gravity, SRP forces and torques, tidal torques, and solar third-body effects outlined in the previous subsections.We add these effects to Equations (15) and (16) to get the following equations of motion: The decoupled equations of motion for the spherical primary are simply The state vector for propagation is defined as X r r , , , , ] , with the corresponding state dynamics being

Methodology
In this study, we investigate the effects of shape changes due to both global deformation and cratering.There is still large uncertainty surrounding Dimorphos's key physical and material properties, so it will be impossible to know exactly how the shape of Dimorphos changed from the DART impact until Hera arrives.Additionally, the exact sensitivity of BYORP coefficients to changes in shape is not well understood (Scheirich et al. 2021).It is possible that small-scale topography could be influential in determining the BYORP coefficient.However, we do not account for local topography given limited observational data and focus solely on global shape changes.Thus, we explore a wide breadth of shape change possibilities to capture shape-change-induced effects on BYORP.In order to quantify the influence of shape on the BYORP effect, we look solely at the change in the dominant along-track 2 0  ( ) coefficient compared to the preimpact shape.From here on out, we will change the notation of the normalized 2 0  ( ) coefficient component to just B and refer to it simply as "the BYORP coefficient."Since we are primarily interested in the change in BYORP coefficient, we ignore any initial BYORP-induced secular drift caused by the actual shape of Dimorphos and use a base oblate spheroid with the extents of Dimorphos listed in Table 1 as a starter shape model.This smooth symmetric shape has a near-zero BYORP coefficient, allowing us to easily capture any ΔB resulting from shape changes.Modifications described in the following subsections are then made to this starter shape model.
Polyhedral shape models are created by first specifying a point cloud of vertices (accounting for any shape modifications) and then using Delaunay triangulation in MATLAB to generate a mesh of triangular facets.For all forms of reshaping considered in this work, we assume rigid bodies and instantaneous reshaping.Work in Raducan & Jutzi (2022) shows that momentum transfer from the DART impact can take up to several hours.However, this process is still shorter than the ∼11 hr orbit period, so instantaneity is a fair assumption.We assume an idealized impact in the sense that the impact direction is aligned with the surface normal at the point of impact.We also assume that the postimpact shape remains fixed for all time following the impact in our dynamical simulations.

Global Deformation Shape Model Generation
We use the methodology employed in Nakano et al. (2022) in order to generate globally deformed shape models.The full methodology can be found in Nakano et al. (2022) but will be summarized here for convenience.This method defines a bodyfixed coordinate system such that the x-, y-, and z-axes are aligned with the oblate spheroid's major, intermediate, and minor semiaxes, respectively.Six surface points are defined as the intersections between the surface of the body and these axes in both the positive and negative directions of the semiaxes.These surface points move locations in the body-fixed coordinate system when global deformation occurs.Reshaping is then quantified as the ratio of the new surface points of the deformed shape to the original surface points of the oblate spheroid.These reshaping parameters are s s s , for positive semiaxis , 32 where s = (x, y, z).
For BYORP coefficient analysis, we define the parameter space of global deformation by varying the a/b and b/c elongation ratios of the deformed body.Thus, assuming an idealized impact along the leading hemisphere, the result is a decreased intermediate axis and increased major axis.We model this by decreasing the +y-axis of the shape model while keeping the −y surface point constant, increasing the ±x-axes by equivalent amounts to capture a specified a/b ratio, and adjusting the ±z-axes such that the volume of the original oblate spheroid is preserved.An estimated (1.3-2.2) × 10 7 kg of ejecta was excavated from the DART impact, corresponding to a lower limit of 0.3%-0.5% of Dimorphos's mass depending on dust size and assuming a system bulk density of 2400 kg m This mass loss is negligible compared to the total mass of Dimorphos, so volume preservation is a fair assumption.Figure 3 shows an example of implementing this methodology of modeling global deformation.

Crater Shape Model Generation
For crater generation, we place the center vertex of the crater at the nominal estimated location of impact: 8.84°± 0.45°S, 264.30°± 0.47°E (Daly et al. 2023a).For BYORP coefficient analysis, we define the parameter space of cratering by varying the depth and diameter of the crater.All craters are either spherical or ellipsoidal in shape and have smooth surfaces.
Additionally, our code also allows for manipulation of the crater rim height relative to the surface normal of the center of the crater.To eliminate the effect of this parameter on BYORP coefficient analysis, we set the rim height to 1 m for all shape models.We assume the semiaxes of the original shape remain unchanged when a crater is added.An example of crater modeling is shown in Figure 4.

Computing BYORP Coefficients
BYORP coefficients are generated for each postimpact scenario using the facets for a given shape model.The methodology in Scheeres (2007) was used for this computation.The coefficients are computed for orders up to but not including n = 3, as it was shown in McMahon & Scheeres (2010b) that orders beyond this do not significantly alter the overall secular trends.The major assumption made for this computation is that the secondary body has zero thermal inertia.However, this assumption does not change the zerothorder coefficients relating to B. The resulting coefficients are in units of area, so they are normalized using Equations (11) and (12) and a volume-equivalent radius of ∼75.69 m (derived from the preimpact shape model) for the results presented in the following sections.

Effect of Shape Changes on BYORP Coefficient
We will assume that the Didymos system is in a tidal-BYORP equilibrium prior to the impact, since the system is likely near this equilibrium configuration (Scheirich & Pravec 2022).Therefore, we can use Equations ( 13) and ( 14) and the physical parameters (separation distance, volumeequivalent diameters of both bodies, and heliocentric orbital elements) provided in Table 1 to compute theoretical BYORP coefficients for the preimpact system.For this computation, we assume equal bulk densities of 2940 kg m −3 for Didymos and Dimorphos, which was computed based off of the 1.24 km separation distance.We use a single system bulk density value, as the bulk density of Dimorphos is essentially unknown.We derive the secondary's mass from the preimpact shape model and compute the primary's mass using the total system mass in Table 1.Since little is known about the tidal parameters of rubble piles, we use three values of k p /Q throughout our study to capture a range of plausible tidal strengths: 4 × 10 −6 , 1 × 10 −5 , and 2.5 × 10 −5 (listed from weakest to strongest).While arbitrarily selected, these values fall within the reasonable range of rubble-pile tidal parameters discussed in the literature (Jacobson & Scheeres 2011a;Nimmo & Matsuyama 2019).Preimpact BYORP coefficients corresponding to these tidal parameters are reported in Table 2.It is important to note that although we are assuming tidal-BYORP equilibrium for our simulations, there is uncertainty surrounding whether or not Didymos is exactly in this equilibrium state (based on the observed ΔM d ).Therefore, the preimpact BYORP coefficients could be slightly different than the ones we have computed.However, since we are also arbitrarily picking several tidal strengths to explore a range of preimpact B values, slight deviations from these values are insignificant compared to the range of B values we are already considering.Thus, choosing to not include the nominal mean anomaly drift in our simulations will not dramatically change the overall behavior of the results in this study.
The first set of results demonstrates how varying the reshaping and crater parameters changes the BYORP coefficient.For global deformation, we vary the a/b elongation ratio from 1.06 (its preimpact value) to 1.5, similar to Nakano et al. (2023b) and Richardson et al. (2023).The upper limit of 1.5 aligns with the lack of observed binary NEAs with elongations >1.5 (Pravec et al. 2016).The ΔB results for these reshaping scenarios are presented in Figure 5.The b/c elongation ratio varies from ∼0.7 to ∼2.1 to bound physically plausible results.Note that the results are bounded by realistic reshaping from the impact geometry (i.e., the intermediate axis must decrease, and the major axis must increase).Any result outside of these bounds would be nonphysical.The ΔB corresponding to the nominal postimpact elongation ratios of a/b = 1.3 and b/c = 1.1 (Table 1) is approximately −3.48 × 10 −2 .
For cratering, we explore diameters in the range of 10-60 m.Simulations performed by the DART Impact Modeling and Simulation Working Group show that depending on the nearsurface cohesive strength of Dimorphos, a defined crater as large as ∼30-60 m in diameter could be formed (Stickle et al. 2023).We vary the depth-to-diameter ratio from 0.1 to 0.5 to fully explore the crater parameter space.For reference, the depth-to-diameter ratio for the crater formed from Hayabusa2ʼs Small Carry-on Impactor on Ryugu was 0.15 ± 0.01 (Arakawa et al. 2020).Natural craters on Ryugu typically range from 0.14 to 0.20 (Sugita et al. 2019).The chosen ratios also align with the average depth-to-diameter ratio of 0.152 ± 0.04 for craters on Dimorphos observed by the DART spacecraft (Barnouin et al. 2023).The ΔB results for the crater cases are shown in Figure 6.
The ΔB results make sense conceptually.Force due to SRP is generated from thermal reemission of sunlight.Thus, a flat surface will be the most efficient at reflecting light.When the secondary is tidally locked, there is a fixed leading and trailing hemisphere of the body.In reality, Dimorphos is not a smooth body, so topography deviating away from a perfect ellipsoid is what causes the reradiated force to differ between the two hemispheres.The fixed hemisphere geometry relative to orbital motion and unequal SRP forces felt on both hemispheres from asymmetric topography is then what enables BYORP and   causes orbital evolution (Ćuk & Nesvorný 2010).Surface area also plays a role in this process, as the more surface area a given hemisphere has, the more reemission of SRP can occur.
For the preimpact oblate spheroid of Dimorphos, the leading hemisphere (particularly the near-equatorial region surrounding the leading point of the impact) has facet normals nearly aligned with the orbital velocity direction.Since the overall body shape is a rounded surface, the normals will not be perfectly aligned, but they will still all point in the same general direction as the velocity vector, and the facets will be effective at reflecting light when lit.The light on the leading hemisphere is reflected in the direction opposite of orbital motion; thus, the SRP force is acting to reduce orbital velocity.For the global deformation cases, we are essentially flattening the leading hemisphere.Compared to the original curved surface of the preimpact shape, the facet normals on the leading hemisphere of the collapsed shape models are now even more aligned with the velocity direction.Since this scenario causes the leading hemisphere to approach the flat surface end case, global deformation increases the effectiveness of light reflection on the leading hemisphere, causing the orbital velocity to be slower than before and the orbit to shrink.This means that the BYORP coefficient must be more negative than the preimpact coefficient, which is what we see in all elongation cases.A smaller b/c ratio for a given a/b corresponds to a shorter y-axis and thus flatter leading hemisphere.As established previously, a flatter surface increases a body's effectiveness at reflecting light and slowing down the orbit rate.Likewise, a larger a/b ratio for a given b/c correlates to a longer x-axis, allowing for more surface area along the leading hemisphere for a given intermediate axis.This larger surface area increases the body's effectiveness at reflecting light on the leading hemisphere, thus resulting in a more negative BYORP coefficient/slower orbital velocity.When looking at Figure 5, we can see that the magnitude in ΔB increases with increasing a/b and decreasing b/c, following this intuition.When a crater is added to the leading hemisphere, surface area that was originally contributing to slowing down the orbit is modified.We can see in Figure 6 that depending on the depth-to-diameter ratio, the effect of the crater on the BYORP coefficient differs.Note that we do not interpolate between crater shape models, as the changing shape models are not as parametric as the global deformation cases.This is likely due to placing the crater at the location of impact and not exactly on the center of the leading hemisphere.However, trends can still be observed from these results.The concavity of the crater causes some of the facet normals to no longer point in the direction near the orbital velocity, resulting in the facets being less effective and light being reflected elsewhere.For deeper craters, the cumulative SRP force generated on the leading hemisphere will be weaker and less efficient at opposing orbital motion, so the BYORP coefficient must be more positive than before.However, for shallower craters, especially with larger diameters, the concavity of the crater is less pronounced, and the surface starts to resemble a flat plate.This provides a similar effect as the global deformation cases, so the BYORP coefficient will become more negative.This shows that a change in B depends strongly on the crater morphology and impact geometry.We do not model morphologies like flatbottom craters, so if a crater does form from the impact, it is difficult to predict how exactly B will change until Hera arrives and characterizes the postimpact shape.
Depending on how much the BYORP coefficient changes due to the impact, the system could experience widely different dynamical results.If B changes but still remains negative, i.e., contractive, then the system should eventually settle at a new theoretical tidal-BYORP equilibrium semimajor axis according to Equation (13).However, if B becomes positive, then the system will experience joint expansive evolution, where both BYORP and tides are working to expand the orbit, and the semimajor axis will secularly expand (Jacobson & Scheeres 2011a).The exact effect on the dynamics of the system depends greatly on the system's tidal parameters k p /Q.In order for a system to be in tidal-BYORP equilibrium, contractive BYORP (and thus a negative BYORP coefficient) are required.If the ΔB caused by reshaping is negative, then B will only become more negative (relative to the nominal B values computed from tidal-BYORP dynamics and reported in Table 2).The orbit will then resettle at a new, closer or smaller, equilibrium semimajor axis, regardless of the tidal strength.This is the case for all global deformation scenarios considered   1).The black star represents the preimpact elongation ratios provided in Table 1.
here, as seen in Figure 7.For the nominal elongation values of a/b = 1.3 and b/c = 1.1, this equates to new equilibrium semimajor axes of 0.793, 0.896, and 1.00 km, corresponding to k p /Q = 4 × 10 −6 , 1 × 10 −5 , and 2.5 × 10 −5 , respectively.This means that the impact could lead to a total decrease in orbit size of approximately 36%, 28%, and 19% relative to the 1.24 km preimpact separation distance, respectively.On the other hand, if ΔB is positive, then this could have different implications for the system evolution depending on the tidal strengths.To demonstrate this, consider a single fixed value of ΔB.If tides are extremely weak, then a negative near-zero B is needed to keep the system in equilibrium.If the positive ΔB is such that it is larger than the magnitude of the preimpact equilibrium B, then this will cause the new B of the system to be positive.If tides are relatively strong (requiring a more negative B to counteract its effects), and the positive ΔB is now less than the magnitude of the original B, although increasing BYORP, the new B will remain negative compared to the weak tidal strength case experiencing the same ΔB.Therefore, the system should resettle at a new higher-equilibrium semimajor axis.Thus, contractive BYORP can turn into expansive BYORP either if tides are weak enough for ΔB to turn B positive or if ΔB is large enough to counteract stronger tides.This is seen in Figure 8, where the same set of ΔBs from crater shape changes is applied to three values of tidal strengths.For the first two tidal parameters (k p /Q = 4 × 10 −6 and 1 × 10 −5 ), we see that some of the ΔBs (typically corresponding to larger depth-todiameter ratios) are large enough relative to the tidal strength to turn the postimpact B positive and expansive.Since k p /Q = 1 × 10 −5 is stronger than k p /Q = 4 × 10 −6 , we see less instances of joint expansive evolution.For the strongest tides (k p /Q = 2.5 × 10 −5 ), we see that none of the ΔBs are large enough to cause expansive BYORP.

Numerical Dynamical Simulations
We simulate the long-term dynamics of the Didymos system the F2BP equations of motion presented in Equation (31).The full set of BYORP coefficients (A n and B n ) at every solar latitude up to n 2 for each shape model is set by adjusting the theoretical equilibrium B values in Table 2 to new values using the ΔBs of every coefficient from the reshaped model.We initialize each simulation with a 1.24 km separation distance between the two components and zero binary orbit inclination.We also account for obliquity between the binary mutual orbit plane and the heliocentric orbit plane and incorporate the heliocentric orbit such that the binary's rotation is retrograde relative to the heliocentric orbit (Richardson et al. 2022).We assume an idealized impact, such that the momentum from the impact is transferred solely in the plane of the mutual orbit and centered on the leading edge (i.e., no instantaneous change in torque).An along-track Δv T of −2.55 mm s −1 is applied to the initial conditions to account for the impact.This value corresponds to a separation distance of 1.24 km and system bulk density of 2940 kg m −3 .We initialize the simulations such that both spin poles of the bodies are  1).The black star represents the preimpact elongation ratios provided in Table 1.Note here that all cases are bounded because ΔB is always 0 in this analysis.aligned with the orbit normal.For each of the three tidal strengths per shape model, we assume the same tidal parameters for both the primary and the secondary.It is important to note that sources such as Jacobson & Scheeres (2011a) and Nimmo & Matsuyama (2019) predict that tidal strengths should scale with the size of the body.There is a lot of debate surrounding whether this scaling is linear or nonlinear in nature, but most sources suggest that Q/k scales with the radius of the body (Jacobson & Scheeres 2011a;Nimmo & Matsuyama 2019).Goldreich & Sari (2009) suggested that Q/k scales inversely with radius, but given the large uncertainty surrounding our knowledge of tidal parameters, we adopt the understanding that Q/k scales proportional to radius (Jacobson & Scheeres 2011a;Nimmo & Matsuyama 2019).This means that when using our convention of k/Q, tidal strength will increase with decreasing radius.Thus, using the same tidal parameters for both the primary and secondary provides us with a conservative prediction on the long-term evolution, as we would expect stronger tides on the secondary due to this scaling law to just damp out any excited libration at a faster rate.
The nominal initial conditions for the simulations are provided in Table 3.The only modifications made from Table 3 are the inertias of the globally deformed shape models and the initial angular velocity of Dimorphos to explore the sensitivity of the dynamics to initial conditions.The inertias were computed using a function in the General Use Binary Asteroid Simulator (GUBAS) that allows for computation of principal moments of inertia of a polyhedral shape model (Davis & Scheeres 2020).GUBAS is a publicly available tool (Davis & Scheeres 2021) based on the F2BP realization from Hou et al. (2017).
The nature of the dynamics was found to be highly dependent on the initial spin rate of Dimorphos.To explore the consequences of various assumptions, we test three different initial spin rates: a preimpact synchronous spin rate (denoted ω s ), a spin rate to conserve the angular momentum of Dimorphos due to the changing inertias from instantaneous reshaping (denoted ω c ), and a spin rate corresponding to a 20% loss in Dimorphos's angular momentum (denoted ω l ).The last option was selected as a test value to both investigate the sensitivity of the dynamics to initial conditions and account for angular momentum loss from ejecta.Since the dynamical and physical effects of Dimorphos's reshaping are not known, we test different spin rates to explore a breadth of possibilities for Hera's observations to compare to.Equations (34)-( 36) are used for computing ω s , ω c , and ω l , respectively: where A 0 , B 0 , and C 0 are the principal inertia values of the preimpact shape model and C reshaped is the inertia value of the postimpact shape model along its z-axis.Note that these initial spin rates are all magnitude quantities as we initialize the spin rate to be solely about Dimorphos's maximum inertia axis.The ω c and ω l spin rates and corresponding maximum inertia values expressed as percentages of their relative preimpact values for  all reshaped models are provided in Table 4.The postimpact elongation ratios reported in this table were selected by starting with the nominal a/b and b/c ratios listed in Table 1 and their respective upper and lower bound uncertainties and then exploring a grid of cases surrounding the predicted postimpact shape.
In actuality, the impact was slightly off-center from the leading edge of Dimorphos, which will change the initial spin rate.That was not accounted for in our simulations, though, as it is expected to be a small effect compared to excitation from the impact itself (Cheng et al. 2023).Investigating a range of spin rates for each shape will also capture the effect from this additional torque.Thus, there are nine simulations for a given modified postimpact shape model (determined by three tidal parameters and three initial spin rates).
The equations of motion are propagated for 10,000 yr following instantaneous reshaping.We use a variable-step, variable-order Runge-Kutta 8(7) integrator.The formulae are outlined in Prince & Dormand (1981), and this integrator is shown to be efficient for nonstiff initial value problems.This integrator was found to be sensitive enough to accurately capture both the long-term secular trends and short-term attitude of the dynamics.
Since the dynamics are non-Keplerian, it is misleading to present classical orbital elements of the mutual orbit, as they do not accurately represent the true dynamics (Meyer et al. 2023a).Thus, for all simulation results, we choose to present solely the separation distance, r, between the primary and secondary and the 1-2-3 Euler angle attitude of the secondary relative to the rotating Hill frame (θ 1 , θ 2 , θ 3 ).The separation distance intuitively is representative of the semimajor axis of the orbit, and the eccentricity of the orbit can be inferred by the magnitude of the variations in r.Larger variations over a given time period correspond to a more eccentric orbit, while smaller variations correspond to a less eccentric orbit.

Simulating Tidal-BYORP Equilibrium
To verify the theoretical equilibrium B values reported in Table 2, we simulate the three different tidal strength cases at equilibrium.Along with the dynamically relaxed orbit initial conditions, we also use the synchronous spin rate of the preimpact shape model, along with its corresponding moments of inertia.The resulting orbital separation distance and Euler angles are shown in Figure 9.
We can see that there is a slight variation in r, but this is small (on the order of ∼2 m), so the mutual orbit remains near zero eccentricity.Roll and pitch are virtually zero, as expected since no spin was initialized in those directions, and there is a small 0.2°planar libration.This is acceptable, as the preimpact synchronous spin rate about the minor axis is computed solely from the gravity model used and does not account for perturbations from the other forces and torques included in the simulations, so small excitation is expected.However, this is still very close to synchronicity.Although small, we also observe the expected outcome of tidal dissipation.We can see that libration and eccentricity are damped out from tides raised on the secondary, with faster damping corresponding to stronger tides.The most significant result from these plots is that tidal-BYORP equilibrium can be simulated even in the presence of additional dynamics.The theory described in Section 2.2 is based solely on BYORP and tides raised on the primary.We show that even when tides on the secondary, YORP torques, and solar third-body effects are included in the dynamics, the equilibrium state is stable and will hold for this system.Simulating this equilibrium is important, as the Didymos system was likely in or near this state prior to impact, so demonstrating we can simulate the equilibrium will allow us to show what happens when the system is knocked from this equilibrium due to the impact.

Simulating Effect of Craters
The dynamics of the crater cases were simulated assuming unchanged inertias and a preimpact synchronous spin rate.As seen in Figure 8, some of the larger, deeper craters can result in expansive BYORP.Figure 10 shows the resulting dynamical evolution for a 50 m wide, 25 m deep crater.For tidal strengths of k/Q = 4 × 10 −6 and 1 × 10 −5 , we see that the BYORP coefficient turns positive, and the mutual orbit should experience expansion.For the k/Q = 2.5 × 10 −5 case, BYORP is no longer expansive.B is still negative but less negative than before, so tides will dominate, and the system should expand and, in theory, settle at a new tidal-BYORP equilibrium semimajor axis of ∼1.4 km.This settling has not completed in the 10,000 yr time frame of this simulation.
A few other of the larger craters' dynamics were propagated, leading to similar results.Tweaking the initial conditions slightly (initial spin rate and inertias) can be enough in some crater cases to knock Dimorphos into a tumbling state.If a crater in fact formed on Dimorphos and elongation did not occur, then it is likely that only planar libration was excited significantly, with tumbling a possibility.Depending on the strength of the tides of Dimorphos, expansive BYORP is also possible, but it is unlikely unless there is a large crater with a larger depth-to-diameter ratio.

Simulating Effect of Global Deformation
The results for global deformation vary widely from those of the crater cases, primarily due to the changing inertias from reshaping.In total, 153 global deformation cases were simulated.In order to visually represent the outcomes of all cases, we have chosen to present the results in a manner similar to Ćuk et al. (2021).For each tidal parameter and spin rate, we plot b/c versus a/b and place a marker at each elongation ratio combination modeled in this study.The color of the marker indicates the initial outcome from the impact (where "initial" is defined as what Hera will observe 4-5 yr following the DART impact), while the shape of the marker represents the long-term behavior.Specifically, blue corresponds to excited planar librations only, orange corresponds to planar and nonplanar excited libration, magenta represents tumbling, and gray represents barrel instability.A square marker indicates eventual recapture to synchronicity, a star marker indicates prolonged chaos through the entirety of the simulation, and a circular marker indicates a case still experiencing significant librations by the end of the simulation (due to expansive BYORP and thus weaker tides/increasing libration).Cases that are marked as recapturing may still experience small librations (<30°) that have not been fully damped out at the conclusion of the simulation.We define tumbling as the roll-pitch-yaw Euler angles circulating about 180°, 90°, and 180°, respectively.Nonplanar libration is defined as bounded 3D oscillations, which differs from the full circulation of tumbling.We define these separately, as BYORP can still act on the system if librating out of plane but will not act on the system if full tumbling is experienced.Barrel instability is a phenomenon where the secondary is primarily aligned with the line of centers (potentially experiencing planar and nonplanar librations) but is slowly rolling along its long axis (Ćuk et al. 2021).While not quite tumbling, this effect will prevent BYORP from acting on the system.Note that for cases denoted as experiencing nonplanar librations and barrel instability from the impact, this does not occur immediately from the impact.We apply a Δv solely in the plane of the orbit, which will only excite planar librations initially.Out-of-plane motion is onset on the order of tens of days (typically <1 yr) following the impact.Therefore, nine different plots accounting for each combination of tidal strength and initial spin rate are presented in Figure 11.The nine plots are organized by increasing tidal strength from the leftmost column to the rightmost column.The top row of plots is for the preimpact synchronous spin rate assumption, the middle row is for the spin rate from conserved angular momentum for the changing inertias, and the bottom row is for the spin rate from the 20% loss in angular momentum.Once again, the postimpact elongation ratios in Figure 11 were chosen by exploring a grid of cases about the nominal postimpact shape with the a/b and b/c ratios listed in Table 1 and their respective uncertainties.
A few observations can be drawn from Figure 11.First, we will compare the effect of tides on the overall trends of the results.Looking from left to right, the results follow the expected trends of tidal dissipation.The weakest tides see the most prolonged attitude instabilities and in some cases are unable to tidally lock the secondary.As tidal strength increases, we see that fewer cases experience chaotic dynamics, and cases that experienced chaos with weaker tides now recapture with stronger tides.This is particularly evident for the middle row of results (spin rate derived from conservation of Dimorphos's angular momentum).Tidal strength does not appear to influence the nature of the dynamics in the direct aftermath of the impact.(This is seen by the fact that the colors remain the same from left to right.) The other observation that can be made from these plots is that both the initial attitude state and long-term behavior of Dimorphos are highly sensitive to the initial spin rate used in the simulations.Any spin rate different from the new postimpact orbital mean motion will induce libration with respect to the Hill frame.For the synchronous spin rate assumption, we see that only one shape experiences tumbling (a/b = 1.1, b/c = 1.4).However, if we move to the spin rate from conserved angular momentum, we see more instances of tumbling.If there is a loss in angular momentum of Dimorphos through events like ejecta loss or nonprincipal axis rotation, we see that for 20% angular momentum loss, all shape models result in tumbling, and most cases recapture to be synchronous.This is due to the fact that 20% loss is sufficient for all cases considered for Dimorphos to spin retrograde relative to the rotating frame, triggering tumbling.While we present the results for 20% rotational angular momentum loss, we also explored other percentages of this loss and found that a threshold of angular momentum loss exists where the spin rate passes through the preimpact synchronous rate.While in some cases this ranges from 5% to 10%, this is highly dependent on the magnitude of inertia changes.This emphasizes that until uncertainty surrounding the postimpact attitude and rotation rate of Dimorphos is resolved by Hera, it is difficult to make a confident prediction as to what is really happening with the short-term dynamics.
We also see that there is a strong dependence on the magnitude of reshaping for the resulting short-term dynamics.Elongation of the shape models (extent of the long axis within the equatorial plane) increases with a/b and typically increases with b/c (depending on whether the minor axis grew or shrank to conserve the total volume through the reshaping process).Looking at the row of synchronous spin rate cases, we see that our short-term dynamics results align closely with the results of Agrusa et al. (2021).In this study, the attitude unstable region for the short-term dynamics of Dimorphos is defined as a function of a/b and b/c for a triaxial ellipsoid for various momentum enhancement factors.Attitude instability is driven by resonances between fundamental frequencies of the system.Looking at the plots for momentum enhancement factors of 3 and 5 in Figure 9 of Agrusa et al. (2021; since the true value of β is likely within this range), we see that all of our cases experiencing barrel instability and tumbling lie within the attitude unstable region, and the cases experiencing excitation of nonplanar libration roughly border this unstable region.We also see more instances of barrel instability rather than pure tumbling, agreeing with the results in Agrusa et al. (2021) that Dimorphos is more prone to rolling about its long axis.Specifically, only the a/b = 1.1, b/c = 1.4 case experiences tumbling, as that is the only shape that falls within the unbounded region of roll and yaw.The cases experiencing barrel instability fall within the region of unbounded roll but a maximum yaw of about 30°-60°, agreeing closely with our results.While the region of instability is not a hard boundary, it is interesting to observe that our results also follow these trends.It is also interesting to note that while attitude instability is predicted in the short term, for the synchronous assumption, we observed that almost all attitude unstable cases eventually recapture a synchronous state.
Looking now to the middle row of results, we observed more cases of prolonged chaotic tumbling.An apparent trend here is that the shapes experiencing prolonged chaos are the shapes experiencing the highest degree of flattening on the leading hemisphere for a given a/b ratio, corresponding to more pronounced asymmetry about the long axis.This results in a higher likelihood of rolling about its long axis and triggering full tumbling/chaos.
Although Figure 11 shows the generalized behavior of each simulation, we observed a variety of results within these categories and choose to highlight a few representative simulations.We will show only cases falling within the uncertainty bounds of the postimpact elongation ratios listed in Table 1.All global deformation simulation results can be found on Zenodo using the following DOI:10.5281/zenodo.10475980(Cueva et al. 2024).Figure 12 shows the result of the nominal postimpact elongations from Table 1 (a/b = 1.3, b/c = 1.1) for the synchronous spin rate assumption.This is an example of the most "relaxed" simulations we observed, where only planar libration is excited (∼20°) and is damped out over time (within ∼500-3500 yr), with stronger tides tidally locking the secondary faster.Since B is now more negative than before due to global deformation, we see the anticipated result of inward migration from stronger BYORP that will eventually settle at a lower tidal-BYORP equilibrium semimajor axis.Figure 13 shows the same nominal shape but with the spin rate from enforcing conservation of angular momentum for the changing inertias.Just by tweaking  one parameter, the system now experiences tumbling from the impact and only recaptures for the two stronger tidal ratios.The weakest tides are unable to recapture the secondary into synchronous lock.Thus, tumbling can be long-lived, even for a tight binary system like Didymos.
We also observe instances of extremely excited dynamics for different shapes falling within the postimpact elongation ratio uncertainty bounds.For a/b = 1.3 and b/c = 1.2, for both the synchronous and conserved angular momentum spin rates, we see a quick onset of nonplanar libration within approximately 30 days following the impact and getting damped out on the order of 10 2 -10 3 yr.The case corresponding to the synchronous spin rate is shown in Figure 14.The case with the updated spin rate follows the same trend but with less libration amplitude.Note that contractive BYORP is still active when Dimorphos is experiencing nonplanar libration.
For the shape of a/b = 1.36 and b/c = 1.2, we observed barrel instability for the synchronous assumption, shown in Figure 15.For this instance, Dimorphos is librating both in and out of the orbit plane, but it is also rolling about its long axis.As we see in the plot of r, the barrel instability terminates BYORP.Once the barrel instability is damped out from tides (in about 10 2 -10 3 yr), BYORP reengages.This is an important result building off of the work from Ćuk et al. (2021), as we show that tidal damping can stop barrel instability.
For the spin rate from conserved angular momentum for this shape, shown in Figure 16, we see that Dimorphos is now tumbling.However, these results demonstrate how Dimorphos can recapture in a different orientation following a period of chaos.For the k/Q = 4 × 10 −6 and 1 × 10 −5 cases, Dimorphos recaptures in a 180°flip about its minor axis (major inertia axis).This behavior has been observed before in Ćuk & Nesvorný (2010).However, this flip has major implications for the long-term evolution, as it reverses the effect of BYORP.Recall that the expected BYORP evolution from global deformation is to shrink the orbit.Now that the body is flipped, BYORP is acting to expand the orbit along with tides.For the k/Q = 2.5 × 10 −5 case, instead of flipping about its minor axis, we see that Dimorphos has instead flipped 180°a bout its major axis (minor inertia axis).This is a new type of observed simulation result and as far as we are aware has not been discussed in the literature.This recaptured orientation has the same effect on BYORP as the 180°flip about the minor axis, as it also reverses the hemispheres and has the reshaped  side of Dimorphos now on the trailing hemisphere relative to the mutual orbit.Note that the orientation of Dimorphos about its major and minor inertia axes from recapture of chaotic motion is random and not possible to predict, agreeing with the conclusions drawn from Ćuk & Nesvorný (2010).
Some of the most interesting observed excited dynamics occur for a/b = 1.24.Figure 17 shows a case where initially only planar libration is excited, but a resonance between the binary orbit and heliocentric orbit (Touma & Wisdom 1998) at ∼1500 yr into the simulation excites eccentricity and more planar libration into the system.This is called the "evection" resonance, which can occur when apsidal precession of the binary mutual orbit is resonant with the mean motion of the heliocentric orbit (Ćuk & Nesvorný 2010).Tides are unable to combat the effect of the resonance and damp out the libration, so the resonances eventually lead to full attitude instability approximately 6500-9000 yr into the simulation.This shows that even if attitude instability is not a direct result of the DART impact, both unstable shapes and other phenomena like these semisecular resonances can drive the system into chaos.This reshaped model also demonstrates that it is possible for the system to settle at the new theoretical tidal-BYORP equilibrium value predicted in Figure 7.This is seen for the k/Q = 2.5 × 10 −5 case, as Dimorphos recaptures at a semimajor axis slightly below its new theoretical equilibrium value at ∼7000 yr, and tidal expansion slowly drives the system to its equilibrium value of ∼1 km.It is possible that the resonance experienced early on in the simulation sped up the system's evolution toward the new equilibrium semimajor axis.This is suggested by the fact that there is a large downward jump in r when Dimorphos recaptures from the chaotic period.However, further studies will have to be conducted to determine if this is a likely outcome from the resonances, or coincidence.
These semisecular resonances are also observed in quite a few of the other simulations.For example, in Figure 18, we see that a resonance at ∼1000 yr excites stronger planar librations and eccentricity for the k/Q = 4 × 10 −6 and 2.5 × 10 −5 cases and nonplanar librations for these cases at ∼4000 yr.Effects from the resonances can be less pronounced as well, leading to significantly smaller excitation of eccentricity and libration, as seen in the k/Q = 1 × 10 −5 case for this shape model.Figure 19 shows a particularly interesting result from the onset of semisecular resonances.After triggering tumbling, Dimorphos recaptures in an orientation where it has flipped 180°about both its major and minor axes.We saw that if flipping about one of these axes occurs, then it will reverse the effect of BYORP.However, in this case, the ends of the long axis have switched locations, but the double flipping keeps the flat reshaped side of Dimorphos in the same location on the leading hemisphere.Thus, BYORP will have the same effect on the orbit and drive inward migration.This orientation change could have interesting implications for other reshaping processes following the initial reshaping from the impact.

Observability of BYORP during Hera
Hera will arrive at Dimorphos in late 2026, with observations beginning in 2027 (Michel et al. 2022).As seen throughout Section 5, different attitude states of Dimorphos have various implications for the BYORP-induced evolution of Dimorphos.Planar and/or nonplanar librations about the 1:1 spin-orbit resonance will still allow for BYORP to act on the system.Attitude instability, such as tumbling or barrel instability, will terminate BYORP.Assuming BYORP is not acting on the system, semimajor axis evolution due solely to tides will be on order of a few centimeters in 5 yr (Richardson et al. 2022).To get an idea of how much secular drift Dimorphos's orbit can experience by the time Hera arrives if both tides and BYORP are acting on the system (i.e., assuming only libration has been excited), we have computed approximate values for the nominal postimpact elongation ratios and their uncertainties listed in Table 1.We observe secular separation distance drift rates ranging from −1.8 to −4.6 cm yr −1 in our simulations, implying an inward drift of ∼9−23 cm in the 5 yr period between DART and Hera.Using Equation (37) from Scheirich & Pravec (2022), we can compute mean motion rates (n ) and mean anomaly drifts (ΔM d ) corresponding to the drift rates pulled from our simulations: We get n  values ranging from 5.5 × 10 −17 to 1.4 × 10 −16 rad s −2 , corresponding to ΔM d ranges of 1.6°-4.0°yr−2 .Due to the orbitattitude coupled nature of this system, variations in the orbit's  separation distance will occur if Dimorphos's attitude is dynamically excited.Based on our predicted secular drift values, it will be challenging to differentiate secular drift from these variations in such a short time span.

BYORP Coefficient of Current Dimorphos Shape Model
While the DART mission was able to provide an extremely high-resolution shape model of Dimorphos preimpact, it is unfortunately only one side of the body.If the back side of Dimorphos was imaged and a full shape model was constructed, then we could compute B for that shape and use our current understanding of tidal-BYORP equilibrium to constrain the strength of tides on Dimorphos.However, we can still compute a BYORP coefficient for the current shape model presented in Daly et al. (2023b) and provided in Table 1.We computed the BYORP coefficient for the shape model generated from a ground sample distance of 972 mm (comprised of 99,846 vertices and 196,608 facets) and found the normalized B value to be −1.069× 10 −2 .This shape model consists of the "front half" topography observed by DRACO stitched onto a smooth oblate spheroid.Thus, the BYORP coefficient we report is for this shape model with a rough front half and smooth back half.By knowing how much the front half of Dimorphos is contributing to BYORP, we can then determine how much the back side of Dimorphos should be contributing depending on the tidal strength, since the smooth back half of the shape model is currently not contributing to BYORP.Future work could use this information to constrain the geography of the back side of Dimorphos.
Hera should be able to provide a full detailed shape model of postimpact Dimorphos.Imaging the back side of Dimorphos, even postimpact, would help better constrain the full preimpact model.Even if elongation and/or surface motion occurred in the time period between DART and Hera, this would at least allow for some estimate of tidal parameters to be obtained from the preimpact equilibrium.Additionally, the full postimpact shape, along with the known attitude state of Dimorphos, will allow us to generate a B value for Dimorphos and fully test the current BYORP theory against how the system is actually evolving.This could also be interesting to help figure out how much material was deposited on the back side of Dimorphos and changed the shape.

Limitations of Current Study
In our study, we assume instantaneous reshaping and a fixed shape of Dimorphos for the entire length of the simulation.In actuality, Dimorphos is likely to experience continued reshaping from other processes throughout its dynamical evolution.For example, surface motion could be triggered from local slope failure caused by the perturbed rotation state and/or tidal acceleration of Dimorphos (Agrusa et al. 2022).The impact also excavated a large amount of ejecta from Dimorphos, approximately (1.3-2.2) × 10 7 kg (Graykowski et al. 2023).The fate of the ejecta could affect the shape of Dimorphos and the dynamics of the system, depending on whether the ejecta redistributes back on Dimorphos, transfers to Didymos, or escapes the system.
Additionally, we assumed an idealized impact such that the momentum is transferred through Dimorphos's center of mass in the plane of the orbit.In reality, the impact imparted a small amount of nonplanar Δv and was off-center from the center of the leading hemisphere by approximately 25 m (Daly et al. 2023a).While not considered in this study and shown to be a relatively small effect (Richardson et al. 2022;Cheng et al. 2023), work in Richardson et al. (2023) shows that the additional torque from the impact introduces a larger degree of libration and attitude instability throughout the a/b and b/c parameter space.
We also assumed that both Didymos and Dimorphos have the same tidal parameters.However, this is not necessarily true, and the rates of tidal dissipation could differ in both bodies, affecting the dynamics.Work performed in Efroimsky (2015) also suggests that tidal parameters are frequency-dependent, whereas we have assumed a constant tidal frequency.It is also possible that tidal dissipation for rubble piles could be driven by viscosity (Efroimsky 2015) rather than the friction-based approach (Goldreich & Sari 2009;Nimmo & Matsuyama 2019) used in this work.
Additionally, our dynamical model does not capture all shape-based gravity effects like GUBAS, so it is possible that some details of the short-term dynamics are inaccurate here; however, the secular trends should be quite similar.Future studies are required in order to investigate this further.

Extensions of Current Study
The methods used in our work can extend beyond this application to Didymos and DART.Other binary asteroid systems may undergo natural processing that can cause reshaping.This includes (but is not limited to) small impacts, planetary flybys, tidal interactions, and thermal cracking or other active asteroid processes (Richardson et al. 1998(Richardson et al. , 2004;;Binzel et al. 2010;Ballouz et al. 2019;Agrusa et al. 2022).In those cases, our methods could be used to predict how other bodies will dynamically evolve.

Conclusions
In this work, we explore the range of possible changes to BYORP due to reshaping of Dimorphos from the DART impact.We consider shape changes from both cratering and global deformation.We leverage the F2BP to propagate the orbit-attitude coupled dynamics of Didymos, evaluating the system's long-term tidal-BYORP evolution.We saw that depending on the initial conditions of the system, degree of reshaping, and tidal parameters, we observe the full range of possible dynamical results.This study has bounded predictions of the long-term secular evolutionary effects of the Didymos system post-DART, which will be later improved upon once uncertainties surrounding the attitude and reshaping of Dimorphos are reduced by Hera.There is not enough time for tides to recapture Dimorphos into the 1:1 spin-orbit resonance by the time Hera arrives (if Dimorphos is in fact dynamically excited), so whichever excited state Dimorphos entered into following the impact will persist throughout Hera's observations.Additionally, detection of BYORP from Hera is unlikely but not impossible.If only libration is excited, BYORP will continue to act, but it will only affect the orbit by a small amount in the few-year gap between DART and Hera.If attitude instability from tumbling or barrel instability is triggered from the impact, then BYORP will be terminated.Tides will still continue to act on the system, possibly allowing for measurement of k/Q if changes in the semimajor axis are observed.Ultimately, Hera will help to improve our knowledge of the state of Dimorphos following the DART impact, allowing us to better refine our predictions for the fate of the secular dynamical evolution of the Didymos system.
A summary of the main conclusions drawn from this work is presented as follows.
1.The formation of a crater results in a net positive ΔB for most cases because the impact location was on the leading hemisphere of Dimorphos, with large/deep crater morphologies leading to a net negative ΔB.All global deformation cases resulted in a net negative ΔB. 2. The preimpact BYORP coefficient of Dimorphos is negative, i.e., causing contractive BYORP.A negative ΔB will drive the system to a new, lower, tidal-BYORP equilibrium semimajor axis.A positive ΔB can drive the system to a new, higher, equilibrium semimajor axis if the overall B remains negative.If ΔB is large enough to make B positive, this will lead to expansive BYORP and an unbounded theoretical semimajor axis.3.If Dimorphos is elongated postimpact, the system evolution can experience the full range of dynamical possibilities depending on the elongation ratios.This includes planar libration, nonplanar libration, barrel instability, and tumbling depending on the particular axis ratios.4. The Didymos system will most likely recapture in a singly synchronous configuration, but it is possible for tumbling or barrel instability to be long-lived.If tides are strong enough, they can damp out these phases.5. Semisecular resonances between the mutual orbit and heliocentric orbit can excite libration and eccentricity in the system, sometimes making the system chaotic.6. Dimorphos can randomly recapture in a different orientation following a period of tumbling.If the body flips 180°about either its minor axis or its major axis, the nature of the BYORP evolution will be reversed.If the body flips 180°about both of these axes, the BYORPinduced evolution will not be affected.

Figure 1 .
Figure 1.Position of the Sun in the secondary's body-fixed frame, described by the vector u ˆ, solar latitude (δ s ), and solar longitude (λ s ).Note that we assume the body-fixed frame x y z , , s s s{ ˆˆˆ} is aligned with the principal moment of inertia axes of the secondary for this work.

Figure 2 .
Figure 2. Depiction of the simplified sphere-ellipsoid dynamical model and relevant frames.x y z , , s s s −3 (Graykowski et al. 2023).Simulations in Raducan et al. (2023) suggest that up to 1% of Dimorphos's mass could have been ejected, and observations in Roth et al. (2023) indicate 0.2%-1.2% of Dimorphos's total mass was ejected.

Figure 4 .
Figure 4. Example of a shape model with a crater at the estimated impact site location.The crater has a diameter of 30 m and depth of 5 m, with a rim height of 1 m.

Figure 6 .
Figure 6.Changes in BYORP coefficient, B, as a function of crater depth and diameter.The white/black crosshairs mark the specific depth and diameter values corresponding to each case.The colored boxes surrounding the crosshairs indicate the magnitude of ΔB.

Figure 5 .
Figure 5. Changes in BYORP coefficient, B, as a function of postimpact elongation ratios a/b and b/c.The solid and dashed black lines represent the nominal postimpact elongation ratios and uncertainty, respectively (Table1).The black star represents the preimpact elongation ratios provided in Table1.

Figure 7 .
Figure 7. New theoretical tidal-BYORP equilibrium semimajor axes for global deformation ΔB values for varying tidal parameters.The solid and dashed gray lines represent the nominal postimpact elongation ratios and uncertainty, respectively (Table1).The black star represents the preimpact elongation ratios provided in Table1.Note here that all cases are bounded because ΔB is always 0 in this analysis.

Figure 8 .
Figure 8. New theoretical tidal-BYORP equilibrium semimajor axes for crater ΔB values for varying tidal parameters.The white/black crosshairs mark the specific depth and diameter values corresponding to each case.The colored boxes surrounding the crosshairs indicate the magnitude of ΔB. White boxes correspond to joint expansive cases where the semimajor axis is theoretically unbounded.

Figure 9 .
Figure 9. Simulation of tidal-BYORP equilibrium for three different tidal strengths.The separation distance is shown on the left, and the attitude of the secondary is shown on the right.

Figure 10 .
Figure 10.Simulation of dynamical evolution for Dimorphos with a 50 m wide, 25 m deep crater and preimpact synchronous spin rate and inertias.The separation distance is shown on the left, and the attitude of the secondary is shown on the right.

Figure 11 .
Figure 11.Overview plots showing dynamical evolution as a function of elongation ratios a/b and b/c.

Figure 12 .
Figure 12.Simulation of dynamical evolution for reshaped model corresponding to a/b = 1.3 and b/c = 1.1 and a preimpact synchronous spin rate.The separation distance is shown on the left, and the attitude of the secondary is shown on the right.

Figure 13 .
Figure 13.Simulation of dynamical evolution for a reshaped model corresponding to a/b = 1.3 and b/c = 1.1 and a spin rate to conserve the angular momentum of Dimorphos for the changing inertias.The separation distance is shown on the left, and the attitude of the secondary is shown on the right.

Figure 14 .
Figure 14.Simulation of dynamical evolution for a reshaped model corresponding to a/b = 1.3 and b/c = 1.2 and a preimpact synchronous spin rate.The separation distance is shown on the left, and the attitude of the secondary is shown on the right.

Figure 15 .
Figure 15.Simulation of dynamical evolution for a reshaped model corresponding to a/b = 1.36 and b/c = 1.2 and a preimpact synchronous spin rate.The separation distance is shown on the left, and the attitude of the secondary is shown on the right.

Figure 17 .
Figure 17.Simulation of dynamical evolution for a reshaped model corresponding to a/b = 1.24 and b/c = 1.1 and a preimpact synchronous spin rate.The separation distance is shown on the left, and the attitude of the secondary is shown on the right.

Figure 16 .
Figure 16.Simulation of dynamical evolution for a reshaped model corresponding to a/b = 1.36 and b/c = 1.2 and a spin rate to conserve the angular momentum of Dimorphos for the changing inertias.The separation distance is shown on the left, and the attitude of the secondary is shown on the right.

Figure 18 .
Figure 18.Simulation of dynamical evolution for a reshaped model corresponding to a/b = 1.24 and b/c = 1.2 and a preimpact synchronous spin rate.The separation distance is shown on the left, and the attitude of the secondary is shown on the right.

Figure 19 .
Figure 19.Simulation of dynamical evolution for a reshaped model corresponding to a/b = 1.24 and b/c = 1.2 and a spin rate to conserve the angular momentum of Dimorphos for the changing inertias for k/Q = 2.5 × 10 −5 only.The separation distance is shown on the left, and the attitude of the secondary is shown on the right.

Table 1
Relevant Physical and Dynamical Parameters of the Didymos System a a Updated with the latest values at the time of this work, taken from the DART mission internal document (DRA v. 5.01).The most updated values are listed in The final equation for this term is

Table 2
Preimpact B Values for Given Tidal Parameters

Table 3
Nominal Initial Conditions of Dynamical Evolution Simulations

Table 4
Updated Inertias and Spin Rates for Globally Deformed Shape Models