Disruption of exo-asteroids around white dwarfs and the release of dust particles in debris rings in co-orbital motion
Kyriaki I. Antoniadou Dimitri Veras
 


   


   

   Department of Mathematics “Tullio Levi-Civita”, University of Padua, 35121, Padua, Italy 
   
 kyriaki.antoniadou@unipd.it 
   Department of Physics, Aristotle University of Thessaloniki, 54124, Thessaloniki, Greece 
   
 
    Centre for Exoplanets and Habitability, University of Warwick, Coventry CV4 7AL, United Kingdom 
   Centre for Space Domain Awareness, University of Warwick, Coventry CV4 7AL, United Kingdom
   Department of Physics, University of Warwick, Coventry CV4 7AL, United Kingdom
         

Co-orbital dynamics in rings

K. I. Antoniadou and D. Veras




  
   Close to the Roche radius of a white dwarf (WD), an asteroid on a circular orbit sheds material that then adopts a very similar orbit. Observations of the resulting debris show a periodic behaviour and changes in flux on short timescales, implying ongoing dynamical activity. Additional encounters from other minor planets may then yield co-orbital rings of debris at different inclinations. The structure, dynamics, and lifetime of these debris discs remains highly uncertain, but is important for understanding WD planetary systems.   
  
   We aim to identify and quantify the locations of co-orbitals in WD–asteroid–dust particle three-body systems by exploring the influence of 1:1 resonant periodic orbits. We begin this exploration with co-planar and inclined orbits in the circular restricted three-body problem (CRTBP) and model the dynamical evolution of these exosystems over observable timescales. The mass ratio parameter for this class of systems (≈ 2 × 10^-11) is one of the lowest ever explored in this dynamical configuration.
  
   We computed the periodic orbits, deduced their linear stability, and suitably seeded the dynamical stability (DS) maps. We carried out a limited suite of N-body simulations to provide direct comparisons with the DS maps.
  
   We derive novel results for this extreme mass ratio in the CRTBP, including new unstable 3D families. We illustrate through the maps and N-body simulations where dust can exist in a stable configuration over observable timescales across a wide expanse of parameter space in the absence of strong external forces.  
  
   Over a timescale of 10 years, the maximum orbital period deviations of stable debris due to the co-orbital perturbations of the asteroid is about a few seconds. Unstable debris in a close encounter with the asteroid typically deviates from the co-orbital configuration by more than about 20 km and is on a near-circular orbit with an eccentricity lower than ≈0.01.

   
   Disruption of exo-asteroids around white dwarfs and the release of dust particles in debris rings in co-orbital motion
    Kyriaki I. Antoniadou<ref>,<ref> Dimitri Veras<ref>,<ref>,<ref>
          
    Received XXXX; Accepted YYYY
======================================================================================================================





§ INTRODUCTION
 

Transiting rocky debris is a signpost of the most dynamically active white dwarf (WD) planetary systems <cit.>. Unlike exoplanet transits around main-sequence stars, which feature a characteristic single solid-body dip in the light curve, minor planets that break up around WDs contain asymmetrical, sharp, shallow, sometimes periodic and sometimes ephemeral transit features on a nightly basis <cit.>.

Although the transit observations are too complex to be explained by a simple disruption model, analytical progress has been made. For WD 1145+017, the first WD discovered with transiting debris <cit.>, signatures corresponding to small periodic dips in the transit curves are assumed to be broken-off fragments and have been leveraged to estimate both the mass (≈ 10^20 kg) and orbital eccentricity (≈ 0) of the progenitor asteroid <cit.>. Furthermore, the frequency and shape of the transit features themselves have been linked to density, composition, and layering of the progenitor asteroid <cit.>.

More generally, the shedding of mass that follows the progenitor orbit has now been investigated extensively <cit.>, as has the process of circularising the debris <cit.>. In contrast, the analytical structure of this phase space has remained relatively unexplored, despite its importance <cit.>. Uncovering the stable and unstable co-orbital regions of these systems can aid in future modelling efforts, and it may help to explain the observations.

