The Yarkovsky effect in REBOUNDx

To more thoroughly study the effects of radiative forces on the orbits of small, astronomical bodies, we introduce the Yarkovsky effect into REBOUNDx, an extensional library for the N-body integrator REBOUND. Two different versions of the Yarkovsky effect (the"Full Version"and the"Simple Version") are available for use, depending on the needs of the user. We provide demonstrations for both versions of the effect and compare their computational efficiency with another previously implemented radiative force. In addition, we show how this effect can be used in tandem with other features in REBOUNDx by simulating the orbits of asteroids during the asymptotic giant branch phase of a 2 $M_{\odot}$ star. This effect is made freely available for use with the latest release of REBOUNDx.


INTRODUCTION
Radiation pressure (Nichols & Hull 1903), Poynting-Robertson drag (Poynting 1904;Robertson 1937), the Yarkovsky effect (Radzievskii 1954;Peterson 1976), and the YORP (Yarkovsky-O'Keefe-Radzievskii-Paddack) effect (Radzievskii 1954;Paddack 1969;O'Keefe 1976) are four radiative forces that can be important for the orbital evolution of small bodies. Both radiation pressure and Poynting-Robertson drag (PR drag) arise from the absorption and scattering of radiation by dust and other particles (Burns et al. 1979). Radiation pressure is generally a force that propels particles away from their host stars, while PR drag is a force that causes particles to lose momentum and spiral inwards (Veras et al. 2015).
Unlike radiation pressure and PR drag, the Yarkovsky and YORP effects appear only when considering objects' angular momentum (Vokrouhlický et al. 2015). Both phenomena occur after a radiation-induced temperature gradient forms on the surface of a body, causing an uneven emission of radiation between its hemispheres. As an object rotates relative to the substellar point, this imbalance in the emission creates a net force on the object. While both effects are based on this principle, the Yarkovsky effect changes the orbit of a body, while the YORP effect creates a torque that changes the body's spin. In addition, the YORP effect only operates on nonspherical, asymmetric objects while the Yarkovsky effect can operate on spherical ones. Even with the differences between them, all of these effects can have important consequences on the dynamical evolution of small bodies; this motivates developing tools to numerically model these phenomena.
is an open-source, N-body integrator for gravitational dynamics (Rein & Liu 2012). Flexible and efficient, can be used in C and is also available as a Python package.
is an extensional library to that contains additional effects and forces (general relativity, tidal dissipation, etc.) that users can add to simulations (Tamayo et al. 2020). This library includes radiative forces (Tamayo et al. 2020), but contained only radiation pressure and PR drag and was lacking both the Yarkovsky and YORP effects. This paper introduces the Yarkovsky effect as another radiative force available in and shows its capabilities and potential for future work.
We present two separate versions of the Yarkovsky effect available for use in : the Full Version and the Simple Version. While both of these versions are based on the Yarkovsky effect model derived in Veras et al. (2015), the Simple Version includes modifications done by Veras et al. (2019) that compromise long-term accuracy for a reduction in computational time and the number of necessary parameters. More detail will be given in later sections on the differences between these two versions and when to appropriately use them. Please note that rotational perturbations created by the YORP effect are not included within this module. YORP requires detailed models of the surface of a body (Nesvorný & Vokrouhlický 2007;Scheeres & Mirrahimi 2008) that are beyond the scope of this paper. For now, we assume that all bodies are spherical and have a uniform density. Adding the YORP effect into will be left for future work.
The Yarkovsky effect has previously been added to other Nbody integrators, including (Chambers 1999) and (Duncan & Levison 1997). While these integrators have their own advantages, contains features that make it better suited for a wide variety of studies. contains a larger array of different integrating schemes, providing nine different options for numerical integrators while and only provide four and five, respectively. In addition, comes standard with features such as collision tracking and the ability to create simulation archives that save particle and orbital parameters at certain times during an ongoing simulation. Features like these would have to be manually added to these other codes by users. will also always yield results that are identical across all platforms and compilers. Finally, gives users greater flexibility when running simulations in , allowing them to easily add multiple effects at once to simulations (see Section 5).
Section 2 describes the main equations that govern the Yarkovsky effect in . In Section 3, we introduce the Full Version of the effect and demonstrate one of its possible applications, and in Section 4, we introduce the Simple Version and show how an analytic equation describing the evolution of an object's semi-major axis can be derived from it under specific circumstances. Finally, we show how the Yarkovsky effect can be combined with other previously implemented effects and test how they perform together in Section 5.

