Free Access
Issue
A&A
Volume 585, January 2016
Article Number A56
Number of page(s) 10
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201526953
Published online 17 December 2015

© ESO, 2015

1. Introduction

A wealth of new data, such as Earth-based broadband photometry and space-based infrared observations, has allowed us to substantially refine our understanding of the inventory of asteroid families in the main belt (see, e.g., Masiero et al. 2015; Nesvorný et al. 2015, and references therein). Previous identification methods were limited in that they only used clustering techniques to find groups of asteroids in the three-dimensional space of proper orbital elements. Recognizing and characterizing individual, overlapping structures, though to be distinct families, in broad asteroid clusters were among the most significant achievements in using the new data.

As an example, it is interesting to discuss the case of the old-recognized Nysa-Polana-Hertha complex (see, e.g., Zappalà et al. 1995, Cellino et al. 2001, Milani et al. 2014, and references therein). Data from Sloan Digital Sky Survey (SDSS; Parker et al. 2008) and results from WISE mission (Masiero et al. 2011, 2013) allowed us to create a “high-definition version” of this orbital region. Previously, the available data hinted that the region was made up of two fundamental components (e.g., Cellino et al. 2001): (i) the high-albedo (S or E spectrally) part at high eccentricity (the Nysa or Hertha segment); and (ii) the low-albedo (B or C spectrally) part at lower eccentricity (the Polana segment). In this paper we deal with (ii) possibly one of the most important source zones of primitive low albedo near-Earth asteroids including 101955 Bennu, which will be visited by NASA’s OSIRIS-Rex mission in the near future (see, e.g., Campins et al. 2010, Bottke et al. 2015, Lauretta et al. 2015).

Further analysis of this region was enabled by including asteroid size data derived from absolute magnitude and albedo estimates. Using this approach, Walsh et al. (2013) and Bottke et al. (2015) argue that the former Polana component consisted of at least two separate families (see also Dykhuis & Greenberg 2015). The most populous and densest cluster has been named the Eulalia family because of its association with the ~40 km size, C-type main belt asteroid 495 Eulalia (e.g., Fieber-Beyer et al. 2008, 2012, for the spectral observation and classification). Interestingly, linking 495 Eulalia, presumably the largest fragment, to the rest of the family members has been challenging; 495 Eulalia lies on the extreme periphery of the family zone identified in proper element space by clustering techniques. This is because 495 Eulalia has a significantly lower proper eccentricity eP value than other family members (Fig. 1) and, at the same time, the high value of the proper semimajor axis aP compared to the family members. Both together prevent direct linkage of 495 Eulalia to its family using a simple clustering method in the proper element space.

thumbnail Fig. 1

Members of the Eulalia family from Walsh et al. (2013) in two-dimensional projections of the proper elements space: (i) semimajor axis aP vs. eccentricity eP (top); and (ii) semimajor axis aP vs. inclination IP (bottom). The star indicates location of 495 Eulalia, the suggested largest fragment in the family. Its very close proximity to the prominent mean motion resonance J3/1 (gray area) makes the proper eccentricity value unstable. On a timescale of tens to hundreds of Myr it can change by ~± 0.02, on occasions be aligned with the mean value in the family defined by members ejected onto safer orbits far from this resonance.

Open with DEXTER

However, Walsh et al. (2013) pointed out that 495 Eulalia resides near the prominent 3:1 mean motion resonance with Jupiter (J3/1) (see already Guillens et al. 2002; and Morbidelli et al. 1995, who discussed the role of mean motion resonances for structure of asteroid families including the Nysa-Polana-Hertha complex). As a result, the orbital structures associated with this resonance make the orbit of 495 Eulalia long-term unstable. On a timescale of tens to hundreds of Myr, Eulalia’s orbital eccentricity may exhibit large excursions. As a result, Walsh et al. (2013) propose that at the moment of collisional origin of the family (possibly some 0.751.1 Gyr ago, see also Bottke et al. 2015), 495 Eulalia was aligned in eP with the other members but subsequently experienced chaotic evolution in this element. Note also that, in spite of its chaotic nature, the dynamical lifetime of 495 Eulalia’s orbit was found to be long enough that it did not conflict with the proposed family’s age. Walsh et al. (2013) mention that only about 15% of close clones of 495 Eulalia fell into the J3/1 resonance over 1 Gyr of their simulation (similarly some 30% were eliminated at 2 Gyr; only at 4 Gyr did more than half of the clones reach the resonance).

While plausible, this scenario may not fully describe the complex nature of the past orbital evolution of 495 Eulalia. An unknown factor is the contribution of Yarkovsky thermal forces (e.g., Bottke et al. 2006). They could have brought this asteroid to its current orbit from some modest location inside the family that was originally more distant from the J3/1 resonance. If so, this might also have implications for the estimated age of the Eulalia family. The crucial missing factor in estimating the magnitude of Yarkovsky drift in semimajor axis is 495 Eulalia’s obliquity value ε. This is because for kilometer-size and larger asteroids the Yarkovsky effect is dominated by the diurnal variant for which da/dt ∝ cosε (e.g., Bottke et al. 2006). So far little was known about ε of Eulalia. Based on a very limited lightcurve dataset from 495 Eulalia apparitions in 1983 and 1984, Binzel (1987) obtained a rotation period of 29.2 ± 0.1 h and suggested it had a near ecliptic pole position at longitude and latitude (224°,2°). This would imply ε ≃ 89°. Admittedly, though, this result is very uncertain and only relies on a qualitative argument.

