© ESO 2021

1 Introduction

Hirayama (1918) was the first to identify concentrations of asteroids with very similar orbital elements that are known today as asteroid families. They are intensively investigated because they provide unique possibilities for studying different processes related to high-energy collisions. According to Nesvorný et al. (2015), more than 120 asteroid families are currently known. A summary of the current state of affairs in the field of asteroid family identification is given in Michel et al. (2015).

Most of the large asteroid families are at least 1 Gyr old (Milani et al. 2017; Carruba et al. 2016). The dynamic history of old asteroid families is affected by a number of limiting factors, such as incomplete data on the Yarkovsky effect, chaos effect, and subsequent collisions (Marzari et al. 1999; Nesvorný et al. 2006), which make reconstructing their dynamical past difficult. The Yarkovsky–O’Keefe–Radzievskii–Paddack (YORP) effect (Carruba et al. 2016) and the chaotic scattering by relatively large trans-Neptunian objects affecting the older families (Galiazzo et al. 2016) may contribute to complicate the already difficult task of estimating their dynamical age even more. The younger an asteroid family, the more straightforward the analysis because the members of such a family only represent the outcome of the process of fragmentation of their parent body prior to the subsequent dynamical and physical changes induced by the orbital and physical evolutionary processes acting over time.

The study of young asteroid families has become a subject of significant interest within the field (see e.g. Pravec et al. 2018). Recent studies have shown that some asteroid families can also be the outcome of the spin-up-induced fission of a critically rotating parent body (fission clusters), see Jacobson & Scheeres (2011). Moreover, cases of subsequent breakups can take place in older families (Carruba et al. 2020). In this paper, we target the effect of the mean motion resonances on changing the semi-major axis of family members and pairs. We focus on the Datura asteroid family because it is well known (Nesvorný & Vokrouhlický 2006) and is very close to the 9:16 mean motion resonance with Mars. Another very interesting young asteroid family is the newly discovered Adelaide family (Vokrouhlický et al. 2021). This family is located near the Datura family region, and future study could be very informative. As a secondary case study, we used the young pair Balam and 2009 BR₆₀.

In all instances we refer to standard notations for orbital elements. a is the semi-major axis in AU, e is the eccentricity, i is the inclination, Ω is the longitude of the ascending node, ω is the argument of perihelion, $\overline{\omega }$ is the longitude of perihelion, which is $\overline{\omega }=\Omega +\omega$ (see e.g. Murray & Dermott 2000). The angular elements are given in degrees. In general, analytic proper orbital elements are based on the solution of averaged equations of motion, and synthetic proper orbital elements are based on the numerical integration results.

We integrated the equations of motion of the systems numerically 800 kyr into the past to study the dynamical evolution of asteroids using the N-body integrator Mercury (Chambers 1999) and the Everhart integration method (Everhart 1985). We used a time-step of one day and an accuracy parameter of 10⁻¹⁸ in our integrations. Considering that the age of the Datura family was estimated in the range of 400–500 kyr (Vokrouhlický et al. 2008), we restricted ourselves to a time interval of integration of 800 kyr.

We used three integration cases in this research. In the first case, we only took into account disturbances from major planets (onlyplanet). In the second case, we added disturbances from Ceres, Vesta, Juno, and Pallas (CVJP). In the third case, we added disturbances from all previous celestial bodies plus Hygiea, Eunomia, Interamnia, and Davida (CVJPHIDE). We retrieved the data used as initial orbital elements from the JPL Horizons website¹. For the masses of large asteroids, which are not published in JPL Horizons, we used the values of masses calculated by the technique ofastrometric mass determination (Baer & Chesley 2008). Table 1 shows the elements of the asteroids discussed in this paper.

We used the preliminary results from Rosaev & Plávalová (2018) as a basis and calculated the values of new numerical integrations.We report on the importance of studying weak effects on the evolution of asteroid orbits, especially the 9:16 resonance ofMars. We organise the sections as follows: The main information about the Datura family is given in Sect. 2. In Sect. 3 we discuss the Lyapunov characteristic exponent (Sect. 3.1) and then determine the resonance position with the Euclid algorithm (Sect. 3.2). We explain our considerations on resonance boundaries in Sect. 3.3 and compare our calculation for the Mars orbit with those of Laskar et al. (2002) results in Sect. 3.4. We propose a simple theoretical model in Sect. 4 and discuss the effect of the large asteroid perturbations on the Datura family members in Sect. 5. We evaluate the degree of the relevance of the Yarkovsky effect on changing the orbits of asteroids (Sect. 6). As another example in which the 9:16 resonance with Mars could be visible, we introduce the pair 3749 Balam and 2009 BR₆₀ in Sect. 7, and finally, we offer our conclusions in Sect. 8.