THE YARKOVSKY EFFECT IN REBOUNDX
Equation (27) from Veras et al. (2015) provides the model for calculating the acceleration on a body due to the Yarkovsky effect. It is implemented by both versions of the Yarkovsky effect in .
where is the luminosity of the object's host star and is the speed of light. In addition, , , and are the bond albedo, radius, and density of the body, and is its distance from the star. The 3 × 3 matrix Y is the rotation matrix that contains the physics of the Yarkovsky effect, and i is the relativistically-corrected direction of the incoming radiation. The parameter is a constant between 0 and 0.25; it was created to avoid the need to model the complex spin behavior of most asteroid-like objects. Realistically, its value will change over time, but in these simulations, it remains constant to reduce complexity. If the target's rotation speed approaches the critical rotation speed at which the target will begin to break up (Walsh et al. 2008;Vokrouhlický et al. 2015), then → 0. If the target's rotation period approaches its orbital period, then → 0.25.
Equation (1) includes the physics for both the diurnal and seasonal components of the Yarkovsky effect. Both of these contributions are based on the same mechanism: the heated side of a body will release photons that are both more abundant and more energetic than those from the dark side. This leads to a net force on the body opposite the emission direction (Bottke et al. 2002). The changes in the orbital parameters of an object depend on how this heated side is oriented when this radiation is emitted. The diurnal contribution comes from the temperature differences created in the body due to its rotation with respect to its orbit. This component can cause an increase or decrease in the semi-major axis of a body, depending on the orientation of its spin axis. The seasonal contribution comes from temperature differences created due to the specific angular momentum of an object in its orbit. This component will always cause a decrease in the semi-major axis of the body and does not depend on the orientation or magnitude of the object's spin. Figure 1 from Bottke et al. (2002) provides a visual diagram of these two contributions, and Section 3.1 provides more information on how they compare and depend on different model parameters. Both of these components are required for a complete model of the Yarkovsky effect.
We will now discuss the equations that comprise both Y and i. The vector i is described by the following: where r is the position vector of the body and v is its velocity vector. If ||v|| << then i ≈ r/ , which is the unit vector in the direction of the incoming radiation. The rotation matrix Y is the product of two matrices: where R (s, ) is the contribution from the diurnal component of the Yarkovsky effect and R (h, ) is the contribution can be constant or time-dependent. Section 5.1 gives an example of how the luminosity can be configured to change over time in a simulation. from the seasonal component of the effect. The vector s is the spin axis of the body, is the thermal lag angle along the equator of the body, h is the object's specific angular momentum axis, and is the thermal lag angle in the orbital plane of the body. The equations for and are based off the specific 1D conduction model found in Brož (2006): where is the Stefan-Boltzmann constant, is the object's emissivity, Σ is its rotational period, Π is its orbital period, and Γ is its thermal inertia. Both and only have meaningful values between 0 • and 45 • . The full equations and derivations for R (s, ) and R (h, ) can be found in Veras et al. (2015). Section 3 and Section 4 will describe how the two versions of the Yarkovsky effect are implemented in .

FULL VERSION
The Full Version of the Yarkovsky effect in is directly based on the model for the Yarkovsky effect found in Veras et al. (2015). Setting a parameter called "ye_flag" to 0 on bodies in simulations will apply this version of the effect to these chosen targets. Users must input all of the necessary parameters contained in the model for this effect to work. As is standard for all effects in -, the parameters must be inputted with the same units as the simulation where the effects are being added. The large number of parameters included in this version of the effect provides a high level of detail for the perturbations created by the Yarkovsky effect on a particular object. However, this means long simulations containing many bodies with the Full Version active will require a large amount of computational time (see Section 5.2).