Table 1

Aspect data for dense photometry observations of 495 Eulalia.

This situation motivated us to obtain new lightcurve observations of 495 Eulalia. Our goal was to better constrain its current pole orientation (Sects. 2 and 3). Additionally, accessing the role of the Yarkovsky effect in the dynamical evolution of 495 Eulalia over many hundreds of Myr to Gyr timescales requires information about its past obliquity state on a comparably long interval of time. For that reason we also take a brief look at Eulalia’s spin dynamics (Sect. 3.1). Finally, in Sect. 4 we use the information about the 495 Eulalia’s spin state and show possible past evolutionary tracks of its orbit with respect to how it may have evolved via the Yarkovsky effect.

2. Observations

We complemented the archival lightcurve observations of 495 Eulalia from Binzel (1987) by obtaining dense photometric datasets in three extensive campaigns in 2012, 2014 and 2015 (Table 1). Additionally, we also made use of sparse photometric data from the US Naval Observatory that were downloaded from AstDyS website1. A total of 118 individual observations between 1998 and 2008 were used in our analysis.

The 2012 observations from the Palmer Divide Observatory were made with a 0.35-m Schmidt-Cassegrain telescope and FLI-1001E CCD camera with a 1024 × 1024 array of 24-micron pixels. Exposures were 240 s with a V filter. A total of 1181 observations were made from 2012 Nov 8 through Dec 23.

The 2014 observations from the Palmer Divide Observatory were made with a 0.30-m Schmidt-Cassegrain telescope ML-1001E CCD camera with a 1024 × 1024 array of 24-micron pixels. Exposures were 120 s with no filter. A total of 1439 observations were made from 2014 Mar 8 to Mar 25.

For both sets of observations, MPO Canopus was used to calibrate and measure the images using differential aperture photometry. The Comparison Star Selector in MPO Canopus was used to find comparison stars of near solar color. The magnitudes for the chosen stars were taken from the MPOSC3 catalog, which is a subset of the 2MASS catalog with BVRI magnitudes derived from formulas by Warner (2007).

The 2014 observations at the Palmer Divide Observatory covered the pre-opposition leg of Eulalia’s observability period. They were complemented by post-opposition observations at the Ondřejov observatory over 7 nights during May and June. The observations were taken with a 0.65-m telescope equipped with a Bessell R filter. They were calibrated in the Johnson-Cousins system using Landolt (1992) standard stars with absolute errors 0.01 mag. Integration times were 180 s, with the telescope tracked at the half-apparent rate of the asteroid. Because of the asteroid’s long rotation period we did not need to take continuous observations. Instead we took a short series of typically four images 13 times per night. The typical separation time between each image was one hour, depending on the scheduling constraints of our other asteroid observations; we usually worked Eulalia as a secondary target on our observation nights. Our observations were reduced using procedures described in Pravec et al. (2006). The resulting data points from each series (typically 4 points) were averaged to form a normal points, which were then used to model the spin and shape.

The 2015 observations from the Blue Mountains Observatory were done due to the favorable positioning of the asteroid in the southern hemisphere during this year’s apparition. The instrument used was a 0.35-m Schmidt Cassegrain telescope operating at f/5.9 and an SBIG ST8XME CCD camera with a 1530 × 1020 array of 9 micron square pixels for most of the observations. Additionally, the observations on June 28 were taken from a location near Perth, West Australia. The instrument used was a 0.30-m Schmidt Cassegrain telescope operating at f/7.2 and an SBIG ST8XME CCD camera. The images from both systems were taken without a filter.

All 2015 images were collected pre-opposition (later observations during this apparition were prevented by asteroid location near the galactic plane). All of our data was reduced using MPO Canopus using deferential aperture photometry. A comparison Star Selector method was used to find solar colored comparison stars. Their magnitudes were derived from the 2MASS catalog with BVRI magnitudes as described in Warner (2007).

thumbnail Fig. 2

Example of photometric data of 495 Eulalia (symbols) fitted with synthetic lightcurves (solid curves) produced by the model in Fig. 4. Because Eulalia’s rotation period is about 29 h, the data from one night cover always only a part of the rotation cycle. The viewing and illumination geometry for the pole solution P1 is given by the aspect angle θ, the solar aspect angle θ0, and the solar phase angle α.

Open with DEXTER

3. Lightcurve inversion

