THE SCHULHOF FAMILY: SOLVING THE AGE PUZZLE

The data set consisted of about 20 light curves (dense photometry) observed during four apparitions and sparse-in-time photometry from the Catalina Sky Survey in the R and V filters that we downloaded from the AstDyS site.¹¹ Table 1 summarizes the available data. Observations in the 2012 apparition were numerous and they extended over several lunations before and after the opposition. This allowed us to unambiguously phase-link the observations from other apparitions and thus obtain an unambiguous solution of the sidereal rotation period for this asteroid.

We took photometric observations of (2384) Schulhof with the 0.65 m telescope in Ondřejov over six nights in 2012 and one night in 2013. The observations were taken with the Bessell R filter and they were calibrated in the Cousins R photometric system using Landolt (1992) standard stars with absolute errors of 0.01 mag. Integration times were 180 s and the telescope was tracked at the half-apparent rate of the asteroid. The observations were reduced using procedures described in Pravec et al. (2006).

Observations at Sugarloaf Mountain Observatory were made using a 0.5 m, f/4.0 reflector. The imaging CCD was a SBIG ST-10XME cooled to −15 C. Images with integration times of 110 seconds were taken through a clear filter. Absolute magnitudes were estimated using a method inherent in the analysis software, which was MPO Canopus. This calibration method is based on referencing a hybrid star catalog consisting mostly of 2MASS stars in the V band. The accuracy is estimated to be 0.07 mag.

The photometric observation of (2384) Schulhof was also obtained on 2012 September 9 with the 0.61 m f/4.3 reflector at the Skalnaté Pleso Observatory through the Cousins R filter and SBIG ST-10XME with 3 × 3 binning with resolution of 1.6 arcsec px⁻¹. 240 s CCD frames were reduced in the standard way using bias, dark, and flat field frames.

We also observed (2384) Schulhof at Abastumani Observatory using the 0.7 m Maksutov telescope over four nights during the 2012–2013 opposition and two nights during the 2013–2014 opposition. Images were taken through a clear filter and reduced using techniques explained in Krugly et al. (2002). The photometric measurements were performed with ASTPHOT software (Mottola et al. 1995). The estimated uncertainty of the photometric points was better than 0.02 mag.

Finally, observations of (2384) Schulhof were also made with the 0.41 m PROMPT-1 telescope at Cerro Tololo on one night in 2015. The observations were taken with a LUM (IR Block) filter with an exposure time of 180 s. The observations were reduced with standard aperture photometry procedures using the photometry program MIRA.

Our well-calibrated observations covered a large enough interval of phase angles to determine Schulhof's absolute magnitude H. From the Ondřejov observations, we derived the mean absolute magnitude in the Cousins R band H_R = 11.64 ± 0.03 and slope factor G = 0.24 ± 0.02. Considering the mean $(V-R)=0.49\pm 0.05$ value for S-type asteroids (e.g., Pravec et al. 2012) and recalling that (2384) Schulhof has measured SDSS colors matching those spectral classifications, we obtain H = 12.13 ± 0.06. This is significantly different from the value of 11.7 quoted above and used by the WISE team to determine Schulhof's size D and geometric albedo p_V.¹² Using the method presented in Section 4 of Pravec et al. (2012), we can improve the WISE solution of D and p_V. In particular, we obtained D = 11.57 ± 0.14 km and ${p}_{V}\;=0.19\pm 0.04$. As explained in Pravec et al. (2012; see also Harris & Harris 1997), a more significant correction concerns the geometric albedo value. We note that both uncertainties are formal. Their realistic values may be somewhat larger (especially as far as the size is concerned).

The size of (2384) Schulhof allows us to estimate two important quantitative parameters needed for the analysis of past convergence of Schulhof family members (Section 4): (i) the escape velocity from this body ≃7 m s⁻¹ and (ii) the Hill radius of its gravitational influence ≃2800 km (both assume 2 g cm⁻³ bulk density). Their values for the parent body of the Schulhof family were only slightly larger. Since we do not have enough information about the smaller family members yet, we shall assume ≃7 m s⁻¹ and ≃3000 km for the parent object escape velocity and Hill radius.