Demonstration
To demonstrate the Full Version of the effect, we perform a simple parameter study that explores how the obliquity, radius, and thermal inertia of an asteroid affect the strength of the Yarkovsky effect. This is similar to a study Carruba et al. (2017) performed to create Figure 15 in their study of the Veritas asteroid family.
We create a 0.01 Myr-long simulation that contains the Sun and an asteroid in a circular, non-inclined orbit at a distance of 3.165802 au, the semi-major axis of 7231 Veritas (Carruba et al. 2017). This value was chosen to place the asteroid at a realistic distance within the Veritas asteroid family. Table 1 lists the unchanging physical properties of the asteroid. The  F . Figure 1 shows how the change in semi-major axis of these simulated asteroids depends on these changing parameters. As expected from Equation (1), the change in semi-major axis increases as the radius of the body decreases. In addition, Equations (4) and (5) describe how bodies with small thermal inertias will radiate away too much heat before they rotate enough for the effect to activate, leading to a nearly negligible change in semi-major axis from the Yarkovsky effect. Figure 1 also shows how the strength of the diurnal and seasonal components of the Yarkovsky effect depend on an object's physical properties. For example, the lines in Figure 1 begin to noticeably trend downwards as thermal inertia reaches values above 1000 . Objects with extremely large thermal inertias will hold onto thermal energy for longer periods of time as they move throughout their orbits, leading to an increased contribution from the seasonal part of the Yarkovsky effect when compared to the diurnal portion. This, combined with the fact that the seasonal component always leads to a decrease in the body's semi-major axis, explains the increased negative change in semi-major axis for bodies with higher thermal inertias.
In addition to thermal inertia, obliquity plays a large role in determining the relative strengths of the seasonal and diurnal components of the Yarkovsky effect. Bodies with obliquities at 0 • and 180 • will feel the strongest perturbations in their semi-major axis from the diurnal component. As the obliquity moves towards 90 • , the diurnal portion's effect on the object's semi-major axis will begin to decrease and its effect on the object's inclination will start to increase. At an obliquity of 90 • , only the seasonal component will affect the body's semi-major axis. This is why objects at 90 • in Figure 1 remain relatively unperturbed until their thermal inertias reach a high enough value where the seasonal component becomes much more significant. In between these extreme values, the body will experience the combined effects of the diurnal and seasonal components. However, because the object's orbital mean motion is small compared to its rotational frequency, the diurnal portion will dominate for most obliquity and thermal inertia values and be the largest determining factor for the object's change in semi-major axis. This is why, in Figure 1, the object's overall change in semi-major axis is greater for obliquity values that are closer to 0 • or 180 • where the diurnal component is strongest in the plane of the object's orbit. Users that want to artificially decrease or eliminate the effects of the diurnal component to make the seasonal component more prominent can simply increase the body's rotation period to an arbitrarily high value in their simulations.
From this model, observational information about changes in the orbital distance of an orbiting body can be used to in SI units. constrain the physical properties of that body, as was done by Vokrouhlický et al. (2008) to constrain several physical parameters of the asteroid 1992 BF and Tardioli et al. (2017) to constrain the obliquities of multiple near-Earth Asteroids. The Full Version of the effect can also be used to approximate the age of asteroid families by determining how the orbits of their constituent asteroids evolved since the family's creation. This information can then be used to determine the time when these asteroids began diverging from the collision that created them. Nesvorný & Bottke (2004) and Carruba et al. (2017) used this technique to decrease the uncertainty in the ages of the Karin and Veritas asteroid families, respectively.

SIMPLE VERSION
The Simple Version of the Yarkovsky effect in is based on modifications made to Equation (1) by Veras et al. (2019). Recall from Equation (1) that Y is the Yarkovsky matrix that describes an object's rotation and thermal emission. Calculating this matrix between each time step of a simulation is time-consuming, and using it requires many parameters. To prevent the need for such a complicated heat conduction model, Veras et al. (2019) placed constant entries into this matrix. While these modifications lead to a less detailed model of the Yarkovsky effect, they decrease the computational time needed for calculations. Two of these prescriptions (models A and B) are available for use in the Simple Version of the Yarkovsky effect.
These matrices have been modified slightly to remove the effects of PR drag by changing the diagonal terms from 1 to 0. PR drag is already a feature in the previously implemented Radiation Forces effect in (Tamayo et al. 2020), so its inclusion in the Simple Version would be superfluous. In addition, while PR drag is important to consider for dustsized particles, the Yarkovsky effect plays a larger role in the orbital evolution of larger bodies. Veras et al. (2015) showed that the Yarkovsky effect should only be active in bodies with diameters greater than 10 m that are spinning moderately fast. In addition, Veras et al. (2015) show that the Yarkovsky effect is proportional to 1/ while PR drag is proportional to 1/ 2 . Therefore, when the Yarkovsky effect is active on a body, it is reasonable to ignore the effects of PR drag.
Model A and Model B are shown in Equations (6) and (7), respectively. Model A is configured to push targets outward (corresponding to prograde rotation), while Model B is configured to drive targets inward (corresponding to retrograde rotation). Setting the parameter "ye_flag" to 1 on an object in the simulation will apply model A to Y, while setting "ye_flag" to -1 will apply model B to Y. In addition to the simplifications made to Y, it's assumed that = 0.25 for this version of the effect. The equation for the Simple Version of the effect is then where Y will either be Model A or Model B depending on what the user chooses. The following section will show how the change in semi-major axis of an asteroid under the influence of this version of the effect can be represented analytically while Section 5.1 gives a demonstration of using this version of the effect in combination with the parameter interpolation feature in .

