Issue 
A&A
Volume 629, September 2019



Article Number  A126  
Number of page(s)  10  
Section  Celestial mechanics and astrometry  
DOI  https://doi.org/10.1051/00046361/201935996  
Published online  16 September 2019 
Driving white dwarf metal pollution through unstable eccentric periodic orbits
^{1}
NaXys, Department of Mathematics, University of Namur, 8 Rempart de la Vierge, 5000 Namur, Belgium
email: kyriaki.antoniadou@unamur.be
^{2}
Department of Physics, University of Warwick, Coventry CV4 7AL, UK
^{3}
Centre for Exoplanets and Habitability, University of Warwick, Coventry CV4 7AL, UK
Received:
31
May
2019
Accepted:
4
August
2019
Context. Planetary debris is observed in the atmospheres of over 1000 white dwarfs, and two white dwarfs are now observed to contain orbiting minor planets. Exoasteroids and planetary core fragments achieve orbits close to the white dwarf through scattering with major planets. However, the architectures that allow for this scattering to take place are timeconsuming to explore with Nbody simulations lasting ∼10^{10} yr; these longrunning simulations restrict the amount of phase space that can be investigated.
Aims. Here we use planar and threedimensional (spatial) elliptic periodic orbits, as well as chaotic indicators through dynamical stability maps, as quick scalefree analytic alternatives to Nbody simulations in order to locate and predict instability in white dwarf planetary systems that consist of one major and one minor planet on very long timescales. We then classify the instability according to ejection versus collisional events.
Methods. We generalized our previous work by allowing eccentricity and inclination of the periodic orbits to increase, thereby adding more realism but also significantly more degrees of freedom to our architectures. We also carried out a suite of computationally expensive 10 Gyr Nbody simulations to provide comparisons with chaotic indicators in a limited region of phase space.
Results. We compute dynamical stability maps that are specific to white dwarf planetary systems and that can be used as tools in future studies to quickly estimate pollution prospects and timescales for oneplanet architectures. We find that these maps also agree well with the outcomes of our Nbody simulations.
Conclusions. As observations of metalpolluted white dwarfs mount exponentially, particularly in the era of Gaia, tools such as periodic orbits can help infer dynamical histories for ensembles of systems.
Key words: celestial mechanics / minor planets / asteroids: general / Kuiper belt: general / white dwarfs / chaos / planets and satellites: dynamical evolution and stability
© ESO 2019
1. Introduction
Observations of white dwarf (WD) exoplanetary systems are abundant and are fast outpacing theoretical efforts to keep up (Veras 2016a). Planetary remnants are now observed in the photospheres of between onequarter and onehalf of Milky Way white dwarfs (Zuckerman et al. 2003, 2010; Koester et al. 2014), over 40 planetary debris disks have been observed to orbit white dwarfs (Zuckerman & Becklin 1987; Gänsicke et al. 2006; Farihi 2016), one exoasteroid is in the process of disintegrating around WD 1145+017 (Vanderburg et al. 2015), and a dense planetary core fragment has been found to orbit SDSS J1228+1040 (Manser et al. 2019). Furthermore, the number of polluted white dwarfs may increase by an order of magnitude over the coming years; this jump will accompany the corresponding increase in known white dwarfs from Gaia in 2018 alone (Gentile Fusillo et al. 2019).
The question then is where the major planets are. Almost certainly, they are hidden from view because white dwarfs are so dim. Nevertheless, the final Gaia data release should reveal at least a dozen planets orbiting white dwarfs through astrometry (Perryman et al. 2014). Moreover, Transiting Exoplanet Survey Satellite (TESS), Large Synoptic Survey Telescope (LSST), and even Laser Interferometer Space Antenna (LISA) have the capability of finding planets that orbit white dwarfs (Cortes & Kipping 2019; Lund et al. 2018; Steffen et al. 2018; Tamanini & Danielski 2018).
The location of these planets is the next question. In our solar system, at least Mars, Jupiter, Saturn, Uranus and Neptune will all survive the Sun’s transition through the giant branch phases into a white dwarf (Schröder & Smith 2008; Veras 2016b). Their semimajor axes will roughly double as a result of mass loss (Omarov 1962; Hadjidemetriou 1963; Veras et al. 2011, 2013a; Dosopoulou & Kalogera 2016a,b). They will accompany a sea of debris produced from supercharged Yarkovsky “O’Keefe” Radzievskii–Paddack (YORP)induced breakup in the asteroid belt (Veras et al. 2014) as well remaining Kuiper belt objects that were orbitally shifted through enhanced Yarkovsky effects (Veras et al. 2015, 2019).
Changes in stability boundaries during the giant branch phase of evolution can trigger gravitational scattering (Debes & Sigurdsson 2002) or delay it for 10^{8}, 10^{9}, or even 10^{10} yr (Veras et al. 2013b, 2016, 2018; Mustill et al. 2014, 2018; Veras & Gänsicke 2015), particularly when secular and mean motion resonance (MMR) sweeping occurs (Bonsor et al. 2011; Debes et al. 2012; Voyatzis et al. 2013; Frewen & Hansen 2014; Smallwood et al. 2018, and in prep.). In at least enough cases for the outcome to be readily observable, the scattering results in small bodies that intersect the white dwarf Roche radius from distances of several au. The subsequent breakup forms a disk (Graham et al. 1990; Jura 2003; Debes et al. 2012; Bear & Soker 2013; Veras et al. 2014, 2017; Brown et al. 2017) that is eventually accreted onto the white dwarf atmosphere and thereby pollutes it.
The authors of several of the references mentioned above performed simulations that ran only on timescales of gigayears or shorter because of natural computational limitations. However, old polluted white dwarfs are common (Hollands et al. 2018), where in this context “old” refers to white dwarfs with cooling ages of at least 5–8 Gyr. Furthermore, fulllifetime simulations that can achieve these timescales (e.g., Veras et al. 2016) have been understandably limited to exploring a narrow band of phase space.
In order to address this difficulty, Antoniadou & Veras (2016) presented an analytical alternative to Nbody simulations by exploring the predictive power of periodic orbits and chaotic indicators through dynamical stability maps (DS maps). They did so in the context of the planar circular restricted threebody problem (CRTBP) with an outer planet (on a circular orbit) and inner asteroid near or inside of the 2:1 MMR. They chose the 2:1 MMR because MMRs have been shown to represent dynamical mechanisms that lead to white dwarf pollution (Bonsor et al. 2011; Voyatzis et al. 2013; Frewen & Hansen 2014; Smallwood et al. 2018, and in prep.), the 2:1 MMR is a strong firstorder resonance, and the 2:1 MMR in particular was a focus of Debes et al. (2012), allowing for comparisons to be made.
Here, we significantly extend Antoniadou & Veras (2016) by instead considering the elliptic restricted threebody problem (ERTBP) as well as, in a more limited capacity, the spatial circular and elliptic restricted threebody problems (3DCRTBP and 3DERTBP).
The ERTBP is more realistic to consider than the CRTBP in white dwarf planetary systems. Even a massive planet on an exactly circular orbit on the main sequence will attain some eccentricity, however small, as a result of giant branch mass loss (Veras et al. 2011). The ERTBP will also allow us to make a direct comparison to previous singleplanet scattering studies, such as that of Debes et al. (2012), which adopted a planetary eccentricity equivalent to the current eccentricity of Jupiter (0.048). The eccentricity of the asteroid through scattering can be excited to the required level (e > 0.99) by a variety of architectures, as demonstrated in the numerical studies listed above, as well as through specific MMRs like the 4:1 (Pichierri et al. 2017). Furthermore, a second major planet in the system can excite the other planet’s eccentricity to any level (Carrera et al. 2019).
In Sect. 2 we provide our model and the basic notions of periodic orbits, linear stability, and the MMR. In Sect. 3 we show the predictive power of the periodic orbits, which can shape the phase space into different domains of stable or chaotic motion. In Sect. 4 we introduce the setup of the Nbody integrations, and in Sect. 5 we present the results of the Nbody simulations. Finally, we summarize our findings in Sect. 6.
2. Model setup, periodic orbits, and MMRs
2.1. Restricted threebody problem
In the framework of the planar restricted threebody problem (2DRTBP), we considered a system consisting of a star, an inner asteroid, and an outer planet of masses m_{S}, m_{A} = 0 and m_{P} = m_{J} = 0.001 m_{⊙}, respectively, and the problem parameter . The asteroid (massless body) does not affect the motions of the star and the planet (primaries) but rather moves under their gravitational attraction. In the 2DRTBP, we defined a suitable planar rotating frame of reference, Oxy, that is centered at the center of mass of the primaries (see, e.g., Hadjidemetriou & Christides 1975).
In the spatial RTBP (3DRTBP), the asteroid was allowed to evolve on an inclined orbit with regard to the plane of the orbits of the primaries. The spatial rotating frame, Gxyz, has an origin that coincides with the center of mass of the massive bodies, its Gzaxis is perpendicular to the Gxyplane, to which the motion of the WD and the planet is confined (see, e.g., Antoniadou & Libert 2019, for more details on dynamics of the 3DRTBP).
Whether in a planar or a spatial rotating frame, the planet was allowed to evolve on a circular e_{P} = 0 (2DCRTBP for a planar or 3DCRTBP for a spatial rotating frame) or an elliptic e_{P} > 0 (2DERTBP for a planar or 3DERTBP for a spatial rotating frame) orbit. More information about the theory of the RTBPs can be found in Szebehely (1967).
We studied periodic orbits under the normalization that sets the sum of the masses and the gravitational constant, G, equal to unity, that is, m_{S} + m_{A} + m_{P} = 1 (or m_{S} = 1 − m_{P} and hence, μ = 0.001) and G = 1.
When viewed in an inertial frame of reference, the orbits of the asteroid and planet correspond to Keplerian ellipses that are characterized by the osculating elements (see, e.g., Murray & Dermott 1999. The notation we use is a (semimajor axis), e (eccentricity), i (inclination), ω (argument of pericenter), M (mean anomaly), and Ω (longitude of ascending node), and ϖ = ω + Ω (longitude of pericenter), Δϖ (apsidal difference), and λ = ϖ + M (mean longitude).
As we study an interior MMR, the massless body (asteroid) evolves in an inner orbit with regard to that of the planet, namely a_{A} < a_{P}. During the computation of the periodic orbits, we always set a_{P} = 1, without loss of generality, so that the period of the planet is T_{0} = 2π.
2.2. Periodic orbits and stability
When we consider a set of positions and velocities, for instance, Q(t), in the phase space of the primaries, then an orbit is periodic if it fulfills the periodic conditions Q(0) = Q(T), where T is the orbital period with t = kT, k ≥ 1 ∈ ℤ. A periodic orbit is symmetric if it is invariant under the fundamental symmetry Σ : (t, x, y) → (−t, x, −y) (Hénon 1973). If the periodic orbit is not invariant, it is asymmetric and Σ will map it to its mirror image. When we take initial conditions on a Poincaré surface of section (Poincaré 1899) of the suitable rotating frame, Oxy or Gxyz, the periodic orbits are the fixed points of the Poincaré map. They can then form a family through different schemes of monoparametric continuation. The origin and continuation of periodic orbits along with the periodicity conditions they should fulfill are described by Antoniadou & Libert (2018a) for the 2DRTBPs and by Antoniadou & Libert (2019) for the 3DRTBPs.
The periodic orbits shape the phase space through their linear stability. We used the conjugate eigenvalues (in reciprocal pairs) of the monodromy matrix of the variational equations of the system in order to classify them as stable or unstable. In general, a periodic orbit is considered stable if and only if all of the eigenvalues lie on the unit circle (see, e.g., Broucke 1969; Marchal 1990, for different types of instabilities). At the orbits where a transition of stability takes place, we have bifurcation points that generate asymmetric periodic orbits. We here plot the stable periodic orbits in blue and the unstable orbits in red. The motion in the domains around the stable orbits is regular, and the longterm stability can be guaranteed thereby. However, in the vicinity of the unstable orbits, chaos emerges and instability events occur, such as collisions or escape.
An intrinsic property of the planar periodic orbits is their vertical stability (Hénon 1973). The vertically stable planar periodic orbits can guarantee the stability of systems with small mutual inclinations around them. The transition of vertical stability to instability along a family is represented by the vertical critical orbit (v.c.o.), which is a bifurcation point that generates spatial periodic orbits. The vertical stability along the planar families is here depicted by solid lines and the vertical instability by dashed lines.
In order to determine the extent of the stable (and/or unstable) regions in the phase space of complex systems, we computed DS maps, which provide a visual representation of the domains. In this respect, we used chaotic indicators, which can help distinguish chaos from order for each of the initial conditions of a DS map. The detrended fast Lyapunov indicator (DFLI; see Froeschlé et al. 1997, for the definition of the FLI), which is the FLI divided by time (Voyatzis 2008), is fast and accurate in tracing chaoticity. When the orbit is regular, the DFLI remains almost constant over time (with values lower than 10). However, when the orbit is chaotic, then the DFLI increases exponentially. We computed the DFLI for a maximum integration time, t_{max} = 250 kyr, which has been proved to be adequate for revealing chaos in the RTBP with the masses we use. In the DS maps, dark colors depict regular orbits, and pale colors highlight traced instability. White shows the failure of the numerical integration at t < t_{max}, which is caused by a small integration step during very close encounters between the bodies. The numerical integration stops when the threshold 10^{30} is reached by the DFLI.
2.3. MMR
In an MMR two proper frequencies of the system become commensurable. An MMR associates the ratio of the mean motions, n_{i}, i standing either for the asteroid, A, or the planet, P, with the rational ratio of two integers, q, p ≠ 0 ∈ ℤ, namely , where q is the order of the resonance. A periodic orbit in the rotating frame shows the exact location of an MMR in phase space.
We here discuss the 2:1 interior MMR, and because the periodic orbits correspond to the stationary solutions of an averaged Hamiltonian, we can define the resonant angles on which they depend as follows:
In the neighborhood of a stable periodic orbit in phase space, the resonant angles and the apsidal difference librate. When all of them librate about 0 or π, the evolution takes place in the vicinity of a symmetric periodic orbit. On the other hand, when they librate about other angles, the periodic orbit is asymmetric. In the neighborhood of unstable periodic orbits, these angles rotate. Given the values of the mean anomalies (0 or π) and the different combinations we can obtain with aligned (Δϖ=0) or antialigned orbits (Δϖ=π), we finally defined four different symmetric configurations: (θ_{1}, θ_{2}) = (0, 0),(0, π),(π, 0), and (π, π). In order to distinguish them on the eccentricities plane, we use the following signconvention. We assign a positive (or negative) value to e_{A} when θ_{1} equals 0 (or π) in the family of each configuration. The same accordingly holds for e_{P} and θ_{2}.
Nonetheless, there also exist other mechanisms that can protect the systems from disruption. For instance, we can observe the secondary resonance inside an MMR, where we consider the average of the frequencies between rotation and libration of θ_{2} and θ_{1}, respectively (Moons & Morbidelli 1993). There also exists the apsidal resonance for the nonresonant orbits, where only the apsidal difference oscillates about 0 or π and the rest of the resonant angles rotate (Murray & Dermott 1999; Morbidelli 2002). For the resonant orbits we can have apsidal resonance only for eccentric orbits that are asymmetric (Michtchenko et al. 2008). Moreover, stability can be attained by apsidal difference oscillation during the transition between two symmetric configurations inside an MMR, where the passage is merely kinematical (Michtchenko et al. 2008). A systematic study of the regions around the periodic orbits in relation to the above mechanisms of protection of the phases can be found in Antoniadou & Libert (2018b) for the 2:1 MMR and Antoniadou & Libert (2018a) for the 3:2, 5:2, 3:1, 4:1, and 5:1 MMRs.
3. Periodic orbits and DS maps as diagnostic tools
3.1. 2DRTBPs
The planar families of periodic orbits in 2:1 MMR when e_{P} > 0 (2DERTBP) for every configuration were computed by Antoniadou & Libert (2018b). We here first employed them to select the symmetric periodic orbits, which then acted as diagnostic tools for the computation of the four DS maps (one for each configuration) on the (e_{A}, e_{P}) plane.
Particularly, as we first aim to explore WD pollution when e_{P} = e_{J} = 0.048, the chosen symmetric periodic orbits correspond to the following initial conditions per symmetric configuration:

(θ_{1}, θ_{2}) = (0, 0) and Δϖ=0: e_{A} = 0.7035, a_{A} = 0.6127, ϖ_{A} = M_{A} = ϖ_{P} = M_{P} = 0°.

(θ_{1}, θ_{2}) = (0, π) and Δϖ=π: e_{A} = 0.7542, a_{A} = 0.6003, ϖ_{A} = M_{P} = 180° and ϖ_{P} = M_{A} = 0°.

(θ_{1}, θ_{2}) = (π, 0) and Δϖ=π: e_{A} = 0.9810, a_{A} = 0.6302, ϖ_{A} = M_{A} = 180° and ϖ_{P} = M_{P} = 0°.

(θ_{1}, θ_{2}) = (π, π) and Δϖ=0: e_{A} = 0.9694, a_{A} = 0.6307, ϖ_{A} = ϖ_{P} = 0° and M_{A} = M_{P} = 180°.
In Fig. 1 we overplot the families of 2:1 resonant symmetric periodic orbits on the DS maps computed on the (e_{A}, e_{P}) plane and highlight the chosen orbits with white crosses. When an island of stability is apparent around the selected periodic orbits (for e_{P} < 0.2), we have a 1:1 secondary resonance where θ_{1} librates around 0 or π according to the symmetric configuration to which the island belongs. More specifically, in the island below the periodic orbit in configuration (0, 0), θ_{1} librates about 0, whereas in the configuration (π, π), it librates about π.
Fig. 1.
Families of symmetric periodic orbits in 2:1 MMR of the 2DERTBP when m_{P} = m_{J} are overplotted on the DS maps of the plane (e_{A}, e_{P}). The white crosses represent the chosen periodic orbits for our study, the orbital elements of which remain fixed for the computation of each map. We use a positive (or negative) value for e_{A} when θ_{1} equals 0 (or π) in the family of each configuration. The same convention applies to e_{P} and θ_{2}. Dark (pale) regions correspond to stable (chaotic) domains in phase space as attributed by the logarithmic values of the DFLI shown by the colored bar. Blue (red) curves stand for the horizontal stability (instability) of the periodic orbits, while the solid (dashed) curves show the vertical stability (instability). The magenta dot corresponds to a transition of the vertical stability. The dashed gray curve depicts the collision line between the planet and the asteroid. 
As the degrees of freedom of the architectures studied here are increased compared to those described in Antoniadou & Veras (2016), the domains in phase space that have to be explored significantly increase as well. Driven by these four periodic orbits, we also computed DS maps on the planes (a_{A}/a_{P}, e_{A}), (ϖ_{A}, M_{A}), and (ϖ_{P}, M_{P}) and show how the domains around the system studied change for each symmetric configuration (see Figs. 2–5).
Fig. 2.
DS maps on the planes (a_{A}/a_{P}, e_{A}) (top), (ϖ_{A}, M_{A}) (middle), and (ϖ_{P}, M_{P}) (bottom) for the symmetric configuration (θ_{1}, θ_{2}) = (0, 0) and Δϖ=0 presented as in Fig. 1. 
Fig. 3.
DS map on the plane (a_{A}/a_{P}, e_{A}) for the symmetric configuration (θ_{1}, θ_{2}) = (0, π) and Δϖ=π. 
Fig. 4.
DS maps for the symmetric configuration (θ_{1}, θ_{2}) = (π, 0) and Δϖ=π presented as in Fig. 2. 
Fig. 5.
DS maps for the symmetric configuration (θ_{1}, θ_{2}) = (π, π) and Δϖ=0 presented as in Fig. 2. 
In particular, we present in Fig. 2 the DS maps for the symmetric configuration (θ_{1}, θ_{2}) = (0, 0) and Δϖ=0. In the small stable domain arising around the periodic orbit, we have 1:1 secondary resonance inside the 2:1 MMR, where θ_{1} librates about 0, and longterm stability is expected therein.
In Fig. 3 we present one DS map for the symmetric configuration (θ_{1}, θ_{2}) = (0, π) and Δϖ=π, that is, only on the plane (a_{A}/a_{P}, e_{A}), because the DS maps on the planes (ϖ_{A}, M_{A}) and (ϖ_{P}, M_{P}) resulted solely in chaotic orbits. The periodic orbit is surrounded by chaotic domains, and instability events are expected to take place in this dynamical neighborhood.
In Fig. 4 we reveal the symmetric configuration (θ_{1}, θ_{2}) = (π, 0) and Δϖ=π. Below the periodic orbit (for e_{A} ≈ 0.8) there exists an island of stability, where the regularity of the orbits is maintained as a result of the 1:1 secondary resonance inside the 2:1 MMR, where the resonant angle θ_{1} librates about π.
In Fig. 5, where the phase space of the configuration (θ_{1}, θ_{2}) = (π, π) and Δϖ=0 is shown, we have an island of stability around the periodic orbit (for e_{A} > 0.8), where the 1:1 secondary resonance takes place with θ_{1} librating about π. This is expected to preserve stability for longtime spans.
Notably, on the plane (a_{A}/a_{P}, e_{A}) of the Figs. 2–5, an instability region is observed for very low values of e_{A} at 2:1 MMR. This region is justified by the Poincaré surface of sections, which have long proved that the chaotic regions therein are bounded and therefore, despite the irregular motion, escape cannot occur. These regions are also sensitive to the eccentricity of Jupiter, which in the DS maps here is equal to 0.048 (ERTBP), a value that is still very close to zero (see the DS maps in the CRTBP studied by Antoniadou & Veras 2016). We refer to Morbidelli & Moons (1993), Michtchenko & FerrazMello (1995), Nesvorný & FerrazMello (1997) and Hadjidemetriou & Voyatzis (2000), for instance, for more details on the dynamics of this particular region.
An additional stable region observed on the plane (a_{A}/a_{P}, e_{A}) for a_{A}/a_{P} ≈ 0.57 and e_{A} > 0.9 in Figs. 3 and 4 is also worth mentioning. This domain is associated with the 7:3 MMR, given the libration of the resonant angles and the apsidal difference. Although this is intriguing, it is beyond the scope of this study.
To our knowledge, asymmetric configurations (Δϖ≠0 or π) linked with WD pollution do not exist in the literature. In order to address the issue of geometric asymmetry (Δϖ=ϖ_{P} − ϖ_{A}) and dynamical asymmetry (ΔM = M_{P} − M_{A}) of the orbits in our Nbody simulations and provide a first insight of how these angles can also cause instability events, we use as initial conditions the following asymmetric periodic orbits, which belong to the Family 1 of Antoniadou & Libert (2018b):

1: e_{A} = 0.0860, e_{P} = 0.5580, a_{A} = 0.6298, ϖ_{A} = M_{A} = 0°, ϖ_{P} = 91.50°, M_{P} = 218.65°.

2: e_{A} = 0.2000, e_{P} = 0.6424, a_{A} = 0.6294, ϖ_{A} = M_{A} = 0°, ϖ_{P} = 322.13°, M_{P} = 294.88°.

3: e_{A} = 0.4009, e_{P} = 0.7481, a_{A} = 0.6293, ϖ_{A} = M_{A} = 0°, ϖ_{P} = 324.51°, M_{P} = 302.12°.
A detailed and exhaustive exploration of the way the asymmetry of orbits affects WD pollution is not the purpose of our present study. By relying on the linear stability of the above asymmetric periodic orbits (all of them are unstable), we directly perform the Nbody simulations in Sect. 5, and provide a neat example of instability events that are linked with them.
3.2. 3DRTBPs
Although in the literature many Nbody investigations of postmainsequence exoplanetary systems have so far assumed coplanarity or nearcoplanarity, it is likely given from our solar system knowledge that none of the asteroids is perfectly coplanar with the planet. In order to address the physical issue of coplanarity, we therefore studied the neighborhood of a spatial periodic orbit in the 3DCRTBP.
More specifically, in Fig. 6 we selected the family of unstable symmetric periodic orbits in the 3DCRTBP (e_{P} = 0) shown by Antoniadou & Libert (2019). For reasons of completeness, it is presented for both prograde (i_{A} < 90°) and retrograde (i_{A} > 90°) orbits. The chosen periodic orbit (white cross) for the computation of the DS map on the (e_{A}, i_{A}) plane has orbital elements:
Fig. 6.
Family of unstable periodic orbits in 2:1 MMR of the 3DCRTBP when m_{P} = m_{J} is overplotted on the DS map of the plane (e_{A}, i_{A}). Presentation as in Fig. 1. 

e_{A} = 0.7922, a_{A} = 0.6295, i_{A} = 15.03°, i_{P} = 0.0°, ϖ_{A} = M_{A} = ϖ_{P} = M_{P} = 0°.
With regard to the 3DERBTP, we only focused on the effect of the vertical stability of the four planar symmetric periodic orbits (white crosses) selected in Fig. 1. We did not study the spatial family of symmetric periodic orbits emanating from the magenta dot (v.c.o. linking the 2DERTBP with the 3DERTBP, see Antoniadou & Libert 2019 for its depiction) in Fig. 1 because it consists of stable spatial periodic orbits.
4. Nbody simulations
4.1. Integrator setup
We used the BulirschStoer integrator from the modified version of MERCURY (Chambers 1999) and imposed an accuracy of 10^{−13}. This provides conservation of energy and angular momentum in the range of 10^{−8}–10^{−12}. We also configured the ejection radius to be at 3 × 10^{5} au.
We ran simulations with sets of 10–50 massless asteroids for different fixed combinations of orbital elements and for up to 14 Gyr, as in Antoniadou & Veras (2016). Many simulations did not run for that long because of the increased parameter space and limited computational resources. Even just 10 Gyr is still older than every known polluted WD (Hollands et al. 2018). The output frequency was 1 Myr.
4.2. Choice of parameters
In addition to the periodic orbits, which expedite the guesswork by providing the initial conditions in order to search the phase space, Nbody integrations were our second means to help us determine the orbital evolution over long time spans. The link between our two tools has to be carefully determined, so that the initial parameters given by the periodic orbits are correctly transformed so that the real system, where the mass of the WD is m_{WD} = 0.6 m_{⊙}, can be accurately studied and represented. In order to convert the data given by the periodic orbits into real units and thus incorporate the periodic orbits into the simulations, we therefore introduce a scaling factor for the sake of the invariance of the equations of motion for the periodic orbits.
We performed the same scaling as was introduced in Antoniadou & Veras (2016) with regard to the semimajor axes used in the simulations,
where is the scaling factor and N represents the orbital elements used in the Nbody simulations. The rest of the orbital elements and the time remained the same.
Furthermore, in all cases we set a_{P} = 10 au (or scaled as au). As the internal properties of the bodies are beyond the scope of our study, we adopted a fictitious WD radius equal to R_{WD} = 10^{6} km ≈1.4378 R_{⊙} and a planetary radius corresponding to R_{P} = 78 000 km or 1.09 R_{J}. The value of R_{WD} conforms well to a conservatively large estimate of the WD Roche radius, see Veras et al. (2017), and during our simulations, we assumed that a collision occurs with the WD when an asteroid intersects this radius.
4.3. Instability timescales and output visualization
We used both colors and shapes in order to visualize different evolution outputs of our Nbody integrations. Namely:

Green circles correspond to stablility

Red triangles represent the ejection of the asteroid

Brown squares indicate the collision between the asteroid and the planet

Purple diamonds signify the collision between the asteroid and the WD
For the stable simulations, we used a colorshading range of [10^{7}, 1.4 × 10^{10}] yr, and for the other simulations, we used a colorshading range of [10^{3}, 1.4 × 10^{10}] yr. The darker the shading, the longer this simulation ran.
We note that in the Nbody integrations stability only reflects the absence of ejection of the asteroid, or collision with the planet or the WD. Even though an evolution of the orbital elements might be irregular (within the neighborhood of unstable periodic orbits), their initial values might not change significantly over time, as is the case of weakly chaotic orbits, which are considered as stable from a physical point of view. This nearinvariance suggests that the output of the Nbody integration would then be classified as stable in the adopted visual representation of the Nbody simulation outcomes.
5. Results of Nbody simulations
In this section we illustrate the agreement between the outcome of the Nbody simulations and the DS maps based on the periodic orbits by presenting a selection of results. Based on the DS maps shown in Sect. 3, we chose the initial conditions for which the DFLI traced chaoticity (DFLI > 15). Then, we plot these initial conditions in the background using a uniform color for all of them and overplot the outcomes of the Nbody simulations in order to compare and contrast the results derived from our tools.
5.1. 2DRTBPs
We here compare the chaotic domains around unstable periodic orbits shown in Sect. 3.1 and the Nbody simulations for the 2DRTBPs for the symmetric and asymmetric cases. More specifically, in Fig. 7 we present the results on the plane (a_{A}/a_{P}, e_{A}) for the configurations (θ_{1}, θ_{2}) = (0, 0) (top panel), (θ_{1}, θ_{2}) = (π, 0) (middle panel), and (θ_{1}, θ_{2}) = (π, π) (bottom panel). Debes et al. (2012) showed that the majority of tidally disrupted asteroids arise from the boundaries of the libration width of 2:1 MMR (Fig. 2 therein). Pichierri et al. (2017) concluded that of the 2:1, 3:1, and 4:1 MMRs, only the 4:1 MMR could excite the eccentricity of the asteroid enough for it to be disrupted by falling onto the central star. By searching for an increased concentration of purple diamonds in each of the panels of Fig. 7, we observe that the configuration (θ_{1}, θ_{2}) = (π, 0) and Δϖ=π (middle panel) can generate WD pollution for the orbits close to the 2:1 MMR.
Fig. 7.
Results of the Nbody simulations on the plane (a_{A}/a_{P}, e_{A}) for the symmetric configuration (θ_{1}, θ_{2}) = (0, 0) and Δϖ=0 (top), (θ_{1}, θ_{2}) = (π, 0) and Δϖ=π (middle) and (θ_{1}, θ_{2}) = (π, π) and Δϖ=0 (bottom). Green circles stand for stable outcome, red triangles for an ejected asteroid, brown squares for an asteroid that hit the planet and purple diamonds for an asteroid that collided with the WD. In the background, the initial conditions with DFLI > 15 are included and uniformly depicted by pale gray. 
In Fig. 8 we present the results for the configuration where (θ_{1}, θ_{2}) = (π, π) and Δϖ=0 on the planes (e_{A}, e_{P}) (top panel) and (ϖ_{A}, M_{A}) (bottom panel). The top panel shows that the tidal disruption of the asteroid can take place for any value of e_{A} when e_{P} > 0.1. Additionally, we observe the remarkable role of the angles in the evolution of the asteroid, as the bottom panel is densely populated by purple diamonds when the stability region is avoided.
Fig. 8.
Results of the Nbody simulations for the symmetric configuration (θ_{1}, θ_{2}) = (π, π) and Δϖ=0 on the planes (e_{A}, e_{P}) (top) and (ϖ_{A}, M_{A}) (bottom). 
Because of the strong dependence of stability on orbital angles, we wished to explore their effect by examining the neighborhood of unstable asymmetric periodic orbits, whose angles vary within the family they belong to and are not constant, like for the symmetric periodic orbits. In Fig. 9 we therefore merged the results of the Nbody simulations for the three asymmetric configurations. In none of them was a stable evolution observed, while in the majority of the simulations, the asteroid either hit the WD or was ejected.
Fig. 9.
Results of the Nbody simulations for the three asymmetric configurations presented on the (e_{A}, e_{P})plane. 
5.2. 3DRTBPs
The aim of this section is to provide clues when architectures with inclined asteroids are considered. First we consider the 3DCRTBP. Based on the DS map of Fig. 6, the Nbody simulations are not encouraged because the region of instability around the unstable family is very small for prograde orbits. This means that the 3DCRTBP, when the asteroid inclination is small, cannot trigger WD pollution, just like the 2DCRTBP that was explored in Antoniadou & Veras (2016).
With regard to the 3DERTBP, we described above that the neighborhood of spatial periodic orbits was not studied because in the elliptic problem, all of them were found to be stable in the 2:1 MMR. From a dynamical point of view, architectures in this problem can therefore only be studied when the vertical stability of the periodic orbits in the 2DERTBP is taken into account. In Fig. 10 we used as initial conditions the four planar symmetric periodic orbits (defined in Sect. 3.1 and delineated by white crosses) that belong to each of the four symmetric configurations (shown in Fig. 1) and added inclination (with the longitudes of ascending node Ω_{A} = Ω_{P} = 0), so that we study the outcome of the simulations in the 3DERTBP. We recall that all of these periodic orbits are vertically unstable (dashed lines in Fig. 1). Figure 10 shows that when the periodic orbit was horizontally unstable (red lines in Fig. 1), instability events were numerous, even for very low inclination values. We may say that WD pollution can be obtained through the configurations (θ_{1}, θ_{2}) = (π, 0) and (θ_{1}, θ_{2}) = (π, π), while in the configuration (θ_{1}, θ_{2}) = (0, 0), an inclination i_{A} > 48° should be achieved first.
Fig. 10.
Results of the Nbody simulations when we considered architectures in the 3DERTBP for each symmetric configuration shown in Fig. 1. The white crosses represent the chosen symmetric periodic orbits with e_{P} = 0.048 from the 2DERTBP, which act as initial conditions for the simulations in the 3DERTBP, as inclination is added to the asteroids thus. 
6. Summary
The old age of white dwarf planetary systems introduces modeling challenges, particularly for longterm numerical simulations. Analytic and semianalytic models can provide timesaving alternatives while simultaneously revealing fundamental insights about structures such as the threebody problem.
Unstable periodic orbits with high eccentricity are here the analytic channel through which we studied the prospect of an eccentric planet that perturbs an asteroid into a white dwarf, thereby polluting it with metals. These periodic orbits enabled us to compute dynamical (in)stability maps for white dwarf planetary systems consisting of one giant planet with the same eccentricity and mass as Jupiter^{1}, plus an internal asteroid (Figs. 1–6). Supplemented with a suite of computationally expensive Nbody simulations, the maps provide a useful diagnostic for identifying initial conditions that can generate white dwarf metal pollution in 2:1 MMR.
In order to link polluted white dwarfs with architectures that arise from planet formation, the initial conditions corresponding to the purple diamonds in Figs. 7–10 could represent finalstate targets for simulations of planetary systems that traverse the mainsequence and giant branch phases. In a comprehensive novel work, Harrison et al. (2018) made the important chemical association between protoplanetarydisk formation locations and condensation temperatures with particular white dwarf pollutants. A next step could be to match the chemical connection with a planetasteroid dynamical association, and the dynamical maps presented here may facilitate this task.
A narrower but still important application of our dynamical maps would be to explore how the accretion rate onto white dwarf systems would change if other massive bodies were present in the system in addition to the one planet considered here. For example, major planets that lie close to the Roche radius of the white dwarf would scatter incoming asteroids, either enhancing or lowering the accretion rate compared to the oneplanet case. The two minor planets already found at or within the Roche radius (Vanderburg et al. 2015; Manser et al. 2019) have prompted additional searches for major planets, and demonstrates that planetary orbits that lie entirely within one solar radius can be achieved.
Finally, the results provided in this work have applications in the broader context of celestial mechanics, and in particular for any architecture that can be efficiently modeled by the 2DERTBP, 3DCRTBP, and 3DERTBP. Examples range from Earthbased architectures (planetspacecraftsatellite, where the spacecraft is massless), to solar systembased architectures (SunEarthJupiter, where Earth is massless) to extrasolar systembased architectures (starstarcircumprimary planet, where the planet is considered massless).
Acknowledgments
We thank the referee for their helpful comments that have improved the manuscript. Computational resources have been provided by the Consortium des Équipements de Calcul Intensif, funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.FNRS) under Grant No. 2.5020.11 and by the Walloon Region. DV gratefully acknowledges the support of the STFC via an Ernest Rutherford Fellowship (grant ST/P003850/1).
References
 Antoniadou, K. I., & Libert, A.S. 2018a, CeMDA, 130, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Antoniadou, K. I., & Libert, A.S. 2018b, A&A, 615, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Antoniadou, K. I., & Libert, A.S. 2019, MNRAS, 483, 2923 [NASA ADS] [CrossRef] [Google Scholar]
 Antoniadou, K. I., & Veras, D. 2016, MNRAS, 463, 4108 [NASA ADS] [CrossRef] [Google Scholar]
 Bear, E., & Soker, N. 2013, New Astron., 19, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Bonsor, A., Mustill, A. J., & Wyatt, M. C. 2011, MNRAS, 414, 930 [NASA ADS] [CrossRef] [Google Scholar]
 Broucke, R. 1969, AIAA J., 7, 1003 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, J. C., Veras, D., & Gänsicke, B. T. 2017, MNRAS, 468, 1575 [NASA ADS] [CrossRef] [Google Scholar]
 Carrera, D., Raymond, S. R., & Davies, M. B. 2019, MNRAS Lett., submitted [arXiv:1903.02564] [Google Scholar]
 Chambers, J. E. 1999, MNRAS, 304, 793 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Cortes, J., & Kipping, D. M. 2019, MNRAS, 488, 1695 [NASA ADS] [CrossRef] [Google Scholar]
 Debes, J. H., & Sigurdsson, S. 2002, ApJ, 572, 556 [NASA ADS] [CrossRef] [Google Scholar]
 Debes, J. H., Walsh, K. J., & Stark, C. 2012, ApJ, 747, 148 [NASA ADS] [CrossRef] [Google Scholar]
 Dosopoulou, F., & Kalogera, V. 2016a, ApJ, 825, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Dosopoulou, F., & Kalogera, V. 2016b, ApJ, 825, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Farihi, J. 2016, New Astron. Rev., 71, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Frewen, S. F. N., & Hansen, B. M. S. 2014, MNRAS, 439, 2442 [NASA ADS] [CrossRef] [Google Scholar]
 Froeschlé, C., Lega, E., & Gonczi, R. 1997, CeMDA, 67, 41 [CrossRef] [Google Scholar]
 Gänsicke, B. T., Marsh, T. R., Southworth, J., & RebassaMansergas, A. 2006, Science, 314, 1908 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Gentile Fusillo, N. P., Tremblay, P.E., Gänsicke, B. T., et al. 2019, MNRAS, 482, 4570 [NASA ADS] [CrossRef] [Google Scholar]
 Graham, J. R., Matthews, K., Neugebauer, G., & Soifer, B. T. 1990, ApJ, 357, 216 [NASA ADS] [CrossRef] [Google Scholar]
 Hadjidemetriou, J. D. 1963, Icarus, 2, 440 [NASA ADS] [CrossRef] [Google Scholar]
 Hadjidemetriou, J. D., & Christides, T. 1975, Celestial Mech., 12, 175 [NASA ADS] [CrossRef] [Google Scholar]
 Hadjidemetriou, J. D., & Voyatzis, G. 2000, CeMDA, 78, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Harrison, J. H. D., Bonsor, A., & Madhusudhan, N. 2018, MNRAS, 479, 3814 [NASA ADS] [CrossRef] [Google Scholar]
 Hénon, M. 1973, A&A, 28, 415 [NASA ADS] [Google Scholar]
 Hollands, M. A., Gänsicke, B. T., & Koester, D. 2018, MNRAS, 477, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Jura, M. 2003, ApJ, 584, L91 [Google Scholar]
 Koester, D., Gänsicke, B. T., & Farihi, J. 2014, A&A, 566, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lund, M. B., Pepper, J. A., Shporer, A., & Stassun, K. G. 2018, AAS J., submitted [arXiv:1809.10900] [Google Scholar]
 Manser, C. J., Gänsicke, B. T., Eggl, S., et al. 2019, Science, 364, 66 [Google Scholar]
 Marchal, C. 1990, The Threebody Problem (Amsterdam: Elsevier) [Google Scholar]
 Michtchenko, T. A., & FerrazMello, S. 1995, A&A, 303, 945 [NASA ADS] [Google Scholar]
 Michtchenko, T. A., Beaugé, C., & FerrazMello, S. 2008, MNRAS, 387, 747 [NASA ADS] [CrossRef] [Google Scholar]
 Moons, M., & Morbidelli, A. 1993, CeMDA, 56, 273 [NASA ADS] [CrossRef] [Google Scholar]
 Morbidelli, A. 2002, Modern Celestial Mechanics: Aspects of Solar System Dynamics (London: Taylor & Francis) [Google Scholar]
 Morbidelli, A., & Moons, M. 1993, Icarus, 102, 316 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge University Press) [Google Scholar]
 Mustill, A. J., Veras, D., & Villaver, E. 2014, MNRAS, 437, 1404 [NASA ADS] [CrossRef] [Google Scholar]
 Mustill, A. J., Villaver, E., Veras, D., Gänsicke, B. T., & Bonsor, A. 2018, MNRAS, 476, 3939 [NASA ADS] [CrossRef] [Google Scholar]
 Nesvorný, D., & FerrazMello, S. 1997, A&A, 320, 672 [NASA ADS] [Google Scholar]
 Omarov, T. B. 1962, Izv. Astrofiz. Inst. Acad. Nauk. KazSSR, 14, 66 [Google Scholar]
 Perryman, M., Hartman, J., Bakos, G. Á., & Lindegren, L. 2014, ApJ, 797, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Pichierri, G., Morbidelli, A., & Lai, D. 2017, A&A, 605, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Poincaré, H. 1899, Les méthodes nouvelles de la méchanique céleste (Paris: GauthierVillars), III [Google Scholar]
 Schröder, K.P., & Smith, R. 2008, MNRAS, 386, 155 [NASA ADS] [CrossRef] [Google Scholar]
 Smallwood, J. L., Martin, R. G., Livio, M., & Lubow, S. H. 2018, MNRAS, 480, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Steffen, J. H., Wu, D.H., & Larson, S. L. 2018, AAS J., submitted [arXiv:1812.03438] [Google Scholar]
 Szebehely, V. 1967, Theory of Orbits. The Restricted Problem of Three Bodies (New York: Academic Press) [Google Scholar]
 Tamanini, N., & Danielski, C. 2018, ArXiv eprints [arXiv:1812.04330] [Google Scholar]
 Vanderburg, A., Johnson, J. A., Rappaport, S., et al. 2015, Nature, 526, 546 [NASA ADS] [CrossRef] [Google Scholar]
 Veras, D. 2016a, R. Soc. Open Sci., 3, 150571 [NASA ADS] [CrossRef] [Google Scholar]
 Veras, D. 2016b, MNRAS, 463, 2958 [NASA ADS] [CrossRef] [Google Scholar]
 Veras, D., & Gänsicke, B. T. 2015, MNRAS, 447, 1049 [NASA ADS] [CrossRef] [Google Scholar]
 Veras, D., Wyatt, M. C., Mustill, A. J., Bonsor, A., & Eldridge, J. J. 2011, MNRAS, 417, 2104 [Google Scholar]
 Veras, D., Hadjidemetriou, J. D., & Tout, C. A. 2013a, MNRAS, 435, 2416 [NASA ADS] [CrossRef] [Google Scholar]
 Veras, D., Mustill, A. J., Bonsor, A., & Wyatt, M. C. 2013b, MNRAS, 431, 1686 [NASA ADS] [CrossRef] [Google Scholar]
 Veras, D., Leinhardt, Z. M., Bonsor, A., & Gänsicke, B. T. 2014, MNRAS, 445, 2244 [NASA ADS] [CrossRef] [Google Scholar]
 Veras, D., Eggl, S., & Gänsicke, B. T. 2015, MNRAS, 451, 2814 [NASA ADS] [CrossRef] [Google Scholar]
 Veras, D., Mustill, A. J., Gänsicke, B. T., et al. 2016, MNRAS, 458, 3942 [NASA ADS] [CrossRef] [Google Scholar]
 Veras, D., Carter, P. J., Leinhardt, Z. M., & Gänsicke, B. T. 2017, MNRAS, 465, 1008 [NASA ADS] [CrossRef] [Google Scholar]
 Veras, D., Georgakarakos, N., Gänsicke, B. T., & DobbsDixon, I. 2018, MNRAS, 481, 2180 [NASA ADS] [CrossRef] [Google Scholar]
 Veras, D., Higuchi, A., & Ida, S. 2019, MNRAS, 485, 708 [Google Scholar]
 Voyatzis, G. 2008, ApJ, 675, 802 [NASA ADS] [CrossRef] [Google Scholar]
 Voyatzis, G., Hadjidemetriou, J. D., Veras, D., & Varvoglis, H. 2013, MNRAS, 430, 3383 [NASA ADS] [CrossRef] [Google Scholar]
 Zuckerman, B., & Becklin, E. E. 1987, Nature, 330, 138 [NASA ADS] [CrossRef] [Google Scholar]
 Zuckerman, B., Koester, D., Reid, I. N., & Hünsch, M. 2003, ApJ, 596, 477 [NASA ADS] [CrossRef] [Google Scholar]
 Zuckerman, B., Melis, C., Klein, B., Koester, D., & Jura, M. 2010, ApJ, 722, 725 [Google Scholar]