Dynamically speaking, WD planetary systems provide extreme examples of the circular restricted three-body problem (CRTBP). The masses of WDs are about half the mass of the Sun, and an orbiting minor planet represents the secondary, leading to extremely low mass ratios. Furthermore, the Roche radius of a WD is located at about 1R_⊙. Hence, the secondary orbits are on the scale of many hours, and spatial distances on the scale of 10^-3 au become important (similar e.g. to Saturn's rings).

The uncertainty about the lifetimes of these debris discs and rings adds to the complexity of these systems <cit.>. If both discs are sufficiently long lived and the flux of asteroids into the WD Roche radius is sufficiently high, then we expect the debris to be replenished anisotropically. Simulations have demonstrated that asteroids (and potentially debris) would enter the Roche radius at a wide variety of inclinations <cit.>. Furthermore, observations do not yet provide a clear picture of the three-dimensional shapes of these debris discs <cit.>. These arguments emphasise the importance of investigating the inclined CRTBP (or 3D-CRTBP) in addition to the co-planar CRTBP (or 2D-CRTBP).

The paper is organised as follows. In Sect. <ref>, a brief literature review of co-orbital motion is provided. In Sect. <ref>, the main aspects of the method are presented, that is, the definition of periodic orbits and mean motion resonances (MMRs), and the linear stability and long-term evolution of systems hosted in specific regions in phase space. In Sect. <ref>, we compute the 2D and 3D families of 1:1 resonant periodic orbits in the CRTBP. In Sect. <ref>, we provide an extended view of the phase space around the WD–asteroid–dust particle system via dynamical stability (DS) maps, which reveal all types of orbits in co-orbital dynamics. In Sects. <ref> and <ref>, we describe the setup and results of the N-body simulations. We discuss our results in Sect. <ref>, and we finally conclude in Sect. <ref>.




§ BRIEF HISTORY AND ASPECTS OF CO-ORBITAL MOTION
 
<cit.> and <cit.> were the first to clearly identify families of what today is referred to as quasi-satellite (QS) orbits in the three-body problem (TBP), based on some preliminary results on retrograde periodic orbits treated by <cit.>, which were then called retrograde satellite orbits. In general, this motion takes place outside of the Hill sphere surrounding the minor primary body of the restricted TBP (RTBP), and these trajectories are much closer to it than to the major primary body <cit.>. These orbits are not related with the tadpole (TP) or horseshoe (HS) orbits and can exist far from the Lagrangian points L_3, L_4, and L_5 <cit.>. 

For clarity, we refer to co-orbital motion as two bodies that share the same orbit, and QS orbits are a special class of TBP orbits that belong to the family of stable symmetric periodic orbits in 1:1 MMR. In the planar circular RTBP (2D-CRTBP), this family of stable periodic orbits was called family f by <cit.>, and <cit.>. The orbits of family f are generated by family E_11^+ <cit.>, when the problem mass parameter μ>0, and they terminate at a collision orbit with the major primary body <cit.>. In Hénon's notation, the generating orbit of family E_11^+ is a third-species orbit (the orbits of the minor primary body and the mass-less body coincide when μ→ 0). The QS orbits of family f are stable for μ<0.0477 <cit.>. 

In the 2D-CRTBP, there exist Lyapounov families with unstable symmetric periodic orbits and stable asymmetric periodic orbits. The former family, called family b by <cit.>, consists of HS orbits that are generated by family E_11^- in Hénon's notation when μ>0, and it includes L_3 (the orbits of the minor primary body and the mass-less body are diametrically opposed when μ→ 0). The latter family, called family E_11^a by Hénon, includes the short-period[The period of short-period orbits tends to 2π, whereas the period of the long-period orbits tends to infinity when μ→ 0.] orbits emanating from L_4 and intersects family E_11^- when the linear stability changes (see Sect. <ref> herein). This family has a mirror image that emanates from L_5. <cit.> called these families of periodic orbits in the rotating frame of the 2D-CRTBP ℒ_3 and short-period ℒ_4^s (or ℒ_5^s). We do not discuss the family of long-period unstable asymmetric periodic orbits, which terminates in the family of short-period asymmetric ones <cit.>, and the families starting from L_1 (prograde satellite motion in family c) and L_2 (family a) <cit.>. 

We instead explore four distinct types of orbits. In our case study, namely WD–asteroid–dust particle, we have the categories listed below. 

        
  * When the motion of the dust particle takes place close to the asteroid and the gravitational perturbation of the WD is assumed to be negligible, the reference system is 'asteroid-centred' and the orbits are almost those of a Keplerian retrograde satellite. The dust particle forms a close binary with the asteroid, and their centre of mass revolves around the WD (Fig. <ref>a). 
        
  * When the motion of the dust particle is quite distant from the asteroid and the gravitational attraction of the WD dominates the orbits, we may use a 'WD-centred' system and compute Keplerian planetary-type orbits (orbits of second kind according to Hénon). We may have three possible configurations:

        (a) The stable configuration of QS_h symmetric orbits is obtained when the bodies are anti-aligned (apsidal difference equal to 180^∘; Fig. <ref>b).

        (b)  The unstable configuration of HS symmetric orbits arises when the bodies are aligned (apsidal difference equal to 0^∘) and move in opposite phases (Fig. <ref>c).

        (c) The stable configuration of TP asymmetric orbits is obtained when the bodies are neither aligned nor anti-aligned (Fig. <ref>d).

        The QS_h orbits and the retrograde satellite orbits are separated by the binary QS_b orbits, which belong to a domain in which none of the primaries influences the mass-less body <cit.>.
          



Osculating orbital elements can be used to describe these orbits, that is, a_i (semi-major axis), e_i (eccentricity), i_i (inclination), ω_i (argument of pericentre), M_i (mean anomaly), and Ω_i (longitude of the ascending node). We also use the notation ϖ_i=ω_i+Ω_i for the longitude of pericentre, Δϖ for the apsidal difference, ΔΩ for the nodal difference, and λ_i=ϖ_i+M_i for the mean longitude. Subscript i= a denotes the asteroid, and i= d signifies the dust particle.
In the CRTBP, these orbits render a system consisting of a dust particle and an asteroid moving around the WD with the same orbital period in elliptic, e_ d>0, and in circular, e_ a=0, orbits, respectively. 

In the Solar System, QS orbits have been computed numerically and analytically for the Mars-Phobos system <cit.>, while the first results were computed for the Phobos spacecraft mission <cit.>. Stable QS orbits around giant planets in our solar neighbourhood were identified by analytical, semi-analytical, or numerical methods <cit.>. Distant retrograde orbits are particularly informative regarding the space mission design to the Moon, Europa, and other satellites and asteroids <cit.>. Asteroid 2003 YN_107 was the first quasi-satellite of Earth that was discovered <cit.>, followed by studies of others in 1:1 resonance with Earth, such as asteroid 2002 AA_29 <cit.>, 2004 GU_9 and 2006 FV_35 <cit.>, 2013 LX_28 <cit.>, and 2014 OL_339 and 2016 HO_3 <cit.>. Other studies involve QS orbits around Venus <cit.>, Jupiter <cit.>, and the main asteroid belt co-orbitals <cit.>. 

<cit.> studied the survival of dust particles in the Didymos-Dimorphos system and their connection with stable QS orbits, while <cit.> explored the QS dust trapping and dynamical evolution in Earth's regime, and <cit.> modelled the evolution of ring particles in the HS regime in the Janus-Epimetheus system. The Yarkovsky-O'Keefe-Radzievskii-Paddack (YORP) effect on a particular class of Earth and Mars co-orbitals was studied for instance by <cit.>.  

The detectability of extrasolar trojans and bodies on QS-like orbits is in general very difficult because they can be discarded as false positives <cit.>. Nevertheless, this has not discouraged theoretical studies <cit.>. 

For reasons of completeness, we refer to the computation of families of QS orbits in the planar general TBP (2D-GTBP). When perturbation was added to the problem and a continuation with respect to the mass of the mass-less body was performed <cit.>, the QS orbits of family f generated stable symmetric periodic orbits in the 1:1 MMR in the so-called family S in the 2D-GTBP computed for various mass ratios by <cit.> and <cit.>. Family S consists of planetary-type (computed in the heliocentric system) and satellite-type orbits (computed in the planetocentric system). This group of families S was later called g(f_1,E_a) by <cit.>. Additionally, numerical, analytical, and semi-analytical treatments of the phase space were conducted for example by <cit.> because planetary configurations in which co-orbital motion could be exhibited are very intriguing <cit.>.   






§ METHOD



 §.§ Notion of periodic orbits related to MMRs and resonant angles

We aim to determine the location of QS, HS, and TP orbits in the 1:1 MMR and explore the influence of the 1:1 resonant symmetric and asymmetric co-planar or inclined periodic orbits on the dynamical behaviour of dust particles sharing the same orbit with an exo-asteroid within a debris disc or ring around a WD. In general, the MMR acts as a phase-protection mechanism safeguarding any planetary system even when the orbits are highly eccentric <cit.>. 

The periodic orbits indicate the exact location of the MMR in phase space and are fundamental for understanding the resonant dynamics of each problem. They coincide with the fixed or periodic points on a Poincaré surface of section and with the stationary equilibrium  points of an appropriately averaged Hamiltonian, as long as the latter is sufficiently accurate <cit.>. The periodic orbits are not isolated in phase space. They are continued mono-parametrically (see Sect. <ref>) and form characteristic curves or families <cit.>. These families of periodic orbits may be either generated by bifurcation points, or are isolated.

We defined the resonant angle θ=λ_ a-λ_ d <cit.> in order to distinguish the various families of periodic orbits in the 2D- and the 3D-RTBPs that dominate and shape each region in phase space. The value of the resonant angle remains constant when computed on the exact periodic orbit. In the following, we discuss the  behaviour of the angles (libration or rotation) that is linked with the linear stability of the periodic orbits and the vicinity of the latter in phase space.  






 §.§ Linear stability and long-term planetary stability

We computed the linear stability <cit.> of the periodic orbits by calculating the eigenvalues of the monodromy matrix of the linearised equations of motion, whose number depends on the degrees of freedom of each problem. If, and only if, all the eigenvalues lie on the unit circle, the periodic orbit is classified as linearly stable. Linearly stable periodic orbits, either planar or spatial, constitute the backbone of stability domains in phase space, where invariant tori exist and the motion is regular and bounded. Therefore, any planetary system therein will survive for long-time spans. In these regions, the resonant angles and the apsidal difference librate about 0^∘ or 180^∘ when the periodic orbit is symmetric, that is, when it is invariant under the fundamental symmetry Σ: (t,x,y) → (-t,x,-y) <cit.>, or about other values when it is asymmetric and Σ maps it to its mirror image <cit.>.

When the periodic orbits are linearly unstable, chaos manifests itself and leads to irregular evolution. This may not change the configuration of the system significantly when the chaos is weak. When the chaos is strong, planetary disruption may occur and lead to close encounters, collisions, or escapes, and hence, to planetary destabilsation. In the vicinities of these chaotic seas, the resonant angles and the apsidal difference rotate. 

We also took an intrinsic property of the planar periodic orbits into account, namely the vertical stability, by computing the vertical stability index <cit.>. A system of celestial bodies moving on different inclined planes can evolve stably when their mutual inclinations are low and when they are located in the neighbourhood of vertically stable 2D periodic orbits. The orbits that are vertically critical (vco) act as bifurcation points that generate 3D periodic orbits. Celestial bodies with a high mutual inclination can only survive when they are located in the neighbourhood of linearly stable 3D resonant periodic orbits <cit.>.

Therefore, the methodical computation of families of periodic orbits and the deduction of their linear stability are precise and rigorous diagnostic tools to guide the exploration of the phase space. In this study, we consider both symmetric and asymmetric  1:1 resonant periodic orbits. Planar stable asymmetric orbits were for example found in the 2D-GTBP by <cit.>. Families of 2D asymmetric periodic orbits were also found to be vertically stable in the 2D-GTBP by <cit.>, and hence, these co-orbital bodies will survive small deviations from co-planar motion.



 §.§ System set-up

We first considered a WD, an asteroid, and a dust particle as point masses, m_ WD, m_ a, and m_ d=0, respectively, and defined the problem parameter as μ=m_ a/m_ WD+m_ a. In this configuration, the motion of the dust particle does not affect the motion of the WD and the asteroid while it moves under their gravitational attraction. In 1:1 resonance, the semi-major axes of the asteroid and the dust particle are almost equal (a_ a≈ a_ d).

In the inertial frames of reference, for instance, OXY or GXYZ, the periodic orbits correspond to almost Keplerian ellipses, which are described by the osculating elements and rotate around the major primary body. In Fig. <ref>, the Keplerian ellipses presented in the 2D inertial frame (left column) are computed for one period, that is, T=2π, and hence, their rotation around the WD is not apparent. 

We introduced suitable rotating frames of reference to reduce the degrees of freedom of the systems and define the periodic orbits. In the 2D-RTBPs, the rotating frame of reference, for instance, Oxy, is centred at the centre of mass of the primaries, that is, the WD (major primary body) and the asteroid  (minor primary body), while the motion of the asteroid is restricted on the Ox-axis <cit.>. The dust particle describes a periodic motion on the Oxy-plane (right column of Fig. <ref>). The exact resonance is defined as the 1:1 resonant periodic motion in this frame.

In the 3D-RTBPs, the dust particle was allowed to evolve on inclined orbits with respect to the orbital plane of the WD and the asteroid (see e.g. the left column in Fig. <ref>). The origin of the 3D rotating frame, Gxyz, called G, coincides with the centre of mass of the primaries, the Gz-axis is perpendicular to the Gxy-plane, and its Gx-axis is directed from the WD to the asteroid <cit.>. 

An orbit, Q(t), is considered as periodic when it satisfies the condition Q(0) = Q(T), where Q(0) is the set of initial conditions (positions and velocities defined in the rotating frame) at t=0, and T is the orbital period. We assumed a solution Q(t)=(x_ d(t),x_ a(t),y_ d(t),ẋ_ a(t),ẋ_ d(t), ẏ_ d(t)) in the 2D-RTBP.  This solution is  periodic with a period T when it satisfies the periodic conditions

    [ ẋ_ a(T)=ẋ_ a(0)=0,                   ;   x_ a(T)=x_ a(0),                   ;  x_ d(T)=x_ d(0) ,   y_ d(T)=y_ d(0),;   ẋ_ d(T)=ẋ_ d(0),   ẏ_ d(T)=ẏ_ d(0). ]
 
Given the Poincaré surface of section at y_ d=0 and the fact that  x_ a=1-m_ a/m_ WD+m_ a is constant and defined by the normalisation we used in the CRTBP, a symmetric periodic orbit that perpendicularly crosses the Ox-axis (ẋ_ d(0)=0) is defined as a point on the plane of initial conditions {(x_ d(0),ẏ_ d(0))}. The asymmetric periodic orbit can for example be defined in the space {(x_ d(0), y_ d(0), ẋ_ d(0), ẏ_ d(0))} when ẋ_ a(0)= 0. 

For the 3D periodic orbits shown in the right column of Fig. <ref>, the xz-symmetric periodic orbits (top panel) can be represented by a point in the 3D space of the initial conditions
{(x_ d(0),z_ d(0),ẏ_ d(0))} and the x-symmetric ones (bottom panel) by {(x_ d(0),ẏ_ d(0),ż_ d(0))}.

The families are formed during the mono-parametric continuation when a parameter, for example, z_ d, is changed for the xz-symmetry, and the rest are differentially corrected in order to satisfy the respective periodicity condition. Continuation with respect to the mass of a body is also possible for example by varying the mass of m_ a and keeping z_ d fixed, for instance, while the respective periodicity condition still holds. The limitations of this continuation method can be found in <cit.>. The periodicity conditions that the periodic orbits must fulfil were described in greater detail for instance in <cit.> for the 2D-RTBPs and in <cit.> for the 3D-RTBPs.

To numerically computate the periodic orbits we considered, we located without loss of generality the asteroid at a_ a=1, and its period was therefore T_0=2π. We also normalised the masses and the gravitational constant, G, so that the sum of the former was equal to unity, namely m_ WD+m_ a+m_ d=1 or m_ WD=1-m_ a and G=1.

In physical units, given the WD 1145+017 system, we adopted a WD mass m_ WD = 0.6m_⊙. Then, we assumed an asteroid with a mass equal to approximately 1:10^ th the mass of Ceres, for instance, 2.5 × 10^19 kg, located at a_ a=0.0054 au around the WD <cit.>. In this configuration, we derived a μ equal to 2.1186× 10^-11. To our knowledge, this value is one of the lowest ever studied for the 2D and 3D CRTBP. A comparably low value of μ (10^-9) for the numerical computation of QS and HS orbits in the CRTBP was adopted by <cit.> for the specific case of retrograde orbits in the 2D-CRTBP and by <cit.> for the Saturn-Janus-Epimetheus system, respectively. <cit.> explored the domains around QS orbits for μ≥ 10^-7. Finally, numerical computations yielding the extent of regular domains around HS and TP orbits in 1:1 MMR as a function of μ were performed by <cit.>.



 §.§ Visualisation of the phase space

In addition to the linear stability of the periodic orbits, which provides a direct insight into the long-term stability of any system in their vicinity, we also computed DS maps with the use of a chaotic indicator to visualise the extent of each domain. We adopted a version of the FLI, called detrended fast Lyapunov indicator (DFLI), that was established as reliable and accurate by <cit.>. We constructed 200 × 100 grid planes and chose a maximum integration time for the computation of the DFLI equal to t_ max=2.5 Myr, that is, almost 4.8 billion orbits of the dust particle, which was deemed adequate for these systems. We halted the integrations either when DFLI(t)>30 or when t_ max was reached. Orbits with DFLI<2 were classified as stable. Extensive numerical computations have shown that a system can be calssified as chaotic when DFLI>15, as the DFLI continues to increase steeply thereafter <cit.>. 





§ FAMILIES OF PERIODIC ORBITS IN THE CRTBP


In this model for the 1:1 MMR, the asteroid moves on a circular orbit (e_ a=0) with a_ a=1, and the dust particle moves on an eccentric orbit (e_ d>0) with a_ d≈ 1. In the following sections, we present the families of 2D and 3D periodic orbits in the CRTBP. 



 §.§ 2D-CRTBP











In Fig. <ref> we present the families of periodic orbits in the 2D-CRTBP on the (x_ d,e_ d) plane, namely the QS and HS orbits. In Fig. <ref>, we present the families I (stable branch) and II (unstable branch) on the (a_ d,e_ d) plane and provide a comparison with μ=0.001 to highlight the fact that we approach the unperturbed case (μ→ 0). Symbol I corresponds to family f which consists of stable (blue coloured) periodic orbits, where the asteroid and the dust particle are in anti-alignment (Δϖ=180^∘ and θ=0^∘), while symbol II denotes family b of unstable (red coloured) periodic orbits, where the bodies are in alignment (Δϖ=0^∘ and θ=180^∘). We note that the two branches shown for each family are equivalent and differ only in the positive or negative direction (namely ẏ_ d>0 or ẏ_ d<0) of the orbit for which the Poincaré surface of section at y_ d=0 was chosen. More particularly, we obtain the configurations shown in Table <ref>, which are equivalent in pairs, that is, I with I' and II with II', at t=0 and at t=T/2. We provide them all here because the continuation process in the 3D-CRTBP, 2D- and 3D-ERTBP (which will be shown in a future study) was computationally easier when we started from the same bifurcation point but with a different cross section, since cusps appear in the rotating frame of reference. 

In Fig. <ref> we show the family of asymmetric periodic orbits, which consists of TP orbits, and starts from L_4 at θ=60^∘ and ends at a symmetric periodic orbit of HS type at θ=180^∘, where the linear stability changes. We also present the mirror family, which starts from L_5 at θ=300^∘. These families are both horizontally and vertically stable (solid blue lines).




 §.§ 3D-CRTBP

We considered the families of QS, HS, and TP orbits presented in the 2D-CRTBP and computed their vertical stability. Two types of vcos generate 3D periodic orbits with two different symmetries <cit.>, namely with respect to the xz-plane or the x-axis of the rotating frame Gxyz. We call the generated families F and G, respectively, and use a superscript to denote the planar family to which the vco (bifurcation point) belongs. We add a hat to represent the respective vcos. For example, vco F^I belongs to the 2D family I and generates the 3D family F^I of xz-symmetric (magenta vco) periodic orbits. In Table <ref>, we provide the eccentricity values for each vco. The family of asymmetric periodic orbits (TP orbits) is vertically stable.




The 3D families of unstable orbits, which stem from the unstable family II' (or II) of HS orbits in the 2D-CRTBP, have not been presented before for any value of μ. The same holds for the vcos that generate them. For the family F^I of stable QS orbits, we find no significant change in the location of the vco as μ→ 0, that is, as we approach the unperturbed case (m_ a=m_ d=0) (see the vcos (magenta dots) in panel a of Fig. <ref> and compare e.g. the F^I (or F^I') with the B_cs points in Table 1 in <cit.>). Nonetheless, we note that as μ→ 0, the QS domain becomes reachable by low-eccentricity orbits and the trajectories approach the minor primary, that is, the asteroid in our case, while for μ=0.001, this domain exists for e_ d>0.18 (see the top panel of Fig. <ref>).