Analytic Comparison
If an object is in a circular orbit with no inclination around a star whose mass and luminosity are constant, then an analytic equation describing the time-evolution of the object's semimajor axis can be derived from the equations of the Simple Version. This derivation is made possible due to the Simple Version's modifications to the Yarkovsky matrix, and it cannot be performed using the equations from the Full Version. While this analytic equation only applies to the specific situation described above, we can use it to show that the Simple Version of the effect properly applies Equation (8) to targets in simulations. Equation (A2) in Veras et al. (2015) shows how the average Yarkovsky-induced time rate of change of a body's semimajor axis depends on the entries in Y: where is the semi-major axis of the body, is the mean motion of the body, is the object's inclination, and Ω is the target's longitude of ascending node. If we assume that the object's inclination and eccentricity are both 0 and that the luminosity of its host star is constant with time, then the equation simplifies to To eliminate all time dependency in variables other than , we must get in terms of : where is the gravitational constant and is the mass of the primary. If we assume that the mass of the body is small compared to the mass of the primary, the body orbits a 1 star, and the units we use are astronomical units, , and years so that G has a value of 4 2 , then Equation (11) simplifies to We can now replace in Equation (10) with Equation (12) to get the following separable differential equation: Solving this differential equation, we get where is time and is a constant of integration. We can transform Equation (14) into a form that more resembles Equation (8) by assuming = 2 and has a value of 0.25: where is the radius of the body. By isolating and solving for when = 0, the equation becomes where 0 is the starting semi-major axis of the body. This is the final form of the equation and can be used to replicate and validate certain results obtained from the Simple Version of the effect. To demonstrate this, we created a simple simulation containing an asteroid in a flat, circular orbit around a 1 star. The asteroid has a density of 3000 kg m −3 , a radius of 1 km, and a semi-major axis of 1 au. The star has a constant luminosity of 100,000 to make the effects of the Yarkovsky effect more noticeable on a shorter timescale. For this simulation, we used Equation (6), which will push the body's orbit outward. The simulation lasts for 1 Myr and records how the semi-major axis of the body changes over time due to the Simple Version of the Yarkovsky effect. We then compared the results from the simulation to the results when using the exact same parameters in Equation (16). Figure 2 shows both the results from the simulation and the curve created by the analytic equation. The two results are virtually identical for at least 1 Myr, showing that the Simple Version will properly apply Equation (8) to targets in simulations with substantial durations.

COMBINING EFFECTS
provides users with the capability to easily use multiple forces, effects, and operators in a single simulation, allowing for detailed studies of situations that require a combination of different physics (giant branch studies, protoplanetary disks, etc.). Section 5.1 shows an example that combines the Yarkovsky effect with parameter interpolation while Section 5.2 shows how the different versions of the Yarkovsky effect perform when combined with the Radiation Forces effect in .

Yarkovsky effect and Parameter Interpolation
effects can be used in conjunction with the parameter interpolator in , a feature that allows one to use imported time-series data for parameters that will change throughout the course of a simulation (Baronett et al. 2022). Based on the cubic spline algorithm from Press et al. (1992), this interpolator takes user-inputted data for a particular parameter and creates a cubic spline that can be called upon at any arbitrary time during a simulation to obtain a value for that parameter. This interpolator is written in C, is machineindependent, and supports both forward and backwards integration.
To demonstrate the combined capabilities of parameter interpolation and the Yarkovsky effect, we run a simple simulation based on studies done by Veras et al. (2019) of asteroids under the influence of the Yarkovsky effect during the latestage phases of stellar evolution. The simulation contains a star with a mass of around 2 that is about to ascend the The actual starting mass is 1.951859 . asymptotic giant branch (AGB) and asteroids placed in circular orbits at 1 au, 3 au, and 10 au. These asteroids each have a density of 3000 kg m −3 and can have radii of 10 m, 100 m, and 1000 m. The Simple Version of the Yarkovsky effect is applied to these nine asteroids for the entirety of the simulation, with the effect set to propel these bodies outward. We also add three control asteroids without the Yarkovsky effect active on them at these three different starting positions. To save computational time and keep the focus of the simulation on the effects of stellar mass loss and the Yarkovsky effect, all asteroids in this simulation are semi-active. The simulation lasts for 1.8 Myr and begins about 1.7 Myr before the star reaches the tip of the AGB. The time-series data for the mass and luminosity of the originally 2 star were obtained from the stellar evolution code (Hurley et al. 2000). During the simulation, the star will lose approximately 67% of its starting mass, which will cause all asteroid orbits to expand adiabatically by roughly a factor of three. The luminosity of the star will reach values thousands of times greater than the luminosity of the Sun, which will cause an increase in the strength of the Yarkovsky effect for the duration of the AGB.
The top two panels of Figure 3 show the mass and luminosity of the star during the course of the AGB phase. The bottom panel of Figure 3 shows the semi-major axes of the asteroids as a function of simulation time and the age of the star. As expected, asteroids with smaller radii are propelled outward much faster by the Yarkovsky effect, and the combination of the Yarkovsky effect and stellar mass loss can push small asteroids to orbits that are over 100 times greater than their starting positions. The three asteroids not under the influence of the Yarkovsky effect remain in relatively unchanging orbits until near the end of the simulation, when the star quickly loses a large portion of its mass.
This simulation and previous studies (Zuckerman et al. 2010;Frewen & Hansen 2014;Veras et al. 2022) have shown that the Yarkovsky effect plays a role in determining the final arrangement of rocky material around white dwarfs and is an important mechanism for phenomena such as white dwarf pollution. The Yarkovsky effect in will be a useful tool to further explore the role this radiative force has on small bodies during the late stages in a star's life.