All Figures
Fig. 1.
Families of symmetric periodic orbits in 2:1 MMR of the 2DERTBP when m_{P} = m_{J} are overplotted on the DS maps of the plane (e_{A}, e_{P}). The white crosses represent the chosen periodic orbits for our study, the orbital elements of which remain fixed for the computation of each map. We use a positive (or negative) value for e_{A} when θ_{1} equals 0 (or π) in the family of each configuration. The same convention applies to e_{P} and θ_{2}. Dark (pale) regions correspond to stable (chaotic) domains in phase space as attributed by the logarithmic values of the DFLI shown by the colored bar. Blue (red) curves stand for the horizontal stability (instability) of the periodic orbits, while the solid (dashed) curves show the vertical stability (instability). The magenta dot corresponds to a transition of the vertical stability. The dashed gray curve depicts the collision line between the planet and the asteroid. 

In the text 
Fig. 2.
DS maps on the planes (a_{A}/a_{P}, e_{A}) (top), (ϖ_{A}, M_{A}) (middle), and (ϖ_{P}, M_{P}) (bottom) for the symmetric configuration (θ_{1}, θ_{2}) = (0, 0) and Δϖ=0 presented as in Fig. 1. 

In the text 
Fig. 3.
DS map on the plane (a_{A}/a_{P}, e_{A}) for the symmetric configuration (θ_{1}, θ_{2}) = (0, π) and Δϖ=π. 