In Fig. <ref> we present the three 3D families of periodic orbits in the 3D-CRTBP on the (e_ d,i_ d) plane. The family F^I consists of stable xz-symmetric periodic orbits (QS orbits) up to i_d, max=30.7^∘. Both F^II' and G^II' have only unstable periodic orbits (HS orbits) with an xz- and x-symmetry, respectively. In Table <ref> we present the configuration of each family. In Fig. <ref>, each family is shown in more detail in (x_ d,ẏ_ d, z_ d) space or (x_ d,ẏ_ d, ż_ d) space, depending on the symmetry of the periodic orbits together with the projections of the planar families and the vcos in the 2D-CRTBP (at z_ d=0 or ż_ d=0).





§ STABLE AND CHAOTIC QS, HS, AND TP ORBIT REGIMES IN THE 3D CRTBP

We identified the regular domains in phase space on the DS maps by monitoring the libration of the resonant angle θ. When it librated about 0^∘, 180^∘ or other values, we labelled the region QS, HS, or TP_4 (or TP_5), respectively. In the pale regions, the resonant angle θ rotates. All types of orbits in co-orbital dynamics are revealed in the DS maps with the boundaries of different domains having been delineated <cit.>. The breakdown of the phase space was organised as follows.

First, we chose three 3D periodic orbits (one from each family, namely F^I, F^II' and G^II') with mutual inclinations, Δ i, equal to 10^∘, 30^∘, and 50^∘ (Figs. <ref>-<ref>). Then, we varied some of the initial conditions to create the grids on the (e_ d,i_ d), (ω_ d,e_ d), and (M_ d,e_ d) plane, while keeping the rest fixed. The values of the angles are given in Table <ref>, and the orbital elements of the selected periodic orbits are also provided in each figure. 