To reconstruct the shape, spin axis orientation, and the sidereal rotation period of 495 Eulalia, we applied the lightcurve inversion method of Kaasalainen et al. (2001) to the lightcurve data listed in Table 1 (see Fig. 2 for examples how the model matches the observations). Note that the nearly 3 year interval covered by our new observations in 2012 and 2015 provides enough baseline to link them in phase with the archival observations in 1983 and 1984. Additionally, the 10 year interval of the sparse photometry provided a valuable phase constraint. As a result, our solution leads to an unambiguous sidereal rotation period of 28.96589 ± 0.00007 hr and two possible pole orientations, (λ,β) = (38°,35°) (pole P1) and (λ,β) = (212°,46°) (pole P2). The obliquity of our two pole solutions is 54° and 45°, respectively.

To estimate the uncertainties of our pole directions, we used the same approach as in Vokrouhlický et al. (2011); we mapped the χ2 for all pole directions and set the boundary of acceptable solutions (Fig. 3). Interestingly, the pole ecliptic longitude was determined with an accuracy of about , while the pole latitude was more poorly constrained with an uncertainty of about 20°. Note the asymmetry of the uncertainty zone in latitude, extending toward higher latitude values. While a significant uncertainty is left in the ecliptic latitude of Eulalia’s rotation pole, we can exclude near ecliptic values, such as suggested by Binzel (1987). The maximum obliquity consistent with the data is 70°, but the true value is likely smaller. The nominal shape model corresponding to the P1 pole is shown in Fig. 4.

thumbnail Fig. 3

Dependence of the level of fit between the data and the model on the pole direction shown in sinusoidal projection of the sky in ecliptic coordinates. The reduced χ2 is color-coded. There are two local minima for poles P1 and P2 (white circles). The estimated uncertainty of the poles is shown as two solid boundaries around P1 and P2. The uncertainty of the pole latitude β is much larger than that of pole longitude λ.

Open with DEXTER

The first pole solution P1 is reminiscent to what is called a Slivan state (see, e.g., Vokrouhlický et al. 2003; Vraštil & Vokrouhlický 2015), namely the spin pole has been captured in a secular spin-orbit resonance related to the s6 planetary frequency. While not impossible, we find it unlikely. Using methods in Vraštil & Vokrouhlický (2015) we found that a stable capture into a Slivan state may only exist for Δ ≤ 0.13, where and (A,B,C) are principal moments of inertia (see Sect. 3.1). The lightcurve inversion technique provides in principle the asteroid shape (such as in Fig. 4) from which Δ could be estimated. However, given the limited amount of data, especially those which are absolutely calibrated, the resulting shape of 495 Eulalia is uncertain. The main degeneracy in our shape solution is in the polar flattening, which directly affects Δ. While the formal best fit solutions with P1 and P2 yield Δ = 0.34 and 0.39, respectively, we estimate (using different shape representations, creating non-convex models, and artificially rescaling the convex models along the rotation axis) that the realistic interval for shape solutions that still match available observational constraints span an interval of ~(0.15−0.42). Lower values of Δ would lead to a signal incompatible with the observations.

thumbnail Fig. 4

Shape model for 495 Eulalia reconstructed from the available lightcurves for the pole P1 solution. Left and middle parts yield two perpendicular equatorial views, right panel shows a view from the rotation pole. The majority of non-calibrated data constrain principally the equatorial axes ratio (which is small in this case because lightcurve amplitude never exceeds 0.5 mag). The polar stretching-factor is less constrained with the present dataset, but our tests indicate the asteroid is not near-spherical.

Open with DEXTER

3.1. Long-term pole evolution of 495 Eulalia

Observations provide the current value of Eulalia’s obliquity. However, in order to assess the role of Yarkovsky-driven migration in the orbit of this asteroid, it is necessary to estimate Eulalia’s obliquity behavior over a very long-timescale. This is not trivial, especially because Eulalia’s rotation state has been found to be prograde, for the following two reasons.

First, the prograde rotators among main-belt asteroids are known to have their obliquities affected by secular spin-orbit resonances (e.g., Vokrouhlický et al. 2006c; Vraštil & Vokrouhlický 2015). We shall not repeat the theoretical background of this phenomenon, and refer the interested reader to the previously mentioned references. Suffice to say that it concerns commensurability between spin axis precession due to solar gravitational torques and one of the modes in which the orbital plane precesses in inertial space. In the case of 495 Eulalia, the orbital precession is dominated by two terms, namely (i) the proper term with frequency s ≃ −46.8 arcsec/yr and inclination 2.5°; and (ii) the forced (planetary) term s6 ≃ −26.3 arcsec/yr and inclination . The spin axis precession may enter into a 1:1 resonance (hereafter called the Cassini resonance) with either of these two terms. If this happens, librations around the exact resonant state (called Cassini state C2) are associated with large-amplitude oscillations in obliquity. Because the polar precession rate is proportional to the dynamical flattening Δ of the body, we can plot obliquity of the Cassini state C2, and the width of the associated resonant zone, as a function of this parameter.