To reconstruct the physical model of Schulhof, we applied the light curve inversion method of Kaasalainen et al. (2001) to the photometric data listed in Table 1; see Figure 1 for the match between the sample of the observations and the model. We found a unique solution for the sidereal rotation period of 3.293677 ± 0.000002 hr for which the pole direction converged to two local minima of the χ² function at ecliptic longitude and latitude (λ, β) = (173°, −81°) (pole P1) and (λ, β) = (65°, −54°) (pole P2). To estimate the uncertainties of these pole directions (Figure 2), we mapped the χ² function for all pole directions on the sky and set the boundary of acceptable solutions in accord with the method described in Vokrouhlický et al. (2011). For both pole values, the region of acceptable solutions was more elongated along the β coordinate, which means that the longitude was determined better than the latitude. Because the available data set is still not very large and the number of adjusted parameters is significant, the uncertainty zone does not have an ellipsoidal shape. Rather, it is irregular and in the case of the P1 pole not centered at the formal best solution. Nevertheless, we note that the maximum southern latitude statistically compatible with the data is ≃ −40°. This corresponds to a minimum obliquity value of ≃135°. The shape model corresponding to the P1 pole is shown in Figure 3. The P2 shape model has a similar pole-on silhouette but is less flattened in the rotation-axis direction. Both shape models, as well as the complete set of dense light curve data used in our analysis, are available at the DAMIT website.¹³

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to assess the structure of the Schulhof family in the space of orbital elements, we computed synthetic proper values for semimajor axis, eccentricity, and sine of inclination. While our method is closely related to the concept introduced by Knežević et al. (2002), we slightly modified the procedure to suit our task. This is because we also want to provide proper elements for less accurate orbits, including some of the single-opposition members in Table 2 (note the AstDyS group computes synthetic proper elements for the numbered and multi-opposition asteroids only). Our procedure is as follows.

For each of the asteroids in Table 2, except for 2008 GW33 whose orbit is the least accurate, we first created 100 clone versions in the six-dimensional space of equinoctical elements ${\boldsymbol{E}}$ at the initial epoch MJD 57200. We obtained them using the probability distribution $p({\boldsymbol{E}})$ ∝ $\mathrm{exp}\left(-\frac{1}{2}\;{\rm{\Delta }}{\boldsymbol{E}}\cdot {\boldsymbol{\Sigma }}\quad \cdot \quad {\rm{\Delta }}{\boldsymbol{E}}\right)$, where (i) Δ${\boldsymbol{E}}$ = ${\boldsymbol{E}}$ − ${{\boldsymbol{E}}}^{\star }$ is the difference with respect to the best-fit orbital values ${\boldsymbol{E}}$^⋆ (Table 2) and (ii) Σ is the covariance matrix of the orbital solution downloaded from the AstDyS Web site. For each of the orbital representations (clones) we also stored its weight given by the $p({\boldsymbol{E}})$ value. Planetary state vectors and velocities for the same epoch were obtained from the JPL DE 405 ephemerides file. Each of the 100 orbital realizations for all asteroids were numerically integrated, together with planets, using a well documented and tested package swift.¹⁵ The integration time step was 5 days and the force model contained gravitational effects from the Sun and planets only. We tracked all orbits for 20 Myr with an output density of 50 years. The output data were then split into five 6.6 Myr long segments uniformly distributed over the integrated time interval (such that the initial time of the first segment was 0 Myr and the last time of the last segment was 20 Myr). The length of the segments was chosen to accommodate 2¹⁷ output points and to be long enough to allow good resolution of the lowest planetary frequencies g₈ and s₈. At each of the segments we (i) computed the mean value of the semimajor axis and (ii) determined the Fourier representation of the non-singular orbital elements $z=e\;\mathrm{exp}({\imath }\varpi )$ and $\zeta =\mathrm{sin}I\;\mathrm{exp}({\imath }{\rm{\Omega }})$, where e is the eccentricity, ϖ the longitude of pericenter, I the inclination, and Ω the longitude of the node. For (ii) we used the frequency modified Fourier transform by Šidlichovský & Nesvorný (1996). We identified planetary terms and the proper term, and computed the amplitude of the latter. For a given asteroid, we performed this computation for each of its clones and each of the five segments in time. We then combined the results using the clone statistical weights given by its $p({\boldsymbol{E}})$ value. This way we obtained mean values of the proper elements a_P, e_P, and sin I_P. Individual values of each clone at each of its segments served to estimate statistical variance that we report as the uncertainty of each of the proper values. We note that typically the uncertainty in a_P is dominated by the uncertainty of the initial osculating semimajor axis, especially for the less constrained orbits, while the uncertainties in e_P and sin I_P are mainly given by the variance in their computation at different time segments.

Figure 4 summarizes the results. We divided the family zone into three zones/boxes A, B, and C. Asteroids in A have orbits very similar to (2384) Schulhof, both in proper and osculating orbital elements. They constitute the immediate vicinity of this largest member in the family and contain multi-opposition objects 2007 EV68, 2008 RA126, 2009 EL11, and a single-opposition member 2015 GH4. Group B is what Pravec & Vokrouhlický (2009) thought to be family around (81337) 2000 GP36. It is composed of a compact cluster of numbered members 81337, 271044, and 286239 and a single-opposition object 2015 FN344. As of today, group B contains the second to fourth largest members in the family (Table 2), though we cannot exclude that other members of a comparable size are yet to be discovered. One of the features that puzzled Vokrouhlický & Nesvorný (2011) was the division of the family into the A and B zones. No objects were known in region C at that time. At this moment, however, we start to detect some Schulhof members even in this zone: a multi-opposition asteroid 2008 EK72 and a single-opposition asteroid 2012 FM46. It is important that 2008 EK72 exhibits the same convergence pattern to (2384) Schulhof as all numbered asteroids in box B (Section 4.1).