In Fig. <ref>,we chose three different 3D periodic orbits in which the dust particle had an inclination equal to 10^∘. In the top row, the periodic orbit that guided the search is stable and belongs to the F^I family with a_ d=1+1.887× 10^-11, e_ d=0.718, and ΔΩ =270^∘. The areas around the stable 3D orbit are populated with regular orbits, and chaoticity becomes apparent for e_ d>0.7 in the top left panel, where the planar family, I, becomes vertically unstable (dashed line in Fig. <ref>). In the middle and bottom rows, the periodic orbits that guided the search are unstable and belong to the F^II' and G^II' families, respectively. The xz-symmetric family has initial conditions a_ d=1-3.4903× 10^-11, e_ d=0.722, and ΔΩ =90^∘, and the x-symmetric family has a_ d=1-2.7865× 10^-11, e_ d=0.923, and ΔΩ =0^∘. The regions around both of them are dominated by irregular orbits. However,  the remaining regular domains QS and TP_4 (or TP_5) become apparent farther away from this periodic orbit with changing angles.














For Fig. <ref>, we performed the same search, but chose three different 3D periodic orbits in which the dust particle had an inclination equal to 30^∘. We only present grids on the (ω_ d,e_ d) and (M_ d,e_ d) plane, as the variation in the (e_ d,i_ d) plane was qualitatively similar with that presented in the left column of Fig. <ref>. In the left column of Fig. <ref>, the stable 3D periodic orbit has a_ d=1+2.366× 10^-11, e_ d=0.881, and ΔΩ =270^∘. This 3D periodic orbit is shown in panel a of Fig. <ref>. In the middle column, the 3D unstable orbit has a_ d=1-3.291× 10^-11, e_ d=0.761, and ΔΩ =90^∘, and in the right column, the initial conditions are a_ d=1-2.449× 10^-11, e_ d=0.904, and ΔΩ =0^∘. 