2 Datura family main facts

The group of minor planets associated with the largest body of (1270) Datura is a young asteroid family. Since its discovery by Nesvorný et al. (2006), both its dynamical evolution and physical properties have been rigorously studied. Vokrouhlický et al. (2008) first postulated that the formation of the Datura asteroid family took place nearly 450 kyr ago. They modelled the dust distributions after the Datura breakup. The Datura family consists of one large 10 km-sized asteroid (parent body) and a few small minor planets, possible fragments resulting from a catastrophic breakup. Until recently, only seven members of the Datura family were known (Nesvorný et al. 2006, 2015; Vokrouhlický et al. 2009), all of them belonging to the so-called S-type of minor planets (Tholen 1984). According to our latest analyses (Rosaev & Plávalová 2017), however, three new members have to be included, and Vokrouhlický et al. (2017) have reported on the existence of another seven. The total sum of known members of the Datura family might therefore be as high as 17 at the time of writing (see Table 1). The obvious question is whether all of these 17 members are the product of a single breakup event or if some of them were brought into existence by subsequent fragmentations. Fatka et al. (2020) introduced the idea of cascade disruptions of some asteroid families.

3 Problem setting

Nesvorný & Vokrouhlický (2006) pointed out that the Datura family is very close to the 9:16M mean motion resonance with Mars. Vokrouhlický et al. (2017) provided a detailed analysis of this resonance in the terms of relative velocity. In particular, they noted that three Datura members, (89309) 2001 VN₃₆, (433382) 2013 ST₇₁, and (485010) 2009 VS₁₁₆, are currently located in this resonance. We consider some other aspects of the interaction of some Datura members with this resonance and offer qualitative explanations of some observed features.

Vokrouhlický et al. (2017) noted that one member, (452713) 2005 YP₁₃₆, even appears to jump over this resonance, accumulating a large perturbation in its proper orbital elements. In our previous work, we outlined this jump in the orbit of 2014 OE₂₀₆ (Rosaev & Plávalová 2017). According to the results of backward integrations when Ceres and Vesta perturbations were taken into account, a further three members (2009 VS₁₁₆, 2001 VN₃₆, and 2011 KT₁₀) have exhibited a similar behaviour (Rosaev & Plávalová 2018). The qualitative explanation of these jumps is the main goal of our present paper.

To study the dependence of the results on initial data variations, we repeated the integration with a limited number of clones within the 1σ interval for three selected members of the Datura family, namely (89309) 2001 VN₃₆, (485010) 2009 VS₁₁₆, and (485010) 2009 VS₁₁₆ (Table 2). We excluded 2014 OE₂₀₆ from our consideration because its orbit is poorly determined. We used six clones for each orbital element: six clones for eccentricity, six clones for orbital inclination, etc. We can accurately describe the most recent (closest in time) encounter and more or less describe the secondary encounter, but beyond this, we can only give a range of possible orbits. In the same way, we can provide a detailed description of the most recent interaction of the orbit of an asteroid with resonance, more or less plot the interval between first and second approach with resonance, but after this, we can only give a range of possible orbits.

Moreover, the question of the interaction between the orbit of an asteroid with resonance cannot be studied by classical statistical methods because the results of the second approach might be higher or lower than the nominal value of resonance with the comparable possibility (Fig. 1). The distribution is therefore expected to be bimodal. Surprisingly, the uncertainties on inclination, eccentricity, node, and perihelion longitudes have a greater effect on the asteroid orbits we studied than the uncertainties on mean anomaly and semi-major axis.

3.1 Lyapunov characteristic exponent

Evidently, the 9:16M is a high-order resonance (7th order) with apparently little strength (Murray & Dermott 2000). However, it has a great influence on the Datura orbits, forcing semi-major axis changes and sometimes moving asteroids from one side of the resonance to the other.

To understand the stability of the Datura family orbits, we employed the Lyapunov characteristic exponent (LCE) on the proper semi-major axis dependence (see e.g. Skokos 2010). LCE measures the mean exponential rate of the divergence of nearby trajectories. The data source that we used for our calculations is from the AstDys website (Knežević & Milani 2003) shown in Table 3².

The distributions of LCE built on analytic and synthetic values of the proper semi-major axis differ (see Fig. 2). Even so, a gap can easily be seen, in which relatively few points lie, and this gap is similar to the values of the 9:16M resonance borders.