Zoom In Zoom Out Reset image size

As for the convergence results reported in the next section, it is interesting to note that asteroids in boxes B and C have opposite initial differences in secular angles with respect to (2384) Schulhof: those in box B have δΩ > 1° and $\delta \varpi \lt -1^\circ$, while those in box C have δΩ < −1° and δ ϖ > 1°. As mentioned above, asteroids in box A have orbits that are very similar to (2384) Schulhof, including the osculating values of the secular angles. In fact, this correlation between the secular angles and position in the family is an indication of their common origin.

We start discussing our results by first considering the case of the most accurately determined orbits, namely numbered asteroids 81337, 271044, and 286239 from zone B and 2008 EK72 from zone C (see Figure 4 and Table 2). All of these objects have a very good convergence to (2384) Schulhof at about the same time, some 800 kyr ago. An example for 81337 is shown in Figures 5 and 6, while the cases for other three asteroids are very similar. This is shown in a less detailed way in the left column of Figures 7 and 8.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Focusing first on our initial method, convergence of secular angles, we constructed a cumulative distribution of successful convergent trials with δV ≤ 5 m s⁻¹, and determined time instants at 10%, 50% (median), and 90% levels. We obtained the following estimates for convergence to (2384) Schulhof: (i) (81337) 2000 GP36 at ${765}_{-80}^{+140}$ kyr, (ii) (271044) 2003 FK6 at ${745}_{-70}^{+140}$ kyr, (iii) (286239) 2001 UR193 at ${810}_{-180}^{+220}$ kyr, and (iv) 2008 EK72 at ${725}_{-140}^{+140}$ kyr. These intervals largely overlap. In order to adopt a conservative approach to combine these estimates, we evaluate the mean value of the medians but keep the extreme minimum age for the 10% level and extreme maximum age for the 90% level over all cases. This provides an age estimate of ${760}_{-190}^{+270}$ kyr, in agreement with findings in Vokrouhlický & Nesvorný (2011).

Next, we examined our second method, the stronger convergence criterion, namely, the true encounter of the clone variants in physical space at a very low relative speed. Figure 6 shows the exemplary case of (81337) 2000 GP36's affinity to (2384) Schulhof. It is important to find that the two bodies have a possibility of very close encounters and that they occur within the well-defined interval of time previously hinted at by the secular angle convergence. While the minimum v of the order of 0.5 m s⁻¹ or less is frequently met, d reaches below the Hill radius only occasionally. This could be because of a slight incompleteness of our force model, the 5 year output frequency (in fact longer than the orbital period of Schulhof members), or a small number of clone variants used. Using conservative criteria d ≤ 4500 km (i.e., 1.5 times the Hill radius) and v ≤ 7 m s⁻¹, we collect a number of successful trials (top panel of Figure 6). Turning the distribution to a cumulative form and determining the 10%, 50% (median), and 90% levels as above, we obtain ${835}_{-123}^{+52}$ kyr. Similar results are also obtained for other asteroids in this class with ages of (i) (271044) 2003 FK6 at ${765}_{-90}^{+150}$ kyr, (ii) (286239) 2001 UR193 at ${910}_{-200}^{+75}$ kyr, and (iii) 2008 EK72 at ${785}_{-150}^{+75}$ kyr. We again collected these results into a conservative estimate of the family's age: ${825}_{-190}^{+160}$ kyr.

Additionally, we probed a simultaneous convergence of the secular angles of all asteroids with well-determined orbits considered in this section. Note that this was the original tool used by Nesvorný et al. (2006) and Nesvorný & Vokrouhlický (2006) to estimate the ages of their discovered young clusters. In this simulation, we used 11 geometric and 21 Yarkovsky clones for (2384) Schulhof, still employing the obliquity constraint from our solution in Section 2. All other smaller members were represented by 21 geometric and 41 Yarkovsky clones, this time allowing both positive and negative drift values in semimajor axis by the thermal accelerations. We probed their orbital evolution to 1.5 Myr ago, storing their state vectors every 5 years. At each of the output epochs we considered 25 million trial identifications between the clones of each of the five asteroids and computed the node and pericenter dispersal δΩ and $\delta \varpi$ using the formula given in Nesvorný & Vokrouhlický (2006). These were inserted into the velocity dispersal estimate by Equation (1). We required δV ≤ 5 m s⁻¹ to record a successful convergence. Figure 9 shows a normalized statistical distribution of these converging solutions binned into 25 kyr intervals of time. As expected, the resulting distribution roughly represents a combined intersection of converging distributions shown in the left column of Figure 7. The median, 10%, and 90% of the distribution provide yet another age estimate for the cluster of ${740}_{-50}^{+110}$ kyr.