For Fig. <ref>, we performed the same search, but chose two unstable 3D periodic orbits in which the dust particle had an inclination equal to 50^∘ with a_ d=1-3.076× 10^-11, e_ d=0.845, and ΔΩ =90^∘ (left column) and a_ d=1-2.685× 10^-11, e_ d=0.858, and ΔΩ =0^∘ (right column). The latter 3D periodic orbit is shown in panel b of Fig. <ref>. The chaotic domains created around the 3D HS orbits span greater areas with increasing inclination and reach even lower eccentricity, e_ d, values for the dust particle.

In Fig. <ref> we created DS maps on the (ω_ d,M_ d) plane around three 2D asymmetric periodic orbits with e_ d=0 (left panel), e_ d=0.5 (middle panel), and e_ d=0.9 (right panel). The regular orbits dominate as the angles vary as long as e_ d does not approximate the bifurcation point on the HS orbit, that is, at 0.917.

In Fig. <ref> we varied the semi-major axis and the eccentricity of the dust particle on the (a_ d,e_ d) plane. In the top and middle rows, we used the same six periodic orbits as in Figs. <ref> and <ref> with Δ i=10^∘ and Δ i=30^∘, respectively. Therefore, islands of stability around the QS orbits were created in the DS maps. In the bottom row, we explored the spatial neighbourhood of the planar vertically stable asymmetric TP orbit close to L_4 (or L_5) for Δ i=10^∘ and three different nodal differences, namely ΔΩ=0^∘, 90^∘ and 180^∘. These numerical results and the boundaries of the domains, namely the libration of the resonant angle, agree with those presented by <cit.>.