Figure 5 shows the result for both contributing terms, with s and s6 frequencies, in orbit precession as well as the current value of the rotation period 29 h. If Eulalia were near spherical, Δ ≤ 0.07, none of the resonances would be possible and the obliquity would exhibit only small oscillations. In the intermediate interval of values of small flattening, 0.07 ≤ Δ ≤ 0.13, only the resonance with the s6 frequency is possible. Specifically, for Δ ≃ 0.1 the resonance obliquity would be close to the current obliquity of Eulalia of 50°. Given the ecliptic longitude of 38°, the pole solution P1 from above would qualify to be a Slivan state (see Vraštil & Vokrouhlický 2015). More likely, though, the dynamical flattening satisfies Δ ≥ 0.13. In this case, the resonant zones of both orbital frequencies exist and overlap in obliquity. The obliquity exhibits large-scale chaotic oscillations within the interval of the overlapping zones. For instance, Fig. 5 makes us predict that if Δ ≃ 0.25, the obliquity would oscillate between 40° and 85°. Obviously, if the initial obliquity was lower than 40° or the dynamical flattening Δ large enough, the resonant phenomena would not be efficient and the obliquity would again have small-amplitude oscillations.

thumbnail Fig. 5

Obliquity of the Cassini equilibrium C2 and the maximum width of the associated resonance zone (gray area) as a function of the Eulalia’s dynamical flattening Δ. The solid line labeled C2(s6) corresponds to the planetary frequency s6, the line labeled C2(s) corresponds to the proper frequency s. For Δ ≥ 0.13 both resonant zones exist and overlap, producing chaotic obliquity evolution within the combined gray area (Fig. 6). These results assume the current rotation period P ≃ 29 h of Eulalia. Our pole solutions yield the current Eulalia’s obliquity ≃ (40°−60°), which could be affected by the chaotic zone unless very low or high value of Δ.

Open with DEXTER

thumbnail Fig. 6

Example of a long-term evolution of Eulalia’s obliquity ε. Here we assumed pole P1 orientation for the initial data and numerically integrated orbit and spin evolution for the next 10 Myr interval of time. The rotation period was assumed constant, and equal to its current value of 29 h, and the dynamical flattening set to Δ = 0.25. Bottom panel: obliquity ε (right ordinate) and/or cosε (left ordinate) vs time in the simulation. The obliquity exhibits large oscillations between ~45° and ~85° in two distinct regimes due to overlap of Cassini resonances related to the s and s6 orbital frequencies. Middle and top panels: Phase angle of the Cassini resonance of the s6 frequency (middle) and the s frequency (top). Intermittent periods of libration correlate with the range of the obliquity oscillations at the bottom panel.

Open with DEXTER

thumbnail Fig. 7

Same as in Fig. 6 but now for the dynamical flattening set to Δ = 0.35.

Open with DEXTER

In order to probe these dynamical effects in more detail, we numerically integrated Eulalia’s spin using the symplectic scheme of Breiter et al. (2005) (for methods see Vraštil & Vokrouhlický 2015). Figure 6 shows an example of Eulalia’s obliquity evolving over 10 Myr. The P1 pole solution was taken as an initial condition and the dynamical flattening parameter was set to Δ = 0.25 (bottom panel; we note that P2 pole solution provides a qualitatively similar result). Our result nicely confirms the behavior expected from Fig. 5, namely chaotic, large-amplitude oscillations due to intermittent interaction with the resonant zones about the Cassini state C2 of the s and s6 frequencies. For a higher value of the dynamical flattening parameter, Δ = 0.35, we expect less violent obliquity evolution. Indeed, Fig. 7 shows results from a 10 Myr simulation for the P1 pole initial state using this value of Δ. Here the oscillations of the obliquity are more regular and only have a small amplitude of ~.

So far we explored Eulalia’s obliquity behavior on a timescale of millions of years. In principle, though, we would like to assess evolution on even longer timescales up to Gyr. In spite of its large size, such timescales are long enough to make Eulalia’s spin state susceptible to tiny perturbations due to the Yarkovsky-O’Keefe-Radzievskii-Paddack (YORP) effect (e.g., Bottke et al. 2006). For example, using a suite of numerical runs for random asteroid shapes in Čapek & Vokrouhlický (2004), and appropriate scaling rules for size and bulk density (see the assumed values cited in Sect. 4.1), we estimate that the YORP effect could produce a long-term change in Eulalia’s rotation rate ω of the order | dω/ dt | ≃ 10-7 rad/s/Myr. Assuming this value, we note that the initial rotation period P of Eulalia could have been ~(10−15) h some Gyr ago, with the value slowly increasing to its current value. Moreover, the results shown in Fig. 5 are invariant with respect to product Δ P. Therefore, if the initial rotation period was smaller than its current value P, and asteroid’s shape is approximately preserved, the initial obliquities of the Cassini state would have been those of an effective Δ value equal to Δ (P′/P). Hence, it is conceivable that the initial Eulalia’s obliquity was once considerably lower. As the asteroid spun down, it could make the obliquity adiabatically evolve together with the Cassini state C2, being first temporarily captured in a Slivan state before entering the large chaotic zone due to overlap of the Cassini resonances of the s and s6 resonances. An example of this can be seen in Fig. 6. Later on, when the spin rate had considerably decreased, the rotation state could have dropped out of the chaotic zone, with the spin state becoming more regular as seen in Fig. 7. Further examples of such evolutionary paths are given in Vokrouhlický et al. (2003) and Vraštil & Vokrouhlický (2015). However, the possibility that Eulalia’s spin state interacted with the chaotic zone of the s and s6 resonances overlap does not permit us to deterministically infer details of this possible scenario.