The differences of the LCE values for the analytic and synthetic semi-major axis may be related to the difference in the time span used inthe calculations, which in synthetic proper elements is 2 Myr. This is greater than the time span used in analytic values. However, we note that for very young asteroid families such as the Datura family, analytic values appear to be more precise; the 9:16M resonance cause the gap in the distribution of the Datura family members by the proper semi-major axis, both in analytic and synthetic values. Some asteroids are on opposite sides of the resonance, meaning that a quantitative analytical study of 9:16M is difficult. Based on the results shown in Fig. 2, we can conclude that the orbits with very close values of their proper semi-major axes (with a difference of about 10⁻⁵ AU) display a large difference in their orbital stability (in Lyapunov time).

Fig. 1 Evolution of the orbit of (89309) 2001 VN36 depending on the initial variations in orbital elements. The dynamical evolutions of the orbit as a function of the eccentricity and inclination are shown in the upper graph, and the node and perihelion longitude variations are shown in the lower graph.

3.2 Determining the resonance position with the Euclid algorithm

To understand these dependences better, we applied the so-called Euclid algorithm of resonance (Subbotin 1968). For a simple mean motion resonance, $$j_1{n}_M+j_2{n}_a=0$$(1)

is valid, where n_M, n_a is the mean motion of Mars and the asteroid,respectively, and j₁, j₂ are the integers.

It has long been established that a long-period term of perturbation (so-called small divisors case) can have a noticeable amplitude only when j₁ and j₂ have low absolutevalues. The most straightforward way to find these values is to use the decomposition of the averaged mean motion ratio n_M ∕n_a into a continued fraction, $$\frac{{\overline{n}}_{\text{M}}}{{\overline{n}}_{\text{a}}}=q_0+\frac{1}{q_1+\frac{1}{q_2+\frac{1}{q_3+\frac{1}{\cdots \frac{1}{q_{k+1.}}}}}}.$$(2)

To determinate values of q_k, we used the Euclid algorithm, $$\begin{array}{ll}\frac{{\overline{n}}_M}{{\overline{n}}_a}\hfill & =\frac{q_0}{1},\hfill \\ \frac{{\overline{n}}_M}{{\overline{n}}_a}\hfill & =q_0+\frac{1}{q_1},\hfill \\ \frac{{\overline{n}}_M}{{\overline{n}}_a}\hfill & =q_0+\frac{1}{q_1+\frac{1}{q2}},\hfill \\ \cdots \hfill & \hfill \end{array}$$(3)

resulting in $$\frac{{\overline{n}}_{\text{M}}}{{\overline{n}}_{\text{a}}}=\frac{q_0}{1}=\frac{q_0{q}_1{q}_2+q_0+q_2}{q_1{q}_2+1}=\cdots =\frac{a_k}{b_k}.$$(4)

The period of perturbation (P) increases due to the approach of resonance, $$P=\frac{2\pi }{|b_k{\overline{n}}_{\text{M}}-a_k{\overline{n}}_{\text{a}}|},$$(5)

and it can be considered as an additional criterion to determine the resonance acting on a minor body. This conclusion is true for the pendulum model of libration at resonance: the period increases as a result of approaching the separatrix. The amplitude of perturbation in linear theory is proportional to $$P{\text{e}}^{|a_k-b_k|},$$(6)

however, it is not true in the non-linear approach because the period tends to infinity when an exact resonance is reached, but in reality, the amplitude always stays finite. Subbotin (1968) showed that when the motion of two planets is studied under the influence of their mutual attraction, it is only necessary to account for one small divisor.

The location of the 9:16M resonance as given by Vokrouhlický et al. (2017) is 2.2359 AU. When we apply the more precise value a_M = 1.523679338 AU for the semi-major axis of Mars, we reach a value of a_r = 2.23603 AU for the centre of this resonance and n_r = 1061.199 arcsec day⁻¹ for the mean motion. According to our numerical integration results, the centre of the chaotic zone (measured as half the sum of the upper and lower boundaries) is about 2.23585 AU. We calculate the distance from resonance (resonance Δ = n_a − n_r arcsec day⁻¹) and resonance period (Table 3).

According to our own findings, the stability can be better characterised by the resonance period and the distance of the body from the resonance than by the LCE in this particular case, as they are clearly affected by the chaotic nature of the dynamical evolution of 2001 VN₃ (Nesvorný et al. 2006). This fact is well supported by the results obtained in the numerical integrations. Otherwise, the LCE leads to the false impression that the dynamical behaviour of 60 151 1999 UZ₆ and (90265) 2003 CL₅ is more unstable than the dynamical evolution of 2001 VN₃₆.