§ SETUP OF THE N-BODY SIMULATIONS


After computing periodic orbits and creating DS maps, we performed N-body simulations. In order to accurately compare these simulations to the DS maps, care must be taken with regard to rescaling <cit.>.

The setup of the N-body simulations requires that the data extracted from the families of periodic orbits and the DS maps are transformed into real units. In this case, the results presented here can be applied to different planetary system configurations. For instance, different units of time, masses, and distances can be used. We needed to set a normalisation with respect to the masses assumed for the computation of the periodic orbits in our case. The details of the invariance of the equations of motion we used to compute the periodic orbits and the appropriate scaling factors were explained by <cit.> for the 2D-CRTBP and by <cit.> for the 2D-ERTBP and the two 3D-RTBPs.

We performed a scaling in our simulations with regard to the semi-major axis of the asteroid and the dust particle, so that 
 
    a_i^ (N) = a_i ζ^1/3,

 where ζ≡m_ WD + m_ a/1m_ is the scaling factor, and N represents the scaled orbital elements used in the N-body simulations. Therefore, the  values of the semi-major axes that were imported to the N-body integrator were multiplied by the factor 0.843432665311676. Time and the remainder of the orbital elements (e.g. the eccentricities and angles) remained the same. 

For all of our N-body simulations, we used the IAS15 integrator in the REBOUND simulation package <cit.>. This integrator uses an adapative time step and conserves energy sufficiently well for our purposes. In order to model a realistically observable timescale, we performed each simulation over 10 years, corresponding to roughly 20,000 orbits of the minor planet. We output data snapshots every 0.05 years, but monitored if and when collisions occurred for every time step. Although our test particles were point masses, our minor planets were not: Their radius was 200 km.