Time Performance
As explained in Section 4, the equations contained in the Simple Version of the effect are less detailed than the ones in the Full Version to save computational time during longer simulations. To demonstrate this difference in time perfor-Setting particles to semi-active in simulations eliminates their gravitational influence on other bodies while retaining the gravitational influence of active particles onto them.  mance, we apply each version of the effect to a 0.01 Myr-long simulation and record its duration. The simulation contains a test particle orbiting the Sun in a circular orbit at 1 au and uses the WHfast integrator (Rein & Tamayo 2015;Rein et al. 2019) with a fixed time step of 0.05 yr. In addition to measuring the two versions of the Yarkovsky effect, we also measure the time performance of the Radiation Forces effect in -, which applies both PR Drag and radiation pressure to chosen particles in a simulation. Measuring this effect will let us compare the performance of our new effect with previously implemented radiative forces.
The following configurations of effects were applied to the simulation: (1) no effects; (2) Radiation Forces; (3) the Simple Version; (4) the Full Version; (5) Radiation Forces and the Simple Version; and (6) Radiation Forces and the Full Version. We performed 10 runs for each combination and calculated their average runtimes and standard deviations. Table 2 shows the results for each configuration. As expected, simulations using the Simple Version require less computational time than ones using the Full Version. On average, we find that the Full Version increases a simulation's duration by a factor of about 2.820 while the Simple Version increases it by a factor of 1.324, which is comparable to the factor of 1.322 from Radiation Forces. We also see a difference in computational times when Radiation Forces is combined with the two different versions of the Yarkovsky F . effect. Radiation Forces combined with the Simple Version increases the length of a simulation by an average factor of 1.831, which is significantly less than the factor of 3.078 from Radiation Forces combined with the Full Version. These results show that the Simple Version of the effect is more useful for long-term simulations containing many bodies, while the Full Version is more useful for short, highly-detailed simulations with fewer bodies.
For the convenience of users, we also provide time performance data for five other integrators available in : IAS15 (Rein & Spiegel 2015), JANUS (Rein & Tamayo 2018), SABA (Rein et al. 2019), EOS (Rein 2020), and Leapfrog. Apart from the integrators, the setup for the simulations and methods for data collection remain the same. Table 3 shows the average durations and standard deviations for simulations using these integrators with the same set of varied combinations of effects active.

CONCLUSION
Two different forms of the Yarkovsky effect have been added to the external library : the Full Version and the Simple Version. The Full Version will be useful for constraining the physical properties of asteroids (Vokrouhlický et al. 2008;Tardioli et al. 2017), estimating the ages of asteroid families (Nesvorný & Bottke 2004;Carruba et al. 2017), and in general, short but detailed simulations containing fewer objects. The Simple Version will be useful for studying the evolution of smaller objects during the post-main sequence (Veras et al. 2019) -which could have important consequences for white dwarf pollution (Zuckerman et al. 2010;Frewen & Hansen 2014;Veras et al. 2022) -and in general, longer simulations with more particles that can tolerate lower levels of detail. We hope these new capabilities for will enable new studies into the role of radiative forces in planetary systems, and we encourage others to make contributions to as well.