Some take away messages from the above discussion are that: (i) Eulalia was likely a prograde rotator over its approximately Gyr-lasting existence; and (ii) the long-term average obliquity of 60° (or equivalently cosε ≃ 0.5) is quite plausible, with a possibility that this value was even lower during the first few hundreds of Myr of its existence.

4. Implications for Eulalia family

In this section we make use of the results above to infer aspects of Eulalia’s orbital history over the past ~1 Gyr.

First, we find that 495 Eulalia is qualitatively reminiscent of Koronis member 2953 Vysheslavia (see, e.g., Milani & Farinella 1995, Vokrouhlický et al. 2001). Recall that Vysheslavia was found residing in a chaotic region very close to the J5/2 mean motion resonance with Jupiter. Unlike the Eulalia case, however, Vysheslavia’s dynamical lifetime was found to be much shorter than the estimated age of the hosting Koronis family (Milani & Farinella 1995). Rather than assuming Vysheslavia is a large fragment from a recent collision, Vokrouhlický et al. (2001) proposed that Vysheslavia was transported to its current orbit by Yarkovsky drift from a stable orbit within the Koronis family. In this case the body originally had a semimajor axis slightly higher than the J5/2 resonance, such that it must have migrated inward. This requires the asteroid to have a retrograde rotation. The prediction ε> 90° for 2953 Vysheslavia has been confirmed by Vokrouhlický et al. (2006b). Even though the exact obliquity evolution of 495 Eulalia is unknown over long timescales, the strongest likelihood is that it was always lower than 90°. As a long term prograde rotator, 495 Eulalia has steadily migrated outward, i.e. towards the J3/1 resonance. Note that if Eulalia a had retrograde rotation, it would evolve away from the J3/1 resonance.

4.1. The expected drift-rate in semimajor axis

Here we estimate how much Eulalia could have migrated in semimajor axis over the past 1 Gyr, the approximate age of the Eulalia family (see Bottke et al. 2015). In order to obtain this result, we combine information about possible past obliquity evolutionary pathways from Sect. 3.1 with additional data about the asteroid.

According to measurements from the WISE spacecraft and subsequent data analysis, 495 Eulalia is a D = 40.0 ± 0.5 km asteroid with a geometric albedo pV = 0.05 ± 0.01 (e.g., Masiero et al. 2011). This solution assumes the asteroid has an absolute magnitude H = 10.8 obtained from astrometric data analysis, a value which is only slightly different from our more accurate result H = 10.95 ± 0.03 (see Warner 2013). According to the methodology in Pravec et al. (2012), such a revision of Eulalia’s absolute magnitude would mainly produce a small change in its albedo. This difference in absolute magnitude H produces a tiny correction that is within the stated uncertainty intervals. In fact, as also discussed by Pravec et al. (2012), the above mentioned 0.5 km uncertainty in Eulalia’s size and 0.01 uncertainty of its geometric albedo are likely formal, with a realistic value as large as 24 km and 0.03. While additional observation techniques, such as polarimetry may help decreasing these realistic uncertainty values (see, e.g., Cellino et al. 2015), these details have luckily no meaningful impact on our general arguments below.

Next, we must assume some value for the thermal inertia Γ of Eulalia’s surface since no direct measurement of this parameter is presently available. According to Delbò et al. (2007) and Delbò & Tanga (2009), a likely value for a main belt body of Eulalia’s size is Γ ≃ 100 J/m2/s1/2/K (SI units), though a lower value may not be excluded as well (see, e.g., Hanuš et al. 2014). Assuming the surface specific heat at constant volume C ≃ 800 J/kg/K and surface density ρs ≃ 1.7 g/cm3, both plausible for Eulalia’s primitive carbonaceous chondrite-like composition, we would obtain a surface thermal conductivity K = Γ2/ (ρsC) ≃ 0.007 W/m/K. As mentioned above, somewhat lower values down to K ≃ 0.001 W/m/K may be also possible. Finally, we assume the bulk density for Eulalia is ρb ≃ 1.3 g/cm3. This value is the same as compositionally similar asteroids like the 60 km diameter main belt asteroid 253 Mathilde (e.g., Yeomans et al. 1997) and the much smaller 0.5 km asteroid 101955 Bennu (e.g., Chesley et al. 2014). Note that Bottke et al. (2015) argue that Bennu originated in the Eulalia family. Should, however, Eulalia’s bulk density be slightly higher one could easily recalibrate our results using the relationship da/d.

