EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 690 (October 2024) A&A, 690 (2024) A249 Full HTML
Open Access
Issue
A&A
Volume 690, October 2024
Article Number A249
Number of page(s) 18
Section Celestial mechanics and astrometry
DOI https://doi.org/10.1051/0004-6361/202451714
Published online 11 October 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

Transiting rocky debris is a signpost of the most dynamically active white dwarf (WD) planetary systems (Vanderburg et al. 2015; Vanderbosch et al. 2020, 2021; Guidry et al. 2021; Farihi et al. 2022). 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 (Gänsicke et al. 2016; Rappaport et al. 2016; Gary et al. 2017; Izquierdo et al. 2018; Rappaport et al. 2018; Aungwerojwit et al. 2024).

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 (Vanderburg et al. 2015), 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 (≈1020 kg) and orbital eccentricity (≈0) of the progenitor asteroid (Rappaport et al. 2016; Gurri et al. 2017). Furthermore, the frequency and shape of the transit features themselves have been linked to density, composition, and layering of the progenitor asteroid (Veras et al. 2017; Duvvuri et al. 2020).

More generally, the shedding of mass that follows the progenitor orbit has now been investigated extensively (Debes et al. 2012; Veras et al. 2014; Malamud & Perets 2020a,b; Li et al. 2021; Brouwers et al. 2022, 2023), as has the process of circularising the debris (Veras et al. 2015, 2022; Nixon et al. 2020; O’Connor & Lai 2020; Malamud et al. 2021). In contrast, the analytical structure of this phase space has remained relatively unexplored, despite its importance (Veras et al. 2016). 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 1 R⊙. 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 (Girven et al. 2012; Veras & Heng 2020). 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 (Veras et al. 2021). Furthermore, observations do not yet provide a clear picture of the three-dimensional shapes of these debris discs (Manser et al. 2016, 2021; Ballering et al. 2022; Goksu et al. 2024). 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. 2, a brief literature review of co-orbital motion is provided. In Sect. 3, 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. 4, we compute the 2D and 3D families of 1:1 resonant periodic orbits in the CRTBP. In Sect. 5, 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. 6 and 7, we describe the setup and results of the N-body simulations. We discuss our results in Sect. 8, and we finally conclude in Sect. 9.

2 Brief history and aspects of co-orbital motion

Brown (1911) and Jackson (1913) 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 Darwin (1897), 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 (Hill 1878; Kogan 1988a, 1990; Mikkola & Innanen 1997). These orbits are not related with the tadpole (TP) or horseshoe (HS) orbits and can exist far from the Lagrangian points L3, L4, and L5 (see e.g. Christou 2000a; Nesvorný et al. 2002; Mikkola et al. 2006, for the behaviour of the phase space around these different types of orbits).

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 2D-CRTBP, this family of stable periodic orbits was called family f by Jackson (1913); Broucke (1968); Hénon (1969); Hénon & Guyot (1970); Benest (1974); Bruno (1994), and Pousse et al. (2017). The orbits of family f are generated by family E11+$E_{11}^ + $ (Hénon 1997), when the problem mass parameter µ > 0, and they terminate at a collision orbit with the major primary body (Hénon 1969). In Hénon’s notation, the generating orbit of family E11+$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 (Hénon & Guyot 1970; Benest 1974).

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 Hénon (1997), consists of HS orbits that are generated by family E11−$E_{11}^ - $ in Hénon’s notation when µ > 0, and it includes L3 (the orbits of the minor primary body and the mass-less body are diametrically opposed when µ → 0). The latter family, called family E11a$E_{11}^a$ by Hénon, includes the short-period1 orbits emanating from L4 and intersects family E11−$E_{11}^ - $ when the linear stability changes (see Sect. 4 herein). This family has a mirror image that emanates from L5. Pousse et al. (2017) called these families of periodic orbits in the rotating frame of the 2D-CRTBP ℒ3 and short-period ℒ4s${\cal L}_4^s$ (or ℒ5s${\cal L}_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 (see e.g. Deprit et al. 1967), and the families starting from L1 (prograde satellite motion in family c) and L2 (family a) (see e.g. Hénon 1969, 1997).

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. 1a).

  • 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 QSh symmetric orbits is obtained when the bodies are anti-aligned (apsidal difference equal to 180°; Fig. 1b).

    • (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. 1c).

    • (c)

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

    The QSh orbits and the retrograde satellite orbits are separated by the binary QSb orbits, which belong to a domain in which none of the primaries influences the mass-less body (see e.g. Pousse et al. 2017; Voyatzis & Antoniadou 2018).

Osculating orbital elements can be used to describe these orbits, that is, ai (semi-major axis), ei (eccentricity), ii (inclination), Ωi (argument of pericentre), Mi (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 + Mi 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, ed > 0, and in circular, ea = 0, orbits, respectively.

In the Solar System, QS orbits have been computed numerically and analytically for the Mars-Phobos system (Lidov & Vashkov’yak 1994), while the first results were computed for the Phobos spacecraft mission (Kogan 1988b; Sagdeev & Zakharov 1989). Stable QS orbits around giant planets in our solar neighbourhood were identified by analytical, semi-analytical, or numerical methods (see e.g. Mikkola & Innanen 1997; Namouni 1999; Namouni et al. 1999; Wiegert et al. 2000; Christou 2000a,b; Nesvorný et al. 2002; Mikkola et al. 2006; Sidorenko et al. 2014; Pousse et al. 2017). Distant retrograde orbits are particularly informative regarding the space mission design to the Moon, Europa, and other satellites and asteroids (see e.g. Minghu et al. 2014; Perozzi et al. 2017; Kimura et al. 2019; Pires & Winter 2020; Li et al. 2022). Asteroid 2003 YN107 was the first quasi-satellite of Earth that was discovered (Connors et al. 2004), followed by studies of others in 1:1 resonance with Earth, such as asteroid 2002 AA29 (Wajer 2009), 2004 GU9 and 2006 FV35 (Wajer 2010), 2013 LX28 (Connors 2014), and 2014 OL339 and 2016 HO3 (de la Fuente Marcos & de la Fuente Marcos 2014, 2016, respectively). Other studies involve QS orbits around Venus (e.g. Mikkola et al. 2004), Jupiter (e.g. Kinoshita & Nakai 2007; Namouni & Morais 2018), and the main asteroid belt co-orbitals (e.g. Christou 2000b).

Madeira et al. (2024) studied the survival of dust particles in the Didymos-Dimorphos system and their connection with stable QS orbits, while Kortenkamp (2013) explored the QS dust trapping and dynamical evolution in Earth’s regime, and Williams & Murray (2011) 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 de la Fuente Marcos & de la Fuente Marcos (2018, 2021).

The detectability of extrasolar trojans and bodies on QS-like orbits is in general very difficult because they can be discarded as false positives (Laughlin & Chambers 2002; Leleu et al. 2019). Nevertheless, this has not discouraged theoretical studies (Raymond et al. 2023a,b).

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 (Hadjidemetriou 1975), 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 Hadjidemetriou et al. (2009) and Hadjidemetriou & Voyatzis (2011). 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 𝑔(f1, Ea) by Voyatzis & Antoniadou (2018). Additionally, numerical, analytical, and semi-analytical treatments of the phase space were conducted for example by Schwarz et al. (2009); Giuppone et al. (2010); Robutel & Pousse (2013); Leleu et al. (2017, 2018); Niederman et al. (2020); Couturier et al. (2022); Sidorenko (2024) because planetary configurations in which co-orbital motion could be exhibited are very intriguing (see e.g. Giuppone et al. 2012; Funk et al. 2013; Lillo-Box et al. 2018; Balsalobre-Ruza et al. 2023).

thumbnail Fig. 1

Periodic orbits in the 2D-CRTBP, where the asteroid, denoted by A (in blue), moves on circular orbits (ea = 0), and the dust particle, D (in red), moves on Keplerian ellipses that rotate around the WD (grey dot), when viewed on an inertial frame of reference (left column). These orbits are exactly periodic in the rotating frame of reference (right column). The orbits were computed for one period, T = 2π, and both bodies orbit the WD with the same period (corresponding in physical space to a distance unit of about 1 R⊙ and to an orbital period of about 4.5 hours). The red dots represent the position of the dust particle at t = 0. In panel a, the stable retrograde satellite orbits that belong to family f are demonstrated, where the dust particle moves on almost circular orbit. In panels b and c, the symmetric stable QSh orbits of family f, and the symmetric unstable horseshoe (HS) orbits are shown, respectively. In panel d, the stable asymmetric tadpole (TP) periodic orbits are shown. The eccentricity of the dust particle chosen for these panels was ed = 0.7, and the initial conditions were taken from the families presented in Fig. 5.

3 Method

3.1 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 (Antoniadou & Voyatzis 2016).

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 (Hadjidemetriou 1993). The periodic orbits are not isolated in phase space. They are continued mono-parametrically (see Sect. 3.3) and form characteristic curves or families (Hadjidemetriou 1975; Hénon 1997). These families of periodic orbits may be either generated by bifurcation points, or are isolated.

We defined the resonant angle θ = λa − λd (Murray & Dermott 1999) 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.

3.2 Linear stability and long-term planetary stability

We computed the linear stability (Marchal 1990) 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) (see Hénon 1973), or about other values when it is asymmetric and Σ maps it to its mirror image (Voyatzis & Hadjidemetriou 2005).

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 (Hénon 1973). 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 (see e.g. Antoniadou & Voyatzis 2014).

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 Giuppone et al. (2010); Hadjidemetriou & Voyatzis (2011). Families of 2D asymmetric periodic orbits were also found to be vertically stable in the 2D-GTBP by Antoniadou et al. (2014a), and hence, these co-orbital bodies will survive small deviations from co-planar motion.