Fig. 2 Lyapunov characteristic exponent (LCE) for the members of the Datura family. The analytic semi-major axis dependence is shown as red dots, and the synthetic semi-major axis dependence is shown as blue squares.

Fig. 3 Evolution of the semi-major axis of 2014 OE206 close to 50 kyr ago.

3.3 Resonance boundaries

The resonance boundaries given by Vokrouhlický et al. (2017) are 2.2355 < a_r < 2.2363 AU. They coincide with the unstable area in the analytic semi-major axis (Fig. 2). For the synthetic semi-major axis, we have the interval value of 2.2350 < a_r < 2.2357 AU (Fig. 2). In our previous work (Rosaev & Plávalová 2018), we gave an estimate of 2.23485 < a_r < 2.23505 AU for the unstable zone. Consistently, the value a_r = 2.23485 AU may be interpreted as the left boundary of the resonance in the synthetic semi-major axis. Rosaev & Plávalová (2018) discussed the jumps of some asteroids from one side of resonance to the other. A detailed picture of one of such jump is given in Fig. 3. Prior to the occurrence of the jump, the oscillation of the semi-major axis acquires a large amplitude. After showing two peaks, the value of the semi-major axis becomes more stable and reaches a new lower level. This effect may be crucial in our understanding of the process to reconstruct the origin of the Datura family. However, we focused on the details of these resonance jumps and tried to explain the reason for them.

It would be natural to assume that the position and width of the resonance changes with time. We tested several possible scenarios to account for this. First of all, we considered the main observational fact: all members of the Datura family were affected by strong perturbations in semi-major axis at the same two epochs, about 50 and 310 kyr ago (Fig. 4). This allows us to assume that the reason for the chaotic dynamical behaviour affecting the Datura family are the orbit variations of Mars. To study this problem, we analysed the orbit of Mars. We observed that the semi-major axis of Mars carries no noticeable secular variations in the considered time interval 0–800 kyr, and so this factor may be excluded. However, the inclination and proper frequency of the node longitude of Mars have strong variations. Because they occurred about 50 kyr and 310 kyr ago, they are the most logical reason for the jump in the values of the semi-major axes of the Datura family members.

During the interval − 800 kyr < t < −310 kyr, the value of the change rate in the node longitude of Mars changed by an average of 0.0048 deg yr⁻¹. However, in the interval −310 kyr < t < −50 kyr, we observe strong chaotic variations in the node rate of Mars; in some intervals, $\dot{\Omega }$ changed sign (Fig. 5). These variations closely correspond with the semi-major axis variations for some members of the Datura family.

Fig. 4 Dynamical evolution of the semi-major axis of Datura family members, illustrating jumps close to 50 and 310 kyr ago.

Fig. 5 Evolution of the node longitude (red line) and inclination (blue line) of the Mars orbit close to 60 kyr ago.

3.4 Comparison of the results of our calculations for the Mars orbit with results of previous calculations

To assess the reliability of the results obtained in this work, we compared our integration results with a classical Mars orbit eccentricity evolution (Laskar et al. 2002). The difference increased with time and reached 10% by the end of the integration interval (800 kyr). This might be explained by the fact that we neglected some perturbations in our calculations, such as those due to relativistic effects (Fig. 6). This discrepancy beyond 450 kyr is noticeable, but considering the age of the Datura family, we can ignore it, although the remaining data until present time (− 450–0 kyr) agree closely with the classical results of Laskar et al. (2002).

4 Simple theoretical model

The effect of splitting the resonance at a high precession rate is well known (see e.g. Morbidelli 2002, and references herein). In general, it leads to chaos due to overlapping neighbouring resonances (Chirikov 1979). However, the outcome is always unique and requires careful consideration in each case. The general form of a resonance argument is (Murray & Dermott 2000) $$\phi =j_1{\lambda }_M+j_2\lambda +j_3{\overline{\omega }}_M+j_4\overline{\omega }+j_5{\Omega }_M+j_6\Omega =0.$$(7)

Here λ_M, λ are the mean longitudes, ${\overline{\omega }}_{\text{M}},\overline{\omega }$ is the perihelion, and Ω_M, Ω are the node longitudes of Mars and the asteroid,respectively. j₁...j₆ are integer number coefficients and valid only if ${\sum }_{i=1}^6{j}_i=0$.