Limiting the integration timescale to 10 years also prevented longer-term forces on the dust particles from becoming important. Several shorter-term forces might also act on the particles. We discuss how these can be contextualized within the outputs of our simulations in Sect. <ref>. These forces typically have a greater effect the smaller the orbital pericentre. In this context, the most consequential results are those with e_ d≈ 0.0, whereas those with more limited applicability have e_ d≳ 0.8. Nevertheless, we sampled e_ d values up to 0.9 in many cases to provide a direct comparison with the DS maps.




§ RESULTS OF THE N-BODY SIMULATIONS


We had a large number of parameters to select in order to report our results. We chose a parameter that has a direct observable link: the maximum variation in the orbital period of the dust particle. The orbital period is directly measured and is independent of other parameters. Transit photometry in these WD systems can be measured on the scale of seconds <cit.>.

We wished to determine the extent of this orbital period variation and how it relates, if it does, to periodic orbits and the structure seen in the DS maps. Because of the computational expense of integrating over a grid of initial conditions with an accurate time-variable integrator such as IAS15, we selected a few parts of a few parameter grids over which to perform simulations. These are outlined below. In all cases, each parameter was sampled uniformly.

For some parameter combinations, some systems become unstable. These instabilities are manifested entirely through collisions between a dust particle and the asteroid. In our simulations, no dust particles escape the system and no dust particles collide with the WD. Unstable systems are indicated by red crosses in the plots.

Within the set of six orbital parameters, the only parameter that determines the location of the dust particle rather than setting the shape of its orbit is M_ d. In some cases, in order to provide a direct comparison to the DS maps, we adopted the value of M_ d that corresponds to a particular family of periodic orbit, or DS map. However, in some other cases, we sampled  a range of M_ d values in order to determine the maximum period deviation independent of M_ d, and we report only those that generated the maximum orbital period deviation.



 §.§ Variation in e_ d versus a_ d
 

Fig. <ref> shows a 91 × 91 grid of a_ d = a_ a∓ 68.135  km and e_ d = 0.0001-0.1000 (left column) and e_ d = 0.0001-0.9000 (right column). For each combination of (a_ d, e_ d), we performed 18 simulations with M_ d uniformly sampled across its entire range. The maximum orbital period deviation occurs along the boundaries of the island of stability that is populated by low-eccentricity orbits with ΔΩ=270^∘ and Δ i=10^∘.




A similar behaviour is shown in Fig. <ref>, but guided by the 3D HS periodic orbit with ΔΩ=90^∘ and Δ i=10^∘. Again, the maximum orbital period deviation takes place along the boundaries of the island where ΔΩ=90^∘. 



Fig. <ref> shows a grid guided by the 3D HS periodic orbit with ΔΩ=0^∘ and Δ i=10^∘. Again, the maximum orbital period deviation occurs along the boundaries of the respective island of stability. 




Fig. <ref> shows a grid guided by a TP periodic orbit with ω_ d=60^∘ and Δ i=10^∘. In contrast to the previous simulations, we performed a single simulation with M_ d = 0^∘. Hence, the value of the maximum orbital period deviation observed in the region populated by TP orbits was lower. 






 §.§ Variation in e_ d versus ω_ d


Fig. <ref> shows a 91 × 91 grid of ω_ d = 0^∘-360^∘ and e_ d = 0.0001-0.9000 guided by the 3D QS periodic orbit with ΔΩ=270^∘ and Δ i=10^∘. For each combination of (ω_ d, e_ d), we performed a single simulation with M_ d = 180^∘ (left column) and 18 simulations with M_ d uniformly sampled across its entire range (right column).  When M_ d = 180^∘, the correspondence between the DS map and the N-body simulations is excellent, although the maximum orbital period deviation observed was negligible in this exploration. This deviation did not increase even when a uniform sample was adopted for M_ d.





 §.§ Variation in i_ d versus e_ d


Fig. <ref> shows a 181 × 181 grid of e_ d = 0.0001-0.8000 and i_ d = 0^∘-90^∘ guided by the HS periodic orbit with Ω_ d=90^∘ and Δ i=50^∘, and we performed a single simulation with M_ d = 180^∘. The  maximum orbital period deviation observed around the 3D HS orbits was negligible, as the semi-major axis of the dust was fixed and equal to that of the periodic orbit. 






§ DISCUSSION


Accurately modelling WD debris discs requires incorporating a vast variety of physical processes, including gas-dust interactions, gas drag, aeolian erosion, Ohmic heating, Lorentz drift, and external perturbations. We have isolated the effect of an important and fundamental force that acts on all the discs: the effect of gravity, and in the likely co-orbital scenario with a disrupted asteroid.

Nevertheless, there are some key physics that we highlight below and relate to our results.



 §.§ General relativity


Similar to Newtonian gravity, general relativity affects all objects orbiting WDs. Furthermore, the effect is greater than that of most known main-sequence extrasolar systems because WD debris discs are so compact. 