Using the analytical formulation provided by Vokrouhlický (1999), we estimated the orbit-averaged change for Eulalia’s semimajor axis (da/dt) due to the Yarkovsky effect for different obliquity and rotation period values. Results shown in Fig. 8 assume Eulalia had its current rotation period of ~29 h for the length of the integration run. Should the rotation period be smaller, which is likely, the results would be simply shifted as indicated by the arrows in Fig. 8. This would make the da/dt value larger than today over the relevant interval of low values in thermal conductivity. Because the diurnal variant of the Yarkovsky effect dominates over the seasonal variant, da/dt scales with cosε, the reference for this evolution being the ε = 0° result (black solid curve). The gray curves show results for ε = 30° and ε = 60°, of which the latter is directly relevant for Eulalia (see Figs. 6 and 7). We conclude that the mean value of da/dt over the past Gyr for Eulalia could have been 5 × 10-6 au/Myr to 7 × 10-6 au/Myr. This suggests the accumulated semimajor axis change was 0.005 au to 0.007 au in 1 Gyr. This value is much smaller than the full spread of the family (Fig. 1), but it is still substantial given the close proximity of 495 Eulalia to the J3/1 resonance.

thumbnail Fig. 8

Estimated mean drift rate da/dt in semimajor axis due to the Yarkovsky effect for 495 Eulalia as a function of the surface thermal conductivity K. We assumed D = 40 km size, P = 29 h rotation period, 1.3 g/cm3 bulk density, 1.7 g/cm3 surface density and 800 J/kg/K. Since the diurnal variant of the Yarkovsky effect dominates, da/dt ∝ cosε, where ε is the obliquity of the spin axis. The upper black curve corresponds to a maximum drift-rate at ε = 0°, the gray curves are for ε = 30° and ε = 60°. Main-belt asteroids of Eulalia size and spectral type are expected to have K in between 0.001 and 0.01 W/m/K. Shorter rotation period would make the curves shifted according to the arrows (examples for 10 and 20 h).

Open with DEXTER

thumbnail Fig. 9

Evolutionary tracks of 100 synthetic D = 10 km asteroids started near the J3/1 resonance in the zone of Eulalia family. Their mean orbital elements were computed in a 5 Myr wide running window with 0.25 Myr shift: mean a vs. mean e (top), mean a vs. mean sinI (bottom). The simulation covered roughly 1 Gyr interval of time until the particles were eliminated from the simulation. Test bodies were given Yarkovsky drift-rate corresponding to the body of the chosen size and an effective obliquity 60° (Sect. 3.1). Initial orbits are highlighted by the upper and left box (the arrow recalls sense of migration in a). The bottom and right box shows current location of 495 Eulalia and its close clones (see the main text).

Open with DEXTER

thumbnail Fig. 10

Same as in Fig. 9 but now a vs time (top), and mean e vs. time (bottom). The mean semimajor axis of 495 Eulalia is 2.487 au, the mean eccentricity of 495 Eulalia is 0.129, both at the center of the respective zones delimited by black lines (their limits computed from excursions of the mean orbital elements of Eulalia’s clones as described in the text). Recall that the particles used in the simulation had 4 times smaller size than 495 Eulalia (10 km vs. 40 km); this means that the abscissa scale needs to be multiplied by a factor of 4 to interpret these results as potential past evolutionary tracks of 495 Eulalia.

Open with DEXTER

4.2. Possible role of the Yarkovsky migration of 495 Eulalia toward the J3/1 resonance

Finally, we conducted numerical simulations that followed potential evolutionary tracks for Eulalia since its formation. As discussed above, our working hypothesis is that Eulalia’s initial orbit was located inside the family and that this largest remnant experienced slow migration toward the J3/1 resonance. Given our unknowns and the chaotic nature of orbital evolution, our runs should be considered as an exploration of this concept rather than a true reconstruction of Eulalia’s past orbital history. Eulalia’s past obliquity evolution, directly linked to the strength of the Yarkovsky effect (da/dt), is intrinsically uncertain due to the chaotic phenomena discussed in Sect. 3.1. Moreover, orbital evolution near the J3/1 resonance is not predictable in a deterministic way over the very long timescale we use here.

Given the expected drift of Eulalia from Sect. 4.1, we constructed the initial orbital data using the following method. First, we decreased Eulalia’s osculating semimajor axis to the interval of values 2.4765 au to 2.4815 au in a random fashion in order to create 100 “clone” possibilities (compare with the current value 2.487 au). Second, we increased the osculating eccentricity by 0.025 to align the initial position of the simulated Eulalia asteroids with Eulalia family members. Indeed, we verified that the proper eccentricity in the first 10 Myr of the simulation was 0.155 for all the clones (compare with Fig. 1). The osculating inclination, as well as the secular angles and longitude in orbit, were not changed. All elements were set to correspond to the initial epoch MJD 57 000. The initial state vector of the planets were taken from the JPL ephemerides file DE405.

Our integrations were performed using the well documented and widely used orbit integration package 2 swift with the Yarkovsky effect included. For the sake of simplicity, we used a simplified variant from Vokrouhlický & Nesvorný (2008) where the effects of the thermal forces were modeled as a transverse acceleration; this produced a secular change in semimajor axis. In order to speed up the calculations, we set the effective size of the simulated particles to D = 10 km. This means that the semimajor axis drift speed was set to be about four times faster than shown in Fig. 8 for Eulalia. While the evolutionary tracks in orbital element space are preserved, the timescale needed to accumulate these effects is 4 times shorter. This computational shortcut has been used often in the previous literature (e.g., Bottke et al. 2001). The propagation timestep was 5 days. We let the integration proceed until all particles were eliminated from the simulation by falling into the J3/1 resonance.