The reach of the exact 9:16M resonance takes place within a very short time interval. In Fig. 7, we present an example of the resonance argument (7) for the evolution of one of the clones of (485010) 2009 VS₁₁₆. In the interval around 320–300 kyr ago, when the jump of the semi-major axis took place, we have the libration of resonance argument. Soon after, in the interval 300–290 kyr ago, we have a chaotic behaviour as opposed to the behaviour found in other intervals, in which a circulation of the resonance argument can be seen.

As is well known (e.g. Murray & Dermott 2000; Morbidelli 2002), each mean motion resonance can be split into several resonances. For the so-called e-resonances and i-resonances, the arguments are $$\begin{array}{ll}\stackrel{.}{{\phi }_1}\hfill & =k{n}_{\text{M}}-pn-\left(k-p\right)\dot{\overline{\omega }},\hfill \\ \stackrel{.}{{\phi }_2}\hfill & =k{n}_{\text{M}}-pn-\left(k-p\right)\dot{\Omega },\hfill \\ \stackrel{.}{{\phi }_3}\hfill & =k{n}_{\text{M}}-pn-\left(k-p\right){\dot{\overline{\omega }}}_M,\hfill \\ \stackrel{.}{{\phi }_4}\hfill & =k{n}_{\text{M}}-pn-\left(k-p\right){\dot{\Omega }}_M,\hfill \end{array}$$(8)

in the general view, $$\begin{array}{ll}\stackrel{.}{{\phi }_i}\hfill & =k{n}_{\text{M}}-p{\alpha }^{-3/2}{n}_{\text{M}}-c\dot{\overline{\omega }}-\text{d}\dot{\Omega }-f{\dot{\Omega }}_{\text{M}}\hfill \\ \hfill & -\left(|k-p|-c-d-f\right){\dot{\overline{\omega }}}_M=0,\hfill \end{array}$$(9)

where k and p are integers and ∣ k − p∣ is the resonance order (in this paper, 16 − 9 = 7), $$\begin{array}{ll}c\hfill & =\mathrm{0...}k-p,\hfill \\ d\hfill & =\mathrm{0...}k-p-c,\hfill \\ f\hfill & =\mathrm{0...}k-p-c-d,\hfill \end{array}$$(10)

and α is represented as $$\begin{array}{ll}{\alpha }_{cdf}\hfill & =(\frac{p}{k}-\hfill \\ \hfill & {-\left[c\dot{\overline{\omega }}+d\dot{\Omega }+f{\dot{\Omega }}_M+\left(|k-p|-c-d-f\right){\dot{\overline{\omega }}}_{\text{M}}\right]/\left(k{n}_{\text{M}}\right))}^{-2/3},\hfill \end{array}$$(11)

and for the 16:9 resonance, $$\begin{array}{ll}{\alpha }_{cdf}\hfill & =(\frac{16}{9}\hfill \\ \hfill & {-\left[c\dot{\overline{\omega }}+d\dot{\Omega }+f{\dot{\Omega }}_{\text{M}}+\left(7-c-d-f\right){\dot{\overline{\omega }}}_{\text{M}}\right]/\left(9n_{\text{M}}\right))}^{-2/3},\hfill \end{array}$$(12)

where n and n_M represent the mean motions of an asteroid and Mars. Equation (11) is indicative in the case of arbitrary resonances. The derivatives in Eq. (11) can be calculated analytically or by numerical integration results. Here we use the second method to compute them.

The distances between these resonances are very small, therefore we can consider them as a single resonance in which the width (and position) varies slowly with time. The positions of resonances a_r = α_cdfa_M calculated by Eq. (12) are given in Table 4, where the common width of the resonance multiplet is given in the last column. When the rate of change of the node longitude of Mars ${\dot{\Omega }}_M$ changes at about 310 kyr and 50 kyr, the commonwidth and position of the centre of the resonance zone also change.

In the considered time interval, the averaged semi-major axis of Mars is 1.523679338 AU, and accordingly, the mean motion is 191.4089028 deg yr⁻¹. The proper frequencies for the Datura asteroid family are $\dot{\Omega }=-0.010272$ deg yr⁻¹, $\overline{\dot{\omega }}=0.009236$ deg yr⁻¹ and remain approximately constant, but the node precession rate of Mars is strongly unstable in the time interval 50–310 kyr ago (Fig. 5). However, in the remaining time interval from 310 kyr, this variation appears to be different, and we can represent it as an increasing line with a value that changes in circulation.