Zoom In Zoom Out Reset image size

Finally, we comment on the Yarkovsky drift values da/dt for secondary components that exhibit successful convergence to (2384) Schulhof. In an ideal case, their distribution would be single-valued, or indicate a strong preference for a small interval of da/dt values. In that case, our results would also predict the sense of rotation of these smaller members in the family. Unfortunately, we find that the da/dt values are flat (randomly distributed). In other words, the secondary clones with different da/dt values have the possibility to converge at different time instants in the past, thus prohibiting us from constraining their rotation parameters. The same conclusion also applies to our runs in the next paragraphs.

We now turn to the case of slightly less accurate, but still reliable, multi-opposition orbits that were not considered in the previous section. These are 2007 EV68, 2008 RA126, and 2009 EL11, all located in box A in Figure 4. This means that all have initial orbits very similar to (2384) Schulhof. We note that the cases of 2007 EV68 and 2009 EL11 were already analyzed by Vokrouhlický & Nesvorný (2011), by then single-opposition asteroids, and they seemed to indicate preferential early convergence to (2384) Schulhof (age <100 kyr).

Results from our first method are shown in the middle column of Figure 7. Both 2007 EV68 and 2009 EL11 still offer young solutions. However, we note that this is very likely a bias of the convergence method that uses just the proximity of secular angles. In fact, their initial values provide an already quite low δV, and their short-period oscillations promptly result in fulfillment of the δV ≤ 7 m s⁻¹ condition. So these early solutions are not fully diagnostic. In this respect, we point to the rising "hump feature" in the number of convergent solutions of 2009 EL11 around 660 kyr (and also to a lesser degree in the case of 2007 EV68). Apparently, the orbit of 2008 RA126 is slightly farther from (2384) Schulhof, and the short-period oscillations of the secular angles do not permit these early solutions. Instead, only solutions from ∼500 kyr appear possible, with many in the formal interval of solutions determined from analysis of the accurate orbits in Section 4.1.

Our conclusions are neatly confirmed when using the second method, thus there is complete convergence in the Cartesian space. In this case both 2007 EV68 and 2009 EL11 show clear preference to ages >500 kyr, with most in the interval derived from the accurate orbits in Section 4.1. This is also confirmed by the solution for 2008 RA126.

We may thus conclude that all multi-opposition asteroids present a consistent picture of the Schulhof family age.

Finally, we consider the least accurate orbits, namely, single-opposition ones. Incidentally, we take one example from each of the boxes A, B, and C from Figure 4 (gray symbols): 2012 FM46 (box C), 2015 FN344 (box B), and 2015 GH4 (box A).

Figure 7 (last column) shows the results of our first method based on the behavior of the secular angles. As expected, the orbits starting from boxes B and C do not exhibit early age solutions. This is because their initial orbits, despite being quite uncertain, are distant enough from (2384) Schulhof that some period of time is needed to build enough changes in Ω and ϖ to approach the values for the largest member in the family. Interestingly, even the nearby (box A) single-opposition orbit of 2015 GH4 currently indicates the same pattern, and the statistical distribution of its successful solutions in time is similar to 2008 RA126. Slightly anomalous to the group is the case of 2012 FM46, which has an earlier peak and a more strongly damped tail than the distribution of successful solutions.

Results from the second convergence method are shown in Figure 8 (last column). The solution for the best constrained orbit in this class, 2015 FN344, fits well the age determined from the best orbits in box B. Data from 2015 GH4 are shifted to a slightly younger age, but this feature is possible to attribute to its currently poor orbital solution. This is because, by taking a constant number of clones for both objects, the method contains a bias toward detecting younger ages as the cloud of clones is quickly dispersing in time. This produces a characteristic $\propto {t}^{-2}$ tail in the number of successful solutions (e.g., J. Žižka et al. 2016, in preparation). While not dwelling on this issue in more detail, we consider solution for 2015 GH4 consistent with a ≃(700–800) kyr age of the family. On the contrary, results from 2012 FM46 remain somewhat anomalous. The formal solution is ${360}_{-100}^{+180}$ kyr and they decay faster than $\propto {t}^{-2}$ toward older age solutions. Only 5% of cases are convergent beyond 700 kyr ago. Given the experience with the previous study and the possibly misleading information from the analysis of the single-opposition objects (note that 2012 FM46 is the least accurate orbit in its class; Table 2), we refrain from drawing bold conclusions from this single case. More details will be obtained if the body is recovered during its next opposition in 2016 February or March.