In the text 
Fig. 4.
DS maps for the symmetric configuration (θ_{1}, θ_{2}) = (π, 0) and Δϖ=π presented as in Fig. 2. 

In the text 
Fig. 5.
DS maps for the symmetric configuration (θ_{1}, θ_{2}) = (π, π) and Δϖ=0 presented as in Fig. 2. 

In the text 
Fig. 6.
Family of unstable periodic orbits in 2:1 MMR of the 3DCRTBP when m_{P} = m_{J} is overplotted on the DS map of the plane (e_{A}, i_{A}). Presentation as in Fig. 1. 

In the text 
Fig. 7.
Results of the Nbody simulations on the plane (a_{A}/a_{P}, e_{A}) for the symmetric configuration (θ_{1}, θ_{2}) = (0, 0) and Δϖ=0 (top), (θ_{1}, θ_{2}) = (π, 0) and Δϖ=π (middle) and (θ_{1}, θ_{2}) = (π, π) and Δϖ=0 (bottom). Green circles stand for stable outcome, red triangles for an ejected asteroid, brown squares for an asteroid that hit the planet and purple diamonds for an asteroid that collided with the WD. In the background, the initial conditions with DFLI > 15 are included and uniformly depicted by pale gray. 

In the text 
Fig. 8.
Results of the Nbody simulations for the symmetric configuration (θ_{1}, θ_{2}) = (π, π) and Δϖ=0 on the planes (e_{A}, e_{P}) (top) and (ϖ_{A}, M_{A}) (bottom). 

In the text 
Fig. 9.
Results of the Nbody simulations for the three asymmetric configurations presented on the (e_{A}, e_{P})plane. 

In the text 
Fig. 10.
Results of the Nbody simulations when we considered architectures in the 3DERTBP for each symmetric configuration shown in Fig. 1. The white crosses represent the chosen symmetric periodic orbits with e_{P} = 0.048 from the 2DERTBP, which act as initial conditions for the simulations in the 3DERTBP, as inclination is added to the asteroids thus. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.