In the time interval −800 kyr < t < −310 kyr, the value of the change rate in the node longitude of Mars is about ${\dot{\Omega }}_{\text{M}}\approx -0.0049$ deg yr⁻¹, and the position of the centre of resonance varies with time until it reaches a value of 2.236029 AU. However, in the time interval −310 kyr < t < −50 kyr, the nodal rate of Mars varies strongly and chaotically. In some intervals where ${\dot{\Omega }}_{\text{M}}$ changes between positive and negative values, the average value of ${\dot{\Omega }}_{\text{M}}\approx 0$ and the resonance position is 2.236052 AU (Fig. 5). According to the results obtained in backward numerical integrations, approximately 50 kyr ago, we calculated the value ${\dot{\Omega }}_{\text{M}}=-0.0049$ deg yr⁻¹ and the centre of the resonance returned to 2.236029 AU. This shows a shift in the resonance position of about 0.00001 AU. Although apparently insignificant, this shift is sufficient to change the conditions of the dynamical interaction of some members of the Datura asteroid family that are in 9:16 M resonance. The centre of the resonance position is always greater than the integration value (2.2359 AU), that is, the chaotic resonance zone shifts towards the perturbing planet.

As a result of these shifts, the boundaries of the resonance change with time (within a range of about 0.00003 AU). This is comparable with the difference in proper semi-major axis between stable and unstable orbits of the Datura family (Fig. 2) and might initiate an approach to the resonance for asteroids of similar semi-major axis. Similar values of semi-major axis are clearly necessary, but this is not the only condition. According to the Yoder diagrammatic model of resonance crossing (Yoder 1979), we can say that because of the short periodic variations of asteroids and Mars semi-major axis, jumps can take place in a narrow range of the resonance argument phase. Simultaneously, a capture in resonance can likewise take place, but this is a rarity.

As an example, approximately 50 kyr ago, the semi-major axis of the centre of the resonance decreased and the orbit of 2014 OE₂₀₆ came close to the resonance. This was captured in resonance, but only for a very short time, about 2 kyr (see Fig. 3). After this, the orbit of 2014 OE₂₀₆ obtained a new value of the semi-major axis beyond that of the resonance; that is, 2014 OE₂₀₆ jumped from one side of the resonance to the other. For the same reason, the orbits of (89309) 2001 VN3 and 2014 OE₂₀₆ had to take time to adjust to the orbital changes of Mars.

In Table 4, the position of the centre of the resonance appears to be able to move away from the Sun as well as towards it. This means that each epoch of resonance position carries an instability in both directions of the jump, depending on the phase of the resonance.

Not all features of the dynamics of the Datura family members can be explained in this way. This is the case of the dynamical behaviour that corresponds to the intervals close to 220 kyr, where no jumps are found. This can be explained by taking into account that other orbital elements of Mars (mainly the orbital eccentricity) play a role in the complete theory. However, despite the inherent complexity of the resonance phenomenon in the complete theory, some relevant features of the orbits in the vicinity of the 9:16M resonance might be clearly understood with the help of the simple model used here.

Fig. 6 Evolution of the Mars orbit eccentricity using our results (red line) compared with the integration (blue line) by Laskar et al. (2002).

Fig. 7 9:16M resonant argument for a representative clone of (485010) 2009 VS116.

5 Effect of the large asteroid perturbations on the Datura family members

In our previous work (Rosaev & Plávalová 2017) devoted to the Datura asteroid family, all the results were obtained by integrating a model solar system that only included the eight major planets and Pluto. However, Carruba et al. (2016) and Novaković et al. (2016) have argued that the dwarf planet Ceres and the large main belt asteroids Vesta, Juno, and Pallas may have a remarkable effect on the gravitational perturbations affecting other smaller members of the main belt, such as the asteroids that belong to the Datura family.

Using these new data, our numerical integrations use a more complex solar system model that includes the dwarf planet Ceres and the large asteroids Vesta, Juno, and Pallas as additional massive perturbers. As a consequence of using this new model, our main results have changed significantly, in particular, we observe a behaviour similar to that of 2014 OE₂₀₆ in the semi-major axis of other members of the Datura family, namely 2001 VN₃₆ and 2011 KT₁₀. The Ceres and Vesta perturbations triggered an effect that destabilized the orbit of 2001 VN₃₆ (Fig. 8) and 2011 KT₁₀ (Fig. 9). We can conclude that the effect of Ceres on the dynamics of the Datura family is quite significant and cannot be neglected.

We note that Ceres and Vesta perturbations extend the boundaries of the resonance-related chaotic zone. As we stated in a previous paper (Rosaev & Plávalová 2017), 2014 OE₂₀₆ can jump from one side of the resonance to the other through the effect of planetary gravitational perturbations alone.