Over the course of a single orbit, general relativity changes each orbital element <cit.>. However, these changes all average out to zero, except for the change in the argument of pericentre. Hence, for circular orbits, the effect of general relativity is negligible. For non-circular orbits, in WD debris discs, ω_ d cycles across its entire range in a time of approximately



    107  yr( M_⋆/0.60 M_⊙)^-3/2( a_ d/1.00R_⊙)^5/2( 1-e_ d^2 ).


As a result, we might imagine a dust particle traversing a horizontal path across the plots in Figs. <ref>-<ref>. However, because the maximum period deviation in the entire phase space is less than a few seconds, this traversal would probably not be observable, regardless of the manner in which it takes place.




 §.§ Sublimation




The creation of gas through dust sublimation generates gas-dust interactions, which are still poorly understood and are difficult to model analytically. Hence, our results are best applied to dust alone.

The location in which dust sublimates depends on the temperature of the WD <cit.>, which in turn depends on how long it has been a WD. This time span is commonly known as its cooling age. WDs have been observed with cooling ages of up to ≈10 Gyr <cit.>. Old WDs hence allow the orbital pericentre of a dust particle to be far inwards of 1R_⊙ without sublimating.

In order to quantify this effect, we computed the maximum value of e_ d allowed for a dust particle to survive as dust before it is sublimated as a function of WD cooling age. To do this, we used the same prescription as in <cit.>, and we report our results in Fig. <ref>. This prescription involves combining equations for the blackbody radiation of the WD, the luminosity of a WD as a function of cooling age, and a relation between effective temperature, distance, and sublimation temperatures of three different materials. The plot considers three different potential compositions of the dust and shows that in many cases, values of e_ d≈ 0.8-0.9 may realistically be adopted.





 §.§ Poynting-Robertson drag


Even when a dust particle is not sublimated, the radiation from the WD will drag the particle inwards by Poynting-Robertson drag. This drag occurs for all dust particles and for pebbles, boulders, and any other objects that would be considered as test particles in this study. Because the speed at which Poynting-Robertson drag acts strongly depends on the size of the particle, the relevance of this effect in our study in turn strongly depends on our assumption of the size of the test particle.

We can estimate the inward drift of a particle through the classic <cit.> formula for the semi-major axis variation due to Poynting-Robertson drag,


    da/dt≈ - 3L_⋆(2 + 3 e_ d^2)/16 π c^2 R_ dρ_ d a_ d(1 - e_ d^2 )^3/2,


and we assumed that the particle density is ρ_ d = 2 g/cm^3.

Figure <ref> illustrates the drag rate for four different particle sizes that are large enough to be more appropriately classed as boulders. Over a 10 yr time span, even the largest boulders would be dragged inwards to a greater extent than their orbit would be gravitationally shifted due to a minor planet in a stable co-orbital configuration. We therefore suggest that observational evidence of a period shift is much more likely to arise from Poynting-Robertson drag than from co-orbital dynamics.




We can further deduce that for e_ d≈ 0, based on the x-axes of the N-body simulation plots for Figs. <ref>-<ref>, a dust particle that starts in a stable co-orbital configuration with a minor planet can be dragged inwards by ≈20-60 km into a region of gravitational instability. In this case, gravity might ultimately generate a collision between the particle and minor planet, even though the trigger was radiative drag.




§ CONCLUSIONS


Motivated by asteroids and dust in co-orbital configurations around WDs, we computed the families of symmetric and asymmetric periodic orbits in the 2D and 3D CRTBP for the WD–asteroid–dust particle configuration with μ≈ 10^-11 in the 1:1 MMR. The 3D families with unstable periodic orbits presented here are novel. In this configuration, QS orbits were additionally found as e_ d→ 0. Nonetheless, we found that the planar families exhibited qualitatively similar attributes with those computed for μ=0.001 by <cit.> and <cit.>, that is, their linear stability and their bifurcation points, even though we approached the unperturbed case. 

We methodically explored the phase space by choosing specific 3D symmetric and 2D asymmetric periodic orbits and monitored the libration of the resonant angle. All types of orbits in co-orbital dynamics were revealed via the DS maps, in which the boundaries of different domains were delineated.

All types of orbits were simulated over a timescale of 10 yr with N-body simulations, namely QS, HS, and TP orbits in the 1:1 MMR. In all cases, an observable orbital period deviation was exhibited in the low-eccentricity regime e_ d<0.1 within each domain and nodal difference, while the maximum (about 3 seconds) was yielded when ΔΩ=90^∘ while drifting away from the periodic orbit (a_ d≠ 1) and being almost at the boundaries of the stability domain. 

Before circularising, an asteroid approaches a WD on an eccentric orbit. Hence, a future useful extension of this work may investigate the three-body problem with an eccentric asteroid at a similarly extreme mass ratio. Even so, whether circular or eccentric, a companion of a WD with transiting dust may have another value of μ when a planetesimal is considered. The 2D and 3D periodic orbits play a crucial role in unravelling and determining possible regimes with expected transits <cit.>[Likewise for possible detections in the GTBP, <cit.>]. How large μ needs to be in order for the co-orbital dynamics to become dominant, that is, with an orbital period variation about 3D periodic orbits, namely a detection, remains to be explored.

 We thank the anonymous reviewer for their astute comments and careful reading which improved the manuscript. The work of KIA was supported by the University of Padua under Grant No. BIRD232319. Results presented in this work have been produced using the Aristotle University of Thessaloniki (AUTh) High Performance Computing Infrastructure and Resources. 

aa