In order to reduce the output dataset, we computed online mean orbital elements – semimajor axis, eccentricity and sine of inclination – using a 5 Myr wide running window. Since our goal was to approximately reconstruct the current orbit of 495 Eulalia, we must adopt some operational way to define it. To that purpose we considered results from the long-term integration of 495 Eulalia and its orbital clones from Walsh et al. (2013). We computed their mean elements in the first 50 Myr interval of our integration time and observed the limits in which they oscillate. We obtained the following intervals of values: (i) mean semimajor axis ranges between 2.4860 au and 2.4875 au; (ii) mean eccentricity ranges between 0.12 and 0.14; (iii) mean sine of inclination ranges between 0.04 and 0.05. The range for eccentricity was found to be rather wide because the bodies interact with the non-linear secular resonance g + g5−2g6. These limits define a box in the space of our mean orbital elements that vaguely characterizes the analogs we have created of current Eulalia’s orbit.

Our results are shown in Figs. 9 and 10. While many of the 100 orbits were scattered by the interaction with the chaotic region near the J3/1 resonance, some made it to the vicinity of the current Eulalia’s location in about 300500 Myr (given the factor of 4 mentioned above, this would translate to a real evolution timescale of ~1.22 Gyr). Obviously, starting the objects slightly closer to the resonance along their path, the final time at which they reach the current orbit of 495 Eulalia would have been smaller. We estimate that about 1015% of our objects eventually evolved into an orbital state similar to that of the observed Eulalia (the box in the mean orbital element space defined above). In this scenario, Eulalia’s orbit is dynamically long-term unstable but could remain near the current state for several hundreds of Myr. Eventually, it will evolve to the J3/1 resonance zone and be ejected from the main belt. This will make it, for a short period of time, a giant member of the near-Earth population (see Guillens et al. 2002, for further examples of large asteroids on unstable orbits near the J3/1 resonance that will evolve into the planet-crossing zone in the next 100 Myr).

5. Conclusions

The orbit of 495 Eulalia is long-term dynamically unstable. This is true in the gravitational-only model (see, e.g., Walsh et al. 2013), but it manifests itself on even shorter timescales when the Yarkovsky effect has been taken into account. Given the results of our simulation in Sect. 4, we estimate that Eulalia will fall into the J3/1 resonance within ~(0.21) Gyr. At that moment, the Eulalia family will become “orphan” it will lose its largest member. This makes the case that this family is rather special. Note there are a number of asteroid families near prominent mean motion resonances (e.g., Morbidelli et al. 1995, Nesvorný et al. 2015), but in these cases the largest fragment is preserved at a safe distance from them; only a portion of the smaller faster drifting members eventually leak out of the main belt through a resonances.

If Eulalia was initially further from the J3/1 resonance, as advocated in Sect. 4, the estimated formation age of the family may need to be revisited. This is because the approach used by Bottke et al. (2015), based on methodology in Vokrouhlický et al. (2006a), requires the center of the family to be defined roughly by the largest fragment. So far the issue of Eulalia’s semimajor axis change by the Yarkovsky effect has not been considered in such studies, presumably because the largest member was assumed to be too large to undergo substantial change in semimajor axis. In the case of ancient families, however, it is possible that this assumption may be invalid. Additionally, this next level generation of approach to study ages of asteroid families may explore effects of other possibly relevant processes such as sub-catastrophic collisions between the asteroids (e.g., Dell’Oro & Cellino 2007). However, an exploration of these interesting topics is beyond the scope of this paper.


2

http:// [2]www.boulder. [2]swri.edu/ [2]h˜al/ [2]swift.html

Acknowledgments

We thank Alberto Cellino whose comments helped to improve the final version of this paper. This work was supported by the Czech Science Foundation (grants P209-12-0229, GA13-01308S and GA15-04816S). Research funds for William Bottke were provided by NASA’s OSIRIS-REx mission; NASA Contract NNM10AA11C (D. S. Lauretta, PI).

References

All Tables

Table 1

Aspect data for dense photometry observations of 495 Eulalia.

All Figures

thumbnail Fig. 1

Members of the Eulalia family from Walsh et al. (2013) in two-dimensional projections of the proper elements space: (i) semimajor axis aP vs. eccentricity eP (top); and (ii) semimajor axis aP vs. inclination IP (bottom). The star indicates location of 495 Eulalia, the suggested largest fragment in the family. Its very close proximity to the prominent mean motion resonance J3/1 (gray area) makes the proper eccentricity value unstable. On a timescale of tens to hundreds of Myr it can change by ~± 0.02, on occasions be aligned with the mean value in the family defined by members ejected onto safer orbits far from this resonance.

Open with DEXTER
In the text
thumbnail Fig. 2