However, when we take the Ceres and Vesta perturbations into account, some other members of the Datura family exhibit a similar behaviour. In particular, 2001 VN₃₆ jumps from one side of 9:16M resonance to the other under perturbations from Ceres (Fig. 8). Similarly, asteroid 2011 KT₁₀ jumps through the resonance when we take perturbations by Ceres and Vesta into account (Fig. 9). We also note a similar effect in the semi-major axis of 2006 SD₃₈₂.

A possible explanation for the Ceres and Vesta effect is that the perturbations from Ceres and Vesta might slightly increase the mean semi-major axis of the Datura family asteroids, moving them closer to resonance. This shift in the values of the semi-major axis is sufficient to initiate strong resonance perturbations in their orbits.

The direct or indirect effect of large asteroids (Juno, Pallas, Hygiea, Eunomia, Interamnia, and Davida) on some Datura family members became evident when we compared the integration results. It is an indirect effect as a consequence of the orbital changes of Mars. The cumulative effect of Ceres, Vesta, Juno, and Pallas on the evolution of the orbital angular elements of Mars crowds out the effect of the other large asteroids (Hygiea, Eunomia, Interamnia, and Davida), making their effect pale into insignificance (see Fig. 10). We recall, however, that during a time of great instability in correlation with the orbital variations of Mars (310, 180, and 50 kyr ago), even the smallest, almost negligible perturbations could alter the character of the orbital evolution of the asteroids mentioned above (Fig. 8).

6 Yarkovsky effect on the resonance

Non-gravitational forces such as the Yarkovsky effect can significantly affect the interaction of asteroids with resonance. The Yarkovsky effect can change the value of the semi-major axis and shift an asteroid away from or towards the Sun. It can increase or decrease the semi-major axis of an asteroid, depending on the axial rotation of the asteroid.

Unfortunately, the parameters of the Yarkovsky effect for the smaller members of the Datura family are poorly understood. We can only estimate the maximum values. Applying the assumption that the rotation axis of all asteroids is colinear and the bulk density is p = 1.1 g cm⁻³, we can normalise the maximum effect using the data of Spoto et al. (2015) for asteroid Bennu. The details of the method we employed are described in our previous study (Kuznetsov et al. 2019), from which we obtained the maximum effect value for semi-major axis ȧ ≈ ± 5 × 10⁻⁴ AU Myr⁻¹ for the bodies we targeted.

These values are enough to move an asteroid inside or outside of the resonance area throughout the estimated lifetime of the Datura asteroid family (450 kyr). Our numerical integrations have shown that the jumps in the semi-major axis of 2001 VN₃₆ and 2009 VS₁₁₆ about 50 kyr ago displayed a wide range of assumed Yarkovsky parameters. We can conclude that studying the Yarkovsky effect is very important (potentially absolutely necessary) to understand the origin of the Datura family, but in order to provide a proper and more accurate study of this effect, new observations are required to improve the knowledge of the Yarkovsky parameters in this particular case.

Fig. 8 Evolution of the semi-major axis of 2001 VN36 in two different integration cases. The upper green line represents only the planetary perturbers (onlyplanet), and the lower red line additionally shows the perturbers Ceres, Vesta, Juno, and Pallas (the CVJP integration case).

Fig. 9 Evolution of the semi-major axis of 2011 KT10 according to three different integration cases. The red line represents only the planetary perturbers (onlyplanet), the blue line shows the CVJP integration case, and the green line illustrates the integration case CVJPHIDE with many perturbers.

Fig. 10 Differences in the orbit perihelion argument and the node longitude of Mars according to two different integration cases. The two lower lines represent the perihelion argument evolution, where the green line is for the CVJP case and the black line shows the CVJPHIDE case. The two upper lines represent the node longitude evolution. The blue line shows the CVJP case, and the red line illustrates the case of CVJPHIDE.

Fig. 11 Dynamical evolution of the semi-major axis of the pair Balam (bottom) and 2009 BR60 (top).

7 Resonance in the case of the pair 3749 Balam and 2009 BR60

Another very interesting example of the 9:16M resonance perturbation is that of the close asteroid pair (3749) Balam and (312497) 2009 BR₆₀ (Naoz et al. 2010). The primary asteroid of this pair is Balam, which has two satellites (Pravec et al. 2019) and a secondary asteroid 2009 BR₆₀ (Vokrouhlický 2009). Pravec et al. (2019) estimated the age of this pair to be about 400 kyr.