thumbnail Fig. 2

Stable (QSh) and unstable (HS) inclined periodic orbits in the 3D-CRTBP with id = 30° (panel a) and id = 50° (panel b). These orbits belong to the 3D families FI with periodic orbits symmetric with respect to the xz-plane and GII′ with periodic orbits symmetric with respect to the x-axis, (both families are shown Fig. 6). These 3D families may either bifurcate from the 2D family I (or I′) or from the 2D family II (or II′) shown in Figs. 3-4. Presentation as in Fig. 1.

3.3 System set-up

We first considered a WD, an asteroid, and a dust particle as point masses, mWD, ma, and md = 0, respectively, and defined the problem parameter as μ=mamWD+ma$\mu = {{{m_{\rm{a}}}} \over {{m_{WD}} + {m_{\rm{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 (aa ≈ ad).

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. 1, 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 (see e.g. Hadjidemetriou & Christides 1975). The dust particle describes a periodic motion on the Oxy-plane (right column of Fig. 1). The exact resonance is defined as the 1:1resonant 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. 2). 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 (see e.g. Szebehely 1967; Murray & Dermott 1999, for more details on various aspects of the RTBPs).

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)=(xd(t),xa(t),yd(t),x˙a(t),x˙d(t),y˙d(t))${\bf{Q}}(t) = \left( {{x_{\rm{d}}}(t),{x_{\rm{a}}}(t),{y_{\rm{d}}}(t),{{\dot x}_{\rm{a}}}(t),{{\dot x}_{\rm{d}}}(t),{{\dot y}_{\rm{d}}}(t)} \right)$ in the 2D-RTBP. This solution is periodic with a period T when it satisfies the periodic conditions x˙a(T)=x˙a(0)=0,xa(T)=xa(0),xd(T)=xd(0),        yd(T)=yd(0),x˙d(T)=x˙d(0),        y˙d(T)=y˙d(0).$\matrix{ {{{\dot x}_{\rm{a}}}(T) = {{\dot x}_{\rm{a}}}(0) = 0,} \hfill \cr {{x_{\rm{a}}}(T) = {x_{\rm{a}}}(0),} \hfill \cr {{x_{\rm{d}}}(T) = {x_{\rm{d}}}(0),\,\,\,\,\,\,\,\,{y_{\rm{d}}}(T) = {y_{\rm{d}}}(0),} \hfill \cr {{{\dot x}_{\rm{d}}}(T) = {{\dot x}_{\rm{d}}}(0),\,\,\,\,\,\,\,\,{{\dot y}_{\rm{d}}}(T) = {{\dot y}_{\rm{d}}}(0).} \hfill \cr } $(1)

Given the Poincaré surface of section at yd = 0 and the fact that xa=1−mamWD+ma${x_{\rm{a}}} = 1 - {{{m_{\rm{a}}}} \over {{m_{WD}} + {m_{\rm{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 { (xd(0),y˙d(0)) }$\left\{ {\left( {{x_{\rm{d}}}(0),{{\dot y}_{\rm{d}}}(0)} \right)} \right\}$. The asymmetric periodic orbit can for example be defined in the space { (xd(0),yd(0),x˙d(0),y˙d(0)) }$\left\{ {\left( {{x_{\rm{d}}}(0),{y_{\rm{d}}}(0),{{\dot x}_{\rm{d}}}(0),{{\dot y}_{\rm{d}}}(0)} \right)} \right\}$ when ẋa(0) = 0.

For the 3D periodic orbits shown in the right column of Fig. 2, the xz-symmetric periodic orbits (top panel) can be represented by a point in the 3D space of the initial conditions { (xd(0),zd(0),y˙d(0)) }$\left\{ {\left( {{x_{\rm{d}}}(0),{z_{\rm{d}}}(0),{{\dot y}_{\rm{d}}}(0)} \right)} \right\}$ and the x-symmetric ones (bottom panel) by { (xd(0),y˙d(0),z˙d(0)) }$\left\{ {\left( {{x_{\rm{d}}}(0),{{\dot y}_{\rm{d}}}(0),{{\dot z}_{\rm{d}}}(0)} \right)} \right\}$.

The families are formed during the mono-parametric continuation when a parameter, for example, zd, 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 ma and keeping zd fixed, for instance, while the respective periodicity condition still holds. The limitations of this continuation method can be found in Hadjidemetriou (1975); Antoniadou & Voyatzis (2013). The periodicity conditions that the periodic orbits must fulfil were described in greater detail for instance in Antoniadou & Libert (2018a) for the 2D-RTBPs and in Antoniadou & Libert (2019) for the 3D-RTBPs.

To numerically computate the periodic orbits we considered, we located without loss of generality the asteroid at aa = 1, and its period was therefore T0 = 2π. We also normalised the masses and the gravitational constant, G, so that the sum of the former was equal to unity, namely mWD + mа + md = 1 or mWD = 1−mа and G = 1.

In physical units, given the WD 1145+017 system, we adopted a WD mass mWD = 0.6m⊙. Then, we assumed an asteroid with a mass equal to approximately 1 : 10th the mass of Ceres, for instance, 2.5 × 1019 kg, located at aa = 0.0054 au around the WD (see e.g. Table 1 in Rappaport et al. 2016). 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 Lidov & Vashkov’yak (1994) for the specific case of retrograde orbits in the 2D-CRTBP and by Llibre & Ollé (2001) for the Saturn-Janus-Epimetheus system, respectively. Pousse et al. (2017) 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 Liberato & Winter (2020).

3.4 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 Voyatzis (2008). We constructed 200 × 100 grid planes and chose a maximum integration time for the computation of the DFLI equal to tmax = 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 tmax 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 (Antoniadou et al. 2014b).

thumbnail Fig. 3

Families of QS and HS orbits in the 2D-CRTBP presented on the (xd, ed) plane, where xd is the position of the periodic orbit (dust particle) in the rotating frame of reference for µ = 2.1186 × 10−11. Families denoted by I correspond to family ƒ, which consists of QS orbits, namely stable (blue) symmetric periodic orbits with ∆ϖ = 180°, and II denotes the family of HS orbits, i.e. unstable (red) symmetric periodic orbits with ∆ϖ = 0°. Primed (or unprimed) symbols indicate the direction (positive or negative crossing) of the orbit on the Poincaré surface of section at yd = 0 and y˙d>0${\dot y_{\rm{d}}} > 0$ (or y˙d<0${\dot y_{\rm{d}}} < 0$), since the two branches of each family are equivalent (see also Table 1). The magenta and green dots correspond to the vcos, and the solid (dashed) lines demonstrate the vertical stability (instability).

4 Families of periodic orbits in the CRTBP

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

4.1 2D-CRTBP

In Fig. 3, we present the families of symmetric periodic orbits in the 2D-CRTBP on the (xd, ed) plane, namely the QS and HS orbits. In Fig. 4, we present the families I (stable branch) and II (unstable branch) on the (ad,ed) 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 y˙d>0${\dot y_{\rm{d}}} > 0$ or y˙d<0${\dot y_{\rm{d}}} < 0$) of the orbit for which the Poincaré surface of section at yd = 0 was chosen. More particularly, we obtain the configurations shown in Table 1, 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. 5, we show the family of asymmetric periodic orbits, which consists of TP orbits, and starts from L4 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 L5 at θ = 300°. These families are both horizontally and vertically stable (solid blue lines).

thumbnail Fig. 4

Families I of QS orbits (panel a) and II of HS orbits (panel b) on the (ad, ed) plane for our model (µ = 2.1186 × 10−11 presented in Fig. 3), together with a comparison with µ = 0.001 and families f and b to highlight the proximity of this model to the unperturbed case.
Table 1

Equivalent configurations of 1:1 resonant periodic orbits in the 2D-CRTBP presented in Fig. 3, namely I with I′ and II with II′, at t = 0 and at t = T/2.

thumbnail Fig. 5

Families of QS, HS, and TP orbits on the (ed, θ) plane. Presentation as in Fig. 3.
thumbnail Fig. 6

Three families, FI, FII′, and GII′ of 3D periodic orbits in the 3D-CRTBP on the (ed, id) plane. The families of QS and HS orbits of the 2D-CRTBP lie on the x-axis at id = 0°. Presentation as in Fig. 3.
thumbnail Fig. 7

Three families of 3D periodic orbits in the 3D-CRTBP of Fig. 6 shown in (xd,y˙d,zd)$\left( {{x_{\rm{d}}},{{\dot y}_{\rm{d}}},{z_{\rm{d}}}} \right)$ space and (xd,y˙d,z˙d)$\left( {{x_{\rm{d}}},{{\dot y}_{\rm{d}}},{{\dot z}_{\rm{d}}}} \right)$ space. The respective families of the 2D-CRTBP together with the bifurcation points (vcos) are also included (at zd = 0 or żd = 0). Presentation as in Fig. 3.
Table 2

Bifurcation points (vcos) from the 2D-CRTBP to the 3D-CRTBP shown in Fig. 3.

4.2 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 (Hénon 1973; Ichtiaroglou & Michalodimitrakis 1980), 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${{{\hat F}^I}}$ belongs to the 2D family and generates the 3D family FI of xz-symmetric (magenta vco) periodic orbits. In Table 2, 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 FI 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 (ma = md = 0) (see the vcos (magenta dots) in panel a of Fig. 4 and compare e.g. the F^I${{{\hat F}^I}}$ (or F^I′${{{\hat F}^{I'}}}$) with the Bcs points in Table 1 in Voyatzis & Antoniadou 2018). 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 ed > 0.18 (see the top panel of Fig. 4).

In Fig. 6, we present the three 3D families of periodic orbits in the 3D-CRTBP on the (ed, id) plane. The family F’ consists of stable xz-symmetric periodic orbits (QS orbits) up to id,max = 30.7°. Both FII′ and GII′ have only unstable periodic orbits (HS orbits) with an xz- and x-symmetry, respectively. In Table 3, we present the configuration of each family. In Fig. 7, each family is shown in more detail in (xd,y˙d,zd)$\left( {{x_{\rm{d}}},{{\dot y}_{\rm{d}}},{z_{\rm{d}}}} \right)$ space or (xd,y˙d,z˙d)$\left( {{x_{\rm{d}}},{{\dot y}_{\rm{d}}},{{\dot z}_{\rm{d}}}} \right)$ 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 zd = 0 or żd = 0).

Table 3

Configuration of the dust particle for each family in the 3D-CRTBP presented in Fig. 6.

5 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 TP4 (or TP5), 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 (see also Namouni 1999; Namouni et al. 1999; Nesvorný et al. 2002; Sidorenko et al. 2014, for more details on the interchange of such regions as the problem’s parameters vary for a different value of µ). The breakdown of the phase space was organised as follows.

First, we chose three 3D periodic orbits (one from each family, namely FI, FII′ and GII′) with mutual inclinations, ∆i, equal to 10°, 30°, and 50° (Figs. 8-10). Then, we varied some of the initial conditions to create the grids on the (ed, id), (Ωd, ed), and (Md, ed) plane, while keeping the rest fixed. The values of the angles are given in Table 3, and the orbital elements of the selected periodic orbits are also provided in each figure.

In Fig. 8, 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 FI family with ad = 1 + 1.887 × 10−11, ed = 0.718, and ∆Ω = 270°. The areas around the stable 3D orbit are populated with regular orbits, and chaoticity becomes apparent for ed > 0.7 in the top left panel, where the planar family, 1, becomes vertically unstable (dashed line in Fig. 3). In the middle and bottom rows, the periodic orbits that guided the search are unstable and belong to the FII′ and GII′ families, respectively. The xz-symmetric family has initial conditions ad = 1−3.4903 × 10−11, ed = 0.722, and ∆Ω = 90°, and the x-symmetric family has ad = 1–2.7865 × 10−11, ed = 0.923, and ∆Ω = 0°. The regions around both of them are dominated by irregular orbits. However, the remaining regular domains QS and TP4 (or TP5) become apparent farther away from this periodic orbit with changing angles.

For Fig. 9, 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, ed) and (Md, ed) plane, as the variation in the (ed, id) plane was qualitatively similar with that presented in the left column of Fig. 8. In the left column of Fig. 9, the stable 3D periodic orbit has ad = 1 + 2.366 × 10−11, ed = 0.881, and ∆Ω = 270°. This 3D periodic orbit is shown in panel a of Fig. 2. In the middle column, the 3D unstable orbit has ad = 1–3.291 × 10−11, ed = 0.761, and ∆Ω = 90°, and in the right column, the initial conditions are ad = 1–2.449 × 10−11, ed = 0.904, and ∆Ω = 0°.

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

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

In Fig. 12, we varied the semi-major axis and the eccentricity of the dust particle on the (ad, ed) plane. In the top and middle rows, we used the same six periodic orbits as in Figs. 8 and 9 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 L4 (or L5) 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 Nesvorný et al. (2002).

thumbnail Fig. 8

DS maps for ∆i = 10° guided by a stable periodic orbit of the FI family (top) and an unstable periodic orbit of the FII′ family (middle) and GII′ family (bottom). The chosen initial conditions, which remained fixed while some varied for the computation of each grid, are reported above each row (see also Table 3). The dashed cyan curves depict the collisions and close encounters between the bodies. In the areas denoted by QS and HS, θ librates about 0 and 180°, respectively. TP4 and TP5 stand for the regular domains created by the asymmetric families generated by L4 and L5, where θ librates accordingly. R denotes areas in which θ rotates. Dark (pale) colours show regular (chaotic) orbits.
thumbnail Fig. 9

DS maps for ∆i = 30° guided by a stable periodic orbit of the FI family (left) and an unstable periodic orbit of the FII′ family (middle) and GII′ family (right). The chosen initial conditions, which remained fixed while some varied for the computation of each grid, are reported above each column (see also Table 3). Presentation as in Fig. 8.
thumbnail Fig. 10

DS maps for ∆i = 50° guided by an unstable periodic orbit of the FII′ family (left) and GII′ family (right). Presentation as in Fig. 9.
thumbnail Fig. 11

DS maps on the (Ωd, Md) plane around asymmetric periodic orbits with ed = 0 (left), ed = 0.5 (middle), and ed = 0.9 (right).
thumbnail Fig. 12

DS maps on the (ad, ed) plane guided by symmetric periodic orbits with ∆i = 10° (top) and ∆i = 30° (middle) and by asymmetric configurations with ∆i = 10° (bottom).
thumbnail Fig. 13

Simulations on the (∆a, ed) plane (top row) and in the (∆a, ed, ∆Pmаx) space (bottom row) guided by the 3D QS periodic orbit with ∆Ω = 270° and ∆i = 10° that formulated the respective DS map in Fig. 12 (top left panel). Here, ad = aа ∓ 68.135 km, and Md is uniformly sampled and incorporates orbits around the 3D QS periodic orbit, and the left column represents a magnification for ed < 0.1. Unstable systems are indicated by red crosses on the (∆a, ed) plane.

6 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 (Antoniadou & Veras 2016, 2019).

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 Antoniadou & Veras (2016) for the 2D-CRTBP and by Antoniadou & Veras (2019) 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 ai(N)=aiζ1/3,$a_i^{({\rm{N}})} = {a_i}{\zeta ^{1/3}},$(2)

where ζ≡mWD+ma1m⊙$\zeta \equiv {{{m_{WD}} + {m_{\rm{a}}}} \over {1{m_ \odot }}}$ 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 (Rein & Liu 2012; Rein & Spiegel 2015). 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. 8. These forces typically have a greater effect the smaller the orbital pericentre. In this context, the most consequential results are those with ed ≈ 0.0, whereas those with more limited applicability have ed ≳ 0.8. Nevertheless, we sampled ed values up to 0.9 in many cases to provide a direct comparison with the DS maps.

7 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 (Gänsicke et al. 2016; Farihi et al. 2022; Aungwerojwit et al. 2024).

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 Md. In some cases, in order to provide a direct comparison to the DS maps, we adopted the value of Md that corresponds to a particular family, or DS map. However, in some other cases, we sampled a range of Md values in order to determine the maximum period deviation independent of Md, and we report only those that generated the maximum orbital period deviation.

7.1 Variation in ed versus ad

Fig. 13 shows a 91 × 91 grid of ad = aa ∓ 68.135 km and ed = 0.0001–0.1000 (left column) and ed = 0.0001–0.9000 (right column). For each combination of (ad, ed), we performed 18 simulations with Md 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. 14, 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. 15 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. 16 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 Md = 0°. Hence, the value of the maximum orbital period deviation observed in the region populated by TP orbits was lower.

thumbnail Fig. 14

Presentation as in Fig. 13, but guided by the 3D HS periodic orbit with ∆Ω = 90° and ∆i = 10° that formulated the respective DS map in Fig. 12 (top middle panel). Unstable systems are indicated by red crosses on the (∆a, ed) plane.

7.2 Variation in ed versus ωd

Fig. 17 shows a 91 × 91 grid of ωd = 0° − 360° and ed = 0.0001–0.9000 guided by the 3D QS periodic orbit with ∆Ω = 270° and ∆i = 10°. For each combination of (ωd, ed), we performed a single simulation with Md = 180° (left column) and 18 simulations with Md uniformly sampled across its entire range (right column). When Md = 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 Md.

thumbnail Fig. 15

Presentation as in Fig. 13, but guided by the 3D HS periodic orbit with ∆Ω = 0° and ∆i = 10° that formulated the respective DS map in Fig. 12 (top right panel). Unstable systems are indicated by red crosses on the (∆a, ed) plane.

7.3 Variation in id versus ed

Fig. 18 shows a 181 × 181 grid of ed = 0.0001 – 0.8000 and id = 0°–90° guided by the HS periodic orbit with Ωd = 90° and ∆i = 50°, and we performed a single simulation with Md = 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.

8 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.

thumbnail Fig. 16

Presentation as in Fig. 13, but guided by a TP periodic orbit with Ωd = 60° and ∆i = 10° that formulated the respective DS map in Fig. 12 (bottom left panel). Md is not uniformly sampled, but equals 0°. Unstable systems are indicated by red crosses on the (∆a, ed) plane.

8.1 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 (Li 2012; Veras 2014). 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(ad1.00 R⊙)5/2(1−ed2).$107{\rm{yr}}{\left( {{{{M_ \star }} \over {0.60{M_ \odot }}}} \right)^{ - 3/2}}{\left( {{{{a_{\rm{d}}}} \over {1.00{R_ \odot }}}} \right)^{5/2}}\left( {1 - e_{\rm{d}}^2} \right).$(3)

As a result, we might imagine a dust particle traversing a horizontal path across the plots in Figs. 8-11. 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.

thumbnail Fig. 17

Simulations on the (Ωd, ed) plane (top row) and in the (Ωd, ed, ∆Pmax) space (bottom row) guided by the 3D QS periodic orbit with ∆Ω = 270° and ∆i = 10° that formulated the respective DS map in Fig. 8 (top middle panel). In the left column, Md = 180°, while it is uniformly sampled and incorporates orbits around the 3D QS periodic orbit on the right.

8.2 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 (Steckloff et al. 2021), 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 (Elms et al. 2022). Old WDs hence allow the orbital pericentre of a dust particle to be far inwards of 1 R⊙ without sublimating.

In order to quantify this effect, we computed the maximum value of ed 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 Veras et al. (2022), and we report our results in Fig. 19. 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 ed ≈ 0.84–0.9 may realistically be adopted.

thumbnail Fig. 18

Simulations on the (id, ed) plane (top row) and in the (id, ed, ∆Pmax) space (bottom row) guided by a HS periodic orbit with Ωd = 90° and ∆i = 50° that formulated the respective DS maps in Fig. 10 (left column). Md is not uniformly sampled, but instead equals 180° in each simulation.

8.3 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 Wyatt & Whipple (1950) formula for the semi-major axis variation due to Poynting-Robertson drag, dadt≈−3L⋆(2+3ed2)16πc2Rdρdad(1−ed2)3/2,${{da} \over {dt}} \approx - {{3{L_ \star }\left( {2 + 3e_{\rm{d}}^2} \right)} \over {16\pi {c^2}{R_{\rm{d}}}{\rho _{\rm{d}}}{a_{\rm{d}}}{{\left( {1 - e_{\rm{d}}^2} \right)}^{3/2}}}},$(4)

and we assumed that the particle density is ρd = 2 g/cm3.

Figure 20 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 ed ≈ 0, based on the x-axes of the N-body simulation plots for Figs. 13–16, 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.

thumbnail Fig. 19

Maximum eccentricity values of the dust particle (right y-axis), consisting of three different materials, before the particles are sublimated, as a function of the WD cooling age. Stronger materials may reach ed ≈ 0.8–0.9 without sublimating when the WD is several billion years old.
thumbnail Fig. 20

Inward drift rate from Poynting–Robertson drag as a function of the WD cooling age. Even the largest boulders may be dragged enough to trigger gravitational close encounters with the co-orbital asteroid.

9 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 ed → 0. Nonetheless, we found that the planar families exhibited qualitatively similar attributes with those computed for µ = 0.001 by Pousse et al. (2017) and Voyatzis & Antoniadou (2018), 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 ed < 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 (ad ≠ 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 (see e.g. the stable domains in the DS maps in Antoniadou & Libert 2018a; Anto-niadou & Libert 2018b; Antoniadou & Libert 2019, for possible detection of terrestrial planets around single-planet systems, i.e. the case of 2D and 3D RTBPs with µ = 0.001)2. 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.

Acknowledgements

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.

References

  1. Antoniadou, K. I., & Libert, A.-S. 2018a, A&A, 615, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Antoniadou, K. I., & Libert, A.-S. 2018b, Celest. Mech. Dyn. Astron., 130, 41 [NASA ADS] [CrossRef] [Google Scholar]
  3. Antoniadou, K. I., & Libert, A.-S. 2019, MNRAS, 483, 2923 [NASA ADS] [CrossRef] [Google Scholar]
  4. Antoniadou, K. I., & Libert, A.-S. 2020, A&A, 640, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Antoniadou, K. I., & Veras, D. 2016, MNRAS, 463, 4108 [NASA ADS] [CrossRef] [Google Scholar]
  6. Antoniadou, K. I., & Veras, D. 2019, A&A, 629, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Antoniadou, K. I., & Voyatzis, G. 2013, Celest. Mech. Dyn. Astron., 115, 161 [Google Scholar]
  8. Antoniadou, K. I., & Voyatzis, G. 2014, Ap&SS, 349, 657 [NASA ADS] [CrossRef] [Google Scholar]
  9. Antoniadou, K. I., & Voyatzis, G. 2016, MNRAS, 461, 3822 [NASA ADS] [CrossRef] [Google Scholar]
  10. Antoniadou, K. I., & Voyatzis, G. 2022, A&A, 661, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Antoniadou, K. I., Voyatzis, G., & Varvoglis, H. 2014a, Proc. Int. Astron. Union, 9, 82 [CrossRef] [Google Scholar]
  12. Antoniadou, K. I., Voyatzis, G., & Varvoglis, H. 2014b, in NumAn2014 Conference on Numerical Analysis, Proceedings of the 6th International Conference on Numerical Analysis (NumAn2014) (AMCL/TUC), 7 [Google Scholar]
  13. Aungwerojwit, A., Gänsicke, B. T., Dhillon, V. S., et al. 2024, MNRAS, 530, 117 [NASA ADS] [CrossRef] [Google Scholar]
  14. Ballering, N. P., Levens, C. I., Su, K. Y. L., & Cleeves, L. I. 2022, ApJ, 939, 108 [NASA ADS] [CrossRef] [Google Scholar]
  15. Balsalobre-Ruza, O., de Gregorio-Monsalvo, I., Lillo-Box, J., et al. 2023, A&A, 675, A172 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Benest, D. 1974, A&A, 32, 39 [NASA ADS] [Google Scholar]
  17. Broucke, R. A. 1968, Periodic Orbits in the Restricted Three-body Problem with Earth-Moon Masses (NASA Jet Propulsion Laboratory) [Google Scholar]
  18. Brouwers, M. G., Bonsor, A., & Malamud, U. 2022, MNRAS, 509, 2404 [NASA ADS] [Google Scholar]
  19. Brouwers, M. G., Bonsor, A., & Malamud, U. 2023, MNRAS, 519, 2646 [Google Scholar]
  20. Brown, E. W. 1911, MNRAS, 71, 438 [NASA ADS] [Google Scholar]
  21. Bruno, A. D. 1994, The Restricted 3-Body Problem: Plane Periodic Orbits (Berlin, New York: Walter de Gruyter) [CrossRef] [Google Scholar]
  22. Christou, A. A. 2000a, Icarus, 144, 1 [NASA ADS] [CrossRef] [Google Scholar]
  23. Christou, A. A. 2000b, A&A, 356, L71 [NASA ADS] [Google Scholar]
  24. Connors, M. 2014, MNRAS, 437, L85 [Google Scholar]
  25. Connors, M., Veillet, C., Brasser, R., et al. 2004, Meteor. Planet. Sci., 39, 1251 [NASA ADS] [CrossRef] [Google Scholar]
  26. Couturier, J., Robutel, P., & Correia, A. C. M. 2022, A&A, 664, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Darwin, G. H. 1897, Acta Math., XXI, 101 [Google Scholar]
  28. de la Fuente Marcos, C., & de la Fuente Marcos, R. 2014, MNRAS, 445, 2985 [CrossRef] [Google Scholar]
  29. de la Fuente Marcos, C., & de la Fuente Marcos, R. 2016, MNRAS, 462, 3441 [CrossRef] [Google Scholar]
  30. de la Fuente Marcos, C., & de la Fuente Marcos, R. 2018, MNRAS, 473, 3434 [Google Scholar]
  31. de la Fuente Marcos, C., & de la Fuente Marcos, R. 2021, MNRAS, 501, 6007 [NASA ADS] [CrossRef] [Google Scholar]
  32. Debes, J. H., Walsh, K. J., & Stark, C. 2012, ApJ, 747, 148 [NASA ADS] [CrossRef] [Google Scholar]
  33. Deprit, A., Jacques, H., Palmore, J., & Price, J. F. 1967, MNRAS, 137, 311 [NASA ADS] [CrossRef] [Google Scholar]
  34. Duvvuri, G. M., Redfield, S., & Veras, D. 2020, ApJ, 893, 166 [NASA ADS] [CrossRef] [Google Scholar]
  35. Elms, A. K., Tremblay, P.-E., Gänsicke, B. T., et al. 2022, MNRAS, 517, 4557 [NASA ADS] [CrossRef] [Google Scholar]
  36. Farihi, J., Hermes, J. J., Marsh, T. R., et al. 2022, MNRAS, 511, 1647 [NASA ADS] [CrossRef] [Google Scholar]
  37. Funk, B., Dvorak, R., & Schwarz, R. 2013, Celest. Mech. Dyn. Astron., 117, 41 [NASA ADS] [CrossRef] [Google Scholar]
  38. Gänsicke, B. T., Aungwerojwit, A., Marsh, T. R., et al. 2016, ApJ, 818, L7 [CrossRef] [Google Scholar]
  39. Gary, B. L., Rappaport, S., Kaye, T. G., Alonso, R., & Hambschs, F. J. 2017, MNRAS, 465, 3267 [Google Scholar]
  40. Girven, J., Brinkworth, C. S., Farihi, J., et al. 2012, ApJ, 749, 154 [NASA ADS] [CrossRef] [Google Scholar]
  41. Giuppone, C. A., Beaugé, C., Michtchenko, T. A., & Ferraz-Mello, S. 2010, MNRAS, 407, 390 [NASA ADS] [CrossRef] [Google Scholar]
  42. Giuppone, C. A., Benítez-Llambay, P., & Beaugé, C. 2012, MNRAS, 421, 356 [Google Scholar]
  43. Goksu, O. A., Kutra, T., & Wu, Y. 2024, AJ, 167, 236 [NASA ADS] [CrossRef] [Google Scholar]
  44. Guidry, J. A., Vanderbosch, Z. P., Hermes, J. J., et al. 2021, ApJ, 912, 125 [NASA ADS] [CrossRef] [Google Scholar]
  45. Gurri, P., Veras, D., & Gänsicke, B. T. 2017, MNRAS, 464, 321 [CrossRef] [Google Scholar]
  46. Hadjidemetriou, J. D. 1975, Celest. Mech., 12, 155 [NASA ADS] [CrossRef] [Google Scholar]
  47. Hadjidemetriou, J. D. 1993, Celest. Mech. Dyn. Astron., 56, 563 [NASA ADS] [CrossRef] [Google Scholar]
  48. Hadjidemetriou, J. D., & Christides, T. 1975, Celest. Mech., 12, 175 [NASA ADS] [CrossRef] [Google Scholar]
  49. Hadjidemetriou, J. D., & Voyatzis, G. 2011, Celest. Mech. Dyn. Astron., 111, 179 [NASA ADS] [CrossRef] [Google Scholar]
  50. Hadjidemetriou, J. D., Psychoyos, D., & Voyatzis, G. 2009, Celest. Mech. Dyn. Astron., 104, 23 [NASA ADS] [CrossRef] [Google Scholar]
  51. Hénon, M. 1969, A&A, 1, 223 [Google Scholar]
  52. Hénon, M. 1973, A&A, 28, 415 [NASA ADS] [Google Scholar]
  53. Hénon, M. 1997, Lecture Notes in Physics Monographs, 52, Generating Families in the Restricted Three-Body Problem (Springer-Verlag Berlin Heidelberg) [Google Scholar]
  54. Hénon, M., & Guyot, M. 1970, in Periodic Orbits Stability and Resonances, ed. G. E. O. Giacaglia, 349 [CrossRef] [Google Scholar]
  55. Hill, G. W. 1878, Am. J. Math., 1, 5 [Google Scholar]
  56. Ichtiaroglou, S., & Michalodimitrakis, M. 1980, A&A, 81, 30 [NASA ADS] [Google Scholar]
  57. Izquierdo, P., Rodríguez-Gil, P., Gänsicke, B. T., et al. 2018, MNRAS, 481, 703 [Google Scholar]
  58. Jackson, J. 1913, MNRAS, 74, 62 [NASA ADS] [Google Scholar]
  59. Kimura, M., Kawamura, M., & Yamada, K. 2019, Adv. Space Res., 63, 1336 [NASA ADS] [CrossRef] [Google Scholar]
  60. Kinoshita, H., & Nakai, H. 2007, Celest. Mech. Dyn. Astron., 98, 181 [NASA ADS] [CrossRef] [Google Scholar]
  61. Kogan, A. Y. 1988a, in Pioneers of Space and the Present Age (Moscow: Izdatel’stvo Nauka), 151 [Google Scholar]
  62. Kogan, A. Y. 1988b, Kosmich. Issled., 26, 813 [NASA ADS] [Google Scholar]
  63. Kogan, A. Y. 1990, in Quasi-satellite Orbits and their Applications, Dresden International Astronautical Federation Congress [Google Scholar]
  64. Kortenkamp, S. J. 2013, Icarus, 226, 1550 [CrossRef] [Google Scholar]
  65. Laughlin, G., & Chambers, J. E. 2002, AJ, 124, 592 [Google Scholar]
  66. Leleu, A., Robutel, P., Correia, A. C. M., & Lillo-Box, J. 2017, A&A, 599, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Leleu, A., Robutel, P., & Correia, A. C. M. 2018, Celest. Mech. Dyn. Astron., 130, 24 [NASA ADS] [CrossRef] [Google Scholar]
  68. Leleu, A., Lillo-Box, J., Sestovic, M., et al. 2019, A&A, 624, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Li, L.-S. 2012, Astrophys. Space Sci., 341, 323 [NASA ADS] [CrossRef] [Google Scholar]
  70. Li, D., Mustill, A. J., & Davies, M. B. 2021, MNRAS, 508, 5671 [NASA ADS] [CrossRef] [Google Scholar]
  71. Li, Q., Tao, Y., & Jiang, F. 2022, Aerospace, 9, 454 [NASA ADS] [CrossRef] [Google Scholar]
  72. Liberato, L., & Winter, O. C. 2020, MNRAS, 496, 3700 [NASA ADS] [CrossRef] [Google Scholar]
  73. Lidov, M. L., & Vashkov’yak, M. A. 1994, Astron. Lett., 20, 676 [NASA ADS] [Google Scholar]
  74. Lillo-Box, J., Barrado, D., Figueira, P., et al. 2018, A&A, 609, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Llibre, J., & Ollé, M. 2001, A&A, 378, 1087 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Madeira, G., Charnoz, S., Rambaux, N., & Robutel, P. 2024, Icarus, 412, 115997 [NASA ADS] [CrossRef] [Google Scholar]
  77. Malamud, U., & Perets, H. B. 2020a, MNRAS, 492, 5561 [Google Scholar]
  78. Malamud, U., & Perets, H. B. 2020b, MNRAS, 493, 698 [NASA ADS] [CrossRef] [Google Scholar]
  79. Malamud, U., Grishin, E., & Brouwers, M. 2021, MNRAS, 501, 3806 [NASA ADS] [Google Scholar]
  80. Manser, C. J., Gänsicke, B. T., Marsh, T. R., et al. 2016, MNRAS, 455, 4467 [Google Scholar]
  81. Manser, C. J., Dennihy, E., Gänsicke, B. T., et al. 2021, MNRAS, 508, 5657 [NASA ADS] [CrossRef] [Google Scholar]
  82. Marchal, C. 1990, The Three-body Problem (Amsterdam: Elsevier) [Google Scholar]
  83. Mikkola, S., & Innanen, K. 1997, in The Dynamical Behaviour of our Planetary System, eds. R. Dvorak, & J. Henrard (Kluwer Academic Publishers), 345 [CrossRef] [Google Scholar]
  84. Mikkola, S., Brasser, R., Wiegert, P., & Innanen, K. 2004, MNRAS, 351, L63 [CrossRef] [Google Scholar]
  85. Mikkola, S., Innanen, K., Wiegert, P., Connors, M., & Brasser, R. 2006, MNRAS, 369, 15 [NASA ADS] [CrossRef] [Google Scholar]
  86. Minghu, T., Ke, Z., L., M., & X., C. 2014, Acta Astronaut., 98, 50 [NASA ADS] [CrossRef] [Google Scholar]
  87. Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge University Press) [Google Scholar]
  88. Namouni, F. 1999, Icarus, 137, 293 [NASA ADS] [CrossRef] [Google Scholar]
  89. Namouni, F., & Morais, M. 2018, MNRAS, 477, L117 [NASA ADS] [CrossRef] [Google Scholar]
  90. Namouni, F., Christou, A. A., & Murray, C. D. 1999, Phys. Rev. Lett., 83, 2506 [NASA ADS] [CrossRef] [Google Scholar]
  91. Nesvorný, D., Thomas, F., Ferraz-Mello, S., & Morbidelli, A. 2002, Celest. Mech. Dyn. Astron., 82, 323 [CrossRef] [Google Scholar]
  92. Niederman, L., Pousse, A., & Robutel, P. 2020, Commun. Math. Phys., 377, 551 [NASA ADS] [CrossRef] [Google Scholar]
  93. Nixon, C. J., Pringle, J. E., Coughlin, E. R., Swan, A., & Farihi, J. 2020, arXiv e-prints [arXiv:2006.07639] [Google Scholar]
  94. O’Connor, C. E., & Lai, D. 2020, MNRAS, 498, 4005 [CrossRef] [Google Scholar]
  95. Perozzi, E., Ceccaroni, M., & B., V. G. 2017, Eur. Phys. J. Plus, 132, 367 [NASA ADS] [CrossRef] [Google Scholar]
  96. Pires, P., & Winter, O. C. 2020, MNRAS, 494, 2727 [CrossRef] [Google Scholar]
  97. Pousse, A., Robutel, P., & Vienne, A. 2017, Celest. Mech. Dyn. Astron., 128, 383 [NASA ADS] [CrossRef] [Google Scholar]
  98. Rappaport, S., Gary, B. L., Kaye, T., et al. 2016, MNRAS, 458, 3904 [Google Scholar]
  99. Rappaport, S., Gary, B. L., Vanderburg, A., et al. 2018, MNRAS, 474, 933 [NASA ADS] [CrossRef] [Google Scholar]
  100. Raymond, S. N., Veras, D., Clement, M. S., et al. 2023a, MNRAS, 521, 2002 [NASA ADS] [CrossRef] [Google Scholar]
  101. Raymond, S. N., Veras, D., Clement, M. S., et al. 2023b, MNRAS, 522, 4875 [NASA ADS] [CrossRef] [Google Scholar]
  102. Rein, H., & Liu, S. F. 2012, A&A, 537, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  103. Rein, H., & Spiegel, D. S. 2015, MNRAS, 446, 1424 [Google Scholar]
  104. Robutel, P., & Pousse, A. 2013, Celest. Mech. Dyn. Astron., 117, 17 [NASA ADS] [CrossRef] [Google Scholar]
  105. Sagdeev, R. Z., & Zakharov, A. V. 1989, Nature, 341, 581 [NASA ADS] [CrossRef] [Google Scholar]
  106. Schwarz, R., Süli, Á., Dvorak, R., & Pilat-Lohinger, E. 2009, Celest. Mech. Dyn. Astron., 104, 69 [NASA ADS] [CrossRef] [Google Scholar]
  107. Sidorenko, V. 2024, Celest. Mech. Dyn. Astron., 136, 25 [NASA ADS] [CrossRef] [Google Scholar]
  108. Sidorenko, V. V., Neishtadt, A. I., Artemyev, A. V., & Zelenyi, L. M. 2014, Celest. Mech. Dyn. Astron., 120, 131 [NASA ADS] [CrossRef] [Google Scholar]
  109. Steckloff, J. K., Debes, J., Steele, A., et al. 2021, ApJ, 913, L31 [NASA ADS] [CrossRef] [Google Scholar]
  110. Szebehely, V. 1967, Theory of Orbits. The Restricted Problem of Three Bodies (Academic Press) [Google Scholar]
  111. Vanderbosch, Z., Hermes, J. J., Dennihy, E., et al. 2020, ApJ, 897, 171 [Google Scholar]
  112. Vanderbosch, Z. P., Rappaport, S., Guidry, J. A., et al. 2021, ApJ, 917, 41 [NASA ADS] [CrossRef] [Google Scholar]
  113. Vanderburg, A., Johnson, J. A., Rappaport, S., et al. 2015, Nature, 526, 546 [Google Scholar]
  114. Veras, D. 2014, MNRAS, 442, L71 [NASA ADS] [CrossRef] [Google Scholar]
  115. Veras, D., & Heng, K. 2020, MNRAS, 496, 2292 [Google Scholar]
  116. Veras, D., Leinhardt, Z. M., Bonsor, A., & Gänsicke, B. T. 2014, MNRAS, 445, 2244 [NASA ADS] [CrossRef] [Google Scholar]
  117. Veras, D., Leinhardt, Z. M., Eggl, S., & Gänsicke, B. T. 2015, MNRAS, 451, 3453 [NASA ADS] [CrossRef] [Google Scholar]
  118. Veras, D., Marsh, T. R., & Gänsicke, B. T. 2016, MNRAS, 461, 1413 [NASA ADS] [CrossRef] [Google Scholar]
  119. Veras, D., Carter, P. J., Leinhardt, Z. M., & Gänsicke, B. T. 2017, MNRAS, 465, 1008 [NASA ADS] [CrossRef] [Google Scholar]
  120. Veras, D., Georgakarakos, N., Mustill, A. J., et al. 2021, MNRAS, 506, 1148 [NASA ADS] [CrossRef] [Google Scholar]
  121. Veras, D., Birader, Y., & Zaman, U. 2022, MNRAS, 510, 3379 [NASA ADS] [CrossRef] [Google Scholar]
  122. Voyatzis, G. 2008, ApJ, 675, 802 [Google Scholar]
  123. Voyatzis, G., & Antoniadou, K. I. 2018, Celest. Mech. Dyn. Astron., 130, 59 [NASA ADS] [CrossRef] [Google Scholar]
  124. Voyatzis, G., & Hadjidemetriou, J. D. 2005, Celest. Mech. Dyn. Astron., 93, 263 [Google Scholar]
  125. Wajer, P. 2009, Icarus, 200, 147 [NASA ADS] [CrossRef] [Google Scholar]
  126. Wajer, P. 2010, Icarus, 209, 488 [NASA ADS] [CrossRef] [Google Scholar]
  127. Wiegert, P., Innanen, K., & Mikkola, S. 2000, AJ, 119, 1978 [NASA ADS] [CrossRef] [Google Scholar]
  128. Williams, G. A., & Murray, C. D. 2011, Icarus, 212, 275 [NASA ADS] [CrossRef] [Google Scholar]
  129. Wyatt, S. P., & Whipple, F. L. 1950, ApJ, 111, 134 [NASA ADS] [CrossRef] [Google Scholar]
1

The period of short-period orbits tends to 2π, whereas the period of the long-period orbits tends to infinity when µ → 0.

2

Likewise for possible detections in the GTBP, (see e.g. Antoniadou & Voyatzis 2016; Antoniadou & Libert 2020; Antoniadou & Voyatzis 2022).

All Tables

Table 1

Equivalent configurations of 1:1 resonant periodic orbits in the 2D-CRTBP presented in Fig. 3, namely I with I′ and II with II′, at t = 0 and at t = T/2.

In the text
Table 2

Bifurcation points (vcos) from the 2D-CRTBP to the 3D-CRTBP shown in Fig. 3.

In the text
Table 3

Configuration of the dust particle for each family in the 3D-CRTBP presented in Fig. 6.

In the text

All Figures

thumbnail Fig. 1

Periodic orbits in the 2D-CRTBP, where the asteroid, denoted by A (in blue), moves on circular orbits (ea = 0), and the dust particle, D (in red), moves on Keplerian ellipses that rotate around the WD (grey dot), when viewed on an inertial frame of reference (left column). These orbits are exactly periodic in the rotating frame of reference (right column). The orbits were computed for one period, T = 2π, and both bodies orbit the WD with the same period (corresponding in physical space to a distance unit of about 1 R⊙ and to an orbital period of about 4.5 hours). The red dots represent the position of the dust particle at t = 0. In panel a, the stable retrograde satellite orbits that belong to family f are demonstrated, where the dust particle moves on almost circular orbit. In panels b and c, the symmetric stable QSh orbits of family f, and the symmetric unstable horseshoe (HS) orbits are shown, respectively. In panel d, the stable asymmetric tadpole (TP) periodic orbits are shown. The eccentricity of the dust particle chosen for these panels was ed = 0.7, and the initial conditions were taken from the families presented in Fig. 5.
In the text
thumbnail Fig. 2

Stable (QSh) and unstable (HS) inclined periodic orbits in the 3D-CRTBP with id = 30° (panel a) and id = 50° (panel b). These orbits belong to the 3D families FI with periodic orbits symmetric with respect to the xz-plane and GII′ with periodic orbits symmetric with respect to the x-axis, (both families are shown Fig. 6). These 3D families may either bifurcate from the 2D family I (or I′) or from the 2D family II (or II′) shown in Figs. 3-4. Presentation as in Fig. 1.
In the text
thumbnail Fig. 3

Families of QS and HS orbits in the 2D-CRTBP presented on the (xd, ed) plane, where xd is the position of the periodic orbit (dust particle) in the rotating frame of reference for µ = 2.1186 × 10−11. Families denoted by I correspond to family ƒ, which consists of QS orbits, namely stable (blue) symmetric periodic orbits with ∆ϖ = 180°, and II denotes the family of HS orbits, i.e. unstable (red) symmetric periodic orbits with ∆ϖ = 0°. Primed (or unprimed) symbols indicate the direction (positive or negative crossing) of the orbit on the Poincaré surface of section at yd = 0 and y˙d>0${\dot y_{\rm{d}}} > 0$ (or y˙d<0${\dot y_{\rm{d}}} < 0$), since the two branches of each family are equivalent (see also Table 1). The magenta and green dots correspond to the vcos, and the solid (dashed) lines demonstrate the vertical stability (instability).
In the text
thumbnail Fig. 4

Families I of QS orbits (panel a) and II of HS orbits (panel b) on the (ad, ed) plane for our model (µ = 2.1186 × 10−11 presented in Fig. 3), together with a comparison with µ = 0.001 and families f and b to highlight the proximity of this model to the unperturbed case.
In the text
thumbnail Fig. 5

Families of QS, HS, and TP orbits on the (ed, θ) plane. Presentation as in Fig. 3.
In the text
thumbnail Fig. 6

Three families, FI, FII′, and GII′ of 3D periodic orbits in the 3D-CRTBP on the (ed, id) plane. The families of QS and HS orbits of the 2D-CRTBP lie on the x-axis at id = 0°. Presentation as in Fig. 3.
In the text
thumbnail Fig. 7

Three families of 3D periodic orbits in the 3D-CRTBP of Fig. 6 shown in (xd,y˙d,zd)$\left( {{x_{\rm{d}}},{{\dot y}_{\rm{d}}},{z_{\rm{d}}}} \right)$ space and (xd,y˙d,z˙d)$\left( {{x_{\rm{d}}},{{\dot y}_{\rm{d}}},{{\dot z}_{\rm{d}}}} \right)$ space. The respective families of the 2D-CRTBP together with the bifurcation points (vcos) are also included (at zd = 0 or żd = 0). Presentation as in Fig. 3.
In the text
thumbnail Fig. 8

DS maps for ∆i = 10° guided by a stable periodic orbit of the FI family (top) and an unstable periodic orbit of the FII′ family (middle) and GII′ family (bottom). The chosen initial conditions, which remained fixed while some varied for the computation of each grid, are reported above each row (see also Table 3). The dashed cyan curves depict the collisions and close encounters between the bodies. In the areas denoted by QS and HS, θ librates about 0 and 180°, respectively. TP4 and TP5 stand for the regular domains created by the asymmetric families generated by L4 and L5, where θ librates accordingly. R denotes areas in which θ rotates. Dark (pale) colours show regular (chaotic) orbits.
In the text
thumbnail Fig. 9

DS maps for ∆i = 30° guided by a stable periodic orbit of the FI family (left) and an unstable periodic orbit of the FII′ family (middle) and GII′ family (right). The chosen initial conditions, which remained fixed while some varied for the computation of each grid, are reported above each column (see also Table 3). Presentation as in Fig. 8.
In the text
thumbnail Fig. 10

DS maps for ∆i = 50° guided by an unstable periodic orbit of the FII′ family (left) and GII′ family (right). Presentation as in Fig. 9.
In the text
thumbnail Fig. 11

DS maps on the (Ωd, Md) plane around asymmetric periodic orbits with ed = 0 (left), ed = 0.5 (middle), and ed = 0.9 (right).
In the text
thumbnail Fig. 12

DS maps on the (ad, ed) plane guided by symmetric periodic orbits with ∆i = 10° (top) and ∆i = 30° (middle) and by asymmetric configurations with ∆i = 10° (bottom).
In the text
thumbnail Fig. 13

Simulations on the (∆a, ed) plane (top row) and in the (∆a, ed, ∆Pmаx) space (bottom row) guided by the 3D QS periodic orbit with ∆Ω = 270° and ∆i = 10° that formulated the respective DS map in Fig. 12 (top left panel). Here, ad = aа ∓ 68.135 km, and Md is uniformly sampled and incorporates orbits around the 3D QS periodic orbit, and the left column represents a magnification for ed < 0.1. Unstable systems are indicated by red crosses on the (∆a, ed) plane.
In the text
thumbnail Fig. 14

Presentation as in Fig. 13, but guided by the 3D HS periodic orbit with ∆Ω = 90° and ∆i = 10° that formulated the respective DS map in Fig. 12 (top middle panel). Unstable systems are indicated by red crosses on the (∆a, ed) plane.
In the text
thumbnail Fig. 15

Presentation as in Fig. 13, but guided by the 3D HS periodic orbit with ∆Ω = 0° and ∆i = 10° that formulated the respective DS map in Fig. 12 (top right panel). Unstable systems are indicated by red crosses on the (∆a, ed) plane.
In the text
thumbnail Fig. 16

Presentation as in Fig. 13, but guided by a TP periodic orbit with Ωd = 60° and ∆i = 10° that formulated the respective DS map in Fig. 12 (bottom left panel). Md is not uniformly sampled, but equals 0°. Unstable systems are indicated by red crosses on the (∆a, ed) plane.
In the text
thumbnail Fig. 17

Simulations on the (Ωd, ed) plane (top row) and in the (Ωd, ed, ∆Pmax) space (bottom row) guided by the 3D QS periodic orbit with ∆Ω = 270° and ∆i = 10° that formulated the respective DS map in Fig. 8 (top middle panel). In the left column, Md = 180°, while it is uniformly sampled and incorporates orbits around the 3D QS periodic orbit on the right.
In the text
thumbnail Fig. 18

Simulations on the (id, ed) plane (top row) and in the (id, ed, ∆Pmax) space (bottom row) guided by a HS periodic orbit with Ωd = 90° and ∆i = 50° that formulated the respective DS maps in Fig. 10 (left column). Md is not uniformly sampled, but instead equals 180° in each simulation.
In the text
thumbnail Fig. 19

Maximum eccentricity values of the dust particle (right y-axis), consisting of three different materials, before the particles are sublimated, as a function of the WD cooling age. Stronger materials may reach ed ≈ 0.8–0.9 without sublimating when the WD is several billion years old.
In the text
thumbnail Fig. 20

Inward drift rate from Poynting–Robertson drag as a function of the WD cooling age. Even the largest boulders may be dragged enough to trigger gravitational close encounters with the co-orbital asteroid.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.

OSZAR »