Example of photometric data of 495 Eulalia (symbols) fitted with synthetic lightcurves (solid curves) produced by the model in Fig. 4. Because Eulalia’s rotation period is about 29 h, the data from one night cover always only a part of the rotation cycle. The viewing and illumination geometry for the pole solution P1 is given by the aspect angle θ, the solar aspect angle θ0, and the solar phase angle α.

Open with DEXTER
In the text
thumbnail Fig. 3

Dependence of the level of fit between the data and the model on the pole direction shown in sinusoidal projection of the sky in ecliptic coordinates. The reduced χ2 is color-coded. There are two local minima for poles P1 and P2 (white circles). The estimated uncertainty of the poles is shown as two solid boundaries around P1 and P2. The uncertainty of the pole latitude β is much larger than that of pole longitude λ.

Open with DEXTER
In the text
thumbnail Fig. 4

Shape model for 495 Eulalia reconstructed from the available lightcurves for the pole P1 solution. Left and middle parts yield two perpendicular equatorial views, right panel shows a view from the rotation pole. The majority of non-calibrated data constrain principally the equatorial axes ratio (which is small in this case because lightcurve amplitude never exceeds 0.5 mag). The polar stretching-factor is less constrained with the present dataset, but our tests indicate the asteroid is not near-spherical.

Open with DEXTER
In the text
thumbnail Fig. 5

Obliquity of the Cassini equilibrium C2 and the maximum width of the associated resonance zone (gray area) as a function of the Eulalia’s dynamical flattening Δ. The solid line labeled C2(s6) corresponds to the planetary frequency s6, the line labeled C2(s) corresponds to the proper frequency s. For Δ ≥ 0.13 both resonant zones exist and overlap, producing chaotic obliquity evolution within the combined gray area (Fig. 6). These results assume the current rotation period P ≃ 29 h of Eulalia. Our pole solutions yield the current Eulalia’s obliquity ≃ (40°−60°), which could be affected by the chaotic zone unless very low or high value of Δ.

Open with DEXTER
In the text
thumbnail Fig. 6

Example of a long-term evolution of Eulalia’s obliquity ε. Here we assumed pole P1 orientation for the initial data and numerically integrated orbit and spin evolution for the next 10 Myr interval of time. The rotation period was assumed constant, and equal to its current value of 29 h, and the dynamical flattening set to Δ = 0.25. Bottom panel: obliquity ε (right ordinate) and/or cosε (left ordinate) vs time in the simulation. The obliquity exhibits large oscillations between ~45° and ~85° in two distinct regimes due to overlap of Cassini resonances related to the s and s6 orbital frequencies. Middle and top panels: Phase angle of the Cassini resonance of the s6 frequency (middle) and the s frequency (top). Intermittent periods of libration correlate with the range of the obliquity oscillations at the bottom panel.

Open with DEXTER
In the text
thumbnail Fig. 7

Same as in Fig. 6 but now for the dynamical flattening set to Δ = 0.35.

Open with DEXTER
In the text
thumbnail Fig. 8

Estimated mean drift rate da/dt in semimajor axis due to the Yarkovsky effect for 495 Eulalia as a function of the surface thermal conductivity K. We assumed D = 40 km size, P = 29 h rotation period, 1.3 g/cm3 bulk density, 1.7 g/cm3 surface density and 800 J/kg/K. Since the diurnal variant of the Yarkovsky effect dominates, da/dt ∝ cosε, where ε is the obliquity of the spin axis. The upper black curve corresponds to a maximum drift-rate at ε = 0°, the gray curves are for ε = 30° and ε = 60°. Main-belt asteroids of Eulalia size and spectral type are expected to have K in between 0.001 and 0.01 W/m/K. Shorter rotation period would make the curves shifted according to the arrows (examples for 10 and 20 h).

Open with DEXTER
In the text
thumbnail Fig. 9

Evolutionary tracks of 100 synthetic D = 10 km asteroids started near the J3/1 resonance in the zone of Eulalia family. Their mean orbital elements were computed in a 5 Myr wide running window with 0.25 Myr shift: mean a vs. mean e (top), mean a vs. mean sinI (bottom). The simulation covered roughly 1 Gyr interval of time until the particles were eliminated from the simulation. Test bodies were given Yarkovsky drift-rate corresponding to the body of the chosen size and an effective obliquity 60° (Sect. 3.1). Initial orbits are highlighted by the upper and left box (the arrow recalls sense of migration in a). The bottom and right box shows current location of 495 Eulalia and its close clones (see the main text).

Open with DEXTER
In the text
thumbnail Fig. 10

Same as in Fig. 9 but now a vs time (top), and mean e vs. time (bottom). The mean semimajor axis of 495 Eulalia is 2.487 au, the mean eccentricity of 495 Eulalia is 0.129, both at the center of the respective zones delimited by black lines (their limits computed from excursions of the mean orbital elements of Eulalia’s clones as described in the text). Recall that the particles used in the simulation had 4 times smaller size than 495 Eulalia (10 km vs. 40 km); this means that the abscissa scale needs to be multiplied by a factor of 4 to interpret these results as potential past evolutionary tracks of 495 Eulalia.

Open with DEXTER
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.