The semi-major axis of Balam (2.23649 AU) is close to the values of the Datura member semi-major axes. For this reason, we can expect an effect, however weak, of the 9:16M resonance on this object. The orbit of the secondary body 2009 BR₆₀ (a = 2.2367 AU) is more distant from this resonance, therefore its perturbations are expected to be weaker. In reality, the amplitude of the semi-major axis variations of Balam is about twice that of 2009 BR₆₀ (Fig. 11). In particular, perturbations that took place about 75 kyr and 310 kyr ago are related to Mars orbit variations. Not all of these features can be explained by Mars resonance perturbations, however, and as yet unknown resonances might act on both the asteroids of this pair, perhaps more so on the primary, Balam. The most probable is the 11:39J resonance with Jupiterat 2.23673 AU.

Figure 11 shows that the dynamical evolution of the semi-major axis of Balam displays two peaks at approximately 310 kyr and 75 kyr ago that have a possible relation to the 9:16M resonance. These peaks are negligible in the case of the secondary body. This can easily be explained by the fact that the location of 2009 BR₆₀ is too distant relative to the 9:16M.

8 Conclusion

Because of the widely accepted value for the age of the Datura family (450 kyr) and the sheer number of its family members, it is difficult to reconstruct a detailed dynamical history of this family. However, the Datura family provides a good case study for the interaction of a compact group of asteroids with a narrow mean motion resonance of high order (9:16). We have provided simple explanations for these interactions in some particular cases using a more accurate model of the solar system that includes not only the gravitational perturbations of the eight major planets of the solar system plus Pluto, but also the perturbations caused by the dwarf planet Ceres and the large asteroid Vesta. Additionally, we note the 9:16M perturbations in the orbit of Balam.

We have highlighted that the orbits with a very close proper semi-major axis (with a difference of about 10⁻⁵ AU) have great differences in their stability in the vicinity of high-order resonance. We further note that the distances between resonance and resonance period are more precise indicators for stability than the LCE, at least in single high-order resonance cases.

Moreover, we studied the significant role of the weaker Ceres and Vesta perturbations on the dynamics of the Datura family members and the pair Balam and 2009 BR₆₀ close to resonance. Although we considered the Yarkovsky effect very briefly in this paper, it is evidently necessary to conduct a thorough study to come to a clear understanding of this very complex phenomenon. The simple theoretical model presented in this paper is able to account for some examples of the most recent chaotic dynamical behaviour, such as resonance jumping, which affect the relatively small main belt asteroids. It can be used as a template on which to build in future work.

References

All Tables

All Figures

	Fig. 1 Evolution of the orbit of (89309) 2001 VN36 depending on the initial variations in orbital elements. The dynamical evolutions of the orbit as a function of the eccentricity and inclination are shown in the upper graph, and the node and perihelion longitude variations are shown in the lower graph.
In the text

	Fig. 2 Lyapunov characteristic exponent (LCE) for the members of the Datura family. The analytic semi-major axis dependence is shown as red dots, and the synthetic semi-major axis dependence is shown as blue squares.
In the text

	Fig. 3 Evolution of the semi-major axis of 2014 OE206 close to 50 kyr ago.
In the text

	Fig. 4 Dynamical evolution of the semi-major axis of Datura family members, illustrating jumps close to 50 and 310 kyr ago.
In the text

	Fig. 5 Evolution of the node longitude (red line) and inclination (blue line) of the Mars orbit close to 60 kyr ago.
In the text

	Fig. 6 Evolution of the Mars orbit eccentricity using our results (red line) compared with the integration (blue line) by Laskar et al. (2002).
In the text

	Fig. 7 9:16M resonant argument for a representative clone of (485010) 2009 VS116.
In the text

	Fig. 8 Evolution of the semi-major axis of 2001 VN36 in two different integration cases. The upper green line represents only the planetary perturbers (onlyplanet), and the lower red line additionally shows the perturbers Ceres, Vesta, Juno, and Pallas (the CVJP integration case).
In the text

	Fig. 9 Evolution of the semi-major axis of 2011 KT10 according to three different integration cases. The red line represents only the planetary perturbers (onlyplanet), the blue line shows the CVJP integration case, and the green line illustrates the integration case CVJPHIDE with many perturbers.
In the text

Fig. 10 Differences in the orbit perihelion argument and the node longitude of Mars according to two different integration cases. The two lower lines represent the perihelion argument evolution, where the green line is for the CVJP case and the black line shows the CVJPHIDE case. The two upper lines represent the node longitude evolution. The blue line shows the CVJP case, and the red line illustrates the case of CVJPHIDE.

In the text

	Fig. 11 Dynamical evolution of the semi-major axis of the pair Balam (bottom) and 2009 BR60 (top).
In the text